Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/CaTS/src/RadiatorSD.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 // ********************************************************************
  2 // * License and Disclaimer                                           *
  3 // *                                                                  *
  4 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  5 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  6 // * conditions of the Geant4 Software License,  included in the file *
  7 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  8 // * include a list of copyright holders.                             *
  9 // *                                                                  *
 10 // * Neither the authors of this software system, nor their employing *
 11 // * institutes,nor the agencies providing financial support for this *
 12 // * work  make  any representation or  warranty, express or implied, *
 13 // * regarding  this  software system or assume any liability for its *
 14 // * use.  Please see the license in the file  LICENSE  and URL above *
 15 // * for the full disclaimer and the limitation of liability.         *
 16 // *                                                                  *
 17 // * This  code  implementation is the result of  the  scientific and *
 18 // * technical work of the GEANT4 collaboration.                      *
 19 // * By using,  copying,  modifying or  distributing the software (or *
 20 // * any work based  on the software)  you  agree  to acknowledge its *
 21 // * use  in  resulting  scientific  publications,  and indicate your *
 22 // * acceptance of all terms of the Geant4 Software license.          *
 23 // ********************************************************************
 24 //
 25 // ********************************************************************
 26 //
 27 //  CaTS (Calorimetry and Tracking Simulation)
 28 //
 29 //  Authors : Hans Wenzel
 30 //            Soon Yung Jun
 31 //            (Fermi National Accelerator Laboratory)
 32 //
 33 // History
 34 //   October 18th, 2021 : first implementation
 35 //
 36 // ********************************************************************
 37 //
 38 /// \file RadiatorSD.cc
 39 /// \brief Implementation of the CaTS::RadiatorSD class
 40 
 41 // Geant4 headers
 42 #include "G4Step.hh"
 43 #include "G4Track.hh"
 44 #ifdef WITH_G4OPTICKS
 45 #  include "G4ThreeVector.hh"
 46 #  include "G4ios.hh"
 47 #  include "G4UnitsTable.hh"
 48 #  include "G4SystemOfUnits.hh"
 49 #  include "G4VProcess.hh"
 50 #  include "G4VRestDiscreteProcess.hh"
 51 #  include "G4SDManager.hh"
 52 #  include "G4HCofThisEvent.hh"
 53 #  include "G4Opticks.hh"
 54 #  include "TrackInfo.hh"
 55 #  include "OpticksGenstep.h"
 56 #  include "OpticksFlags.hh"
 57 #  include "G4OpticksHit.hh"
 58 #  include "G4EventManager.hh"
 59 #  include "G4Event.hh"
 60 #  include "G4RunManager.hh"
 61 #  include "G4Version.hh"
 62 #  include "PhotonSD.hh"
 63 #  include "G4Cerenkov.hh"
 64 #  include "G4Scintillation.hh"
 65 #  include <string>
 66 #endif
 67 // project headers
 68 #include "RadiatorSD.hh"
 69 #include "ConfigurationManager.hh"
 70 
 71 RadiatorSD::RadiatorSD(G4String name)
 72   : G4VSensitiveDetector(name)
 73 {
 74   verbose = ConfigurationManager::getInstance()->isEnable_verbose();
 75 }
 76 
 77 void RadiatorSD::Initialize(G4HCofThisEvent*) {}
 78 
 79 G4bool RadiatorSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
 80 {
 81   G4double edep = aStep->GetTotalEnergyDeposit();
 82   if(edep == 0.)
 83     return false;
 84   // only deal with charged particles
 85   G4Track* aTrack = aStep->GetTrack();
 86   G4double charge = aTrack->GetDynamicParticle()->GetCharge();
 87   if(charge == 0)
 88     return false;
 89 #ifdef WITH_G4OPTICKS
 90   if(ConfigurationManager::getInstance()->isEnable_opticks())
 91   {
 92     G4int materialIndex = 0;
 93     if(first)
 94     {
 95       aMaterial     = aTrack->GetMaterial();
 96       materialIndex = aMaterial->GetIndex();
 97       if(verbose)
 98       {
 99         G4cout << "*******************************" << G4endl;
100         G4cout << "RadiatorSD::ProcessHits initializing Material:  "
101                << aMaterial->GetName() << " " << G4endl;
102         G4cout << "RadiatorSD::ProcessHits: Name "
103                << aStep->GetPreStepPoint()
104                     ->GetPhysicalVolume()
105                     ->GetLogicalVolume()
106                     ->GetName()
107                << G4endl;
108       }
109       aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
110       if(verbose)
111       {
112         aMaterialPropertiesTable->DumpTable();
113       }
114       //
115       // properties related to Scintillation
116       //
117 #  if(G4VERSION_NUMBER > 1072)
118       YieldRatio =
119         aMaterialPropertiesTable->GetConstProperty(kSCINTILLATIONYIELD1) /
120         aMaterialPropertiesTable->GetConstProperty(
121           kSCINTILLATIONYIELD2);  // slowerRatio,
122       FastTimeConstant = aMaterialPropertiesTable->GetConstProperty(
123         kSCINTILLATIONTIMECONSTANT1);  // TimeConstant,
124       SlowTimeConstant = aMaterialPropertiesTable->GetConstProperty(
125         kSCINTILLATIONTIMECONSTANT2);  // slowerTimeConstant,
126 #  else
127       Fast_Intensity = aMaterialPropertiesTable->GetProperty(kFASTCOMPONENT);
128       Slow_Intensity = aMaterialPropertiesTable->GetProperty(kSLOWCOMPONENT);
129       YieldRatio     = aMaterialPropertiesTable->GetConstProperty(kYIELDRATIO);
130 #  endif
131       ScintillationType = Slow;
132       //
133       // properties related to Cerenkov
134       //
135       Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
136 #  if(G4VERSION_NUMBER > 1072)
137       Pmin = Rindex->GetMinEnergy();
138       Pmax = Rindex->GetMaxEnergy();
139 #  else
140       Pmin           = Rindex->GetMinLowEdgeEnergy();
141       Pmax           = Rindex->GetMaxLowEdgeEnergy();
142 #  endif
143       dp   = Pmax - Pmin;
144       nMax = Rindex->GetMaxValue();
145       if(verbose)
146       {
147         G4cout << "nMax: " << nMax << "  Pmin: " << Pmin << "  Pmax: " << Pmax
148                << "  dp: " << dp << G4endl;
149         Rindex->DumpValues();
150       }
151       //
152       first = false;
153     }                    // end if first
154     G4int Sphotons = 0;  // number of scintillation photons this step
155     G4int Cphotons = 0;  // number of Cerenkov photons this step
156     //
157     // info needed for generating Cerenkov photons on the GPU;
158     //
159     G4double maxCos                      = 0.0;
160     G4double maxSin2                     = 0.0;
161     G4double beta                        = 0.0;
162     G4double beta1                       = 0.0;
163     G4double beta2                       = 0.0;
164     G4double BetaInverse                 = 0.0;
165     G4double MeanNumberOfPhotons1        = 0.0;
166     G4double MeanNumberOfPhotons2        = 0.0;
167     G4SteppingManager* fpSteppingManager = G4EventManager::GetEventManager()
168                                              ->GetTrackingManager()
169                                              ->GetSteppingManager();
170     G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
171     if(stepStatus != fAtRestDoItProc)
172     {
173       G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector();
174       size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops();
175       for(size_t i3 = 0; i3 < MAXofPostStepLoops; i3++)
176       {
177         if((*procPost)[i3]->GetProcessName() == "Cerenkov")
178         {
179           G4Cerenkov* proc = (G4Cerenkov*) (*procPost)[i3];
180           thePhysicsTable  = proc->GetPhysicsTable();
181           CerenkovAngleIntegrals =
182             (G4PhysicsOrderedFreeVector*) ((*thePhysicsTable)(materialIndex));
183           Cphotons = proc->GetNumPhotons();
184           if(Cphotons > 0)
185           {
186             beta1       = aStep->GetPreStepPoint()->GetBeta();
187             beta2       = aStep->GetPostStepPoint()->GetBeta();
188             beta        = (beta1 + beta2) * 0.5;
189             BetaInverse = 1. / beta;
190             maxCos      = BetaInverse / nMax;
191             maxSin2     = (1.0 - maxCos) * (1.0 + maxCos);
192             MeanNumberOfPhotons1 =
193               proc->GetAverageNumberOfPhotons(charge, beta1, aMaterial, Rindex);
194             MeanNumberOfPhotons2 =
195               proc->GetAverageNumberOfPhotons(charge, beta2, aMaterial, Rindex);
196           }
197         }
198         if((*procPost)[i3]->GetProcessName() == "Scintillation")
199         {
200           G4Scintillation* proc1 = (G4Scintillation*) (*procPost)[i3];
201           Sphotons               = proc1->GetNumPhotons();
202         }
203       }
204     }
205     tSphotons += Sphotons;
206     tCphotons += Cphotons;
207     G4ThreeVector deltaPosition = aStep->GetDeltaPosition();
208     G4double ScintillationTime  = 0. * ns;
209     G4int scntId                = 1;
210     G4StepPoint* pPreStepPoint  = aStep->GetPreStepPoint();
211     G4ThreeVector x0            = pPreStepPoint->GetPosition();
212     G4ThreeVector p0            = aStep->GetDeltaPosition().unit();
213     //
214     // harvest the Scintillation photon gensteps:
215     //
216     if(Sphotons > 0)
217     {
218       G4double ScintillationRiseTime = 0.0;
219       G4Opticks::Get()->collectGenstep_G4Scintillation_1042(
220         aTrack, aStep, Sphotons, scntId, ScintillationTime,
221         ScintillationRiseTime);
222     }
223     //
224     // harvest the Cerenkov photon gensteps:
225     //
226     if(Cphotons > 0)
227     {
228       G4Opticks::Get()->collectGenstep_G4Cerenkov_1042(
229         aTrack, aStep, Cphotons, BetaInverse, Pmin, Pmax, maxCos, maxSin2,
230         MeanNumberOfPhotons1, MeanNumberOfPhotons2);
231     }
232     G4Opticks* g4ok      = G4Opticks::Get();
233     G4RunManager* rm     = G4RunManager::GetRunManager();
234     const G4Event* event = rm->GetCurrentEvent();
235     G4int eventid        = event->GetEventID();
236     G4OpticksHit hit;
237     unsigned num_photons = g4ok->getNumPhotons();
238     if(num_photons > ConfigurationManager::getInstance()->getMaxPhotons())
239     {
240       g4ok->propagateOpticalPhotons(eventid);
241       G4HCtable* hctable = G4SDManager::GetSDMpointer()->GetHCtable();
242       for(G4int i = 0; i < hctable->entries(); ++i)
243       {
244         std::string sdn   = hctable->GetSDname(i);
245         std::size_t found = sdn.find("Photondetector");
246         if(found != std::string::npos)
247         {
248           PhotonSD* aSD =
249             (PhotonSD*) G4SDManager::GetSDMpointer()->FindSensitiveDetector(
250               sdn);
251           aSD->AddOpticksHits();
252         }
253       }
254       g4ok->reset();
255     }
256   }
257 #endif
258   return true;
259 }
260 
261 void RadiatorSD::EndOfEvent(G4HCofThisEvent*)
262 {
263   tSphotons = 0;
264   tCphotons = 0;
265 }
266