Skip to content

Commit 9e00cf5

Browse files
dsekihatdsekihat
authored
PWGEM/PhotonMeson: update dilepton QC MC (#3906)
1 parent 496e128 commit 9e00cf5

File tree

5 files changed

+170
-20
lines changed

5 files changed

+170
-20
lines changed

‎PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx‎

Lines changed: 9 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -202,11 +202,6 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
202202
list->Add(new TH1F("hNpair_lspp", "Number of LS++ pairs per collision", 101, -0.5f, 100.5f));
203203
list->Add(new TH1F("hNpair_lsmm", "Number of LS-- pairs per collision", 101, -0.5f, 100.5f));
204204

205-
if (TString(histClass) == "Generated") {
206-
TH2F* hMvsPt = new TH2F("hMvsPt", "m_{ll} vs. p_{T,ll};m_{ll} (GeV/c^{2});p_{T,ll}", 110, 0, 1.1f, 100, 0, 10.f);
207-
hMvsPt->Sumw2();
208-
list->Add(hMvsPt);
209-
}
210205
} // end of Dalitz
211206
if (TString(histClass) == "Track") {
212207
float maxP = 10.f;
@@ -387,14 +382,12 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
387382
list->Add(new TH2F("hPhotonPhivsRxy", "conversion point of #varphi vs. R_{xy} MC;#varphi (rad.);R_{xy} (cm);N_{e}", 360, 0.0f, TMath::TwoPi(), 200, 0, 200));
388383
}
389384

390-
// Generated, particles
391385
if (TString(subGroup) == "Photon") {
392386
list->Add(new TH1F("hPt_Photon", "photon pT;p_{T} (GeV/c)", 1000, 0.0f, 10));
393387
list->Add(new TH1F("hY_Photon", "photon y;rapidity y", 40, -2.0f, 2.0f));
394388
list->Add(new TH1F("hPhi_Photon", "photon #varphi;#varphi (rad.)", 180, 0, TMath::TwoPi()));
395389
}
396390

397-
// Generated, particles
398391
if (TString(subGroup) == "Pi0Eta") {
399392
static constexpr std::string_view parnames[2] = {"Pi0", "Eta"};
400393
for (int i = 0; i < 2; i++) {
@@ -403,6 +396,15 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
403396
list->Add(new TH1F(Form("hPhi_%s", parnames[i].data()), Form("%s #varphi;#varphi (rad.)", parnames[i].data()), 180, 0, TMath::TwoPi()));
404397
}
405398
}
399+
if (TString(subGroup) == "dielectron") {
400+
TH2F* hMvsPt = new TH2F("hMvsPt", "m_{ee} vs. p_{T,ee};m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c)", 110, 0, 1.1f, 100, 0, 10.f);
401+
hMvsPt->Sumw2();
402+
list->Add(hMvsPt);
403+
} else if (TString(subGroup) == "dimuon") {
404+
TH2F* hMvsPt = new TH2F("hMvsPt", "m_{#mu#mu} vs. p_{T,#mu#mu};m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c)", 90, 0.2, 1.1f, 20, 0, 2.f);
405+
hMvsPt->Sumw2();
406+
list->Add(hMvsPt);
407+
}
406408
}
407409

408410
const int nmgg04 = 201;

‎PWGEM/PhotonMeson/TableProducer/associateMCinfo.cxx‎

