Skip to content

Commit 70158dd

Browse files
dsekihatdsekihat
authored
[PWGEM/Dilepton] update studyMCTruth.cxx (#10111)
1 parent b177f17 commit 70158dd

File tree

2 files changed

+118
-53
lines changed

2 files changed

+118
-53
lines changed

‎PWGEM/Dilepton/Tasks/associateMCcollision.cxx‎

Lines changed: 31 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,6 @@
1717
#include "Framework/runDataProcessing.h"
1818
#include "Framework/AnalysisTask.h"
1919
#include "Framework/ASoAHelpers.h"
20-
2120
#include "PWGEM/Dilepton/DataModel/dileptonTables.h"
2221

2322
using namespace o2;
@@ -33,24 +32,25 @@ struct associateMCcollision {
3332

3433
void init(InitContext&)
3534
{
35+
if (doprocessNcontrib && doprocessNcontrib_Derived) {
36+
LOGF(fatal, "Please select only 1 process function.");
37+
}
3638
addhistograms();
3739
}
3840

3941
~associateMCcollision() {}
4042

4143
void addhistograms()
4244
{
43-
fRegistry.add("hReccollsPerMCcoll", "Rec. colls per MC coll;Rec. colls per MC coll;Number of MC collisions", kTH1F, {{21, -0.5, 20.5}}, false);
45+
fRegistry.add("hReccollsPerMCcoll", "Rec. colls per MC coll;Rec. colls per MC coll;Number of MC collisions", kTH1D, {{21, -0.5, 20.5}}, false);
4446
}
4547

46-
using MyCollisions = soa::Join<aod::EMEvents, aod::EMMCEventLabels>;
47-
using MyCollision = MyCollisions::iterator;
48-
PresliceUnsorted<MyCollisions> recColperMcCollision = aod::emmceventlabel::emmceventId;
49-
50-
void processNcontrib(aod::EMMCEvents const& mccollisions, MyCollisions const& collisions)
48+
template <typename TMCCollisions, typename TCollisions, typename TPreslice>
49+
void runMC(TMCCollisions const& mcCollisions, TCollisions const& collisions, TPreslice const& perMCCollision)
5150
{
52-
for (auto& mccollision : mccollisions) {
53-
auto rec_colls_per_mccoll = collisions.sliceBy(recColperMcCollision, mccollision.globalIndex());
51+
52+
for (auto& mcCollision : mcCollisions) {
53+
auto rec_colls_per_mccoll = collisions.sliceBy(perMCCollision, mcCollision.globalIndex());
5454
fRegistry.fill(HIST("hReccollsPerMCcoll"), rec_colls_per_mccoll.size());
5555
uint32_t maxNumContrib = 0;
5656
int rec_col_globalIndex = -999;
@@ -63,8 +63,28 @@ struct associateMCcollision {
6363
// LOGF(info, "rec_col_globalIndex = %d", rec_col_globalIndex);
6464
mpemeventIds(rec_col_globalIndex);
6565
} // end of mc collision
66-
} // end of process
67-
PROCESS_SWITCH(associateMCcollision, processNcontrib, "produce most probable emeventId based on Ncontrib to PV", true);
66+
67+
} // end of runMC
68+
69+
using MyCollisions = soa::Join<aod::Collisions, aod::McCollisionLabels>;
70+
using MyCollision = MyCollisions::iterator;
71+
PresliceUnsorted<aod::McCollisionLabels> recColperMcCollision = aod::mccollisionlabel::mcCollisionId;
72+
73+
using MyEMCollisions = soa::Join<aod::EMEvents, aod::EMMCEventLabels>;
74+
using MyEMCollision = MyEMCollisions::iterator;
75+
PresliceUnsorted<aod::EMMCEventLabels> recColperMcCollision_derived = aod::emmceventlabel::emmceventId;
76+
77+
void processNcontrib_Derived(aod::EMMCEvents const& mcCollisions, MyEMCollisions const& collisions)
78+
{
79+
runMC(mcCollisions, collisions, recColperMcCollision_derived);
80+
}
81+
PROCESS_SWITCH(associateMCcollision, processNcontrib_Derived, "produce most probable emeventId based on Ncontrib to PV for derived AOD", true);
82+
83+
void processNcontrib(aod::McCollisions const& mcCollisions, MyCollisions const& collisions)
84+
{
85+
runMC(mcCollisions, collisions, recColperMcCollision);
86+
}
87+
PROCESS_SWITCH(associateMCcollision, processNcontrib, "produce most probable emeventId based on Ncontrib to PV for original AOD", false);
6888
};
6989

7090
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

‎PWGEM/Dilepton/Tasks/studyMCTruth.cxx‎

Lines changed: 87 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
#include "ReconstructionDataFormats/Track.h"
2626
#include "Common/Core/TableHelper.h"
2727
#include "Common/DataModel/EventSelection.h"
28+
#include "PWGEM/Dilepton/DataModel/dileptonTables.h"
2829
#include "PWGEM/Dilepton/Utils/MCUtilities.h"
2930

3031
using namespace o2;
@@ -74,29 +75,39 @@ struct studyMCTruth {
7475
}
7576

7677
static constexpr std::string_view dileptonSigns[3] = {"uls/", "lspp/", "lsmm/"};
77-
static constexpr std::string_view evNames[4] = {"before/", "after/"};
78+
static constexpr std::string_view evNames[4] = {"allMC/", "selectedMC/", "selectedMC_and_Rec/", "selectedMC_and_selectedRec/"};
7879

7980
void addHistograms()
8081
{
8182
const AxisSpec axis_mll{ConfMllBins, "m_{ll} (GeV/c^{2})"};
8283
const AxisSpec axis_ptll{ConfPtllBins, "p_{T,ll} (GeV/c)"};
8384

84-
fRegistry.add("Event/before/hZvtx", "MC Zvtx;Z_{vtx} (cm)", kTH1D, {{100, -50, +50}}, false);
85-
fRegistry.add("Event/before/hImpactParameter", "impact parameter;impact parameter b (fm)", kTH1D, {{200, 0, 20}}, false);
86-
fRegistry.addClone("Event/before/", "Event/after/");
87-
88-
fRegistry.add("Pair/before/Pi0/uls/hMvsPt", "m_{ll} vs. p_{T,ll}", kTH2D, {axis_mll, axis_ptll}, true);
89-
fRegistry.addClone("Pair/before/Pi0/uls/", "Pair/before/Pi0/lspp/");
90-
fRegistry.addClone("Pair/before/Pi0/uls/", "Pair/before/Pi0/lsmm/");
91-
fRegistry.addClone("Pair/before/Pi0/", "Pair/before/Eta/");
92-
fRegistry.addClone("Pair/before/Pi0/", "Pair/before/EtaPrime/");
93-
fRegistry.addClone("Pair/before/Pi0/", "Pair/before/Rho/");
94-
fRegistry.addClone("Pair/before/Pi0/", "Pair/before/Omega/");
95-
fRegistry.addClone("Pair/before/Pi0/", "Pair/before/Phi/");
96-
fRegistry.addClone("Pair/before/Pi0/", "Pair/before/JPsi/");
97-
fRegistry.addClone("Pair/before/Pi0/", "Pair/before/ccbar/");
98-
fRegistry.addClone("Pair/before/Pi0/", "Pair/before/bbbar/");
99-
fRegistry.addClone("Pair/before/", "Pair/after/");
85+
fRegistry.add("Event/hDiffBC", "diffrence in BC;BC_{rec. coll.} - BC_{mc coll.}", kTH1D, {{101, -50.5, +50.5}}, false);
86+
fRegistry.add("Event/allMC/hZvtx", "MC Zvtx;Z_{vtx} (cm)", kTH1D, {{100, -50, +50}}, false);
87+
fRegistry.add("Event/allMC/hImpactParameter", "impact parameter;impact parameter b (fm)", kTH1D, {{200, 0, 20}}, false);
88+
fRegistry.addClone("Event/allMC/", "Event/selectedMC/");
89+
fRegistry.addClone("Event/allMC/", "Event/selectedMC_and_Rec/");
90+
fRegistry.addClone("Event/allMC/", "Event/selectedMC_and_selectedRec/");
91+
92+
fRegistry.add("Pair/allMC/Pi0/uls/hMvsPt", "m_{ll} vs. p_{T,ll}", kTH2D, {axis_mll, axis_ptll}, true);
93+
fRegistry.addClone("Pair/allMC/Pi0/uls/", "Pair/allMC/Pi0/lspp/");
94+
fRegistry.addClone("Pair/allMC/Pi0/uls/", "Pair/allMC/Pi0/lsmm/");
95+
fRegistry.addClone("Pair/allMC/Pi0/", "Pair/allMC/Eta/");
96+
fRegistry.addClone("Pair/allMC/Pi0/", "Pair/allMC/EtaPrime/");
97+
fRegistry.addClone("Pair/allMC/Pi0/", "Pair/allMC/Rho/");
98+
fRegistry.addClone("Pair/allMC/Pi0/", "Pair/allMC/Omega/");
99+
fRegistry.addClone("Pair/allMC/Pi0/", "Pair/allMC/Phi/");
100+
fRegistry.addClone("Pair/allMC/Pi0/", "Pair/allMC/JPsi/");
101+
fRegistry.addClone("Pair/allMC/Pi0/", "Pair/allMC/ccbar/");
102+
fRegistry.addClone("Pair/allMC/Pi0/", "Pair/allMC/bbbar/");
103+
fRegistry.addClone("Pair/allMC/", "Pair/selectedMC/");
104+
fRegistry.addClone("Pair/allMC/", "Pair/selectedMC_and_Rec/");
105+
fRegistry.addClone("Pair/allMC/", "Pair/selectedMC_and_selectedRec/");
106+
107+
// allMC = all mc collisions
108+
// selectedMC = mc collisions selected by your event cuts
109+
// selectedMC_and_Rec = mc collisions selected by your event cuts and at least 1 reconstructed collision is associated.
110+
// selectedMC_and_selectedRec = mc collisions selected by your event cuts and at least 1 reconstructed collision is associated, and the associated rec. collision is selected by your cuts.
100111
}
101112

102113
template <typename TTrack, typename TMCParticles>
@@ -131,27 +142,22 @@ struct studyMCTruth {
131142
}
132143
}
133144

134-
template <typename TBCs, typename TMCCollision>
135-
bool isSelectedCollision(TMCCollision const& mcCollision)
145+
template <typename TCollision, typename TBC>
146+
bool isSelectedCollision(TCollision const& collision, TBC const& bc)
136147
{
137-
const auto& bc = mcCollision.template bc_as<TBCs>();
138-
139-
if (mcCollision.posZ() < eventcuts.cfgZvtxMin || eventcuts.cfgZvtxMax < mcCollision.posZ()) {
148+
if (collision.posZ() < eventcuts.cfgZvtxMin || eventcuts.cfgZvtxMax < collision.posZ()) {
140149
return false;
141150
}
142151

143152
if (eventcuts.cfgRequireFT0AND && !bc.selection_bit(o2::aod::evsel::kIsTriggerTVX)) {
144153
return false;
145154
}
146-
147155
if (eventcuts.cfgRequireNoTFB && !bc.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) {
148156
return false;
149157
}
150-
151158
if (eventcuts.cfgRequireNoITSROFB && !bc.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
152159
return false;
153160
}
154-
155161
return true;
156162
}
157163

@@ -225,41 +231,78 @@ struct studyMCTruth {
225231
}
226232
}
227233

228-
template <typename TMCCollisions, typename TMCParticles, typename TBCs>
229-
void runMC(TMCCollisions const& mcCollisions, TMCParticles const& mcParticles, TBCs const&)
234+
template <typename TMCCollisions, typename TMCParticles, typename TBCs, typename TCollisions>
235+
void runMC(TMCCollisions const& mcCollisions, TMCParticles const& mcParticles, TBCs const&, TCollisions const& collisions)
230236
{
231237
for (const auto& mcCollision : mcCollisions) {
232-
fRegistry.fill(HIST("Event/before/hZvtx"), mcCollision.posZ());
233-
fRegistry.fill(HIST("Event/before/hImpactParameter"), mcCollision.impactParameter());
238+
const auto& bc_from_mcCollision = mcCollision.template bc_as<TBCs>();
239+
bool isSelectedMC = isSelectedCollision(mcCollision, bc_from_mcCollision);
240+
241+
bool isSelectedRec = false;
242+
bool hasRecCollision = false;
243+
if (mcCollision.mpemeventId() >= 0) {
244+
hasRecCollision = true;
245+
const auto& collision = collisions.rawIteratorAt(mcCollision.mpemeventId()); // most probable reconstructed collision
246+
const auto& bc_from_collision = collision.template foundBC_as<TBCs>();
247+
isSelectedRec = isSelectedCollision(collision, bc_from_collision);
248+
fRegistry.fill(HIST("Event/hDiffBC"), bc_from_collision.globalBC() - bc_from_mcCollision.globalBC());
249+
}
250+
251+
fRegistry.fill(HIST("Event/allMC/hZvtx"), mcCollision.posZ());
252+
fRegistry.fill(HIST("Event/allMC/hImpactParameter"), mcCollision.impactParameter());
253+
if (isSelectedMC) {
254+
fRegistry.fill(HIST("Event/selectedMC/hZvtx"), mcCollision.posZ());
255+
fRegistry.fill(HIST("Event/selectedMC/hImpactParameter"), mcCollision.impactParameter());
256+
if (hasRecCollision) {
257+
fRegistry.fill(HIST("Event/selectedMC_and_Rec/hZvtx"), mcCollision.posZ());
258+
fRegistry.fill(HIST("Event/selectedMC_and_Rec/hImpactParameter"), mcCollision.impactParameter());
259+
if (isSelectedRec) {
260+
fRegistry.fill(HIST("Event/selectedMC_and_selectedRec/hZvtx"), mcCollision.posZ());
261+
fRegistry.fill(HIST("Event/selectedMC_and_selectedRec/hImpactParameter"), mcCollision.impactParameter());
262+
}
263+
}
264+
}
234265

235266
// store MC true information
236267
auto posLeptons_per_mccollision = mcPosLeptons.sliceBy(perMcCollision, mcCollision.globalIndex());
237268
auto negLeptons_per_mccollision = mcNegLeptons.sliceBy(perMcCollision, mcCollision.globalIndex());
238269

239-
bool isSelected = isSelectedCollision<TBCs>(mcCollision);
240-
if (isSelected) {
241-
fRegistry.fill(HIST("Event/after/hZvtx"), mcCollision.posZ());
242-
fRegistry.fill(HIST("Event/after/hImpactParameter"), mcCollision.impactParameter());
243-
}
244-
245270
for (const auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posLeptons_per_mccollision, negLeptons_per_mccollision))) { // ULS
246271
fillTrueInfo<0, 0>(pos, neg, mcParticles);
247-
if (isSelected) {
272+
if (isSelectedMC) {
248273
fillTrueInfo<1, 0>(pos, neg, mcParticles);
274+
if (hasRecCollision) {
275+
fillTrueInfo<2, 0>(pos, neg, mcParticles);
276+
if (isSelectedRec) {
277+
fillTrueInfo<3, 0>(pos, neg, mcParticles);
278+
}
279+
}
249280
}
250281
} // end of ULS pair loop
251282

252283
for (auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posLeptons_per_mccollision, posLeptons_per_mccollision))) { // LS++
253284
fillTrueInfo<0, 1>(pos1, pos2, mcParticles);
254-
if (isSelected) {
285+
if (isSelectedMC) {
255286
fillTrueInfo<1, 1>(pos1, pos2, mcParticles);
287+
if (hasRecCollision) {
288+
fillTrueInfo<2, 1>(pos1, pos2, mcParticles);
289+
if (isSelectedRec) {
290+
fillTrueInfo<3, 1>(pos1, pos2, mcParticles);
291+
}
292+
}
256293
}
257294
} // end of LS++ pair loop
258295

259296
for (auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negLeptons_per_mccollision, negLeptons_per_mccollision))) { // LS--
260297
fillTrueInfo<0, 2>(neg1, neg2, mcParticles);
261-
if (isSelected) {
298+
if (isSelectedMC) {
262299
fillTrueInfo<1, 2>(neg1, neg2, mcParticles);
300+
if (hasRecCollision) {
301+
fillTrueInfo<2, 2>(neg1, neg2, mcParticles);
302+
if (isSelectedRec) {
303+
fillTrueInfo<3, 2>(neg1, neg2, mcParticles);
304+
}
305+
}
263306
}
264307
} // end of LS-- pair loop
265308

@@ -270,19 +313,21 @@ struct studyMCTruth {
270313
SliceCache cache;
271314
Preslice<aod::McParticles> perMcCollision = aod::mcparticle::mcCollisionId;
272315

316+
using MyMcCollisions = soa::Join<aod::McCollisions, aod::MostProbableEMEventIdsInMC>;
317+
273318
Filter collisionFilter = eventcuts.cfgMinImpPar < o2::aod::mccollision::impactParameter && o2::aod::mccollision::impactParameter < eventcuts.cfgMaxImpPar;
274-
using FilteredMcCollisions = soa::Filtered<aod::McCollisions>;
319+
using FilteredMyMcCollisions = soa::Filtered<MyMcCollisions>;
275320

276321
Partition<aod::McParticles> mcPosLeptons = o2::aod::mcparticle::pdgCode == -mccuts.cfgPdgCodeLepton && (mccuts.cfgMinPtGen < o2::aod::mcparticle::pt && o2::aod::mcparticle::pt < mccuts.cfgMaxPtGen) && (mccuts.cfgMinEtaGen < o2::aod::mcparticle::eta && o2::aod::mcparticle::eta < mccuts.cfgMaxEtaGen);
277322
Partition<aod::McParticles> mcNegLeptons = o2::aod::mcparticle::pdgCode == mccuts.cfgPdgCodeLepton && (mccuts.cfgMinPtGen < o2::aod::mcparticle::pt && o2::aod::mcparticle::pt < mccuts.cfgMaxPtGen) && (mccuts.cfgMinEtaGen < o2::aod::mcparticle::eta && o2::aod::mcparticle::eta < mccuts.cfgMaxEtaGen);
278323

279-
void processMC(FilteredMcCollisions const& mcCollisions, aod::McParticles const& mcParticles, soa::Join<aod::BCsWithTimestamps, aod::BcSels> const& bcs)
324+
void processMC(FilteredMyMcCollisions const& mcCollisions, aod::McParticles const& mcParticles, soa::Join<aod::BCsWithTimestamps, aod::BcSels> const& bcs, soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels> const& collisions)
280325
{
281-
runMC(mcCollisions, mcParticles, bcs);
326+
runMC(mcCollisions, mcParticles, bcs, collisions);
282327
}
283328
PROCESS_SWITCH(studyMCTruth, processMC, "process", true);
284329

285-
void processDummy(FilteredMcCollisions const&) {}
330+
void processDummy(FilteredMyMcCollisions const&) {}
286331
PROCESS_SWITCH(studyMCTruth, processDummy, "process Dummy", false);
287332
};
288333

0 commit comments

Comments
 (0)