Skip to content

Commit e8827fe

Browse files
rolavickrolavick
committed
MC tau tables ready
1 parent 0bb6fb9 commit e8827fe

File tree

2 files changed

+103
-62
lines changed

2 files changed

+103
-62
lines changed

‎PWGUD/DataModel/TauEventTables.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,8 @@ DECLARE_SOA_COLUMN(TrueDaugPx, trueDaugPx, float[2]);
9494
DECLARE_SOA_COLUMN(TrueDaugPy, trueDaugPy, float[2]);
9595
DECLARE_SOA_COLUMN(TrueDaugPz, trueDaugPz, float[2]);
9696
DECLARE_SOA_COLUMN(TrueDaugPdgCode, trueDaugPdgCode, int[2]);
97+
// additional info
98+
DECLARE_SOA_COLUMN(ProblematicEvent, problematicEvent, bool);
9799

98100
} // namespace tau_tree
99101
DECLARE_SOA_TABLE(TauTwoTracks, "AOD", "TAUTWOTRACK",
@@ -215,7 +217,8 @@ DECLARE_SOA_TABLE(TrueTauTwoTracks, "AOD", "TRUETAUTWOTRACK",
215217
tau_tree::TrueDaugPx,
216218
tau_tree::TrueDaugPy,
217219
tau_tree::TrueDaugPz,
218-
tau_tree::TrueDaugPdgCode);
220+
tau_tree::TrueDaugPdgCode,
221+
tau_tree::ProblematicEvent);
219222

220223
} // namespace o2::aod
221224

‎PWGUD/TableProducer/tauEventTableProducer.cxx‎

Lines changed: 99 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,6 @@
4040
#include "Common/DataModel/TrackSelectionTables.h"
4141
#include "PWGUD/Core/UPCTauCentralBarrelHelperRL.h"
4242
#include "PWGUD/DataModel/UDTables.h"
43-
#include "PWGUD/DataModel/UDIndex.h" // for UDMcParticles2UDTracks table
4443
#include "PWGUD/DataModel/TauEventTables.h"
4544
#include "PWGUD/Core/SGSelector.h"
4645

