From e4fda256ecc383189e582dee6060034ad8a953e9 Mon Sep 17 00:00:00 2001 From: Roman Lavicka Date: Wed, 18 Mar 2026 17:36:49 +0100 Subject: [PATCH] Updating the table producer for VM->twotracks Monte Carlo --- .../twoTracksEventTableProducer.cxx | 107 ++++++++++-------- 1 file changed, 61 insertions(+), 46 deletions(-) diff --git a/PWGUD/TableProducer/twoTracksEventTableProducer.cxx b/PWGUD/TableProducer/twoTracksEventTableProducer.cxx index 874761db585..3af6a964021 100644 --- a/PWGUD/TableProducer/twoTracksEventTableProducer.cxx +++ b/PWGUD/TableProducer/twoTracksEventTableProducer.cxx @@ -626,6 +626,7 @@ struct TwoTracksEventTableProducer { // get particles associated to generated collision auto const& partsFromMcColl = parts.sliceBy(partPerMcCollision, mccoll.globalIndex()); int countMothers = 0; + int totalChargedDaughters = 0; for (const auto& particle : partsFromMcColl) { // select only mothers with checking if particle has no mother if (particle.has_mothers()) @@ -644,74 +645,85 @@ struct TwoTracksEventTableProducer { trueMotherY[countMothers - 1] = particle.py(); trueMotherZ[countMothers - 1] = particle.pz(); - // get daughters of the tau + // get daughters of the mother const auto& daughters = particle.daughters_as(); - int countDaughters = 0; for (const auto& daughter : daughters) { // check if it is the charged particle (= no pi0 or neutrino) if (enumMyParticle(daughter.pdgCode()) == -1) continue; - countDaughters++; - // check there is only 1 charged daughter related to 1 tau - if (countDaughters > 1) { + // check we do not have more than 2 charged daughters in total + if (totalChargedDaughters >= 2) { if (verboseInfo) - printLargeMessage("Truth collision has more than 1 charged daughters of no mother particles. Breaking the daughter loop."); + printLargeMessage("Truth collision has more than 2 total charged daughters. Breaking the daughter loop."); histos.get(HIST("Truth/hTroubles"))->Fill(3); problem = true; break; } // fill info for each daughter - trueDaugX[countMothers - 1] = daughter.px(); - trueDaugY[countMothers - 1] = daughter.py(); - trueDaugZ[countMothers - 1] = daughter.pz(); - trueDaugPdgCode[countMothers - 1] = daughter.pdgCode(); + trueDaugX[totalChargedDaughters] = daughter.px(); + trueDaugY[totalChargedDaughters] = daughter.py(); + trueDaugZ[totalChargedDaughters] = daughter.pz(); + trueDaugPdgCode[totalChargedDaughters] = daughter.pdgCode(); // get tracks associated to MC daughter (how well the daughter was reconstructed) auto const& tracksFromDaughter = trks.sliceBy(trackPerMcParticle, daughter.globalIndex()); // check there is exactly 1 track per 1 particle + if (tracksFromDaughter.size() == 0) { + if (verboseInfo) + printLargeMessage("Daughter has no associated track. Skipping this daughter."); + histos.get(HIST("Truth/hTroubles"))->Fill(5); + problem = true; + totalChargedDaughters++; + continue; + } if (tracksFromDaughter.size() > 1) { if (verboseInfo) printLargeMessage("Daughter has more than 1 associated track. Skipping this daughter."); histos.get(HIST("Truth/hTroubles"))->Fill(4); problem = true; + totalChargedDaughters++; continue; } - // grab the track and fill info for reconstructed track (should be done twice) + // grab the track and fill info for reconstructed track const auto& trk = tracksFromDaughter.iteratorAt(0); - px[countMothers - 1] = trk.px(); - py[countMothers - 1] = trk.py(); - pz[countMothers - 1] = trk.pz(); - sign[countMothers - 1] = trk.sign(); - dcaxy[countMothers - 1] = trk.dcaXY(); - dcaz[countMothers - 1] = trk.dcaZ(); - trkTimeRes[countMothers - 1] = trk.trackTimeRes(); - if (countMothers == 1) { + px[totalChargedDaughters] = trk.px(); + py[totalChargedDaughters] = trk.py(); + pz[totalChargedDaughters] = trk.pz(); + sign[totalChargedDaughters] = trk.sign(); + dcaxy[totalChargedDaughters] = trk.dcaXY(); + dcaz[totalChargedDaughters] = trk.dcaZ(); + trkTimeRes[totalChargedDaughters] = trk.trackTimeRes(); + if (totalChargedDaughters == 0) { itsClusterSizesTrk1 = trk.itsClusterSizes(); } else { itsClusterSizesTrk2 = trk.itsClusterSizes(); } - tpcSignal[countMothers - 1] = trk.tpcSignal(); - tpcEl[countMothers - 1] = trk.tpcNSigmaEl(); - tpcMu[countMothers - 1] = trk.tpcNSigmaMu(); - tpcPi[countMothers - 1] = trk.tpcNSigmaPi(); - tpcKa[countMothers - 1] = trk.tpcNSigmaKa(); - tpcPr[countMothers - 1] = trk.tpcNSigmaPr(); - tpcIP[countMothers - 1] = trk.tpcInnerParam(); - tofSignal[countMothers - 1] = trk.tofSignal(); - tofEl[countMothers - 1] = trk.tofNSigmaEl(); - tofMu[countMothers - 1] = trk.tofNSigmaMu(); - tofPi[countMothers - 1] = trk.tofNSigmaPi(); - tofKa[countMothers - 1] = trk.tofNSigmaKa(); - tofPr[countMothers - 1] = trk.tofNSigmaPr(); - tofEP[countMothers - 1] = trk.tofExpMom(); + tpcSignal[totalChargedDaughters] = trk.tpcSignal(); + tpcEl[totalChargedDaughters] = trk.tpcNSigmaEl(); + tpcMu[totalChargedDaughters] = trk.tpcNSigmaMu(); + tpcPi[totalChargedDaughters] = trk.tpcNSigmaPi(); + tpcKa[totalChargedDaughters] = trk.tpcNSigmaKa(); + tpcPr[totalChargedDaughters] = trk.tpcNSigmaPr(); + tpcIP[totalChargedDaughters] = trk.tpcInnerParam(); + tofSignal[totalChargedDaughters] = trk.tofSignal(); + tofEl[totalChargedDaughters] = trk.tofNSigmaEl(); + tofMu[totalChargedDaughters] = trk.tofNSigmaMu(); + tofPi[totalChargedDaughters] = trk.tofNSigmaPi(); + tofKa[totalChargedDaughters] = trk.tofNSigmaKa(); + tofPr[totalChargedDaughters] = trk.tofNSigmaPr(); + tofEP[totalChargedDaughters] = trk.tofExpMom(); + totalChargedDaughters++; } // daughters + if (problem) + break; } // particles } else { // get only the truth information. The reco-level info is left on default // get particles associated to generated collision auto const& partsFromMcColl = parts.sliceBy(partPerMcCollision, mccoll.globalIndex()); int countMothers = 0; + int totalChargedDaughters = 0; for (const auto& particle : partsFromMcColl) { - // select only tauons with checking if particle has no mother + // select only motherless particles if (particle.has_mothers()) continue; countMothers++; @@ -723,39 +735,42 @@ struct TwoTracksEventTableProducer { problem = true; break; } - // fill info for each tau + // fill info for each mother trueMotherX[countMothers - 1] = particle.px(); trueMotherY[countMothers - 1] = particle.py(); trueMotherZ[countMothers - 1] = particle.pz(); - // get daughters of the tau + // get daughters of the mother const auto& daughters = particle.daughters_as(); - int countDaughters = 0; for (const auto& daughter : daughters) { // select only the charged particle (= no pi0 or neutrino) if (enumMyParticle(daughter.pdgCode()) == -1) continue; - countDaughters++; - // check there is only 1 charged daughter related to 1 tau - if (countDaughters > 1) { + // check we do not have more than 2 charged daughters in total + if (totalChargedDaughters >= 2) { if (verboseInfo) - printLargeMessage("Truth collision has more than 1 charged daughters of no mother particles. Breaking the daughter loop."); + printLargeMessage("Truth collision has more than 2 total charged daughters. Breaking the daughter loop."); histos.get(HIST("Truth/hTroubles"))->Fill(13); problem = true; break; } // fill info for each daughter - trueDaugX[countMothers - 1] = daughter.px(); - trueDaugY[countMothers - 1] = daughter.py(); - trueDaugZ[countMothers - 1] = daughter.pz(); - trueDaugPdgCode[countMothers - 1] = daughter.pdgCode(); + trueDaugX[totalChargedDaughters] = daughter.px(); + trueDaugY[totalChargedDaughters] = daughter.py(); + trueDaugZ[totalChargedDaughters] = daughter.pz(); + trueDaugPdgCode[totalChargedDaughters] = daughter.pdgCode(); + totalChargedDaughters++; } // daughters + if (problem) + break; } // particles } // collisions - // decide the channel and set the variable. Only two cahnnels suported now. + // decide the channel and set the variable if ((enumMyParticle(trueDaugPdgCode[0]) == P_ELECTRON) && (enumMyParticle(trueDaugPdgCode[1]) == P_ELECTRON)) trueChannel = CH_EE; + if ((enumMyParticle(trueDaugPdgCode[0]) == P_MUON) && (enumMyParticle(trueDaugPdgCode[1]) == P_MUON)) + trueChannel = CH_MUMU; if ((enumMyParticle(trueDaugPdgCode[0]) == P_ELECTRON) && ((enumMyParticle(trueDaugPdgCode[1]) == P_PION) || (enumMyParticle(trueDaugPdgCode[1]) == P_MUON))) trueChannel = CH_EMUPI; if ((enumMyParticle(trueDaugPdgCode[1]) == P_ELECTRON) && ((enumMyParticle(trueDaugPdgCode[0]) == P_PION) || (enumMyParticle(trueDaugPdgCode[0]) == P_MUON)))