Skip to content

Commit 8f7a995

Browse files
dsekihatdsekihat
authored
PWGEM/PhotonMeson: fix MC indexing for daughters (#5600)
1 parent 575f012 commit 8f7a995

File tree

6 files changed

+110
-116
lines changed

6 files changed

+110
-116
lines changed

‎PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -403,12 +403,12 @@ void o2::aod::pwgem::photon::histogram::DefineHistograms(THashList* list, const
403403
}
404404

405405
const int nmgg = 401;
406-
float mgg[nmgg] = {};
406+
double mgg[nmgg] = {};
407407
for (int i = 0; i < nmgg; i++) {
408408
mgg[i] = 0.002 * i;
409409
}
410410
const int npTgg = 71;
411-
float pTgg[npTgg] = {};
411+
double pTgg[npTgg] = {};
412412
for (int i = 0; i < 50; i++) {
413413
pTgg[i] = 0.1 * (i - 0) + 0.0; // from 0 to 5 GeV/c, every 0.1 GeV/c
414414
}
@@ -527,12 +527,12 @@ void o2::aod::pwgem::photon::histogram::DefineHistograms(THashList* list, const
527527
}
528528

529529
const int nmgg04 = 201;
530-
float mgg04[nmgg04] = {};
530+
double mgg04[nmgg04] = {};
531531
for (int i = 0; i < nmgg04; i++) {
532532
mgg04[i] = 0.002 * i;
533533
}
534534
const int npTgg10 = 61;
535-
float pTgg10[npTgg10] = {};
535+
double pTgg10[npTgg10] = {};
536536
for (int i = 0; i < 50; i++) {
537537
pTgg10[i] = 0.1 * (i - 0) + 0.0; // from 0 to 5 GeV/c, every 0.1 GeV/c
538538
}

‎PWGEM/PhotonMeson/DataModel/gammaTables.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -141,7 +141,7 @@ namespace emmcparticle
141141
{
142142
DECLARE_SOA_INDEX_COLUMN(EMMCEvent, emmcevent);
143143
DECLARE_SOA_SELF_ARRAY_INDEX_COLUMN(Mothers, mothers); //! Mother tracks (possible empty) array. Iterate over mcParticle.mothers_as<aod::McParticles>())
144-
DECLARE_SOA_SELF_SLICE_INDEX_COLUMN(Daughters, daughters); //! Daughter tracks (possibly empty) slice. Check for non-zero with mcParticle.has_daughters(). Iterate over mcParticle.daughters_as<aod::McParticles>())
144+
DECLARE_SOA_SELF_ARRAY_INDEX_COLUMN(Daughters, daughters); //! Daughter tracks (possibly empty) array. Check for non-zero with mcParticle.has_daughters(). Iterate over mcParticle.daughters_as<aod::McParticles>())
145145
DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float px, float py) -> float { return RecoDecay::sqrtSumOfSquares(px, py); });
146146
DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, [](float px, float py, float pz) -> float { return RecoDecay::eta(std::array{px, py, pz}); });
147147
DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, [](float px, float py) -> float { return RecoDecay::phi(px, py); });
@@ -160,7 +160,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(Y, y, //! Particle rapidity
160160
DECLARE_SOA_TABLE_FULL(EMMCParticles, "EMMCParticles", "AOD", "EMMCPARTICLE", //! MC track information (on disk)
161161
o2::soa::Index<>, emmcparticle::EMMCEventId,
162162
mcparticle::PdgCode, mcparticle::Flags,
163-
emmcparticle::MothersIds, emmcparticle::DaughtersIdSlice,
163+
emmcparticle::MothersIds, emmcparticle::DaughtersIds,
164164
mcparticle::Px, mcparticle::Py, mcparticle::Pz, mcparticle::E,
165165
mcparticle::Vx, mcparticle::Vy, mcparticle::Vz, mcparticle::Vt,
166166

‎PWGEM/PhotonMeson/TableProducer/associateMCinfo.cxx‎

Lines changed: 5 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -433,13 +433,15 @@ struct AssociateMCInfo {
433433
}
434434
}
435435

436-
// TODO: Check that the daughter slice in the skimmed table works as expected
437-
// Note that not all daughters from the original table are preserved in the skimmed MC stack
436+
// Note that not all daughters from the original table are preserved in the skimmed MC stack
438437
std::vector<int> daughters;
439438
if (mctrack.has_daughters()) {
439+
// LOGF(info, "daughter range in original MC stack pdg = %d | %d - %d , n dau = %d", mctrack.pdgCode(), mctrack.daughtersIds()[0], mctrack.daughtersIds()[1], mctrack.daughtersIds()[1] -mctrack.daughtersIds()[0] +1);
440440
for (int d = mctrack.daughtersIds()[0]; d <= mctrack.daughtersIds()[1]; ++d) {
441441
// TODO: remove this check as soon as issues with MC production are fixed
442442
if (d < mcTracks.size()) { // protect against bad daughter indices
443+
// auto dau_tmp = mcTracks.iteratorAt(d);
444+
// LOGF(info, "daughter pdg = %d", dau_tmp.pdgCode());
443445
if (fNewLabels.find(d) != fNewLabels.end()) {
444446
daughters.push_back(fNewLabels.find(d)->second);
445447
}
@@ -449,14 +451,9 @@ struct AssociateMCInfo {
449451
}
450452
}
451453
}
452-
int daughterRange[2] = {-1, -1};
453-
if (daughters.size() > 0) {
454-
daughterRange[0] = daughters[0];
455-
daughterRange[1] = daughters[daughters.size() - 1];
456-
}
457454

458455
emmcparticles(fEventIdx.find(oldLabel)->second, mctrack.pdgCode(), mctrack.flags(),
459-
mothers, daughterRange,
456+
mothers, daughters,
460457
mctrack.px(), mctrack.py(), mctrack.pz(), mctrack.e(),
461458
mctrack.vx(), mctrack.vy(), mctrack.vz(), mctrack.vt());
462459
} // end loop over labels

‎PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx‎

Lines changed: 37 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -466,22 +466,12 @@ struct Pi0EtaToGammaGammaMC {
466466
auto g1mc = mcparticles.iteratorAt(photonid1);
467467
auto g2mc = mcparticles.iteratorAt(photonid2);
468468

469-
if constexpr (pairtype == PairType::kPCMPCM) {
470-
if (!IsConversionPointInAcceptance(g1mc, maxRgen, maxY_track, margin_z_mc, mcparticles)) {
471-
continue;
472-
}
473-
if (!IsConversionPointInAcceptance(g2mc, maxRgen, maxY_track, margin_z_mc, mcparticles)) {
474-
continue;
475-
}
476-
}
477-
478469
pi0id = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 111, mcparticles);
479470
etaid = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 221, mcparticles);
480471

