include/npdet/OpticalPhotonEfficiencyStackingAction.h

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013 #ifndef OPTICALPHOTONEFFICIENCYSTACKINGACTION_H
0014 #define OPTICALPHOTONEFFICIENCYSTACKINGACTION_H
0015
0016 #include "DDG4/Geant4Random.h"
0017 #include "DDG4/Geant4StackingAction.h"
0018 #include "G4OpticalPhoton.hh"
0019 #include "G4Track.hh"
0020
0021 /// Namespace for the AIDA detector description toolkit
0022 namespace dd4hep {
0023
0024   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0025   namespace sim {
0026
0027     class OpticalPhotonEfficiencyStackingAction: public Geant4StackingAction {
0028     public:
0029       /// Standard constructor with initializing arguments
0030       OpticalPhotonEfficiencyStackingAction(Geant4Context* c, const std::string& n)
0031       : Geant4StackingAction(c, n) {
0032         declareProperty("LambdaMin", m_lambda_min);
0033         declareProperty("LambdaMax", m_lambda_max);
0034         declareProperty("Efficiency", m_efficiency);
0035         declareProperty("LogicalVolume", m_logical_volume);
0036       };
0037       /// Default destructor
0038       virtual ~OpticalPhotonEfficiencyStackingAction() {
0039         printout(DEBUG, name(), "Suppressed %d of %d photons in lv %s",
0040           m_killed_photons, m_total_photons, m_logical_volume.c_str());
0041         printout(DEBUG, name(), "lambda range: [%f,%f] nm",
0042           m_lambda_min / CLHEP::nm, m_lambda_max / CLHEP::nm);
0043         std::ostringstream oss_efficiency;
0044         std::copy(m_efficiency.begin(), m_efficiency.end(),
0045           std::ostream_iterator<double>(oss_efficiency, " "));
0046         std::string str_efficiency = oss_efficiency.str();
0047         printout(DEBUG, name(), "efficiency: %s", str_efficiency.c_str());
0048       };
0049       /// New-stage callback
0050       virtual void newStage(G4StackManager*) override { };
0051       /// Preparation callback
0052       virtual void prepare(G4StackManager*) override { };
0053       /// Return TrackClassification with enum G4ClassificationOfNewTrack or NoTrackClassification
0054       virtual TrackClassification classifyNewTrack(G4StackManager*, const G4Track* aTrack) override {
0055         // Only apply to optical photons
0056         if (aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) {
0057           auto* pv = aTrack->GetVolume();
0058           auto* lv = pv->GetLogicalVolume();
0059           printout(VERBOSE, name(), "photon in pv %s lv %s",
0060             pv->GetName().c_str(), lv->GetName().c_str());
0061           // Only apply to specified logical volume
0062           if (lv->GetName() == m_logical_volume) {
0063             double mom = aTrack->GetMomentum().mag();
0064             double lambda = CLHEP::hbarc * CLHEP::twopi / mom;
0065             printout(VERBOSE, name(), "with mom = %f eV, lambda = %f nm",
0066               mom / CLHEP::eV, lambda / CLHEP::nm);
0067
0068             m_total_photons++;
0069             if (m_lambda_min < lambda && lambda < m_lambda_max) {
0070               double efficiency{0.};
0071               if (m_efficiency.size() == 0) {
0072                 // No efficiency specified, assume zero
0073                 efficiency = 0.;
0074                 // which means kill
0075                 ++m_killed_photons;
0076                 return TrackClassification(fKill);
0077
0078               } else if (m_efficiency.size() == 1) {
0079                 // Single constant value over lambda range
0080                 efficiency = m_efficiency.front();
0081
0082               } else {
0083                 // Linear interpolation on lambda grid
0084                 double lambda_step = (m_lambda_max - m_lambda_min) / (m_efficiency.size() - 1);
0085                 double div = (lambda - m_lambda_min) / lambda_step;
0086                 auto i = std::llround(std::floor(div));
0087                 double t = div - i;
0088                 double a = m_efficiency[i];
0089                 double b = m_efficiency[i+1];
0090                 efficiency = a + t * (b - a);
0091                 printout(VERBOSE, name(), "a = %f, b = %f, t = %f", a, b, t);
0092                 printout(VERBOSE, name(), "efficiency %f", efficiency);
0093               }
0094
0095               // Edge cases
0096               if (efficiency == 0.0) {
0097                 ++m_killed_photons;
0098                 return TrackClassification(fKill);
0099               }
0100               if (efficiency == 1.0) return TrackClassification();
0101
0102               // Throw random value
0103               Geant4Event&  evt = context()->event();
0104               Geant4Random& rnd = evt.random();
0105               double random = rnd.uniform();
0106               if (random > efficiency) {
0107                 printout(VERBOSE, name(), "photon killed");
0108                 ++m_killed_photons;
0109                 return TrackClassification(fKill);
0110               }
0111             } else {
0112               printout(VERBOSE, name(), "outside lambda range [%f,%f] nm", m_lambda_min / CLHEP::nm, m_lambda_max / CLHEP::nm);
0113             }
0114           } else {
0115             printout(VERBOSE, name(), "not in volume %s", m_logical_volume.c_str());
0116           }
0117         }
0118         return TrackClassification();
0119       };
0120     private:
0121       double m_lambda_min{0.}, m_lambda_max{0.};
0122       std::vector<double> m_efficiency;
0123       std::string m_logical_volume;
0124       std::size_t m_total_photons{0}, m_killed_photons{0};
0125     };
0126   }    // End namespace sim
0127 }      // End namespace dd4hep
0128
0129 #endif // OPTICALPHOTONEFFICIENCYSTACKINGACTION_H