Geant4 Cross Reference

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