481472
if (pi0id < 0 && etaid < 0) {
482473
continue;
483474
}
484-
485475
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
486476
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
487477
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
@@ -492,12 +482,25 @@ struct Pi0EtaToGammaGammaMC {
492482
if (pi0id > 0) {
493483
auto pi0mc = mcparticles.iteratorAt(pi0id);
494484
if (pi0mc.isPhysicalPrimary() || pi0mc.producedByGenerator()) {
485+
if constexpr (pairtype == PairType::kPCMPCM) {
486+
if (!IsConversionPointInAcceptance(g1mc, maxRgen, maxY_track, margin_z_mc, mcparticles) || !IsConversionPointInAcceptance(g2mc, maxRgen, maxY_track, margin_z_mc, mcparticles)) {
487+
continue;
488+
}
489+
}
495490
reinterpret_cast<TH2F*>(list_pair_ss->FindObject(Form("%s_%s", cut.GetName(), cut.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Pi0_Primary"))->Fill(v12.M(), v12.Pt());
496491
} else if (IsFromWD(pi0mc.emmcevent(), pi0mc, mcparticles)) {
497492
reinterpret_cast<TH2F*>(list_pair_ss->FindObject(Form("%s_%s", cut.GetName(), cut.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Pi0_FromWD"))->Fill(v12.M(), v12.Pt());
498493
}
499494
} else if (etaid > 0) {
500-
reinterpret_cast<TH2F*>(list_pair_ss->FindObject(Form("%s_%s", cut.GetName(), cut.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Eta_Primary"))->Fill(v12.M(), v12.Pt());
495+
auto etamc = mcparticles.iteratorAt(etaid);
496+
if (etamc.isPhysicalPrimary() || etamc.producedByGenerator()) {
497+
if constexpr (pairtype == PairType::kPCMPCM) {
498+
if (!IsConversionPointInAcceptance(g1mc, maxRgen, maxY_track, margin_z_mc, mcparticles) || !IsConversionPointInAcceptance(g2mc, maxRgen, maxY_track, margin_z_mc, mcparticles)) {
499+
continue;
500+
}
501+
}
502+
reinterpret_cast<TH2F*>(list_pair_ss->FindObject(Form("%s_%s", cut.GetName(), cut.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Eta_Primary"))->Fill(v12.M(), v12.Pt());
503+
}
501504
}
502505
} // end of combination
503506
} // end of paircutloop
@@ -515,6 +518,9 @@ struct Pi0EtaToGammaGammaMC {
515518
continue;
516519
}
517520

521+
int photonid1 = -1;
522+
int photonid2 = -1;
523+
518524
if constexpr (pairtype == PairType::kPCMPHOS || pairtype == PairType::kPCMEMC) {
519525
auto pos = g1.template posTrack_as<MyMCV0Legs>();
520526
auto ele = g1.template negTrack_as<MyMCV0Legs>();
@@ -542,20 +548,14 @@ struct Pi0EtaToGammaGammaMC {
542548
auto ele2mc = ele2.template emmcparticle_as<aod::EMMCParticles>();
543549
// LOGF(info,"pos1mc.globalIndex() = %d , ele1mc.globalIndex() = %d , pos2mc.globalIndex() = %d , ele2mc.globalIndex() = %d", pos1mc.globalIndex(), ele1mc.globalIndex(), pos2mc.globalIndex(), ele2mc.globalIndex());
544550

545-
int photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); // real photon
551+
photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); // real photon
546552
if (photonid1 < 0) {
547553
continue;
548554
}
549555
auto g1mc = mcparticles.iteratorAt(photonid1);
550556

551-
if constexpr (pairtype == PairType::kPCMDalitzEE || pairtype == kPCMDalitzMuMu) {
552-
if (!IsConversionPointInAcceptance(g1mc, maxRgen, maxY_track, margin_z_mc, mcparticles)) {
553-
continue;
554-
}
555-
}
556-
557-
if (cut2.IsPhotonConversionSelected()) { // v0photon + photon conversion on ITSib stored in dielectron table. pi0 -> gamma gamma
558-
int photonid2 = FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles); // photon conversion stored in dielectron table
557+
if (cut2.IsPhotonConversionSelected()) { // v0photon + photon conversion on ITSib stored in dielectron table. pi0 -> gamma gamma
558+
photonid2 = FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles); // photon conversion stored in dielectron table
559559
if (photonid2 < 0) {
560560
continue;
561561
}
@@ -581,7 +581,7 @@ struct Pi0EtaToGammaGammaMC {
581581
auto ele2mc = ele2.template emmcparticle_as<aod::EMMCParticles>();
582582
// LOGF(info,"pos1mc.globalIndex() = %d , ele1mc.globalIndex() = %d , pos2mc.globalIndex() = %d , ele2mc.globalIndex() = %d", pos1mc.globalIndex(), ele1mc.globalIndex(), pos2mc.globalIndex(), ele2mc.globalIndex());
583583

584-
int photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); // real photon
584+
photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); // real photon
585585
if (photonid1 < 0) {
586586
continue;
587587
}
@@ -606,12 +606,27 @@ struct Pi0EtaToGammaGammaMC {
606606
if (pi0id > 0) {
607607
auto pi0mc = mcparticles.iteratorAt(pi0id);
608608
if (pi0mc.isPhysicalPrimary() || pi0mc.producedByGenerator()) {
609+
if constexpr (pairtype == PairType::kPCMDalitzEE || pairtype == kPCMDalitzMuMu) {
610+
auto g1mc = mcparticles.iteratorAt(photonid1);
611+
if (!IsConversionPointInAcceptance(g1mc, maxRgen, maxY_track, margin_z_mc, mcparticles)) {
612+
continue;
613+
}
614+
}
609615
reinterpret_cast<TH2F*>(list_pair_ss->FindObject(Form("%s_%s", cut1.GetName(), cut2.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Pi0_Primary"))->Fill(v12.M(), v12.Pt());
610616
} else if (IsFromWD(pi0mc.emmcevent(), pi0mc, mcparticles)) {
611617
reinterpret_cast<TH2F*>(list_pair_ss->FindObject(Form("%s_%s", cut1.GetName(), cut2.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Pi0_FromWD"))->Fill(v12.M(), v12.Pt());
612618
}
613619
} else if (etaid > 0) {
614-
reinterpret_cast<TH2F*>(list_pair_ss->FindObject(Form("%s_%s", cut1.GetName(), cut2.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Eta_Primary"))->Fill(v12.M(), v12.Pt());
620+
auto etamc = mcparticles.iteratorAt(etaid);
621+
if (etamc.isPhysicalPrimary() || etamc.producedByGenerator()) {
622+
if constexpr (pairtype == PairType::kPCMDalitzEE || pairtype == kPCMDalitzMuMu) {
623+
auto g1mc = mcparticles.iteratorAt(photonid1);
624+
if (!IsConversionPointInAcceptance(g1mc, maxRgen, maxY_track, margin_z_mc, mcparticles)) {
625+
continue;
626+
}
627+
}
628+
reinterpret_cast<TH2F*>(list_pair_ss->FindObject(Form("%s_%s", cut1.GetName(), cut2.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Eta_Primary"))->Fill(v12.M(), v12.Pt());
629+
}
615630
}
616631

617632
} // end of combination
@@ -656,25 +671,6 @@ struct Pi0EtaToGammaGammaMC {
656671
}
657672

658673
auto mccollision = collision.emmcevent();
659-
reinterpret_cast<TH1F*>(list_gen_pair->FindObject("hZvtx_before"))->Fill(mccollision.posZ());
660-
reinterpret_cast<TH1F*>(list_gen_pair->FindObject("hCollisionCounter"))->Fill(1.0); // all
661-
if (!collision.sel8()) {
662-
continue;
663-
}
664-
reinterpret_cast<TH1F*>(list_gen_pair->FindObject("hCollisionCounter"))->Fill(2.0); // FT0VX i.e. FT0and
665-
666-
if (collision.numContrib() < 0.5) {
667-
continue;
668-
}
669-
reinterpret_cast<TH1F*>(list_gen_pair->FindObject("hCollisionCounter"))->Fill(3.0); // Ncontrib > 0
670-
671-
if (abs(collision.posZ()) > 10.0) {
672-
continue;
673-
}
674-
675-
reinterpret_cast<TH1F*>(list_gen_pair->FindObject("hZvtx_after"))->Fill(mccollision.posZ());
676-
reinterpret_cast<TH1F*>(list_gen_pair->FindObject("hCollisionCounter"))->Fill(4.0); // |Zvtx| < 10 cm
677-
678674
auto mctracks_coll = mcparticles.sliceBy(perMcCollision, mccollision.globalIndex());
679675
for (auto& mctrack : mctracks_coll) {
680676
if (abs(mctrack.y()) > maxY_track) {

0 commit comments

Comments
 (0)