| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // *******************************************
27 //
28 // CaTS (Calorimetry and Tracking Simulation)
29 //
30 // Authors : Hans Wenzel
31 // Soon Yung Jun
32 // (Fermi National Accelerator Labo
33 //
34 // History
35 // October 18th, 2021 : first implementation
36 //
37 // *******************************************
38 //
39 /// \file lArTPCSD.cc
40 /// \brief Implementation of the CaTS::lArTPCS
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(
80 if(verbose)
81 {
82 G4cout << collectionName.size() << " lAr
83 << " collection Name: " << HCname <
84 }
85 fHCID = -1;
86 first = true;
87 }
88
89 void lArTPCSD::Initialize(G4HCofThisEvent* hce
90 {
91 flArTPCHitsCollection =
92 new lArTPCHitsCollection(SensitiveDetector
93 if(fHCID < 0)
94 {
95 if(verbose)
96 {
97 G4cout << "lArTPCSD::Initialize: " << S
98 << collectionName[0] << G4endl;
99 }
100 fHCID = G4SDManager::GetSDMpointer()->GetC
101 }
102 hce->AddHitsCollection(fHCID, flArTPCHitsCol
103 }
104
105 G4bool lArTPCSD::ProcessHits(G4Step* aStep, G4
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
113 if(charge == 0)
114 return false;
115 G4double ds = aStep->GetStepLength();
116 lArTPCHit* newHit = new lArTPCHit(
117 NumElectrons(edep, ds), aStep->GetPostStep
118 aStep->GetPostStepPoint()->GetPosition().g
119 aStep->GetPostStepPoint()->GetPosition().g
120 flArTPCHitsCollection->insert(newHit);
121 #ifdef WITH_G4OPTICKS
122 if(ConfigurationManager::getInstance()->isEn
123 {
124 if(first)
125 {
126 aMaterial = aTrack->GetMaterial();
127 materialIndex = aMaterial->GetIndex();
128 if(verbose)
129 {
130 G4cout << "***************************
131 G4cout << "RadiatorSD::ProcessHits ini
132 << aMaterial->GetName() << " "
133 G4cout << "RadiatorSD::ProcessHits: Na
134 << aStep->GetPreStepPoint()
135 ->GetPhysicalVolume()
136 ->GetLogicalVolume()
137 ->GetName()
138 << G4endl;
139 }
140 aMaterialPropertiesTable = aMaterial->Ge
141 if(verbose)
142 {
143 aMaterialPropertiesTable->DumpTable();
144 }
145 //
146 // properties related to Scintillation
147 //
148 # if(G4VERSION_NUMBER > 1072)
149 YieldRatio =
150 aMaterialPropertiesTable->GetConstProp
151 aMaterialPropertiesTable->GetConstProp
152 kSCINTILLATIONYIELD2); // slowerRat
153 FastTimeConstant = aMaterialPropertiesTa
154 kSCINTILLATIONTIMECONSTANT1); // Time
155 SlowTimeConstant = aMaterialPropertiesTa
156 kSCINTILLATIONTIMECONSTANT2); // slow
157 # else
158 Fast_Intensity = aMaterialPropertiesTabl
159 Slow_Intensity = aMaterialPropertiesTabl
160 YieldRatio = aMaterialPropertiesTabl
161 # endif
162 ScintillationType = Slow;
163 //
164 // properties related to Cerenkov
165 //
166 Rindex = aMaterialPropertiesTable->GetPr
167 # if(G4VERSION_NUMBER > 1072)
168 Pmin = Rindex->GetMinEnergy();
169 Pmax = Rindex->GetMaxEnergy();
170 # else
171 Pmin = Rindex->GetMinLowEdgeEn
172 Pmax = Rindex->GetMaxLowEdgeEn
173 # endif
174 dp = Pmax - Pmin;
175 if(verbose)
176 {
177 G4cout << "nMax: " << nMax << "Pmin: "
178 << "dp: " << dp << G4endl;
179 Rindex->DumpValues();
180 }
181 //
182 first = false;
183 }
184 G4int Sphotons = 0; // number of scintill
185 G4int Cphotons = 0; // number of Cerenkov
186 //
187 // info needed for generating Cerenkov pho
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 = G4E
198 -
199 -
200 G4StepStatus stepStatus = fpSteppingManage
201 if(stepStatus != fAtRestDoItProc)
202 {
203 G4ProcessVector* procPost = fpSteppingMa
204 size_t MAXofPostStepLoops = fpSteppingMa
205 for(size_t i3 = 0; i3 < MAXofPostStepLoo
206 {
207 if((*procPost)[i3]->GetProcessName() =
208 {
209 G4Cerenkov* proc = (G4Cerenkov*) (*p
210 thePhysicsTable = proc->GetPhysicsT
211 CerenkovAngleIntegrals =
212 (G4PhysicsOrderedFreeVector*) ((*t
213 Cphotons = proc->GetNumPhotons();
214 if(Cphotons > 0)
215 {
216 beta1 = aStep->GetPreStepPoi
217 beta2 = aStep->GetPostStepPo
218 beta = (beta1 + beta2) * 0.
219 BetaInverse = 1. / beta;
220 maxCos = BetaInverse / nMax;
221 maxSin2 = (1.0 - maxCos) * (1.
222 MeanNumberOfPhotons1 =
223 proc->GetAverageNumberOfPhotons(
224 MeanNumberOfPhotons2 =
225 proc->GetAverageNumberOfPhotons(
226 }
227 }
228 if((*procPost)[i3]->GetProcessName() =
229 {
230 G4Scintillation* proc1 = (G4Scintill
231 Sphotons = proc1->GetN
232 }
233 }
234 }
235 tSphotons += Sphotons;
236 tCphotons += Cphotons;
237 G4ThreeVector deltaPosition = aStep->GetDe
238 G4double ScintillationTime = 0. * ns;
239 G4int scntId = 1;
240 G4StepPoint* pPreStepPoint = aStep->GetPr
241 G4ThreeVector x0 = pPreStepPoin
242 G4ThreeVector p0 = aStep->GetDe
243 //
244 // harvest the Scintillation photon genste
245 //
246 if(Sphotons > 0)
247 {
248 G4double ScintillationRiseTime = 0.0;
249 G4Opticks::Get()->collectGenstep_G4Scint
250 aTrack, aStep, Sphotons, scntId, Scint
251 ScintillationRiseTime);
252 }
253 //
254 // harvest the Cerenkov photon gensteps:
255 //
256 if(Cphotons > 0)
257 {
258 G4Opticks::Get()->collectGenstep_G4Ceren
259 aTrack, aStep, Cphotons, BetaInverse,
260 MeanNumberOfPhotons1, MeanNumberOfPhot
261 }
262 G4Opticks* g4ok = G4Opticks::Get();
263 G4RunManager* rm = G4RunManager::GetRu
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::get
269 {
270 g4ok->propagateOpticalPhotons(eventid);
271 G4HCtable* hctable = G4SDManager::GetSDM
272 for(G4int i = 0; i < hctable->entries();
273 {
274 std::string sdn = hctable->GetSDname
275 std::size_t found = sdn.find("Photonde
276 if(found != std::string::npos)
277 {
278 PhotonSD* aSD =
279 (PhotonSD*) G4SDManager::GetSDMpoi
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->entrie
296 if(verbose)
297 {
298 G4cout << " Number of lArTPCHits: " << Nb
299 }
300 }
301
302 G4double lArTPCSD::NumElectrons(G4double edep,
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
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
323 return fNumIonElectrons;
324 }
325