Skip to content

Commit 4ae5e28

Browse files
spucilloalibuildromainschotter
authored
[PWGLF] Fix MC association and processMonteCarloRec function in cascadeAnalysisLightIonsDerivedData.cxx (#14843)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch> Co-authored-by: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com>
1 parent 6c3659b commit 4ae5e28

File tree

1 file changed

+147
-133
lines changed

1 file changed

+147
-133
lines changed

‎PWGLF/Tasks/Strangeness/cascadeAnalysisLightIonsDerivedData.cxx‎

Lines changed: 147 additions & 133 deletions
Original file line numberDiff line numberDiff line change
@@ -906,174 +906,188 @@ struct CascadeAnalysisLightIonsDerivedData {
906906

907907
PROCESS_SWITCH(CascadeAnalysisLightIonsDerivedData, processData, "Process data", true);
908908

909-
void processMonteCarloRec(SimCollisions const& RecCols, CascadeMCCandidates const& fullCascades, DaughterTracks const&, CascadeMCCores const&)
909+
void processMonteCarloRec(SimCollisions::iterator const& RecCol, CascadeMCCandidates const& fullCascades, DaughterTracks const&, CascadeMCCores const&)
910910
{
911-
for (const auto& RecCol : RecCols) {
912-
// Fill event counter before event selection
913-
registryMC.fill(HIST("number_of_events_mc_rec"), 0);
911+
// Fill event counter before event selection
912+
registryMC.fill(HIST("number_of_events_mc_rec"), 0);
914913

915-
// Initialize CCDB objects using the BC info
916-
initCCDB(RecCol);
914+
// Initialize CCDB objects using the BC info
915+
initCCDB(RecCol);
917916

918-
// Define the event centrality using different estimators
919-
float centralityMcRec = -1;
920-
float multiplicityMcRec = -1;
917+
// Define the event centrality using different estimators
918+
float centralityMcRec = -1;
919+
float multiplicityMcRec = -1;
921920

922-
if (centralityEstimator == Option::kFT0C) {
923-
centralityMcRec = RecCol.centFT0C();
924-
multiplicityMcRec = RecCol.multFT0C();
925-
}
926-
if (centralityEstimator == Option::kFT0M) {
927-
centralityMcRec = RecCol.centFT0M();
928-
multiplicityMcRec = RecCol.multFT0C() + RecCol.multFT0A();
929-
}
930-
if (centralityEstimator == Option::kFV0A) {
931-
centralityMcRec = RecCol.centFV0A();
932-
multiplicityMcRec = RecCol.multFV0A();
933-
}
934-
if (centralityEstimator == Option::kNGlobal) {
935-
centralityMcRec = RecCol.centNGlobal();
936-
multiplicityMcRec = RecCol.multNTracksGlobal();
937-
}
921+
if (centralityEstimator == Option::kFT0C) {
922+
centralityMcRec = RecCol.centFT0C();
923+
multiplicityMcRec = RecCol.multFT0C();
924+
}
925+
if (centralityEstimator == Option::kFT0M) {
926+
centralityMcRec = RecCol.centFT0M();
927+
multiplicityMcRec = RecCol.multFT0C() + RecCol.multFT0A();
928+
}
929+
if (centralityEstimator == Option::kFV0A) {
930+
centralityMcRec = RecCol.centFV0A();
931+
multiplicityMcRec = RecCol.multFV0A();
932+
}
933+
if (centralityEstimator == Option::kNGlobal) {
934+
centralityMcRec = RecCol.centNGlobal();
935+
multiplicityMcRec = RecCol.multNTracksGlobal();
936+
}
938937

939-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 0, centralityMcRec);
938+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 0, centralityMcRec);
940939

941-
// event selections
942-
if (applySel8 && !RecCol.sel8())
943-
continue;
944-
registryMC.fill(HIST("number_of_events_mc_rec"), 1);
945-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 1, centralityMcRec);
940+
// event selections
941+
if (applySel8 && !RecCol.sel8())
942+
return;
943+
registryMC.fill(HIST("number_of_events_mc_rec"), 1);
944+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 1, centralityMcRec);
946945

