[Common] Fix matching to MC for daughters interacting with material (#9391)

Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
Co-authored-by: Vít Kučera <vit.kucera@cern.ch>

Author: 3 people

Common/Core/RecoDecay.h

@@ -668,9 +668,11 @@ struct RecoDecay {
   }
 
   /// Checks whether the reconstructed decay candidate is the expected decay.
-  /// \param checkProcess  switch to accept only decay daughters by checking the production process of MC particles
-  /// \param acceptIncompleteReco  switch to accept candidates with only part of the daughters reconstructed
+  /// \tparam acceptFlavourOscillation  switch to accept flavour oscillastion (i.e. B0 -> B0bar -> D+pi-)
+  /// \tparam checkProcess  switch to accept only decay daughters by checking the production process of MC particles
+  /// \tparam acceptIncompleteReco  switch to accept candidates with only part of the daughters reconstructed
   /// \tparam acceptTrackDecay  switch to accept candidates with daughter tracks of pions and kaons which decayed
+  /// \tparam acceptTrackIntWithMaterial switch to accept candidates with final (i.e. p, K, pi) daughter tracks interacting with material
   /// \param particlesMC  table with MC particles
   /// \param arrDaughters  array of candidate daughters
   /// \param PDGMother  expected mother PDG code
@@ -680,8 +682,9 @@ struct RecoDecay {
   /// \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,16 +693,18 @@ struct RecoDecay {
                              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
-    int8_t sgn = 0;                        // 1 if the expected mother is particle, -1 if antiparticle (w.r.t. PDGMother)
-    int8_t nPiToMuLocal = 0;               // number of pion prongs decayed to a muon
-    int8_t nKaToPiLocal = 0;               // number of kaon prongs decayed to a pion
-    int indexMother = -1;                  // index of the mother particle
-    std::vector<int> arrAllDaughtersIndex; // vector of indices of all daughters of the mother of the first provided daughter
-    std::array<int, N> arrDaughtersIndex;  // array of indices of provided daughters
+    int8_t coefFlavourOscillation = 1;         // 1 if no B0(s) flavour oscillation occured, -1 else
+    int8_t sgn = 0;                            // 1 if the expected mother is particle, -1 if antiparticle (w.r.t. PDGMother)
+    int8_t nPiToMuLocal = 0;                   // number of pion prongs decayed to a muon
+    int8_t nKaToPiLocal = 0;                   // number of kaon prongs decayed to a pion
+    int8_t nInteractionsWithMaterialLocal = 0; // number of interactions with material
+    int indexMother = -1;                      // index of the mother particle
+    std::vector<int> arrAllDaughtersIndex;     // vector of indices of all daughters of the mother of the first provided daughter
+    std::array<int, N> arrDaughtersIndex;      // array of indices of provided daughters
     if (sign) {
       *sign = sgn;
     }
@@ -737,6 +742,28 @@ struct RecoDecay {
           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 != TMCProcess::kPDecay && process != TMCProcess::kPPrimary && pdgI == pdgMotherI) {
+          if (!particleI.has_mothers()) {
+            break;
+          }
+          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 == TMCProcess::kPDecay || process == TMCProcess::kPPrimary) { // we found the original daughter that interacted with material
+              nInteractionsWithMaterialLocal++;
+            }
+          }
+        }
+      }
       arrDaughtersIndex[iProng] = particleI.globalIndex();
       // Get the list of daughter indices from the mother of the first prong.
       if (iProng == 0) {
@@ -817,6 +844,11 @@ struct RecoDecay {
         *nKaToPi = nKaToPiLocal;
       }
     }
+    if constexpr (acceptTrackIntWithMaterial) {
+      if (nInteractionsWithMaterial) {
+        *nInteractionsWithMaterial = nInteractionsWithMaterialLocal;
+      }
+    }
     return indexMother;
   }