@@ -125,8 +124,6 @@ struct TauEventTableProducer {
125124
using FullMCUDTracks = soa::Join<aod::UDTracks, aod::UDTracksExtra, aod::UDTracksDCA, aod::UDTracksPID, aod::UDTracksFlags, aod::UDMcTrackLabels>;
126125
using FullMCSGUDCollisions = soa::Join<aod::UDCollisions, aod::UDCollisionsSels, aod::UDCollisionSelExtras, aod::SGCollisions, aod::UDMcCollsLabels>;
127126
using FullMCSGUDCollision = FullMCSGUDCollisions::iterator;
128-
using UDMcParticlesWithUDTracks = soa::Join<aod::UDMcParticles, aod::UDMcParticlesToUDTracks>;
129-
using UDMcCollisionsWithUDCollisions = soa::Join<aod::UDMcCollisions, aod::UDMcCollisionsToUDCollisions>;
130127

131128
// init
132129
void init(InitContext&)
@@ -136,6 +133,8 @@ struct TauEventTableProducer {
136133

137134
mySetITShitsRule(cutGlobalTrack.cutITShitsRule);
138135

136+
histos.add("Truth/hTroubles", "Counter of unwanted issues;;Number of troubles (-)", HistType::kTH1D, {{15, 0.5, 15.5}});
137+
139138
} // end init
140139

141140
template <typename C>
@@ -408,19 +407,19 @@ struct TauEventTableProducer {
408407
PROCESS_SWITCH(TauEventTableProducer, processDataSG, "Iterate UD tables with measured data created by SG-Candidate-Producer.", false);
409408

410409
PresliceUnsorted<aod::UDMcParticles> partPerMcCollision = aod::udmcparticle::udMcCollisionId;
411-
// PresliceUnsorted<UDMcParticlesWithUDTracks> partWtrkPerMcCollision = aod::udmcparticle::udMcCollisionId;
412410
PresliceUnsorted<FullMCSGUDCollisions> colPerMcCollision = aod::udcollision::udMcCollisionId;
413411
PresliceUnsorted<FullMCUDTracks> trackPerMcParticle = aod::udmctracklabel::udMcParticleId;
414-
PresliceUnsorted<FullMCUDTracks> trackPerCollision = aod::udtrack::udCollisionId;
412+
Preslice<FullMCUDTracks> trackPerCollision = aod::udtrack::udCollisionId; // sorted preslice used because the pair track-collision is already sorted in processDataSG function
415413

416414
void processMonteCarlo(aod::UDMcCollisions const& mccollisions,
417415
aod::UDMcParticles const& parts,
418416
FullMCSGUDCollisions const& recolls,
419417
FullMCUDTracks const& trks)
420418
{
421-
419+
// start loop over generated collisions
422420
for(const auto& mccoll : mccollisions){
423421

422+
// prepare local variables for output table
424423
int32_t runNumber = -999;
425424
int bc = -999;
426425
int nTrks[3] = {-999,-999,-999};//totalTracks, numContrib, globalNonPVtracks
@@ -465,17 +464,24 @@ struct TauEventTableProducer {
465464
float trueDaugY[2] = {-999., -999.};
466465
float trueDaugZ[2] = {-999., -999.};
467466
int trueDaugPdgCode[2] = {-999,-999};
467+
bool problem = false;
468468

469+
// find reconstructed collisions associated to the generated collision
469470
auto const& collFromMcColls = recolls.sliceBy(colPerMcCollision, mccoll.globalIndex());
470-
// if (mccoll.has_udcollisions()) {
471-
if (collFromMcColls.size() > 0) {
471+
// check the generated collision was reconstructed
472+
if (collFromMcColls.size() > 0) { // get the truth and reco-level info
472473
trueHasRecoColl = true;
473-
// auto const& collFromMcColls = mccoll.udcollisions_as<FullMCSGUDCollisions>();
474+
// check there is exactly one reco-level collision associated to generated collision
474475
if (collFromMcColls.size() > 1) {
475-
printLargeMessage("Truth collision has more than 1 reco collision. Skipping this event.");
476+
if (verboseInfo)
477+
printLargeMessage("Truth collision has more than 1 reco collision. Skipping this event.");
478+
histos.get<TH1>(HIST("Truth/hTroubles"))->Fill(1);
479+
problem = true;
476480
continue;
477481
}
482+
// grap reco-level collision
478483
auto const& collFromMcColl = collFromMcColls.iteratorAt(0);
484+
// grab tracks from the reco-level collision to get info to match measured data tables (processDataSG function)
479485
auto const& trksFromColl = trks.sliceBy(trackPerCollision, collFromMcColl.globalIndex());
480486
int countTracksPerCollision = 0;
481487
int countGoodNonPVtracks = 0;
@@ -487,6 +493,7 @@ struct TauEventTableProducer {
487493
}
488494
}
489495

496+
// fill info for reconstructed collision
490497
runNumber = collFromMcColl.runNumber();
491498
bc = collFromMcColl.globalBC();
492499
nTrks[0] = countTracksPerCollision;
@@ -513,125 +520,156 @@ struct TauEventTableProducer {
513520
timesFIT[1] = collFromMcColl.timeFT0C();
514521
timesFIT[2] = collFromMcColl.timeFV0A();
515522

523+
// get particles associated to generated collision
516524
auto const& partsFromMcColl = parts.sliceBy(partPerMcCollision, mccoll.globalIndex());
517525
int countMothers = 0;
518526
for (const auto& particle : partsFromMcColl) {
527+
// select only tauons with checking if particle has no mother
519528
if (particle.has_mothers())
520529
continue;
521530
countMothers++;
531+
// check the generated collision does not have more than 2 tauons
522532
if (countMothers > 2) {
523-
printLargeMessage("Truth collision has more than 2 no mother particles. Breaking the particle loop.");
533+
if (verboseInfo)
534+
printLargeMessage("Truth collision has more than 2 no mother particles. Breaking the particle loop.");
535+
histos.get<TH1>(HIST("Truth/hTroubles"))->Fill(2);
536+
problem = true;
524537
break;
525538
}
539+
// fill info for each tau
526540
trueTauX[countMothers-1] = particle.px();
527541
trueTauY[countMothers-1] = particle.py();
528542
trueTauZ[countMothers-1] = particle.pz();
529543

544+
// get daughters of the tau
530545
const auto& daughters = particle.daughters_as<aod::UDMcParticles>();
531546
int countDaughters = 0;
532547
for (const auto& daughter : daughters){
533548
// check if it is the charged particle (= no pi0 or neutrino)
534549
if (enumMyParticle(daughter.pdgCode()) == -1)
535550
continue;
536551
countDaughters++;
537-
if (countDaughters > 2) {
538-
printLargeMessage("Truth collision has more than 2 charged daughters of no mother particles. Breaking the daughter loop.");
552+
// check there is only 1 charged daughter related to 1 tau
553+
if (countDaughters > 1) {
554+
if (verboseInfo)
555+
printLargeMessage("Truth collision has more than 1 charged daughters of no mother particles. Breaking the daughter loop.");
556+
histos.get<TH1>(HIST("Truth/hTroubles"))->Fill(3);
557+
problem = true;
539558
break;
540559
}
541-
trueDaugX[countDaughters-1] = daughter.px();
542-
trueDaugY[countDaughters-1] = daughter.py();
543-
trueDaugZ[countDaughters-1] = daughter.pz();
544-
trueDaugPdgCode[countDaughters-1] = daughter.pdgCode();
560+
// fill info for each daughter
561+
trueDaugX[countMothers-1] = daughter.px();
562+
trueDaugY[countMothers-1] = daughter.py();
563+
trueDaugZ[countMothers-1] = daughter.pz();
564+
trueDaugPdgCode[countMothers-1] = daughter.pdgCode();
545565

546-
// const auto& tracksFromDaughter = daughter.udtracks_as<FullMCUDTracks>();
566+
// get tracks associated to MC daughter (how well the daughter was reconstructed)
547567
auto const& tracksFromDaughter = trks.sliceBy(trackPerMcParticle, daughter.globalIndex());
568+
// check there is exactly 1 track per 1 particle
548569
if (tracksFromDaughter.size() > 1) {
549-
printLargeMessage("Daughter has more than 1 associated track. Skipping this daughter.");
570+
if (verboseInfo)
571+
printLargeMessage("Daughter has more than 1 associated track. Skipping this daughter.");
572+
histos.get<TH1>(HIST("Truth/hTroubles"))->Fill(4);
573+
problem = true;
550574
continue;
551575
}
576+
// grab the track and fill info for reconstructed track (should be done twice)
552577
const auto& trk = tracksFromDaughter.iteratorAt(0);
553-
px[countDaughters-1] = trk.px();
554-
py[countDaughters-1] = trk.py();
555-
pz[countDaughters-1] = trk.pz();
556-
sign[countDaughters-1] = trk.sign();
557-
dcaxy[countDaughters-1] = trk.dcaXY();
558-
dcaz[countDaughters-1] = trk.dcaZ();
559-
trkTimeRes[countDaughters-1] = trk.trackTimeRes();
560-
if (countDaughters == 1) {
578+
px[countMothers-1] = trk.px();
579+
py[countMothers-1] = trk.py();
580+
pz[countMothers-1] = trk.pz();
581+
sign[countMothers-1] = trk.sign();
582+
dcaxy[countMothers-1] = trk.dcaXY();
583+
dcaz[countMothers-1] = trk.dcaZ();
584+
trkTimeRes[countMothers-1] = trk.trackTimeRes();
585+
if (countMothers == 1) {
561586
itsClusterSizesTrk1 = trk.itsClusterSizes();
562587
} else {
563588
itsClusterSizesTrk2 = trk.itsClusterSizes();
564589
}
565-
tpcSignal[countDaughters-1] = trk.tpcSignal();
566-
tpcEl[countDaughters-1] = trk.tpcNSigmaEl();
567-
tpcMu[countDaughters-1] = trk.tpcNSigmaMu();
568-
tpcPi[countDaughters-1] = trk.tpcNSigmaPi();
569-
tpcKa[countDaughters-1] = trk.tpcNSigmaKa();
570-
tpcPr[countDaughters-1] = trk.tpcNSigmaPr();
571-
tpcIP[countDaughters-1] = trk.tpcInnerParam();
572-
tofSignal[countDaughters-1] = trk.tofSignal();
573-
tofEl[countDaughters-1] = trk.tofNSigmaEl();
574-
tofMu[countDaughters-1] = trk.tofNSigmaMu();
575-
tofPi[countDaughters-1] = trk.tofNSigmaPi();
576-
tofKa[countDaughters-1] = trk.tofNSigmaKa();
577-
tofPr[countDaughters-1] = trk.tofNSigmaPr();
578-
tofEP[countDaughters-1] = trk.tofExpMom();
590+
tpcSignal[countMothers-1] = trk.tpcSignal();
591+
tpcEl[countMothers-1] = trk.tpcNSigmaEl();
592+
tpcMu[countMothers-1] = trk.tpcNSigmaMu();
593+
tpcPi[countMothers-1] = trk.tpcNSigmaPi();
594+
tpcKa[countMothers-1] = trk.tpcNSigmaKa();
595+
tpcPr[countMothers-1] = trk.tpcNSigmaPr();
596+
tpcIP[countMothers-1] = trk.tpcInnerParam();
597+
tofSignal[countMothers-1] = trk.tofSignal();
598+
tofEl[countMothers-1] = trk.tofNSigmaEl();
599+
tofMu[countMothers-1] = trk.tofNSigmaMu();
600+
tofPi[countMothers-1] = trk.tofNSigmaPi();
601+
tofKa[countMothers-1] = trk.tofNSigmaKa();
602+
tofPr[countMothers-1] = trk.tofNSigmaPr();
603+
tofEP[countMothers-1] = trk.tofExpMom();
579604
}// daughters
580605
}// particles
581-
} else {
606+
} else { // get only the truth information. The reco-level info is left on default
607+
// get particles associated to generated collision
582608
auto const& partsFromMcColl = parts.sliceBy(partPerMcCollision, mccoll.globalIndex());
583609
int countMothers = 0;
584610
for (const auto& particle : partsFromMcColl) {
611+
// select only tauons with checking if particle has no mother
585612
if (particle.has_mothers())
586613
continue;
587614
countMothers++;
615+
// check the generated collision does not have more than 2 tauons
588616
if (countMothers > 2) {
589-
printLargeMessage("Truth collision has more than 2 no mother particles. Breaking the particle loop.");
617+
if (verboseInfo)
618+
printLargeMessage("Truth collision has more than 2 no mother particles. Breaking the particle loop.");
619+
histos.get<TH1>(HIST("Truth/hTroubles"))->Fill(12);
620+
problem = true;
590621
break;
591622
}
623+
// fill info for each tau
592624
trueTauX[countMothers-1] = particle.px();
593625
trueTauY[countMothers-1] = particle.py();
594626
trueTauZ[countMothers-1] = particle.pz();
595627

628+
// get daughters of the tau
596629
const auto& daughters = particle.daughters_as<aod::UDMcParticles>();
597630
int countDaughters = 0;
598631
for (const auto& daughter : daughters){
599-
// check if it is the charged particle (= no pi0 or neutrino)
632+
// select only the charged particle (= no pi0 or neutrino)
600633
if (enumMyParticle(daughter.pdgCode()) == -1)
601634
continue;
602635
countDaughters++;
603-
if (countDaughters > 2) {
604-
printLargeMessage("Truth collision has more than 2 charged daughters of no mother particles. Breaking the daughter loop.");
636+
// check there is only 1 charged daughter related to 1 tau
637+
if (countDaughters > 1) {
638+
if (verboseInfo)
639+
printLargeMessage("Truth collision has more than 1 charged daughters of no mother particles. Breaking the daughter loop.");
640+
histos.get<TH1>(HIST("Truth/hTroubles"))->Fill(13);
641+
problem = true;
605642
break;
606643
}
607-
trueDaugX[countDaughters-1] = daughter.px();
608-
trueDaugY[countDaughters-1] = daughter.py();
609-
trueDaugZ[countDaughters-1] = daughter.pz();
610-
trueDaugPdgCode[countDaughters-1] = daughter.pdgCode();
644+
// fill info for each daughter
645+
trueDaugX[countMothers-1] = daughter.px();
646+
trueDaugY[countMothers-1] = daughter.py();
647+
trueDaugZ[countMothers-1] = daughter.pz();
648+
trueDaugPdgCode[countMothers-1] = daughter.pdgCode();
611649
}// daughters
612650
}// particles
613651
}// collisions
614652

653+
// decide the channel and set the variable. Only two cahnnels suported now.
654+
if ((enumMyParticle(trueDaugPdgCode[0]) == P_ELECTRON) && (enumMyParticle(trueDaugPdgCode[1]) == P_ELECTRON))
655+
trueChannel = CH_EE;
656+
if ((enumMyParticle(trueDaugPdgCode[0]) == P_ELECTRON) && ((enumMyParticle(trueDaugPdgCode[1]) == P_PION) || (enumMyParticle(trueDaugPdgCode[1]) == P_MUON)))
657+
trueChannel = CH_EMUPI;
658+
if ((enumMyParticle(trueDaugPdgCode[1]) == P_ELECTRON) && ((enumMyParticle(trueDaugPdgCode[0]) == P_PION) || (enumMyParticle(trueDaugPdgCode[0]) == P_MUON)))
659+
trueChannel = CH_EMUPI;
660+
615661
trueTauTwoTracks(runNumber, bc, nTrks[0], nTrks[1], nTrks[2], vtxPos[0], vtxPos[1], vtxPos[2],
616662
recoMode, occupancy, hadronicRate, bcSels[0], bcSels[1], bcSels[2],
617663
bcSels[3], bcSels[4], bcSels[5], bcSels[6], bcSels[7],
618-
amplitudesFIT[0], amplitudesFIT[1], amplitudesFIT[2], -999.,-999.,
619-
timesFIT[0], timesFIT[1], timesFIT[2], -999.,-999.,
664+
amplitudesFIT[0], amplitudesFIT[1], amplitudesFIT[2], -999.,-999.,// no ZDC info in MC
665+
timesFIT[0], timesFIT[1], timesFIT[2], -999.,-999.,// no ZDC info in MC
620666
px, py, pz, sign, dcaxy, dcaz, trkTimeRes,
621667
itsClusterSizesTrk1, itsClusterSizesTrk2,
622668
tpcSignal, tpcEl, tpcMu, tpcPi, tpcKa, tpcPr, tpcIP,
623669
tofSignal, tofEl, tofMu, tofPi, tofKa, tofPr, tofEP,
624670
trueChannel, trueHasRecoColl, mccoll.posX(), mccoll.posY(), mccoll.posZ(),
625-
trueTauX, trueTauY, trueTauZ, trueDaugX, trueDaugY, trueDaugZ, trueDaugPdgCode);
626-
627-
628-
// auto colSlice = recolls.sliceBy(colPerMcCollision, mccoll.globalIndex());
629-
// LOGF(info, "collision slice size %i ", colSlice.size());
671+
trueTauX, trueTauY, trueTauZ, trueDaugX, trueDaugY, trueDaugZ, trueDaugPdgCode, problem);
630672
}// mccollisions
631-
632-
633-
// const auto& mccollision = mccollisions.iteratorAt(0);
634-
635673
}
636674
PROCESS_SWITCH(TauEventTableProducer, processMonteCarlo, "Iterate UD tables with simulated data created by SG-Candidate-Producer.", false);
637675

0 commit comments

Comments
 (0)