947-
if (applyVtxZ && std::fabs(RecCol.posZ()) > zVtx)
948-
continue;
949-
registryMC.fill(HIST("number_of_events_mc_rec"), 2);
950-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 2, centralityMcRec);
946+
if (applyVtxZ && std::fabs(RecCol.posZ()) > zVtx)
947+
return;
948+
registryMC.fill(HIST("number_of_events_mc_rec"), 2);
949+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 2, centralityMcRec);
951950

952-
if (rejectITSROFBorder && !RecCol.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
953-
continue;
954-
}
955-
registryMC.fill(HIST("number_of_events_mc_rec"), 3 /* Not at ITS ROF border */);
956-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 3, centralityMcRec);
951+
if (rejectITSROFBorder && !RecCol.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
952+
return;
953+
}
954+
registryMC.fill(HIST("number_of_events_mc_rec"), 3 /* Not at ITS ROF border */);
955+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 3, centralityMcRec);
957956

958-
if (rejectTFBorder && !RecCol.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) {
959-
continue;
960-
}
961-
registryMC.fill(HIST("number_of_events_mc_rec"), 4 /* Not at TF border */);
962-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 4, centralityMcRec);
957+
if (rejectTFBorder && !RecCol.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) {
958+
return;
959+
}
960+
registryMC.fill(HIST("number_of_events_mc_rec"), 4 /* Not at TF border */);
961+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 4, centralityMcRec);
963962

964-
if (requireVertexITSTPC && !RecCol.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) {
965-
continue;
966-
}
967-
registryMC.fill(HIST("number_of_events_mc_rec"), 5 /* Contains at least one ITS-TPC track */);
968-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 5, centralityMcRec);
963+
if (requireVertexITSTPC && !RecCol.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) {
964+
return;
965+
}
966+
registryMC.fill(HIST("number_of_events_mc_rec"), 5 /* Contains at least one ITS-TPC track */);
967+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 5, centralityMcRec);
969968

970-
if (requireIsGoodZvtxFT0VsPV && !RecCol.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) {
971-
continue;
972-
}
973-
registryMC.fill(HIST("number_of_events_mc_rec"), 6 /* PV position consistency check */);
974-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 6, centralityMcRec);
969+
if (requireIsGoodZvtxFT0VsPV && !RecCol.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) {
970+
return;
971+
}
972+
registryMC.fill(HIST("number_of_events_mc_rec"), 6 /* PV position consistency check */);
973+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 6, centralityMcRec);
975974

976-
if (requireIsVertexTOFmatched && !RecCol.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) {
977-
continue;
978-
}
979-
registryMC.fill(HIST("number_of_events_mc_rec"), 7 /* PV with at least one contributor matched with TOF */);
980-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 7, centralityMcRec);
975+
if (requireIsVertexTOFmatched && !RecCol.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) {
976+
return;
977+
}
978+
registryMC.fill(HIST("number_of_events_mc_rec"), 7 /* PV with at least one contributor matched with TOF */);
979+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 7, centralityMcRec);
981980

982-
if (requireIsVertexTRDmatched && !RecCol.selection_bit(o2::aod::evsel::kIsVertexTRDmatched)) {
983-
continue;
984-
}
985-
registryMC.fill(HIST("number_of_events_mc_rec"), 8 /* PV with at least one contributor matched with TRD */);
986-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 8, centralityMcRec);
981+
if (requireIsVertexTRDmatched && !RecCol.selection_bit(o2::aod::evsel::kIsVertexTRDmatched)) {
982+
return;
983+
}
984+
registryMC.fill(HIST("number_of_events_mc_rec"), 8 /* PV with at least one contributor matched with TRD */);
985+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 8, centralityMcRec);
987986

