diff --git a/src/dd4pod/python/npsim.py b/src/dd4pod/python/npsim.py index 2c35e91..c3e8d77 100755 --- a/src/dd4pod/python/npsim.py +++ b/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": { diff --git a/src/plugins/include/npdet/OpticalPhotonEfficiencyStackingAction.h b/src/plugins/include/npdet/OpticalPhotonEfficiencyStackingAction.h index e902f4a..4a57f6d 100644 --- a/src/plugins/include/npdet/OpticalPhotonEfficiencyStackingAction.h +++ b/src/plugins/include/npdet/OpticalPhotonEfficiencyStackingAction.h @@ -17,6 +17,12 @@ #include "DDG4/Geant4StackingAction.h" #include "G4OpticalPhoton.hh" #include "G4Track.hh" +#include +#include +#include +#include +#include +#include /// Namespace for the AIDA detector description toolkit namespace dd4hep { @@ -31,20 +37,15 @@ 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(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 { }; @@ -52,6 +53,7 @@ namespace dd4hep { 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(); @@ -59,40 +61,18 @@ namespace dd4hep { 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 m_lambda_values; + std::vector m_interp_lambda_values; std::vector 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