diff --git a/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx b/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx index fba1ce3cf58..eeceadca07e 100644 --- a/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx +++ b/PWGCF/MultiparticleCorrelations/Tasks/multiparticleCumulants.cxx @@ -13,6 +13,7 @@ /// \brief Task for producing multiparticle cumulants /// \author Pei-Ying Kuan, TU München, pei-ying.kuan@cern.ch +#include "Common/CCDB/EventSelectionParams.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" @@ -31,6 +32,7 @@ #include #include +#include #include #include #include @@ -129,8 +131,7 @@ enum EnCorrHistograms { static constexpr std::array CorrHistNames = { "Centrality", - "Multiplicity", -}; + "Multiplicity"}; enum EnCentEstm { eCentFT0C, @@ -156,22 +157,37 @@ static constexpr std::array MultEstmNames = { "FT0M", "FV0A"}; -enum EnCutBeforeAfter { +enum EnEstmCorr { + eFT0CFT0M, + eFT0CFV0A, + eFT0MFV0A, + eEstmCorr_N +}; + +static constexpr std::array EstmCorrName = { + "FT0CFT0M", + "FT0CFV0A", + "FT0MFV0A"}; + +enum EnBeforeAfter { eBefore, eAfter, - eCutBeforeAfter_N + eBeforeAfter_N }; -static constexpr std::array CutBeforeAfterNames = { +static constexpr std::array BeforeAfterNames = { "before", "after", -}; + "Before", + "After"}; + +static constexpr int NumHarmonics = 3; // *) Main task: struct MultiparticleCumulants { // this name is used in lower-case format to name the TDirectoryFile in AnalysisResults.root // *) Base TList to hold all output objects: - TString sBaseListName = "Default list name"; // yes, I declare it separately, because I need it also later in BailOut() function + TString sBaseListName = "Default list name"; OutputObj fBaseList{sBaseListName.Data(), OutputObjHandlingPolicy::AnalysisObject, OutputObjSourceType::OutputObjSource}; // *) Service: @@ -180,14 +196,15 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // *) Define configurables: Configurable cfDryRun{"cfDryRun", false, "book all histos and run without filling and calculating anything"}; - Configurable cfCentEstm{"cfCentEstm", "FT0M", "centrality estimator: TBD"}; - Configurable cfMultEstm{"cfMultEstm", "FV0A", "multiplicity estimator: TBD"}; + Configurable cfCentEstm{"cfCentEstm", "FT0M", "centrality estimator: FT0C, FT0M, FV0A"}; + Configurable cfMultEstm{"cfMultEstm", "FV0A", "multiplicity estimator: FT0C, FT0M, FV0A"}; Configurable cfQASwitch{"cfQASwitch", true, "quality assurance switch"}; Configurable cfWeightSwitch{"cfWeightSwitch", true, "weight switch"}; Configurable cfPrintSwitch{"cfPrintSwitch", true, "printing result switch"}; // *) Event cut switches + Configurable cfGeneralCutSwitch{"cfGeneralCutSwitch", true, "NoCollInTimeRangeStandard, NoCollInRofStandard, NoSameBunchPileup, IsVertexITSTPC, IsGoodITSLayersAll, IsGoodZvtxFT0vsPV, NoHighMultCollInPrevRof"}; Configurable cfVertexZCutSwitch{"cfVertexZCutSwitch", true, "vertex z cut switch"}; Configurable cfSel8CutSwitch{"cfSel8CutSwitch", true, "Sel8 cut switch"}; Configurable cfCentCutSwitch{"cfCentCutSwitch", true, "centrality cut switch"}; @@ -202,13 +219,16 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam Configurable cfTpcNClsFoundCutSwitch{"cfTpcNClsFoundCutSwitch", true, ""}; Configurable cfDCAXYCutSwitch{"cfDCAXYCutSwitch", true, ""}; Configurable cfDCAZCutSwitch{"cfDCAZCutSwitch", true, ""}; + Configurable cfEtaGapSwitch{"cfEtaGapSwitch", true, "Eta gap switch"}; // *) Event cut Configurable> cfVertexZCut{"cfVertexZCut", {-10., 10.}, "vertex z position range: {min, max}[cm]"}; Configurable> cfCentCut{"cfCentCut", {10., 20.}, "centrality range: {min, max}[%]"}; Configurable> cfNumContribCut{"cfNumContribCut", {0, 3000.}, "NContribution range: {min, max}"}; - Configurable cfCentCorrCut{"cfCentCorrCut", 5, "centrality difference maximum"}; - Configurable cfMultCorrCut{"cfMultCorrCut", 0.2, "multiplicity difference maximum"}; + Configurable> cfCentCorrCut{"cfCentCorrCut", {1., 10.}, "parameters of limits in centralities 2D histograms, (x-t)/m < y < mx+t: {m, t}"}; + Configurable> cfMultABCorrCut{"cfMultABCorrCut", {1.4, 300., 1.4, 300.}, "parameters of limits in FT0C vs. FT0M multiplicities 2D histograms, ax+b < y < cx+d: {a,b,c,d}"}; + Configurable> cfMultACCorrCut{"cfMultACCorrCut", {1.4, 300., 1.4, 300.}, "parameters of limits in FT0C vs. FV0A multiplicities 2D histograms, ax+b < y < cx+d: {a,b,c,d}"}; + Configurable> cfMultBCCorrCut{"cfMultBCCorrCut", {1.4, 300., 1.4, 300.}, "parameters of limits in FT0M vs. FV0A multiplicities 2D histograms, ax+b < y < cx+d: {a,b,c,d}"}; // *) Particle cut Configurable> cfPtCut{"cfPtCut", {0.2, 5.0}, "Pt range: {min, max}[GeV], with convention: min <= Pt < max"}; @@ -217,6 +237,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam Configurable> cfTpcNClsFoundCut{"cfTpcNClsFoundCut", {70., 160.}, "range of found TPC clusters for this track geometry: {min, max}"}; Configurable> cfDCAXYCut{"cfDCAXYCut", {-3.2, 3.2}, "range of distance-of-closest-approach (DCA) of the extrapolated track to the primary position in XY-direction: {min, max}[cm]"}; Configurable> cfDCAZCut{"cfDCAZCut", {-2.4, 2.4}, "range of distance-of-closest-approach (DCA) of the extrapolated track to the primary position in Z-direction: {min, max}[cm]"}; + Configurable cfEtaGap{"cfEtaGap", 1., "|dEta| > gap"}; // *) Others Configurable cfFileWithWeights{"cfFileWithWeights", "/scratch3/go52dab/O2tutorial/tutorial3-6/weights.root", "path to external ROOT file which holds all particle weights in O2 format"}; @@ -225,11 +246,15 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam Configurable> cfPtBins{"cfPtBins", {1000, 0., 100.}, "nPtBins, ptMin, ptMax"}; Configurable> cfPhiBins{"cfPhiBins", {1000, 0., o2::constants::math::TwoPI}, "nPhiBins, phiMin, phiMax"}; Configurable> cfCentBins{"cfCentBins", {100, 0., 100.}, "nCenBins, cenMin, cenMax"}; - Configurable> cfMultBins{"cfMultBins", {100, 0., 5000.}, "nMultBins, MultMin, MultMax"}; + Configurable> cfFT0CMultBins{"cfFT0CMultBins", {100, 0., 60000.}, "nMultBins, MultMin, MultMax"}; + Configurable> cfFT0MMultBins{"cfFT0MMultBins", {100, 0., 250000.}, "nMultBins, MultMin, MultMax"}; + Configurable> cfFV0AMultBins{"cfFV0AMultBins", {100, 0., 200000.}, "nMultBins, MultMin, MultMax"}; Configurable> cfVerXBins{"cfVerXBins", {100, -0.05, 0.05}, "nVerXBins, VerXMin, VerXMax"}; Configurable> cfVerYBins{"cfVerYBins", {100, -0.05, 0.05}, "nVerYBins, VerYMin, VerYMax"}; Configurable> cfVerZBins{"cfVerZBins", {100, -50., 50.}, "nVerZBins, VerZMin, VerZMax"}; Configurable> cfNumContribBins{"cfNumContribBins", {100, 0., 5000.}, "nNumContribBins, NumContribMin, NumContribMax"}; + Configurable> cfTwoParticleCorrBins{"cfTwoParticleCorrBins", {100, 0., 2e-3}, "nTwoParticleCorrBins, TwoParticleCorrMin, TwoParticleCorrMax"}; + Configurable> cfFourParticleCorrBins{"cfFourParticleCorrBins", {100, 0., 1e-7}, "nFourParticleCorrBins, FourParticleCorrMin, FourParticleCorrMax"}; // *) Define and initialize all data members to be called in the main process* functions: // **) Task configuration: @@ -242,6 +267,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam bool fPrintSwitch = true; + bool fGeneralCutSwitch = true; bool fVertexZCutSwitch = true; bool fSel8CutSwitch = true; bool fCentCutSwitch = true; @@ -255,29 +281,36 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam bool fTpcNClsFoundCutSwitch = true; bool fDCAXYCutSwitch = true; bool fDCAZCutSwitch = true; + bool fEtaGapSwitch = true; std::vector fVertexZCut = {-10., 10.}; std::vector fCentCut = {10., 20.}; std::vector fNumContribCut = {0, 3000.}; - float fCentCorrCut = 5; - float fMultCorrCut = 0.2; - + std::vector fCentCorrCut = {1., 10.}; + std::vector>> fMultAllCorrCut = {{{0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}}, {{0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}}, {{0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}}}; std::vector fPtCut = {0.2, 5.0}; std::vector fEtaCut = {-0.8, 0.8}; std::vector fSignCut = {1, 0, 1}; std::vector fTpcNClsFoundCut = {70., 160.}; std::vector fDCAXYCut = {-3.2, 3.2}; std::vector fDCAZCut = {-2.4, 2.4}; + float fEtaGap = 1.; - std::vector fPtBins = {0., 100.}; - std::vector fPhiBins = {0., o2::constants::math::TwoPI}; + std::vector fPtBins = {1000, 0., 100.}; + std::vector fPhiBins = {1000, 0., o2::constants::math::TwoPI}; - std::vector fCentBins = {0., 100.}; - std::vector fMultBins = {0, 5000}; - std::vector fVerXBins = {-0.05, 0.05}; - std::vector fVerYBins = {-0.05, 0.05}; - std::vector fVerZBins = {-50., 50.}; - std::vector fNumContribBins = {0, 5000}; + std::vector fCentBins = {100, 0., 100.}; + std::vector fMultBins = {100, 0, 5000}; + std::vector fFT0CMultBins = {100, 0., 60000.}; + std::vector fFT0MMultBins = {100, 0., 250000.}; + std::vector fFV0AMultBins = {100, 0., 200000.}; + + std::vector fVerXBins = {100, -0.05, 0.05}; + std::vector fVerYBins = {100, -0.05, 0.05}; + std::vector fVerZBins = {100, -50., 50.}; + std::vector fNumContribBins = {100, 0, 5000}; + std::vector fTwoParticleCorrBins = {100, 0, 2e-3}; + std::vector fFourParticleCorrBins = {100, 0, 1e-7}; std::string fFileWithWeights = "/scratch3/go52dab/O2tutorial/analysis_code/weights.root"; @@ -285,12 +318,12 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam struct ParticleHistograms { TList* fParticleHistogramsList = nullptr; - std::array, 2>, eParticleHistograms_N> fParticleHistograms{}; + std::array, eBeforeAfter_N>, eParticleHistograms_N> fParticleHistograms{}; // [type][before/after cut][rec/sim] } pc; struct EventHistograms { TList* fEventHistogramsList = nullptr; - std::array, 2>, eEventHistograms_N> fEventHistograms{}; + std::array, eBeforeAfter_N>, eEventHistograms_N> fEventHistograms{}; // [type][before/after cut][rec/sim] } ev; struct QAHistograms { @@ -301,7 +334,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam struct CorrHistograms { TList* fCorrHistogramsList = nullptr; - std::array, eMultEstm_N>, eMultEstm_N>, eCorrHistograms_N> fCorrHistograms{}; + std::array, eMultEstm_N>, eMultEstm_N>, eCorrHistograms_N> fCorrHistograms{}; // [mult/cent][type][type][before/after cut] + std::array, eCorrHistograms_N + 2> fCorrBounds{}; //[cent/multAB/multAC/multBC][upper/lower] } cr; struct WeightHistograms { @@ -320,24 +354,6 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam TH1F* fDummyPtWeightHistogram = nullptr; } wt; - struct EventByEventQuantities { - int fRunNumber = 0; - float fReferenceMultiplicity = 0.; - float fCentrality = 0.; - float fCentralitySim = 0.; - float fImpactParameter = 0.; - float fNumContrib = 0.; - 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 = nullptr; - std::map fMultiparticleCorrelationByRunMap; - std::array fTwoParticleCorrelationProfiles = {{nullptr}}; - std::array fFourParticleCorrelationProfiles = {{nullptr}}; - } mc; - struct MultiparticleCorrelationCalculation { int h1 = 0; int h2 = 0; @@ -349,18 +365,47 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam 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 MaxHarmonic = 17; // need 4 + 2 + 2 + 4 + 1 at least for SC(4,2) static constexpr int MaxPower = MaxCorrelator + 1; - static constexpr int NumSC = 2; // need SC(3,2) and SC(4,2) std::array, MaxHarmonic> fQvectorBefore; - std::array, MaxHarmonic> fQvectorAfter; // All needed Q-vector components + std::array, MaxHarmonic> fQvectorAfter; + std::array, MaxHarmonic> fQvectorAfterA; // Q-vector with eta gap + std::array, MaxHarmonic> fQvectorAfterB; // Q-vector with eta gap } mcc; - template + struct MultiparticleCorrelationProfile { + TList* fMultiparticleCorrelationProfilesList = nullptr; + TList* fMultiparticleCorrelationBackupList = nullptr; + std::map fMultiparticleCorrelationByRunMap; + std::array fTwoParticleCorrelationProfiles{}; // [cut] + std::array fFourParticleCorrelationProfiles{}; + std::array, eBeforeAfter_N>, eBeforeAfter_N> fTwoParticleCorrelationHistograms{}; // [cut][event weight][n] + std::array, eBeforeAfter_N>, eBeforeAfter_N> fFourParticleCorrelationHistograms{}; + std::array fTwoParticleCorrelationGapHistograms{}; + } mc; + + struct EventByEventQuantities { + int fRunNumber = 0; + float fReferenceMultiplicity = 0.; + std::array fMultiplicityMaxEbye{}; + float fCentrality = 0.; + float fCentralitySim = 0.; + float fImpactParameter = 0.; + float fNumContrib = 0.; + std::array, eBeforeAfter_N> fTwoParticleCorrelationEbye{}; //<2> [before, after][v2^2, v3^2, v4^2] + std::array, eBeforeAfter_N> fTwoParticleCorrelationMinEbye{}; + std::array, eBeforeAfter_N> fTwoParticleCorrelationMaxEbye{}; + std::array, eBeforeAfter_N> fFourParticleCorrelationEbye{}; //<4> [before, after][v2^4, v3^2 v2^2, v4^2 v2^2] + std::array, eBeforeAfter_N> fFourParticleCorrelationMinEbye{}; + std::array, eBeforeAfter_N> fFourParticleCorrelationMaxEbye{}; + } ebye; + + template // rm should be rec and sim bool ctEventCuts(T1 const& collision, T2 const& rlCollisionCentAll, T3 const& rlCollisionMultAll) { bool pass = true; + bool bGeneralCut = true; bool bVertexZCut = true; bool bSel8Cut = true; bool bCentCut = true; @@ -368,12 +413,20 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam bool bCentCorrCut = true; bool bMultCorrCut = true; - // *) For real event and MC event + // *) For rec event and sim event bVertexZCut = collision.posZ() < tc.fVertexZCut[1] && collision.posZ() > tc.fVertexZCut[0]; - // *) For real event only - if constexpr (rm == eRl) { + // *) For rec event only + if constexpr (rs == eRec) { + // *) General cuts + bGeneralCut = collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard) && + collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard) && + collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup) && + collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC) && + collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll) && + collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV) && + collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof); // *) Sel8Cut bSel8Cut = collision.sel8(); // *) CentCut @@ -389,29 +442,38 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam iCent = rlCollisionCentAll[i]; for (int j = i + 1; j < eCentEstm_N; j++) { jCent = rlCollisionCentAll[j]; - bCentCorrCut = std::abs((iCent - jCent)) < tc.fCentCorrCut; + float upper = cr.fCorrBounds[eCorrCent][0]->Eval(iCent); + float lower = cr.fCorrBounds[eCorrCent][1]->Eval(iCent); + bCentCorrCut &= jCent >= lower && jCent <= upper; } } // *) MultCorrCut float iMult = 0.; float jMult = 0.; + int nComb = 0; for (int i = 0; i < eMultEstm_N; i++) { iMult = rlCollisionMultAll[i]; for (int j = i + 1; j < eMultEstm_N; j++) { jMult = rlCollisionMultAll[j]; - bMultCorrCut = std::abs((iMult - jMult) / iMult) < tc.fMultCorrCut; + float upper = cr.fCorrBounds[eCorrMult + nComb][0]->Eval(iMult); + float lower = cr.fCorrBounds[eCorrMult + nComb][1]->Eval(iMult); + bMultCorrCut &= jMult >= lower && jMult <= upper; + nComb += 1; } } } - // *) For MC event only - if constexpr (rm == eMc) { + // *) For sim event only + if constexpr (rs == eSim) { // *) CentCut bCentCut = ebye.fCentralitySim < tc.fCentCut[1] && ebye.fCentralitySim > tc.fCentCut[0]; } // *) Combine all switches + if (tc.fGeneralCutSwitch) { + pass &= bGeneralCut; + } if (tc.fVertexZCutSwitch) { pass &= bVertexZCut; } @@ -434,7 +496,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam return pass; } - template + template bool ctParticleCuts(T1 const& track) { bool pass = true; @@ -446,12 +508,12 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam bool bDCAXYCut = true; bool bDCAZCut = true; - // *) For real event and MC event + // *) For rec event and sim event bPtCut = track.pt() < tc.fPtCut[1] && track.pt() > tc.fPtCut[0]; bEtaCut = track.eta() < tc.fEtaCut[1] && track.eta() > tc.fEtaCut[0]; - // *) For real event only - if constexpr (rm == eRl) { + // *) For rec event only + if constexpr (rs == eRec) { bSignCut = (track.sign() == -1 && tc.fSignCut[0]) || (track.sign() == 0 && tc.fSignCut[1]) || (track.sign() == 1 && tc.fSignCut[2]); @@ -463,8 +525,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam track.dcaZ() > tc.fDCAZCut[0]; } - // *) For mc event only - if constexpr (rm == eMc) { + // *) For sim event only + if constexpr (rs == eSim) { TParticlePDG* particle = pdg->GetParticle(track.pdgCode()); if (!particle) { // LOGF(warning, "PDG code %d not found", track.pdgCode()); @@ -500,33 +562,53 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam return pass; } - TComplex mccQ(int n, int p, EnCutBeforeAfter eba) + TComplex mccQ(int n, int p, EnBeforeAfter eba, int egap) { // Using the fact that Q{-n,p} = Q{n,p}^*. - if (eba == eBefore) { + // Q-vector before cut: if (n >= 0) { return mcc.fQvectorBefore[n][p]; } return TComplex::Conjugate(mcc.fQvectorBefore[-n][p]); } + if (egap == 0) { + // Q-vector after cut: + if (n >= 0) { + return mcc.fQvectorAfter[n][p]; + } + return TComplex::Conjugate(mcc.fQvectorAfter[-n][p]); + } + if (egap == 1) { + // Q-vector after cut, eta < gap: + if (n >= 0) { + return mcc.fQvectorAfterA[n][p]; + } + return TComplex::Conjugate(mcc.fQvectorAfterA[-n][p]); + } + // Q-vector after cut, eta > gap: if (n >= 0) { - return mcc.fQvectorAfter[n][p]; + return mcc.fQvectorAfterB[n][p]; } - return TComplex::Conjugate(mcc.fQvectorAfter[-n][p]); + return TComplex::Conjugate(mcc.fQvectorAfterB[-n][p]); + } + + TComplex mccTwo(int n1, int n2, EnBeforeAfter eba, int egap) + { + return mccQ(n1, 1, eba, egap) * mccQ(n2, 1, eba, egap) - mccQ(n1 + n2, 2, eba, egap); } template - TComplex mccRecursion(int n, std::array harmonic, EnCutBeforeAfter eba, int mult = 1, int skip = 0) + TComplex mccRecursion(int n, std::array harmonic, EnBeforeAfter eba, int egap, 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)); + TComplex c(mccQ(harmonic[nm1], mult, eba, egap)); if (nm1 == 0) { return c; } - c *= mccRecursion(nm1, harmonic, eba); + c *= mccRecursion(nm1, harmonic, eba, egap); if (nm1 == skip) { return c; } @@ -537,7 +619,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam int hhold = harmonic[counter1]; harmonic[counter1] = harmonic[nm2]; harmonic[nm2] = hhold + harmonic[nm1]; - TComplex c2(mccRecursion(nm1, harmonic, eba, multp1, nm2)); + TComplex c2(mccRecursion(nm1, harmonic, eba, egap, multp1, nm2)); int counter2 = n - 3; while (counter2 >= skip) { harmonic[nm2] = harmonic[counter1]; @@ -546,7 +628,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam hhold = harmonic[counter1]; harmonic[counter1] = harmonic[nm2]; harmonic[nm2] = hhold + harmonic[nm1]; - c2 += mccRecursion(nm1, harmonic, eba, multp1, counter2); + c2 += mccRecursion(nm1, harmonic, eba, egap, multp1, counter2); --counter2; } harmonic[nm2] = harmonic[counter1]; @@ -563,8 +645,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // Get TObject pointer from TList, even if it's in some nested TList. // Foreseen to be used to fetch histograms or profiles from files directly. // Some ideas taken from TCollection::ls() - // If you have added histograms directly to files (without TList's), then - // you can fetch them directly with file->Get("hist-name"). + // If you have added histograms directly to files (without TList's), then you can fetch them directly with file->Get("hist-name"). // Usage: TH1D* hist = (TH1D*) // getObjectFromList("some-valid-TList-pointer","some-object-name"); @@ -676,9 +757,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } } - // OK, we got the desired TList with efficiency corrections, after that we - // can use the common code for all 3 cases (local, AliEn, CCDB, that - // common code is below) + // OK, we got the desired TList with efficiency corrections, after that we can use the common code for all 3 cases (local, AliEn, CCDB, that common code is below) } else { // this is the local case, please handle this one now: // Check if the external ROOT file exists at specified path: @@ -756,6 +835,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam for (int p = 0; p < mcc.MaxPower; p++) { mcc.fQvectorBefore[h][p] = TComplex(0., 0.); mcc.fQvectorAfter[h][p] = TComplex(0., 0.); + mcc.fQvectorAfterA[h][p] = TComplex(0., 0.); + mcc.fQvectorAfterB[h][p] = TComplex(0., 0.); } } @@ -765,9 +846,6 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // 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); @@ -776,9 +854,6 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // 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); @@ -787,9 +862,6 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // 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); @@ -845,24 +917,66 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam 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]); + mc.fMultiparticleCorrelationBackupList = nullptr; + mc.fMultiparticleCorrelationBackupList = new TList(); + mc.fMultiparticleCorrelationBackupList->SetName("backup"); + mc.fMultiparticleCorrelationBackupList->SetOwner(kTRUE); + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fMultiparticleCorrelationBackupList); + + for (int i = 0; i < eBeforeAfter_N; i++) { // before/after cut + mc.fTwoParticleCorrelationProfiles[i] = nullptr; + mc.fFourParticleCorrelationProfiles[i] = nullptr; + + for (int j = 0; j < eBeforeAfter_N; j++) { // before/after weight + mc.fTwoParticleCorrelationHistograms[i][j].fill(nullptr); + mc.fFourParticleCorrelationHistograms[i][j].fill(nullptr); + } + } + + for (int i = 0; i < eBeforeAfter_N; i++) { // before/after cut + mc.fTwoParticleCorrelationProfiles[i] = new TProfile(Form("prof2%sCut", BeforeAfterNames[i + 2]), Form("2-p correlation %s cut", BeforeAfterNames[i]), 3, 2., 5.); + mc.fTwoParticleCorrelationProfiles[i]->Sumw2(); + + for (int k = 0; k < NumHarmonics; k++) { // v2, v3, v4 + mc.fTwoParticleCorrelationProfiles[i]->GetXaxis()->SetBinLabel(k + 1, Form("v%d^2", k + 2)); + if (static_cast(i)) { + mc.fTwoParticleCorrelationGapHistograms[k] = new TH1D(Form("hist2v%dAfterCutWithGap", k + 2), Form("2-p correlation v%d^2 after cut with gap", k + 2), static_cast(tc.fTwoParticleCorrBins[0]), tc.fTwoParticleCorrBins[1], tc.fTwoParticleCorrBins[2]); + mc.fTwoParticleCorrelationGapHistograms[k]->Sumw2(); + mc.fTwoParticleCorrelationGapHistograms[k]->SetOption("HIST"); + } + for (int j = 0; j < eBeforeAfter_N; j++) { // before/after weight + mc.fTwoParticleCorrelationHistograms[i][j][k] = new TH1D(Form("hist2v%d%sCut%sWeight", k + 2, BeforeAfterNames[i + 2], BeforeAfterNames[j + 2]), Form("2-p correlation v%d^2 %s cut %s weight", k + 2, BeforeAfterNames[i], BeforeAfterNames[j]), static_cast(tc.fTwoParticleCorrBins[0]), tc.fTwoParticleCorrBins[1], tc.fTwoParticleCorrBins[2]); + if (static_cast(j)) { + mc.fTwoParticleCorrelationHistograms[i][j][k]->Sumw2(); + mc.fTwoParticleCorrelationHistograms[i][j][k]->SetOption("HIST"); + } + mc.fMultiparticleCorrelationBackupList->Add(mc.fTwoParticleCorrelationHistograms[i][j][k]); + } + if (static_cast(i)) { + mc.fMultiparticleCorrelationBackupList->Add(mc.fTwoParticleCorrelationGapHistograms[k]); + } + } + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fTwoParticleCorrelationProfiles[i]); + } + + for (int i = 0; i < eBeforeAfter_N; i++) { // before/after cut + mc.fFourParticleCorrelationProfiles[i] = new TProfile(Form("prof4%sCut", BeforeAfterNames[i + 2]), Form("4-p correlation %s cut", BeforeAfterNames[i]), NumHarmonics, 2., 5.); + mc.fFourParticleCorrelationProfiles[i]->Sumw2(); + + for (int k = 0; k < NumHarmonics; k++) { // v2v2, v2v3, v2v4 + mc.fFourParticleCorrelationProfiles[i]->GetXaxis()->SetBinLabel(k + 1, Form("v2^2v%d^2", k + 2)); + + for (int j = 0; j < eBeforeAfter_N; j++) { // before/after weight + mc.fFourParticleCorrelationHistograms[i][j][k] = new TH1D(Form("hist4v2v%d%sCut%sWeight", k + 2, BeforeAfterNames[i + 2], BeforeAfterNames[j + 2]), Form("4-p correlation v2^2v%d^2 %s cut %s weight", k + 2, BeforeAfterNames[i], BeforeAfterNames[j]), static_cast(tc.fFourParticleCorrBins[0]), tc.fFourParticleCorrBins[1], tc.fFourParticleCorrBins[2]); + if (static_cast(j)) { + mc.fFourParticleCorrelationHistograms[i][j][k]->Sumw2(); + mc.fFourParticleCorrelationHistograms[i][j][k]->SetOption("HIST"); + } + mc.fMultiparticleCorrelationBackupList->Add(mc.fFourParticleCorrelationHistograms[i][j][k]); + } + } + mc.fMultiparticleCorrelationByRunMap[ebye.fRunNumber]->Add(mc.fFourParticleCorrelationProfiles[i]); + } } // Real data centrality: @@ -897,6 +1011,12 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam if (tc.fMultEstm == MultEstmNames[i]) { rlCollisionMult = rlCollisionMultAll[i]; } + if (rlCollisionMultAll[i] >= ebye.fMultiplicityMaxEbye[i]) { + ebye.fMultiplicityMaxEbye[i] = rlCollisionMultAll[i]; + if (tc.fPrintSwitch) { + LOGF(info, "max %s multiplicity so far = %f", MultEstmNames[i], ebye.fMultiplicityMaxEbye[i]); + } + } } // Real data nContrib: @@ -935,33 +1055,46 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam 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]; - if (!h) { - LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrCent][%d][%d][eBefore]", i, j); + if (tc.fCentCorrCutSwitch) { + 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]; + if (!h) { + LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrCent][%d][%d][eBefore]", i, j); + } + if (rlCollisionCentAll[i] >= 0. && rlCollisionCentAll[j] >= 0.) { + h->Fill(rlCollisionCentAll[i], rlCollisionCentAll[j]); + } } - h->Fill(rlCollisionCentAll[i], rlCollisionCentAll[j]); } } // 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]; - if (!h) { - LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrMult][%d][%d][eBefore]", i, j); + 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][eBefore]; + if (!h) { + LOGF(fatal, "Missing histogram cr.fCorrHistograms[eCorrMult][%d][%d][eBefore]", i, j); + } + if (rlCollisionMultAll[i] >= 0. && rlCollisionMultAll[j] >= 0.) { + h->Fill(rlCollisionMultAll[i], rlCollisionMultAll[j]); + } } - h->Fill(rlCollisionMultAll[i], rlCollisionMultAll[j]); } } + // Fill NContrib vs. multiplicity QA histogram before cut: + if (qa.fQASwitch) { + qa.fQAHistograms[eQAMultNumContrib][eBefore]->Fill(rlCollisionMult, rlCollisionNumContrib); + } + // If RecAndSim: if constexpr (rs == eRecAndSim) { if (!collision.has_mcCollision()) { if (tc.fPrintSwitch) { - LOGF(warning, " No MC collision for this collision, skip..."); + LOGF(warning, "No MC collision for this collision, skip..."); } } else { // Define MC collision: @@ -969,9 +1102,11 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam // 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; + double mcCollisionCentDouble = 0.; + double b = mccollision.impactParameter() * 1e-15; // convert fm to m + double xs = 7.71e-28; // convert barn to m^2 + mcCollisionCentDouble = o2::constants::math::PI * b * b / xs * 100; + mcCollisionCent = static_cast(mcCollisionCentDouble); // Event-by-event quantity: ebye.fCentralitySim = mcCollisionCent; @@ -989,27 +1124,26 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam ev.fEventHistograms[eVertexZ][eSim][eBefore]->Fill(mccollision.posZ()); // Fill MC event histograms after cut: - if (ctEventCuts(mccollision, rlCollisionCentAll, rlCollisionMultAll)) { + 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 centrality QA histogram before/after cut: 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)) { + 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)) { + if (ctEventCuts(collision, rlCollisionCentAll, rlCollisionMultAll)) { ev.fEventHistograms[eCent][eRec][eAfter]->Fill(rlCollisionCent); ev.fEventHistograms[eMult][eRec][eAfter]->Fill(rlCollisionMult); ev.fEventHistograms[eVertexX][eRec][eAfter]->Fill(collision.posX()); @@ -1043,16 +1177,24 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } } - // Fail the event cut, skip this collision: + // Fill NContrib vs. multiplicity QA histogram after cut: + if (qa.fQASwitch) { + qa.fQAHistograms[eQAMultNumContrib][eAfter]->Fill(rlCollisionMult, rlCollisionNumContrib); + } + } else { + // Fail the event cut, skip this collision: return; } } int nTracksBefore = tracks.size(); int nTracksAfter = 0; + int nTracksAfterA = 0; + int nTracksAfterB = 0; // Calculate Q-vectors for available angles and weights: + double dEta = 0.; double dPhi = 0.; // particle angle double dPt = 0.; double wPhi = 1.; // particle weight @@ -1075,6 +1217,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam pc.fParticleHistograms[ePhi][eRec][eBefore]->Fill(track.phi()); // Calculating Q-vector before cut: + dEta = track.eta(); dPhi = track.phi(); dPt = track.pt(); if (wt.fWeightSwitch) { @@ -1089,6 +1232,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam auto* histPt = wt.fPtWeightHistogramsMap.at(ebye.fRunNumber); wPhi = histPhi->GetBinContent(histPhi->GetXaxis()->FindBin(dPhi)); wPt = histPt->GetBinContent(histPt->GetXaxis()->FindBin(dPt)); + wPhi *= wPt; } @@ -1101,7 +1245,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } } - if (ctParticleCuts(track)) { + if (ctParticleCuts(track)) { // Fill particle histograms after cut: pc.fParticleHistograms[ePt][eRec][eAfter]->Fill(track.pt()); @@ -1114,10 +1258,22 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam wPhiToPowerP = std::pow(wPhi, p); } mcc.fQvectorAfter[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); + if (tc.fEtaGapSwitch) { + if (dEta < -tc.fEtaGap / 2.) { + mcc.fQvectorAfterA[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); + } else if (dEta > tc.fEtaGap / 2.) { + mcc.fQvectorAfterB[h][p] += TComplex(wPhiToPowerP * std::cos(h * dPhi), wPhiToPowerP * std::sin(h * dPhi)); + } + } } } nTracksAfter += 1; + if (dEta < -tc.fEtaGap / 2) { + nTracksAfterA += 1; + } else if (dEta > tc.fEtaGap / 2) { + nTracksAfterB += 1; + } } // ... @@ -1142,7 +1298,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam pc.fParticleHistograms[ePt][eSim][eBefore]->Fill(mcparticle.pt()); pc.fParticleHistograms[ePhi][eSim][eBefore]->Fill(mcparticle.phi()); - if (ctParticleCuts(mcparticle)) { + 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()); @@ -1152,50 +1308,137 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } // 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++) { + for (int i = 0; i < NumHarmonics; i++) { + // 2-p correlations: mcc.h1 = -(i + 2); mcc.h2 = i + 2; 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(); + // Before cut: + TComplex twoRecursionBefore = mccRecursion(2, harmonicsTwoNum, eBefore, 0) / mccRecursion(2, harmonicsTwoDen, eBefore, 0).Re(); + double wTwoRecursionBefore = mccRecursion(2, harmonicsTwoDen, eBefore, 0).Re(); ebye.fTwoParticleCorrelationEbye[eBefore][i] = twoRecursionBefore.Re(); - TComplex twoRecursionAfter = mccRecursion(2, harmonicsTwoNum, eAfter) / mccRecursion(2, harmonicsTwoDen, eAfter).Re(); - double wTwoRecursionAfter = mccRecursion(2, harmonicsTwoDen, eAfter).Re(); - ebye.fTwoParticleCorrelationEbye[eAfter][i] = twoRecursionAfter.Re(); - - if (nTracksBefore > 1) { + if (nTracksBefore > 1 && wTwoRecursionBefore > 0.) { + // Fill vn^2 profile and histogram: mc.fTwoParticleCorrelationProfiles[eBefore]->Fill(i + 2.5, ebye.fTwoParticleCorrelationEbye[eBefore][i], wTwoRecursionBefore); + mc.fTwoParticleCorrelationHistograms[eBefore][eBefore][i]->Fill(ebye.fTwoParticleCorrelationEbye[eBefore][i]); + mc.fTwoParticleCorrelationHistograms[eBefore][eAfter][i]->Fill(ebye.fTwoParticleCorrelationEbye[eBefore][i], wTwoRecursionBefore); + + // Find max and min: + if (ebye.fTwoParticleCorrelationEbye[eBefore][i] > ebye.fTwoParticleCorrelationMaxEbye[eBefore][i]) { + ebye.fTwoParticleCorrelationMaxEbye[eBefore][i] = ebye.fTwoParticleCorrelationEbye[eBefore][i]; + LOGF(info, "max v%d^2 so far = %e", i + 2, ebye.fTwoParticleCorrelationEbye[eBefore][i]); + } + if (ebye.fTwoParticleCorrelationEbye[eBefore][i] < ebye.fTwoParticleCorrelationMinEbye[eBefore][i]) { + ebye.fTwoParticleCorrelationMinEbye[eBefore][i] = ebye.fTwoParticleCorrelationEbye[eBefore][i]; + LOGF(info, "min v%d^2 so far = %e", i + 2, ebye.fTwoParticleCorrelationEbye[eBefore][i]); + } + } else { + LOGF(warning, "cent=%f, nTracksBefore = %d, wTwoRecursionBefore = %e", rlCollisionCent, nTracksBefore, wTwoRecursionBefore); } - if (nTracksAfter > 1) { + + // After cut: + TComplex twoRecursionAfter = mccRecursion(2, harmonicsTwoNum, eAfter, 0) / mccRecursion(2, harmonicsTwoDen, eAfter, 0).Re(); + double wTwoRecursionAfter = mccRecursion(2, harmonicsTwoDen, eAfter, 0).Re(); + ebye.fTwoParticleCorrelationEbye[eAfter][i] = twoRecursionAfter.Re(); + + if (nTracksAfter > 1 && wTwoRecursionAfter > 0.) { + // Fill vn^2 profile and histogram: mc.fTwoParticleCorrelationProfiles[eAfter]->Fill(i + 2.5, ebye.fTwoParticleCorrelationEbye[eAfter][i], wTwoRecursionAfter); + mc.fTwoParticleCorrelationHistograms[eAfter][eBefore][i]->Fill(ebye.fTwoParticleCorrelationEbye[eAfter][i]); + mc.fTwoParticleCorrelationHistograms[eAfter][eAfter][i]->Fill(ebye.fTwoParticleCorrelationEbye[eAfter][i], wTwoRecursionAfter); + + // Find max and min: + if (ebye.fTwoParticleCorrelationEbye[eAfter][i] > ebye.fTwoParticleCorrelationMaxEbye[eAfter][i]) { + ebye.fTwoParticleCorrelationMaxEbye[eAfter][i] = ebye.fTwoParticleCorrelationEbye[eAfter][i]; + LOGF(info, "max v%d^2 so far = %e", i + 2, ebye.fTwoParticleCorrelationEbye[eAfter][i]); + } + if (ebye.fTwoParticleCorrelationEbye[eAfter][i] < ebye.fTwoParticleCorrelationMinEbye[eAfter][i]) { + ebye.fTwoParticleCorrelationMinEbye[eAfter][i] = ebye.fTwoParticleCorrelationEbye[eAfter][i]; + LOGF(info, "min v%d^2 so far = %e", i + 2, ebye.fTwoParticleCorrelationEbye[eAfter][i]); + } + } else { + LOGF(warning, "cent=%f, nTracksAfter = %d, wTwoRecursionAfter = %e", rlCollisionCent, nTracksAfter, wTwoRecursionAfter); } + + // After cut, with gap: + TComplex qva = mcc.fQvectorAfterA[mcc.h2][0]; + TComplex qvb = mcc.fQvectorAfterB[mcc.h2][0]; + double gap = 0.; + if (nTracksAfterA * nTracksAfterB > 0) { + gap = (qva * TComplex::Conjugate(qvb)).Re() / (nTracksAfterA * nTracksAfterB); + mc.fTwoParticleCorrelationGapHistograms[i]->Fill(gap, nTracksAfterA * nTracksAfterB); + } else { + LOGF(warning, "cent=%f, nTracksAfterA = %d, nTracksAfterB = %d", rlCollisionCent, nTracksAfterA, nTracksAfterB); + } + + // // Check if results are consistent: + // double direct = (mccTwo(mcc.h1, mcc.h2, eBefore, 0) / mccTwo(0, 0, eBefore, 0).Re()).Re(); + // double recur = twoRecursionBefore.Re(); + // double diff = std::abs((direct - recur) / recur); + + // if(mccTwo(0, 0, eBefore, 0).Re() < 0. || wTwoRecursionBefore < 0. || diff > 0.01) { + // LOGF(warning, "Compare <2>: cent=%f, n=%d, direct=%e, recursion=%e, diff=%e, denDirect=%e, denRec=%e", rlCollisionCent, mcc.h2, direct, recur, diff, mccTwo(0, 0, eBefore, 0).Re(), wTwoRecursionBefore); + // } } - for (int i = 0; i < mcc.NumSC; i++) { + for (int i = 0; i < NumHarmonics; i++) { // 4-p correlations: - mcc.h1 = -(i + 3); + mcc.h1 = -(i + 2); mcc.h2 = -2; mcc.h3 = 2; - mcc.h4 = i + 3; + mcc.h4 = i + 2; 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(); + // Before cut: + TComplex fourRecursionBefore = mccRecursion(4, harmonicsFourNum, eBefore, 0) / mccRecursion(4, harmonicsFourDen, eBefore, 0).Re(); + double wFourRecursionBefore = mccRecursion(4, harmonicsFourDen, eBefore, 0).Re(); ebye.fFourParticleCorrelationEbye[eBefore][i] = fourRecursionBefore.Re(); - TComplex fourRecursionAfter = mccRecursion(4, harmonicsFourNum, eAfter) / mccRecursion(4, harmonicsFourDen, eAfter).Re(); - double wFourRecursionAfter = mccRecursion(4, harmonicsFourDen, eAfter).Re(); - ebye.fFourParticleCorrelationEbye[eAfter][i] = fourRecursionAfter.Re(); + if (nTracksBefore > mcc.MaxCorrelator && wFourRecursionBefore > 0.) { + // Fill v2^2vn^2 profile and histogram: + mc.fFourParticleCorrelationProfiles[eBefore]->Fill(i + 2.5, ebye.fFourParticleCorrelationEbye[eBefore][i], wFourRecursionBefore); + mc.fFourParticleCorrelationHistograms[eBefore][eBefore][i]->Fill(ebye.fFourParticleCorrelationEbye[eBefore][i]); + mc.fFourParticleCorrelationHistograms[eBefore][eAfter][i]->Fill(ebye.fFourParticleCorrelationEbye[eBefore][i], wFourRecursionBefore); - if (nTracksBefore > mcc.MaxCorrelator) { - mc.fFourParticleCorrelationProfiles[eBefore]->Fill(i + 3.5, ebye.fFourParticleCorrelationEbye[eBefore][i], wFourRecursionBefore); + // Find max and min: + if (ebye.fFourParticleCorrelationEbye[eBefore][i] > ebye.fFourParticleCorrelationMaxEbye[eBefore][i]) { + ebye.fFourParticleCorrelationMaxEbye[eBefore][i] = ebye.fFourParticleCorrelationEbye[eBefore][i]; + LOGF(info, "max v2^2v%d^2 so far = %e", i + 2, ebye.fFourParticleCorrelationEbye[eBefore][i]); + } + if (ebye.fFourParticleCorrelationEbye[eBefore][i] < ebye.fFourParticleCorrelationMinEbye[eBefore][i]) { + ebye.fFourParticleCorrelationMinEbye[eBefore][i] = ebye.fFourParticleCorrelationEbye[eBefore][i]; + LOGF(info, "min v2^2v%d^2 so far = %e", i + 2, ebye.fFourParticleCorrelationEbye[eBefore][i]); + } + } else { + LOGF(warning, "cent=%f, nTracksBefore = %d, wFourRecursionBefore = %e", rlCollisionCent, nTracksBefore, wFourRecursionBefore); } - if (nTracksAfter > mcc.MaxCorrelator) { - mc.fFourParticleCorrelationProfiles[eAfter]->Fill(i + 3.5, ebye.fFourParticleCorrelationEbye[eAfter][i], wFourRecursionAfter); + + // After cut: + TComplex fourRecursionAfter = mccRecursion(4, harmonicsFourNum, eAfter, 0) / mccRecursion(4, harmonicsFourDen, eAfter, 0).Re(); + double wFourRecursionAfter = mccRecursion(4, harmonicsFourDen, eAfter, 0).Re(); + ebye.fFourParticleCorrelationEbye[eAfter][i] = fourRecursionAfter.Re(); + + if (nTracksAfter > mcc.MaxCorrelator && wFourRecursionAfter > 0.) { + // Fill v2^2vn^2 profile and histogram: + mc.fFourParticleCorrelationProfiles[eAfter]->Fill(i + 2.5, ebye.fFourParticleCorrelationEbye[eAfter][i], wFourRecursionAfter); + mc.fFourParticleCorrelationHistograms[eAfter][eBefore][i]->Fill(ebye.fFourParticleCorrelationEbye[eAfter][i]); + mc.fFourParticleCorrelationHistograms[eAfter][eAfter][i]->Fill(ebye.fFourParticleCorrelationEbye[eAfter][i], wFourRecursionAfter); + + // Find max and min: + if (ebye.fFourParticleCorrelationEbye[eAfter][i] > ebye.fFourParticleCorrelationMaxEbye[eAfter][i]) { + ebye.fFourParticleCorrelationMaxEbye[eAfter][i] = ebye.fFourParticleCorrelationEbye[eAfter][i]; + LOGF(info, "max v2^2v%d^2 so far = %e", i + 2, ebye.fFourParticleCorrelationEbye[eAfter][i]); + } + if (ebye.fFourParticleCorrelationEbye[eAfter][i] < ebye.fFourParticleCorrelationMinEbye[eAfter][i]) { + ebye.fFourParticleCorrelationMinEbye[eAfter][i] = ebye.fFourParticleCorrelationEbye[eAfter][i]; + LOGF(info, "min v2^2v%d^2 so far = %e", i + 2, ebye.fFourParticleCorrelationEbye[eAfter][i]); + } + } else { + LOGF(warning, "cent=%f, nTracksAfter = %d, wFourRecursionAfter = %e", rlCollisionCent, nTracksAfter, wFourRecursionAfter); } } } @@ -1208,10 +1451,10 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam float minPt = lPtBins[1]; float maxPt = lPtBins[2]; - for (int ba = 0; ba < eCutBeforeAfter_N; ba++) { + for (int ba = 0; ba < eBeforeAfter_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 nameRec = Form("fHist%s[eRec][%s cut]", ParticleHistNames[histType], BeforeAfterNames[ba]); + std::string nameSim = Form("fHist%s[eSim][%s cut]", ParticleHistNames[histType], BeforeAfterNames[ba]); std::string nameRecfull = Form("%s distribution for reconstructed particles", ParticleHistNames[histType]); std::string nameSimfull = Form("%s distribution for simulated particles", ParticleHistNames[histType]); @@ -1237,10 +1480,10 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam float minCent = lCentBins[1]; float maxCent = lCentBins[2]; - for (int ba = 0; ba < eCutBeforeAfter_N; ba++) { + for (int ba = 0; ba < eBeforeAfter_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], BeforeAfterNames[ba]); + std::string nameSim = Form("fHist%s[eSim][%s cut]", EventHistNames[histType], BeforeAfterNames[ba]); std::string nameRecfull; std::string nameSimfull; @@ -1300,8 +1543,8 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam maxCent = lCentBins[2]; } - for (int ba = 0; ba < eCutBeforeAfter_N; ba++) { - std::string name = Form("fHist%s[%s cut]", EventHistNames[histType], CutBeforeAfterNames[ba]); + for (int ba = 0; ba < eBeforeAfter_N; ba++) { + std::string name = Form("fHist%s[%s cut]", EventHistNames[histType], BeforeAfterNames[ba]); std::string namefull; if constexpr (histType == eCent) { @@ -1329,14 +1572,56 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam void bookCorrHistograms(T1 const& lCrBins) { - 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]; + // book limit functions: + if constexpr (histType == eCorrCent) { + cr.fCorrBounds[histType][0] = new TF1(Form("fCorrUpperBound%s", CorrHistNames[histType]), "[0] * x + [1]"); + cr.fCorrBounds[histType][1] = new TF1(Form("fCorrLowerBound%s", CorrHistNames[histType]), "(x - [1]) / [0]"); + cr.fCorrBounds[histType][0]->SetTitle(Form("%s upper bound", CorrHistNames[histType])); + cr.fCorrBounds[histType][1]->SetTitle(Form("%s lower bound", CorrHistNames[histType])); + } else if constexpr (histType == eCorrMult) { + for (int i = 0; i < eEstmCorr_N; i++) { + cr.fCorrBounds[histType + i][0] = new TF1(Form("fCorrUpperBound%s%s", EstmCorrName[i], CorrHistNames[histType]), "[0] * x + [1]"); + cr.fCorrBounds[histType + i][1] = new TF1(Form("fCorrLowerBound%s%s", EstmCorrName[i], CorrHistNames[histType]), "[0] * x + [1]"); + cr.fCorrBounds[histType + i][0]->SetTitle(Form("%s %s upper bound", EstmCorrName[i], CorrHistNames[histType])); + cr.fCorrBounds[histType + i][1]->SetTitle(Form("%s %s lower bound", EstmCorrName[i], CorrHistNames[histType])); + } + } + + int nBinsCent = 0; + float minCent = 0.; + float maxCent = 0.; + int nBinsXMult = 0; + float minXMult = 0.; + float maxXMult = 0.; + int nBinsYMult = 0; + float minYMult = 0.; + float maxYMult = 0.; + + if constexpr (histType == eCorrCent) { - for (int ba = 0; ba < eCutBeforeAfter_N; ba++) { + // (cent) get bins: + const auto& lCentBins = lCrBins[histType]; + nBinsCent = static_cast(lCentBins[0]); + minCent = lCentBins[1]; + maxCent = lCentBins[2]; + + // (cent) set functions and add to list: + cr.fCorrBounds[histType][0]->SetRange(minCent, maxCent); + cr.fCorrBounds[histType][1]->SetRange(minCent, maxCent); + cr.fCorrBounds[histType][0]->SetMinimum(0); + cr.fCorrBounds[histType][1]->SetMinimum(0); + cr.fCorrBounds[histType][0]->SetMaximum(100); + cr.fCorrBounds[histType][1]->SetMaximum(100); + cr.fCorrBounds[histType][0]->SetParameters(tc.fCentCorrCut[0], tc.fCentCorrCut[1]); + cr.fCorrBounds[histType][1]->SetParameters(tc.fCentCorrCut[0], tc.fCentCorrCut[1]); + cr.fCorrHistogramsList->Add(cr.fCorrBounds[histType][0]); + cr.fCorrHistogramsList->Add(cr.fCorrBounds[histType][1]); + + } // else if constexpr (histType == eCorrMult) {} + + for (int ba = 0; ba < eBeforeAfter_N; ba++) { std::string name; - std::string namefull = Form("%s correlation %s cut", CorrHistNames[histType], CutBeforeAfterNames[ba]); + std::string namefull = Form("%s correlation %s cut", CorrHistNames[histType], BeforeAfterNames[ba]); int nEstm = 0; if constexpr (histType == eCorrCent) { @@ -1344,33 +1629,67 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam } else if constexpr (histType == eCorrMult) { nEstm = eMultEstm_N; } + + int nComb = 0; // count combination AB, AC, BC + for (int i = 0; i < nEstm; i++) { std::string titleX; - if (histType == eCorrCent) { + if constexpr (histType == eCorrCent) { titleX = Form("%s %s", CentEstmNames[i], CorrHistNames[histType]); - } else if (histType == eCorrMult) { + + } else if constexpr (histType == eCorrMult) { titleX = Form("%s %s", MultEstmNames[i], CorrHistNames[histType]); - // x axis bin... } for (int j = i + 1; j < nEstm; j++) { std::string titleY; - if (histType == eCorrCent) { - name = Form("fHist%s[%s vs. %s][%s cut]", CorrHistNames[histType], CentEstmNames[i], CentEstmNames[j], CutBeforeAfterNames[ba]); + if constexpr (histType == eCorrCent) { + name = Form("fHist%s[%s vs. %s][%s cut]", CorrHistNames[histType], CentEstmNames[i], CentEstmNames[j], BeforeAfterNames[ba]); titleY = Form("%s %s", CentEstmNames[j], CorrHistNames[histType]); + + // (cent) book corr 2D histogram: 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]); + + } else if constexpr (histType == eCorrMult) { + name = Form("fHist%s[%s vs. %s][%s cut]", CorrHistNames[histType], MultEstmNames[i], MultEstmNames[j], BeforeAfterNames[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); + + // (mult) get bins: + const auto& lMultXBins = lCrBins[i + 1]; + nBinsXMult = static_cast(lMultXBins[0]); + minXMult = lMultXBins[1]; + maxXMult = lMultXBins[2]; + const auto& lMultYBins = lCrBins[j + 1]; + nBinsYMult = static_cast(lMultYBins[0]); + minYMult = lMultYBins[1]; + maxYMult = lMultYBins[2]; + + // (mult) set functions and add to list: + if (static_cast(ba)) { + cr.fCorrBounds[histType + nComb][0]->SetRange(minXMult, maxXMult); + cr.fCorrBounds[histType + nComb][1]->SetRange(minXMult, maxXMult); + cr.fCorrBounds[histType + nComb][0]->SetParameters(tc.fMultAllCorrCut[i][j][0], tc.fMultAllCorrCut[i][j][1]); + cr.fCorrBounds[histType + nComb][1]->SetParameters(tc.fMultAllCorrCut[i][j][2], tc.fMultAllCorrCut[i][j][3]); + cr.fCorrBounds[histType + nComb][0]->SetMinimum(0); + cr.fCorrBounds[histType + nComb][1]->SetMinimum(0); + cr.fCorrBounds[histType + nComb][0]->SetMaximum(maxYMult); + cr.fCorrBounds[histType + nComb][1]->SetMaximum(maxYMult); + cr.fCorrHistogramsList->Add(cr.fCorrBounds[histType + nComb][0]); + cr.fCorrHistogramsList->Add(cr.fCorrBounds[histType + nComb][1]); + } + + // (mult) book corr 2D histogram: + cr.fCorrHistograms[histType][i][j][ba] = new TH2F(name.c_str(), namefull.c_str(), nBinsXMult, minXMult, maxXMult, nBinsYMult, minYMult, maxYMult); } + // (cent/mult) set corr histogram cr.fCorrHistograms[histType][i][j][ba]->GetYaxis()->SetTitle(titleY.c_str()); cr.fCorrHistograms[histType][i][j][ba]->GetXaxis()->SetTitle(titleX.c_str()); cr.fCorrHistogramsList->Add(cr.fCorrHistograms[histType][i][j][ba]); + + nComb += 1; } } } @@ -1393,6 +1712,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam tc.fCentEstm = cfCentEstm; tc.fMultEstm = cfMultEstm; + tc.fGeneralCutSwitch = cfGeneralCutSwitch; tc.fVertexZCutSwitch = cfVertexZCutSwitch; tc.fSel8CutSwitch = cfSel8CutSwitch; tc.fCentCutSwitch = cfCentCutSwitch; @@ -1406,12 +1726,15 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam tc.fTpcNClsFoundCutSwitch = cfTpcNClsFoundCutSwitch; tc.fDCAXYCutSwitch = cfDCAXYCutSwitch; tc.fDCAZCutSwitch = cfDCAZCutSwitch; + tc.fEtaGapSwitch = cfEtaGapSwitch; tc.fVertexZCut = cfVertexZCut; tc.fCentCut = cfCentCut; tc.fNumContribCut = cfNumContribCut; tc.fCentCorrCut = cfCentCorrCut; - tc.fMultCorrCut = cfMultCorrCut; + tc.fMultAllCorrCut[0][1] = cfMultABCorrCut; + tc.fMultAllCorrCut[0][2] = cfMultACCorrCut; + tc.fMultAllCorrCut[1][2] = cfMultBCCorrCut; tc.fPtCut = cfPtCut; tc.fEtaCut = cfEtaCut; @@ -1419,16 +1742,22 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam tc.fTpcNClsFoundCut = cfTpcNClsFoundCut; tc.fDCAXYCut = cfDCAXYCut; tc.fDCAZCut = cfDCAZCut; + tc.fEtaGap = cfEtaGap; tc.fPtBins = cfPtBins; tc.fPhiBins = cfPhiBins; tc.fCentBins = cfCentBins; - tc.fMultBins = cfMultBins; + tc.fFT0CMultBins = cfFT0CMultBins; + tc.fFT0MMultBins = cfFT0MMultBins; + tc.fFV0AMultBins = cfFV0AMultBins; + tc.fVerXBins = cfVerXBins; tc.fVerYBins = cfVerYBins; tc.fVerZBins = cfVerZBins; tc.fNumContribBins = cfNumContribBins; + tc.fTwoParticleCorrBins = cfTwoParticleCorrBins; + tc.fFourParticleCorrBins = cfFourParticleCorrBins; tc.fPrintSwitch = cfPrintSwitch; tc.fFileWithWeights = cfFileWithWeights; @@ -1436,6 +1765,16 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam qa.fQASwitch = cfQASwitch; wt.fWeightSwitch = cfWeightSwitch; + if (tc.fMultEstm == "FT0C") { + tc.fMultBins = tc.fFT0CMultBins; + } else if (tc.fMultEstm == "FT0M") { + tc.fMultBins = tc.fFT0MMultBins; + } else if (tc.fMultEstm == "FV0A") { + tc.fMultBins = tc.fFV0AMultBins; + } else { + LOGF(fatal, "Invalid multiplicity estimator: %s", tc.fMultEstm.c_str()); + } + // *) Book base list: auto* temp = new TList(); temp->SetOwner(kTRUE); @@ -1475,7 +1814,7 @@ struct MultiparticleCumulants { // this name is used in lower-case format to nam std::vector> lPcBins = {tc.fPtBins, tc.fPhiBins}; std::vector> lEvBins = {tc.fCentBins, tc.fMultBins, tc.fVerXBins, tc.fVerYBins, tc.fVerZBins, tc.fNumContribBins}; std::vector> lQABins = {tc.fCentBins, tc.fMultBins, tc.fNumContribBins}; - std::vector> lCrBins = {tc.fCentBins, tc.fMultBins}; + std::vector> lCrBins = {tc.fCentBins, tc.fFT0CMultBins, tc.fFT0MMultBins, tc.fFV0AMultBins}; bookParticleHistograms(lPcBins); bookParticleHistograms(lPcBins);