eic-opticks/dd4hepplugins/OpticsSteppingAction.cc

0001 #include "OpticsSteppingAction.hh"
0002
0003 #include <DD4hep/InstanceCount.h>
0004 #include <DDG4/Factories.h>
0005
0006 #include <G4AutoLock.hh>
0007 #include <G4Cerenkov.hh>
0008 #include <G4Material.hh>
0009 #include <G4MaterialPropertiesTable.hh>
0010 #include <G4OpticalPhoton.hh>
0011 #include <G4ProcessManager.hh>
0012 #include <G4Scintillation.hh>
0013 #include <G4Step.hh>
0014 #include <G4SteppingManager.hh>
0015 #include <G4Track.hh>
0016
0017 #include <U4.hh>
0018
0019 namespace
0020 {
0021 G4Mutex genstep_mutex = G4MUTEX_INITIALIZER;
0022 }
0023
0024 namespace ddeicopticks
0025 {
0026 //---------------------------------------------------------------------------//
0027 OpticsSteppingAction::OpticsSteppingAction(dd4hep::sim::Geant4Context *ctxt, std::string const &name)
0028     : dd4hep::sim::Geant4SteppingAction(ctxt, name)
0029 {
0030     dd4hep::InstanceCount::increment(this);
0031     declareProperty("Verbose", verbose_);
0032 }
0033
0034 //---------------------------------------------------------------------------//
0035 OpticsSteppingAction::~OpticsSteppingAction()
0036 {
0037     dd4hep::InstanceCount::decrement(this);
0038 }
0039
0040 //---------------------------------------------------------------------------//
0041 void OpticsSteppingAction::operator()(const G4Step *step, G4SteppingManager *mgr)
0042 {
0043     G4SteppingManager *fpSteppingManager = mgr;
0044
0045     const G4Track *track = step->GetTrack();
0046
0047     G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
0048     if (stepStatus == fAtRestDoItProc)
0049         return;
0050
0051     // Skip optical photons — they don't produce gensteps
0052     if (track->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
0053         return;
0054
0055     G4ProcessVector *procPost = fpSteppingManager->GetfPostStepDoItVector();
0056     size_t nproc = fpSteppingManager->GetMAXofPostStepLoops();
0057
0058     if (verbose_ > 1)
0059     {
0060         std::string procs;
0061         for (size_t i = 0; i < nproc; i++)
0062             if ((*procPost)[i])
0063                 procs += (*procPost)[i]->GetProcessName() + " ";
0064         info("  nproc=%zu processes=[%s]", nproc, procs.c_str());
0065     }
0066
0067     for (size_t i = 0; i < nproc; i++)
0068     {
0069         if (!(*procPost)[i])
0070             continue;
0071
0072         const G4String &procName = (*procPost)[i]->GetProcessName();
0073
0074         if (procName == "Cerenkov")
0075         {
0076             const G4DynamicParticle *aParticle = track->GetDynamicParticle();
0077             G4double charge = aParticle->GetDefinition()->GetPDGCharge();
0078             const G4Material *mat = track->GetMaterial();
0079             G4MaterialPropertiesTable *mpt = mat->GetMaterialPropertiesTable();
0080             if (!mpt)
0081                 continue;
0082
0083             G4MaterialPropertyVector *Rindex = mpt->GetProperty(kRINDEX);
0084             if (!Rindex || Rindex->GetVectorLength() == 0)
0085                 continue;
0086
0087             G4Cerenkov *cer = static_cast<G4Cerenkov *>((*procPost)[i]);
0088             G4int numPhotons = cer->GetNumPhotons();
0089             if (numPhotons <= 0)
0090                 continue;
0091
0092             G4double Pmin = Rindex->Energy(0);
0093             G4double Pmax = Rindex->GetMaxEnergy();
0094             G4double nMax = Rindex->GetMaxValue();
0095             G4double beta1 = step->GetPreStepPoint()->GetBeta();
0096             G4double beta2 = step->GetPostStepPoint()->GetBeta();
0097             G4double beta = (beta1 + beta2) * 0.5;
0098             G4double BetaInverse = 1.0 / beta;
0099             G4double maxCos = BetaInverse / nMax;
0100             G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
0101             G4double MeanNumberOfPhotons1 = cer->GetAverageNumberOfPhotons(charge, beta1, mat, Rindex);
0102             G4double MeanNumberOfPhotons2 = cer->GetAverageNumberOfPhotons(charge, beta2, mat, Rindex);
0103
0104             {
0105                 G4AutoLock lock(&genstep_mutex);
0106                 U4::CollectGenstep_G4Cerenkov_modified(track, step, numPhotons, BetaInverse, Pmin, Pmax, maxCos,
0107                                                        maxSin2, MeanNumberOfPhotons1, MeanNumberOfPhotons2);
0108             }
0109
0110             if (verbose_ > 0)
0111                 info("Cerenkov genstep: %d photons", numPhotons);
0112         }
0113         else if (procName == "Scintillation")
0114         {
0115             G4Scintillation *scint = static_cast<G4Scintillation *>((*procPost)[i]);
0116             G4int numPhotons = scint->GetNumPhotons();
0117             if (numPhotons <= 0)
0118                 continue;
0119
0120             const G4Material *mat = track->GetMaterial();
0121             G4MaterialPropertiesTable *mpt = mat->GetMaterialPropertiesTable();
0122             if (!mpt)
0123                 continue;
0124
0125             G4double scintTime = 0.0;
0126             if (mpt->ConstPropertyExists(kSCINTILLATIONTIMECONSTANT1))
0127                 scintTime = mpt->GetConstProperty(kSCINTILLATIONTIMECONSTANT1);
0128
0129             {
0130                 G4AutoLock lock(&genstep_mutex);
0131                 U4::CollectGenstep_DsG4Scintillation_r4695(track, step, numPhotons,
0132                                                            /*scnt=*/1, scintTime);
0133             }
0134
0135             if (verbose_ > 0)
0136                 info("Scintillation genstep: %d photons", numPhotons);
0137         }
0138     }
0139 }
0140
0141 //---------------------------------------------------------------------------//
0142 } // namespace ddeicopticks
0143
0144 DECLARE_GEANT4ACTION_NS(ddeicopticks, OpticsSteppingAction)