@@ -668,9 +668,11 @@ struct RecoDecay {
668668 }
669669
670670 // / Checks whether the reconstructed decay candidate is the expected decay.
671- // / \param checkProcess switch to accept only decay daughters by checking the production process of MC particles
672- // / \param acceptIncompleteReco switch to accept candidates with only part of the daughters reconstructed
671+ // / \tparam acceptFlavourOscillation switch to accept flavour oscillastion (i.e. B0 -> B0bar -> D+pi-)
672+ // / \tparam checkProcess switch to accept only decay daughters by checking the production process of MC particles
673+ // / \tparam acceptIncompleteReco switch to accept candidates with only part of the daughters reconstructed
673674 // / \tparam acceptTrackDecay switch to accept candidates with daughter tracks of pions and kaons which decayed
675+ // / \tparam acceptTrackIntWithMaterial switch to accept candidates with final (i.e. p, K, pi) daughter tracks interacting with material
674676 // / \param particlesMC table with MC particles
675677 // / \param arrDaughters array of candidate daughters
676678 // / \param PDGMother expected mother PDG code
@@ -680,8 +682,9 @@ struct RecoDecay {
680682 // / \param depthMax maximum decay tree level to check; Daughters up to this level will be considered. If -1, all levels are considered.
681683 // / \param nPiToMu number of pion prongs decayed to a muon
682684 // / \param nKaToPi number of kaon prongs decayed to a pion
685+ // / \param nInteractionsWithMaterial number of daughter particles that interacted with material
683686 // / \return index of the mother particle if the mother and daughters are correct, -1 otherwise
684- template <bool acceptFlavourOscillation = false , bool checkProcess = false , bool acceptIncompleteReco = false , bool acceptTrackDecay = false , std::size_t N, typename T, typename U>
687+ template <bool acceptFlavourOscillation = false , bool checkProcess = false , bool acceptIncompleteReco = false , bool acceptTrackDecay = false , bool acceptTrackIntWithMaterial = false , std::size_t N, typename T, typename U>
685688 static int getMatchedMCRec (const T& particlesMC,
686689 const std::array<U, N>& arrDaughters,
687690 int PDGMother,
@@ -690,16 +693,18 @@ struct RecoDecay {
690693 int8_t * sign = nullptr ,
691694 int depthMax = 1 ,
692695 int8_t * nPiToMu = nullptr ,
693- int8_t * nKaToPi = nullptr )
696+ int8_t * nKaToPi = nullptr ,
697+ int8_t * nInteractionsWithMaterial = nullptr )
694698 {
695699 // Printf("MC Rec: Expected mother PDG: %d", PDGMother);
696- int8_t coefFlavourOscillation = 1 ; // 1 if no B0(s) flavour oscillation occured, -1 else
697- int8_t sgn = 0 ; // 1 if the expected mother is particle, -1 if antiparticle (w.r.t. PDGMother)
698- int8_t nPiToMuLocal = 0 ; // number of pion prongs decayed to a muon
699- int8_t nKaToPiLocal = 0 ; // number of kaon prongs decayed to a pion
700- int indexMother = -1 ; // index of the mother particle
701- std::vector<int > arrAllDaughtersIndex; // vector of indices of all daughters of the mother of the first provided daughter
702- std::array<int , N> arrDaughtersIndex; // array of indices of provided daughters
700+ int8_t coefFlavourOscillation = 1 ; // 1 if no B0(s) flavour oscillation occured, -1 else
701+ int8_t sgn = 0 ; // 1 if the expected mother is particle, -1 if antiparticle (w.r.t. PDGMother)
702+ int8_t nPiToMuLocal = 0 ; // number of pion prongs decayed to a muon
703+ int8_t nKaToPiLocal = 0 ; // number of kaon prongs decayed to a pion
704+ int8_t nInteractionsWithMaterialLocal = 0 ; // number of interactions with material
705+ int indexMother = -1 ; // index of the mother particle
706+ std::vector<int > arrAllDaughtersIndex; // vector of indices of all daughters of the mother of the first provided daughter
707+ std::array<int , N> arrDaughtersIndex; // array of indices of provided daughters
703708 if (sign) {
704709 *sign = sgn;
705710 }
@@ -737,6 +742,28 @@ struct RecoDecay {
737742 particleI = motherI;
738743 }
739744 }
745+ if constexpr (acceptTrackIntWithMaterial) {
746+ // Replace the MC particle associated with the prong by its mother for part → part due to material interactions.
747+ // It keeps looking at the mother iteratively, until it finds a particle from decay or primary
748+ auto process = particleI.getProcess ();
749+ auto pdgI = std::abs (particleI.pdgCode ());
750+ auto pdgMotherI = std::abs (particleI.pdgCode ());
751+ while (process != TMCProcess::kPDecay && process != TMCProcess::kPPrimary && pdgI == pdgMotherI) {
752+ if (!particleI.has_mothers ()) {
753+ break ;
754+ }
755+ auto motherI = particleI.template mothers_first_as <T>();
756+ pdgI = std::abs (particleI.pdgCode ());
757+ pdgMotherI = std::abs (motherI.pdgCode ());
758+ if (pdgI == pdgMotherI) {
759+ particleI = motherI;
760+ process = particleI.getProcess ();
761+ if (process == TMCProcess::kPDecay || process == TMCProcess::kPPrimary ) { // we found the original daughter that interacted with material
762+ nInteractionsWithMaterialLocal++;
763+ }
764+ }
765+ }
766+ }
740767 arrDaughtersIndex[iProng] = particleI.globalIndex ();
741768 // Get the list of daughter indices from the mother of the first prong.
742769 if (iProng == 0 ) {
@@ -817,6 +844,11 @@ struct RecoDecay {
817844 *nKaToPi = nKaToPiLocal;
818845 }
819846 }
847+ if constexpr (acceptTrackIntWithMaterial) {
848+ if (nInteractionsWithMaterial) {
849+ *nInteractionsWithMaterial = nInteractionsWithMaterialLocal;
850+ }
851+ }
820852 return indexMother;
821853 }
822854
0 commit comments