988-
if (rejectSameBunchPileup && !RecCol.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
989-
continue;
990-
}
991-
registryMC.fill(HIST("number_of_events_mc_rec"), 9 /* Not at same bunch pile-up */);
992-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 9, centralityMcRec);
987+
if (rejectSameBunchPileup && !RecCol.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
988+
return;
989+
}
990+
registryMC.fill(HIST("number_of_events_mc_rec"), 9 /* Not at same bunch pile-up */);
991+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 9, centralityMcRec);
993992

994-
if (requireInel0 && RecCol.multNTracksPVeta1() < 1) {
995-
continue;
996-
}
997-
registryMC.fill(HIST("number_of_events_mc_rec"), 10 /* INEL > 0 */);
998-
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 10, centralityMcRec);
993+
if (requireInel0 && RecCol.multNTracksPVeta1() < 1) {
994+
return;
995+
}
996+
registryMC.fill(HIST("number_of_events_mc_rec"), 10 /* INEL > 0 */);
997+
registryMC.fill(HIST("number_of_events_mc_rec_vs_centrality"), 10, centralityMcRec);
999998

1000-
// Store the Zvtx
1001-
registryQC.fill(HIST("hVertexZRec"), RecCol.posZ());
999+
// Store the Zvtx
1000+
registryQC.fill(HIST("hVertexZRec"), RecCol.posZ());
1001+
1002+
// Store the event centrality
1003+
registryMC.fill(HIST("hCentEstimator_truerec"), centralityMcRec);
1004+
registryMC.fill(HIST("hCentralityVsNch_truerec"), centralityMcRec, RecCol.multNTracksPVeta1());
1005+
registryMC.fill(HIST("hCentralityVsMultiplicity_truerec"), centralityMcRec, multiplicityMcRec);
10021006

