New method to select primaries

Author: rliotino99

DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx

@@ -34,6 +34,7 @@
 #include <Framework/runDataProcessing.h>
 
 #include <TGeoManager.h>
+#include <TPDGCode.h>
 
 #include <HMPIDBase/Param.h>
 
@@ -73,6 +74,12 @@ struct HmpidTableProducer {
   Configurable<bool> requireTPC{"requireTPC", true, "Require TPC track"};
   Configurable<bool> requireTOF{"requireTOF", true, "Require TOF track"};
 
+  Configurable<bool> useInAbsorberGeomMethod{"useInAbsorberGeomMethod", false, "Use geometrical method to check if daughters are born in absorber"};
+
+  // (reference) 473 cm - was the legacy value in run2 simulation
+  Configurable<float> survivalThresholdRich2{"survivalThresholdRich2", 437.5f, "survivalThresholdRich2"};
+  Configurable<float> survivalThresholdRich4{"survivalThresholdRich4", 439.0f, "survivalThresholdRich4"};
+
   using CollisionCandidates = o2::soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFV0As>;
 
   using TrackCandidates = soa::Join<aod::Tracks, aod::TracksExtra,
@@ -94,27 +101,7 @@ struct HmpidTableProducer {
 
   static constexpr int Rich2 = 2, Rich4 = 4;
 
-  // -----------------------------------------------------------------------
-  // HMPID absorber geometry (hardcoded from HMPIDSimulation/Detector.cxx,
-  // Detector::ConstructGeometry / Detector::createAbsorber).
-  // The experiment is finalised and this geometry will not change, so the
-  // values are copied here instead of being re-derived from TGeoManager/CCDB
-  // at runtime. If the detector geometry code is ever revisited, these
-  // constants must be updated accordingly.
-  //
-  // Each absorber is a box (TGeoBBox) whose LOCAL->GLOBAL transform is built as:
-  //   pMatrix->SetTranslation(T);
-  //   pMatrix->RotateZ(theta);
-  // which yields, for a local point p: p_glob = Rz(theta) * p_loc + T.
-  // In particular the box CENTER in global coordinates is exactly T (the
-  // rotation does not affect T, since it is applied to p_loc only, not to
-  // the already-set translation). Only the box AXES are rotated by theta
-  // with respect to the global x,y axes.
-  //
-  // To test whether a global point lies inside the box we invert the
-  // transform: p_loc = Rz(-theta) * (p_glob - T), then compare component-wise
-  // against the box half-widths.
-  // -----------------------------------------------------------------------
+  // (reference) HMPID Detector class in O2
   static constexpr double AbsThetaDeg = 33.5;
   const double mAbsCosT = std::cos(AbsThetaDeg * TMath::DegToRad());
   const double mAbsSinT = std::sin(AbsThetaDeg * TMath::DegToRad());
@@ -159,6 +146,10 @@ struct HmpidTableProducer {
                kTH1F, {{4, -0.5, 3.5, ""}});
 
     histos.add("hProdVertex", ";X (cm);Y (cm);Z (cm)", HistType::kTH3F, {{500, -500., 500.}, {500, -500., 500.}, {500, -500., 500.}});
+    histos.add("hDaughterRCyl_Rich2", "hDaughterRCyl_Rich2", kTH1F, {{600, 0., 600.}});
+    histos.add("hDaughterRCyl_Rich4", "hDaughterRCyl_Rich4", kTH1F, {{600, 0., 600.}});
+    histos.add("hDaughterRSph_Rich2", "hDaughterRSph_Rich2", kTH1F, {{600, 0., 600.}});
+    histos.add("hDaughterRSph_Rich4", "hDaughterRSph_Rich4", kTH1F, {{600, 0., 600.}});
   }
 
   // -----------------------------------------------------------------------
@@ -361,17 +352,32 @@ struct HmpidTableProducer {
       return false;
     }
 
-    // subtract the box center (translation is not rotated, see geometry block above)
-    const double rx = vx * mAbsCosT + vy * mAbsSinT;
-    const double ry = -vx * mAbsSinT + vy * mAbsCosT;
-    const double rz = vz;
+    // translate to box center
+    const double lx = vx - centerX;
+    const double ly = vy; // centerY = 0
+    const double lz = vz - centerZ;
 
     // rotate by -theta into the box local frame
-    const double lx = rx - centerX;
-    const double ly = ry; // centerY = 0
-    const double lz = rz - centerZ;
+    const double rx = lx * mAbsCosT + ly * mAbsSinT;
+    const double ry = -lx * mAbsSinT + ly * mAbsCosT;
+    const double rz = lz;
+
+    return std::abs(rx) <= halfX && std::abs(ry) <= AbsHalfY && std::abs(rz) <= AbsHalfZ;
+  }
+
+  bool survivedAbsorber(double vx, double vy, int chamber)
+  {
+    float thresholdR = 0.;
+    if (chamber == Rich2) {
+      thresholdR = survivalThresholdRich2;
+    } else if (chamber == Rich4) {
+      thresholdR = survivalThresholdRich4;
+    } else {
+      return false;
+    }
 
-    return std::abs(lx) <= halfX && std::abs(ly) <= AbsHalfY && std::abs(lz) <= AbsHalfZ;
+    const float r = std::hypot(vx, vy);
+    return r > thresholdR;
   }
 
   void processEvent(CollisionCandidates::iterator const& col,
@@ -514,23 +520,76 @@ struct HmpidTableProducer {
 
           if ((chamberM3 == Rich2 || chamberM3 == Rich4) && mc.has_daughters()) {
             auto dIds = mc.daughtersIds();
-
-            for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) {
-              auto daughter = mcParticles.rawIteratorAt(idx);
-
-              histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz());
-
-              if (isInAbsorber(daughter.vx(), daughter.vy(), daughter.vz(), chamberM3)) {
-                interactionInAbsorber = true;
-                break;
-              }
-            } // end loop daughters
+            bool foundRelevantDaughter = false; // true if at least one non-delta/photon daughter was examined
+
+            if (useInAbsorberGeomMethod) {
+              for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) {
+                auto daughter = mcParticles.rawIteratorAt(idx);
+
+                int absPdg = std::abs(daughter.pdgCode());
+                if (absPdg == kElectron || absPdg == kGamma)
+                  continue;
+
+                foundRelevantDaughter = true;
+
+                // diagnostics on daughters distribution
+                histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz());
+
+                double rCyl = std::hypot(daughter.vx(), daughter.vy());
+                double rSph = std::hypot(daughter.vx(), daughter.vy(), daughter.vz());
+                if (chamberM3 == Rich2) {
+                  histos.fill(HIST("hDaughterRCyl_Rich2"), rCyl);
+                  histos.fill(HIST("hDaughterRSph_Rich2"), rSph);
+                } else {
+                  histos.fill(HIST("hDaughterRCyl_Rich4"), rCyl);
+                  histos.fill(HIST("hDaughterRSph_Rich4"), rSph);
+                }
+
+                if (isInAbsorber(daughter.vx(), daughter.vy(), daughter.vz(), chamberM3)) {
+                  interactionInAbsorber = true;
+                }
+              } // end loop daughters
+            } else {
+              bool survived = false;
+              for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) {
+                auto daughter = mcParticles.rawIteratorAt(idx);
+
+                // skip delta rays (e-/e+), photons, and HMPID Cherenkov/feedback
+                int absPdg = std::abs(daughter.pdgCode());
+                if (absPdg == kElectron || absPdg == kGamma)
+                  continue;
+
+                foundRelevantDaughter = true;
+
+                histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz());
+
+                double rCyl = std::hypot(daughter.vx(), daughter.vy());
+                double rSph = std::hypot(daughter.vx(), daughter.vy(), daughter.vz());
+                if (chamberM3 == Rich2) {
+                  histos.fill(HIST("hDaughterRCyl_Rich2"), rCyl);
+                  histos.fill(HIST("hDaughterRSph_Rich2"), rSph);
+                } else {
+                  histos.fill(HIST("hDaughterRCyl_Rich4"), rCyl);
+                  histos.fill(HIST("hDaughterRSph_Rich4"), rSph);
+                }
+
+                if (survivedAbsorber(daughter.vx(), daughter.vy(), chamberM3)) {
+                  survived = true;
+                }
+              } // end loop daughters
+
+              // No relevant daughter found (only delta rays/photons, or no
+              // daughters at all): no evidence of a genuine interaction -> treat
+              // as primary/survived, consistent with the "no daughters" case.
+              interactionInAbsorber = foundRelevantDaughter ? !survived : false;
+            }
           } // end if has_daughters
 
           hmpidAnalysisMC(mc.pdgCode(), mc.vx(), mc.vy(), mc.vz(),
                           mc.isPhysicalPrimary(), mc.getProcess(), interactionInAbsorber);
         } else {
-          hmpidAnalysisMC(-1, 0.f, 0.f, 0.f, false, -100, false);
+          // No MC truth associated to this track
+          hmpidAnalysisMC(-999, -999.f, -999.f, -999.f, false, -100, false);
         }
       } // end if constexpr (isMC)