@@ -626,6 +626,7 @@ struct TwoTracksEventTableProducer {
626626 // get particles associated to generated collision
627627 auto const & partsFromMcColl = parts.sliceBy (partPerMcCollision, mccoll.globalIndex ());
628628 int countMothers = 0 ;
629+ int totalChargedDaughters = 0 ;
629630 for (const auto & particle : partsFromMcColl) {
630631 // select only mothers with checking if particle has no mother
631632 if (particle.has_mothers ())
@@ -644,74 +645,85 @@ struct TwoTracksEventTableProducer {
644645 trueMotherY[countMothers - 1 ] = particle.py ();
645646 trueMotherZ[countMothers - 1 ] = particle.pz ();
646647
647- // get daughters of the tau
648+ // get daughters of the mother
648649 const auto & daughters = particle.daughters_as <aod::UDMcParticles>();
649- int countDaughters = 0 ;
650650 for (const auto & daughter : daughters) {
651651 // check if it is the charged particle (= no pi0 or neutrino)
652652 if (enumMyParticle (daughter.pdgCode ()) == -1 )
653653 continue ;
654- countDaughters++;
655- // check there is only 1 charged daughter related to 1 tau
656- if (countDaughters > 1 ) {
654+ // check we do not have more than 2 charged daughters in total
655+ if (totalChargedDaughters >= 2 ) {
657656 if (verboseInfo)
658- printLargeMessage (" Truth collision has more than 1 charged daughters of no mother particles . Breaking the daughter loop." );
657+ printLargeMessage (" Truth collision has more than 2 total charged daughters. Breaking the daughter loop." );
659658 histos.get <TH1>(HIST (" Truth/hTroubles" ))->Fill (3 );
660659 problem = true ;
661660 break ;
662661 }
663662 // fill info for each daughter
664- trueDaugX[countMothers - 1 ] = daughter.px ();
665- trueDaugY[countMothers - 1 ] = daughter.py ();
666- trueDaugZ[countMothers - 1 ] = daughter.pz ();
667- trueDaugPdgCode[countMothers - 1 ] = daughter.pdgCode ();
663+ trueDaugX[totalChargedDaughters ] = daughter.px ();
664+ trueDaugY[totalChargedDaughters ] = daughter.py ();
665+ trueDaugZ[totalChargedDaughters ] = daughter.pz ();
666+ trueDaugPdgCode[totalChargedDaughters ] = daughter.pdgCode ();
668667
669668 // get tracks associated to MC daughter (how well the daughter was reconstructed)
670669 auto const & tracksFromDaughter = trks.sliceBy (trackPerMcParticle, daughter.globalIndex ());
671670 // check there is exactly 1 track per 1 particle
671+ if (tracksFromDaughter.size () == 0 ) {
672+ if (verboseInfo)
673+ printLargeMessage (" Daughter has no associated track. Skipping this daughter." );
674+ histos.get <TH1>(HIST (" Truth/hTroubles" ))->Fill (5 );
675+ problem = true ;
676+ totalChargedDaughters++;
677+ continue ;
678+ }
672679 if (tracksFromDaughter.size () > 1 ) {
673680 if (verboseInfo)
674681 printLargeMessage (" Daughter has more than 1 associated track. Skipping this daughter." );
675682 histos.get <TH1>(HIST (" Truth/hTroubles" ))->Fill (4 );
676683 problem = true ;
684+ totalChargedDaughters++;
677685 continue ;
678686 }
679- // grab the track and fill info for reconstructed track (should be done twice)
687+ // grab the track and fill info for reconstructed track
680688 const auto & trk = tracksFromDaughter.iteratorAt (0 );
681- px[countMothers - 1 ] = trk.px ();
682- py[countMothers - 1 ] = trk.py ();
683- pz[countMothers - 1 ] = trk.pz ();
684- sign[countMothers - 1 ] = trk.sign ();
685- dcaxy[countMothers - 1 ] = trk.dcaXY ();
686- dcaz[countMothers - 1 ] = trk.dcaZ ();
687- trkTimeRes[countMothers - 1 ] = trk.trackTimeRes ();
688- if (countMothers == 1 ) {
689+ px[totalChargedDaughters ] = trk.px ();
690+ py[totalChargedDaughters ] = trk.py ();
691+ pz[totalChargedDaughters ] = trk.pz ();
692+ sign[totalChargedDaughters ] = trk.sign ();
693+ dcaxy[totalChargedDaughters ] = trk.dcaXY ();
694+ dcaz[totalChargedDaughters ] = trk.dcaZ ();
695+ trkTimeRes[totalChargedDaughters ] = trk.trackTimeRes ();
696+ if (totalChargedDaughters == 0 ) {
689697 itsClusterSizesTrk1 = trk.itsClusterSizes ();
690698 } else {
691699 itsClusterSizesTrk2 = trk.itsClusterSizes ();
692700 }
693- tpcSignal[countMothers - 1 ] = trk.tpcSignal ();
694- tpcEl[countMothers - 1 ] = trk.tpcNSigmaEl ();
695- tpcMu[countMothers - 1 ] = trk.tpcNSigmaMu ();
696- tpcPi[countMothers - 1 ] = trk.tpcNSigmaPi ();
697- tpcKa[countMothers - 1 ] = trk.tpcNSigmaKa ();
698- tpcPr[countMothers - 1 ] = trk.tpcNSigmaPr ();
699- tpcIP[countMothers - 1 ] = trk.tpcInnerParam ();
700- tofSignal[countMothers - 1 ] = trk.tofSignal ();
701- tofEl[countMothers - 1 ] = trk.tofNSigmaEl ();
702- tofMu[countMothers - 1 ] = trk.tofNSigmaMu ();
703- tofPi[countMothers - 1 ] = trk.tofNSigmaPi ();
704- tofKa[countMothers - 1 ] = trk.tofNSigmaKa ();
705- tofPr[countMothers - 1 ] = trk.tofNSigmaPr ();
706- tofEP[countMothers - 1 ] = trk.tofExpMom ();
701+ tpcSignal[totalChargedDaughters] = trk.tpcSignal ();
702+ tpcEl[totalChargedDaughters] = trk.tpcNSigmaEl ();
703+ tpcMu[totalChargedDaughters] = trk.tpcNSigmaMu ();
704+ tpcPi[totalChargedDaughters] = trk.tpcNSigmaPi ();
705+ tpcKa[totalChargedDaughters] = trk.tpcNSigmaKa ();
706+ tpcPr[totalChargedDaughters] = trk.tpcNSigmaPr ();
707+ tpcIP[totalChargedDaughters] = trk.tpcInnerParam ();
708+ tofSignal[totalChargedDaughters] = trk.tofSignal ();
709+ tofEl[totalChargedDaughters] = trk.tofNSigmaEl ();
710+ tofMu[totalChargedDaughters] = trk.tofNSigmaMu ();
711+ tofPi[totalChargedDaughters] = trk.tofNSigmaPi ();
712+ tofKa[totalChargedDaughters] = trk.tofNSigmaKa ();
713+ tofPr[totalChargedDaughters] = trk.tofNSigmaPr ();
714+ tofEP[totalChargedDaughters] = trk.tofExpMom ();
715+ totalChargedDaughters++;
707716 } // daughters
717+ if (problem)
718+ break ;
708719 } // particles
709720 } else { // get only the truth information. The reco-level info is left on default
710721 // get particles associated to generated collision
711722 auto const & partsFromMcColl = parts.sliceBy (partPerMcCollision, mccoll.globalIndex ());
712723 int countMothers = 0 ;
724+ int totalChargedDaughters = 0 ;
713725 for (const auto & particle : partsFromMcColl) {
714- // select only tauons with checking if particle has no mother
726+ // select only motherless particles
715727 if (particle.has_mothers ())
716728 continue ;
717729 countMothers++;
@@ -723,39 +735,42 @@ struct TwoTracksEventTableProducer {
723735 problem = true ;
724736 break ;
725737 }
726- // fill info for each tau
738+ // fill info for each mother
727739 trueMotherX[countMothers - 1 ] = particle.px ();
728740 trueMotherY[countMothers - 1 ] = particle.py ();
729741 trueMotherZ[countMothers - 1 ] = particle.pz ();
730742
731- // get daughters of the tau
743+ // get daughters of the mother
732744 const auto & daughters = particle.daughters_as <aod::UDMcParticles>();
733- int countDaughters = 0 ;
734745 for (const auto & daughter : daughters) {
735746 // select only the charged particle (= no pi0 or neutrino)
736747 if (enumMyParticle (daughter.pdgCode ()) == -1 )
737748 continue ;
738- countDaughters++;
739- // check there is only 1 charged daughter related to 1 tau
740- if (countDaughters > 1 ) {
749+ // check we do not have more than 2 charged daughters in total
750+ if (totalChargedDaughters >= 2 ) {
741751 if (verboseInfo)
742- printLargeMessage (" Truth collision has more than 1 charged daughters of no mother particles . Breaking the daughter loop." );
752+ printLargeMessage (" Truth collision has more than 2 total charged daughters. Breaking the daughter loop." );
743753 histos.get <TH1>(HIST (" Truth/hTroubles" ))->Fill (13 );
744754 problem = true ;
745755 break ;
746756 }
747757 // fill info for each daughter
748- trueDaugX[countMothers - 1 ] = daughter.px ();
749- trueDaugY[countMothers - 1 ] = daughter.py ();
750- trueDaugZ[countMothers - 1 ] = daughter.pz ();
751- trueDaugPdgCode[countMothers - 1 ] = daughter.pdgCode ();
758+ trueDaugX[totalChargedDaughters] = daughter.px ();
759+ trueDaugY[totalChargedDaughters] = daughter.py ();
760+ trueDaugZ[totalChargedDaughters] = daughter.pz ();
761+ trueDaugPdgCode[totalChargedDaughters] = daughter.pdgCode ();
762+ totalChargedDaughters++;
752763 } // daughters
764+ if (problem)
765+ break ;
753766 } // particles
754767 } // collisions
755768
756- // decide the channel and set the variable. Only two cahnnels suported now.
769+ // decide the channel and set the variable
757770 if ((enumMyParticle (trueDaugPdgCode[0 ]) == P_ELECTRON) && (enumMyParticle (trueDaugPdgCode[1 ]) == P_ELECTRON))
758771 trueChannel = CH_EE;
772+ if ((enumMyParticle (trueDaugPdgCode[0 ]) == P_MUON) && (enumMyParticle (trueDaugPdgCode[1 ]) == P_MUON))
773+ trueChannel = CH_MUMU;
759774 if ((enumMyParticle (trueDaugPdgCode[0 ]) == P_ELECTRON) && ((enumMyParticle (trueDaugPdgCode[1 ]) == P_PION) || (enumMyParticle (trueDaugPdgCode[1 ]) == P_MUON)))
760775 trueChannel = CH_EMUPI;
761776 if ((enumMyParticle (trueDaugPdgCode[1 ]) == P_ELECTRON) && ((enumMyParticle (trueDaugPdgCode[0 ]) == P_PION) || (enumMyParticle (trueDaugPdgCode[0 ]) == P_MUON)))
0 commit comments