diff --git a/PWGLF/DataModel/LFResonanceTables.h b/PWGLF/DataModel/LFResonanceTables.h index fcbe3b6378a..8a04eb4ffdd 100644 --- a/PWGLF/DataModel/LFResonanceTables.h +++ b/PWGLF/DataModel/LFResonanceTables.h @@ -68,6 +68,7 @@ DECLARE_SOA_COLUMN(EvtPlResAB, evtPlResAB, float); // DECLARE_SOA_COLUMN(EvtPlResAC, evtPlResAC, float); //! Second harmonic event plane resolution of A-C sub events DECLARE_SOA_COLUMN(EvtPlResBC, evtPlResBC, float); //! Second harmonic event plane resolution of B-C sub events DECLARE_SOA_COLUMN(BMagField, bMagField, float); //! Magnetic field +DECLARE_SOA_COLUMN(IsRecINELgt0, isRecINELgt0, bool); //! Is reconstructed INEL>0 event // MC DECLARE_SOA_COLUMN(IsVtxIn10, isVtxIn10, bool); //! Vtx10 DECLARE_SOA_COLUMN(IsINELgt0, isINELgt0, bool); //! INEL>0 @@ -75,16 +76,20 @@ DECLARE_SOA_COLUMN(IsTriggerTVX, isTriggerTVX, bool); //! TriggerTVX DECLARE_SOA_COLUMN(IsInSel8, isInSel8, bool); //! InSel8 DECLARE_SOA_COLUMN(IsInAfterAllCuts, isInAfterAllCuts, bool); //! InAfterAllCuts DECLARE_SOA_COLUMN(ImpactParameter, impactParameter, float); //! ImpactParameter +DECLARE_SOA_COLUMN(MCMultiplicity, mcMultiplicity, float); //! MC Multiplicity } // namespace resocollision DECLARE_SOA_TABLE(ResoCollisions, "AOD", "RESOCOLLISION", o2::soa::Index<>, o2::aod::mult::MultNTracksPV, + o2::aod::mult::MultNTracksPVeta1, + o2::aod::mult::MultNTracksPVetaHalf, collision::PosX, collision::PosY, collision::PosZ, resocollision::Cent, - resocollision::BMagField); + resocollision::BMagField, + resocollision::IsRecINELgt0); using ResoCollision = ResoCollisions::iterator; DECLARE_SOA_TABLE(ResoCollisionColls, "AOD", "RESOCOLLISIONCOL", @@ -98,7 +103,8 @@ DECLARE_SOA_TABLE(ResoMCCollisions, "AOD", "RESOMCCOLLISION", resocollision::IsTriggerTVX, resocollision::IsInSel8, resocollision::IsInAfterAllCuts, - resocollision::ImpactParameter); + resocollision::ImpactParameter, + resocollision::MCMultiplicity); using ResoMCCollision = ResoMCCollisions::iterator; DECLARE_SOA_TABLE(ResoSpheroCollisions, "AOD", "RESOSPHEROCOLLISION", @@ -229,6 +235,8 @@ DECLARE_SOA_COLUMN(IsPhysicalPrimary, isPhysicalPrimary, bool); DECLARE_SOA_COLUMN(ProducedByGenerator, producedByGenerator, bool); DECLARE_SOA_COLUMN(MotherId, motherId, int); //! Id of the mother particle DECLARE_SOA_COLUMN(MotherPDG, motherPDG, int); //! PDG code of the mother particle +DECLARE_SOA_COLUMN(MotherPt, motherPt, float); //! Pt of the mother particle +DECLARE_SOA_COLUMN(MotherRap, motherRap, float); //! Rapidity of the mother particle DECLARE_SOA_COLUMN(DaughterPDG1, daughterPDG1, int); //! PDG code of the first Daughter particle DECLARE_SOA_COLUMN(DaughterPDG2, daughterPDG2, int); //! PDG code of the second Daughter particle DECLARE_SOA_COLUMN(DaughterID1, daughterID1, int); //! Id of the first Daughter particle @@ -812,6 +820,8 @@ DECLARE_SOA_TABLE(ResoMCV0s, "AOD", "RESOMCV0", mcparticle::PdgCode, resodaughter::MotherId, resodaughter::MotherPDG, + resodaughter::MotherPt, + resodaughter::MotherRap, resodaughter::DaughterID1, resodaughter::DaughterID2, resodaughter::DaughterPDG1, @@ -824,6 +834,8 @@ DECLARE_SOA_TABLE(ResoMCCascades, "AOD", "RESOMCCASCADE", mcparticle::PdgCode, resodaughter::MotherId, resodaughter::MotherPDG, + resodaughter::MotherPt, + resodaughter::MotherRap, resodaughter::BachTrkID, resodaughter::V0ID, resodaughter::DaughterPDG1, diff --git a/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx b/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx index 79514728d10..e1c4b3cc5d8 100644 --- a/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx +++ b/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx @@ -11,11 +11,12 @@ /// /// \file resonanceInitializer.cxx /// \brief Initializes variables for the resonance candidate producers -/// \author Bong-Hwi Lim +/// \author Bong-Hwi Lim , Minjae Kim /// #include "PWGLF/DataModel/LFResonanceTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/mcCentrality.h" #include "PWGLF/Utils/collisionCuts.h" #include "Common/Core/EventPlaneHelper.h" @@ -24,6 +25,7 @@ #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/McCollisionExtra.h" #include "Common/DataModel/Qvectors.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -205,6 +207,27 @@ struct ResonanceInitializer { Configurable cfgSecondaryCrossMassCutWindow{"cfgSecondaryCrossMassCutWindow", 0.05, "Secondary inv mass selection window with (anti)lambda hypothesis"}; } SecondaryCuts; + struct : ConfigurableGroup { + Configurable cfgGenMult05{"cfgGenMult05", true, "GenEvent: multiplicity in |eta| < 0.5"}; + Configurable cfgGenMult10{"cfgGenMult10", false, "GenEvent: multiplicity in |eta| < 1.0"}; + Configurable cfgGenMultPercentile{"cfgGenMultPercentile", false, "Inherit Centrality(Multiplicity) percentile from MC collision only using LF-mc-centrality task"}; + + Configurable isZvtxcutGen{"isZvtxcutGen", true, "z-vertex cut for the GenCollision"}; + Configurable cutzvertexGen{"cutzvertexGen", 10.0f, "z-vertex cut for the GenCollision"}; + Configurable checkIsTrueINELgt0{"checkIsTrueINELgt0", true, "Check true INEL>0 for the Gen. Collision"}; + + ConfigurableAxis ptAxisGen{"ptAxisGen", {400, 0.0f, 20.0f}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis multNTracksAxis{"multNTracksAxis", {500, 0.0f, +5000.0f}, "Number of charged particles"}; + ConfigurableAxis impactParameterAxis{"impactParameterAxis", {500, 0, 50}, "IP (fm)"}; + + Configurable isDaughterCheck{"isDaughterCheck", 1, "Check if the mother has the correct daughters when it is considered"}; + Configurable cfgRapidityCutGen{"cfgRapidityCutGen", 0.5, "Rapidity cut for the truth particle"}; + Configurable pdgTruthMother{"pdgTruthMother", 3324, "pdgcode for the truth mother e.g. Xi(1530) (3324)"}; + Configurable pdgTruthDaughter1{"pdgTruthDaughter1", 3312, "pdgcode for the daughter 1, e.g. Xi- 3312"}; + Configurable pdgTruthDaughter2{"pdgTruthDaughter2", 211, "pdgcode for the daughter 2, e.g. pi+ 211"}; + } GenCuts; + Configurable checkIsRecINELgt0{"checkIsRecINELgt0", true, "Check rec INEL>0 for the Rec. Collision"}; + HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::AnalysisObject}; // Pre-filters for efficient process @@ -245,6 +268,7 @@ struct ResonanceInitializer { || (nabs(aod::mcparticle::pdgCode) == 123324); // Xi(1820)-0 using ResoEvents = soa::Join; + using ResoEvents001 = soa::Join; using ResoRun2Events = soa::Join; using ResoEventsMC = soa::Join; using ResoRun2EventsMC = soa::Join; @@ -931,6 +955,22 @@ struct ResonanceInitializer { } return lMothersPDGs; }; + auto getMothersPt = [&](auto const& theMcParticle) { + std::vector lMothersPts{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother pdgcode lMother: %f", lMother.pt()); + lMothersPts.push_back(lMother.pt()); + } + return lMothersPts; + }; + auto getMothersRap = [&](auto const& theMcParticle) { + std::vector lMothersRaps{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother rap lMother: %f", lMother.y()); + lMothersRaps.push_back(lMother.y()); + } + return lMothersRaps; + }; auto getDaughtersIndeces = [&](auto const& theMcParticle) { std::vector lDaughtersIndeces{}; for (auto const& lDaughter : theMcParticle.template daughters_as()) { @@ -954,6 +994,8 @@ struct ResonanceInitializer { // ------ std::vector mothers = {-1, -1}; std::vector motherPDGs = {-1, -1}; + std::vector mothersPts = {-1.0f, -1.0f}; + std::vector mothersRaps = {-1.0f, -1.0f}; std::vector daughters = {-1, -1}; std::vector daughterPDGs = {-1, -1}; if (v0.has_mcParticle()) { @@ -961,10 +1003,13 @@ struct ResonanceInitializer { if (v0mc.has_mothers()) { mothers = getMothersIndeces(v0mc); motherPDGs = getMothersPDGCodes(v0mc); + mothersPts = getMothersPt(v0mc); + mothersRaps = getMothersRap(v0mc); } while (mothers.size() > 2) { mothers.pop_back(); motherPDGs.pop_back(); + mothersPts.pop_back(); } if (v0mc.has_daughters()) { daughters = getDaughtersIndeces(v0mc); @@ -978,6 +1023,8 @@ struct ResonanceInitializer { reso2mcv0s(v0mc.pdgCode(), mothers[0], motherPDGs[0], + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -988,6 +1035,8 @@ struct ResonanceInitializer { reso2mcv0s(0, mothers[0], motherPDGs[0], + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1017,6 +1066,22 @@ struct ResonanceInitializer { } return lMothersPDGs; }; + auto getMothersPt = [&](auto const& theMcParticle) { + std::vector lMothersPts{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother pdgcode lMother: %f", lMother.pt()); + lMothersPts.push_back(lMother.pt()); + } + return lMothersPts; + }; + auto getMothersRap = [&](auto const& theMcParticle) { + std::vector lMothersRaps{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother rap lMother: %f", lMother.y()); + lMothersRaps.push_back(lMother.y()); + } + return lMothersRaps; + }; auto getDaughtersIndeces = [&](auto const& theMcParticle) { std::vector lDaughtersIndeces{}; for (auto const& lDaughter : theMcParticle.template daughters_as()) { @@ -1042,15 +1107,20 @@ struct ResonanceInitializer { std::vector motherPDGs = {-1, -1}; std::vector daughters = {-1, -1}; std::vector daughterPDGs = {-1, -1}; + std::vector mothersPts = {-1.0f, -1.0f}; + std::vector mothersRaps = {-1.0f, -1.0f}; if (casc.has_mcParticle()) { auto cascmc = casc.mcParticle(); if (cascmc.has_mothers()) { mothers = getMothersIndeces(cascmc); + mothersPts = getMothersPt(cascmc); + mothersRaps = getMothersRap(cascmc); motherPDGs = getMothersPDGCodes(cascmc); } while (mothers.size() > 2) { mothers.pop_back(); motherPDGs.pop_back(); + mothersPts.pop_back(); } if (cascmc.has_daughters()) { daughters = getDaughtersIndeces(cascmc); @@ -1064,6 +1134,8 @@ struct ResonanceInitializer { reso2mccascades(cascmc.pdgCode(), mothers[0], motherPDGs[0], + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1074,6 +1146,8 @@ struct ResonanceInitializer { reso2mccascades(0, mothers[0], motherPDGs[0], + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1112,20 +1186,53 @@ struct ResonanceInitializer { } } + template + void fillMCGenParticles(TotalMCParts const& mcParticles, MCCentGen const& Cent, MCMultGen const& MCMult, MCIPGen const& IP, evtType const& eventType) + { + for (auto const& mcPart : mcParticles) { + + if (std::abs(mcPart.pdgCode()) != GenCuts.pdgTruthMother || std::abs(mcPart.y()) >= GenCuts.cfgRapidityCutGen) + continue; + std::vector daughterPDGs; + if (mcPart.has_daughters()) { + auto daughter01 = mcParticles.rawIteratorAt(mcPart.daughtersIds()[0] - mcParticles.offset()); + auto daughter02 = mcParticles.rawIteratorAt(mcPart.daughtersIds()[1] - mcParticles.offset()); + daughterPDGs = {daughter01.pdgCode(), daughter02.pdgCode()}; + } else { + daughterPDGs = {-1, -1}; + } + + if (GenCuts.isDaughterCheck) { + bool pass1 = std::abs(daughterPDGs[0]) == GenCuts.pdgTruthDaughter1 || std::abs(daughterPDGs[1]) == GenCuts.pdgTruthDaughter1; + bool pass2 = std::abs(daughterPDGs[0]) == GenCuts.pdgTruthDaughter2 || std::abs(daughterPDGs[1]) == GenCuts.pdgTruthDaughter2; + if (!pass1 || !pass2) + continue; + } + if (mcPart.pdgCode() > 0) // Consider INELt0 or INEL + qaRegistry.fill(HIST("EventGen/h5ResonanceTruth"), eventType, mcPart.pt(), Cent, MCMult, IP); + else + qaRegistry.fill(HIST("EventGen/h5ResonanceTruthAnti"), eventType, mcPart.pt(), Cent, MCMult, IP); + + daughterPDGs.clear(); + } + } + template - void fillMCCollision(MCCol const& mccol, MCPart const& mcparts, float impactpar = -999.0) + void fillMCCollision(MCCol const& mccol, MCPart const& mcparts, float impactpar = -999.0, float mult = -1.0) { auto centrality = 0.0; if constexpr (!isRun2) centrality = centEst(mccol); else centrality = mccol.centRun2V0M(); - bool inVtx10 = (std::abs(mccol.mcCollision().posZ()) > 10.) ? false : true; + + // bool inVtx10 = (std::abs(mccol.mcCollision().posZ()) > 10.) ? false : true; -> Gen. level informations will be processed in processMCGen + bool inVtx10 = (std::abs(mccol.posZ()) > 10.) ? false : true; bool isTrueINELgt0 = isTrueINEL0(mccol, mcparts); bool isTriggerTVX = mccol.selection_bit(aod::evsel::kIsTriggerTVX); bool isSel8 = mccol.sel8(); bool isSelected = colCuts.isSelected(mccol); - resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, impactpar); + resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, impactpar, mult); // QA for Trigger efficiency qaRegistry.fill(HIST("Event/hMCEventIndices"), centrality, aod::resocollision::kINEL); @@ -1239,6 +1346,22 @@ struct ResonanceInitializer { qaRegistry.add("hGoodMCCascIndices", "hGoodMCCascIndices", kTH1F, {idxAxis}); qaRegistry.add("Phi", "#phi distribution", kTH1F, {{65, -0.1, 6.4}}); } + + TString hNEventsMCLabels[4] = {"All", "z vrtx", "INEL", "INEL>0"}; + if (doprocessMCgen) { + AxisSpec centAxisGen = {binsCent, "Centrality (%)"}; + qaRegistry.add("EventGen/hNEventsMC", "EventGen/hNEventsMC", kTH1D, {{4, 0.0f, 4.0f}}); + for (int n = 1; n <= qaRegistry.get(HIST("EventGen/hNEventsMC"))->GetNbinsX(); n++) { + qaRegistry.get(HIST("EventGen/hNEventsMC"))->GetXaxis()->SetBinLabel(n, hNEventsMCLabels[n - 1]); + } + qaRegistry.add("EventGen/h5ResonanceTruth", "EventGen/h5ResonanceTruth", kTHnSparseD, {{2, 0.0f, 2.0f}, GenCuts.ptAxisGen, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/h5ResonanceTruthAnti", "EventGen/h5ResonanceTruthAnti", kTHnSparseD, {{2, 0.0f, 2.0f}, GenCuts.ptAxisGen, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/hZCollisionGen", "EventGen/hZCollisionGen", kTH1D, {{100, -20.0f, 20.0f}}); + + qaRegistry.add("EventGen/h4MultCent_genMC", "EventGen/h4MultCent_genMC", kTHnSparseD, {{2, 0.0f, 2.0f}, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/h4MultCent_recMC", "EventGen/h4MultCent_recMC", kTHnSparseD, {{2, 0.0f, 2.0f}, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/h2CentralityVsMultMC", "EventGen/h2CentralityVsMultMC", kTH2D, {centAxisGen, GenCuts.multNTracksAxis}); + } } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) // Simple copy from LambdaKzeroFinder.cxx @@ -1303,7 +1426,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1329,7 +1452,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1356,7 +1479,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1383,7 +1506,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1409,7 +1532,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1424,7 +1547,7 @@ struct ResonanceInitializer { } PROCESS_SWITCH(ResonanceInitializer, processTrackV0DataRun2, "Process for data", false); - void processTrackV0CascData(ResoEvents::iterator const& collision, + void processTrackV0CascData(ResoEvents001::iterator const& collision, soa::Filtered const& tracks, ResoV0s const& V0s, ResoCascades const& Cascades, @@ -1438,8 +1561,11 @@ struct ResonanceInitializer { if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) return; colCuts.fillQA(collision); + bool isRecINELgt0 = 0; + if (checkIsRecINELgt0) + isRecINELgt0 = collision.isInelGt0(); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(collision.multNTracksPV(), collision.multNTracksPVeta1(), collision.multNTracksPVetaHalf(), collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, isRecINELgt0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1466,7 +1592,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1493,7 +1619,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1525,7 +1651,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1552,7 +1678,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1582,7 +1708,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1611,7 +1737,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1632,25 +1758,47 @@ struct ResonanceInitializer { } PROCESS_SWITCH(ResonanceInitializer, processTrackV0MCRun2, "Process for MC", false); - void processTrackV0CascMC(soa::Join::iterator const& collision, - aod::McCollisions const&, soa::Filtered const& tracks, + void processTrackV0CascMC(soa::Join::iterator const& collision, + soa::Join const& mcCollisions, soa::Filtered const& tracks, ResoV0sMC const& V0s, ResoCascadesMC const& Cascades, aod::McParticles const& mcParticles, aod::BCsWithTimestamps const&) { + auto bc = collision.bc_as(); /// adding timestamp to access magnetic field later initCCDB(bc); if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + float Cent = 100.5f; + + const auto mcId = collision.mcCollisionId(); + auto mcCollision = mcCollisions.iteratorAt(mcId); + float impactpar = mcCollision.impactParameter(); + float mult = -1.0f; + if (GenCuts.cfgGenMultPercentile) + Cent = mcCollision.centFT0M(); + else + Cent = centEst(collision); + + bool isRecINELgt0 = 0; + if (checkIsRecINELgt0) + isRecINELgt0 = collision.isInelGt0(); + + resoCollisions(collision.multNTracksPV(), collision.multNTracksPVeta1(), collision.multNTracksPVetaHalf(), collision.posX(), collision.posY(), collision.posZ(), Cent, dBz, isRecINELgt0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } resoSpheroCollisions(computeSpherocity(tracks, trackSphMin, trackSphDef)); resoEvtPlCollisions(0, 0, 0, 0); - fillMCCollision(collision, mcParticles); + + if (GenCuts.cfgGenMult05) + mult = mcCollision.multMCNParticlesEta05(); + else if (GenCuts.cfgGenMult10) + mult = mcCollision.multMCNParticlesEta10(); + + fillMCCollision(collision, mcParticles, impactpar, mult); // Loop over tracks fillTracks(collision, tracks); @@ -1661,13 +1809,98 @@ struct ResonanceInitializer { fillCascades(collision, Cascades, tracks); // Loop over all MC particles - auto mcParts = selectedMCParticles->sliceBy(perMcCollision, collision.mcCollision().globalIndex()); + auto mcParts = selectedMCParticles->sliceBy(perMcCollision, mcId); fillMCParticles(mcParts, mcParticles); } PROCESS_SWITCH(ResonanceInitializer, processTrackV0CascMC, "Process for MC", false); + // Following the discussions at the PAG meeting (https://indico.cern.ch/event/1583408/) + // we have introduced an auxiliary task that, when the resonanceInitializer.cxx is used, + // Only consider N_rec / N_gen i.e. not consider level of N_gen at least once + void processMCgen(soa::Join::iterator const& mcCollision, + aod::McParticles const& mcParticles, + const soa::SmallGroups>& collisions, + aod::BCsWithTimestamps const&) + { + auto bc = mcCollision.bc_as(); /// adding timestamp to access magnetic field later + initCCDB(bc); + + auto getCentGen = [&]() { + if (cfgMultName.value == "FT0M") { // FT0A,C results wiill be updated later when CCDB is available + return mcCollision.centFT0M(); + } + return 100.5f; + }; + auto getCentReco = [&](auto const& col) { + if (cfgMultName.value == "FT0M") { + return col.centFT0M(); + } else if (cfgMultName.value == "FT0C") { + return col.centFT0C(); + } else if (cfgMultName.value == "FT0A") { + return col.centFT0A(); + } + return 100.5f; + }; + + float cent = getCentGen(); + float IP = mcCollision.impactParameter(); + float mult = -1; + if (GenCuts.cfgGenMult05) { + mult = mcCollision.multMCNParticlesEta05(); + } else if (GenCuts.cfgGenMult10) { + mult = mcCollision.multMCNParticlesEta10(); + } + + qaRegistry.fill(HIST("EventGen/hNEventsMC"), 0.5); + + if (GenCuts.isZvtxcutGen && std::fabs(mcCollision.posZ()) > GenCuts.cutzvertexGen) { + return; + } + qaRegistry.fill(HIST("EventGen/hZCollisionGen"), mcCollision.posZ()); + qaRegistry.fill(HIST("EventGen/hNEventsMC"), 1.5); + + int evType = 0; + + qaRegistry.fill(HIST("EventGen/hNEventsMC"), 2.5); + if (GenCuts.checkIsTrueINELgt0 && mcCollision.isInelGt0()) { + evType++; + qaRegistry.fill(HIST("EventGen/hNEventsMC"), 3.5); + } + + bool atLeastOne = false; + int biggestNContribs = -1; + + float centReco = 100.5f; + for (const auto& collision : collisions) { + if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) + continue; + if (!colCuts.isSelected(collision)) // Bug is appeared in colCuts-> double counting in event QA histo, will be fixed later + continue; + if (biggestNContribs < collision.multPVTotalContributors()) { + biggestNContribs = collision.multPVTotalContributors(); + centReco = getCentReco(collision); + } + + atLeastOne = true; + } + + if (GenCuts.cfgGenMultPercentile) { + fillMCGenParticles(mcParticles, cent, mult, IP, evType); + qaRegistry.fill(HIST("EventGen/h4MultCent_genMC"), evType, cent, mult, IP); + } else { + fillMCGenParticles(mcParticles, centReco, mult, IP, evType); + qaRegistry.fill(HIST("EventGen/h4MultCent_genMC"), evType, centReco, mult, IP); + qaRegistry.fill(HIST("EventGen/h2CentralityVsMultMC"), centReco, mult); + } + + if (atLeastOne) { + qaRegistry.fill(HIST("EventGen/h4MultCent_recMC"), evType, centReco, mult, IP); + } + } + PROCESS_SWITCH(ResonanceInitializer, processMCgen, "Process for MCGen", true); + void processTrackV0CascMCRun2(soa::Join::iterator const& collision, - aod::McCollisions const&, soa::Filtered const& tracks, + aod::McCollisions const&, ResoTracksMC const& tracks, ResoV0sMC const& V0s, ResoCascadesMC const& Cascades, aod::McParticles const& mcParticles, BCsWithRun2Info const&) @@ -1675,7 +1908,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } diff --git a/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx b/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx index a1ade00c627..5a4bc9f2be1 100644 --- a/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx +++ b/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx @@ -475,7 +475,7 @@ struct ResonanceModuleInitializer { bool isTriggerTVX = mccol.selection_bit(aod::evsel::kIsTriggerTVX); bool isSel8 = mccol.sel8(); bool isSelected = colCuts.isSelected(mccol); - resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, mcCent); + resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, mcCent, -1.0f); // QA for Trigger efficiency qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINEL); @@ -546,7 +546,7 @@ struct ResonanceModuleInitializer { colCuts.fillQA(collision); centrality = centEst(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -569,7 +569,7 @@ struct ResonanceModuleInitializer { colCuts.fillQARun2(collision); centrality = collision.centRun2V0M(); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1123,6 +1123,8 @@ struct ResonanceDaughterInitializer { reso2mcv0s(v0mc.pdgCode(), mothers[0], motherPDGs[0], + 0, + 0, daughters[0], daughters[1], daughterPDGs[0], @@ -1133,6 +1135,8 @@ struct ResonanceDaughterInitializer { reso2mcv0s(0, mothers[0], motherPDGs[0], + 0, + 0, daughters[0], daughters[1], daughterPDGs[0], @@ -1289,6 +1293,8 @@ struct ResonanceDaughterInitializer { reso2mccascades(cascmc.pdgCode(), mothers[0], motherPDGs[0], + 0, + 0, daughters[0], daughters[1], daughterPDGs[0], @@ -1299,6 +1305,8 @@ struct ResonanceDaughterInitializer { reso2mccascades(0, mothers[0], motherPDGs[0], + 0, + 0, daughters[0], daughters[1], daughterPDGs[0],