diff --git a/PWGLF/Tasks/Strangeness/cascadeAnalysisLightIonsDerivedData.cxx b/PWGLF/Tasks/Strangeness/cascadeAnalysisLightIonsDerivedData.cxx index 627a602a7d8..e0f2198d899 100644 --- a/PWGLF/Tasks/Strangeness/cascadeAnalysisLightIonsDerivedData.cxx +++ b/PWGLF/Tasks/Strangeness/cascadeAnalysisLightIonsDerivedData.cxx @@ -906,174 +906,188 @@ struct CascadeAnalysisLightIonsDerivedData { PROCESS_SWITCH(CascadeAnalysisLightIonsDerivedData, processData, "Process data", true); - void processMonteCarloRec(SimCollisions const& RecCols, CascadeMCCandidates const& fullCascades, DaughterTracks const&, CascadeMCCores const&) + void processMonteCarloRec(SimCollisions::iterator const& RecCol, CascadeMCCandidates const& fullCascades, DaughterTracks const&, CascadeMCCores const&) { - for (const auto& RecCol : RecCols) { - // Fill event counter before event selection - registryMC.fill(HIST("number_of_events_mc_rec"), 0); + // Fill event counter before event selection + registryMC.fill(HIST("number_of_events_mc_rec"), 0); - // Initialize CCDB objects using the BC info - initCCDB(RecCol); + // Initialize CCDB objects using the BC info + initCCDB(RecCol); - // Define the event centrality using different estimators - float centralityMcRec = -1; - float multiplicityMcRec = -1; + // Define the event centrality using different estimators + float centralityMcRec = -1; + float multiplicityMcRec = -1; - if (centralityEstimator == Option::kFT0C) { - centralityMcRec = RecCol.centFT0C(); - multiplicityMcRec = RecCol.multFT0C(); - } - if (centralityEstimator == Option::kFT0M) { - centralityMcRec = RecCol.centFT0M(); - multiplicityMcRec = RecCol.multFT0C() + RecCol.multFT0A(); - } - if (centralityEstimator == Option::kFV0A) { - centralityMcRec = RecCol.centFV0A(); - multiplicityMcRec = RecCol.multFV0A(); - } - if (centralityEstimator == Option::kNGlobal) { - centralityMcRec = RecCol.centNGlobal(); - multiplicityMcRec = RecCol.multNTracksGlobal(); - } + if (centralityEstimator == Option::kFT0C) { + centralityMcRec = RecCol.centFT0C(); + multiplicityMcRec = RecCol.multFT0C(); + } + if (centralityEstimator == Option::kFT0M) { + centralityMcRec = RecCol.centFT0M(); + multiplicityMcRec = RecCol.multFT0C() + RecCol.multFT0A(); + } + if (centralityEstimator == Option::kFV0A) { + centralityMcRec = RecCol.centFV0A(); + multiplicityMcRec = RecCol.multFV0A(); + } + if (centralityEstimator == Option::kNGlobal) { + centralityMcRec = RecCol.centNGlobal(); + multiplicityMcRec = RecCol.multNTracksGlobal(); + } - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 0, centralityMcRec); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 0, centralityMcRec); - // event selections - if (applySel8 && !RecCol.sel8()) - continue; - registryMC.fill(HIST("number_of_events_mc_rec"), 1); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 1, centralityMcRec); + // event selections + if (applySel8 && !RecCol.sel8()) + return; + registryMC.fill(HIST("number_of_events_mc_rec"), 1); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 1, centralityMcRec); - if (applyVtxZ && std::fabs(RecCol.posZ()) > zVtx) - continue; - registryMC.fill(HIST("number_of_events_mc_rec"), 2); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 2, centralityMcRec); + if (applyVtxZ && std::fabs(RecCol.posZ()) > zVtx) + return; + registryMC.fill(HIST("number_of_events_mc_rec"), 2); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 2, centralityMcRec); - if (rejectITSROFBorder && !RecCol.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { - continue; - } - registryMC.fill(HIST("number_of_events_mc_rec"), 3 /* Not at ITS ROF border */); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 3, centralityMcRec); + if (rejectITSROFBorder && !RecCol.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + return; + } + registryMC.fill(HIST("number_of_events_mc_rec"), 3 /* Not at ITS ROF border */); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 3, centralityMcRec); - if (rejectTFBorder && !RecCol.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { - continue; - } - registryMC.fill(HIST("number_of_events_mc_rec"), 4 /* Not at TF border */); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 4, centralityMcRec); + if (rejectTFBorder && !RecCol.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { + return; + } + registryMC.fill(HIST("number_of_events_mc_rec"), 4 /* Not at TF border */); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 4, centralityMcRec); - if (requireVertexITSTPC && !RecCol.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { - continue; - } - registryMC.fill(HIST("number_of_events_mc_rec"), 5 /* Contains at least one ITS-TPC track */); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 5, centralityMcRec); + if (requireVertexITSTPC && !RecCol.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + return; + } + registryMC.fill(HIST("number_of_events_mc_rec"), 5 /* Contains at least one ITS-TPC track */); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 5, centralityMcRec); - if (requireIsGoodZvtxFT0VsPV && !RecCol.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { - continue; - } - registryMC.fill(HIST("number_of_events_mc_rec"), 6 /* PV position consistency check */); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 6, centralityMcRec); + if (requireIsGoodZvtxFT0VsPV && !RecCol.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + return; + } + registryMC.fill(HIST("number_of_events_mc_rec"), 6 /* PV position consistency check */); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 6, centralityMcRec); - if (requireIsVertexTOFmatched && !RecCol.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) { - continue; - } - registryMC.fill(HIST("number_of_events_mc_rec"), 7 /* PV with at least one contributor matched with TOF */); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 7, centralityMcRec); + if (requireIsVertexTOFmatched && !RecCol.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) { + return; + } + registryMC.fill(HIST("number_of_events_mc_rec"), 7 /* PV with at least one contributor matched with TOF */); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 7, centralityMcRec); - if (requireIsVertexTRDmatched && !RecCol.selection_bit(o2::aod::evsel::kIsVertexTRDmatched)) { - continue; - } - registryMC.fill(HIST("number_of_events_mc_rec"), 8 /* PV with at least one contributor matched with TRD */); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 8, centralityMcRec); + if (requireIsVertexTRDmatched && !RecCol.selection_bit(o2::aod::evsel::kIsVertexTRDmatched)) { + return; + } + registryMC.fill(HIST("number_of_events_mc_rec"), 8 /* PV with at least one contributor matched with TRD */); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 8, centralityMcRec); - if (rejectSameBunchPileup && !RecCol.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { - continue; - } - registryMC.fill(HIST("number_of_events_mc_rec"), 9 /* Not at same bunch pile-up */); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 9, centralityMcRec); + if (rejectSameBunchPileup && !RecCol.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + return; + } + registryMC.fill(HIST("number_of_events_mc_rec"), 9 /* Not at same bunch pile-up */); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 9, centralityMcRec); - if (requireInel0 && RecCol.multNTracksPVeta1() < 1) { - continue; - } - registryMC.fill(HIST("number_of_events_mc_rec"), 10 /* INEL > 0 */); - registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 10, centralityMcRec); + if (requireInel0 && RecCol.multNTracksPVeta1() < 1) { + return; + } + registryMC.fill(HIST("number_of_events_mc_rec"), 10 /* INEL > 0 */); + registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 10, centralityMcRec); - // Store the Zvtx - registryQC.fill(HIST("hVertexZRec"), RecCol.posZ()); + // Store the Zvtx + registryQC.fill(HIST("hVertexZRec"), RecCol.posZ()); + + // Store the event centrality + registryMC.fill(HIST("hCentEstimator_truerec"), centralityMcRec); + registryMC.fill(HIST("hCentralityVsNch_truerec"), centralityMcRec, RecCol.multNTracksPVeta1()); + registryMC.fill(HIST("hCentralityVsMultiplicity_truerec"), centralityMcRec, multiplicityMcRec); - // Store the event centrality - registryMC.fill(HIST("hCentEstimator_truerec"), centralityMcRec); - registryMC.fill(HIST("hCentralityVsNch_truerec"), centralityMcRec, RecCol.multNTracksPVeta1()); - registryMC.fill(HIST("hCentralityVsMultiplicity_truerec"), centralityMcRec, multiplicityMcRec); + for (const auto& casc : fullCascades) { + if (etaMin > casc.bacheloreta() || casc.bacheloreta() > etaMax || + etaMin > casc.negativeeta() || casc.negativeeta() > etaMax || + etaMin > casc.positiveeta() || casc.positiveeta() > etaMax) + continue; // remove acceptance that's badly reproduced by MC / superfluous in future - for (const auto& casc : fullCascades) { - if (etaMin > casc.bacheloreta() || casc.bacheloreta() > etaMax || - etaMin > casc.negativeeta() || casc.negativeeta() > etaMax || - etaMin > casc.positiveeta() || casc.positiveeta() > etaMax) - continue; // remove acceptance that's badly reproduced by MC / superfluous in future + if (!casc.has_cascMCCore()) + continue; - if (!casc.has_cascMCCore()) - continue; + auto cascMC = casc.template cascMCCore_as(); - auto cascMC = casc.template cascMCCore_as(); + auto bach = casc.bachTrackExtra_as(); + auto pos = casc.posTrackExtra_as(); + auto neg = casc.negTrackExtra_as(); - auto bach = casc.bachTrackExtra_as(); - auto pos = casc.posTrackExtra_as(); - auto neg = casc.negTrackExtra_as(); + int pdgParent = cascMC.pdgCode(); + bool isPhysPrim = cascMC.isPhysicalPrimary(); + if (pdgParent == 0) + continue; + if (!isPhysPrim) + continue; - int pdgParent = cascMC.pdgCode(); - bool isPhysPrim = cascMC.isPhysicalPrimary(); - if (pdgParent == 0) - continue; - if (!isPhysPrim) - continue; + float ptmc = RecoDecay::sqrtSumOfSquares(cascMC.pxMC(), cascMC.pyMC()); - float ptmc = RecoDecay::sqrtSumOfSquares(cascMC.pxMC(), cascMC.pyMC()); - - // ------------------------------------- Store selctions distribution for QC - registryQC.fill(HIST("hv0cosPARec"), casc.v0cosPA(RecCol.posX(), RecCol.posY(), RecCol.posZ())); - registryQC.fill(HIST("hcasccosPARec"), casc.casccosPA(RecCol.posX(), RecCol.posY(), RecCol.posZ())); - registryQC.fill(HIST("hv0radiusRec"), casc.v0radius()); - registryQC.fill(HIST("hcascradiusRec"), casc.cascradius()); - registryQC.fill(HIST("hdcaV0daughtersRec"), casc.dcaV0daughters()); - registryQC.fill(HIST("hdcacascdaughtersRec"), casc.dcacascdaughters()); - registryQC.fill(HIST("hdcapostopvRec"), casc.dcapostopv()); - registryQC.fill(HIST("hdcanegtopvRec"), casc.dcanegtopv()); - registryQC.fill(HIST("hdcabachtopvRec"), casc.dcabachtopv()); - registryQC.fill(HIST("hdcav0topvRec"), casc.dcav0topv(RecCol.posX(), RecCol.posY(), RecCol.posZ())); - - // ------------------------------------- Store selctions distribution for analysis - if (casc.sign() < 0) { - if (pdgParent == kXiMinus) { - registryMC.fill(HIST("hMassXineg_truerec"), centralityMcRec, ptmc, casc.mXi()); - } - if (pdgParent == kOmegaMinus) { - registryMC.fill(HIST("hMassOmeganeg_truerec"), centralityMcRec, ptmc, casc.mOmega()); - } + // ------------------------------------- Store selctions distribution for QC + registryQC.fill(HIST("hv0cosPARec"), casc.v0cosPA(RecCol.posX(), RecCol.posY(), RecCol.posZ())); + registryQC.fill(HIST("hcasccosPARec"), casc.casccosPA(RecCol.posX(), RecCol.posY(), RecCol.posZ())); + registryQC.fill(HIST("hv0radiusRec"), casc.v0radius()); + registryQC.fill(HIST("hcascradiusRec"), casc.cascradius()); + registryQC.fill(HIST("hdcaV0daughtersRec"), casc.dcaV0daughters()); + registryQC.fill(HIST("hdcacascdaughtersRec"), casc.dcacascdaughters()); + registryQC.fill(HIST("hdcapostopvRec"), casc.dcapostopv()); + registryQC.fill(HIST("hdcanegtopvRec"), casc.dcanegtopv()); + registryQC.fill(HIST("hdcabachtopvRec"), casc.dcabachtopv()); + registryQC.fill(HIST("hdcav0topvRec"), casc.dcav0topv(RecCol.posX(), RecCol.posY(), RecCol.posZ())); + + // ------------------------------------- Store selctions distribution for analysis + if (casc.sign() < 0) { + if (pdgParent == kXiMinus) { + registryMC.fill(HIST("hMassXineg_truerec"), centralityMcRec, ptmc, casc.mXi()); + } + if (pdgParent == kOmegaMinus) { + registryMC.fill(HIST("hMassOmeganeg_truerec"), centralityMcRec, ptmc, casc.mOmega()); } + } - if (casc.sign() > 0) { - if (pdgParent == kXiPlusBar) { - registryMC.fill(HIST("hMassXipos_truerec"), centralityMcRec, ptmc, casc.mXi()); - } - if (pdgParent == kOmegaPlusBar) { - registryMC.fill(HIST("hMassOmegapos_truerec"), centralityMcRec, ptmc, casc.mOmega()); - } + if (casc.sign() > 0) { + if (pdgParent == kXiPlusBar) { + registryMC.fill(HIST("hMassXipos_truerec"), centralityMcRec, ptmc, casc.mXi()); + } + if (pdgParent == kOmegaPlusBar) { + registryMC.fill(HIST("hMassOmegapos_truerec"), centralityMcRec, ptmc, casc.mOmega()); } + } - if (casc.sign() < 0 && pdgParent == kXiMinus && passedXiSelection(casc, pos, neg, bach, RecCol)) { + // compute MC association + const bool isXiMC = (std::abs(pdgParent) == PDG_t::kXiMinus); + const bool isOmegaMC = (std::abs(pdgParent) == PDG_t::kOmegaMinus); + bool isTrueMCCascade = false; + bool isTrueMCCascadeDecay = false; + bool isCorrectLambdaDecay = false; + + if (isPhysPrim && (isXiMC || isOmegaMC)) + isTrueMCCascade = true; + if (isTrueMCCascade && ((casc.sign() > 0 && cascMC.pdgCodePositive() == PDG_t::kPiPlus && cascMC.pdgCodeNegative() == PDG_t::kProtonBar) || (casc.sign() < 0 && cascMC.pdgCodePositive() == PDG_t::kProton && cascMC.pdgCodeNegative() == PDG_t::kPiMinus))) + isCorrectLambdaDecay = true; + if (isTrueMCCascade && isCorrectLambdaDecay && ((isXiMC && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kPiPlus) || (isOmegaMC && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kKPlus))) + isTrueMCCascadeDecay = true; + + if (isTrueMCCascadeDecay) { + if (pdgParent == kXiMinus && passedXiSelection(casc, pos, neg, bach, RecCol)) { registryMC.fill(HIST("hMassXinegSelected_truerec"), centralityMcRec, ptmc, casc.mXi()); } - if (casc.sign() < 0 && pdgParent == kOmegaMinus && passedOmegaSelection(casc, pos, neg, bach, RecCol)) { + if (pdgParent == kOmegaMinus && passedOmegaSelection(casc, pos, neg, bach, RecCol)) { registryMC.fill(HIST("hMassOmeganegSelected_truerec"), centralityMcRec, ptmc, casc.mOmega()); } - if (casc.sign() > 0 && pdgParent == kXiPlusBar && passedXiSelection(casc, pos, neg, bach, RecCol)) { + if (pdgParent == kXiPlusBar && passedXiSelection(casc, pos, neg, bach, RecCol)) { registryMC.fill(HIST("hMassXiposSelected_truerec"), centralityMcRec, ptmc, casc.mXi()); } - if (casc.sign() > 0 && pdgParent == kOmegaPlusBar && passedOmegaSelection(casc, pos, neg, bach, RecCol)) { + if (pdgParent == kOmegaPlusBar && passedOmegaSelection(casc, pos, neg, bach, RecCol)) { registryMC.fill(HIST("hMassOmegaposSelected_truerec"), centralityMcRec, ptmc, casc.mOmega()); } - } // casc loop - } // rec.collision loop + } + } // casc loop } PROCESS_SWITCH(CascadeAnalysisLightIonsDerivedData, processMonteCarloRec, "Process MC Rec", false);