diff --git a/PWGJE/DataModel/PhotonChargedTriggerCorrelation.h b/PWGJE/DataModel/PhotonChargedTriggerCorrelation.h index 9b3eef803a8..637a7f18c35 100644 --- a/PWGJE/DataModel/PhotonChargedTriggerCorrelation.h +++ b/PWGJE/DataModel/PhotonChargedTriggerCorrelation.h @@ -43,11 +43,13 @@ DECLARE_SOA_COLUMN(Eta, eta, float); namespace collision_extra_corr { DECLARE_SOA_COLUMN(SelEv, selEv, bool); +DECLARE_SOA_COLUMN(ThirdDecimalTruePosZ, thirdDecimalTruePosZ, int); DECLARE_SOA_COLUMN(TrigEv, trigEv, bool); +DECLARE_SOA_COLUMN(PtMax, ptMax, float); DECLARE_SOA_COLUMN(NGlobalTracks, nGlobalTracks, int); } // namespace collision_extra_corr DECLARE_SOA_TABLE(CollisionsExtraCorr, "AOD", "COLLISIONSEXTRACORR", - collision_extra_corr::SelEv, collision_extra_corr::TrigEv, collision_extra_corr::NGlobalTracks); + collision_extra_corr::SelEv, collision_extra_corr::ThirdDecimalTruePosZ, collision_extra_corr::TrigEv, collision_extra_corr::PtMax, collision_extra_corr::NGlobalTracks); // trigger namespace trigger @@ -83,12 +85,12 @@ using Pipm = Pipms::iterator; namespace photon_pcm { DECLARE_SOA_INDEX_COLUMN_FULL(V0PhotonKF, v0PhotonKF, int, V0PhotonsKF, ""); -DECLARE_SOA_COLUMN(PosTrackId, posTrackId, int); -DECLARE_SOA_COLUMN(NegTrackId, negTrackId, int); +DECLARE_SOA_INDEX_COLUMN_FULL(PosJetTrack, posJetTrack, int, JetTracks, "_pos"); +DECLARE_SOA_INDEX_COLUMN_FULL(NegJetTrack, negJetTrack, int, JetTracks, "_neg"); } // namespace photon_pcm DECLARE_SOA_TABLE(PhotonPCMs, "AOD", "PHOTONPCMS", o2::soa::Index<>, corr_particle::JetCollisionId, photon_pcm::V0PhotonKFId, - photon_pcm::PosTrackId, photon_pcm::NegTrackId, + photon_pcm::PosJetTrackId, photon_pcm::NegJetTrackId, corr_particle::Pt, corr_particle::Phi, corr_particle::Eta); using PhotonPCM = PhotonPCMs::iterator; @@ -97,15 +99,15 @@ namespace photon_pcm_pair { DECLARE_SOA_INDEX_COLUMN_FULL(V0PhotonKF1, v0PhotonKF1, int, V0PhotonsKF, "_1"); DECLARE_SOA_INDEX_COLUMN_FULL(V0PhotonKF2, v0PhotonKF2, int, V0PhotonsKF, "_2"); -DECLARE_SOA_COLUMN(PosTrack1Id, posTrack1Id, int); -DECLARE_SOA_COLUMN(NegTrack1Id, negTrack1Id, int); -DECLARE_SOA_COLUMN(PosTrack2Id, posTrack2Id, int); -DECLARE_SOA_COLUMN(NegTrack2Id, negTrack2Id, int); +DECLARE_SOA_INDEX_COLUMN_FULL(PosJetTrack1, posJetTrack1, int, JetTracks, "_pos1"); +DECLARE_SOA_INDEX_COLUMN_FULL(NegJetTrack1, negJetTrack1, int, JetTracks, "_neg1"); +DECLARE_SOA_INDEX_COLUMN_FULL(PosJetTrack2, posJetTrack2, int, JetTracks, "_pos2"); +DECLARE_SOA_INDEX_COLUMN_FULL(NegJetTrack2, negJetTrack2, int, JetTracks, "_neg2"); DECLARE_SOA_COLUMN(Mgg, mgg, float); } // namespace photon_pcm_pair DECLARE_SOA_TABLE(PhotonPCMPairs, "AOD", "PHOTONPCMPAIRS", o2::soa::Index<>, corr_particle::JetCollisionId, photon_pcm_pair::V0PhotonKF1Id, photon_pcm_pair::V0PhotonKF2Id, - photon_pcm_pair::PosTrack1Id, photon_pcm_pair::NegTrack1Id, photon_pcm_pair::PosTrack2Id, photon_pcm_pair::NegTrack2Id, + photon_pcm_pair::PosJetTrack1Id, photon_pcm_pair::NegJetTrack1Id, photon_pcm_pair::PosJetTrack2Id, photon_pcm_pair::NegJetTrack2Id, corr_particle::Pt, corr_particle::Phi, corr_particle::Eta, photon_pcm_pair::Mgg); using PhotonPCMPair = PhotonPCMPairs::iterator; @@ -114,11 +116,13 @@ using PhotonPCMPair = PhotonPCMPairs::iterator; // mcCollision extension namespace mc_collision_extra_corr { +DECLARE_SOA_COLUMN(ThirdDecimalTruePosZ, thirdDecimalTruePosZ, int); DECLARE_SOA_COLUMN(TrigEv, trigEv, bool); +DECLARE_SOA_COLUMN(PtMax, ptMax, float); DECLARE_SOA_COLUMN(NChargedInEtaRange, nChargedInEtaRange, int); } // namespace mc_collision_extra_corr DECLARE_SOA_TABLE(McCollisionsExtraCorr, "AOD", "MCCOLLISIONSEXTRACORR", - mc_collision_extra_corr::TrigEv, mc_collision_extra_corr::NChargedInEtaRange); + mc_collision_extra_corr::ThirdDecimalTruePosZ, mc_collision_extra_corr::TrigEv, mc_collision_extra_corr::PtMax, mc_collision_extra_corr::NChargedInEtaRange); // trigger namespace trigger_particle diff --git a/PWGJE/Tasks/photonChargedTriggerCorrelation.cxx b/PWGJE/Tasks/photonChargedTriggerCorrelation.cxx index 81cd439194e..57d8f1c2b47 100644 --- a/PWGJE/Tasks/photonChargedTriggerCorrelation.cxx +++ b/PWGJE/Tasks/photonChargedTriggerCorrelation.cxx @@ -30,6 +30,7 @@ #include #include +#include #include #include #include @@ -52,13 +53,18 @@ #include #include #include +#include #include #include #include #include #include +#include #include +#include #include +#include +#include #include const double absEtaMaxDefault = 0.8; @@ -77,108 +83,109 @@ using CorrMcCollision = CorrMcCollisions::iterator; using BinningZPvMult = ColumnBinningPolicy; -// correlation analysis ======================================================================================================================================================================= +// correlation analysis /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// =/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/ struct PhotonChargedTriggerCorrelation { // configurables // general (kenobi) - Configurable pathCcdbEff{"pathCcdbEff", "Users/j/jkinner/efficiency/set_in_config", "base path to the ccdb efficiencies"}; + Configurable pathCcdbCorrection{"pathCcdbCorrection", "Users/j/jkinner/correction/set_in_config", "base path to the ccdb corrections"}; Configurable urlCcdb{"urlCcdb", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable noLaterThanCcdbConfig{"noLaterThanCcdbConfig", -1, "latest acceptable timestamp of creation for the object (-1 for task start time)"}; + Configurable splitMcEvN{"splitMcEvN", 1, "number of equal event fractions for selection in MC (max 10) | (select 1 for data)"}; + Configurable splitMcEvSelect{"splitMcEvSelect", 1, "selected event fraction number in MC (starts at 1) | (select 1 for data)"}; // analysis - Configurable doEffCorrectionTrigger{"doEffCorrectionTrigger", false, "whether to do on-the-fly mixing correction for triggers"}; - Configurable doEffCorrectionHadron{"doEffCorrectionHadron", false, "whether to do on-the-fly mixing correction for hadrons"}; - Configurable doEffCorrectionPipm{"doEffCorrectionPipm", false, "whether to do on-the-fly mixing correction for pipm"}; - Configurable doEffCorrectionPhotonPCM{"doEffCorrectionPhotonPCM", false, "whether to do on-the-fly mixing correction for photonPCM"}; - + Configurable applyCTotTrigger{"applyCTotTrigger", false, "whether to apply on-the-fly total correction factor for triggers"}; Configurable doTrigEvMixing{"doTrigEvMixing", false, "whether to use trigger events for trigger mixing"}; - Configurable nTriggerSavedForMixing{"nTriggerSavedForMixing", 2048, "number of triggers that are saved for mixing with other events"}; - Configurable nTriggerMixingMcTrue{"nTriggerMixingMcTrue", 8, "number of triggers that are used for mc true mixing"}; - Configurable nTriggerMixingHadron{"nTriggerMixingHadron", 64, "number of triggers that are used for hadron mixing"}; - Configurable nTriggerMixingPipm{"nTriggerMixingPipm", 64, "number of triggers that are used for pipm mixing"}; - Configurable nTriggerMixingPhotonPCM{"nTriggerMixingPhotonPCM", 256, "number of triggers that are saved for photonPCM mixing"}; - Configurable nTriggerMixingH0PCM{"nTriggerMixingH0PCM", 256, "number of triggers that are saved for h0PCM (pi0, eta) mixing"}; + + Configurable nTriggerSavedForMixing{"nTriggerSavedForMixing", 8192, "number of triggers that are saved for mixing with other events"}; + Configurable nTriggerBinMinThreshold{"nTriggerBinMinThreshold", 128, "threshold minimum number of triggers in mixing bin required"}; + Configurable iTriggerStartFraction{"iTriggerStartFraction", 0.5, "ratio of bin size for randomised iTrigger start"}; + Configurable nTriggerScaleCoefficient{"nTriggerScaleCoefficient", 2, "coeffcient a in mixing-number pt power scaling"}; + Configurable nTriggerScaleExponent{"nTriggerScaleExponent", 2, "exponent a in mixing-number pt power scaling"}; + Configurable nMixingAt0McTrue{"nMixingAt0McTrue", 4, "number of triggers that are used for mc true mixing"}; + Configurable nMixingAt0Hadron{"nMixingAt0Hadron", 16, "number of triggers that are used for hadron mixing"}; + Configurable nMixingAt0Pipm{"nMixingAt0Pipm", 16, "number of triggers that are used for pipm mixing"}; + Configurable nMixingAt0PhotonPCM{"nMixingAt0PhotonPCM", 256, "number of triggers that are saved for photonPCM mixing"}; + Configurable nMixingAt0H0PCM{"nMixingAt0H0PCM", 1024, "number of triggers that are saved for h0PCM (pi0, eta) mixing"}; + Configurable nNeighboursMixingPhotonPCMPair{"nNeighboursMixingPhotonPCMPair", 32, "number neighbours used for for photonPCM pair mixing"}; Configurable> pi0PCMPeakMassRange{"pi0PCMPeakMassRange", {0.10, 0.15}, "photon-pair mass integration range for pi0PCM"}; - Configurable> pi0PCMSideMassRange{"pi0PCMSideMassRange", {0.16, 0.24}, "photon-pair mass integration range outside pi0PCM region"}; + Configurable> pi0PCMSideMassRange{"pi0PCMSideMassRange", {0.16, 0.22}, "photon-pair mass integration range outside pi0PCM region"}; Configurable> etaPCMPeakMassRange{"etaPCMPeakMassRange", {0.51, 0.56}, "photon-pair mass integration range for etaPCM"}; - Configurable> etaPCMLowSideMassRange{"etaPCMLowSideMassRange", {0.45, 0.50}, "photon-pair mass integration range below etaPCM region"}; - Configurable> etaPCMHighSideMassRange{"etaPCMHighSideMassRange", {0.56, 0.65}, "photon-pair mass integration range above etaPCM region"}; - - Configurable doTrigEvEff{"doTrigEvEff", false, "whether to use trigger events for efficiency histograms"}; - Configurable ptCutTrigEvEff{"ptCutTrigEvEff", 4, "pT cut for efficieny calculation in trigger events (to avoid trigger bias)"}; - Configurable requireSingleCollisionPurity{"requireSingleCollisionPurity", true, "whether particle from single chosen MC-col associated to reco-col (else just type/kin match)"}; + Configurable> etaPCMLowSideMassRange{"etaPCMLowSideMassRange", {0.44, 0.48}, "photon-pair mass integration range below etaPCM region"}; + Configurable> etaPCMHighSideMassRange{"etaPCMHighSideMassRange", {0.58, 0.62}, "photon-pair mass integration range above etaPCM region"}; // for histograms - Configurable nBinsZPv{"nBinsZPv", 100, "number zPv bins in histos for QA"}; - Configurable nBinsZPvSmol{"nBinsZPvSmol", 28, "number zPv bins but smaller"}; - Configurable nBinsMult{"nBinsMult", 200, "number multiplicity bins in histos for QA"}; - Configurable nBinsMultSmol{"nBinsMultSmol", 20, "number multiplicity bins but smaller"}; + Configurable nBinsZPv{"nBinsZPv", 28, "number zPv bins in histos"}; + Configurable nBinsMult{"nBinsMult", 64, "number multiplicity bins in histos"}; Configurable nBinsOccupancy{"nBinsOccupancy", 2000, "number occupancy bins in histos for QA"}; Configurable nBinsPhi{"nBinsPhi", 72, "number phi bins"}; - Configurable nBinsEta{"nBinsEta", 40, "number eta bins"}; - Configurable nBinsMgg{"nBinsMgg", 160, "number mass-photon-pair bins"}; + Configurable nBinsEta{"nBinsEta", 32, "number eta bins"}; + Configurable nBinsDCAz{"nBinsDCAz", 100, "number DCAz bins"}; + Configurable nBinsMgg{"nBinsMgg", 240, "number mass-photon-pair bins"}; - Configurable> binsPtTrig{"binsPtTrig", {5, 10, 25, 50}, "correlation ptTrig bins"}; + Configurable> binsPtTrig{"binsPtTrig", {0, 5, 10, 25, 50, 100}, "correlation ptTrig bins"}; Configurable> binsPtAssoc{"binsPtAssoc", - {0.2, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10, 12.5, 15, 17.5, 20, 30, 40}, + {0, 0.2, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10, + 12.5, 15, 17.5, 20, 25, 35, 50, 100}, "correlation ptAssoc bins"}; Configurable> binsDPhi{"binsDPhi", {0.00 * DPHI_SCALE, - 0.04 * DPHI_SCALE, 0.08 * DPHI_SCALE, 0.11 * DPHI_SCALE, 0.14 * DPHI_SCALE, - 0.16 * DPHI_SCALE, 0.18 * DPHI_SCALE, 0.20 * DPHI_SCALE, 0.22 * DPHI_SCALE, - 0.23 * DPHI_SCALE, 0.24 * DPHI_SCALE, 0.25 * DPHI_SCALE, 0.26 * DPHI_SCALE, 0.27 * DPHI_SCALE, 0.28 * DPHI_SCALE, - 0.30 * DPHI_SCALE, 0.32 * DPHI_SCALE, 0.34 * DPHI_SCALE, 0.36 * DPHI_SCALE, - 0.39 * DPHI_SCALE, 0.42 * DPHI_SCALE, 0.46 * DPHI_SCALE, 0.50 * DPHI_SCALE, - 0.54 * DPHI_SCALE, 0.58 * DPHI_SCALE, 0.61 * DPHI_SCALE, 0.64 * DPHI_SCALE, - 0.66 * DPHI_SCALE, 0.68 * DPHI_SCALE, 0.70 * DPHI_SCALE, 0.72 * DPHI_SCALE, - 0.74 * DPHI_SCALE, 0.76 * DPHI_SCALE, 0.78 * DPHI_SCALE, - 0.80 * DPHI_SCALE, 0.82 * DPHI_SCALE, 0.84 * DPHI_SCALE, 0.86 * DPHI_SCALE, - 0.89 * DPHI_SCALE, 0.92 * DPHI_SCALE, 0.96 * DPHI_SCALE, 1.00 * DPHI_SCALE}, + 0.05 * DPHI_SCALE, 0.10 * DPHI_SCALE, 0.14 * DPHI_SCALE, 0.17 * DPHI_SCALE, 0.20 * DPHI_SCALE, + 0.22 * DPHI_SCALE, 0.24 * DPHI_SCALE, 0.26 * DPHI_SCALE, 0.28 * DPHI_SCALE, 0.30 * DPHI_SCALE, + 0.33 * DPHI_SCALE, 0.36 * DPHI_SCALE, 0.40 * DPHI_SCALE, 0.45 * DPHI_SCALE, 0.50 * DPHI_SCALE, + 0.55 * DPHI_SCALE, 0.60 * DPHI_SCALE, 0.64 * DPHI_SCALE, 0.68 * DPHI_SCALE, 0.71 * DPHI_SCALE, + 0.74 * DPHI_SCALE, 0.76 * DPHI_SCALE, 0.79 * DPHI_SCALE, + 0.82 * DPHI_SCALE, 0.86 * DPHI_SCALE, 0.90 * DPHI_SCALE, 0.95 * DPHI_SCALE, 1.00 * DPHI_SCALE}, "correlation bins DeltaPhi"}; Configurable> binsDEta{"binsDEta", {0 / 32. * DETA_SCALE, - 1 / 32. * DETA_SCALE, 2 / 32. * DETA_SCALE, 3 / 32. * DETA_SCALE, 4 / 32. * DETA_SCALE, - 5 / 32. * DETA_SCALE, 6 / 32. * DETA_SCALE, 7 / 32. * DETA_SCALE, 8 / 32. * DETA_SCALE, - 9 / 32. * DETA_SCALE, 10 / 32. * DETA_SCALE, 11 / 32. * DETA_SCALE, 12 / 32. * DETA_SCALE, 13 / 32. * DETA_SCALE, 14 / 32. * DETA_SCALE, - 59 / 128. * DETA_SCALE, 62 / 128. * DETA_SCALE, 64 / 128. * DETA_SCALE, 66 / 128. * DETA_SCALE, 69 / 128. * DETA_SCALE, 18 / 32. * DETA_SCALE, - 19 / 32. * DETA_SCALE, 20 / 32. * DETA_SCALE, 21 / 32. * DETA_SCALE, 22 / 32. * DETA_SCALE, 23 / 32. * DETA_SCALE, 24 / 32. * DETA_SCALE, - 25 / 32. * DETA_SCALE, 26 / 32. * DETA_SCALE, 27 / 32. * DETA_SCALE, 28 / 32. * DETA_SCALE, - 29 / 32. * DETA_SCALE, 30 / 32. * DETA_SCALE, 31 / 32. * DETA_SCALE, 32 / 32. * DETA_SCALE}, + 2 / 32. * DETA_SCALE, 4 / 32. * DETA_SCALE, 6 / 32. * DETA_SCALE, 8 / 32. * DETA_SCALE, + 9.5 / 32. * DETA_SCALE, 11 / 32. * DETA_SCALE, 12.5 / 32. * DETA_SCALE, 14 / 32. * DETA_SCALE, + 15.5 / 32. * DETA_SCALE, 16.5 / 32. * DETA_SCALE, 18 / 32. * DETA_SCALE, + 19.5 / 32. * DETA_SCALE, 21 / 32. * DETA_SCALE, 22.5 / 32. * DETA_SCALE, 24 / 32. * DETA_SCALE, + 26 / 32. * DETA_SCALE, 28 / 32. * DETA_SCALE, 30 / 32. * DETA_SCALE, 32 / 32. * DETA_SCALE}, "correlation bins DeltaEta"}; - Configurable> binsZPv{"binsZPv", - {-7, -5, -3, -1, 1, 3, 5, 7}, - "zPv mixing bins"}; - Configurable> binsMult{"binsMult", - {-0.5, 9.5, 14.5, 19.5, 25.5, 32}, - "multiplicity mixing bins for mc true"}; - Configurable> binsZPvMcTrue{"binsZPvMcTrue", - {-10000, 10000}, - "zPv mixing bins"}; - Configurable> binsMultMcTrue{"binsMultMcTrue", - {-0.5, 10.5, 15.5, 20.5, 25.5, 30.5, 35.5, 40.5, 50.5, 64}, - "multiplicity mixing bins for mc true"}; + Configurable> binsZPvBinning{"binsZPvBinning", + {-7, -5, -3, -1, 1, 3, 5, 7}, + "zPv mixing bins"}; + Configurable> binsMultBinning{"binsMultBinning", + {-0.5, 10.5, 15.5, 20.5, 27.5, 42.5}, + "multiplicity mixing bins"}; + Configurable> binsZPvBinningMcTrue{"binsZPvBinningMcTrue", + {-10000, 10000}, + "zPv mixing bins for mc true"}; + Configurable> binsMultBinningMcTrue{"binsMultBinningMcTrue", + {-0.5, 16.5, 24.5, 31.5, 39.5, 64.5}, + "multiplicity mixing bins for mc true"}; // configurables from other tasks - double etaMax; + double absEtaMax; + + // further variables + + std::mt19937 randomMt{3141u}; + int randomIntInInterval(int const lowEdgeInclusive, int const upEdgeInclusive) + { + return lowEdgeInclusive + (randomMt() % (upEdgeInclusive - lowEdgeInclusive + 1)); // speeds up task by factor > 2; modulo bias should not be relevant here + // return std::uniform_int_distribution(lowEdgeInclusive, upEdgeInclusive)(randomMt); + } // objects to hold histograms HistogramRegistry histos{"histogramRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; // ccdb calls - const int64_t noLaterThanCcdb = noLaterThanCcdbConfig == -1 ? std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count() : noLaterThanCcdbConfig; + const int64_t noLaterThanCcdb = + noLaterThanCcdbConfig == -1 ? std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count() : noLaterThanCcdbConfig; Service ccdb; // for mc Service pdg; - // random number generation - static constexpr unsigned int SeedRandomEngine = 12345; - std::mt19937 randomEngine{SeedRandomEngine}; - // partitions SliceCache cache; @@ -194,21 +201,18 @@ struct PhotonChargedTriggerCorrelation { // combinations binning // cumbersome, but still better than having extra configurable or figuring out how to init binningZPvMult later while declaring it here - std::function(std::vector const&, double)> prependValueToVector = + std::function(std::vector const&, double const)> const prependValueToVector = [](std::vector const& vec, double const value) { std::vector resultVec = {value}; resultVec.insert(resultVec.end(), vec.begin(), vec.end()); return resultVec; }; - BinningZPvMult binningZPvMult{{prependValueToVector(binsZPv.value, VARIABLE_WIDTH), prependValueToVector(binsMult.value, VARIABLE_WIDTH)}, true}; + BinningZPvMult binningZPvMult{{prependValueToVector(binsZPvBinning.value, VARIABLE_WIDTH), prependValueToVector(binsMultBinning.value, VARIABLE_WIDTH)}, true}; // declare analysis variables - // efficiency histograms - TH1D* h1PtInvEffTrigger; - TH1D* h1PtInvEffHadron; - TH1D* h1PtInvEffPipm; - TH1D* h1PtInvEffPhotonPCM; + // correction histograms + TH1D* h1PtCorrectionTrigger; // mixing trigger memory struct MixingTrigger { @@ -219,131 +223,149 @@ struct PhotonChargedTriggerCorrelation { }; // class to handle trigger info from previous collisions (beyond single dataframe) // organised as zPv- and mult-bin matrix of deque (pt, phi, eta) to save trigger info beyond single dataframe - // extra bin for mult overflow - // with adjusted zVtx (see triggerBinValuesZPv in init) and mult overflow -> all events accounted for + // last mult bin open to account for overflow and adjusted zVtx limits to account for rounding errors -> all events accounted for // (possibly replace by some advanced derived data method and O2 event mixing in future?) class MixingTriggerMemory { public: - // finds bin that value belongs to (assumes ordered bins) (starts at 0; includes underflow (return -1) and overlflow (return bins.size() - 1)) + // finds bin that value belongs to (assumes ordered bin edges) (starts at 0; includes underflow (return -1) and overflow (return binEdges.size() - 1)) // should be faster than some std binary search due to small number of bins (zPv, mult) - static int findIntervalBin(double value, const std::vector& bins) + static int findIntervalBin(double value, const std::vector& binEdges) { - const int n = bins.size() - 1; - if (value < bins[0]) + if (binEdges.empty()) { + throw std::invalid_argument("binEdges in findIntervalBin"); + } + const int n = binEdges.size() - 1; + if (value < binEdges[0]) return -1; // underflow for (int i_bin = 0; i_bin < n; i_bin++) - if (value < bins[i_bin + 1]) + if (value < binEdges[i_bin + 1]) return i_bin; return n; // overflow } - MixingTriggerMemory(int const nTriggerSavedForMixingIn, std::vector binsZPv, std::vector binsMult) + MixingTriggerMemory(std::vector const& binEdgesZPvIn, std::vector const& binEdgesMultIn, size_t const nTriggerPerBinLimitIn) + : binEdgesZPv(binEdgesZPvIn), + nBinsZPv(binEdgesZPvIn.size() - 1), + binEdgesMult(binEdgesMultIn), + nBinsMult(binEdgesMultIn.size() - 1), + nTriggerPerBinLimit(nTriggerPerBinLimitIn) { - nTriggerSavedForMixing = nTriggerSavedForMixingIn; - triggerBinValuesZPv = binsZPv; - triggerBinValuesMult = binsMult; - // prevent rounding errors in bin finding (multiplicity accounted for by it going to 0 and already considering overflow separately) - triggerBinValuesZPv.front() *= zPvRoundingErrorAdjust; - triggerBinValuesZPv.back() *= zPvRoundingErrorAdjust; + int const nEdgesMin = 2; + if (binEdgesZPv.size() < nEdgesMin || binEdgesMult.size() < nEdgesMin) { + throw std::invalid_argument("too few bin edges to define zPv or mult bins"); + } + if (!std::is_sorted(binEdgesZPv.begin(), binEdgesZPv.end()) || !std::is_sorted(binEdgesMult.begin(), binEdgesMult.end())) { + throw std::invalid_argument("bin edges must be sorted"); + } + // prevent zPv over/underflow due to rounding errors + binEdgesZPv.front() = binEdgesZPv.front() < 0 ? binEdgesZPv.front() * zPvRoundingErrorAdjust : binEdgesZPv.front() / zPvRoundingErrorAdjust; + binEdgesZPv.back() = binEdgesZPv.back() > 0 ? binEdgesZPv.back() * zPvRoundingErrorAdjust : binEdgesZPv.back() / zPvRoundingErrorAdjust; // init correct size of zPv-mult matrix - savedTriggersZPvMult.resize(binsZPv.size() - 1); - for (size_t i_zPv = 0; i_zPv < binsZPv.size() - 1; i_zPv++) { - savedTriggersZPvMult[i_zPv].resize(binsMult.size()); + savedTriggersZPvMult.resize(nBinsZPv); + for (size_t i_zPv = 0; i_zPv < nBinsZPv; i_zPv++) { + savedTriggersZPvMult[i_zPv].resize(nBinsMult); } } // save trigger for mixing - // up to nTriggerSavedForMixing stored (LIFO) + // up to nTriggerPerBinLimit stored in rolling buffer for each zPv-mult bin void saveTrigger(float const pt, float const phi, float const eta, double const zPv, double const mult) { - int const iBinCorrZPv = findIntervalBin(zPv, triggerBinValuesZPv); - int const iBinCorrMult = findIntervalBin(mult, triggerBinValuesMult); - // special cases (floating point precision errors, mult overflow) should be taken care of by triggerBinValuesZPv and triggerBinValuesMult - savedTriggersZPvMult[iBinCorrZPv][iBinCorrMult].push_front(MixingTrigger{pt, phi, eta}); - if (static_cast(savedTriggersZPvMult[iBinCorrZPv][iBinCorrMult].size()) > nTriggerSavedForMixing) { - savedTriggersZPvMult[iBinCorrZPv][iBinCorrMult].pop_back(); + int const iBinCorrZPv = getZPvBin(zPv); + int const iBinCorrMult = getMultBin(mult); + // fill + std::deque& triggerBin = savedTriggersZPvMult[iBinCorrZPv][iBinCorrMult]; + triggerBin.push_front(MixingTrigger{pt, phi, eta}); + if (triggerBin.size() > nTriggerPerBinLimit) { + triggerBin.pop_back(); + return; // skip nTriggerBinMin update when unnecessary + } + // update nTriggerBinMin + nTriggerBinMin = nTriggerPerBinLimit; + for (size_t i_zPv = 0; i_zPv < nBinsZPv; i_zPv++) { + for (size_t i_mult = 0; i_mult < nBinsMult; i_mult++) { + nTriggerBinMin = std::min(nTriggerBinMin, savedTriggersZPvMult[i_zPv][i_mult].size()); + } } } // return deques of trigger pt, phi, eta in the given zPv/mult bin std::deque const& getTriggers(double const zPv, double const mult) const { - int const iBinCorrZPv = findIntervalBin(zPv, triggerBinValuesZPv); - int const iBinCorrMult = findIntervalBin(mult, triggerBinValuesMult); + int const iBinCorrZPv = getZPvBin(zPv); + int const iBinCorrMult = getMultBin(mult); return savedTriggersZPvMult[iBinCorrZPv][iBinCorrMult]; } + size_t getBinSizeMin() const + { + return nTriggerBinMin; + } + private: double const zPvRoundingErrorAdjust = 1.0001; - int nTriggerSavedForMixing; - std::vector triggerBinValuesZPv; - std::vector triggerBinValuesMult; + std::vector binEdgesZPv; + size_t const nBinsZPv; + std::vector binEdgesMult; + size_t const nBinsMult; std::vector>> savedTriggersZPvMult; + size_t const nTriggerPerBinLimit; + size_t nTriggerBinMin = 0; + + size_t getZPvBin(double const zPv) const + { + int const iBinInit = findIntervalBin(zPv, binEdgesZPv); + if (iBinInit == -1 || iBinInit == static_cast(nBinsZPv)) { + throw std::runtime_error("zPv underflow or overflow in MixingTriggerMemory"); + } + return iBinInit; + } + size_t getMultBin(double const mult) const + { + int const iBinInit = findIntervalBin(mult, binEdgesMult); + if (iBinInit == static_cast(nBinsMult)) { + return iBinInit - 1; + } + if (iBinInit == -1) { + throw std::runtime_error("mult underflow in MixingTriggerMemory"); + } + return iBinInit; + } }; - MixingTriggerMemory mixingTriggerMemoryReco{nTriggerSavedForMixing.value, binsZPv.value, binsMult.value}; - MixingTriggerMemory mixingTriggerMemoryTrue{nTriggerSavedForMixing.value, binsZPvMcTrue.value, binsMult.value}; - MixingTriggerMemory mixingTriggerMemoryRecoColTrue{nTriggerSavedForMixing.value, binsZPvMcTrue.value, binsMult.value}; + // different buffers for mixing modes + enum class TriggerMemoryMode { Reco, + True, + TrueAssocEv, + TrueAssocEvRecoPtTrig }; + + std::unique_ptr triggerMemoryReco; + std::unique_ptr triggerMemoryTrue; + std::unique_ptr triggerMemoryTrueAssocEv; + std::unique_ptr triggerMemoryTrueAssocEvRecoPtTrig; - // functions ================================================================================================================================================================================ + // functions //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ========================================================================================================================================================================================== // general (kenobi) ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ // get histograms from ccdb - // save efficiencies from ccdb in histogram registry + // save corrections from ccdb in histogram registry void initCcdbHistograms() { // trigger - h1PtInvEffTrigger = nullptr; - if (doEffCorrectionTrigger) { - h1PtInvEffTrigger = ccdb->getForTimeStamp(pathCcdbEff.value + "/trigger", noLaterThanCcdb); - - const double* effBinsTrigger = h1PtInvEffTrigger->GetXaxis()->GetXbins()->GetArray(); - const AxisSpec axisPtEffTrigger{std::vector(effBinsTrigger, effBinsTrigger + h1PtInvEffTrigger->GetNbinsX() + 1), "#it{p}_{T}"}; - histos.add("usedEff/h1_pt_invEff_trigger_ccdb", "h1_pt_invEff_trigger_ccdb", kTH1D, {axisPtEffTrigger}, true); - for (int iBin = 1; iBin <= h1PtInvEffTrigger->GetNbinsX(); iBin++) { - histos.get(HIST("usedEff/h1_pt_invEff_trigger_ccdb"))->SetBinContent(iBin, h1PtInvEffTrigger->GetBinContent(iBin)); - histos.get(HIST("usedEff/h1_pt_invEff_trigger_ccdb"))->SetBinError(iBin, h1PtInvEffTrigger->GetBinError(iBin)); - } - } - // hadron - h1PtInvEffHadron = nullptr; - if (doEffCorrectionHadron) { - h1PtInvEffHadron = ccdb->getForTimeStamp(pathCcdbEff.value + "/hadron", noLaterThanCcdb); - - const double* effBinsHadron = h1PtInvEffHadron->GetXaxis()->GetXbins()->GetArray(); - const AxisSpec axisPtEffHadron{std::vector(effBinsHadron, effBinsHadron + h1PtInvEffHadron->GetNbinsX() + 1), "#it{p}_{T}"}; - histos.add("usedEff/h1_pt_invEff_hadron_ccdb", "h1_pt_invEff_hadron_ccdb", kTH1D, {axisPtEffHadron}, true); - for (int iBin = 1; iBin <= h1PtInvEffHadron->GetNbinsX(); iBin++) { - histos.get(HIST("usedEff/h1_pt_invEff_hadron_ccdb"))->SetBinContent(iBin, h1PtInvEffHadron->GetBinContent(iBin)); - histos.get(HIST("usedEff/h1_pt_invEff_hadron_ccdb"))->SetBinError(iBin, h1PtInvEffHadron->GetBinError(iBin)); - } - } - // pipm - h1PtInvEffPipm = nullptr; - if (doEffCorrectionPipm) { - h1PtInvEffPipm = ccdb->getForTimeStamp(pathCcdbEff.value + "/pipm", noLaterThanCcdb); - - const double* effBinsPipm = h1PtInvEffPipm->GetXaxis()->GetXbins()->GetArray(); - const AxisSpec axisPtEffPipm{std::vector(effBinsPipm, effBinsPipm + h1PtInvEffPipm->GetNbinsX() + 1), "#it{p}_{T}"}; - histos.add("usedEff/h1_pt_invEff_pipm_ccdb", "h1_pt_invEff_pipm_ccdb", kTH1D, {axisPtEffPipm}, true); - for (int iBin = 1; iBin <= h1PtInvEffPipm->GetNbinsX(); iBin++) { - histos.get(HIST("usedEff/h1_pt_invEff_pipm_ccdb"))->SetBinContent(iBin, h1PtInvEffPipm->GetBinContent(iBin)); - histos.get(HIST("usedEff/h1_pt_invEff_pipm_ccdb"))->SetBinError(iBin, h1PtInvEffPipm->GetBinError(iBin)); - } - } - // photonPCM - h1PtInvEffPhotonPCM = nullptr; - if (doEffCorrectionPhotonPCM) { - h1PtInvEffPhotonPCM = ccdb->getForTimeStamp(pathCcdbEff.value + "/photonPCM", noLaterThanCcdb); - - const double* effBinsPhotonPCM = h1PtInvEffPhotonPCM->GetXaxis()->GetXbins()->GetArray(); - const AxisSpec axisPtEffPhotonPCM{std::vector(effBinsPhotonPCM, effBinsPhotonPCM + h1PtInvEffPhotonPCM->GetNbinsX() + 1), "#it{p}_{T}"}; - histos.add("usedEff/h1_pt_invEff_photonPCM_ccdb", "h1_pt_invEff_photonPCM_ccdb", kTH1D, {axisPtEffPhotonPCM}, true); - for (int iBin = 1; iBin <= h1PtInvEffPhotonPCM->GetNbinsX(); iBin++) { - histos.get(HIST("usedEff/h1_pt_invEff_photonPCM_ccdb"))->SetBinContent(iBin, h1PtInvEffPhotonPCM->GetBinContent(iBin)); - histos.get(HIST("usedEff/h1_pt_invEff_photonPCM_ccdb"))->SetBinError(iBin, h1PtInvEffPhotonPCM->GetBinError(iBin)); + h1PtCorrectionTrigger = nullptr; + if (applyCTotTrigger) { + h1PtCorrectionTrigger = ccdb->getForTimeStamp(pathCcdbCorrection.value + "/trigger", noLaterThanCcdb); + + const double* correctionBinsTrigger = h1PtCorrectionTrigger->GetXaxis()->GetXbins()->GetArray(); + const AxisSpec axisPtCorrectionTrigger{std::vector(correctionBinsTrigger, correctionBinsTrigger + h1PtCorrectionTrigger->GetNbinsX() + 1), "#it{p}_{T}"}; + histos.add("usedCorrection/h1_pt_correction_trigger_ccdb", "h1_pt_correction_trigger_ccdb", kTH1D, {axisPtCorrectionTrigger}, true); + for (int iBin = 1; iBin <= h1PtCorrectionTrigger->GetNbinsX(); iBin++) { + histos.get(HIST("usedCorrection/h1_pt_correction_trigger_ccdb"))->SetBinContent(iBin, h1PtCorrectionTrigger->GetBinContent(iBin)); + histos.get(HIST("usedCorrection/h1_pt_correction_trigger_ccdb"))->SetBinError(iBin, h1PtCorrectionTrigger->GetBinError(iBin)); } } } @@ -355,43 +377,36 @@ struct PhotonChargedTriggerCorrelation { const AxisSpec axisN{1, 0., 1., "#it{N}_{something}"}; const AxisSpec axisCategories{16, 0., 16., "categories"}; - const AxisSpec axisZPv{nBinsZPv, -10, 10, "#it{z}_{pv}"}; - const AxisSpec axisZPvSmol{nBinsZPvSmol, -7, 7, "#it{z}_{pv}"}; + const AxisSpec axisZPv{nBinsZPv, -7, 7, "#it{z}_{pv}"}; const AxisSpec axisMult{nBinsMult + 1, -0.5, nBinsMult + 0.5, "multiplicity"}; - const AxisSpec axisMultSmol{nBinsMultSmol + 1, -0.5, nBinsMultSmol + 0.5, "multiplicity"}; const AxisSpec axisOccupancy{nBinsOccupancy + 1, -0.5, nBinsOccupancy + 0.5, "occupancy"}; const AxisSpec axisPhi{nBinsPhi, 0, constants::math::TwoPI, "#it{#varphi}"}; - const AxisSpec axisEta{nBinsEta, -etaMax, etaMax, "#it{#eta}"}; + const AxisSpec axisEta{nBinsEta, -absEtaMax, absEtaMax, "#it{#eta}"}; + const AxisSpec axisDCAz{nBinsDCAz, -5, 5, "DCA_{z}"}; const AxisSpec axisMgg{nBinsMgg, 0, 0.8, "#it{m}_{#gamma#gamma}"}; const AxisSpec axisPtTrig{binsPtTrig, "#it{p}_{T}^{trig}"}; const AxisSpec axisPtAssoc{binsPtAssoc, "#it{p}_{T}^{assoc}"}; const AxisSpec axisDPhi{binsDPhi, "#Delta#it{#varphi}"}; const AxisSpec axisDEta{binsDEta, "#Delta#it{#eta}"}; - const AxisSpec axisZPvBinning{binsZPv, "#it{z}_{pv} correlation binning"}; - const AxisSpec axisMultBinning{binsMult, "multiplicity correlation binning"}; - const AxisSpec axisZPvBinningMcTrue{binsZPvMcTrue, "#it{z}_{pv} correlation binning for mc true"}; - const AxisSpec axisMultBinningMcTrue{binsMultMcTrue, "multiplicity correlation binning for mc true"}; + const AxisSpec axisZPvBinning{binsZPvBinning, "#it{z}_{pv} correlation binning"}; + const AxisSpec axisMultBinning{binsMultBinning, "multiplicity correlation binning"}; + const AxisSpec axisZPvBinningMcTrue{binsZPvBinningMcTrue, "#it{z}_{pv} correlation binning for mc true"}; + const AxisSpec axisMultBinningMcTrue{binsMultBinningMcTrue, "multiplicity correlation binning for mc true"}; // reco info - histos.add("reco/info/h1_nEvents", "h1_nEvents", kTH1D, {axisCategories}); - histos.get(HIST("reco/info/h1_nEvents"))->GetXaxis()->SetBinLabel(1, "#it{N}_{ev}^{sel}"); - histos.get(HIST("reco/info/h1_nEvents"))->GetXaxis()->SetBinLabel(2, "#it{N}_{ev}"); - histos.get(HIST("reco/info/h1_nEvents"))->GetXaxis()->SetBinLabel(3, "#it{N}_{ev}^{trig}"); - - histos.add("reco/info/h2_zPvMult", "h2_zPvMult", kTHnSparseD, {axisZPv, axisMult}, true); - histos.add("reco/info/h1_occupancy", "h1_occupancy", kTH1D, {axisOccupancy}, true); - histos.add("reco/info/h2_zPvMult_trigEv", "h2_zPvMult_trigEv", kTHnSparseD, {axisZPv, axisMult}, true); - histos.add("reco/info/h1_occupancy_trigEv", "h1_occupancy_trigEv", kTH1D, {axisOccupancy}, true); + histos.add("reco/info/h1_nEvents", "h1_nEvents", kTH1D, {axisCategories}, true); + histos.add("reco/info/h3_ptTrigZPvMult", "h3_ptTrigZPvMult", kTHnSparseD, {axisPtTrig, axisZPv, axisMult}, true); + histos.add("reco/info/h2_ptTrigOccupancy", "h2_ptTrigOccupancy", kTHnSparseD, {axisPtTrig, axisOccupancy}, true); // reco (correlation) analysis histos.add("reco/corr/h3_ptPhiEta_trig", "h3_ptPhiEta_trig", kTHnSparseD, {axisPtAssoc, axisPhi, axisEta}, true); - std::function add_corrHists = - [&](std::string const name_id) { - histos.add(std::format("reco/corr/h3_ptPhiEta_assoc_{}", name_id).data(), std::format("h3_ptPhiEta_assoc_{}", name_id).data(), - kTHnSparseD, {axisPtAssoc, axisPhi, axisEta}, true); + auto const add_corrHists = + [&](std::string const& name_id) { + histos.add(std::format("reco/corr/h4_ptTrigPtAssocPhiEta_assoc_{}", name_id).data(), std::format("h4_ptTrigPtAssocPhiEta_assoc_{}", name_id).data(), + kTHnSparseD, {axisPtTrig, axisPtAssoc, axisPhi, axisEta}, true); histos.add(std::format("reco/corr/h6_corr_{}", name_id).data(), std::format("h6_corr_{}", name_id).data(), kTHnSparseF, {axisDPhi, axisDEta, axisPtTrig, axisPtAssoc, axisZPvBinning, axisMultBinning}, true); histos.add(std::format("reco/corr/h6_mix_{}", name_id).data(), std::format("h6_mix_{}", name_id).data(), @@ -407,9 +422,10 @@ struct PhotonChargedTriggerCorrelation { histos.add("reco/plain/h3_ptPhiEta_photonPCM", "h3_ptPhiEta_photonPCM", kTHnSparseD, {axisPtAssoc, axisPhi, axisEta}, true); add_corrHists("photonPCM"); // photonPCM pairs - histos.add("reco/plain/h4_ptMggZPvMult_photonPCMPair", "h4_ptMggZPvMult_photonPCMPair", kTHnSparseD, {axisPtAssoc, axisMgg, axisZPvBinning, axisMultBinning}, true); - histos.add("reco/plain/h4_ptMggZPvMult_trigEv_photonPCMPair", "h4_ptMggZPvMult_trigEv_photonPCMPair", kTHnSparseD, {axisPtAssoc, axisMgg, axisZPvBinning, axisMultBinning}, true); - histos.add("reco/corr/h5_ptTrigPtAssocMggZPvMult_assoc_photonPCMPair", "h5_ptTrigPtAssocMggZPvMult_assoc_photonPCMPair", kTHnSparseD, {axisPtTrig, axisPtAssoc, axisMgg, axisZPvBinning, axisMultBinning}, true); + histos.add("reco/plain/h5_ptTrigPtAssocMggZPvMult_photonPCMPair", "h5_ptTrigPtAssocMggZPvMult_photonPCMPair", + kTHnSparseD, {axisPtTrig, axisPtAssoc, axisMgg, axisZPvBinning, axisMultBinning}, true); + histos.add("reco/corr/h5_ptTrigPtAssocMggZPvMult_assoc_photonPCMPair", "h5_ptTrigPtAssocMggZPvMult_assoc_photonPCMPair", + kTHnSparseD, {axisPtTrig, axisPtAssoc, axisMgg, axisZPvBinning, axisMultBinning}, true); // pi0PCM add_corrHists("pi0PCMPeak"); add_corrHists("pi0PCMSide"); @@ -420,60 +436,153 @@ struct PhotonChargedTriggerCorrelation { // event mixing for photon pairs histos.add("reco/plain/h2_zPvMult_photonPCMPair_evMix", "h2_zPvMult_photonPCMPair_evMix", kTHnSparseD, {axisZPv, axisMult}, true); histos.add("reco/plain/h4_ptMggZPvMult_photonPCMPair_evMix", "h4_ptMggZPvMult_photonPCMPair_evMix", kTHnSparseD, {axisPtAssoc, axisMgg, axisZPvBinning, axisMultBinning}, true); - histos.add("reco/plain/h4_ptMggZPvMult_trigEv_photonPCMPair_evMix", "h4_ptMggZPvMult_trigEv_photonPCMPair_evMix", kTHnSparseD, {axisPtAssoc, axisMgg, axisZPvBinning, axisMultBinning}, true); // mc info - histos.add("mc/info/h1_nEvents_mcTrue", "h1_nEvents_mcTrue", kTH1D, {axisN}); - histos.add("mc/info/h2_zPvMult_mcTrue", "h2_zPvMult_mcTrue", kTHnSparseD, {axisZPv, axisMult}, true); - histos.add("mc/info/h1_nTrigEv_mcTrue", "h1_nTrigEv_mcTrue", kTH1D, {axisN}); - histos.add("mc/info/h2_zPvMult_trigEv_mcTrue", "h2_zPvMult_trigEv_mcTrue", kTHnSparseD, {axisZPv, axisMult}, true); - histos.add("mc/info/h1_nRecoCol_mcTrue", "h1_nRecoCol_mcTrue", kTH1D, {axisN}); - histos.add("mc/info/h2_zPvMult_recoCol_mcTrue", "h2_zPvMult_recoCol_mcTrue", kTHnSparseD, {axisZPv, axisMult}, true); + histos.add("mc/info/h3_ptTrigZPvMult_true", "h3_ptTrigZPvMult_true", kTHnSparseD, {axisPtTrig, axisZPv, axisMult}, true); + histos.add("mc/info/h3_ptTrigZPvMult_trueAssocEv", "h3_ptTrigZPvMult_trueAssocEv", kTHnSparseD, {axisPtTrig, axisZPv, axisMult}, true); + histos.add("mc/info/h3_ptTrigZPvMult_trueAssocEvRecoPtTrig", "h3_ptTrigZPvMult_trueAssocEvRecoPtTrig", kTHnSparseD, {axisPtTrig, axisZPv, axisMult}, true); // reco and true collision correlations - const std::vector assocMcCorrHistNames = {"hadron", "pipm", "photon", "pi0", "eta"}; - for (auto const& collision_type : {"true", "recoCol_true"}) { - histos.add(std::format("mc/{}/corr/h3_ptPhiEta_trig", collision_type).data(), "h3_ptPhiEta_trig", kTHnSparseD, {axisPtAssoc, axisPhi, axisEta}, true); - for (auto const& assocName : assocMcCorrHistNames) { - histos.add(std::format("mc/{}/corr/h6_corr_{}", collision_type, assocName).data(), std::format("h6_corr_{}", assocName).data(), + const std::vector assocMcCorrNamesMcAll = {"hadron", "pipm", "photon", "pi0", "eta"}; + for (auto const& correlationType : {"true", "trueAssocEv", "trueAssocEvRecoPtTrig"}) { + histos.add(std::format("mc/corr/h3_ptPhiEta_trig_{}", correlationType).data(), std::format("h3_ptPhiEta_trig_{}", correlationType).data(), + kTHnSparseD, {axisPtAssoc, axisPhi, axisEta}, true); + for (auto const& assocName : assocMcCorrNamesMcAll) { + histos.add(std::format("mc/corr/h6_corr_{}_{}", correlationType, assocName).data(), std::format("h6_corr_{}_{}", correlationType, assocName).data(), kTHnSparseD, {axisDPhi, axisDEta, axisPtTrig, axisPtAssoc, axisZPvBinningMcTrue, axisMultBinningMcTrue}, true); - histos.add(std::format("mc/{}/corr/h6_mix_{}", collision_type, assocName).data(), std::format("h6_mix_{}", assocName).data(), + histos.add(std::format("mc/corr/h6_mix_{}_{}", correlationType, assocName).data(), std::format("h6_mix_{}_{}", correlationType, assocName).data(), kTHnSparseD, {axisDPhi, axisDEta, axisPtTrig, axisPtAssoc, axisZPvBinningMcTrue, axisMultBinningMcTrue}, true); } } - // mc efficiency/purity - std::function add_effHists = - [&](std::string const name_id) { - histos.add(std::format("mc/eff/h3_ptPhiEta_{}", name_id).data(), std::format("h3_ptPhiEta_{}", name_id).data(), - kTHnSparseD, {axisPtAssoc, axisPhi, axisEta}, true); - histos.add(std::format("mc/eff/h3_ptZPvMult_{}", name_id).data(), std::format("h3_ptZPvMult_{}", name_id).data(), - kTHnSparseD, {axisPtAssoc, axisZPvSmol, axisMultSmol}, true); + // decay correlation extra info (just true level) + const std::vector assocMcCorrNamesMcDecayAddition = {"photonDecay", "photonDirect", "photonPi0", "photonEta", + "omega", "photonOmega", "etaPrime", "photonEtaPrime", "photonOtherMother"}; + for (auto const& assocName : assocMcCorrNamesMcDecayAddition) { + histos.add(std::format("mc/corr/h6_corr_true_{}", assocName).data(), std::format("h6_corr_true_{}", assocName).data(), + kTHnSparseD, {axisDPhi, axisDEta, axisPtTrig, axisPtAssoc, axisZPvBinningMcTrue, axisMultBinningMcTrue}, true); + histos.add(std::format("mc/corr/h6_mix_true_{}", assocName).data(), std::format("h6_mix_true_{}", assocName).data(), + kTHnSparseD, {axisDPhi, axisDEta, axisPtTrig, axisPtAssoc, axisZPvBinningMcTrue, axisMultBinningMcTrue}, true); + } + + // extra reco correlations with MC info + auto const addMcRecoCorrHists = + [&](std::string const& correlationType, std::string const& assocName) { + histos.add(std::format("mc/corr/h6_corr_{}_{}", correlationType, assocName).data(), std::format("h6_corr_{}_{}", correlationType, assocName).data(), + kTHnSparseD, {axisDPhi, axisDEta, axisPtTrig, axisPtAssoc, axisZPvBinning, axisMultBinning}, true); + histos.add(std::format("mc/corr/h6_mix_{}_{}", correlationType, assocName).data(), std::format("h6_mix_{}_{}", correlationType, assocName).data(), + kTHnSparseD, {axisDPhi, axisDEta, axisPtTrig, axisPtAssoc, axisZPvBinning, axisMultBinning}, true); + }; + // matchable associated particles + const std::vector assocCorrNamesMcExtra = {"hadron", "pipm", "photonPCM"}; + for (auto const& correlationType : {"recoAssocEv", "recoPure", "recoPureTruePtAssoc"}) { + for (auto const& assocName : assocCorrNamesMcExtra) { + addMcRecoCorrHists(correlationType, assocName); + } + } + // invariant mass associated particles + // mc pseudo yield + addMcRecoCorrHists("pseudoReco", "pi0PCM"); + addMcRecoCorrHists("pseudoRecoAssocEv", "pi0PCM"); + addMcRecoCorrHists("pseudoRecoPure", "pi0PCM"); + addMcRecoCorrHists("pseudoRecoPureTruePtAssoc", "pi0PCM"); + // efficiency advanced info + addMcRecoCorrHists("trueGeoAcc", "pi0"); + addMcRecoCorrHists("trueMeasDecay", "pi0"); + + // mc correction + auto const addCorrectionHistsResolved = + [&](std::string const& nameId) { + histos.add(std::format("mc/correction/resol/h4_ptTrigPtAssocPhiEta_{}", nameId).data(), std::format("h4_ptTrigPtAssocPhiEta_{}", nameId).data(), + kTHnSparseD, {axisPtTrig, axisPtAssoc, axisPhi, axisEta}, true); + histos.add(std::format("mc/correction/resol/h4_ptTrigPtAssocZPvMult_{}", nameId).data(), std::format("h4_ptTrigPtAssocZPvMult_{}", nameId).data(), + kTHnSparseD, {axisPtTrig, axisPtAssoc, axisZPv, axisMult}, true); + }; + auto const addCorrectionHistsResolvedBasicTypes = + [&](std::string const& nameReco, std::string const& nameTrue) { + std::vector const recoTypes = {"reco", "recoAssocEv", "recoPure"}; + std::vector const trueTypes = {"true", "trueAssocEv", "trueAssocEvRecoPtTrig"}; + for (std::string const& recoType : recoTypes) { + addCorrectionHistsResolved(std::format("{}_{}", recoType, nameReco)); + } + for (std::string const& trueType : trueTypes) { + addCorrectionHistsResolved(std::format("{}_{}", trueType, nameTrue)); + } + }; + + auto const addCorrectionHistPt = + [&](std::string const& nameId) { + histos.add(std::format("mc/correction/h2_ptTrigPtAssoc_{}", nameId).data(), std::format("h2_ptTrigPtAssoc_{}", nameId).data(), + kTHnSparseD, {axisPtTrig, axisPtAssoc}, true); + }; + auto const addCorrectionHistPtCategories = + [&](std::string const& nameId) { + histos.add(std::format("mc/correction/h3_ptTrigPtAssocCategory_{}", nameId).data(), std::format("h3_ptTrigPtAssocCategory_{}", nameId).data(), + kTHnSparseD, {axisPtTrig, axisPtAssoc, axisCategories}, true); + }; + auto const addCorrectionHistsV0 = + [&](std::string const& nameReco, std::string const& nameTrue) { + addCorrectionHistPt(std::format("true_{}", nameTrue).data()); + addCorrectionHistPt(std::format("trueAssocEv_{}", nameTrue).data()); + addCorrectionHistPt(std::format("trueAssocEvRecoPtTrig_{}", nameTrue).data()); + // mc pseudo yield + addCorrectionHistPt(std::format("pseudoReco_{}", nameReco).data()); + addCorrectionHistPt(std::format("pseudoRecoAssocEv_{}", nameReco).data()); + addCorrectionHistPt(std::format("pseudoRecoPure_{}", nameReco).data()); + // efficiency advanced info + addCorrectionHistPt(std::format("trueGeoAcc_{}", nameTrue).data()); + addCorrectionHistPt(std::format("trueMeasDecay_{}", nameTrue).data()); }; - // mc tracks - add_effHists("mcReco_hadron"); - add_effHists("mcReco_hasCorrectMc_hadron"); - add_effHists("mcTrue_recoCol_hadron"); - // mc pipm PID - add_effHists("mcReco_pipm"); - add_effHists("mcReco_hasCorrectMc_pipm"); - add_effHists("mcTrue_recoCol_pipm"); - // mc photonPCM - add_effHists("mcReco_photonPCM"); - add_effHists("mcReco_hasCorrectMc_photonPCM"); - add_effHists("mcTrue_recoCol_photon"); - // mc pi0 - add_effHists("mcTrue_recoCol_pi0"); - // mc eta - add_effHists("mcTrue_recoCol_eta"); - - // test of the test while testing another test. featuring a test - histos.add("test/h2_mult_comp", "h2_mult_comp", kTH2D, {axisMult, axisMult}, true); - histos.add("test/h2_tracks_zPvMultDep", "h2_tracks_zPvMultDep", kTH2D, {axisZPv, axisMult}, true); - histos.add("test/h2_globalTracks_zPvMultDep", "h2_globalTracks_zPvMultDep", kTH2D, {axisZPv, axisMult}, true); + + // matchables + // tracks + addCorrectionHistsResolvedBasicTypes("hadron", "hadron"); + // pipm PID + addCorrectionHistsResolvedBasicTypes("pipm", "pipm"); + // photonPCM + addCorrectionHistsResolvedBasicTypes("photonPCM", "photon"); + // impurities + addCorrectionHistPtCategories("recoImpurities_photonPCM"); + + // invariant mass reconstruction + // pi0 + addCorrectionHistsV0("pi0PCM", "pi0"); + // impurities + addCorrectionHistPtCategories("pseudoRecoImpurities_pi0PCM"); + // eta + addCorrectionHistsV0("etaPCM", "eta"); + + // tests + + // DCAz + histos.add("reco/plain/h5_ptTrigPtAssocDCAzZPvMult_photonPCM", "h5_ptTrigPtAssocDCAzZPvMult_photonPCM", + kTHnSparseD, {axisPtTrig, axisPtAssoc, axisDCAz, axisZPvBinning, axisMultBinning}, true); + histos.add("mc/plain/h5_ptTrigPtAssocDCAzZPvMult_photonPCM", "h5_ptTrigPtAssocDCAzZPvMult_photonPCM", + kTHnSparseD, {axisPtTrig, axisPtAssoc, axisDCAz, axisZPvBinning, axisMultBinning}, true); } // selections /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ + + // mc split event selection based in the thrid decimal of mc-true posZ + template + bool checkSplitMcEventSelection(T_collision const& collision) + { + // select based on configurables + return collision.thirdDecimalTruePosZ() % splitMcEvN == splitMcEvSelect - 1; + } + + // total event selection (basic selections from producer and analysis level selections) + template + bool totalEvSel(T_collision const& collision) + { + bool isSelEv = true; + + if constexpr (requires { collision.selEv(); }) { + isSelEv = collision.selEv(); + } + return isSelEv && checkSplitMcEventSelection(collision); + } // checks if mcParticle is charged template @@ -484,52 +593,69 @@ struct PhotonChargedTriggerCorrelation { return false; return true; } - // checks if mcParticle should be detected (physicalPrimary, |eta|) + // checks if fast decaying mcParticle is 'primary' template - bool checkPrimaryEtaMc(T_mcParticle const& mcParticle) + bool checkDecayPrimary(T_mcParticle const& mcParticle) { - if (!mcParticle.isPhysicalPrimary()) - return false; - if (std::abs(mcParticle.eta()) > etaMax) - return false; - return true; + // identify decaying primary + return mcParticle.producedByGenerator(); } - // checks if mcParticle should be detected as primary track (physicalPrimary, charge, |eta|) + // checks if mcParticle daughters are two photons template - bool checkPrimaryTrackMc(T_mcParticle const& mcParticle) + bool checkToGG(T_mcParticle const& mcParticle) { - if (!checkPrimaryEtaMc(mcParticle)) - return false; - if (!checkChargedMc(mcParticle)) + auto const& daughters = mcParticle.template daughters_as(); + constexpr int NDaughtersToGG = 2; + if (daughters.size() != NDaughtersToGG) return false; + auto daughter = daughters.begin(); + auto daughterEnd = daughters.end(); + while (daughter != daughterEnd) { + if ((*daughter).pdgCode() != PDG_t::kGamma) + return false; + daughter++; + } return true; } - // checks if mcParticle should be detected as 'primary' (|eta| not checked) + // check if particle has mother in parent tree template - bool checkH0Primary(T_mcParticle const& mcParticle, int const pdg) + bool checkForMother(T_mcParticle mcParticle, int const pdgCode, bool const checkAntiParticle) { - if (mcParticle.pdgCode() != pdg) + if (!mcParticle.has_mothers()) return false; - const auto& h0Daughters = mcParticle.template daughters_as(); - // identify primary h0 (account for 0 daughters for some reason) - if (h0Daughters.size() == 0) - return false; - for (auto const& h0_daughter : h0Daughters) { - if (!h0_daughter.isPhysicalPrimary()) - return false; + auto const mothers = mcParticle.template mothers_as(); + auto mother = mothers.begin(); + auto motherEnd = mothers.end(); + while (mother != motherEnd) { + + // LOGF(info, "searchPdgCode: %i, current[ pdgCode: %i, status: %i, primary: %i ], mother[ pdgCode: %i, status: %i, primary: %i ]", + // pdgCode, + // mcParticle.pdgCode(), mcParticle.getGenStatusCode(), mcParticle.isPhysicalPrimary(), + // (*mother).pdgCode(), (*mother).getGenStatusCode(), (*mother).isPhysicalPrimary()); + + if ((*mother).pdgCode() == pdgCode || (checkAntiParticle && (*mother).pdgCode() == -pdgCode)) + return true; + if (checkForMother(*mother, pdgCode, checkAntiParticle)) + return true; + mother++; } - return true; + return false; } - // checks if mcParticle should be detected as 'primary' pi0->gammagamma (|eta| not checked) + + // checks if daughters are in acceptance template - bool checkH0ToGG(T_mcParticle const& mcParticle, int const pdg) + bool checkDaughtersInAcceptance(T_mcParticle const mcParticle, std::pair const etaRange, std::pair const phiRange = {0, constants::math::TwoPI}) { - if (!checkH0Primary(mcParticle, pdg)) - return false; - // select h0 -> gg - constexpr int NDaughtersH0ToGG = 2; - if (mcParticle.template daughters_as().size() != NDaughtersH0ToGG) - return false; + auto const& daughterPhotons = mcParticle.template daughters_as(); + auto daughterPhoton = daughterPhotons.begin(); + auto daughterPhotonEnd = daughterPhotons.end(); + while (daughterPhoton != daughterPhotonEnd) { + if ((*daughterPhoton).eta() < etaRange.first || (*daughterPhoton).eta() > etaRange.second) + return false; + if ((*daughterPhoton).phi() < phiRange.first || (*daughterPhoton).phi() > phiRange.second) + return false; + daughterPhoton++; + } return true; } @@ -550,51 +676,46 @@ struct PhotonChargedTriggerCorrelation { return true; }; - // checks if tracks come from pi0 double conversion + // checks if tracks come from double conversion template - bool isGGFromPi0(T_track const& posTrack1, T_track const& negTrack1, T_track const& posTrack2, T_track const& negTrack2) + bool isGGFromDoubleConversion(T_track const& posTrack1, T_track const& negTrack1, T_track const& posTrack2, T_track const& negTrack2, int const pdgCode) { if (!isConversionPhoton(posTrack1, negTrack1) || !isConversionPhoton(posTrack2, negTrack2)) return false; // check same mother auto const& mothers1 = (*(posTrack1.mcParticle().template mothers_as().begin())).template mothers_as(); auto const& mothers2 = (*(posTrack2.mcParticle().template mothers_as().begin())).template mothers_as(); - constexpr int NMothersPhotonFromPi0 = 2; // for some reason two mothers (same particle) for pi0 decays (contradicts PYTHIA documentation, but whatever) - if (mothers1.size() != NMothersPhotonFromPi0 || mothers2.size() != NMothersPhotonFromPi0) + constexpr int NMothersPhotonFromH0 = 2; // for some reason two mothers (same particle) for h0 decays (contradicts PYTHIA documentation, but whatever) + if (mothers1.size() != NMothersPhotonFromH0 || mothers2.size() != NMothersPhotonFromH0) return false; if (mothers1.begin()->globalIndex() != mothers2.begin()->globalIndex()) return false; - // check pi0 - if (mothers1.begin()->pdgCode() != PDG_t::kPi0) + // check particle type + if (mothers1.begin()->pdgCode() != pdgCode) return false; return true; }; // analysis helpers ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ template double getH1ValueAt(T_h1 const* const h1, double const value) { return h1->GetBinContent(h1->FindFixBin(value)); } - // efficiency helpers - enum class EffParticleType { Trigger, - Hadron, - Pipm, - PhotonPCM }; - // efficiency function - template - double getInvEff(double const value) + // correction helpers + enum class CorrectionParticleType { Trigger, + Hadron, + Pipm, + PhotonPCM }; + + template + double getPtCorrection(double const value) { - if constexpr (T_effParticleType == EffParticleType::Trigger) { - return doEffCorrectionTrigger ? getH1ValueAt(h1PtInvEffTrigger, value) : 1; - } else if constexpr (T_effParticleType == EffParticleType::Hadron) { - return doEffCorrectionHadron ? getH1ValueAt(h1PtInvEffHadron, value) : 1; - } else if constexpr (T_effParticleType == EffParticleType::Pipm) { - return doEffCorrectionPipm ? getH1ValueAt(h1PtInvEffPipm, value) : 1; - } else if constexpr (T_effParticleType == EffParticleType::PhotonPCM) { - return doEffCorrectionPhotonPCM ? getH1ValueAt(h1PtInvEffPhotonPCM, value) : 1; + if constexpr (T_correctionParticleType == CorrectionParticleType::Trigger) { + return applyCTotTrigger ? getH1ValueAt(h1PtCorrectionTrigger, value) : 1; } else { return 1; } @@ -641,21 +762,21 @@ struct PhotonChargedTriggerCorrelation { } // check if invariant mass range - enum class MassRange { pi0PCMPeak, - pi0PCMSide, - etaPCMPeak, - etaPCMSide }; + enum class MassRange { Pi0PCMPeak, + Pi0PCMSide, + EtaPCMPeak, + EtaPCMSide }; - template + template bool checkMassRange(double const mgg) { - if constexpr (T_massRange == MassRange::pi0PCMPeak) { + if constexpr (range == MassRange::Pi0PCMPeak) { return mgg > pi0PCMPeakMassRange.value[0] && mgg < pi0PCMPeakMassRange.value[1]; - } else if constexpr (T_massRange == MassRange::pi0PCMSide) { + } else if constexpr (range == MassRange::Pi0PCMSide) { return mgg > pi0PCMSideMassRange.value[0] && mgg < pi0PCMSideMassRange.value[1]; - } else if constexpr (T_massRange == MassRange::etaPCMPeak) { + } else if constexpr (range == MassRange::EtaPCMPeak) { return mgg > etaPCMPeakMassRange.value[0] && mgg < etaPCMPeakMassRange.value[1]; - } else if constexpr (T_massRange == MassRange::etaPCMSide) { + } else if constexpr (range == MassRange::EtaPCMSide) { return (mgg > etaPCMLowSideMassRange.value[0] && mgg < etaPCMLowSideMassRange.value[1]) || (mgg > etaPCMHighSideMassRange.value[0] && mgg < etaPCMHighSideMassRange.value[1]); } @@ -663,6 +784,7 @@ struct PhotonChargedTriggerCorrelation { } // analysis ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ // generalised correlation functions // per collision @@ -691,34 +813,67 @@ struct PhotonChargedTriggerCorrelation { } } + // mixing-number pt scaling with power law + size_t nMixingPtPowerScaling(double const pt, size_t const n0) + { + double const rawScale = 1 + nTriggerScaleCoefficient * std::pow(pt, nTriggerScaleExponent); + return n0 * static_cast(rawScale); + } + // mixing - template void corrProcessMixing(T_collision const& collision, T_associatedThisEvent const& associatedThisEvent, T_funcMixing&& funcMixing, - size_t const nTriggerMixing, size_t const nTriggersThisDataFrame) + size_t const nTriggerMixingAt0) { - // skip if event does not contain valid trigger - if (doTrigEvMixing && !collision.trigEv()) + size_t triggerMemoryBinSizeMin = 0; + + // mixing loops (more dynamic than O2 mixing) + std::deque const& savedTriggers = + [&]() -> std::deque const& { + if constexpr (triggerMemoryMode == TriggerMemoryMode::Reco) { + triggerMemoryBinSizeMin = triggerMemoryReco->getBinSizeMin(); + return triggerMemoryReco->getTriggers(collision.posZ(), collision.nGlobalTracks()); + } else if constexpr (triggerMemoryMode == TriggerMemoryMode::True) { + triggerMemoryBinSizeMin = triggerMemoryTrue->getBinSizeMin(); + return triggerMemoryTrue->getTriggers(collision.posZ(), collision.nChargedInEtaRange()); + } else if constexpr (triggerMemoryMode == TriggerMemoryMode::TrueAssocEv) { + triggerMemoryBinSizeMin = triggerMemoryTrueAssocEv->getBinSizeMin(); + return triggerMemoryTrueAssocEv->getTriggers(collision.posZ(), collision.nChargedInEtaRange()); + } else if constexpr (triggerMemoryMode == TriggerMemoryMode::TrueAssocEvRecoPtTrig) { + triggerMemoryBinSizeMin = triggerMemoryTrueAssocEvRecoPtTrig->getBinSizeMin(); + return triggerMemoryTrueAssocEvRecoPtTrig->getTriggers(collision.posZ(), collision.nChargedInEtaRange()); + } + }(); + // require enough mixing statistics + if (triggerMemoryBinSizeMin < static_cast(nTriggerBinMinThreshold)) { return; - - // mixing loops (more efficient than O2 mixing (for now)) - auto savedTriggers = mixingTriggerMemoryReco.getTriggers(collision.posZ(), collision.nGlobalTracks()); - // number of triggers - const size_t mixUpToTriggerN = std::min(savedTriggers.size(), nTriggerMixing + nTriggersThisDataFrame); - const float perTriggerWeight = 1. / (mixUpToTriggerN - nTriggersThisDataFrame); // mixUpToTriggerN <= nTriggersThisDataFrame not problematic since no loop then - // mixing loops - for (size_t i_mixingTrigger = nTriggersThisDataFrame; i_mixingTrigger < mixUpToTriggerN; i_mixingTrigger++) { - for (auto const& associated : associatedThisEvent) { - funcMixing(collision, savedTriggers[i_mixingTrigger].pt(), savedTriggers[i_mixingTrigger].phi(), savedTriggers[i_mixingTrigger].eta(), associated, perTriggerWeight); + } + const size_t nSavedTriggers = savedTriggers.size(); + const size_t iStartTriggerParticle = randomIntInInterval(0, iTriggerStartFraction * nSavedTriggers); + // associated loop + for (auto const& associated : associatedThisEvent) { + // trigger numbers + const size_t nTriggerMixing = nMixingPtPowerScaling(associated.pt(), nTriggerMixingAt0); + const size_t mixUpToTriggerN = std::min(nSavedTriggers, iStartTriggerParticle + nTriggerMixing); + const float perTriggerWeight = 1. / (mixUpToTriggerN - iStartTriggerParticle); // mixUpToTriggerN <= iStartTriggerParticle caught by break statement of loop + // trigger loop + for (size_t i_mixingTrigger = iStartTriggerParticle; i_mixingTrigger < mixUpToTriggerN; i_mixingTrigger++) { + MixingTrigger const& mixingTrigger = savedTriggers[i_mixingTrigger]; + funcMixing(collision, mixingTrigger, associated, perTriggerWeight); } } } + // analysis tasks /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // =/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/=/ + void init(InitContext& initContext) { // analysis info - ccdb->setURL(urlCcdb.value); + ccdb->setURL(urlCcdb); // enabling object caching (otherwise each call goes to CCDB server) ccdb->setCaching(true); // ccdb->setLocalObjectValidityChecking(); @@ -728,75 +883,29 @@ struct PhotonChargedTriggerCorrelation { // init analysis variables // get variables from other tasks - o2::common::core::getTaskOptionValue(initContext, "photon-charged-trigger-producer", "etaMax", etaMax, false); + o2::common::core::getTaskOptionValue(initContext, "photon-charged-trigger-producer", "absEtaMax", absEtaMax, false); // histograms from ccdb initCcdbHistograms(); // create analysis histograms initHistograms(); - } - - // reconstructed //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - void processInfo(CorrCollision const& collision) - { - // all events - histos.fill(HIST("reco/info/h1_nEvents"), 1.5); - - // event selection - if (!collision.selEv()) - return; - histos.fill(HIST("reco/info/h1_nEvents"), 0.5); - - histos.fill(HIST("reco/info/h2_zPvMult"), collision.posZ(), collision.nGlobalTracks()); - histos.fill(HIST("reco/info/h1_occupancy"), collision.trackOccupancyInTimeRange()); - - // trigger events - if (!collision.trigEv()) - return; - histos.fill(HIST("reco/info/h1_nEvents"), 2.5); - histos.fill(HIST("reco/info/h2_zPvMult_trigEv"), collision.posZ(), collision.nGlobalTracks()); - histos.fill(HIST("reco/info/h1_occupancy_trigEv"), collision.trackOccupancyInTimeRange()); + // mixing trigger memory + triggerMemoryReco = std::make_unique(binsZPvBinning, binsMultBinning, nTriggerSavedForMixing); + triggerMemoryTrue = std::make_unique(binsZPvBinningMcTrue, binsMultBinningMcTrue, nTriggerSavedForMixing); + triggerMemoryTrueAssocEv = std::make_unique(binsZPvBinningMcTrue, binsMultBinningMcTrue, nTriggerSavedForMixing); + triggerMemoryTrueAssocEvRecoPtTrig = std::make_unique(binsZPvBinningMcTrue, binsMultBinningMcTrue, nTriggerSavedForMixing); } - PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processInfo, "process general info on collisions and tracks for analysis and qa", false); - void processCorrFirst(CorrCollisions const& collisions, aod::Triggers const& triggers) - { - // do at beginning of each data frame (before other reco correlation process functions) - // (PROCESS_SWITCH of this process has to be declared first) - - // [wow, such empty] - - for (auto const& collision : collisions) { - // event selection - if (!collision.selEv()) - continue; - - // group collision - auto const triggersThisEvent = triggers.sliceBy(perColTriggers, collision.globalIndex()); - - // trigger loop - for (auto const& trigger : triggersThisEvent) { - // trigger info - histos.fill(HIST("reco/corr/h3_ptPhiEta_trig"), trigger.pt(), trigger.phi(), trigger.eta(), - getInvEff(trigger.pt())); - - // save trigger for mixing - mixingTriggerMemoryReco.saveTrigger(trigger.pt(), trigger.phi(), trigger.eta(), collision.posZ(), collision.nGlobalTracks()); - } - } - } - PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processCorrFirst, "process to gather info before correlation processes", false); + // reconstructed //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ========================================================================================================================================================================================== void processCorrHadron(CorrCollisions const& collisions, aod::Triggers const& triggers, aod::Hadrons const& hadrons) { - size_t const nTriggersThisDataFrame = triggers.size(); - for (auto const& collision : collisions) { // event selection - if (!collision.selEv()) + if (!totalEvSel(collision)) continue; // group collision @@ -805,8 +914,7 @@ struct PhotonChargedTriggerCorrelation { auto const funcPlain = [this]([[maybe_unused]] auto const& collision, auto const& associated) { histos.fill(HIST("reco/plain/h3_ptPhiEta_hadron"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(associated.pt())); + associated.pt(), associated.phi(), associated.eta()); }; corrProcessPlain(collision, hadronsThisEvent, funcPlain); @@ -815,37 +923,38 @@ struct PhotonChargedTriggerCorrelation { if (trigger.jetTrackId() == associated.jetTrackId()) return; - histos.fill(HIST("reco/corr/h3_ptPhiEta_assoc_hadron"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(trigger.pt()) * getInvEff(associated.pt())); + histos.fill(HIST("reco/corr/h4_ptTrigPtAssocPhiEta_assoc_hadron"), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta(), + getPtCorrection(trigger.pt())); histos.fill(HIST("reco/corr/h6_corr_hadron"), getDeltaPhi(trigger.phi(), associated.phi()), trigger.eta() - associated.eta(), trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), - getInvEff(trigger.pt()) * getInvEff(associated.pt())); + getPtCorrection(trigger.pt())); }; corrProcessCorrelation(collision, triggersThisEvent, hadronsThisEvent, funcCorrelation); - auto const funcMixing = [this](auto const& collision, - float const mixingTriggerPt, float const mixingTriggerPhi, float const mixingTriggerEta, auto const& associated, auto const perTriggerWeight) { + // select mixing events + if (doTrigEvMixing && !collision.trigEv()) + continue; + + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { histos.fill(HIST("reco/corr/h6_mix_hadron"), - getDeltaPhi(mixingTriggerPhi, associated.phi()), - mixingTriggerEta - associated.eta(), - mixingTriggerPt, associated.pt(), collision.posZ(), collision.nGlobalTracks(), - perTriggerWeight * getInvEff(mixingTriggerPt) * getInvEff(associated.pt())); + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); }; - corrProcessMixing(collision, hadronsThisEvent, funcMixing, nTriggerMixingHadron, nTriggersThisDataFrame); + corrProcessMixing(collision, hadronsThisEvent, funcMixing, nMixingAt0Hadron); } } PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processCorrHadron, "process standard correlation for associated hardons", false); void processCorrPipm(CorrCollisions const& collisions, aod::Triggers const& triggers, aod::Pipms const& pipms) { - size_t const nTriggersThisDataFrame = triggers.size(); - for (auto const& collision : collisions) { // event selection - if (!collision.selEv()) + if (!totalEvSel(collision)) continue; // group collision @@ -854,8 +963,7 @@ struct PhotonChargedTriggerCorrelation { auto const funcPlain = [this]([[maybe_unused]] auto const& collision, auto const& associated) { histos.fill(HIST("reco/plain/h3_ptPhiEta_pipm"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(associated.pt())); + associated.pt(), associated.phi(), associated.eta()); }; corrProcessPlain(collision, pipmsThisEvent, funcPlain); @@ -864,37 +972,38 @@ struct PhotonChargedTriggerCorrelation { if (trigger.jetTrackId() == associated.jetTrackId()) return; - histos.fill(HIST("reco/corr/h3_ptPhiEta_assoc_pipm"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(trigger.pt()) * getInvEff(associated.pt())); + histos.fill(HIST("reco/corr/h4_ptTrigPtAssocPhiEta_assoc_pipm"), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta(), + getPtCorrection(trigger.pt())); histos.fill(HIST("reco/corr/h6_corr_pipm"), getDeltaPhi(trigger.phi(), associated.phi()), trigger.eta() - associated.eta(), trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), - getInvEff(trigger.pt()) * getInvEff(associated.pt())); + getPtCorrection(trigger.pt())); }; corrProcessCorrelation(collision, triggersThisEvent, pipmsThisEvent, funcCorrelation); - auto const funcMixing = [this](auto const& collision, - float const mixingTriggerPt, float const mixingTriggerPhi, float const mixingTriggerEta, auto const& associated, auto const perTriggerWeight) { + // select mixing events + if (doTrigEvMixing && !collision.trigEv()) + continue; + + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { histos.fill(HIST("reco/corr/h6_mix_pipm"), - getDeltaPhi(mixingTriggerPhi, associated.phi()), - mixingTriggerEta - associated.eta(), - mixingTriggerPt, associated.pt(), collision.posZ(), collision.nGlobalTracks(), - perTriggerWeight * getInvEff(mixingTriggerPt) * getInvEff(associated.pt())); + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); }; - corrProcessMixing(collision, pipmsThisEvent, funcMixing, nTriggerMixingPipm, nTriggersThisDataFrame); + corrProcessMixing(collision, pipmsThisEvent, funcMixing, nMixingAt0Pipm); } } PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processCorrPipm, "process standard correlation for associated pipm", false); void processCorrPhotonPCM(CorrCollisions const& collisions, aod::Triggers const& triggers, aod::PhotonPCMs const& photonPCMs) { - size_t const nTriggersThisDataFrame = triggers.size(); - for (auto const& collision : collisions) { // event selection - if (!collision.selEv()) + if (!totalEvSel(collision)) continue; // group collision @@ -903,47 +1012,47 @@ struct PhotonChargedTriggerCorrelation { auto const funcPlain = [this]([[maybe_unused]] auto const& collision, auto const& associated) { histos.fill(HIST("reco/plain/h3_ptPhiEta_photonPCM"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(associated.pt())); + associated.pt(), associated.phi(), associated.eta()); }; corrProcessPlain(collision, photonPCMsThisEvent, funcPlain); auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { // exclude self correlation - if (trigger.jetTrackId() == associated.posTrackId() || trigger.jetTrackId() == associated.negTrackId()) + if (trigger.jetTrackId() == associated.posJetTrackId() || trigger.jetTrackId() == associated.negJetTrackId()) return; - histos.fill(HIST("reco/corr/h3_ptPhiEta_assoc_photonPCM"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(trigger.pt()) * getInvEff(associated.pt())); + histos.fill(HIST("reco/corr/h4_ptTrigPtAssocPhiEta_assoc_photonPCM"), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta(), + getPtCorrection(trigger.pt())); histos.fill(HIST("reco/corr/h6_corr_photonPCM"), getDeltaPhi(trigger.phi(), associated.phi()), trigger.eta() - associated.eta(), trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), - getInvEff(trigger.pt()) * getInvEff(associated.pt())); + getPtCorrection(trigger.pt())); }; corrProcessCorrelation(collision, triggersThisEvent, photonPCMsThisEvent, funcCorrelation); - auto const funcMixing = [this](auto const& collision, - float const mixingTriggerPt, float const mixingTriggerPhi, float const mixingTriggerEta, auto const& associated, auto const perTriggerWeight) { + // select mixing events + if (doTrigEvMixing && !collision.trigEv()) + continue; + + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { histos.fill(HIST("reco/corr/h6_mix_photonPCM"), - getDeltaPhi(mixingTriggerPhi, associated.phi()), - mixingTriggerEta - associated.eta(), - mixingTriggerPt, associated.pt(), collision.posZ(), collision.nGlobalTracks(), - perTriggerWeight * getInvEff(mixingTriggerPt) * getInvEff(associated.pt())); + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); }; - corrProcessMixing(collision, photonPCMsThisEvent, funcMixing, nTriggerMixingPhotonPCM, nTriggersThisDataFrame); + corrProcessMixing(collision, photonPCMsThisEvent, funcMixing, nMixingAt0PhotonPCM); } } PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processCorrPhotonPCM, "process standard correlation for associated photonPCM", false); void processCorrPhotonPCMPair(CorrCollisions const& collisions, aod::Triggers const& triggers, aod::PhotonPCMPairs const& photonPCMPairs) { - size_t const nTriggersThisDataFrame = triggers.size(); - for (auto const& collision : collisions) { // event selection - if (!collision.selEv()) + if (!totalEvSel(collision)) continue; // group collision @@ -951,113 +1060,110 @@ struct PhotonChargedTriggerCorrelation { auto const photonPCMPairsThisEvent = photonPCMPairs.sliceBy(perColPhotonPCMPairs, collision.globalIndex()); auto const funcPlain = [this](auto const& collision, auto const& associated) { - histos.fill(HIST("reco/plain/h4_ptMggZPvMult_photonPCMPair"), associated.pt(), associated.mgg(), collision.posZ(), collision.nGlobalTracks()); + histos.fill(HIST("reco/plain/h5_ptTrigPtAssocMggZPvMult_photonPCMPair"), collision.ptMax(), associated.pt(), associated.mgg(), collision.posZ(), collision.nGlobalTracks()); }; corrProcessPlain(collision, photonPCMPairsThisEvent, funcPlain); - auto const funcPlainTrigEv = [this](auto const& collision, auto const& associated) { - histos.fill(HIST("reco/plain/h4_ptMggZPvMult_trigEv_photonPCMPair"), associated.pt(), associated.mgg(), collision.posZ(), collision.nGlobalTracks()); - }; - if (collision.trigEv()) - corrProcessPlain(collision, photonPCMPairsThisEvent, funcPlainTrigEv); - auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { // exclude self correlation - if (trigger.jetTrackId() == associated.posTrack1Id() || trigger.jetTrackId() == associated.negTrack1Id() || - trigger.jetTrackId() == associated.negTrack2Id() || trigger.jetTrackId() == associated.posTrack2Id()) + if (trigger.jetTrackId() == associated.posJetTrack1Id() || trigger.jetTrackId() == associated.negJetTrack1Id() || + trigger.jetTrackId() == associated.negJetTrack2Id() || trigger.jetTrackId() == associated.posJetTrack2Id()) return; histos.fill(HIST("reco/corr/h5_ptTrigPtAssocMggZPvMult_assoc_photonPCMPair"), trigger.pt(), associated.pt(), associated.mgg(), collision.posZ(), collision.nGlobalTracks(), - getInvEff(trigger.pt())); + getPtCorrection(trigger.pt())); // pi0 - if (checkMassRange(associated.mgg())) { - histos.fill(HIST("reco/corr/h3_ptPhiEta_assoc_pi0PCMPeak"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(trigger.pt())); + if (checkMassRange(associated.mgg())) { + histos.fill(HIST("reco/corr/h4_ptTrigPtAssocPhiEta_assoc_pi0PCMPeak"), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta(), + getPtCorrection(trigger.pt())); histos.fill(HIST("reco/corr/h6_corr_pi0PCMPeak"), getDeltaPhi(trigger.phi(), associated.phi()), trigger.eta() - associated.eta(), trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), - getInvEff(trigger.pt())); + getPtCorrection(trigger.pt())); return; } - if (checkMassRange(associated.mgg())) { - histos.fill(HIST("reco/corr/h3_ptPhiEta_assoc_pi0PCMSide"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(trigger.pt())); + if (checkMassRange(associated.mgg())) { + histos.fill(HIST("reco/corr/h4_ptTrigPtAssocPhiEta_assoc_pi0PCMSide"), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta(), + getPtCorrection(trigger.pt())); histos.fill(HIST("reco/corr/h6_corr_pi0PCMSide"), getDeltaPhi(trigger.phi(), associated.phi()), trigger.eta() - associated.eta(), trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), - getInvEff(trigger.pt())); + getPtCorrection(trigger.pt())); return; } // eta - if (checkMassRange(associated.mgg())) { - histos.fill(HIST("reco/corr/h3_ptPhiEta_assoc_etaPCMPeak"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(trigger.pt())); + if (checkMassRange(associated.mgg())) { + histos.fill(HIST("reco/corr/h4_ptTrigPtAssocPhiEta_assoc_etaPCMPeak"), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta(), + getPtCorrection(trigger.pt())); histos.fill(HIST("reco/corr/h6_corr_etaPCMPeak"), getDeltaPhi(trigger.phi(), associated.phi()), trigger.eta() - associated.eta(), trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), - getInvEff(trigger.pt())); + getPtCorrection(trigger.pt())); return; } - if (checkMassRange(associated.mgg())) { - histos.fill(HIST("reco/corr/h3_ptPhiEta_assoc_etaPCMSide"), - associated.pt(), associated.phi(), associated.eta(), - getInvEff(trigger.pt())); + if (checkMassRange(associated.mgg())) { + histos.fill(HIST("reco/corr/h4_ptTrigPtAssocPhiEta_assoc_etaPCMSide"), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta(), + getPtCorrection(trigger.pt())); histos.fill(HIST("reco/corr/h6_corr_etaPCMSide"), getDeltaPhi(trigger.phi(), associated.phi()), trigger.eta() - associated.eta(), trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), - getInvEff(trigger.pt())); + getPtCorrection(trigger.pt())); return; } }; corrProcessCorrelation(collision, triggersThisEvent, photonPCMPairsThisEvent, funcCorrelation); - auto const funcMixing = [this](auto const& collision, - float const mixingTriggerPt, float const mixingTriggerPhi, float const mixingTriggerEta, auto const& associated, auto const perTriggerWeight) { + // select mixing events + if (doTrigEvMixing && !collision.trigEv()) + continue; + + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { // pi0 - if (checkMassRange(associated.mgg())) { + if (checkMassRange(associated.mgg())) { histos.fill(HIST("reco/corr/h6_mix_pi0PCMPeak"), - getDeltaPhi(mixingTriggerPhi, associated.phi()), - mixingTriggerEta - associated.eta(), - mixingTriggerPt, associated.pt(), collision.posZ(), collision.nGlobalTracks(), - perTriggerWeight * getInvEff(mixingTriggerPt)); + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); return; } - if (checkMassRange(associated.mgg())) { + if (checkMassRange(associated.mgg())) { histos.fill(HIST("reco/corr/h6_mix_pi0PCMSide"), - getDeltaPhi(mixingTriggerPhi, associated.phi()), - mixingTriggerEta - associated.eta(), - mixingTriggerPt, associated.pt(), collision.posZ(), collision.nGlobalTracks(), - perTriggerWeight * getInvEff(mixingTriggerPt)); + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); return; } // eta - if (checkMassRange(associated.mgg())) { + if (checkMassRange(associated.mgg())) { histos.fill(HIST("reco/corr/h6_mix_etaPCMPeak"), - getDeltaPhi(mixingTriggerPhi, associated.phi()), - mixingTriggerEta - associated.eta(), - mixingTriggerPt, associated.pt(), collision.posZ(), collision.nGlobalTracks(), - perTriggerWeight * getInvEff(mixingTriggerPt)); + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); return; } - if (checkMassRange(associated.mgg())) { + if (checkMassRange(associated.mgg())) { histos.fill(HIST("reco/corr/h6_mix_etaPCMSide"), - getDeltaPhi(mixingTriggerPhi, associated.phi()), - mixingTriggerEta - associated.eta(), - mixingTriggerPt, associated.pt(), collision.posZ(), collision.nGlobalTracks(), - perTriggerWeight * getInvEff(mixingTriggerPt)); + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); return; } }; - corrProcessMixing(collision, photonPCMPairsThisEvent, funcMixing, nTriggerMixingH0PCM, nTriggersThisDataFrame); + corrProcessMixing(collision, photonPCMPairsThisEvent, funcMixing, nMixingAt0H0PCM); } } PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processCorrPhotonPCMPair, "process standard correlation for associated pi0PCM", false); @@ -1072,15 +1178,15 @@ struct PhotonChargedTriggerCorrelation { auto const& [collision1, photonPCMs1, collision2, photonPCMs2] = *pair; // // check that current und mixing-trigger event are from the same zPv/mult bins - // if (checkSameBin(collision1.posZ(), collision2.posZ(), binsZPv) == -1) { + // if (checkSameBin(collision1.posZ(), collision2.posZ(), binsZPvBinning) == -1) { // std::printf("ERROR: zPv bins do not match\n"); continue; // } - // if (checkSameBin(collision1.nGlobalTracks(), collision2.nGlobalTracks(), binsMult) == -1) { + // if (checkSameBin(collision1.nGlobalTracks(), collision2.nGlobalTracks(), binsMultBinning) == -1) { // std::printf("ERROR: multiplicity bins do not match\n"); continue; // } // event selection - if (!collision1.selEv() || !collision2.selEv()) + if (!totalEvSel(collision1) || !totalEvSel(collision2)) continue; // event info histos.fill(HIST("reco/plain/h2_zPvMult_photonPCMPair_evMix"), collision1.posZ(), collision1.nGlobalTracks()); @@ -1092,393 +1198,1200 @@ struct PhotonChargedTriggerCorrelation { histos.fill(HIST("reco/plain/h4_ptMggZPvMult_photonPCMPair_evMix"), p4photonPCMPair.pt(), p4photonPCMPair.M(), collision1.posZ(), collision1.nGlobalTracks()); } - - // trigger events - if (!collision1.trigEv() || !collision2.trigEv()) - continue; - // mixing loop - for (auto const& [photonPCM1, photonPCM2] : soa::combinations(soa::CombinationsFullIndexPolicy(photonPCMs1, photonPCMs2))) { - ROOT::Math::PtEtaPhiMVector const p4photonPCM1(photonPCM1.pt(), photonPCM1.eta(), photonPCM1.phi(), 0.); - ROOT::Math::PtEtaPhiMVector const p4photonPCM2(photonPCM2.pt(), photonPCM2.eta(), photonPCM2.phi(), 0.); - ROOT::Math::PtEtaPhiMVector const p4photonPCMPair = p4photonPCM1 + p4photonPCM2; - - histos.fill(HIST("reco/plain/h4_ptMggZPvMult_trigEv_photonPCMPair_evMix"), p4photonPCMPair.pt(), p4photonPCMPair.M(), collision1.posZ(), collision1.nGlobalTracks()); - } } } PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processCorrPhotonPCMPairMix, "process gamma-gamma mixing for photonPCM", false); - // mc /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - void processMcInfo(CorrMcCollisions const& mcCollisions, CorrMcDCollisions const& collisions) + void processPCMDCAz(CorrCollision const& collision, aod::PhotonPCMs const& photonPCMs, aod::V0PhotonsKF const&) { - for (auto const& mcCollision : mcCollisions) { - // all events - histos.fill(HIST("mc/info/h1_nEvents_mcTrue"), 0.5); - histos.fill(HIST("mc/info/h2_zPvMult_mcTrue"), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); + // event selection + if (!totalEvSel(collision)) + return; - // trigger events - if (!mcCollision.trigEv()) - continue; - histos.fill(HIST("mc/info/h1_nTrigEv_mcTrue"), 0.5); - histos.fill(HIST("mc/info/h2_zPvMult_trigEv_mcTrue"), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); - } - for (auto const& collision : collisions) { - // event selection - if (!collision.selEv()) - continue; - histos.fill(HIST("mc/info/h1_nRecoCol_mcTrue"), 0.5); - histos.fill(HIST("mc/info/h2_zPvMult_recoCol_mcTrue"), collision.mcCollision_as().posZ(), collision.mcCollision_as().nChargedInEtaRange()); + for (auto const& photonPCM : photonPCMs) { + histos.fill(HIST("reco/plain/h5_ptTrigPtAssocDCAzZPvMult_photonPCM"), + collision.ptMax(), photonPCM.pt(), photonPCM.v0PhotonKF().dcaZtopv(), collision.posZ(), collision.nGlobalTracks()); } } - PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcInfo, "process general info on mc collisions and tracks for analysis and qa", false); + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processPCMDCAz, "process to test DCAz distribution of photons in trigger collisions", false); + + // MC /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ========================================================================================================================================================================================== + + // correlations ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ // (sad) attempt at reducing code duplication - enum class McCorrEventType : int { True = 0, - RecoColTrue = 1 }; + enum class McCorrTruthLevel : int { True = 0, + TrueAssocEv = 1, + TrueAssocEvRecoPtTrig = 2, + TrueMeasDecay = 3, + TrueGeoAcc = 4 }; enum class McCorrCorrelationType : int { Correlation = 0, Mixing = 1 }; - enum class McCorrAssociatedType : int { Hadron = 0, - Pipm = 1, - Photon = 2, - Pi0 = 3, - Eta = 4 }; - static constexpr const char* McHistPaths[2][2][5] = { - {{"mc/true/corr/h6_corr_hadron", "mc/true/corr/h6_corr_pipm", "mc/true/corr/h6_corr_photon", - "mc/true/corr/h6_corr_pi0", "mc/true/corr/h6_corr_eta"}, - {"mc/true/corr/h6_mix_hadron", "mc/true/corr/h6_mix_pipm", "mc/true/corr/h6_mix_photon", - "mc/true/corr/h6_mix_pi0", "mc/true/corr/h6_mix_eta"}}, - {{"mc/recoCol_true/corr/h6_corr_hadron", "mc/recoCol_true/corr/h6_corr_pipm", "mc/recoCol_true/corr/h6_corr_photon", - "mc/recoCol_true/corr/h6_corr_pi0", "mc/recoCol_true/corr/h6_corr_eta"}, - {"mc/recoCol_true/corr/h6_mix_hadron", "mc/recoCol_true/corr/h6_mix_pipm", "mc/recoCol_true/corr/h6_mix_photon", - "mc/recoCol_true/corr/h6_mix_pi0", "mc/recoCol_true/corr/h6_mix_eta"}}}; - static constexpr const char* getMcHistPath(McCorrEventType eventType, McCorrCorrelationType correlationType, McCorrAssociatedType associatedType) + enum class McCorrAssocType : int { Hadron = 0, + Pipm = 1, + Photon = 2, + PhotonDecay = 3, + PhotonDirect = 4, + Pi0 = 5, + PhotonPi0 = 6, + Eta = 7, + PhotonEta = 8, + Omega = 9, + PhotonOmega = 10, + EtaPrime = 11, + PhotonEtaPrime = 12, + PhotonOtherMother = 13 }; + static constexpr const char* McTrueCorrHistPaths[5][2][14] = { + {{"mc/corr/h6_corr_true_hadron", "mc/corr/h6_corr_true_pipm", + "mc/corr/h6_corr_true_photon", "mc/corr/h6_corr_true_photonDecay", "mc/corr/h6_corr_true_photonDirect", + "mc/corr/h6_corr_true_pi0", "mc/corr/h6_corr_true_photonPi0", "mc/corr/h6_corr_true_eta", "mc/corr/h6_corr_true_photonEta", + "mc/corr/h6_corr_true_omega", "mc/corr/h6_corr_true_photonOmega", "mc/corr/h6_corr_true_etaPrime", "mc/corr/h6_corr_true_photonEtaPrime", "mc/corr/h6_corr_true_photonOtherMother"}, + {"mc/corr/h6_mix_true_hadron", "mc/corr/h6_mix_true_pipm", + "mc/corr/h6_mix_true_photon", "mc/corr/h6_mix_true_photonDecay", "mc/corr/h6_mix_true_photonDirect", + "mc/corr/h6_mix_true_pi0", "mc/corr/h6_mix_true_photonPi0", "mc/corr/h6_mix_true_eta", "mc/corr/h6_mix_true_photonEta", + "mc/corr/h6_mix_true_omega", "mc/corr/h6_mix_true_photonOmega", "mc/corr/h6_mix_true_etaPrime", "mc/corr/h6_mix_true_photonEtaPrime", "mc/corr/h6_mix_true_photonOtherMother"}}, + {{"mc/corr/h6_corr_trueAssocEv_hadron", "mc/corr/h6_corr_trueAssocEv_pipm", + "mc/corr/h6_corr_trueAssocEv_photon", "", "", + "mc/corr/h6_corr_trueAssocEv_pi0", "", "mc/corr/h6_corr_trueAssocEv_eta", "", + "", "", "", "", ""}, + {"mc/corr/h6_mix_trueAssocEv_hadron", "mc/corr/h6_mix_trueAssocEv_pipm", + "mc/corr/h6_mix_trueAssocEv_photon", "", "", + "mc/corr/h6_mix_trueAssocEv_pi0", "", "mc/corr/h6_mix_trueAssocEv_eta", "", + "", "", "", "", ""}}, + {{"mc/corr/h6_corr_trueAssocEvRecoPtTrig_hadron", "mc/corr/h6_corr_trueAssocEvRecoPtTrig_pipm", + "mc/corr/h6_corr_trueAssocEvRecoPtTrig_photon", "", "", + "mc/corr/h6_corr_trueAssocEvRecoPtTrig_pi0", "", "mc/corr/h6_corr_trueAssocEvRecoPtTrig_eta", "", + "", "", "", "", ""}, + {"mc/corr/h6_mix_trueAssocEvRecoPtTrig_hadron", "mc/corr/h6_mix_trueAssocEvRecoPtTrig_pipm", + "mc/corr/h6_mix_trueAssocEvRecoPtTrig_photon", "", "", + "mc/corr/h6_mix_trueAssocEvRecoPtTrig_pi0", "", "mc/corr/h6_mix_trueAssocEvRecoPtTrig_eta", "", + "", "", "", "", ""}}, + {{"", "", + "", "", "", + "mc/corr/h6_corr_trueMeasDecay_pi0", "", "", "", + "", "", "", "", ""}, + {"", "", + "", "", "", + "mc/corr/h6_mix_trueMeasDecay_pi0", "", "", "", + "", "", "", "", ""}}, + {{"", "", + "", "", "", + "mc/corr/h6_corr_trueGeoAcc_pi0", "", "", "", + "", "", "", "", ""}, + {"", "", + "", "", "", + "mc/corr/h6_mix_trueGeoAcc_pi0", "", "", "", + "", "", "", "", ""}}}; + static constexpr const char* getMcTrueCorrHistPath(McCorrTruthLevel truthLevel, McCorrCorrelationType correlationType, McCorrAssocType assocType) { - return McHistPaths[static_cast(eventType)][static_cast(correlationType)][static_cast(associatedType)]; + return McTrueCorrHistPaths[static_cast(truthLevel)][static_cast(correlationType)][static_cast(assocType)]; } - // fill mc correaltion histograms based on given associated mc particle - template - void fillMcCorrHists(auto const& mcCollision, auto const& trigger, auto const& associated, double const weight) + // fill mc correlation histograms based on given associated mc particle + template + void fillMcCorrHists(auto const& mcCollision, auto const& triggerIn, auto const& associated, double const weight) { + // check mc + if constexpr (requires { triggerIn.template jetTrack_as(); }) { + if (!triggerIn.template jetTrack_as().has_mcParticle()) + return; + } + // trigger info separation + auto const triggerParticle = [&]() { + if constexpr (correlationType == McCorrCorrelationType::Correlation) { + if constexpr (truthLevel == McCorrTruthLevel::TrueAssocEvRecoPtTrig) { + static_assert(requires { triggerIn.template jetTrack_as(); }, "triggerIn must be reco level for TrueAssocEvRecoPtTrig"); + return triggerIn.template jetTrack_as().mcParticle(); + } else if constexpr (truthLevel == McCorrTruthLevel::TrueAssocEv || truthLevel == McCorrTruthLevel::True) { + static_assert(requires { triggerIn.template jetMcParticle_as(); }, "triggerIn must be true level for TrueAssocEv and True"); + return triggerIn; + } + } else if constexpr (correlationType == McCorrCorrelationType::Mixing) { + static_assert(std::same_as, MixingTrigger>, "triggerIn must be MixingTrigger for mixing"); + return triggerIn; + } + }(); + auto const& triggerForPt = triggerIn; + + // exclude self correlation + if constexpr (correlationType == McCorrCorrelationType::Correlation) { + if constexpr (truthLevel == McCorrTruthLevel::TrueAssocEvRecoPtTrig) { + if (triggerParticle.globalIndex() == associated.globalIndex()) + return; + } else if constexpr (truthLevel == McCorrTruthLevel::TrueAssocEv || truthLevel == McCorrTruthLevel::True) { + if (triggerParticle.jetMcParticleId() == associated.globalIndex()) + return; + } + } + + double const dPhi = getDeltaPhi(triggerParticle.phi(), associated.phi()); + double const dEta = triggerParticle.eta() - associated.eta(); + double const ptTrig = triggerForPt.pt(); + double const ptAssoc = associated.pt(); + double const posZ = mcCollision.posZ(); + double const mult = mcCollision.nChargedInEtaRange(); + + if (std::abs(associated.eta()) > absEtaMax) + return; + // standard particles (marked physical primary) - if (checkPrimaryEtaMc(associated)) { - // charged primary ('hadron') selection + if (associated.isPhysicalPrimary()) { + // charged primary ('hadron') if (checkChargedMc(associated)) { - histos.fill(HIST(getMcHistPath(eventType, correlationType, McCorrAssociatedType::Hadron)), - getDeltaPhi(trigger.phi(), associated.phi()), - trigger.eta() - associated.eta(), - trigger.pt(), associated.pt(), mcCollision.posZ(), mcCollision.nChargedInEtaRange(), - weight); + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::Hadron)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); } - // pipm selection + // pipm if (std::abs(associated.pdgCode()) == PDG_t::kPiPlus) { - histos.fill(HIST(getMcHistPath(eventType, correlationType, McCorrAssociatedType::Pipm)), - getDeltaPhi(trigger.phi(), associated.phi()), - trigger.eta() - associated.eta(), - trigger.pt(), associated.pt(), mcCollision.posZ(), mcCollision.nChargedInEtaRange(), - weight); + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::Pipm)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); return; } - // photon selection + // photon if (associated.pdgCode() == PDG_t::kGamma) { - histos.fill(HIST(getMcHistPath(eventType, correlationType, McCorrAssociatedType::Photon)), - getDeltaPhi(trigger.phi(), associated.phi()), - trigger.eta() - associated.eta(), - trigger.pt(), associated.pt(), mcCollision.posZ(), mcCollision.nChargedInEtaRange(), - weight); + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::Photon)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + + // extra info for decay correlation only for total true level + if constexpr (truthLevel != McCorrTruthLevel::True) { + return; + } + + // decay and direct + int const statusDecayLow = 91; + int const statusDecayUp = 99; + int const statusDirect = 62; + + if (associated.getGenStatusCode() >= statusDecayLow && associated.getGenStatusCode() <= statusDecayUp) { + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::PhotonDecay)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, + weight); + // decays from different mothers + int const pdgMother = associated.template mothers_as().begin()->pdgCode(); + switch (pdgMother) { + case PDG_t::kPi0: + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::PhotonPi0)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + break; + case constants::physics::Pdg::kEta: + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::PhotonEta)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + break; + case constants::physics::Pdg::kOmega: + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::PhotonOmega)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + break; + case constants::physics::Pdg::kEtaPrime: + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::PhotonEtaPrime)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + break; + default: + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::PhotonOtherMother)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + break; + } + } else { + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::PhotonDirect)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, + weight); + + if (associated.getGenStatusCode() != statusDirect) { + LOGF(info, "filled primary photon with status: %i", associated.getGenStatusCode()); + } + } return; } return; } // decaying particles (not marked physical primary) - if ((std::abs(associated.eta()) < etaMax)) { - // pi0 selection - if (checkH0Primary(associated, PDG_t::kPi0)) { - histos.fill(HIST(getMcHistPath(eventType, correlationType, McCorrAssociatedType::Pi0)), - getDeltaPhi(trigger.phi(), associated.phi()), - trigger.eta() - associated.eta(), - trigger.pt(), associated.pt(), mcCollision.posZ(), mcCollision.nChargedInEtaRange(), - weight); + if (!checkDecayPrimary(associated)) + return; + // pi0 + if (associated.pdgCode() == PDG_t::kPi0) { + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::Pi0)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + // efficincy extra info + if constexpr (truthLevel != McCorrTruthLevel::TrueAssocEvRecoPtTrig) return; - } - // eta selection - if (checkH0Primary(associated, 221)) { - histos.fill(HIST(getMcHistPath(eventType, correlationType, McCorrAssociatedType::Eta)), - getDeltaPhi(trigger.phi(), associated.phi()), - trigger.eta() - associated.eta(), - trigger.pt(), associated.pt(), mcCollision.posZ(), mcCollision.nChargedInEtaRange(), - weight); + // chosen decay + if (!checkToGG(associated)) return; - } + histos.fill(HIST(getMcTrueCorrHistPath(McCorrTruthLevel::TrueMeasDecay, correlationType, McCorrAssocType::Pi0)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + // daughters in acceptance + if (!checkDaughtersInAcceptance(associated, {-absEtaMax, absEtaMax})) + return; + histos.fill(HIST(getMcTrueCorrHistPath(McCorrTruthLevel::TrueGeoAcc, correlationType, McCorrAssocType::Pi0)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + return; + } + // eta + if (associated.pdgCode() == constants::physics::Pdg::kEta) { + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::Eta)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + return; + } + // extra info for decay correlation only for total true level + if constexpr (truthLevel != McCorrTruthLevel::True) { + return; + } + // omega + if (associated.pdgCode() == constants::physics::Pdg::kOmega) { + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::Omega)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + return; + } + // etaPrime + if (associated.pdgCode() == constants::physics::Pdg::kEtaPrime) { + histos.fill(HIST(getMcTrueCorrHistPath(truthLevel, correlationType, McCorrAssocType::EtaPrime)), + dPhi, dEta, ptTrig, ptAssoc, posZ, mult, weight); + return; } } void processMcTrueCorr(CorrMcCollisions const& mcCollisions, aod::TriggerParticles const& triggerParticles, aod::JetParticles const& mcParticles) { for (auto const& mcCollision : mcCollisions) { + // event selection + if (!totalEvSel(mcCollision)) + continue; + // group collision auto const triggerParticlesThisEvent = triggerParticles.sliceBy(perColTriggerParticles, mcCollision.globalIndex()); auto const mcParticlesThisEvent = mcParticles.sliceBy(perColMcParticles, mcCollision.globalIndex()); - // trigger pairing loop - for (auto const& trigger : triggerParticlesThisEvent) { - // trigger info - histos.fill(HIST("mc/true/corr/h3_ptPhiEta_trig"), trigger.pt(), trigger.phi(), trigger.eta()); - - // save trigger for mixing - mixingTriggerMemoryTrue.saveTrigger(trigger.pt(), trigger.phi(), trigger.eta(), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); + auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { + fillMcCorrHists(collision, trigger, associated, 1); + }; + corrProcessCorrelation(mcCollision, triggerParticlesThisEvent, mcParticlesThisEvent, funcCorrelation); - for (auto const& associated : mcParticlesThisEvent) { - // exclude self correlation - if (trigger.jetMcParticleId() == associated.globalIndex()) - continue; + // select mixing events + if (doTrigEvMixing && !mcCollision.trigEv()) + continue; - fillMcCorrHists(mcCollision, trigger, associated, 1); - } - } + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { + fillMcCorrHists(collision, trigger, associated, perTriggerWeight); + }; + corrProcessMixing(mcCollision, mcParticlesThisEvent, funcMixing, nMixingAt0McTrue); } } PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcTrueCorr, "process mc-true (all collisions) correlation for multiple associated particles", false); - void processMcTrueMix(CorrMcCollisions const& mcCollisions, aod::TriggerParticles const& triggerParticles, aod::JetParticles const& mcParticles) + void processMcTrueAssocEvCorr(CorrMcDCollisions const& collisions, CorrMcCollisions const&, aod::TriggerParticles const& triggerParticles, aod::JetParticles const& mcParticles) { - for (auto const& mcCollision : mcCollisions) { + for (auto const& collision : collisions) { + // event selection + if (!totalEvSel(collision)) + continue; + + auto const& mcCollision = collision.mcCollision_as(); + // group collision - auto const triggerParticlesThisEvent = triggerParticles.sliceBy(perColTriggerParticles, mcCollision.globalIndex()); - auto const mcParticlesThisEvent = mcParticles.sliceBy(perColMcParticles, mcCollision.globalIndex()); + auto const triggerParticlesThisEvent = triggerParticles.sliceBy(perColTriggerParticles, collision.mcCollisionId()); + auto const mcParticlesThisEvent = mcParticles.sliceBy(perColMcParticles, collision.mcCollisionId()); - const size_t nTriggerParticlesThisDataFrame = triggerParticles.size(); - auto savedTriggers = mixingTriggerMemoryTrue.getTriggers(mcCollision.posZ(), mcCollision.nChargedInEtaRange()); - const size_t mixUpToTriggerN = std::min(savedTriggers.size(), static_cast(nTriggerMixingMcTrue) + nTriggerParticlesThisDataFrame); - const float perTriggerWeight = 1. / (mixUpToTriggerN - nTriggerParticlesThisDataFrame); + auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { + fillMcCorrHists(collision, trigger, associated, 1); + }; + corrProcessCorrelation(mcCollision, triggerParticlesThisEvent, mcParticlesThisEvent, funcCorrelation); - // trigger loop - for (size_t i_mixingTrigger = nTriggerParticlesThisDataFrame; i_mixingTrigger < mixUpToTriggerN; i_mixingTrigger++) { - MixingTrigger const& mixingTrigger = savedTriggers[i_mixingTrigger]; - for (auto const& associated : mcParticlesThisEvent) { - fillMcCorrHists(mcCollision, mixingTrigger, associated, perTriggerWeight); - } - } + // select mixing events + if (doTrigEvMixing && !mcCollision.trigEv()) + continue; + + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { + fillMcCorrHists(collision, trigger, associated, perTriggerWeight); + }; + corrProcessMixing(mcCollision, mcParticlesThisEvent, funcMixing, nMixingAt0McTrue); } } - PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcTrueMix, "process mc-true (all collisions) correlation mixing for multiple associated particles", false); + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcTrueAssocEvCorr, "process mc-true (reco collisions) correlation for multiple associated particles", false); - void processMcRecoColTrueCorr(CorrMcDCollisions const& collisions, CorrMcCollisions const&, aod::TriggerParticles const& triggerParticles, aod::JetParticles const& mcParticles) + void processMcTrueAssocEvRecoPtTrigCorr(CorrMcDCollisions const& collisions, aod::Triggers const& triggers, aod::JetTracksMCD const&, CorrMcCollisions const&, aod::JetParticles const& mcParticles) { for (auto const& collision : collisions) { // event selection - if (!collision.selEv()) + if (!totalEvSel(collision)) continue; - // group collision - auto const triggerParticlesThisEvent = triggerParticles.sliceBy(perColTriggerParticles, collision.mcCollisionId()); + auto const& mcCollision = collision.mcCollision_as(); + + // group collision + auto const triggersThisEvent = triggers.sliceBy(perColTriggers, collision.globalIndex()); auto const mcParticlesThisEvent = mcParticles.sliceBy(perColMcParticles, collision.mcCollisionId()); - auto const& mcCollision = collision.mcCollision_as(); + auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { + fillMcCorrHists(collision, trigger, associated, 1); + }; + corrProcessCorrelation(mcCollision, triggersThisEvent, mcParticlesThisEvent, funcCorrelation); - // trigger pairing loop - for (auto const& trigger : triggerParticlesThisEvent) { - // trigger info - histos.fill(HIST("mc/recoCol_true/corr/h3_ptPhiEta_trig"), trigger.pt(), trigger.phi(), trigger.eta()); + // select mixing events + if (doTrigEvMixing && !mcCollision.trigEv()) + continue; - // save trigger for mixing - mixingTriggerMemoryRecoColTrue.saveTrigger(trigger.pt(), trigger.phi(), trigger.eta(), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { + fillMcCorrHists(collision, trigger, associated, perTriggerWeight); + }; + corrProcessMixing(mcCollision, mcParticlesThisEvent, funcMixing, nMixingAt0McTrue); + } + } + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcTrueAssocEvRecoPtTrigCorr, "process mc-true (reco collisions with reco ptTrig) correlation for multiple associated particles", false); - // hadrons (tracks) and pipm - for (auto const& associated : mcParticlesThisEvent) { - // exclude self correlation - if (trigger.jetMcParticleId() == associated.globalIndex()) - continue; + void processMcRecoCorrHadron(CorrMcDCollisions const& collisions, aod::Triggers const& triggers, aod::Hadrons const& hadrons, + aod::JetTracksMCD const&, aod::JetParticles const&) + { + for (auto const& collision : collisions) { + // event selection + if (!totalEvSel(collision)) + continue; - fillMcCorrHists(mcCollision, trigger, associated, 1); - } - } + // group collision + auto const triggersThisEvent = triggers.sliceBy(perColTriggers, collision.globalIndex()); + auto const hadronsThisEvent = hadrons.sliceBy(perColHadrons, collision.globalIndex()); + + auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { + // exclude self correlation + if (trigger.jetTrackId() == associated.jetTrackId()) + return; + + // check mc + if (!associated.template jetTrack_as().has_mcParticle()) + return; + auto const& associatedMcParticle = associated.template jetTrack_as().mcParticle(); + // collision association + if (associatedMcParticle.mcCollisionId() != collision.mcCollisionId()) + return; + histos.fill(HIST("mc/corr/h6_corr_recoAssocEv_hadron"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // purity + if (!checkChargedMc(associatedMcParticle) || !associatedMcParticle.isPhysicalPrimary()) + return; + histos.fill(HIST("mc/corr/h6_corr_recoPure_hadron"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // true pt + histos.fill(HIST("mc/corr/h6_corr_recoPureTruePtAssoc_hadron"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associatedMcParticle.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + }; + corrProcessCorrelation(collision, triggersThisEvent, hadronsThisEvent, funcCorrelation); + + // select mixing events + if (doTrigEvMixing && !collision.trigEv()) + continue; + + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { + // check mc + if (!associated.template jetTrack_as().has_mcParticle()) + return; + auto const& associatedMcParticle = associated.template jetTrack_as().mcParticle(); + // collision association + if (associatedMcParticle.mcCollisionId() != collision.mcCollisionId()) + return; + histos.fill(HIST("mc/corr/h6_mix_recoAssocEv_hadron"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // purity + if (!checkChargedMc(associatedMcParticle) || !associatedMcParticle.isPhysicalPrimary()) + return; + histos.fill(HIST("mc/corr/h6_mix_recoPure_hadron"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // true pt + histos.fill(HIST("mc/corr/h6_mix_recoPureTruePtAssoc_hadron"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associatedMcParticle.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + }; + corrProcessMixing(collision, hadronsThisEvent, funcMixing, nMixingAt0Hadron); } } - PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcRecoColTrueCorr, "process mc-true (reco collisions) correlation for multiple associated particles", false); + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcRecoCorrHadron, "process correlation for associated hardons with additional mc information", false); - void processMcRecoColTrueMix(CorrMcDCollisions const& collisions, CorrMcCollisions const&, aod::TriggerParticles const& triggerParticles, aod::JetParticles const& mcParticles) + void processMcRecoCorrPipm(CorrMcDCollisions const& collisions, aod::Triggers const& triggers, aod::Pipms const& pipms, + aod::JetTracksMCD const&, aod::JetParticles const&) { for (auto const& collision : collisions) { // event selection - if (!collision.selEv()) + if (!totalEvSel(collision)) continue; // group collision - auto const triggerParticlesThisEvent = triggerParticles.sliceBy(perColTriggerParticles, collision.mcCollisionId()); - auto const mcParticlesThisEvent = mcParticles.sliceBy(perColMcParticles, collision.mcCollisionId()); + auto const triggersThisEvent = triggers.sliceBy(perColTriggers, collision.globalIndex()); + auto const pipmsThisEvent = pipms.sliceBy(perColPipms, collision.globalIndex()); - auto const& mcCollision = collision.mcCollision_as(); + auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { + // exclude self correlation + if (trigger.jetTrackId() == associated.jetTrackId()) + return; + + // check mc + if (!associated.template jetTrack_as().has_mcParticle()) + return; + auto const& associatedMcParticle = associated.template jetTrack_as().mcParticle(); + // collision association + if (associatedMcParticle.mcCollisionId() != collision.mcCollisionId()) + return; + histos.fill(HIST("mc/corr/h6_corr_recoAssocEv_pipm"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // purity + if (std::abs(associatedMcParticle.pdgCode()) != PDG_t::kPiPlus || !associatedMcParticle.isPhysicalPrimary()) + return; + histos.fill(HIST("mc/corr/h6_corr_recoPure_pipm"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // true pt + histos.fill(HIST("mc/corr/h6_corr_recoPureTruePtAssoc_pipm"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associatedMcParticle.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + }; + corrProcessCorrelation(collision, triggersThisEvent, pipmsThisEvent, funcCorrelation); - const size_t nTriggerParticlesThisDataFrame = triggerParticles.size(); - auto savedTriggers = mixingTriggerMemoryRecoColTrue.getTriggers(mcCollision.posZ(), mcCollision.nChargedInEtaRange()); - const size_t mixUpToTriggerN = std::min(savedTriggers.size(), static_cast(nTriggerMixingMcTrue) + nTriggerParticlesThisDataFrame); - const float perTriggerWeight = 1. / (mixUpToTriggerN - nTriggerParticlesThisDataFrame); + // select mixing events + if (doTrigEvMixing && !collision.trigEv()) + continue; - // trigger loop - for (size_t i_mixingTrigger = nTriggerParticlesThisDataFrame; i_mixingTrigger < mixUpToTriggerN; i_mixingTrigger++) { - MixingTrigger const& mixingTrigger = savedTriggers[i_mixingTrigger]; - for (auto const& associated : mcParticlesThisEvent) { - fillMcCorrHists(mcCollision, mixingTrigger, associated, perTriggerWeight); - } - } + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { + // check mc + if (!associated.template jetTrack_as().has_mcParticle()) + return; + auto const& associatedMcParticle = associated.template jetTrack_as().mcParticle(); + // collision association + if (associatedMcParticle.mcCollisionId() != collision.mcCollisionId()) + return; + histos.fill(HIST("mc/corr/h6_mix_recoAssocEv_pipm"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // purity + if (std::abs(associatedMcParticle.pdgCode()) != PDG_t::kPiPlus || !associatedMcParticle.isPhysicalPrimary()) + return; + histos.fill(HIST("mc/corr/h6_mix_recoPure_pipm"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // true pt + histos.fill(HIST("mc/corr/h6_mix_recoPureTruePtAssoc_pipm"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associatedMcParticle.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + }; + corrProcessMixing(collision, pipmsThisEvent, funcMixing, nMixingAt0Pipm); + } + } + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcRecoCorrPipm, "process correlation for associated pipms with additional mc information", false); + + void processMcRecoCorrPhotonPCM(CorrMcDCollisions const& collisions, aod::Triggers const& triggers, aod::PhotonPCMs const& photonPCMs, + aod::JetTracksMCD const&, aod::JetParticles const&) + { + for (auto const& collision : collisions) { + // event selection + if (!totalEvSel(collision)) + continue; + + // group collision + auto const triggersThisEvent = triggers.sliceBy(perColTriggers, collision.globalIndex()); + auto const photonPCMsThisEvent = photonPCMs.sliceBy(perColPhotonPCMs, collision.globalIndex()); + + auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { + // exclude self correlation + if (trigger.jetTrackId() == associated.posJetTrackId() || trigger.jetTrackId() == associated.negJetTrackId()) + return; + + // check mc + auto const& posTrack = associated.template posJetTrack_as(); + auto const& negTrack = associated.template negJetTrack_as(); + if (!posTrack.has_mcParticle() || !negTrack.has_mcParticle()) + return; + // collision association + if (posTrack.mcParticle().mcCollisionId() != collision.mcCollisionId() || negTrack.mcParticle().mcCollisionId() != collision.mcCollisionId()) + return; + histos.fill(HIST("mc/corr/h6_corr_recoAssocEv_photonPCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // purity + auto const& photons = posTrack.mcParticle().template mothers_as(); + if (!isConversionPhoton(posTrack, negTrack) || !photons.begin()->isPhysicalPrimary()) + return; + histos.fill(HIST("mc/corr/h6_corr_recoPure_photonPCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // true pt + histos.fill(HIST("mc/corr/h6_corr_recoPureTruePtAssoc_photonPCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), photons.begin()->pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + }; + corrProcessCorrelation(collision, triggersThisEvent, photonPCMsThisEvent, funcCorrelation); + + // select mixing events + if (doTrigEvMixing && !collision.trigEv()) + continue; + + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { + // check mc + auto const& posTrack = associated.template posJetTrack_as(); + auto const& negTrack = associated.template negJetTrack_as(); + if (!posTrack.has_mcParticle() || !negTrack.has_mcParticle()) + return; + // collision association + if (posTrack.mcParticle().mcCollisionId() != collision.mcCollisionId() || negTrack.mcParticle().mcCollisionId() != collision.mcCollisionId()) + return; + histos.fill(HIST("mc/corr/h6_mix_recoAssocEv_photonPCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // purity + auto const& photons = posTrack.mcParticle().template mothers_as(); + if (!isConversionPhoton(posTrack, negTrack) || !photons.begin()->isPhysicalPrimary()) + return; + histos.fill(HIST("mc/corr/h6_mix_recoPure_photonPCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // true pt + histos.fill(HIST("mc/corr/h6_mix_recoPureTruePtAssoc_photonPCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), photons.begin()->pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + }; + corrProcessMixing(collision, photonPCMsThisEvent, funcMixing, nMixingAt0PhotonPCM); } } - PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcRecoColTrueMix, "process mc-true (reco collisions) correlation mixing for multiple associated particles", false); + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcRecoCorrPhotonPCM, "process correlation for associated photonPCMs with additional mc information", false); - void processMcRecoColEff(CorrMcDCollision const& collision, aod::JetTracksMCD const& tracks, - aod::Hadrons const& hadrons, aod::Pipms const& pipms, aod::PhotonPCMs const& photonPCMs, - CorrMcCollisions const&, aod::JetParticles const& mcParticles, aod::TriggerParticles const& triggerParticles) + void processMcRecoCorrPi0PCM(CorrMcDCollisions const& collisions, aod::Triggers const& triggers, aod::PhotonPCMPairs const& photonPCMPairs, + aod::JetTracksMCD const&, aod::JetParticles const&) + { + for (auto const& collision : collisions) { + // event selection + if (!totalEvSel(collision)) + continue; + + // group collision + auto const triggersThisEvent = triggers.sliceBy(perColTriggers, collision.globalIndex()); + auto const photonPCMPairsThisEvent = photonPCMPairs.sliceBy(perColPhotonPCMPairs, collision.globalIndex()); + + auto const funcCorrelation = [this](auto const& collision, auto const& trigger, auto const& associated) { + // exclude self correlation + if (trigger.jetTrackId() == associated.posJetTrack1Id() || trigger.jetTrackId() == associated.negJetTrack1Id() || + trigger.jetTrackId() == associated.posJetTrack2Id() || trigger.jetTrackId() == associated.negJetTrack2Id()) + return; + + // check mc + auto const& posTrack1 = associated.template posJetTrack1_as(); + auto const& negTrack1 = associated.template negJetTrack1_as(); + auto const& posTrack2 = associated.template posJetTrack2_as(); + auto const& negTrack2 = associated.template negJetTrack2_as(); + if (!posTrack1.has_mcParticle() || !negTrack1.has_mcParticle() || !posTrack2.has_mcParticle() || !negTrack2.has_mcParticle()) + return; + // pseudo yield + if (!isGGFromDoubleConversion(posTrack1, negTrack1, posTrack2, negTrack2, PDG_t::kPi0)) + return; + histos.fill(HIST("mc/corr/h6_corr_pseudoReco_pi0PCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // collision association + if (posTrack1.mcParticle().mcCollisionId() != collision.mcCollisionId()) + return; + histos.fill(HIST("mc/corr/h6_corr_pseudoRecoAssocEv_pi0PCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // purity (just secondaries) + auto const& photons1 = posTrack1.mcParticle().template mothers_as(); + auto const& mothersOfPhoton = photons1.begin()->template mothers_as(); + if (!checkDecayPrimary(*(mothersOfPhoton.begin()))) + return; + histos.fill(HIST("mc/corr/h6_corr_pseudoRecoPure_pi0PCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + // true pt + histos.fill(HIST("mc/corr/h6_corr_pseudoRecoPureTruePtAssoc_pi0PCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), mothersOfPhoton.begin()->pt(), collision.posZ(), collision.nGlobalTracks(), + getPtCorrection(trigger.pt())); + }; + corrProcessCorrelation(collision, triggersThisEvent, photonPCMPairsThisEvent, funcCorrelation); + + // select mixing events + if (doTrigEvMixing && !collision.trigEv()) + continue; + + auto const funcMixing = [this](auto const& collision, auto const& trigger, auto const& associated, auto const perTriggerWeight) { + // check mc + auto const& posTrack1 = associated.template posJetTrack1_as(); + auto const& negTrack1 = associated.template negJetTrack1_as(); + auto const& posTrack2 = associated.template posJetTrack2_as(); + auto const& negTrack2 = associated.template negJetTrack2_as(); + if (!posTrack1.has_mcParticle() || !negTrack1.has_mcParticle() || !posTrack2.has_mcParticle() || !negTrack2.has_mcParticle()) + return; + // pseudo yield + if (!isGGFromDoubleConversion(posTrack1, negTrack1, posTrack2, negTrack2, PDG_t::kPi0)) + return; + histos.fill(HIST("mc/corr/h6_mix_pseudoReco_pi0PCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // collision association + if (posTrack1.mcParticle().mcCollisionId() != collision.mcCollisionId()) + return; + histos.fill(HIST("mc/corr/h6_mix_pseudoRecoAssocEv_pi0PCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // purity (just secondaries) + auto const& photons1 = posTrack1.mcParticle().template mothers_as(); + auto const& mothersOfPhoton = photons1.begin()->template mothers_as(); + if (!checkDecayPrimary(*(mothersOfPhoton.begin()))) + return; + histos.fill(HIST("mc/corr/h6_mix_pseudoRecoPure_pi0PCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), associated.pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + // true pt + histos.fill(HIST("mc/corr/h6_mix_pseudoRecoPureTruePtAssoc_pi0PCM"), + getDeltaPhi(trigger.phi(), associated.phi()), + trigger.eta() - associated.eta(), + trigger.pt(), mothersOfPhoton.begin()->pt(), collision.posZ(), collision.nGlobalTracks(), + perTriggerWeight * getPtCorrection(trigger.pt())); + }; + corrProcessMixing(collision, photonPCMPairsThisEvent, funcMixing, nMixingAt0PhotonPCM); + } + } + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcRecoCorrPi0PCM, "process correlation for associated pi0PCMs with additional mc information", false); + + // other //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ + + enum class McCorrectionObservable : int { PhiEta = 0, + ZPvMult = 1 }; + enum class McResolvedTruthLevel : int { Reco = 0, + RecoAssocEv = 1, + RecoPure = 2, + TrueAssocEvRecoPtTrig = 3, + TrueAssocEv = 4, + True = 5 }; + enum class CorrAssocType : int { Hadron = 0, + Pipm = 1, + PhotonPCM = 2 }; + static constexpr char const* McResolvedCorrectionHistPaths[2][6][3] = { + {{"mc/correction/resol/h4_ptTrigPtAssocPhiEta_reco_hadron", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_reco_pipm", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_reco_photonPCM"}, + {"mc/correction/resol/h4_ptTrigPtAssocPhiEta_recoAssocEv_hadron", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_recoAssocEv_pipm", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_recoAssocEv_photonPCM"}, + {"mc/correction/resol/h4_ptTrigPtAssocPhiEta_recoPure_hadron", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_recoPure_pipm", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_recoPure_photonPCM"}, + {"mc/correction/resol/h4_ptTrigPtAssocPhiEta_trueAssocEvRecoPtTrig_hadron", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_trueAssocEvRecoPtTrig_pipm", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_trueAssocEvRecoPtTrig_photon"}, + {"mc/correction/resol/h4_ptTrigPtAssocPhiEta_trueAssocEv_hadron", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_trueAssocEv_pipm", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_trueAssocEv_photon"}, + {"mc/correction/resol/h4_ptTrigPtAssocPhiEta_true_hadron", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_true_pipm", + "mc/correction/resol/h4_ptTrigPtAssocPhiEta_true_photon"}}, + {{"mc/correction/resol/h4_ptTrigPtAssocZPvMult_reco_hadron", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_reco_pipm", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_reco_photonPCM"}, + {"mc/correction/resol/h4_ptTrigPtAssocZPvMult_recoAssocEv_hadron", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_recoAssocEv_pipm", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_recoAssocEv_photonPCM"}, + {"mc/correction/resol/h4_ptTrigPtAssocZPvMult_recoPure_hadron", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_recoPure_pipm", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_recoPure_photonPCM"}, + {"mc/correction/resol/h4_ptTrigPtAssocZPvMult_trueAssocEvRecoPtTrig_hadron", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_trueAssocEvRecoPtTrig_pipm", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_trueAssocEvRecoPtTrig_photon"}, + {"mc/correction/resol/h4_ptTrigPtAssocZPvMult_trueAssocEv_hadron", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_trueAssocEv_pipm", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_trueAssocEv_photon"}, + {"mc/correction/resol/h4_ptTrigPtAssocZPvMult_true_hadron", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_true_pipm", + "mc/correction/resol/h4_ptTrigPtAssocZPvMult_true_photon"}}}; + static constexpr char const* getMcResolvedCorrectionHistPath(McCorrectionObservable observable, McResolvedTruthLevel truthLevel, CorrAssocType assocType) + { + return McResolvedCorrectionHistPaths[static_cast(observable)][static_cast(truthLevel)][static_cast(assocType)]; + } + template + void fillResolvedCorrectionHists(auto const& collision, auto const& associated) + { + if constexpr (truthLevel == McResolvedTruthLevel::Reco || truthLevel == McResolvedTruthLevel::RecoAssocEv || truthLevel == McResolvedTruthLevel::RecoPure) { + histos.fill(HIST(getMcResolvedCorrectionHistPath(McCorrectionObservable::PhiEta, truthLevel, assocType)), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta()); + histos.fill(HIST(getMcResolvedCorrectionHistPath(McCorrectionObservable::ZPvMult, truthLevel, assocType)), + collision.ptMax(), associated.pt(), collision.posZ(), collision.nGlobalTracks()); + } + if constexpr (truthLevel == McResolvedTruthLevel::TrueAssocEvRecoPtTrig) { + histos.fill(HIST(getMcResolvedCorrectionHistPath(McCorrectionObservable::PhiEta, truthLevel, assocType)), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta()); + histos.fill(HIST(getMcResolvedCorrectionHistPath(McCorrectionObservable::ZPvMult, truthLevel, assocType)), + collision.ptMax(), associated.pt(), collision.template mcCollision_as().posZ(), collision.template mcCollision_as().nChargedInEtaRange()); + } + if constexpr (truthLevel == McResolvedTruthLevel::True || truthLevel == McResolvedTruthLevel::TrueAssocEv) { + histos.fill(HIST(getMcResolvedCorrectionHistPath(McCorrectionObservable::PhiEta, truthLevel, assocType)), + collision.ptMax(), associated.pt(), associated.phi(), associated.eta()); + histos.fill(HIST(getMcResolvedCorrectionHistPath(McCorrectionObservable::ZPvMult, truthLevel, assocType)), + collision.ptMax(), associated.pt(), collision.posZ(), collision.nChargedInEtaRange()); + } + } + + enum class PhotonImpurity { WrongEv, + Misid, + K0Short, + K0Long, + Lambda0, + Other }; + static constexpr double binValuePi0Impurities(PhotonImpurity const impurity) + { + return static_cast(impurity) + 0.5; + } + + void processMcRecoCorrection(CorrMcDCollision const& collision, aod::JetTracksMCD const&, + aod::Hadrons const& hadrons, aod::Pipms const& pipms, aod::PhotonPCMs const& photonPCMs, aod::PhotonPCMPairs const& photonPCMPairs, + CorrMcCollisions const&, aod::JetParticles const& mcParticles) { // event selection - if (!collision.selEv()) + if (!totalEvSel(collision)) return; + auto const& mcCollision = collision.mcCollision_as(); + + // group collision auto const mcParticlesThisEvent = mcParticles.sliceBy(perColMcParticles, collision.mcCollisionId()); - auto const triggerParticlesThisEvent = triggerParticles.sliceBy(perColTriggerParticles, collision.mcCollisionId()); // hadrons for (auto const& hadron : hadrons) { - if (doTrigEvEff && !collision.trigEv() && hadron.pt() < ptCutTrigEvEff) - continue; - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcReco_hadron"), hadron.pt(), hadron.phi(), hadron.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcReco_hadron"), hadron.pt(), collision.posZ(), collision.nGlobalTracks()); - // purity + // reconstructed + fillResolvedCorrectionHists(collision, hadron); + // check mc if (!hadron.jetTrack_as().has_mcParticle()) continue; - auto const hadronParticle = hadron.jetTrack_as().mcParticle(); - if (!checkPrimaryTrackMc(hadronParticle)) + auto const& hadronParticle = hadron.jetTrack_as().mcParticle(); + // collision association + if (hadronParticle.mcCollisionId() != collision.mcCollisionId()) continue; - if (requireSingleCollisionPurity && hadronParticle.mcCollisionId() != collision.mcCollisionId()) + fillResolvedCorrectionHists(collision, hadron); + // purity + if (!checkChargedMc(hadronParticle) || !hadronParticle.isPhysicalPrimary()) continue; - - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcReco_hasCorrectMc_hadron"), hadron.pt(), hadron.phi(), hadron.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcReco_hasCorrectMc_hadron"), hadron.pt(), collision.posZ(), collision.nGlobalTracks()); + fillResolvedCorrectionHists(collision, hadron); } // pipm for (auto const& pipm : pipms) { - if (doTrigEvEff && !collision.trigEv() && pipm.pt() < ptCutTrigEvEff) - continue; - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcReco_pipm"), pipm.pt(), pipm.phi(), pipm.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcReco_pipm"), pipm.pt(), collision.posZ(), collision.nGlobalTracks()); - // purity + // reconstructed + fillResolvedCorrectionHists(collision, pipm); + // check mc if (!pipm.jetTrack_as().has_mcParticle()) continue; - auto const pipmParticle = pipm.jetTrack_as().mcParticle(); - if (std::abs(pipmParticle.pdgCode()) != PDG_t::kPiPlus || !checkPrimaryEtaMc(pipmParticle)) + auto const& pipmParticle = pipm.jetTrack_as().mcParticle(); + // collision association + if (pipmParticle.mcCollisionId() != collision.mcCollisionId()) continue; - if (requireSingleCollisionPurity && pipmParticle.mcCollisionId() != collision.mcCollisionId()) + fillResolvedCorrectionHists(collision, pipm); + // purity + if (std::abs(pipmParticle.pdgCode()) != PDG_t::kPiPlus || !pipmParticle.isPhysicalPrimary()) continue; - - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcReco_hasCorrectMc_pipm"), pipm.pt(), pipm.phi(), pipm.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcReco_hasCorrectMc_pipm"), pipm.pt(), collision.posZ(), collision.nGlobalTracks()); + fillResolvedCorrectionHists(collision, pipm); } // photonPCM for (auto const& photonPCM : photonPCMs) { - if (doTrigEvEff && !collision.trigEv()) + // reconstructed + fillResolvedCorrectionHists(collision, photonPCM); + // check mc + auto const& posTrack = photonPCM.posJetTrack_as(); + auto const& negTrack = photonPCM.negJetTrack_as(); + if (!posTrack.has_mcParticle() || !negTrack.has_mcParticle()) continue; - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcReco_photonPCM"), photonPCM.pt(), photonPCM.phi(), photonPCM.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcReco_photonPCM"), photonPCM.pt(), collision.posZ(), collision.nGlobalTracks()); - + // collision association + if (posTrack.mcParticle().mcCollisionId() != collision.mcCollisionId() || negTrack.mcParticle().mcCollisionId() != collision.mcCollisionId()) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_recoImpurities_photonPCM"), + collision.ptMax(), photonPCM.pt(), binValuePi0Impurities(PhotonImpurity::WrongEv)); + continue; + } + fillResolvedCorrectionHists(collision, photonPCM); // purity - // (V0Legs does not have the tracks reference as index column (just int)??) - auto const& posTrack = tracks.rawIteratorAt(photonPCM.posTrackId() - tracks.offset()); - auto const& negTrack = tracks.rawIteratorAt(photonPCM.negTrackId() - tracks.offset()); - if (!posTrack.has_mcParticle() || !negTrack.has_mcParticle()) + if (!isConversionPhoton(posTrack, negTrack)) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_recoImpurities_photonPCM"), + collision.ptMax(), photonPCM.pt(), binValuePi0Impurities(PhotonImpurity::Misid)); continue; - if (!isConversionPhoton(posTrack, negTrack) || !checkPrimaryEtaMc(*(posTrack.mcParticle().mothers_as().begin()))) + } + auto const& photons = posTrack.mcParticle().mothers_as(); + if (!photons.begin()->isPhysicalPrimary()) { + if (checkForMother(*(photons.begin()), PDG_t::kK0Short, true)) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_recoImpurities_photonPCM"), + collision.ptMax(), photonPCM.pt(), binValuePi0Impurities(PhotonImpurity::K0Short)); + continue; + } + if (checkForMother(*(photons.begin()), PDG_t::kK0Long, true)) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_recoImpurities_photonPCM"), + collision.ptMax(), photonPCM.pt(), binValuePi0Impurities(PhotonImpurity::K0Long)); + continue; + } + if (checkForMother(*(photons.begin()), PDG_t::kLambda0, true)) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_recoImpurities_photonPCM"), + collision.ptMax(), photonPCM.pt(), binValuePi0Impurities(PhotonImpurity::Lambda0)); + continue; + } + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_recoImpurities_photonPCM"), + collision.ptMax(), photonPCM.pt(), binValuePi0Impurities(PhotonImpurity::Other)); continue; - if (requireSingleCollisionPurity && posTrack.mcParticle().mcCollisionId() != collision.mcCollisionId()) + } + fillResolvedCorrectionHists(collision, photonPCM); + } + + // h0 PCM + for (auto const& photonPCMPair : photonPCMPairs) { + // check mc + auto const& posTrack1 = photonPCMPair.posJetTrack1_as(); + auto const& negTrack1 = photonPCMPair.negJetTrack1_as(); + auto const& posTrack2 = photonPCMPair.posJetTrack2_as(); + auto const& negTrack2 = photonPCMPair.negJetTrack2_as(); + if (!posTrack1.has_mcParticle() || !negTrack1.has_mcParticle() || !posTrack2.has_mcParticle() || !negTrack2.has_mcParticle()) continue; - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcReco_hasCorrectMc_photonPCM"), photonPCM.pt(), photonPCM.phi(), photonPCM.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcReco_hasCorrectMc_photonPCM"), photonPCM.pt(), collision.posZ(), collision.nGlobalTracks()); + // pi0PCM + if (isGGFromDoubleConversion(posTrack1, negTrack1, posTrack2, negTrack2, PDG_t::kPi0)) { + // pseudo reconstructed + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_pseudoReco_pi0PCM"), collision.ptMax(), photonPCMPair.pt()); + // collision association + if (posTrack1.mcParticle().mcCollisionId() != collision.mcCollisionId()) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_pseudoRecoImpurities_pi0PCM"), + collision.ptMax(), photonPCMPair.pt(), binValuePi0Impurities(PhotonImpurity::WrongEv)); + continue; + } + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_pseudoRecoAssocEv_pi0PCM"), collision.ptMax(), photonPCMPair.pt()); + // purity + // note: do not dereference 'mothers' and store element in 'const&' (leads to seg fault) + auto const& photons1 = posTrack1.mcParticle().mothers_as(); + auto const& mothersOfPhoton = photons1.begin()->mothers_as(); + auto const& grandmothersOfPhoton = mothersOfPhoton.begin()->mothers_as(); + if (!checkDecayPrimary(*(mothersOfPhoton.begin()))) { + if (std::abs(grandmothersOfPhoton.begin()->pdgCode()) == PDG_t::kK0Short) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_pseudoRecoImpurities_pi0PCM"), + collision.ptMax(), photonPCMPair.pt(), binValuePi0Impurities(PhotonImpurity::K0Short)); + continue; + } + if (std::abs(grandmothersOfPhoton.begin()->pdgCode()) == PDG_t::kK0Long) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_pseudoRecoImpurities_pi0PCM"), + collision.ptMax(), photonPCMPair.pt(), binValuePi0Impurities(PhotonImpurity::K0Long)); + continue; + } + if (std::abs(grandmothersOfPhoton.begin()->pdgCode()) == PDG_t::kLambda0) { + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_pseudoRecoImpurities_pi0PCM"), + collision.ptMax(), photonPCMPair.pt(), binValuePi0Impurities(PhotonImpurity::Lambda0)); + continue; + } + histos.fill(HIST("mc/correction/h3_ptTrigPtAssocCategory_pseudoRecoImpurities_pi0PCM"), + collision.ptMax(), photonPCMPair.pt(), binValuePi0Impurities(PhotonImpurity::Other)); + continue; + } + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_pseudoRecoPure_pi0PCM"), + collision.ptMax(), photonPCMPair.pt()); + continue; + } + // etaPCM + if (isGGFromDoubleConversion(posTrack1, negTrack1, posTrack2, negTrack2, constants::physics::Pdg::kEta)) { + // pseudo reconstructed + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_pseudoReco_etaPCM"), collision.ptMax(), photonPCMPair.pt()); + // collision association + if (posTrack1.mcParticle().mcCollisionId() != collision.mcCollisionId()) + continue; + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_pseudoRecoAssocEv_etaPCM"), collision.ptMax(), photonPCMPair.pt()); + // purity + auto const& photons1 = posTrack1.mcParticle().mothers_as(); + auto const& mothersOfPhoton = photons1.begin()->mothers_as(); + if (!checkDecayPrimary(*(mothersOfPhoton.begin()))) + continue; + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_pseudoRecoPure_etaPCM"), collision.ptMax(), photonPCMPair.pt()); + continue; + } } // mcParticle loop for (auto const& mcParticle : mcParticlesThisEvent) { - bool const countChargedTrigEvEff = !doTrigEvEff || collision.trigEv() || mcParticle.pt() > ptCutTrigEvEff; - bool const countOtherTrigEvEff = !doTrigEvEff || collision.trigEv(); + if (std::abs(mcParticle.eta()) > absEtaMax) + continue; // standard particles (marked physical primary) - if (checkPrimaryEtaMc(mcParticle)) { + if (mcParticle.isPhysicalPrimary()) { + // charged // hadrons - if (checkChargedMc(mcParticle) && countChargedTrigEvEff) { - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcTrue_recoCol_hadron"), mcParticle.pt(), mcParticle.phi(), mcParticle.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcTrue_recoCol_hadron"), mcParticle.pt(), collision.mcCollision_as().posZ(), collision.nGlobalTracks()); + if (checkChargedMc(mcParticle)) { + fillResolvedCorrectionHists(collision, mcParticle); + fillResolvedCorrectionHists(mcCollision, mcParticle); } // pipm - if (std::abs(mcParticle.pdgCode()) == PDG_t::kPiPlus && countChargedTrigEvEff) { - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcTrue_recoCol_pipm"), mcParticle.pt(), mcParticle.phi(), mcParticle.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcTrue_recoCol_pipm"), mcParticle.pt(), collision.mcCollision_as().posZ(), collision.nGlobalTracks()); + if (std::abs(mcParticle.pdgCode()) == PDG_t::kPiPlus) { + fillResolvedCorrectionHists(collision, mcParticle); + fillResolvedCorrectionHists(mcCollision, mcParticle); + continue; } + // non charged // photons - if (mcParticle.pdgCode() == PDG_t::kGamma && countOtherTrigEvEff) { - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcTrue_recoCol_photon"), mcParticle.pt(), mcParticle.phi(), mcParticle.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcTrue_recoCol_photon"), mcParticle.pt(), collision.mcCollision_as().posZ(), collision.nGlobalTracks()); + if (mcParticle.pdgCode() == PDG_t::kGamma) { + fillResolvedCorrectionHists(collision, mcParticle); + fillResolvedCorrectionHists(mcCollision, mcParticle); + continue; } + continue; } // decaying particles (not marked physical primary) - if ((std::abs(mcParticle.eta()) < etaMax)) { - // pi0 - if (checkH0ToGG(mcParticle, PDG_t::kPi0) && countOtherTrigEvEff) { - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcTrue_recoCol_pi0"), mcParticle.pt(), mcParticle.phi(), mcParticle.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcTrue_recoCol_pi0"), mcParticle.pt(), collision.mcCollision_as().posZ(), collision.nGlobalTracks()); + if (!checkDecayPrimary(mcParticle)) + continue; + // pi0 + if (mcParticle.pdgCode() == PDG_t::kPi0) { + // true + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_trueAssocEv_pi0"), mcCollision.ptMax(), mcParticle.pt()); + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_trueAssocEvRecoPtTrig_pi0"), collision.ptMax(), mcParticle.pt()); + // chosen decay + if (!checkToGG(mcParticle)) + continue; + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_trueMeasDecay_pi0"), collision.ptMax(), mcParticle.pt()); + // daughters in acceptance + if (!checkDaughtersInAcceptance(mcParticle, {-absEtaMax, absEtaMax})) + continue; + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_trueGeoAcc_pi0"), collision.ptMax(), mcParticle.pt()); + continue; + } + // eta + if (mcParticle.pdgCode() == constants::physics::Pdg::kEta) { + // true + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_trueAssocEv_eta"), mcCollision.ptMax(), mcParticle.pt()); + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_trueAssocEvRecoPtTrig_eta"), collision.ptMax(), mcParticle.pt()); + // chosen decay + if (!checkToGG(mcParticle)) + continue; + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_trueMeasDecay_eta"), collision.ptMax(), mcParticle.pt()); + // daughters in acceptance + if (!checkDaughtersInAcceptance(mcParticle, {-absEtaMax, absEtaMax})) + continue; + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_trueGeoAcc_eta"), collision.ptMax(), mcParticle.pt()); + continue; + } + } + } + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcRecoCorrection, "process mc reconstruction to calculate corrections", false); + + void processMcTrueCorrection(CorrMcCollision const& mcCollision, aod::JetParticles const& mcParticles) + { + // event selection + if (!totalEvSel(mcCollision)) + return; + + for (auto const& mcParticle : mcParticles) { + if (std::abs(mcParticle.eta()) > absEtaMax) + continue; + + // standard particles (marked physical primary) + if (mcParticle.isPhysicalPrimary()) { + // charged + // hadrons + if (checkChargedMc(mcParticle)) { + fillResolvedCorrectionHists(mcCollision, mcParticle); } - // eta - if (checkH0ToGG(mcParticle, 221) && countOtherTrigEvEff) { - histos.fill(HIST("mc/eff/h3_ptPhiEta_mcTrue_recoCol_eta"), mcParticle.pt(), mcParticle.phi(), mcParticle.eta()); - histos.fill(HIST("mc/eff/h3_ptZPvMult_mcTrue_recoCol_eta"), mcParticle.pt(), collision.mcCollision_as().posZ(), collision.nGlobalTracks()); + // pipm + if (std::abs(mcParticle.pdgCode()) == PDG_t::kPiPlus) { + fillResolvedCorrectionHists(mcCollision, mcParticle); + continue; } + // non charged + // photons + if (mcParticle.pdgCode() == PDG_t::kGamma) { + fillResolvedCorrectionHists(mcCollision, mcParticle); + continue; + } + continue; + } + + // decaying particles (not marked physical primary) + if (!checkDecayPrimary(mcParticle)) + continue; + // pi0 + if (mcParticle.pdgCode() == PDG_t::kPi0) { + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_true_pi0"), mcCollision.ptMax(), mcParticle.pt()); + continue; + } + // eta + if (mcParticle.pdgCode() == constants::physics::Pdg::kEta) { + histos.fill(HIST("mc/correction/h2_ptTrigPtAssoc_true_eta"), mcCollision.ptMax(), mcParticle.pt()); + continue; } } } - PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcRecoColEff, "process MC data to calculate efficiencies and purities", false); - - // test ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcTrueCorrection, "process mc true data to calculate corrections", false); - void processTest(CorrCollision const& collision, - soa::Join const& tracks, soa::Join const&, - aod::Hadrons const& hadrons) + void processPCMDCAzTrue(CorrMcDCollision const& collision, aod::PhotonPCMs const& photonPCMs, aod::V0PhotonsKF const&, aod::JetTracksMCD const&, + CorrMcCollisions const&, aod::JetParticles const&) { // event selection - if (!collision.selEv()) + if (!totalEvSel(collision)) return; - histos.fill(HIST("test/h2_mult_comp"), collision.nGlobalTracks(), hadrons.size()); + for (auto const& photonPCM : photonPCMs) { + // check mc + auto const& posTrack = photonPCM.posJetTrack_as(); + auto const& negTrack = photonPCM.negJetTrack_as(); + if (!posTrack.has_mcParticle() || !negTrack.has_mcParticle()) + continue; + // collision association + if (posTrack.mcParticle().mcCollisionId() != collision.mcCollisionId() || negTrack.mcParticle().mcCollisionId() != collision.mcCollisionId()) { + continue; + } + + histos.fill(HIST("mc/plain/h5_ptTrigPtAssocDCAzZPvMult_photonPCM"), + collision.ptMax(), photonPCM.pt(), photonPCM.v0PhotonKF().dcaZtopv(), collision.posZ(), collision.nGlobalTracks()); + } + } + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processPCMDCAzTrue, "process to test true DCAz distribution of photons in trigger collisions", false); - for (auto const& track : tracks) { - auto const fullTrack = track.track_as>(); + // last ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ - constexpr float Mincrossedrows = 40; - constexpr float Maxchi2tpc = 5.0; - constexpr float Maxchi2its = 6.0; - constexpr float MaxR = 83.1; - constexpr float MinPtTrackiu = 0.1; + void processRecoLast(CorrCollisions const& collisions, aod::Triggers const& triggers) + { + // do at end of each data frame (after other reco correlation process functions) + // (PROCESS_SWITCH of this process has to be declared last) - if (!fullTrack.hasITS() && !fullTrack.hasTPC()) + for (auto const& collision : collisions) { + // all events + histos.fill(HIST("reco/info/h1_nEvents"), 0.5); + + // mc split + if (!checkSplitMcEventSelection(collision)) continue; - if (fullTrack.x() * fullTrack.x() + fullTrack.y() * fullTrack.y() > MaxR * MaxR || fullTrack.pt() < MinPtTrackiu) + histos.fill(HIST("reco/info/h1_nEvents"), 1.5); + + // standard event selection + if (!collision.selEv()) continue; - if (fullTrack.hasTPC()) { - if (fullTrack.tpcNClsCrossedRows() < Mincrossedrows || fullTrack.tpcChi2NCl() > Maxchi2tpc) - continue; - } - if (fullTrack.hasITS()) { - if (fullTrack.itsChi2NCl() > Maxchi2its) - continue; + histos.fill(HIST("reco/info/h1_nEvents"), 2.5); + + histos.fill(HIST("reco/info/h3_ptTrigZPvMult"), collision.ptMax(), collision.posZ(), collision.nGlobalTracks()); + histos.fill(HIST("reco/info/h2_ptTrigOccupancy"), collision.ptMax(), collision.trackOccupancyInTimeRange()); + + // group collision + auto const triggersThisEvent = triggers.sliceBy(perColTriggers, collision.globalIndex()); + + // trigger loop + for (auto const& trigger : triggersThisEvent) { + // trigger info + histos.fill(HIST("reco/corr/h3_ptPhiEta_trig"), trigger.pt(), trigger.phi(), trigger.eta(), + getPtCorrection(trigger.pt())); + + // save trigger for mixing + triggerMemoryReco->saveTrigger(trigger.pt(), trigger.phi(), trigger.eta(), collision.posZ(), collision.nGlobalTracks()); } + } + } + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processRecoLast, "process general info for reconstructed collisions and later correlations last", true); + + void processMcLast(CorrMcDCollisions const& collisions, aod::Triggers const& triggers, aod::JetTracksMCD const&, + CorrMcCollisions const& mcCollisions, aod::TriggerParticles const& triggerParticles, aod::JetParticles const&) + { + for (auto const& mcCollision : mcCollisions) { + // event selection + if (!totalEvSel(mcCollision)) + continue; + + // group collision + auto const triggerParticlesThisEvent = triggerParticles.sliceBy(perColTriggerParticles, mcCollision.globalIndex()); + + histos.fill(HIST("mc/info/h3_ptTrigZPvMult_true"), mcCollision.ptMax(), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); - histos.fill(HIST("test/h2_tracks_zPvMultDep"), collision.posZ(), collision.nGlobalTracks()); + // trigger loop + for (auto const& trigger : triggerParticlesThisEvent) { + // trigger info + histos.fill(HIST("mc/corr/h3_ptPhiEta_trig_true"), trigger.pt(), trigger.phi(), trigger.eta()); + + // save trigger for mixing + triggerMemoryTrue->saveTrigger(trigger.pt(), trigger.phi(), trigger.eta(), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); + } } + for (auto const& collision : collisions) { + // event selection + if (!totalEvSel(collision)) + continue; + + auto const& mcCollision = collision.mcCollision_as(); - histos.fill(HIST("test/h2_globalTracks_zPvMultDep"), collision.posZ(), collision.nGlobalTracks(), hadrons.size()); + // group collision + auto const triggerParticlesThisEvent = triggerParticles.sliceBy(perColTriggerParticles, mcCollision.globalIndex()); + auto const triggersThisEvent = triggers.sliceBy(perColTriggers, collision.globalIndex()); + + histos.fill(HIST("mc/info/h3_ptTrigZPvMult_trueAssocEv"), mcCollision.ptMax(), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); + histos.fill(HIST("mc/info/h3_ptTrigZPvMult_trueAssocEvRecoPtTrig"), collision.ptMax(), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); + + // trigger loops + for (auto const& triggerParticle : triggerParticlesThisEvent) { + // trigger info + histos.fill(HIST("mc/corr/h3_ptPhiEta_trig_trueAssocEv"), triggerParticle.pt(), triggerParticle.phi(), triggerParticle.eta()); + + // save trigger for mixing + triggerMemoryTrueAssocEv->saveTrigger(triggerParticle.pt(), triggerParticle.phi(), triggerParticle.eta(), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); + } + for (auto const& trigger : triggersThisEvent) { + if (!trigger.jetTrack_as().has_mcParticle()) + continue; + + auto const& triggerParticle = trigger.jetTrack_as().mcParticle(); + + // trigger info + histos.fill(HIST("mc/corr/h3_ptPhiEta_trig_trueAssocEvRecoPtTrig"), trigger.pt(), triggerParticle.phi(), triggerParticle.eta()); + + // save trigger for mixing + triggerMemoryTrueAssocEvRecoPtTrig->saveTrigger(trigger.pt(), triggerParticle.phi(), triggerParticle.eta(), mcCollision.posZ(), mcCollision.nChargedInEtaRange()); + } + } } - PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processTest, "process just to test things", false); + PROCESS_SWITCH(PhotonChargedTriggerCorrelation, processMcLast, "process general MC info for collisions and later correlations last", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& configContext) diff --git a/PWGJE/Tasks/photonChargedTriggerProducer.cxx b/PWGJE/Tasks/photonChargedTriggerProducer.cxx index 4ffc372ad67..e8507e86fa1 100644 --- a/PWGJE/Tasks/photonChargedTriggerProducer.cxx +++ b/PWGJE/Tasks/photonChargedTriggerProducer.cxx @@ -26,20 +26,26 @@ #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" +#include #include #include #include #include #include +#include +#include #include #include +#include #include #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include #include +#include #include +#include #include using namespace o2; @@ -47,24 +53,36 @@ using namespace o2::framework; // correlation derived data =================================================================================================================================================================== +template +concept IsJetMcCollision = std::same_as, aod::JetMcCollision>; +template +concept IsJetCollisionMCD = std::same_as, aod::JetCollisionMCD>; +template +concept HasTrueCollision = IsJetMcCollision || IsJetCollisionMCD; + struct PhotonChargedTriggerProducer { - // reco - Produces collisionExtraCorrTable; - Produces triggerTable; - Produces hadronTable; - Produces pipmTable; - Produces photonPCMTable; - Produces photonPCMPairTable; - // mc - Produces mcCollisionExtraCorrTable; - Produces triggerParticleTable; + struct : ProducesGroup { + // reco + Produces collisionExtraCorr; + Produces trigger; + Produces hadron; + Produces pipm; + Produces photonPCM; + Produces photonPCMPair; + // mc + Produces mcCollisionExtraCorr; + Produces triggerParticle; + } tableProduction; Configurable zPvMax{"zPvMax", 7, "maximum absZ primary-vertex cut"}; Configurable occupancyMin{"occupancyMin", 0, "minimum occupancy cut"}; Configurable occupancyMax{"occupancyMax", 2000, "maximum occupancy cut"}; - Configurable etaMax{"etaMax", 0.8, "maximum absEta cut"}; + Configurable absEtaMax{"absEtaMax", 0.8, "maximum absEta cut"}; Configurable eventSelections{"eventSelections", "sel8", "JE framework - event selection"}; + Configurable skipMBGapEvents{"skipMBGapEvents", true, "skip MB gap events"}; + Configurable applyRctSelection{"applyRctSelection", true, "apply RCT selection"}; + Configurable rctLabel{"rctLabel", "CBT_hadronPID", "RCT selection label"}; Configurable trackSelections{"trackSelections", "globalTracks", "JE framework - track selections"}; Configurable triggerMasks{"triggerMasks", "", "JE framework - skimmed data trigger masks (relevent for correlation: fTrackLowPt,fTrackHighPt)"}; @@ -75,7 +93,10 @@ struct PhotonChargedTriggerProducer { Configurable> nSigmaPiTof{"nSigmaPiTof", {-1, 2}, "minimum-maximum nSigma for pipm in tof"}; Configurable> nSigmaPiRelRise{"nSigmaPiRelRise", {0, 2}, "minimum-maximum nSigma pipm tpc at high pt"}; + Configurable> pcmAbsDcaZLim{"pcmAbsDcaZLim", {0, 5}, "absolute values of DCAz limits for PCM"}; + Configurable ptTrigMin{"ptTrigMin", 5, "minimum pT of triggers"}; + Configurable oneTriggerPerEvent{"oneTriggerPerEvent", false, "use only one trigger (of highest pt) in an event"}; // derivatives of configurables @@ -83,6 +104,8 @@ struct PhotonChargedTriggerProducer { int trackSelection = -1; std::vector triggerMaskBits; + HistogramRegistry histos{"histogramRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + // for mc Service pdg; @@ -104,7 +127,7 @@ struct PhotonChargedTriggerProducer { { if (!jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) return false; - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents, applyRctSelection, rctLabel)) return false; if (std::abs(collision.posZ()) > zPvMax) return false; @@ -113,13 +136,28 @@ struct PhotonChargedTriggerProducer { return true; } + // mc split event selection based in the thrid decimal of mc-true posZ + template + int getThirdDecimalTruePosZ(T_collision const& collision) + { + // get third decimal + double absPosZ; + if constexpr (IsJetCollisionMCD) { + absPosZ = std::abs(collision.template mcCollision_as().posZ()); + } else { + absPosZ = std::abs(collision.posZ()); + } + int const thirdDecimalPosZ = static_cast(1000 * absPosZ) % 10; + return thirdDecimalPosZ; + } + // checks global track cuts template bool checkGlobalTrackEta(T_track const& track) { if (!jetderiveddatautilities::selectTrack(track, trackSelection)) return false; - if (!jetderiveddatautilities::applyTrackKinematics(track, 0.1, 1000, -1 * etaMax, etaMax)) + if (!jetderiveddatautilities::applyTrackKinematics(track, 0.1, 1000, -1 * absEtaMax, absEtaMax)) return false; return true; } @@ -172,14 +210,44 @@ struct PhotonChargedTriggerProducer { eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(triggerMasks); + + // histograms + AxisSpec const axisZPv{28, -7, 7, "#it{z}_{pv}"}; + AxisSpec const axisMult{65, -0.5, 64.5, "multiplicity"}; + AxisSpec const axisHadronicRate{1001, -0.5, 1000.5, "hadronicRate"}; + AxisSpec const axisPt{50, 0, 25, "#it{p}_{T}"}; + AxisSpec const axisPhi{72, 0, constants::math::TwoPI, "#varphi"}; + AxisSpec const axisEta{32, -absEtaMax, absEtaMax, "#eta"}; + // reco + histos.add("reco/h2_zPvMult", "h2_zPvMult", kTHnSparseD, {axisZPv, axisMult}, true); + histos.add("reco/h1_hadronicRate", "h1_hadronicRate", kTH1D, {axisHadronicRate}, true); + histos.add("reco/h3_ptPhiEta_trigger", "h3_ptPhiEta_trigger", kTHnSparseD, {axisPt, axisPhi, axisEta}, true); + histos.add("reco/h3_ptPhiEta_hadron", "h3_ptPhiEta_hadron", kTHnSparseD, {axisPt, axisPhi, axisEta}, true); + histos.add("reco/h3_ptPhiEta_pipm", "h3_ptPhiEta_pipm", kTHnSparseD, {axisPt, axisPhi, axisEta}, true); + histos.add("reco/h3_ptPhiEta_photonPCM", "h3_ptPhiEta_photonPCM", kTHnSparseD, {axisPt, axisPhi, axisEta}, true); + histos.add("reco/h3_ptPhiEta_photonPCMPair", "h3_ptPhiEta_photonPCMPair", kTHnSparseD, {axisPt, axisPhi, axisEta}, true); + // true + histos.add("true/h2_zPvMult", "h2_zPvMult", kTHnSparseD, {axisZPv, axisMult}, true); + histos.add("true/h3_ptPhiEta_trigger", "h3_ptPhiEta_trigger", kTHnSparseD, {axisPt, axisPhi, axisEta}, true); } - void processRecoCollisionTrigger(aod::JetCollision const& collision, aod::JetTracks const& tracks) + template + void fillRecoCollisionTrigger(T_collision const& collision, T_tracks const& tracks) { // event selection - const bool isSelectedEvent = checkEventSelection(collision); - // trigger event check - bool isTriggerEvent = false; + bool isSelectedEvent = checkEventSelection(collision); + // mc split event selection + // default -1 for data reconstruction + int thirdDecimalTruePosZ = -1; + if constexpr (IsJetCollisionMCD) { + thirdDecimalTruePosZ = getThirdDecimalTruePosZ(collision); + } + // maximum track pt + float maxTrackPt = -1; + int maxTriggerCollisionId = -1; + int maxTriggerGlobalIndex = -1; + float maxTriggerPhi = 0; + float maxTriggerEta = 0; // number global tracks int nGlobalTracks = 0; @@ -190,23 +258,49 @@ struct PhotonChargedTriggerProducer { continue; nGlobalTracks++; + if (track.pt() > maxTrackPt) { + maxTrackPt = track.pt(); + maxTriggerCollisionId = track.collisionId(); + maxTriggerGlobalIndex = track.globalIndex(); + maxTriggerPhi = track.phi(); + maxTriggerEta = track.eta(); + } if (!isSelectedEvent) continue; if (track.pt() < ptTrigMin) continue; - - isTriggerEvent = true; + if (oneTriggerPerEvent) + continue; // trigger info - triggerTable(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); + histos.fill(HIST("reco/h3_ptPhiEta_trigger"), track.pt(), track.phi(), track.eta()); + tableProduction.trigger(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); + } + bool const isTriggerEvent = maxTrackPt > ptTrigMin; + if (oneTriggerPerEvent && isSelectedEvent && isTriggerEvent) { + histos.fill(HIST("reco/h3_ptPhiEta_trigger"), maxTrackPt, maxTriggerPhi, maxTriggerEta); + tableProduction.trigger(maxTriggerCollisionId, maxTriggerGlobalIndex, maxTrackPt, maxTriggerPhi, maxTriggerEta); } // collision info - collisionExtraCorrTable(isSelectedEvent, isTriggerEvent, nGlobalTracks); + histos.fill(HIST("reco/h2_zPvMult"), collision.posZ(), nGlobalTracks); + histos.fill(HIST("reco/h1_hadronicRate"), collision.hadronicRate()); + tableProduction.collisionExtraCorr(isSelectedEvent, thirdDecimalTruePosZ, isTriggerEvent, maxTrackPt, nGlobalTracks); + } + + void processRecoCollisionTrigger(aod::JetCollision const& collision, aod::JetTracks const& tracks) + { + fillRecoCollisionTrigger(collision, tracks); } PROCESS_SWITCH(PhotonChargedTriggerProducer, processRecoCollisionTrigger, "process correlation collision_extra and trigger table (reconstructed)", false); + void processMcRecoCollisionTrigger(aod::JetCollisionMCD const& collision, aod::JetTracks const& tracks, aod::JetMcCollisions const&) + { + fillRecoCollisionTrigger(collision, tracks); + } + PROCESS_SWITCH(PhotonChargedTriggerProducer, processMcRecoCollisionTrigger, "process correlation collision_extra and trigger table (reconstructed MC)", false); + void processRecoPipmTPCTOF(aod::JetCollision const& collision, soa::Join const& tracks, soa::Join const&) { @@ -221,7 +315,8 @@ struct PhotonChargedTriggerProducer { continue; // hadron - hadronTable(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); + histos.fill(HIST("reco/h3_ptPhiEta_hadron"), track.pt(), track.phi(), track.eta()); + tableProduction.hadron(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); // pipm selection auto const& trackPID = track.track_as>(); @@ -229,7 +324,8 @@ struct PhotonChargedTriggerProducer { continue; // pipm - pipmTable(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); + histos.fill(HIST("reco/h3_ptPhiEta_pipm"), track.pt(), track.phi(), track.eta()); + tableProduction.pipm(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); } } PROCESS_SWITCH(PhotonChargedTriggerProducer, processRecoPipmTPCTOF, "process pipm (TPC-TOF) table (reconstructed)", false); @@ -248,7 +344,8 @@ struct PhotonChargedTriggerProducer { continue; // hadron - hadronTable(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); + histos.fill(HIST("reco/h3_ptPhiEta_hadron"), track.pt(), track.phi(), track.eta()); + tableProduction.hadron(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); // pipm selection auto const& trackPID = track.track_as>(); @@ -256,7 +353,8 @@ struct PhotonChargedTriggerProducer { continue; // pipm - pipmTable(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); + histos.fill(HIST("reco/h3_ptPhiEta_pipm"), track.pt(), track.phi(), track.eta()); + tableProduction.pipm(track.collisionId(), track.globalIndex(), track.pt(), track.phi(), track.eta()); } } PROCESS_SWITCH(PhotonChargedTriggerProducer, processRecoPipmTPC, "process pipm (TPC) table (reconstructed)", false); @@ -274,12 +372,16 @@ struct PhotonChargedTriggerProducer { // photonPCM for (auto const& v0Photon : v0PhotonsThisEvent) { // photon selection - if (std::abs(v0Photon.eta()) > etaMax) + if (std::abs(v0Photon.eta()) > absEtaMax) continue; + if (std::abs(v0Photon.dcaZtopv()) < pcmAbsDcaZLim.value[0] || std::abs(v0Photon.dcaZtopv()) > pcmAbsDcaZLim.value[1]) { + continue; + } // photon PCM - photonPCMTable(v0Photon.collisionId(), v0Photon.globalIndex(), - v0Photon.posTrack().trackId(), v0Photon.negTrack().trackId(), v0Photon.pt(), v0Photon.phi(), v0Photon.eta()); + histos.fill(HIST("reco/h3_ptPhiEta_photonPCM"), v0Photon.pt(), v0Photon.phi(), v0Photon.eta()); + tableProduction.photonPCM(v0Photon.collisionId(), v0Photon.globalIndex(), + v0Photon.posTrack().trackId(), v0Photon.negTrack().trackId(), v0Photon.pt(), v0Photon.phi(), v0Photon.eta()); } // photonPCm pairs @@ -290,21 +392,28 @@ struct PhotonChargedTriggerProducer { ROOT::Math::PtEtaPhiMVector const p4V0PCMPair = p4V0PCM1 + p4V0PCM2; // pi0 selection - if (std::abs(p4V0PCMPair.Eta()) > etaMax) + if (std::abs(p4V0PCMPair.Eta()) > absEtaMax) continue; // save info - photonPCMPairTable(v0Photon1.collisionId(), v0Photon1.globalIndex(), v0Photon2.globalIndex(), - v0Photon1.posTrack().trackId(), v0Photon1.negTrack().trackId(), v0Photon2.posTrack().trackId(), v0Photon2.negTrack().trackId(), - p4V0PCMPair.Pt(), RecoDecay::constrainAngle(p4V0PCMPair.Phi(), 0), p4V0PCMPair.Eta(), p4V0PCMPair.M()); + histos.fill(HIST("reco/h3_ptPhiEta_photonPCMPair"), p4V0PCMPair.Pt(), RecoDecay::constrainAngle(p4V0PCMPair.Phi(), 0), p4V0PCMPair.Eta()); + tableProduction.photonPCMPair(v0Photon1.collisionId(), v0Photon1.globalIndex(), v0Photon2.globalIndex(), + v0Photon1.posTrack().trackId(), v0Photon1.negTrack().trackId(), v0Photon2.posTrack().trackId(), v0Photon2.negTrack().trackId(), + p4V0PCMPair.Pt(), RecoDecay::constrainAngle(p4V0PCMPair.Phi(), 0), p4V0PCMPair.Eta(), p4V0PCMPair.M()); } } PROCESS_SWITCH(PhotonChargedTriggerProducer, processRecoPhotonPCM, "process photonPCM table (reconstructed)", false); - void processMcCorrTables(aod::JetMcCollision const&, aod::JetParticles const& mcParticles) + void processMcCorrTables(aod::JetMcCollision const& mcCollision, aod::JetParticles const& mcParticles) { - // trigger event check - bool isTriggerEvent = false; + // mc split event selection + int const thirdDecimalTruePosZ = getThirdDecimalTruePosZ(mcCollision); + // maximum charged particle pt + float maxChargedPt = -1; + int maxTriggerCollisionId = -1; + int maxTriggerGlobalIndex = -1; + float maxTriggerPhi = 0; + float maxTriggerEta = 0; // number charged particles in eta range int nCharged = 0; @@ -316,23 +425,37 @@ struct PhotonChargedTriggerProducer { continue; if (!mcParticle.isPhysicalPrimary()) continue; - if (std::abs(mcParticle.eta()) > etaMax) + if (std::abs(mcParticle.eta()) > absEtaMax) continue; nCharged++; + if (mcParticle.pt() > maxChargedPt) { + maxChargedPt = mcParticle.pt(); + maxTriggerCollisionId = mcParticle.mcCollisionId(); + maxTriggerGlobalIndex = mcParticle.globalIndex(); + maxTriggerPhi = mcParticle.phi(); + maxTriggerEta = mcParticle.eta(); + } // trigger selection if (mcParticle.pt() < ptTrigMin) continue; - - isTriggerEvent = true; + if (oneTriggerPerEvent) + continue; // trigger info - triggerParticleTable(mcParticle.mcCollisionId(), mcParticle.globalIndex(), mcParticle.pt(), mcParticle.phi(), mcParticle.eta()); + histos.fill(HIST("true/h3_ptPhiEta_trigger"), mcParticle.pt(), mcParticle.phi(), mcParticle.eta()); + tableProduction.triggerParticle(mcParticle.mcCollisionId(), mcParticle.globalIndex(), mcParticle.pt(), mcParticle.phi(), mcParticle.eta()); + } + bool const isTriggerEvent = maxChargedPt > ptTrigMin; + if (oneTriggerPerEvent && isTriggerEvent) { + histos.fill(HIST("true/h3_ptPhiEta_trigger"), maxChargedPt, maxTriggerPhi, maxTriggerEta); + tableProduction.triggerParticle(maxTriggerCollisionId, maxTriggerGlobalIndex, maxChargedPt, maxTriggerPhi, maxTriggerEta); } // collision info - mcCollisionExtraCorrTable(isTriggerEvent, nCharged); + histos.fill(HIST("true/h2_zPvMult"), mcCollision.posZ(), nCharged); + tableProduction.mcCollisionExtraCorr(thirdDecimalTruePosZ, isTriggerEvent, maxChargedPt, nCharged); } PROCESS_SWITCH(PhotonChargedTriggerProducer, processMcCorrTables, "process table production (mc)", false); };