| 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 /// \file optical/OpNovice2/src/SteppingAction
28 /// \brief Implementation of the SteppingActio
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31
32 #include "SteppingAction.hh"
33
34 #include "HistoManager.hh"
35 #include "Run.hh"
36 #include "SteppingMessenger.hh"
37 #include "TrackInformation.hh"
38
39 #include "G4Cerenkov.hh"
40 #include "G4Event.hh"
41 #include "G4EventManager.hh"
42 #include "G4OpBoundaryProcess.hh"
43 #include "G4OpticalPhoton.hh"
44 #include "G4ProcessManager.hh"
45 #include "G4RunManager.hh"
46 #include "G4Scintillation.hh"
47 #include "G4Step.hh"
48 #include "G4SteppingManager.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4Track.hh"
51
52 //....oooOO0OOooo........oooOO0OOooo........oo
53 SteppingAction::SteppingAction() : G4UserStepp
54 {
55 fSteppingMessenger = new SteppingMessenger(t
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oo
59 SteppingAction::~SteppingAction()
60 {
61 delete fSteppingMessenger;
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oo
65 void SteppingAction::UserSteppingAction(const
66 {
67 static G4ParticleDefinition* opticalphoton =
68
69 G4AnalysisManager* analysisMan = G4AnalysisM
70 Run* run = static_cast<Run*>(G4RunManager::G
71
72 G4Track* track = step->GetTrack();
73 G4StepPoint* endPoint = step->GetPostStepPoi
74 G4StepPoint* startPoint = step->GetPreStepPo
75
76 const G4DynamicParticle* theParticle = track
77 const G4ParticleDefinition* particleDef = th
78
79 auto trackInfo = (TrackInformation*)(track->
80
81 if (particleDef == opticalphoton) {
82 const G4VProcess* pds = endPoint->GetProce
83 G4String procname = pds->GetProcessName();
84 if (procname == "OpAbsorption") {
85 run->AddOpAbsorption();
86 if (trackInfo->GetIsFirstTankX()) {
87 run->AddOpAbsorptionPrior();
88 }
89 }
90 else if (procname == "OpRayleigh") {
91 run->AddRayleigh();
92 }
93 else if (procname == "OpWLS") {
94 G4double en = track->GetKineticEnergy();
95 run->AddWLSAbsorption();
96 run->AddWLSAbsorptionEnergy(en);
97 analysisMan->FillH1(4, en / eV); // abs
98 // loop over secondaries, create statist
99 // const std::vector<const G4Track*>* se
100 auto secondaries = step->GetSecondaryInC
101 for (auto sec : *secondaries) {
102 en = sec->GetKineticEnergy();
103 run->AddWLSEmission();
104 run->AddWLSEmissionEnergy(en);
105 analysisMan->FillH1(5, en / eV); // e
106 G4double time = sec->GetGlobalTime();
107 analysisMan->FillH1(6, time / ns);
108 }
109 }
110 else if (procname == "OpWLS2") {
111 G4double en = track->GetKineticEnergy();
112 run->AddWLS2Absorption();
113 run->AddWLS2AbsorptionEnergy(en);
114 analysisMan->FillH1(7, en / eV); // abs
115 // loop over secondaries, create statist
116 // const std::vector<const G4Track*>* se
117 auto secondaries = step->GetSecondaryInC
118 for (auto sec : *secondaries) {
119 en = sec->GetKineticEnergy();
120 run->AddWLS2Emission();
121 run->AddWLS2EmissionEnergy(en);
122 analysisMan->FillH1(8, en / eV); // e
123 G4double time = sec->GetGlobalTime();
124 analysisMan->FillH1(9, time / ns);
125 }
126 }
127
128 // optical process has endpt on bdry,
129 if (endPoint->GetStepStatus() == fGeomBoun
130 G4ThreeVector p0 = startPoint->GetMoment
131 G4ThreeVector p1 = endPoint->GetMomentum
132
133 G4OpBoundaryProcessStatus theStatus = Un
134
135 G4ProcessManager* OpManager = opticalpho
136 G4ProcessVector* postStepDoItVector = Op
137 G4int n_proc = postStepDoItVector->entri
138
139 if (trackInfo->GetIsFirstTankX()) {
140 G4double px1 = p1.x();
141 G4double py1 = p1.y();
142 G4double pz1 = p1.z();
143 // do not count Absorbed or Detected p
144 if (track->GetTrackStatus() != fStopAn
145 if (px1 < 0.) {
146 analysisMan->FillH1(11, px1);
147 analysisMan->FillH1(12, py1);
148 analysisMan->FillH1(13, pz1);
149 }
150 else {
151 analysisMan->FillH1(14, px1);
152 analysisMan->FillH1(15, py1);
153 analysisMan->FillH1(16, pz1);
154 }
155 }
156
157 trackInfo->SetIsFirstTankX(false);
158 run->AddTotalSurface();
159
160 for (G4int i = 0; i < n_proc; ++i) {
161 G4VProcess* currentProcess = (*postS
162
163 auto opProc = dynamic_cast<G4OpBound
164 if (opProc) {
165 G4double angle = std::acos(p0.x())
166 theStatus = opProc->GetStatus();
167 analysisMan->FillH1(10, theStatus)
168 switch (theStatus) {
169 case Transmission:
170 run->AddTransmission();
171 analysisMan->FillH1(25, angle
172 break;
173 case FresnelRefraction:
174 run->AddFresnelRefraction();
175 analysisMan->FillH1(17, px1);
176 analysisMan->FillH1(18, py1);
177 analysisMan->FillH1(19, pz1);
178 analysisMan->FillH1(20, angle
179 break;
180 case FresnelReflection:
181 run->AddFresnelReflection();
182 analysisMan->FillH1(21, angle
183 analysisMan->FillH1(23, angle
184 break;
185 case TotalInternalReflection:
186 run->AddTotalInternalReflectio
187 analysisMan->FillH1(22, angle
188 analysisMan->FillH1(23, angle
189 break;
190 case LambertianReflection:
191 run->AddLambertianReflection()
192 break;
193 case LobeReflection:
194 run->AddLobeReflection();
195 break;
196 case SpikeReflection:
197 run->AddSpikeReflection();
198 analysisMan->FillH1(26, angle
199 break;
200 case BackScattering:
201 run->AddBackScattering();
202 break;
203 case Absorption:
204 analysisMan->FillH1(24, angle
205 run->AddAbsorption();
206 break;
207 case Detection:
208 run->AddDetection();
209 break;
210 case NotAtBoundary:
211 run->AddNotAtBoundary();
212 break;
213 case SameMaterial:
214 run->AddSameMaterial();
215 break;
216 case StepTooSmall:
217 run->AddStepTooSmall();
218 break;
219 case NoRINDEX:
220 run->AddNoRINDEX();
221 break;
222 case PolishedLumirrorAirReflecti
223 run->AddPolishedLumirrorAirRef
224 break;
225 case PolishedLumirrorGlueReflect
226 run->AddPolishedLumirrorGlueRe
227 break;
228 case PolishedAirReflection:
229 run->AddPolishedAirReflection(
230 break;
231 case PolishedTeflonAirReflection
232 run->AddPolishedTeflonAirRefle
233 break;
234 case PolishedTiOAirReflection:
235 run->AddPolishedTiOAirReflecti
236 break;
237 case PolishedTyvekAirReflection:
238 run->AddPolishedTyvekAirReflec
239 break;
240 case PolishedVM2000AirReflection
241 run->AddPolishedVM2000AirRefle
242 break;
243 case PolishedVM2000GlueReflectio
244 run->AddPolishedVM2000AirRefle
245 break;
246 case EtchedLumirrorAirReflection
247 run->AddEtchedLumirrorAirRefle
248 break;
249 case EtchedLumirrorGlueReflectio
250 run->AddEtchedLumirrorGlueRefl
251 break;
252 case EtchedAirReflection:
253 run->AddEtchedAirReflection();
254 break;
255 case EtchedTeflonAirReflection:
256 run->AddEtchedTeflonAirReflect
257 break;
258 case EtchedTiOAirReflection:
259 run->AddEtchedTiOAirReflection
260 break;
261 case EtchedTyvekAirReflection:
262 run->AddEtchedTyvekAirReflecti
263 break;
264 case EtchedVM2000AirReflection:
265 run->AddEtchedVM2000AirReflect
266 break;
267 case EtchedVM2000GlueReflection:
268 run->AddEtchedVM2000AirReflect
269 break;
270 case GroundLumirrorAirReflection
271 run->AddGroundLumirrorAirRefle
272 break;
273 case GroundLumirrorGlueReflectio
274 run->AddGroundLumirrorGlueRefl
275 break;
276 case GroundAirReflection:
277 run->AddGroundAirReflection();
278 break;
279 case GroundTeflonAirReflection:
280 run->AddGroundTeflonAirReflect
281 break;
282 case GroundTiOAirReflection:
283 run->AddGroundTiOAirReflection
284 break;
285 case GroundTyvekAirReflection:
286 run->AddGroundTyvekAirReflecti
287 break;
288 case GroundVM2000AirReflection:
289 run->AddGroundVM2000AirReflect
290 break;
291 case GroundVM2000GlueReflection:
292 run->AddGroundVM2000AirReflect
293 break;
294 case Dichroic:
295 run->AddDichroic();
296 break;
297 case CoatedDielectricReflection:
298 run->AddCoatedDielectricReflec
299 break;
300 case CoatedDielectricRefraction:
301 run->AddCoatedDielectricRefrac
302 break;
303 case CoatedDielectricFrustratedT
304 run->AddCoatedDielectricFrustr
305 break;
306 default:
307 G4cout << "theStatus: " << the
308 break;
309 }
310 }
311 }
312 }
313 // when studying boundary scattering, it
314 // visualize the photon before and after
315 // selected, kill the photon when reachi
316 // (note that there are 2 steps at the b
317 // equals 0 and 1 on the first surface)
318 if (fKillOnSecondSurface) {
319 if (trackInfo->GetReflectionNumber() >
320 track->SetTrackStatus(fStopAndKill);
321 }
322 }
323 trackInfo->IncrementReflectionNumber();
324 }
325
326 // This block serves to test that G4OpBoun
327 // velocity correctly. It is not necessary
328 // Only steps where pre- and post- are the
329 // incorrect checks (so, in practice, set
330 // for particles to step in the interior o
331 if (endPoint->GetMaterial() == startPoint-
332 G4double trackVelocity = track->GetVeloc
333 G4double materialVelocity = CLHEP::c_lig
334 G4MaterialPropertyVector* velVector =
335 endPoint->GetMaterial()->GetMaterialPr
336 if (velVector) {
337 materialVelocity = velVector->Value(th
338 }
339
340 if (std::abs(trackVelocity - materialVel
341 G4ExceptionDescription ed;
342 ed << "Optical photon group velocity:
343 << " cm/ns is not what is expected
344 << materialVelocity / (cm / ns) <<
345 G4Exception("OpNovice2 SteppingAction"
346 }
347 }
348 // end of group velocity test
349 }
350
351 else { // particle != opticalphoton
352 // print how many Cerenkov and scint photo
353 // this demonstrates use of GetNumPhotons(
354 auto proc_man = track->GetDynamicParticle(
355 G4ProcessVector* proc_vec = proc_man->GetP
356 G4int n_proc = proc_vec->entries();
357
358 G4int n_scint = 0;
359 G4int n_cer = 0;
360 for (G4int i = 0; i < n_proc; ++i) {
361 G4String proc_name = (*proc_vec)[i]->Get
362 if (proc_name == "Cerenkov") {
363 auto cer = (G4Cerenkov*)(*proc_vec)[i]
364 n_cer = cer->GetNumPhotons();
365 }
366 else if (proc_name == "Scintillation") {
367 auto scint = (G4Scintillation*)(*proc_
368 n_scint = scint->GetNumPhotons();
369 }
370 }
371 if (fVerbose > 0) {
372 if (n_cer > 0 || n_scint > 0) {
373 G4cout << "In this step, " << n_cer <<
374 << " scintillation photons were
375 }
376 }
377
378 // loop over secondaries, create statistic
379 const std::vector<const G4Track*>* seconda
380
381 for (auto sec : *secondaries) {
382 if (sec->GetDynamicParticle()->GetPartic
383 G4String creator_process = sec->GetCre
384 if (creator_process == "Cerenkov") {
385 G4double en = sec->GetKineticEnergy(
386 run->AddCerenkovEnergy(en);
387 run->AddCerenkov();
388 analysisMan->FillH1(1, en / eV);
389 }
390 else if (creator_process == "Scintilla
391 G4double en = sec->GetKineticEnergy(
392 run->AddScintillationEnergy(en);
393 run->AddScintillation();
394 analysisMan->FillH1(2, en / eV);
395
396 G4double time = sec->GetGlobalTime()
397 analysisMan->FillH1(3, time / ns);
398 }
399 }
400 }
401 }
402
403 return;
404 }
405
406 //....oooOO0OOooo........oooOO0OOooo........oo
407