diff --git a/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx b/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx index 4abfd65310a..6f8e60f981c 100644 --- a/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx +++ b/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx @@ -38,13 +38,16 @@ #include #include #include +#include #include #include #include +#include #include #include +#include #include #include @@ -94,7 +97,7 @@ enum EnEventHistograms { eEventHistograms_N }; -const char* eventHistNames[eEventHistograms_N] = { +static constexpr std::array EventHistNames = { "Centrality", "Multiplicity", "VertexX", @@ -108,7 +111,7 @@ enum EnParticleHistograms { eParticleHistograms_N }; -const char* particleHistNames[eParticleHistograms_N] = { +static constexpr std::array ParticleHistNames = { "Pt", "Phi"}; @@ -124,7 +127,7 @@ enum EnCorrHistograms { eCorrHistograms_N }; -const char* corrHistNames[eCorrHistograms_N] = { +static constexpr std::array CorrHistNames = { "Centrality", "Multiplicity", }; @@ -136,7 +139,7 @@ enum EnCentEstm { eCentEstm_N }; -const char* centEstmNames[eCentEstm_N] = { +static constexpr std::array CentEstmNames = { "FT0C", "FT0M", "FV0A"}; @@ -148,7 +151,7 @@ enum EnMultEstm { eMultEstm_N }; -const char* multEstmNames[eMultEstm_N] = { +static constexpr std::array MultEstmNames = { "FT0C", "FT0M", "FV0A"}; @@ -159,7 +162,7 @@ enum EnCutBeforeAfter { eCutBeforeAfter_N }; -const char* cutBeforeAfterNames[eCutBeforeAfter_N] = { +static constexpr std::array CutBeforeAfterNames = { "before", "after", }; @@ -172,8 +175,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam OutputObj fBaseList{sBaseListName.Data(), OutputObjHandlingPolicy::AnalysisObject, OutputObjSourceType::OutputObjSource}; // *) Service: - Service ccdb; // support for offline callibration data base - Service pdg; + Service ccdb{}; // support for offline callibration data base + Service pdg{}; // *) Define configurables: Configurable cfDryRun{"cfDryRun", false, "book all histos and run without filling and calculating anything"}; @@ -217,7 +220,6 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // *) Others Configurable cfFileWithWeights{"cfFileWithWeights", "/scratch3/go52dab/O2tutorial/tutorial3-6/weights.root", "path to external ROOT file which holds all particle weights in O2 format"}; - Configurable cfRunNumber{"cfRunNumber", "000123456", "run number"}; // *) Bins Configurable> cfPtBins{"cfPtBins", {1000, 0., 100.}, "nPtBins, ptMin, ptMax"}; @@ -232,8 +234,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // *) Define and initialize all data members to be called in the main process* functions: // **) Task configuration: struct TaskConfiguration { - bool fProcess[eProcess_N] = {false}; // Set what to process. See enum EnProcess for full description. Set via implicit variables within a PROCESS_SWITCH clause. - bool fDryRun = false; // book all histos and run without filling and calculating anything + std::array fProcess = {{false}}; // Set what to process. See enum EnProcess for full description. Set via implicit variables within a PROCESS_SWITCH clause. + bool fDryRun = false; std::string fCentEstm = "FT0M"; std::string fMultEstm = "FV0A"; @@ -277,68 +279,85 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam std::vector fVerZBins = {-50., 50.}; std::vector fNumContribBins = {0, 5000}; - std::string fFileWithWeights = "/scratch3/go52dab/O2tutorial/tutorial3-6/weights.root"; - std::string fRunNumber = "000123456"; + std::string fFileWithWeights = "/scratch3/go52dab/O2tutorial/analysis_code/weights.root"; } tc; struct ParticleHistograms { - TList* fParticleHistogramsList = NULL; - TH1F* fParticleHistograms[eParticleHistograms_N][2][2] = {{{NULL}}}; + TList* fParticleHistogramsList = nullptr; + std::array, 2>, eParticleHistograms_N> fParticleHistograms = {{{{nullptr}}}}; } pc; struct EventHistograms { - TList* fEventHistogramsList = NULL; - TH1F* fEventHistograms[eEventHistograms_N][2][2] = {{{NULL}}}; // [ type - see enum EnEventHistograms ][reco,sim][before,after event cuts] + TList* fEventHistogramsList = nullptr; + std::array, 2>, eEventHistograms_N> fEventHistograms = {{{{nullptr}}}}; } ev; struct QAHistograms { bool fQASwitch = kTRUE; - TList* fQAHistogramsList = NULL; - TH2F* fQAHistograms[eQAHistograms_N][2] = {{NULL}}; // [type][before/after cut] + TList* fQAHistogramsList = nullptr; + std::array, eQAHistograms_N> fQAHistograms = {{{nullptr}}}; // [type][before/after cut] } qa; struct CorrHistograms { - TList* fCorrHistogramsList = NULL; - TH2F* fCorrHistograms[eCorrHistograms_N][eMultEstm_N][eMultEstm_N][2] = {{{{NULL}}}}; + TList* fCorrHistogramsList = nullptr; + std::array, eMultEstm_N>, eMultEstm_N>, eCorrHistograms_N> fCorrHistograms = {{{{nullptr}}}}; } cr; struct WeightHistograms { bool fWeightSwitch = kTRUE; - TList* fWeightHistogramsList = NULL; - std::vector fWeightHistograms; + TList* fWeightHistogramsList = nullptr; + // Fill phi/pt histograms with that run number: + std::map fPtRealByRunMap; // MC rec (too lazy to change name) data pt histograms, valid if processMonteCarlo + std::map fPtMCByRunMap; // MC sim (too lazy to change name) data pt histograms, valid if processMonteCarlo + std::map fPhiByRunMap; // Phi histograms, valid if processRealData + // Make weight histograms locally. Upload weight histograms to CCDB: + std::vector fWeightHistograms; // Get all weight histograms with that run number + std::map fPhiWeightHistogramsMap; // Get phi weight histograms + std::map fPtWeightHistogramsMap; // Get pt weight histograms + // Null weight histograms. Use them when no weight histograms found in the given run number: + TH1F* fDummyPhiWeightHistogram = nullptr; + TH1F* fDummyPtWeightHistogram = nullptr; } wt; struct EventByEventQuantities { + int fRunNumber = 0; float fReferenceMultiplicity = 0.; float fCentrality = 0.; float fCentralitySim = 0.; float fImpactParameter = 0.; float fNumContrib = 0.; - - float fTwoParticleCorrelationEbye[eCutBeforeAfter_N][3] = {{0., 0., 0.}, {0., 0., 0.}}; //<2> [before, after][v2^2, v3^2, v4^2] - float fFourParticleCorrelationEbye[eCutBeforeAfter_N][2] = {{0., 0.}, {0., 0.}}; //<4> [before, after][v3^2 v2^2, v4^2 v2^2] + std::array, eCutBeforeAfter_N> fTwoParticleCorrelationEbye = {{{0., 0., 0.}, {0., 0., 0.}}}; //<2> [before, after][v2^2, v3^2, v4^2] + std::array, eCutBeforeAfter_N> fFourParticleCorrelationEbye = {{{0., 0.}, {0., 0.}}}; //<4> [before, after][v3^2 v2^2, v4^2 v2^2] } ebye; struct MultiparticleCorrelationProfile { - TList* fMultiparticleCorrelationProfilesList = NULL; - TProfile* fTwoParticleCorrelationProfiles[eCutBeforeAfter_N] = {NULL}; - TProfile* fFourParticleCorrelationProfiles[eCutBeforeAfter_N] = {NULL}; + TList* fMultiparticleCorrelationProfilesList = nullptr; + std::map fMultiparticleCorrelationByRunMap; + std::array fTwoParticleCorrelationProfiles = {{nullptr}}; + std::array fFourParticleCorrelationProfiles = {{nullptr}}; } mc; struct MultiparticleCorrelationCalculation { - int h1, h2, h3, h4, h5, h6, h7, h8; + int h1 = 0; + int h2 = 0; + int h3 = 0; + int h4 = 0; + int h5 = 0; + int h6 = 0; + int h7 = 0; + int h8 = 0; // Book Q-vector components: static constexpr int MaxCorrelator = 4; // <> static constexpr int MaxHarmonic = 5; // n+1 in vn, in this case n=4 as we need v2, v3, v4 static constexpr int MaxPower = MaxCorrelator + 1; static constexpr int NumSC = 2; // need SC(3,2) and SC(4,2) - TComplex fQvectorBefore[MaxHarmonic][MaxPower]; - TComplex fQvectorAfter[MaxHarmonic][MaxPower]; // All needed Q-vector components + std::array, MaxHarmonic> fQvectorBefore; + std::array, MaxHarmonic> fQvectorAfter; // All needed Q-vector components } mcc; template - bool ctEventCuts(T1 const& collision, T2 rlCollisionCentAll, T3 rlCollisionMultAll) + bool ctEventCuts(T1 const& collision, T2 const& rlCollisionCentAll, T3 const& rlCollisionMultAll) { bool pass = true; @@ -394,22 +413,22 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // *) Combine all switches if (tc.fVertexZCutSwitch) { - pass = pass && bVertexZCut; + pass &= bVertexZCut; } if (tc.fSel8CutSwitch) { - pass = pass && bSel8Cut; + pass &= bSel8Cut; } if (tc.fCentCutSwitch) { - pass = pass && bCentCut; + pass &= bCentCut; } if (tc.fNumContribCutSwitch) { - pass = pass && bNumContribCut; + pass &= bNumContribCut; } if (tc.fCentCorrCutSwitch) { - pass = pass && bCentCorrCut; + pass &= bCentCorrCut; } if (tc.fMultCorrCutSwitch) { - pass = pass && bMultCorrCut; + pass &= bMultCorrCut; } return pass; @@ -460,22 +479,22 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } if (tc.fPtCutSwitch) { - pass = pass && bPtCut; + pass &= bPtCut; } if (tc.fEtaCutSwitch) { - pass = pass && bEtaCut; + pass &= bEtaCut; } if (tc.fSignCutSwitch) { - pass = pass && bSignCut; + pass &= bSignCut; } if (tc.fTpcNClsFoundCutSwitch) { - pass = pass && bTpcNClsFoundCut; + pass &= bTpcNClsFoundCut; } if (tc.fDCAXYCutSwitch) { - pass = pass && bDCAXYCut; + pass &= bDCAXYCut; } if (tc.fDCAZCutSwitch) { - pass = pass && bDCAZCut; + pass &= bDCAZCut; } return pass; @@ -488,29 +507,30 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam if (eba == eBefore) { if (n >= 0) { return mcc.fQvectorBefore[n][p]; - } else { - return TComplex::Conjugate(mcc.fQvectorBefore[-n][p]); } + return TComplex::Conjugate(mcc.fQvectorBefore[-n][p]); } else { if (n >= 0) { return mcc.fQvectorAfter[n][p]; - } else { - return TComplex::Conjugate(mcc.fQvectorAfter[-n][p]); } + return TComplex::Conjugate(mcc.fQvectorAfter[-n][p]); } } - TComplex mccRecursion(int n, int* harmonic, EnCutBeforeAfter eba, int mult = 1, int skip = 0) + template + TComplex mccRecursion(int n, std::array harmonic, EnCutBeforeAfter eba, int mult = 1, int skip = 0) { // Calculate multi-particle correlators by using recursion (an improved faster version) originally developed by Kristjan Gulbrandsen (gulbrand@nbi.dk). int nm1 = n - 1; TComplex c(mccQ(harmonic[nm1], mult, eba)); - if (nm1 == 0) + if (nm1 == 0) { return c; + } c *= mccRecursion(nm1, harmonic, eba); - if (nm1 == skip) + if (nm1 == skip) { return c; + } int multp1 = mult + 1; int nm2 = n - 2; @@ -533,8 +553,9 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam harmonic[nm2] = harmonic[counter1]; harmonic[counter1] = hhold; - if (mult == 1) + if (mult == 1) { return c - c2; + } return c - static_cast(mult) * c2; } @@ -557,7 +578,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__); } if (0 == list->GetEntries()) { - return NULL; + return nullptr; } // The object is in the current base list: @@ -567,18 +588,19 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } // Otherwise, search for the object recursively in the nested lists: - TObject* objectIter; // iterator object in the loop below + TObject* objectIter = nullptr; // iterator object in the loop below TIter next(list); while ((objectIter = next())) // double round braces are to silence the warnings { if (TString(objectIter->ClassName()).EqualTo("TList")) { - objectFinal = getObjectFromList(reinterpret_cast(objectIter), objectName); - if (objectFinal) + objectFinal = getObjectFromList(dynamic_cast(objectIter), objectName); + if (objectFinal) { return objectFinal; + } } } // while(objectIter = next()) - return NULL; + return nullptr; } // TObject* getObjectFromList(TList* list, char* objectName) @@ -586,8 +608,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam { // a) Return value: std::vector histograms; - TList* baseList = NULL; // base top-level list in the TFile, e.g. named "ccdb_object" - TList* listWithRuns = NULL; // nested list with run-wise TList's holding run-specific weights + TList* baseList = nullptr; // base top-level list in the TFile, e.g. named "ccdb_object" + TList* listWithRuns = nullptr; // nested list with run-wise TList's holding run-specific weights // c) Determine from filePath if the file in on a local machine, or in home dir AliEn, or in CCDB: // Algorithm: @@ -607,7 +629,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam if (bFileIsInAliEn) { // File you want to access is in your home dir in AliEn: - TGrid* alien = TGrid::Connect("alien", gSystem->Getenv("USER"), "", ""); // do not forget to add #include to the preamble of your analysis task + TGrid* const alien = TGrid::Connect("alien", gSystem->Getenv("USER"), "", ""); // do not forget to add #include to the preamble of your analysis task if (!alien) { LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__); } @@ -617,19 +639,20 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } weightsFile->GetObject("ccdb_object", baseList); if (!baseList) { - weightsFile->ls(); LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__); } // Finally, from the top-level TList, get the desired nested TList => the technical problem here is that it can be nested at any level, // for thare there is a helper utility function GetObjectFromList(...) , see its implementation further below - listWithRuns = reinterpret_cast(getObjectFromList(baseList, runNumber)); + listWithRuns = dynamic_cast(getObjectFromList(baseList, runNumber)); if (!listWithRuns) { TString runNumberWithLeadingZeroes = "000"; runNumberWithLeadingZeroes += runNumber; // another try, with "000" prepended to run number - listWithRuns = reinterpret_cast(getObjectFromList(baseList, runNumberWithLeadingZeroes.Data())); + listWithRuns = dynamic_cast(getObjectFromList(baseList, runNumberWithLeadingZeroes.Data())); if (!listWithRuns) { - LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__); + LOGF(warning, "\033[1;31m%s at line %d : this crash can happen if in the output file there is no list with weights for the current runnumber = %s\033[0m", __FUNCTION__, __LINE__, runNumber); + histograms = {nullptr}; + return histograms; } } } else if (bFileIsInCCDB) { @@ -637,19 +660,20 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // Remember that here I do not access the file; instead, I directly access the object in that file. // My home dir in CCDB: https://alice-ccdb.cern.ch/browse/Users/a/abilandz/ => adapt for your case ccdb->setURL("https://alice-ccdb.cern.ch"); // to be able to use "ccdb" this object in your analysis task, see 4b/ below - baseList = reinterpret_cast(ccdb->get(TString(filePath).ReplaceAll("/alice-ccdb.cern.ch/", "").Data())); - baseList->ls(); + baseList = dynamic_cast(ccdb->get(TString(filePath).ReplaceAll("/alice-ccdb.cern.ch/", "").Data())); if (!baseList) { LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__); } - listWithRuns = reinterpret_cast(getObjectFromList(baseList, runNumber)); + listWithRuns = dynamic_cast(getObjectFromList(baseList, runNumber)); if (!listWithRuns) { TString runNumberWithLeadingZeroes = "000"; runNumberWithLeadingZeroes += runNumber; // another try, with "000" prepended to run number - listWithRuns = reinterpret_cast(getObjectFromList(baseList, runNumberWithLeadingZeroes.Data())); + listWithRuns = dynamic_cast(getObjectFromList(baseList, runNumberWithLeadingZeroes.Data())); if (!listWithRuns) { - LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__); + LOGF(warning, "\033[1;31m%s at line %d : this crash can happen if in the output file there is no list with weights for the current runnumber = %s\033[0m", __FUNCTION__, __LINE__, runNumber); + histograms = {nullptr}; + return histograms; } } @@ -661,8 +685,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // Check if the external ROOT file exists at specified path: if (gSystem->AccessPathName(filePath, kFileExists)) { - LOGF(info, - "\033[1;33m if(gSystem->AccessPathName(filePath,kFileExists)), filePath = %s \033[0m", filePath); + LOGF(info, "\033[1;33m if(gSystem->AccessPathName(filePath, kFileExists)), filePath = %s \033[0m", filePath); LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__); } @@ -674,18 +697,18 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam weightsFile->GetObject("ccdb_object", baseList); if (!baseList) { - weightsFile->ls(); LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__); } - listWithRuns = reinterpret_cast(getObjectFromList(baseList, runNumber)); + listWithRuns = dynamic_cast(getObjectFromList(baseList, runNumber)); if (!listWithRuns) { TString runNumberWithLeadingZeroes = "000"; runNumberWithLeadingZeroes += runNumber; // another try, with "000" prepended to run number - listWithRuns = reinterpret_cast(getObjectFromList(baseList, runNumberWithLeadingZeroes.Data())); + listWithRuns = dynamic_cast(getObjectFromList(baseList, runNumberWithLeadingZeroes.Data())); if (!listWithRuns) { - baseList->ls(); - LOGF(fatal, "\033[1;31m%s at line %d : this crash can happen if in the output file there is no list with weights for the current runnumber = %s\033[0m", __FUNCTION__, __LINE__, runNumber); + LOGF(warning, "\033[1;31m%s at line %d : this crash can happen if in the output file there is no list with weights for the current runnumber = %s\033[0m", __FUNCTION__, __LINE__, runNumber); + histograms = {nullptr}; + return histograms; } } } @@ -708,7 +731,9 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam hist->SetDirectory(0); auto* histClone = dynamic_cast(hist->Clone()); if (!histClone) { - LOGF(fatal, "Failed to clone histogram %s", hist->GetName()); + LOGF(warning, "Failed to clone histogram %s", hist->GetName()); + histograms = {nullptr}; + return histograms; } histClone->SetTitle(Form("%s:%s", filePath, histClone->GetName())); histograms.push_back(histClone); @@ -727,8 +752,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam return; } - // TH1F* histAcceptanceWeight = wt.fWeightHistograms[1]; - + // Book Q-vector arrays: for (int h = 0; h < mcc.MaxHarmonic; h++) { for (int p = 0; p < mcc.MaxPower; p++) { mcc.fQvectorBefore[h][p] = TComplex(0., 0.); @@ -736,12 +760,114 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } } - if (tc.fPrintSwitch) { - // Print current run number: - LOGF(info, "Run number: %d", collision.bc().runNumber()); + // Get run number: + ebye.fRunNumber = collision.bc().runNumber(); + std::string stringRunNumber = std::to_string(ebye.fRunNumber); + + // Book phi histogram with this run number: + if (!wt.fPhiByRunMap.contains(ebye.fRunNumber)) { + // wt.fPhiByRunMap[ebye.fRunNumber] = new TH1F(Form("hPhi_run%d", ebye.fRunNumber), Form("phi distribution for run %d", ebye.fRunNumber), static_cast(tc.fPhiBins[0]), tc.fPhiBins[1], tc.fPhiBins[2]); + // wt.fPhiByRunMap[ebye.fRunNumber]->SetDirectory(nullptr); + // wt.fWeightHistogramsList->Add(wt.fPhiByRunMap[ebye.fRunNumber]); + auto* hPhi = new TH1F(Form("hPhi_run%d", ebye.fRunNumber), Form("phi distribution for run %d", ebye.fRunNumber), static_cast(tc.fPhiBins[0]), tc.fPhiBins[1], tc.fPhiBins[2]); + hPhi->SetDirectory(nullptr); + wt.fPhiByRunMap.try_emplace(ebye.fRunNumber, hPhi); + wt.fWeightHistogramsList->Add(wt.fPhiByRunMap[ebye.fRunNumber]); + } + + // Book pt MC rec histogram with this run number: + if (!wt.fPtRealByRunMap.contains(ebye.fRunNumber)) { + // wt.fPtRealByRunMap[ebye.fRunNumber] = new TH1F(Form("hPtReal_run%d", ebye.fRunNumber), Form("pt MC rec distribution for run %d", ebye.fRunNumber), static_cast(tc.fPtBins[0]), tc.fPtBins[1], tc.fPtBins[2]); + // wt.fPtRealByRunMap[ebye.fRunNumber]->SetDirectory(nullptr); + // wt.fWeightHistogramsList->Add(wt.fPtRealByRunMap[ebye.fRunNumber]); + auto* hPtReal = new TH1F(Form("hPtReal_run%d", ebye.fRunNumber), Form("pt MC rec distribution for run %d", ebye.fRunNumber), static_cast(tc.fPtBins[0]), tc.fPtBins[1], tc.fPtBins[2]); + hPtReal->SetDirectory(nullptr); + wt.fPtRealByRunMap.try_emplace(ebye.fRunNumber, hPtReal); + wt.fWeightHistogramsList->Add(wt.fPtRealByRunMap[ebye.fRunNumber]); + } + + // Book pt MC sim histogram with this run number: + if (!wt.fPtMCByRunMap.contains(ebye.fRunNumber)) { + // wt.fPtMCByRunMap[ebye.fRunNumber] = new TH1F(Form("hPtMC_run%d", ebye.fRunNumber), Form("pt MC sim distribution for run %d", ebye.fRunNumber), static_cast(tc.fPtBins[0]), tc.fPtBins[1], tc.fPtBins[2]); + // wt.fPtMCByRunMap[ebye.fRunNumber]->SetDirectory(nullptr); + // wt.fWeightHistogramsList->Add(wt.fPtMCByRunMap[ebye.fRunNumber]); + auto* hPtMC = new TH1F(Form("hPtMC_run%d", ebye.fRunNumber), Form("pt MC sim distribution for run %d", ebye.fRunNumber), static_cast(tc.fPtBins[0]), tc.fPtBins[1], tc.fPtBins[2]); + hPtMC->SetDirectory(nullptr); + wt.fPtMCByRunMap.try_emplace(ebye.fRunNumber, hPtMC); + wt.fWeightHistogramsList->Add(wt.fPtMCByRunMap[ebye.fRunNumber]); + } + + // Get phi and pt weight histogram with this run number: + if (wt.fWeightSwitch && !wt.fPhiWeightHistogramsMap.contains(ebye.fRunNumber)) { + + TH1F* phiWeightHist = dynamic_cast(wt.fDummyPhiWeightHistogram->Clone(Form("wPhi_run%d", ebye.fRunNumber))); + TH1F* ptWeightHist = dynamic_cast(wt.fDummyPtWeightHistogram->Clone(Form("wPt_run%d", ebye.fRunNumber))); + + wt.fWeightHistograms = getHistogramsWithWeights(tc.fFileWithWeights.c_str(), stringRunNumber.c_str()); + + for (auto const& hist : wt.fWeightHistograms) { + if (!hist) { + LOGF(warning, "Fail to loop weight histograms"); + continue; + } + TString histName = hist->GetName(); + if (histName.BeginsWith("wPhi")) { + delete phiWeightHist; + phiWeightHist = dynamic_cast(hist->Clone(Form("wPhi_run%d", ebye.fRunNumber))); + break; + } + } + + for (auto const& hist : wt.fWeightHistograms) { + if (!hist) { + LOGF(warning, "Fail to loop weight histograms"); + continue; + } + TString histName = hist->GetName(); + if (histName.BeginsWith("wPt")) { + delete ptWeightHist; + ptWeightHist = dynamic_cast(hist->Clone(Form("wPt_run%d", ebye.fRunNumber))); + break; + } + } + + phiWeightHist->SetDirectory(nullptr); + wt.fPhiWeightHistogramsMap.try_emplace(ebye.fRunNumber, phiWeightHist); + wt.fWeightHistogramsList->Add(wt.fPhiWeightHistogramsMap[ebye.fRunNumber]); + + ptWeightHist->SetDirectory(nullptr); + wt.fPtWeightHistogramsMap.try_emplace(ebye.fRunNumber, ptWeightHist); + wt.fWeightHistogramsList->Add(wt.fPtWeightHistogramsMap[ebye.fRunNumber]); + } + + // Multiparticle correlation list with this run number: + if (!mc.fMultiparticleCorrelationByRunMap.contains(ebye.fRunNumber)) { + mc.fMultiparticleCorrelationByRunMap.try_emplace(ebye.fRunNumber, new TList()); + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->SetName(Form("mcc_run%d", ebye.fRunNumber)); + mc.fMultiparticleCorrelationProfilesList->Add(mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]); + + mc.fTwoParticleCorrelationProfiles[eBefore] = nullptr; + mc.fTwoParticleCorrelationProfiles[eAfter] = nullptr; + mc.fFourParticleCorrelationProfiles[eBefore] = nullptr; + mc.fFourParticleCorrelationProfiles[eAfter] = nullptr; + + mc.fTwoParticleCorrelationProfiles[eBefore] = new TProfile("prof2Before", "2-p correlation before cut", 3, 2., 5.); + mc.fTwoParticleCorrelationProfiles[eAfter] = new TProfile("prof2After", "2-p correlation after cut", 3, 2., 5.); + mc.fTwoParticleCorrelationProfiles[eBefore]->Sumw2(); + mc.fTwoParticleCorrelationProfiles[eAfter]->Sumw2(); + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fTwoParticleCorrelationProfiles[eBefore]); + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fTwoParticleCorrelationProfiles[eAfter]); + + mc.fFourParticleCorrelationProfiles[eBefore] = new TProfile("prof4Before", "4-p correlation before cut", 2, 3., 5.); + mc.fFourParticleCorrelationProfiles[eAfter] = new TProfile("prof4After", "4-p correlation after cut", 2, 3., 5.); + mc.fFourParticleCorrelationProfiles[eBefore]->Sumw2(); + mc.fFourParticleCorrelationProfiles[eAfter]->Sumw2(); + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fFourParticleCorrelationProfiles[eBefore]); + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fFourParticleCorrelationProfiles[eAfter]); } - float rlCollisionCentAll[eCentEstm_N] = { + // Real data centrality: + std::array rlCollisionCentAll = { collision.centFT0C(), collision.centFT0M(), collision.centFV0A()}; @@ -750,15 +876,16 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam for (int i = 0; i < eCentEstm_N; i++) { if (tc.fPrintSwitch) { - LOGF(info, "%s Centrality: %f", centEstmNames[i], + LOGF(info, "%s Centrality: %f", CentEstmNames[i], rlCollisionCentAll[i]); } - if (tc.fCentEstm == centEstmNames[i]) { + if (tc.fCentEstm == CentEstmNames[i]) { rlCollisionCent = rlCollisionCentAll[i]; } } - float rlCollisionMultAll[eMultEstm_N] = { + // Real data multiplicity: + std::array rlCollisionMultAll = { static_cast(collision.multFT0C()), static_cast(collision.multFT0M()), static_cast(collision.multFV0A())}; @@ -766,37 +893,41 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam for (int i = 0; i < eMultEstm_N; i++) { if (tc.fPrintSwitch) { - LOGF(info, "%s Multiplicity: %f", multEstmNames[i], rlCollisionMultAll[i]); + LOGF(info, "%s Multiplicity: %f", MultEstmNames[i], rlCollisionMultAll[i]); } - if (tc.fMultEstm == multEstmNames[i]) { + if (tc.fMultEstm == MultEstmNames[i]) { rlCollisionMult = rlCollisionMultAll[i]; } } + // Real data nContrib: float rlCollisionNumContrib = 0.; rlCollisionNumContrib = static_cast(collision.numContrib()); + // Event-by-event quantity: + ebye.fCentrality = rlCollisionCent; + ebye.fReferenceMultiplicity = rlCollisionMult; + ebye.fNumContrib = rlCollisionNumContrib; + + // Print... if (tc.fPrintSwitch) { - // Print centrality estimated with "FT0M" estimator: - LOGF(info, "Centrality: %f", rlCollisionCent); - // Print multiplicity: + LOGF(info, "Run number: %d", ebye.fRunNumber); + + LOGF(info, "Centrality: %f", rlCollisionCent); LOGF(info, "Multiplicity: %f", static_cast(rlCollisionMult)); - // Print vertex position: LOGF(info, "Vertex X position: %f", collision.posX()); LOGF(info, "Vertex Y position: %f", collision.posY()); LOGF(info, "Vertex Z position: %f", collision.posZ()); - // Print NContributors - LOGF(info, "NContributors: %f", - static_cast(rlCollisionNumContrib)); + LOGF(info, "NContributors: %f", static_cast(rlCollisionNumContrib)); } - ebye.fCentrality = rlCollisionCent; - ebye.fReferenceMultiplicity = rlCollisionMult; - ebye.fNumContrib = rlCollisionNumContrib; + // If Rec or RecAndSim: if constexpr (rs == eRec || rs == eRecAndSim) { + + // Fill real event histograms before cut: ev.fEventHistograms[eCent][eRec][eBefore]->Fill(rlCollisionCent); ev.fEventHistograms[eMult][eRec][eBefore]->Fill(rlCollisionMult); ev.fEventHistograms[eVertexX][eRec][eBefore]->Fill(collision.posX()); @@ -804,6 +935,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam ev.fEventHistograms[eVertexZ][eRec][eBefore]->Fill(collision.posZ()); ev.fEventHistograms[eNumContrib][eRec][eBefore]->Fill(rlCollisionNumContrib); + // Fill centrality correlation histograms before cut: for (int i = 0; i < eCentEstm_N; i++) { for (int j = i + 1; j < eCentEstm_N; j++) { auto* h = cr.fCorrHistograms[eCorrCent][i][j][eBefore]; @@ -814,6 +946,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } } + // Fill multiplicity correlation histograms before cut: for (int i = 0; i < eMultEstm_N; i++) { for (int j = i + 1; j < eMultEstm_N; j++) { auto* h = cr.fCorrHistograms[eCorrMult][i][j][eBefore]; @@ -824,53 +957,59 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } } + // If RecAndSim: if constexpr (rs == eRecAndSim) { + if (!collision.has_mcCollision()) { if (tc.fPrintSwitch) { LOGF(warning, " No MC collision for this collision, skip..."); } - return; - } - - auto mccollision = collision.mcCollision(); + } else { + // Define MC collision: + auto mccollision = collision.mcCollision(); - float mcCollisionCent = 0.; + // Define MC centrality: + float mcCollisionCent = 0.; + float b = mccollision.impactParameter() * std::pow(10, -15); // convert fm to m + float xs = 7.71 * std::pow(10, -28); // convert barn to m^2 + mcCollisionCent = o2::constants::math::PI * b * b / xs * 100; - float b = mccollision.impactParameter() * std::pow(10, -15); // convert fm to m - float xs = 7.71 * std::pow(10, -28); // convert barn to m^2 - mcCollisionCent = o2::constants::math::PI * b * b / xs * 100; + // Event-by-event quantity: + ebye.fCentralitySim = mcCollisionCent; + ebye.fImpactParameter = b; - ebye.fCentralitySim = mcCollisionCent; - ebye.fImpactParameter = b; - - if (tc.fPrintSwitch) { - LOGF(info, "mc impact param (fm): %f", mccollision.impactParameter()); - LOGF(info, "mc centrality: %f", mcCollisionCent); - } - - ev.fEventHistograms[eCent][eSim][eBefore]->Fill(mcCollisionCent); - ev.fEventHistograms[eVertexX][eSim][eBefore]->Fill(mccollision.posX()); - ev.fEventHistograms[eVertexY][eSim][eBefore]->Fill(mccollision.posY()); - ev.fEventHistograms[eVertexZ][eSim][eBefore]->Fill(mccollision.posZ()); + if (tc.fPrintSwitch) { + LOGF(info, "mc impact param (fm): %f", mccollision.impactParameter()); + LOGF(info, "mc centrality: %f", mcCollisionCent); + } - if (ctEventCuts(mccollision, rlCollisionCentAll, rlCollisionMultAll)) { - ev.fEventHistograms[eCent][eSim][eAfter]->Fill(mcCollisionCent); - ev.fEventHistograms[eVertexX][eSim][eAfter]->Fill(mccollision.posX()); - ev.fEventHistograms[eVertexY][eSim][eAfter]->Fill(mccollision.posY()); - ev.fEventHistograms[eVertexZ][eSim][eAfter]->Fill(mccollision.posZ()); - } + // Fill MC event histograms before cut: + ev.fEventHistograms[eCent][eSim][eBefore]->Fill(mcCollisionCent); + ev.fEventHistograms[eVertexX][eSim][eBefore]->Fill(mccollision.posX()); + ev.fEventHistograms[eVertexY][eSim][eBefore]->Fill(mccollision.posY()); + ev.fEventHistograms[eVertexZ][eSim][eBefore]->Fill(mccollision.posZ()); + + // Fill MC event histograms after cut: + if (ctEventCuts(mccollision, rlCollisionCentAll, rlCollisionMultAll)) { + ev.fEventHistograms[eCent][eSim][eAfter]->Fill(mcCollisionCent); + ev.fEventHistograms[eVertexX][eSim][eAfter]->Fill(mccollision.posX()); + ev.fEventHistograms[eVertexY][eSim][eAfter]->Fill(mccollision.posY()); + ev.fEventHistograms[eVertexZ][eSim][eAfter]->Fill(mccollision.posZ()); + } - if (qa.fQASwitch) { - qa.fQAHistograms[eQACent][eBefore]->Fill(rlCollisionCent, mcCollisionCent); - qa.fQAHistograms[eQAMultNumContrib][eBefore]->Fill(rlCollisionMult, rlCollisionNumContrib); - if (ctEventCuts(collision, rlCollisionCentAll, rlCollisionMultAll) && - ctEventCuts(mccollision, rlCollisionCentAll, rlCollisionMultAll)) { - qa.fQAHistograms[eQACent][eAfter]->Fill(rlCollisionCent, mcCollisionCent); - qa.fQAHistograms[eQAMultNumContrib][eAfter]->Fill(rlCollisionMult, rlCollisionNumContrib); + if (qa.fQASwitch) { + qa.fQAHistograms[eQACent][eBefore]->Fill(rlCollisionCent, mcCollisionCent); + qa.fQAHistograms[eQAMultNumContrib][eBefore]->Fill(rlCollisionMult, rlCollisionNumContrib); + if (ctEventCuts(collision, rlCollisionCentAll, rlCollisionMultAll) && + ctEventCuts(mccollision, rlCollisionCentAll, rlCollisionMultAll)) { + qa.fQAHistograms[eQACent][eAfter]->Fill(rlCollisionCent, mcCollisionCent); + qa.fQAHistograms[eQAMultNumContrib][eAfter]->Fill(rlCollisionMult, rlCollisionNumContrib); + } } } } + // Fill real event histograms after cut if (ctEventCuts(collision, rlCollisionCentAll, rlCollisionMultAll)) { ev.fEventHistograms[eCent][eRec][eAfter]->Fill(rlCollisionCent); ev.fEventHistograms[eMult][eRec][eAfter]->Fill(rlCollisionMult); @@ -879,6 +1018,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam ev.fEventHistograms[eVertexZ][eRec][eAfter]->Fill(collision.posZ()); ev.fEventHistograms[eNumContrib][eRec][eAfter]->Fill(rlCollisionNumContrib); + // Fill centrality correlation histograms after cut: if (tc.fCentCorrCutSwitch) { for (int i = 0; i < eCentEstm_N; i++) { for (int j = i + 1; j < eCentEstm_N; j++) { @@ -891,8 +1031,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } } + // Fill multiplicity correlation histograms after cut: if (tc.fMultCorrCutSwitch) { - for (int i = 0; i < eMultEstm_N; i++) { for (int j = i + 1; j < eMultEstm_N; j++) { auto* h = cr.fCorrHistograms[eCorrMult][i][j][eAfter]; @@ -904,6 +1044,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } } + // Fail the event cut, skip this collision: } else { return; } @@ -914,25 +1055,42 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // Calculate Q-vectors for available angles and weights: double dPhi = 0.; // particle angle + double dPt = 0.; double wPhi = 1.; // particle weight - double wPhiToPowerP = 1.; // particle weight raised to power p + double wPt = 1.; + double wPhiToPowerP = 1.; // particle weight raised to power p. wPhi is actually wPt*wPhi but I'm too lazy to change the name. // Main loop over particles: for (auto const& track : tracks) { // LOGF(info, "Track azimuthal angle: %f", track.phi()); // LOGF(info, "Transverse momentum: %f", track.pt()); - // Fill reconstructed ...: if constexpr (rs == eRec || rs == eRecAndSim) { - // Fill track pt distribution: + // Fill phi/pt real histogram with this run number: + wt.fPhiByRunMap.at(ebye.fRunNumber)->Fill(track.phi()); + wt.fPtRealByRunMap.at(ebye.fRunNumber)->Fill(track.pt()); + + // Fill track histograms before cut: pc.fParticleHistograms[ePt][eRec][eBefore]->Fill(track.pt()); pc.fParticleHistograms[ePhi][eRec][eBefore]->Fill(track.phi()); + // Calculating Q-vector before cut: dPhi = track.phi(); + dPt = track.pt(); if (wt.fWeightSwitch) { - auto* hist = wt.fWeightHistograms[1]; - wPhi = hist->GetBinContent(wt.fWeightHistograms[1]->GetXaxis()->FindBin(dPhi)); + if (!wt.fPhiWeightHistogramsMap.contains(ebye.fRunNumber) || !wt.fPhiWeightHistogramsMap.at(ebye.fRunNumber)) { + LOGF(fatal, "Missing Phi weight histogram for run %d", ebye.fRunNumber); + } + if (!wt.fPtWeightHistogramsMap.contains(ebye.fRunNumber) || !wt.fPtWeightHistogramsMap.at(ebye.fRunNumber)) { + LOGF(fatal, "Missing Pt weight histogram for run %d", ebye.fRunNumber); + } + + auto* histPhi = wt.fPhiWeightHistogramsMap.at(ebye.fRunNumber); + auto* histPt = wt.fPtWeightHistogramsMap.at(ebye.fRunNumber); + wPhi = histPhi->GetBinContent(histPhi->GetXaxis()->FindBin(dPhi)); + wPt = histPt->GetBinContent(histPt->GetXaxis()->FindBin(dPt)); + wPhi *= wPt; } for (int h = 0; h < mcc.MaxHarmonic; h++) { @@ -945,9 +1103,12 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } if (ctParticleCuts(track)) { + + // Fill particle histograms after cut: pc.fParticleHistograms[ePt][eRec][eAfter]->Fill(track.pt()); pc.fParticleHistograms[ePhi][eRec][eAfter]->Fill(track.phi()); + // Calculating Q-vector after cut: for (int h = 0; h < mcc.MaxHarmonic; h++) { for (int p = 0; p < mcc.MaxPower; p++) { if (wt.fWeightSwitch) { @@ -971,26 +1132,32 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam if (tc.fPrintSwitch) { LOGF(warning, " No MC particle for this track, skip..."); } - return; - } - auto mcparticle = track.mcParticle(); // corresponding MC truth simulated particle - pc.fParticleHistograms[ePt][eSim][eBefore]->Fill(mcparticle.pt()); - pc.fParticleHistograms[ePhi][eSim][eBefore]->Fill(mcparticle.phi()); - if (ctParticleCuts(mcparticle)) { - pc.fParticleHistograms[ePt][eSim][eAfter]->Fill(mcparticle.pt()); - pc.fParticleHistograms[ePhi][eSim][eAfter]->Fill(mcparticle.phi()); + } else { + // Corresponding MC truth simulated particle + auto mcparticle = track.mcParticle(); + + // Fill pt MC sim histogram with this run number: + wt.fPtMCByRunMap.at(ebye.fRunNumber)->Fill(mcparticle.pt()); + + // Fill MC particle histograms before cut: + pc.fParticleHistograms[ePt][eSim][eBefore]->Fill(mcparticle.pt()); + pc.fParticleHistograms[ePhi][eSim][eBefore]->Fill(mcparticle.phi()); + + if (ctParticleCuts(mcparticle)) { + // Fill MC particle histograms after cut: + pc.fParticleHistograms[ePt][eSim][eAfter]->Fill(mcparticle.pt()); + pc.fParticleHistograms[ePhi][eSim][eAfter]->Fill(mcparticle.phi()); + } } } // end of if constexpr (rs == eRecAndSim) { - } // if constexpr (rs == eRec || rs == eRecAndSim) { - } // end of for (int64_t i = 0; i < tracks.size(); i++) { for (int i = 0; i < mcc.MaxHarmonic - 2; i++) { mcc.h1 = -(i + 2); mcc.h2 = i + 2; - int harmonicsTwoNum[2] = {mcc.h1, mcc.h2}; - int harmonicsTwoDen[2] = {0, 0}; + std::array harmonicsTwoNum = {mcc.h1, mcc.h2}; + std::array harmonicsTwoDen = {0, 0}; TComplex twoRecursionBefore = mccRecursion(2, harmonicsTwoNum, eBefore) / mccRecursion(2, harmonicsTwoDen, eBefore).Re(); double wTwoRecursionBefore = mccRecursion(2, harmonicsTwoDen, eBefore).Re(); @@ -1014,8 +1181,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam mcc.h2 = -2; mcc.h3 = 2; mcc.h4 = i + 3; - int harmonicsFourNum[4] = {mcc.h1, mcc.h2, mcc.h3, mcc.h4}; - int harmonicsFourDen[4] = {0, 0, 0, 0}; + std::array harmonicsFourNum = {mcc.h1, mcc.h2, mcc.h3, mcc.h4}; + std::array harmonicsFourDen = {0, 0, 0, 0}; TComplex fourRecursionBefore = mccRecursion(4, harmonicsFourNum, eBefore) / mccRecursion(4, harmonicsFourDen, eBefore).Re(); double wFourRecursionBefore = mccRecursion(4, harmonicsFourDen, eBefore).Re(); @@ -1035,70 +1202,70 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } template - void bookParticleHistograms(T1 const& lPcBins, ParticleHistograms& pc) + void bookParticleHistograms(T1 const& lPcBins/* , ParticleHistograms& pc */) { - std::vector lPtBins = lPcBins[histType]; // define local array and initialize it from an array set in the configurables + const auto& lPtBins = lPcBins[histType]; // define local array and initialize it from an array set in the configurables int nBinsPt = static_cast(lPtBins[0]); float minPt = lPtBins[1]; float maxPt = lPtBins[2]; for (int ba = 0; ba < eCutBeforeAfter_N; ba++) { - std::string nameRec = Form("fHist%s[eRec][%s cut]", particleHistNames[histType], cutBeforeAfterNames[ba]); - std::string nameSim = Form("fHist%s[eSim][%s cut]", particleHistNames[histType], cutBeforeAfterNames[ba]); - std::string nameRecfull = Form("%s distribution for reconstructed particles", particleHistNames[histType]); - std::string nameSimfull = Form("%s distribution for simulated particles", particleHistNames[histType]); + std::string nameRec = Form("fHist%s[eRec][%s cut]", ParticleHistNames[histType], CutBeforeAfterNames[ba]); + std::string nameSim = Form("fHist%s[eSim][%s cut]", ParticleHistNames[histType], CutBeforeAfterNames[ba]); + std::string nameRecfull = Form("%s distribution for reconstructed particles", ParticleHistNames[histType]); + std::string nameSimfull = Form("%s distribution for simulated particles", ParticleHistNames[histType]); if (doprocessRec || doprocessRecSim) { pc.fParticleHistograms[histType][eRec][ba] = new TH1F(nameRec.c_str(), nameRecfull.c_str(), nBinsPt, minPt, maxPt); - pc.fParticleHistograms[histType][eRec][ba]->GetXaxis()->SetTitle(particleHistNames[histType]); + pc.fParticleHistograms[histType][eRec][ba]->GetXaxis()->SetTitle(ParticleHistNames[histType]); pc.fParticleHistogramsList->Add(pc.fParticleHistograms[histType][eRec][ba]); } if (doprocessSim || doprocessRecSim) { pc.fParticleHistograms[histType][eSim][ba] = new TH1F(nameSim.c_str(), nameSimfull.c_str(), nBinsPt, minPt, maxPt); - pc.fParticleHistograms[histType][eSim][ba]->GetXaxis()->SetTitle(particleHistNames[histType]); + pc.fParticleHistograms[histType][eSim][ba]->GetXaxis()->SetTitle(ParticleHistNames[histType]); pc.fParticleHistogramsList->Add(pc.fParticleHistograms[histType][eSim][ba]); } } } template - void bookEventHistograms(T1 const& lEvBins, EventHistograms& ev) + void bookEventHistograms(T1 const& lEvBins/* , EventHistograms& ev */) { - std::vector lCentBins = lEvBins[histType]; // define local array and initialize it from an array set in the configurables + const auto& lCentBins = lEvBins[histType]; // define local array and initialize it from an array set in the configurables int nBinsCent = static_cast(lCentBins[0]); float minCent = lCentBins[1]; float maxCent = lCentBins[2]; for (int ba = 0; ba < eCutBeforeAfter_N; ba++) { - std::string nameRec = Form("fHist%s[eRec][%s cut]", eventHistNames[histType], cutBeforeAfterNames[ba]); - std::string nameSim = Form("fHist%s[eSim][%s cut]", eventHistNames[histType], cutBeforeAfterNames[ba]); + std::string nameRec = Form("fHist%s[eRec][%s cut]", EventHistNames[histType], CutBeforeAfterNames[ba]); + std::string nameSim = Form("fHist%s[eSim][%s cut]", EventHistNames[histType], CutBeforeAfterNames[ba]); std::string nameRecfull; std::string nameSimfull; if constexpr (histType == eCent) { - std::string nameRecfull = Form("%s %s distribution for reconstructed events", tc.fCentEstm.c_str(), eventHistNames[histType]); - std::string nameSimfull = Form("%s %s distribution for simulated events", tc.fCentEstm.c_str(), eventHistNames[histType]); + nameRecfull = Form("%s %s distribution for reconstructed events", tc.fCentEstm.c_str(), EventHistNames[histType]); + nameSimfull = Form("%s %s distribution for simulated events", tc.fCentEstm.c_str(), EventHistNames[histType]); } else if constexpr (histType == eMult) { - std::string nameRecfull = Form("%s %s distribution for reconstructed events", tc.fMultEstm.c_str(), eventHistNames[histType]); - std::string nameSimfull = Form("%s %s distribution for simulated events", tc.fMultEstm.c_str(), eventHistNames[histType]); + nameRecfull = Form("%s %s distribution for reconstructed events", tc.fMultEstm.c_str(), EventHistNames[histType]); + nameSimfull = Form("%s %s distribution for simulated events", tc.fMultEstm.c_str(), EventHistNames[histType]); } else { - std::string nameRecfull = Form("%s distribution for reconstructed events", eventHistNames[histType]); - std::string nameSimfull = Form("%s distribution for simulated events", eventHistNames[histType]); + nameRecfull = Form("%s distribution for reconstructed events", EventHistNames[histType]); + nameSimfull = Form("%s distribution for simulated events", EventHistNames[histType]); } if (doprocessRec || doprocessRecSim) { ev.fEventHistograms[histType][eRec][ba] = new TH1F(nameRec.c_str(), nameRecfull.c_str(), nBinsCent, minCent, maxCent); - ev.fEventHistograms[histType][eRec][ba]->GetXaxis()->SetTitle(eventHistNames[histType]); + ev.fEventHistograms[histType][eRec][ba]->GetXaxis()->SetTitle(EventHistNames[histType]); ev.fEventHistogramsList->Add(ev.fEventHistograms[histType][eRec][ba]); } if (doprocessSim || doprocessRecSim) { if constexpr (histType != eNumContrib && histType != eMult) { ev.fEventHistograms[histType][eSim][ba] = new TH1F(nameSim.c_str(), nameSimfull.c_str(), nBinsCent, minCent, maxCent); - ev.fEventHistograms[histType][eSim][ba]->GetXaxis()->SetTitle(eventHistNames[histType]); + ev.fEventHistograms[histType][eSim][ba]->GetXaxis()->SetTitle(EventHistNames[histType]); ev.fEventHistogramsList->Add(ev.fEventHistograms[histType][eSim][ba]); } // No nContrib and multiplicity for processSim } @@ -1106,7 +1273,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } template - void bookQAHistograms(T1 const& lQABins, QAHistograms& qa) + void bookQAHistograms(T1 const& lQABins/* , QAHistograms& qa */) { int nBinsCentX = 0; float minCentX = 0.; @@ -1119,31 +1286,31 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam float maxCent = 0.; if (histType == eMult) { - std::vector lCentBinsX = lQABins[1]; // MultBins + const auto& lCentBinsX = lQABins[1]; // MultBins nBinsCentX = static_cast(lCentBinsX[0]); minCentX = lCentBinsX[1]; maxCentX = lCentBinsX[2]; - std::vector lCentBinsY = lQABins[2]; // nContribBins + const auto& lCentBinsY = lQABins[2]; // nContribBins nBinsCentY = static_cast(lCentBinsY[0]); minCentY = lCentBinsY[1]; maxCentY = lCentBinsY[2]; } else { - std::vector lCentBins = lQABins[histType]; + const auto& lCentBins = lQABins[histType]; nBinsCent = static_cast(lCentBins[0]); minCent = lCentBins[1]; maxCent = lCentBins[2]; } for (int ba = 0; ba < eCutBeforeAfter_N; ba++) { - std::string name = Form("fHist%s[%s cut]", eventHistNames[histType], cutBeforeAfterNames[ba]); + std::string name = Form("fHist%s[%s cut]", EventHistNames[histType], CutBeforeAfterNames[ba]); std::string namefull; if constexpr (histType == eCent) { - std::string namefull = Form("Quality assurance of %s %s", tc.fCentEstm.c_str(), eventHistNames[histType]); + namefull = Form("Quality assurance of %s %s", tc.fCentEstm.c_str(), EventHistNames[histType]); } else if constexpr (histType == eMult) { - std::string namefull = Form("Quality assurance of %s %s vs. NContributors", tc.fMultEstm.c_str(), eventHistNames[histType]); + namefull = Form("Quality assurance of %s %s vs. NContributors", tc.fMultEstm.c_str(), EventHistNames[histType]); } else { - std::string namefull = Form("Quality assurance of %s", eventHistNames[histType]); + namefull = Form("Quality assurance of %s", EventHistNames[histType]); } if constexpr (histType == eMult) { @@ -1152,25 +1319,25 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam qa.fQAHistograms[histType][ba]->GetXaxis()->SetTitle("Reference multiplicity"); } else { qa.fQAHistograms[histType][ba] = new TH2F(name.c_str(), namefull.c_str(), nBinsCent, minCent, maxCent, nBinsCent, minCent, maxCent); - qa.fQAHistograms[histType][ba]->GetYaxis()->SetTitle(Form("Simulated %s", eventHistNames[histType])); - qa.fQAHistograms[histType][ba]->GetXaxis()->SetTitle(Form("Reconstructed %s", eventHistNames[histType])); + qa.fQAHistograms[histType][ba]->GetYaxis()->SetTitle(Form("Simulated %s", EventHistNames[histType])); + qa.fQAHistograms[histType][ba]->GetXaxis()->SetTitle(Form("Reconstructed %s", EventHistNames[histType])); } qa.fQAHistogramsList->Add(qa.fQAHistograms[histType][ba]); } } template - void bookCorrHistograms(T1 const& lCrBins, CorrHistograms& cr) + void bookCorrHistograms(T1 const& lCrBins/* , CorrHistograms& cr */) { - std::vector lCentBins = lCrBins[histType]; // define local array and initialize it from an array set in the configurables + const auto& lCentBins = lCrBins[histType]; // define local array and initialize it from an array set in the configurables int nBinsCent = static_cast(lCentBins[0]); float minCent = lCentBins[1]; float maxCent = lCentBins[2]; for (int ba = 0; ba < eCutBeforeAfter_N; ba++) { - std::string name = ""; - std::string namefull = Form("%s correlation %s cut", corrHistNames[histType], cutBeforeAfterNames[ba]); + std::string name; + std::string namefull = Form("%s correlation %s cut", CorrHistNames[histType], CutBeforeAfterNames[ba]); int nEstm = 0; if constexpr (histType == eCorrCent) { @@ -1180,24 +1347,24 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } for (int i = 0; i < nEstm; i++) { - std::string titleX = ""; + std::string titleX; if (histType == eCorrCent) { - titleX = Form("%s %s", centEstmNames[i], corrHistNames[histType]); + titleX = Form("%s %s", CentEstmNames[i], CorrHistNames[histType]); } else if (histType == eCorrMult) { - titleX = Form("%s %s", multEstmNames[i], corrHistNames[histType]); + titleX = Form("%s %s", MultEstmNames[i], CorrHistNames[histType]); // x axis bin... } for (int j = i + 1; j < nEstm; j++) { - std::string titleY = ""; + std::string titleY; if (histType == eCorrCent) { - name = Form("fHist%s[%s vs. %s][%s cut]", corrHistNames[histType], centEstmNames[i], centEstmNames[j], cutBeforeAfterNames[ba]); - titleY = Form("%s %s", centEstmNames[j], corrHistNames[histType]); + name = Form("fHist%s[%s vs. %s][%s cut]", CorrHistNames[histType], CentEstmNames[i], CentEstmNames[j], CutBeforeAfterNames[ba]); + titleY = Form("%s %s", CentEstmNames[j], CorrHistNames[histType]); cr.fCorrHistograms[histType][i][j][ba] = new TH2F(name.c_str(), namefull.c_str(), nBinsCent, minCent, maxCent, nBinsCent, minCent, maxCent); } else if (histType == eCorrMult) { - name = Form("fHist%s[%s vs. %s][%s cut]", corrHistNames[histType], multEstmNames[i], multEstmNames[j], cutBeforeAfterNames[ba]); - titleY = Form("%s %s", multEstmNames[j], corrHistNames[histType]); + name = Form("fHist%s[%s vs. %s][%s cut]", CorrHistNames[histType], MultEstmNames[i], MultEstmNames[j], CutBeforeAfterNames[ba]); + titleY = Form("%s %s", MultEstmNames[j], CorrHistNames[histType]); // y axis bins... cr.fCorrHistograms[histType][i][j][ba] = new TH2F(name.c_str(), namefull.c_str(), nBinsCent, minCent, maxCent, nBinsCent, minCent, maxCent); } @@ -1266,7 +1433,6 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam tc.fPrintSwitch = cfPrintSwitch; tc.fFileWithWeights = cfFileWithWeights; - tc.fRunNumber = cfRunNumber; qa.fQASwitch = cfQASwitch; wt.fWeightSwitch = cfWeightSwitch; @@ -1279,32 +1445,32 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // *) Book and nest all other TLists: pc.fParticleHistogramsList = new TList(); pc.fParticleHistogramsList->SetName("ParticleHistograms"); - pc.fParticleHistogramsList->SetOwner(kFALSE); + pc.fParticleHistogramsList->SetOwner(kTRUE); fBaseList->Add(pc.fParticleHistogramsList); // any nested TList in the base TList appears as a subdir in the output ROOT file ev.fEventHistogramsList = new TList(); ev.fEventHistogramsList->SetName("EventHistograms"); - ev.fEventHistogramsList->SetOwner(kFALSE); + ev.fEventHistogramsList->SetOwner(kTRUE); fBaseList->Add(ev.fEventHistogramsList); qa.fQAHistogramsList = new TList(); qa.fQAHistogramsList->SetName("QualityAssuranceHistograms"); - qa.fQAHistogramsList->SetOwner(kFALSE); + qa.fQAHistogramsList->SetOwner(kTRUE); fBaseList->Add(qa.fQAHistogramsList); wt.fWeightHistogramsList = new TList(); wt.fWeightHistogramsList->SetName("WeightHistograms"); - wt.fWeightHistogramsList->SetOwner(kFALSE); + wt.fWeightHistogramsList->SetOwner(kTRUE); fBaseList->Add(wt.fWeightHistogramsList); cr.fCorrHistogramsList = new TList(); cr.fCorrHistogramsList->SetName("CorrelationHistograms"); - cr.fCorrHistogramsList->SetOwner(kFALSE); + cr.fCorrHistogramsList->SetOwner(kTRUE); fBaseList->Add(cr.fCorrHistogramsList); mc.fMultiparticleCorrelationProfilesList = new TList(); mc.fMultiparticleCorrelationProfilesList->SetName("MultiparticleCorrelationProfiles"); - mc.fMultiparticleCorrelationProfilesList->SetOwner(kFALSE); + mc.fMultiparticleCorrelationProfilesList->SetOwner(kTRUE); fBaseList->Add(mc.fMultiparticleCorrelationProfilesList); std::vector> lPcBins = {tc.fPtBins, tc.fPhiBins}; @@ -1312,39 +1478,27 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam std::vector> lQABins = {tc.fCentBins, tc.fMultBins, tc.fNumContribBins}; std::vector> lCrBins = {tc.fCentBins, tc.fMultBins}; - bookParticleHistograms(lPcBins, pc); - bookParticleHistograms(lPcBins, pc); - bookEventHistograms(lEvBins, ev); - bookEventHistograms(lEvBins, ev); - bookEventHistograms(lEvBins, ev); - bookEventHistograms(lEvBins, ev); - bookEventHistograms(lEvBins, ev); - bookEventHistograms(lEvBins, ev); - bookQAHistograms(lQABins, qa); - bookQAHistograms(lQABins, qa); - bookCorrHistograms(lCrBins, cr); // if switch on ... - bookCorrHistograms(lCrBins, cr); - - if (wt.fWeightSwitch) { - wt.fWeightHistograms = getHistogramsWithWeights(tc.fFileWithWeights.c_str(), tc.fRunNumber.c_str()); - for (auto* const& hist : wt.fWeightHistograms) { // o2-linter: disable=const-ref-in-for-loop - wt.fWeightHistogramsList->Add(hist); - } + bookParticleHistograms(lPcBins/* , pc */); + bookParticleHistograms(lPcBins/* , pc */); + bookEventHistograms(lEvBins/* , ev */); + bookEventHistograms(lEvBins/* , ev */); + bookEventHistograms(lEvBins/* , ev */); + bookEventHistograms(lEvBins/* , ev */); + bookEventHistograms(lEvBins/* , ev */); + bookEventHistograms(lEvBins/* , ev */); + bookQAHistograms(lQABins/* , qa */); + bookQAHistograms(lQABins/* , qa */); + bookCorrHistograms(lCrBins/* , cr */); // if switch on ... + bookCorrHistograms(lCrBins/* , cr */); + + wt.fDummyPhiWeightHistogram = new TH1F("fDummyPhiWeightHistogram", "Dummy phi weight histogram", static_cast(tc.fPhiBins[0]), tc.fPhiBins[1], tc.fPhiBins[2]); + for (int i = 1; i <= wt.fDummyPhiWeightHistogram->GetNbinsX(); i++) { + wt.fDummyPhiWeightHistogram->SetBinContent(i, 1.); + } + wt.fDummyPtWeightHistogram = new TH1F("fDummyPtWeightHistogram", "Dummy pt weight histogram", static_cast(tc.fPtBins[0]), tc.fPtBins[1], tc.fPtBins[2]); + for (int i = 1; i <= wt.fDummyPtWeightHistogram->GetNbinsX(); i++) { + wt.fDummyPtWeightHistogram->SetBinContent(i, 1.); } - - mc.fTwoParticleCorrelationProfiles[eBefore] = new TProfile("prof2Before", "2-p correlation before cut", 3, 2., 5.); - mc.fTwoParticleCorrelationProfiles[eAfter] = new TProfile("prof2After", "2-p correlation after cut", 3, 2., 5.); - mc.fTwoParticleCorrelationProfiles[eBefore]->Sumw2(); - mc.fTwoParticleCorrelationProfiles[eAfter]->Sumw2(); - mc.fMultiparticleCorrelationProfilesList->Add(mc.fTwoParticleCorrelationProfiles[eBefore]); - mc.fMultiparticleCorrelationProfilesList->Add(mc.fTwoParticleCorrelationProfiles[eAfter]); - - mc.fFourParticleCorrelationProfiles[eBefore] = new TProfile("prof4Before", "4-p correlation before cut", 2, 3., 5.); - mc.fFourParticleCorrelationProfiles[eAfter] = new TProfile("prof4After", "4-p correlation after cut", 2, 3., 5.); - mc.fFourParticleCorrelationProfiles[eBefore]->Sumw2(); - mc.fFourParticleCorrelationProfiles[eAfter]->Sumw2(); - mc.fMultiparticleCorrelationProfilesList->Add(mc.fFourParticleCorrelationProfiles[eBefore]); - mc.fMultiparticleCorrelationProfilesList->Add(mc.fFourParticleCorrelationProfiles[eAfter]); } // end of void init(InitContext&) {