Lines changed: 2 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -99,7 +99,7 @@ struct AssociateMCInfo {
9999
auto groupedMcTracks = mcTracks.sliceBy(perMcCollision, mcCollision.globalIndex());
100100

101101
for (auto& mctrack : groupedMcTracks) {
102-
if (mctrack.pt() < 1e-2 || abs(mctrack.y()) > 1.5 || abs(mctrack.vz()) > 250 || sqrt(pow(mctrack.vx(), 2) + pow(mctrack.vy(), 2)) > max_rxy_gen) {
102+
if (mctrack.pt() < 1e-3 || abs(mctrack.y()) > 1.5 || abs(mctrack.vz()) > 250 || sqrt(pow(mctrack.vx(), 2) + pow(mctrack.vy(), 2)) > max_rxy_gen) {
103103
continue;
104104
}
105105
int pdg = mctrack.pdgCode();
@@ -143,9 +143,6 @@ struct AssociateMCInfo {
143143
auto ele = v0.template negTrack_as<aod::V0Legs>();
144144
auto pos = v0.template posTrack_as<aod::V0Legs>();
145145

146-
// auto o2track_ele = ele.template track_as<TracksMC>();
147-
// auto o2track_pos = pos.template track_as<TracksMC>();
148-
149146
auto o2track_ele = o2tracks.iteratorAt(pos.trackId());
150147
auto o2track_pos = o2tracks.iteratorAt(ele.trackId());
151148

@@ -192,7 +189,7 @@ struct AssociateMCInfo {
192189
} // end of em primary track loop
193190
}
194191
if constexpr (static_cast<bool>(system & kDalitzMuMu)) {
195-
// for dalitz ee
192+
// for dalitz mumu
196193
auto emprimarymuons_coll = emprimarymuons.sliceBy(perCollision_mu, collision.globalIndex());
197194
for (auto& emprimarymuon : emprimarymuons_coll) {
198195
auto o2track = o2tracks.iteratorAt(emprimarymuon.trackId());

‎PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx‎

Lines changed: 25 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -347,13 +347,6 @@ struct Pi0EtaToGammaGammaMC {
347347
int photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles);
348348
int photonid2 = FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles);
349349

350-
// if (photonid1 < 0) { // check swap, true electron is reconstructed as positron and vice versa.
351-
// photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, 11, -11, 22, mcparticles);
352-
// }
353-
// if (photonid2 < 0) { // check swap, true electron is reconstructed as positron and vice versa.
354-
// photonid2 = FindCommonMotherFrom2Prongs(pos2mc, ele2mc, 11, -11, 22, mcparticles);
355-
// }
356-
357350
// LOGF(info,"photonid1 = %d , photonid2 = %d", photonid1, photonid2);
358351
if (photonid1 < 0 || photonid2 < 0) {
359352
continue;
@@ -435,14 +428,38 @@ struct Pi0EtaToGammaGammaMC {
435428
auto g1mc = mcparticles.iteratorAt(photonid1);
436429
pi0id = FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -11, 11, 111, mcparticles);
437430
etaid = FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -11, 11, 221, mcparticles);
431+
432+
} else if constexpr (pairtype == PairType::kPCMDalitzMuMu) { // check 4 legs
433+
auto pos1 = g1.template posTrack_as<MyMCV0Legs>();
434+
auto ele1 = g1.template negTrack_as<MyMCV0Legs>();
435+
auto pos2 = g2.template posTrack_as<MyMCMuons>();
436+
auto ele2 = g2.template negTrack_as<MyMCMuons>();
437+
if (pos1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) {
438+
continue;
439+
}
440+
441+
auto pos1mc = pos1.template emmcparticle_as<aod::EMMCParticles>();
442+
auto ele1mc = ele1.template emmcparticle_as<aod::EMMCParticles>();
443+
auto pos2mc = pos2.template emmcparticle_as<aod::EMMCParticles>();
444+
auto ele2mc = ele2.template emmcparticle_as<aod::EMMCParticles>();
445+
// LOGF(info,"pos1mc.globalIndex() = %d , ele1mc.globalIndex() = %d , pos2mc.globalIndex() = %d , ele2mc.globalIndex() = %d", pos1mc.globalIndex(), ele1mc.globalIndex(), pos2mc.globalIndex(), ele2mc.globalIndex());
446+
447+
int photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); // real photon
448+
if (photonid1 < 0) {
449+
continue;
450+
}
451+
auto g1mc = mcparticles.iteratorAt(photonid1);
452+
pi0id = -1;
453+
etaid = FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -13, 13, 221, mcparticles);
438454
}
455+
439456
if (pi0id < 0 && etaid < 0) {
440457
continue;
441458
}
442459

443460
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
444461
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
445-
if constexpr (pairtype == PairType::kPCMDalitzEE) {
462+
if constexpr (pairtype == PairType::kPCMDalitzEE || pairtype == PairType::kPCMDalitzMuMu) {
446463
v2.SetM(g2.mass());
447464
}
448465
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;

‎PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx‎

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@ struct DalitzEEQCMC {
5353
OutputObj<THashList> fOutputEvent{"Event"};
5454
OutputObj<THashList> fOutputTrack{"Track"};
5555
OutputObj<THashList> fOutputDalitzEE{"DalitzEE"};
56+
OutputObj<THashList> fOutputGen{"Generated"};
5657
THashList* fMainList = new THashList();
5758

5859
void addhistograms()
@@ -71,6 +72,10 @@ struct DalitzEEQCMC {
7172
o2::aod::emphotonhistograms::AddHistClass(fMainList, "DalitzEE");
7273
THashList* list_dalitzee = reinterpret_cast<THashList*>(fMainList->FindObject("DalitzEE"));
7374

75+
o2::aod::emphotonhistograms::AddHistClass(fMainList, "Generated");
76+
THashList* list_gen = reinterpret_cast<THashList*>(fMainList->FindObject("Generated"));
77+
o2::aod::emphotonhistograms::DefineHistograms(list_gen, "Generated", "dielectron");
78+
7479
for (const auto& cut : fDalitzEECuts) {
7580
const char* cutname = cut.GetName();
7681
o2::aod::emphotonhistograms::AddHistClass(list_track, cutname);
@@ -114,6 +119,7 @@ struct DalitzEEQCMC {
114119
fOutputEvent.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Event")));
115120
fOutputTrack.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Track")));
116121
fOutputDalitzEE.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("DalitzEE")));
122+
fOutputGen.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Generated")));
117123
}
118124

119125
template <typename TTrack, typename TMCTracks>
@@ -220,6 +226,66 @@ struct DalitzEEQCMC {
220226
} // end of process
221227
PROCESS_SWITCH(DalitzEEQCMC, processQCMC, "run Dalitz QC", true);
222228

229+
Configurable<float> min_mcPt{"min_mcPt", 0.05, "min. MC pT"};
230+
Configurable<float> max_mcPt{"max_mcPt", 1e+10, "max. MC pT"};
231+
Configurable<float> max_mcEta{"max_mcEta", 0.9, "max. MC eta"};
232+
Partition<aod::EMMCParticles> posTracks = o2::aod::mcparticle::pdgCode == -11; // e+
233+
Partition<aod::EMMCParticles> negTracks = o2::aod::mcparticle::pdgCode == +11; // e-
234+
PresliceUnsorted<aod::EMMCParticles> perMcCollision = aod::emmcparticle::emreducedmceventId;
235+
void processGen(MyCollisions const& collisions, aod::EMReducedMCEvents const&, aod::EMMCParticles const& mcparticles)
236+
{
237+
// loop over mc stack and fill histograms for pure MC truth signals
238+
// all MC tracks which belong to the MC event corresponding to the current reconstructed event
239+
for (auto& collision : collisions) {
240+
auto mccollision = collision.emreducedmcevent();
241+
// LOGF(info, "mccollision.globalIndex() = %d", mccollision.globalIndex());
242+
// auto mctracks_coll = mcparticles.sliceBy(perMcCollision, mccollision.globalIndex());
243+
244+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(1.0);
245+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hZvtx_before"))->Fill(mccollision.posZ());
246+
if (!collision.sel8()) {
247+
continue;
248+
}
249+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(2.0);
250+
251+
if (collision.numContrib() < 0.5) {
252+
continue;
253+
}
254+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(3.0);
255+
256+
if (abs(collision.posZ()) > 10.0) {
257+
continue;
258+
}
259+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(4.0);
260+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hZvtx_after"))->Fill(mccollision.posZ());
261+
auto posTracks_per_coll = posTracks->sliceByCachedUnsorted(o2::aod::emmcparticle::emreducedmceventId, mccollision.globalIndex(), cache);
262+
auto negTracks_per_coll = negTracks->sliceByCachedUnsorted(o2::aod::emmcparticle::emreducedmceventId, mccollision.globalIndex(), cache);
263+
for (auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) {
264+
// LOGF(info, "pdg1 = %d, pdg2 = %d", t1.pdgCode(), t2.pdgCode());
265+
266+
if (t1.pt() < min_mcPt || max_mcPt < t1.pt() || abs(t1.eta()) > max_mcEta) {
267+
continue;
268+
}
269+
if (t2.pt() < min_mcPt || max_mcPt < t2.pt() || abs(t2.eta()) > max_mcEta) {
270+
continue;
271+
}
272+
273+
int mother_id = FindLF(t1, t2, mcparticles);
274+
if (mother_id < 0) {
275+
continue;
276+
}
277+
auto mcmother = mcparticles.iteratorAt(mother_id);
278+
if (IsPhysicalPrimary(mcmother.emreducedmcevent(), mcmother, mcparticles)) {
279+
ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron);
280+
ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron);
281+
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
282+
reinterpret_cast<TH2F*>(fMainList->FindObject("Generated")->FindObject("hMvsPt"))->Fill(v12.M(), v12.Pt());
283+
} // end of LF
284+
} // end of true ULS pair loop
285+
} // end of collision loop
286+
}
287+
PROCESS_SWITCH(DalitzEEQCMC, processGen, "run genrated info", true);
288+
223289
void processDummy(MyCollisions const& collisions) {}
224290
PROCESS_SWITCH(DalitzEEQCMC, processDummy, "Dummy function", false);
225291
};

