feat: apply dRICH and pfRICH quantum efficiency at photon generation (stacking) stage

src/dd4pod/python/npsim.py

@@ -13,6 +13,7 @@
 
 from DDSim.DD4hepSimulation import DD4hepSimulation
 
+from g4units import mm
 
 if __name__ == "__main__":
   logging.basicConfig(format='%(name)-16s %(levelname)s %(message)s', level=logging.INFO, stream=sys.stdout)
@@ -69,8 +70,37 @@ def setupCerenkov(kernel):
   RUNNER.action.mapActions['PFRICH'] = 'Geant4OpticalTrackerAction'
   RUNNER.action.mapActions['DIRC'] = 'Geant4OpticalTrackerAction'
 
-  # Use the optical photon efficiency stacking action for hpDIRC
+  # Use the optical photon efficiency stacking action
+  nm = 1e-6 * mm  # dd4hep length unit
+  drich_lambda_values = [l * nm for l in [ # [nm]
+    315, 325, 340, 350, 370, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 1000
+  ]]
+  drich_efficiency = [e/100. for e in [ # [%]
+    0.00, 4.00, 10.0, 20.0, 30.0, 35.0, 40.0, 38.0, 35.0, 27.0, 20.0, 15.0, 12.0, 8.00, 6.00, 4.00, 0.00
+  ]]
+  pfrich_lambda_values = [l * nm for l in [ # [nm]
+    315, 325, 340, 350, 370, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 1000
+  ]]
+  pfrich_efficiency = [e/100. for e in [ # [%]
+    0.00, 4.00, 10.0, 20.0, 30.0, 35.0, 40.0, 38.0, 35.0, 27.0, 20.0, 15.0, 12.0, 8.00, 6.00, 4.00, 0.00
+  ]]
   RUNNER.action.stack = [
+    {
+      "name": "OpticalPhotonEfficiencyStackingAction",
+      "parameter": {
+        "LambdaValues": drich_lambda_values,
+        "LogicalVolume": "DRICH_(gas|aerogel)",
+        "Efficiency": drich_efficiency,
+      }
+    },
+    {
+      "name": "OpticalPhotonEfficiencyStackingAction",
+      "parameter": {
+        "LambdaValues": pfrich_lambda_values,
+        "LogicalVolume": "PFRICH(_GasVol|-aerogel-[0-9]+-[0-9]+)",
+        "Efficiency": pfrich_efficiency,
+      }
+    },
     {
       "name": "OpticalPhotonEfficiencyStackingAction",
       "parameter": {

src/plugins/include/npdet/OpticalPhotonEfficiencyStackingAction.h

@@ -17,6 +17,12 @@
 #include "DDG4/Geant4StackingAction.h"
 #include "G4OpticalPhoton.hh"
 #include "G4Track.hh"
+#include <algorithm>
+#include <cmath>
+#include <iterator>
+#include <mutex>
+#include <regex>
+#include <stdexcept>
 
 /// Namespace for the AIDA detector description toolkit
 namespace dd4hep {
@@ -31,68 +37,42 @@ namespace dd4hep {
       : Geant4StackingAction(c, n) {
         declareProperty("LambdaMin", m_lambda_min);
         declareProperty("LambdaMax", m_lambda_max);
+        declareProperty("LambdaValues", m_lambda_values);
         declareProperty("Efficiency", m_efficiency);
         declareProperty("LogicalVolume", m_logical_volume);
       };
       /// Default destructor
       virtual ~OpticalPhotonEfficiencyStackingAction() {
-        printout(DEBUG, name(), "Suppressed %d of %d photons in lv %s",
-          m_killed_photons, m_total_photons, m_logical_volume.c_str());
-        printout(DEBUG, name(), "lambda range: [%f,%f] nm",
-          m_lambda_min / CLHEP::nm, m_lambda_max / CLHEP::nm);
-        std::ostringstream oss_efficiency;
-        std::copy(m_efficiency.begin(), m_efficiency.end(),
-          std::ostream_iterator<double>(oss_efficiency, " "));
-        std::string str_efficiency = oss_efficiency.str();
-        printout(DEBUG, name(), "efficiency: %s", str_efficiency.c_str());
+        double pct = m_total_photons > 0 ? 100.0 * m_killed_photons / m_total_photons : 0.0;
+        printout(WARNING, name(), "Suppressed %zu of %zu photons (%.1f%%) in lv %s",
+          m_killed_photons, m_total_photons, pct, m_logical_volume.c_str());
       };
       /// New-stage callback
       virtual void newStage(G4StackManager*) override { };
       /// Preparation callback
       virtual void prepare(G4StackManager*) override { };
       /// Return TrackClassification with enum G4ClassificationOfNewTrack or NoTrackClassification
       virtual TrackClassification classifyNewTrack(G4StackManager*, const G4Track* aTrack) override {
+        std::call_once(m_init_once, [this]() { this->initialize(); });
         // Only apply to optical photons
         if (aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) {
           auto* pv = aTrack->GetVolume();
           if (pv == nullptr) return TrackClassification();
           auto* lv = pv->GetLogicalVolume();
           printout(VERBOSE, name(), "photon in pv %s lv %s",
             pv->GetName().c_str(), lv->GetName().c_str());
-          // Only apply to specified logical volume
-          if (lv->GetName() == m_logical_volume) {
+          // Only apply to specified logical volume regex
+          if (std::regex_match(lv->GetName().c_str(), m_logical_volume_regex)) {
             double mom = aTrack->GetMomentum().mag();
             double lambda = CLHEP::hbarc * CLHEP::twopi / mom;
             printout(VERBOSE, name(), "with mom = %f eV, lambda = %f nm",
               mom / CLHEP::eV, lambda / CLHEP::nm);
 
             m_total_photons++;
-            if (m_lambda_min < lambda && lambda < m_lambda_max) {
-              double efficiency{0.};
-              if (m_efficiency.size() == 0) {
-                // No efficiency specified, assume zero
-                efficiency = 0.;
-                // which means kill
-                ++m_killed_photons;
-                return TrackClassification(fKill);
-
-              } else if (m_efficiency.size() == 1) {
-                // Single constant value over lambda range
-                efficiency = m_efficiency.front();
-
-              } else {
-                // Linear interpolation on lambda grid
-                double lambda_step = (m_lambda_max - m_lambda_min) / (m_efficiency.size() - 1);
-                double div = (lambda - m_lambda_min) / lambda_step;
-                auto i = std::llround(std::floor(div));
-                double t = div - i;
-                double a = m_efficiency[i];
-                double b = m_efficiency[i+1];
-                efficiency = a + t * (b - a);
-                printout(VERBOSE, name(), "a = %f, b = %f, t = %f", a, b, t);
-                printout(VERBOSE, name(), "efficiency %f", efficiency);
-              }
-
+            auto lambda_min = !m_interp_lambda_values.empty() ? m_interp_lambda_values.front() : m_lambda_min;
+            auto lambda_max = !m_interp_lambda_values.empty() ? m_interp_lambda_values.back() : m_lambda_max;
+            if (lambda_min < lambda && lambda < lambda_max) {
+              double efficiency = interpolate(lambda);
               // Edge cases
               if (efficiency == 0.0) {
                 ++m_killed_photons;
@@ -110,7 +90,10 @@ namespace dd4hep {
                 return TrackClassification(fKill);
               }
             } else {
-              printout(VERBOSE, name(), "outside lambda range [%f,%f] nm", m_lambda_min / CLHEP::nm, m_lambda_max / CLHEP::nm);
+              // Outside the QE-specified wavelength range: efficiency is 0
+              printout(VERBOSE, name(), "outside lambda range [%f,%f] nm", lambda_min / CLHEP::nm, lambda_max / CLHEP::nm);
+              ++m_killed_photons;
+              return TrackClassification(fKill);
             }
           } else {
             printout(VERBOSE, name(), "not in volume %s", m_logical_volume.c_str());
@@ -119,9 +102,58 @@ namespace dd4hep {
         return TrackClassification();
       };
     private:
+      void initialize() {
+        try {
+          m_logical_volume_regex = std::regex(m_logical_volume);
+        } catch (const std::regex_error&) {
+          throw std::runtime_error("Invalid LogicalVolume regex");
+        }
+        if (m_efficiency.size() > 1) {
+          if (m_lambda_values.size() == m_efficiency.size()) {
+            m_interp_lambda_values = m_lambda_values;
+          } else {
+            auto lambda_min = m_lambda_min;
+            auto lambda_max = m_lambda_max;
+            m_interp_lambda_values.reserve(m_efficiency.size());
+            auto lambda_step = (lambda_max - lambda_min) / (m_efficiency.size() - 1);
+            for (std::size_t i = 0; i < m_efficiency.size(); ++i) {
+              m_interp_lambda_values.push_back(lambda_min + i * lambda_step);
+            }
+          }
+        } else if (m_lambda_values.size() > 1) {
+          m_interp_lambda_values = m_lambda_values;
+        }
+      }
+      inline double interpolate(double lambda) const {
+        if (m_efficiency.size() == 1) {
+          return m_efficiency.front();
+        }
+        if (m_efficiency.size() < 2 || m_interp_lambda_values.size() < 2) {
+          return 0.0;
+        }
+        auto upper = std::upper_bound(m_interp_lambda_values.begin(), m_interp_lambda_values.end(), lambda);
+        auto i = std::distance(m_interp_lambda_values.begin(), upper) - 1;
+        double a_lambda = m_interp_lambda_values[i];
+        double b_lambda = m_interp_lambda_values[i+1];
+        double t = (lambda - a_lambda) / (b_lambda - a_lambda);
+        double a = m_efficiency[i];
+        double b = m_efficiency[i+1];
+        printout(VERBOSE, name(), "a = %f, b = %f, t = %f", a, b, t);
+#if defined(__cpp_lib_interpolate) && __cpp_lib_interpolate >= 201902L
+        double efficiency = std::lerp(a, b, t);
+#else
+        double efficiency = a + t * (b - a);
+#endif
+        printout(VERBOSE, name(), "efficiency %f", efficiency);
+        return efficiency;
+      }
       double m_lambda_min{0.}, m_lambda_max{0.};
+      std::vector<double> m_lambda_values;
+      std::vector<double> m_interp_lambda_values;
       std::vector<double> m_efficiency;
       std::string m_logical_volume;
+      std::regex m_logical_volume_regex{};
+      std::once_flag m_init_once{};
       std::size_t m_total_photons{0}, m_killed_photons{0};
     };
   }    // End namespace sim