1003-
// Store the event centrality
1004-
registryMC.fill(HIST("hCentEstimator_truerec"), centralityMcRec);
1005-
registryMC.fill(HIST("hCentralityVsNch_truerec"), centralityMcRec, RecCol.multNTracksPVeta1());
1006-
registryMC.fill(HIST("hCentralityVsMultiplicity_truerec"), centralityMcRec, multiplicityMcRec);
1007+
for (const auto& casc : fullCascades) {
1008+
if (etaMin > casc.bacheloreta() || casc.bacheloreta() > etaMax ||
1009+
etaMin > casc.negativeeta() || casc.negativeeta() > etaMax ||
1010+
etaMin > casc.positiveeta() || casc.positiveeta() > etaMax)
1011+
continue; // remove acceptance that's badly reproduced by MC / superfluous in future
10071012

1008-
for (const auto& casc : fullCascades) {
1009-
if (etaMin > casc.bacheloreta() || casc.bacheloreta() > etaMax ||
1010-
etaMin > casc.negativeeta() || casc.negativeeta() > etaMax ||
1011-
etaMin > casc.positiveeta() || casc.positiveeta() > etaMax)
1012-
continue; // remove acceptance that's badly reproduced by MC / superfluous in future
1013+
if (!casc.has_cascMCCore())
1014+
continue;
10131015

1014-
if (!casc.has_cascMCCore())
1015-
continue;
1016+
auto cascMC = casc.template cascMCCore_as<CascadeMCCores>();
10161017

1017-
auto cascMC = casc.template cascMCCore_as<CascadeMCCores>();
1018+
auto bach = casc.bachTrackExtra_as<DaughterTracks>();
1019+
auto pos = casc.posTrackExtra_as<DaughterTracks>();
1020+
auto neg = casc.negTrackExtra_as<DaughterTracks>();
10181021

1019-
auto bach = casc.bachTrackExtra_as<DaughterTracks>();
1020-
auto pos = casc.posTrackExtra_as<DaughterTracks>();
1021-
auto neg = casc.negTrackExtra_as<DaughterTracks>();
1022+
int pdgParent = cascMC.pdgCode();
1023+
bool isPhysPrim = cascMC.isPhysicalPrimary();
1024+
if (pdgParent == 0)
1025+
continue;
1026+
if (!isPhysPrim)
1027+
continue;
10221028

1023-
int pdgParent = cascMC.pdgCode();
1024-
bool isPhysPrim = cascMC.isPhysicalPrimary();
1025-
if (pdgParent == 0)
1026-
continue;
1027-
if (!isPhysPrim)
1028-
continue;
1029+
float ptmc = RecoDecay::sqrtSumOfSquares(cascMC.pxMC(), cascMC.pyMC());
10291030

1030-
float ptmc = RecoDecay::sqrtSumOfSquares(cascMC.pxMC(), cascMC.pyMC());
1031-
1032-
// ------------------------------------- Store selctions distribution for QC
1033-
registryQC.fill(HIST("hv0cosPARec"), casc.v0cosPA(RecCol.posX(), RecCol.posY(), RecCol.posZ()));
1034-
registryQC.fill(HIST("hcasccosPARec"), casc.casccosPA(RecCol.posX(), RecCol.posY(), RecCol.posZ()));
1035-
registryQC.fill(HIST("hv0radiusRec"), casc.v0radius());
1036-
registryQC.fill(HIST("hcascradiusRec"), casc.cascradius());
1037-
registryQC.fill(HIST("hdcaV0daughtersRec"), casc.dcaV0daughters());
1038-
registryQC.fill(HIST("hdcacascdaughtersRec"), casc.dcacascdaughters());
1039-
registryQC.fill(HIST("hdcapostopvRec"), casc.dcapostopv());
1040-
registryQC.fill(HIST("hdcanegtopvRec"), casc.dcanegtopv());
1041-
registryQC.fill(HIST("hdcabachtopvRec"), casc.dcabachtopv());
1042-
registryQC.fill(HIST("hdcav0topvRec"), casc.dcav0topv(RecCol.posX(), RecCol.posY(), RecCol.posZ()));
1043-
1044-
// ------------------------------------- Store selctions distribution for analysis
1045-
if (casc.sign() < 0) {
1046-
if (pdgParent == kXiMinus) {
1047-
registryMC.fill(HIST("hMassXineg_truerec"), centralityMcRec, ptmc, casc.mXi());
1048-
}
1049-
if (pdgParent == kOmegaMinus) {
1050-
registryMC.fill(HIST("hMassOmeganeg_truerec"), centralityMcRec, ptmc, casc.mOmega());
1051-
}
1031+
// ------------------------------------- Store selctions distribution for QC
1032+
registryQC.fill(HIST("hv0cosPARec"), casc.v0cosPA(RecCol.posX(), RecCol.posY(), RecCol.posZ()));
1033+
registryQC.fill(HIST("hcasccosPARec"), casc.casccosPA(RecCol.posX(), RecCol.posY(), RecCol.posZ()));
1034+
registryQC.fill(HIST("hv0radiusRec"), casc.v0radius());
1035+
registryQC.fill(HIST("hcascradiusRec"), casc.cascradius());
1036+
registryQC.fill(HIST("hdcaV0daughtersRec"), casc.dcaV0daughters());
1037+
registryQC.fill(HIST("hdcacascdaughtersRec"), casc.dcacascdaughters());
1038+
registryQC.fill(HIST("hdcapostopvRec"), casc.dcapostopv());
1039+
registryQC.fill(HIST("hdcanegtopvRec"), casc.dcanegtopv());
1040+
registryQC.fill(HIST("hdcabachtopvRec"), casc.dcabachtopv());
1041+
registryQC.fill(HIST("hdcav0topvRec"), casc.dcav0topv(RecCol.posX(), RecCol.posY(), RecCol.posZ()));
1042+
1043+
// ------------------------------------- Store selctions distribution for analysis
1044+
if (casc.sign() < 0) {
1045+
if (pdgParent == kXiMinus) {
1046+
registryMC.fill(HIST("hMassXineg_truerec"), centralityMcRec, ptmc, casc.mXi());
1047+
}
1048+
if (pdgParent == kOmegaMinus) {
1049+
registryMC.fill(HIST("hMassOmeganeg_truerec"), centralityMcRec, ptmc, casc.mOmega());
10521050
}
1051+
}
10531052

1054-
if (casc.sign() > 0) {
1055-
if (pdgParent == kXiPlusBar) {
1056-
registryMC.fill(HIST("hMassXipos_truerec"), centralityMcRec, ptmc, casc.mXi());
1057-
}
1058-
if (pdgParent == kOmegaPlusBar) {
1059-
registryMC.fill(HIST("hMassOmegapos_truerec"), centralityMcRec, ptmc, casc.mOmega());
1060-
}
1053+
if (casc.sign() > 0) {
1054+
if (pdgParent == kXiPlusBar) {
1055+
registryMC.fill(HIST("hMassXipos_truerec"), centralityMcRec, ptmc, casc.mXi());
1056+
}
1057+
if (pdgParent == kOmegaPlusBar) {
1058+
registryMC.fill(HIST("hMassOmegapos_truerec"), centralityMcRec, ptmc, casc.mOmega());
10611059
}
1060+
}
10621061

1063-
if (casc.sign() < 0 && pdgParent == kXiMinus && passedXiSelection(casc, pos, neg, bach, RecCol)) {
1062+
// compute MC association
1063+
const bool isXiMC = (std::abs(pdgParent) == PDG_t::kXiMinus);
1064+
const bool isOmegaMC = (std::abs(pdgParent) == PDG_t::kOmegaMinus);
1065+
bool isTrueMCCascade = false;
1066+
bool isTrueMCCascadeDecay = false;
1067+
bool isCorrectLambdaDecay = false;
1068+
1069+
if (isPhysPrim && (isXiMC || isOmegaMC))
1070+
isTrueMCCascade = true;
1071+
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)))
1072+
isCorrectLambdaDecay = true;
1073+
if (isTrueMCCascade && isCorrectLambdaDecay && ((isXiMC && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kPiPlus) || (isOmegaMC && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kKPlus)))
1074+
isTrueMCCascadeDecay = true;
1075+
1076+
if (isTrueMCCascadeDecay) {
1077+
if (pdgParent == kXiMinus && passedXiSelection(casc, pos, neg, bach, RecCol)) {
10641078
registryMC.fill(HIST("hMassXinegSelected_truerec"), centralityMcRec, ptmc, casc.mXi());
10651079
}
1066-
if (casc.sign() < 0 && pdgParent == kOmegaMinus && passedOmegaSelection(casc, pos, neg, bach, RecCol)) {
1080+
if (pdgParent == kOmegaMinus && passedOmegaSelection(casc, pos, neg, bach, RecCol)) {
10671081
registryMC.fill(HIST("hMassOmeganegSelected_truerec"), centralityMcRec, ptmc, casc.mOmega());
10681082
}
1069-
if (casc.sign() > 0 && pdgParent == kXiPlusBar && passedXiSelection(casc, pos, neg, bach, RecCol)) {
1083+
if (pdgParent == kXiPlusBar && passedXiSelection(casc, pos, neg, bach, RecCol)) {
10701084
registryMC.fill(HIST("hMassXiposSelected_truerec"), centralityMcRec, ptmc, casc.mXi());
10711085
}
1072-
if (casc.sign() > 0 && pdgParent == kOmegaPlusBar && passedOmegaSelection(casc, pos, neg, bach, RecCol)) {
1086+
if (pdgParent == kOmegaPlusBar && passedOmegaSelection(casc, pos, neg, bach, RecCol)) {
10731087
registryMC.fill(HIST("hMassOmegaposSelected_truerec"), centralityMcRec, ptmc, casc.mOmega());
10741088
}
1075-
} // casc loop
1076-
} // rec.collision loop
1089+
}
1090+
} // casc loop
10771091
}
10781092

10791093
PROCESS_SWITCH(CascadeAnalysisLightIonsDerivedData, processMonteCarloRec, "Process MC Rec", false);

0 commit comments

Comments
 (0)