@@ -562,62 +562,75 @@ class RecoDecay
562562 }
563563
564564 // / Finds the mother of an MC particle by looking for the expected PDG code in the mother chain.
565+ // / \param particlesMC table with MC particles
565566 // / \param particle MC particle
566567 // / \param PDGMother expected mother PDG code
567568 // / \param acceptAntiParticles switch to accept the antiparticle of the expected mother
568569 // / \param sign antiparticle indicator of the found mother w.r.t. PDGMother; 1 if particle, -1 if antiparticle, 0 if mother not found
569570 // / \param depthMax maximum decay tree level to check; Mothers up to this level will be considered. If -1, all levels are considered.
570571 // / \return index of the mother particle if found, -1 otherwise
571572 template <typename T>
572- static int getMother (const T& particle,
573+ static int getMother (const T& particlesMC,
574+ const typename T::iterator& particle,
573575 int PDGMother,
574576 bool acceptAntiParticles = false ,
575577 int8_t * sign = nullptr ,
576578 int8_t depthMax = -1 )
577579 {
578- int8_t sgn = 0 ; // 1 if the expected mother is particle, -1 if antiparticle (w.r.t. PDGMother)
579- int indexMother = -1 ; // index of the final matched mother, if found
580- auto particleMother = particle; // Initialise loop over mothers.
581- int stage = 0 ; // mother tree level (just for debugging)
580+ int8_t sgn = 0 ; // 1 if the expected mother is particle, -1 if antiparticle (w.r.t. PDGMother)
581+ int indexMother = -1 ; // index of the final matched mother, if found
582+ int stage = 0 ; // mother tree level (just for debugging)
583+ bool motherFound = false ; // true when the desired mother particle is found in the kine tree
582584 if (sign) {
583585 *sign = sgn;
584586 }
585- std::vector<int > listId;
586- listId.push_back (particleMother.globalIndex ());
587- while (particleMother.has_mothers ()) {
588- if (depthMax > -1 && -stage >= depthMax) { // Maximum depth has been reached.
589- return -1 ;
590- }
591- if (particleMother.mothersIds ()[0 ] == -1 ) { // Additional check for converted Run 2 MC
592- return -1 ;
593- }
594- particleMother = particleMother.template mothers_first_as <typename std::decay_t <T>::parent_t >(); // get mother 0
595- auto indexMotherTmp = particleMother.globalIndex ();
596- // Check mother's PDG code.
597- auto PDGParticleIMother = particleMother.pdgCode (); // PDG code of the mother
598- // printf("getMother: ");
599- // for (int i = stage; i < 0; i++) // Indent to make the tree look nice.
600- // printf(" ");
601- // printf("Stage %d: Mother PDG: %d, Index: %d\n", stage, PDGParticleIMother, indexMotherTmp);
602- if (std::find (listId.begin (), listId.end (), indexMotherTmp) != listId.end ()) {
603- LOGF (fatal, " Circular mothership at particle index %d"
604- return -1 ;
605- }
606- listId.push_back (indexMotherTmp);
607- if (PDGParticleIMother == PDGMother) { // exact PDG match
608- sgn = 1 ;
609- indexMother = indexMotherTmp;
610- break ;
611- } else if (acceptAntiParticles && PDGParticleIMother == -PDGMother) { // antiparticle PDG match
612- sgn = -1 ;
613- indexMother = indexMotherTmp;
614- break ;
587+
588+ // vector of vectors with mother indices; each line corresponds to a "stage"
589+ std::vector<std::vector<long int >> arrayIds{};
590+ std::vector<long int > initVec{particle.globalIndex ()};
591+ arrayIds.push_back (initVec); // the first vector contains the index of the original particle
592+
593+ while (!motherFound && arrayIds[-stage].size () > 0 && (depthMax < 0 || -stage < depthMax)) {
594+ // vector of mother indices for the current stage
595+ std::vector<long int > arrayIdsStage{};
596+ for (auto & iPart : arrayIds[-stage]) { // check all the particles that were the mothers at the previous stage
597+ auto particleMother = particlesMC.rawIteratorAt (iPart - particlesMC.offset ());
598+ if (particleMother.has_mothers ()) {
599+ for (auto iMother = particleMother.mothersIds ().front (); iMother <= particleMother.mothersIds ().back (); ++iMother) { // loop over the mother particles of the analysed particle
600+ if (std::find (arrayIdsStage.begin (), arrayIdsStage.end (), iMother) != arrayIdsStage.end ()) { // if a mother is still present in the vector, do not check it again
601+ continue ;
602+ }
603+ auto mother = particlesMC.rawIteratorAt (iMother - particlesMC.offset ());
604+ // Check mother's PDG code.
605+ auto PDGParticleIMother = mother.pdgCode (); // PDG code of the mother
606+ // printf("getMother: ");
607+ // for (int i = stage; i < 0; i++) // Indent to make the tree look nice.
608+ // printf(" ");
609+ // printf("Stage %d: Mother PDG: %d, Index: %d\n", stage, PDGParticleIMother, iMother);
610+ if (PDGParticleIMother == PDGMother) { // exact PDG match
611+ sgn = 1 ;
612+ indexMother = iMother;
613+ motherFound = true ;
614+ break ;
615+ } else if (acceptAntiParticles && PDGParticleIMother == -PDGMother) { // antiparticle PDG match
616+ sgn = -1 ;
617+ indexMother = iMother;
618+ motherFound = true ;
619+ break ;
620+ }
621+ // add mother index in the vector for the current stage
622+ arrayIdsStage.push_back (iMother);
623+ }
624+ }
615625 }
626+ // add vector of mother indices for the current stage
627+ arrayIds.push_back (arrayIdsStage);
616628 stage--;
617629 }
618630 if (sign) {
619631 *sign = sgn;
620632 }
633+
621634 return indexMother;
622635 }
623636
@@ -723,7 +736,7 @@ class RecoDecay
723736 if (iProng == 0 ) {
724737 // Get the mother index and its sign.
725738 // PDG code of the first daughter's mother determines whether the expected mother is a particle or antiparticle.
726- indexMother = getMother (particleI, PDGMother, acceptAntiParticles, &sgn, depthMax);
739+ indexMother = getMother (particlesMC, particleI, PDGMother, acceptAntiParticles, &sgn, depthMax);
727740 // Check whether mother was found.
728741 if (indexMother <= -1 ) {
729742 // Printf("MC Rec: Rejected: bad mother index or PDG");
0 commit comments