[Common] Fix matching to MC for daughters interacting with material

Common/Core/RecoDecay.h

@@ -42,7 +42,7 @@
                     Prompt,
                     NonPrompt };
 
   static constexpr int8_t PdgStatusCodeAfterFlavourOscillation = 92; // decay products after B0(s) flavour oscillation
 
   // Auxiliary functions
 
@@ -680,8 +680,9 @@
   /// \param depthMax  maximum decay tree level to check; Daughters up to this level will be considered. If -1, all levels are considered.
   /// \param nPiToMu  number of pion prongs decayed to a muon
   /// \param nKaToPi  number of kaon prongs decayed to a pion
+  /// \param nInteractionsWithMaterial  number of daughter particles that interacted with material
   /// \return index of the mother particle if the mother and daughters are correct, -1 otherwise
-  template <bool acceptFlavourOscillation = false, bool checkProcess = false, bool acceptIncompleteReco = false, bool acceptTrackDecay = false, std::size_t N, typename T, typename U>
+  template <bool acceptFlavourOscillation = false, bool checkProcess = false, bool acceptIncompleteReco = false, bool acceptTrackDecay = false, bool acceptTrackIntWithMaterial = false, std::size_t N, typename T, typename U>
   static int getMatchedMCRec(const T& particlesMC,
                              const std::array<U, N>& arrDaughters,
                              int PDGMother,
@@ -690,7 +691,8 @@
                              int8_t* sign = nullptr,
                              int depthMax = 1,
                              int8_t* nPiToMu = nullptr,
-                             int8_t* nKaToPi = nullptr)
+                             int8_t* nKaToPi = nullptr,
+                             int8_t* nInteractionsWithMaterial = nullptr)
   {
     // Printf("MC Rec: Expected mother PDG: %d", PDGMother);
     int8_t coefFlavourOscillation = 1;     // 1 if no B0(s) flavour oscillation occured, -1 else
@@ -737,6 +739,25 @@
           particleI = motherI;
         }
       }
+      if constexpr (acceptTrackIntWithMaterial) {
+        // Replace the MC particle associated with the prong by its mother for part → part due to material interactions.
+        // It keeps looking at the mother iteratively, until it finds a particle from decay or primary
+        auto process = particleI.getProcess();
+        auto pdgI = std::abs(particleI.pdgCode());
+        auto pdgMotherI = std::abs(particleI.pdgCode());
+        while (process != kPDecay && process != kPPrimary && pdgI == pdgMotherI) {
+          auto motherI = particleI.template mothers_first_as<T>();
+          pdgI = std::abs(particleI.pdgCode());
+          pdgMotherI = std::abs(motherI.pdgCode());
+          if (pdgI == pdgMotherI) {
+            particleI = motherI;
+            process = particleI.getProcess();
+            if (process == kPDecay || kPPrimary) { // we found the original daughter that interacted with material
+              nInteractionsWithMaterial++;
+            }
+          }
+        }
+      }
       arrDaughtersIndex[iProng] = particleI.globalIndex();
       // Get the list of daughter indices from the mother of the first prong.
       if (iProng == 0) {
@@ -1380,15 +1401,15 @@
     auto vecXYZFromVec = pVector(vecPtEtaPhi);
     auto vecXYZFromVars = pVector(pt(vecPtEtaPhi), eta(vecPtEtaPhi), phi(vecPtEtaPhi));
     // Test px, py, pz, pt, p, eta, phi, e, y, m
     printf("RecoDecay test\n");
     printf("px: In: %g, XYZ: %g, XYZ from vec: %g, XYZ from vars: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", pxIn, vecXYZ[0], vecXYZFromVec[0], vecXYZFromVars[0], px(vecPtEtaPhi), px(vecPtEtaPhiM));
     printf("py: In: %g, XYZ: %g, XYZ from vec: %g, XYZ from vars: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", pyIn, vecXYZ[1], vecXYZFromVec[1], vecXYZFromVars[1], py(vecPtEtaPhi), py(vecPtEtaPhiM));
     printf("pz: In: %g, XYZ: %g, XYZ from vec: %g, XYZ from vars: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", pzIn, vecXYZ[2], vecXYZFromVec[2], vecXYZFromVars[2], pz(vecPtEtaPhi), pz(vecPtEtaPhiM));
     printf("pt: XYZ: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", RecoDecay::pt(vecXYZ), pt(vecPtEtaPhi), pt(vecPtEtaPhiM));
     printf("p: XYZ: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", RecoDecay::p(vecXYZ), p(vecPtEtaPhi), p(vecPtEtaPhiM));
     printf("eta: XYZ: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", RecoDecay::eta(vecXYZ), eta(vecPtEtaPhi), eta(vecPtEtaPhiM));
     printf("phi: XYZ: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", RecoDecay::phi(vecXYZ), phi(vecPtEtaPhi), phi(vecPtEtaPhiM));
     printf("e: XYZ: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", RecoDecay::e(vecXYZ, mIn), e(vecPtEtaPhi, mIn), e(vecPtEtaPhiM));
     printf("y: XYZ: %g, PtEtaPhi: %g, PtEtaPhiM: %g\n", RecoDecay::y(vecXYZ, mIn), y(vecPtEtaPhi, mIn), y(vecPtEtaPhiM));
     printf("m: In: %g, XYZ(p, E): %g, XYZ(pVec, E): %g, XYZ(arr): %g, PtEtaPhiM: %g\n",
            mIn,