‎PWGEM/PhotonMeson/Tasks/dalitzMuMuQCMC.cxx‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@ struct DalitzMuMuQCMC {
5353
OutputObj<THashList> fOutputEvent{"Event"};
5454
OutputObj<THashList> fOutputTrack{"Track"};
5555
OutputObj<THashList> fOutputDalitzMuMu{"DalitzMuMu"};
56+
OutputObj<THashList> fOutputGen{"Generated"};
5657
THashList* fMainList = new THashList();
5758

5859
void addhistograms()
@@ -71,6 +72,10 @@ struct DalitzMuMuQCMC {
7172
o2::aod::emphotonhistograms::AddHistClass(fMainList, "DalitzMuMu");
7273
THashList* list_dalitzmumu = reinterpret_cast<THashList*>(fMainList->FindObject("DalitzMuMu"));
7374

75+
o2::aod::emphotonhistograms::AddHistClass(fMainList, "Generated");
76+
THashList* list_gen = reinterpret_cast<THashList*>(fMainList->FindObject("Generated"));
77+
o2::aod::emphotonhistograms::DefineHistograms(list_gen, "Generated", "dimuon");
78+
7479
for (const auto& cut : fDalitzMuMuCuts) {
7580
const char* cutname = cut.GetName();
7681
o2::aod::emphotonhistograms::AddHistClass(list_track, cutname);
@@ -114,6 +119,7 @@ struct DalitzMuMuQCMC {
114119
fOutputEvent.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Event")));
115120
fOutputTrack.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Track")));
116121
fOutputDalitzMuMu.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("DalitzMuMu")));
122+
fOutputGen.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Generated")));
117123
}
118124

119125
template <typename TTrack, typename TMCTracks>
@@ -177,6 +183,8 @@ struct DalitzMuMuQCMC {
177183
auto posmc = pos.template emmcparticle_as<aod::EMMCParticles>();
178184
auto elemc = ele.template emmcparticle_as<aod::EMMCParticles>();
179185

186+
// LOGF(info, "posmc.pdgCode() = %d, , posmc.has_mothers()() = %d | elemc.pdgCode() = %d, elemc.has_mothers() = %d", posmc.pdgCode(), posmc.has_mothers(), elemc.pdgCode(), elemc.has_mothers());
187+
180188
int mother_id = FindLF(posmc, elemc, mcparticles);
181189
if (mother_id < 0) {
182190
continue;
@@ -209,6 +217,66 @@ struct DalitzMuMuQCMC {
209217
} // end of process
210218
PROCESS_SWITCH(DalitzMuMuQCMC, processQCMC, "run Dalitz QC", true);
211219

220+
Configurable<float> min_mcPt{"min_mcPt", 0.05, "min. MC pT"};
221+
Configurable<float> max_mcPt{"max_mcPt", 0.6, "max. MC pT"};
222+
Configurable<float> max_mcEta{"max_mcEta", 0.9, "max. MC eta"};
223+
Partition<aod::EMMCParticles> posTracks = o2::aod::mcparticle::pdgCode == -13; // mu+
224+
Partition<aod::EMMCParticles> negTracks = o2::aod::mcparticle::pdgCode == +13; // mu-
225+
PresliceUnsorted<aod::EMMCParticles> perMcCollision = aod::emmcparticle::emreducedmceventId;
226+
void processGen(MyCollisions const& collisions, aod::EMReducedMCEvents const&, aod::EMMCParticles const& mcparticles)
227+
{
228+
// loop over mc stack and fill histograms for pure MC truth signals
229+
// all MC tracks which belong to the MC event corresponding to the current reconstructed event
230+
for (auto& collision : collisions) {
231+
auto mccollision = collision.emreducedmcevent();
232+
// LOGF(info, "mccollision.globalIndex() = %d", mccollision.globalIndex());
233+
// auto mctracks_coll = mcparticles.sliceBy(perMcCollision, mccollision.globalIndex());
234+
235+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(1.0);
236+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hZvtx_before"))->Fill(mccollision.posZ());
237+
if (!collision.sel8()) {
238+
continue;
239+
}
240+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(2.0);
241+
242+
if (collision.numContrib() < 0.5) {
243+
continue;
244+
}
245+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(3.0);
246+
247+
if (abs(collision.posZ()) > 10.0) {
248+
continue;
249+
}
250+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(4.0);
251+
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hZvtx_after"))->Fill(mccollision.posZ());
252+
auto posTracks_per_coll = posTracks->sliceByCachedUnsorted(o2::aod::emmcparticle::emreducedmceventId, mccollision.globalIndex(), cache);
253+
auto negTracks_per_coll = negTracks->sliceByCachedUnsorted(o2::aod::emmcparticle::emreducedmceventId, mccollision.globalIndex(), cache);
254+
for (auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) {
255+
// LOGF(info, "pdg1 = %d, pdg2 = %d", t1.pdgCode(), t2.pdgCode());
256+
257+
if (t1.pt() < min_mcPt || max_mcPt < t1.pt() || abs(t1.eta()) > max_mcEta) {
258+
continue;
259+
}
260+
if (t2.pt() < min_mcPt || max_mcPt < t2.pt() || abs(t2.eta()) > max_mcEta) {
261+
continue;
262+
}
263+
264+
int mother_id = FindLF(t1, t2, mcparticles);
265+
if (mother_id < 0) {
266+
continue;
267+
}
268+
auto mcmother = mcparticles.iteratorAt(mother_id);
269+
if (IsPhysicalPrimary(mcmother.emreducedmcevent(), mcmother, mcparticles)) {
270+
ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassMuon);
271+
ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassMuon);
272+
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
273+
reinterpret_cast<TH2F*>(fMainList->FindObject("Generated")->FindObject("hMvsPt"))->Fill(v12.M(), v12.Pt());
274+
} // end of LF
275+
} // end of true ULS pair loop
276+
} // end of collision loop
277+
}
278+
PROCESS_SWITCH(DalitzMuMuQCMC, processGen, "run genrated info", true);
279+
212280
void processDummy(MyCollisions const& collisions) {}
213281
PROCESS_SWITCH(DalitzMuMuQCMC, processDummy, "Dummy function", false);
214282
};

0 commit comments

Comments
 (0)