| Geant4 Cross Reference |
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 /// \file optical/OpNovice2/src/SteppingAction.cc
28 /// \brief Implementation of the SteppingAction class
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
53 SteppingAction::SteppingAction() : G4UserSteppingAction()
54 {
55 fSteppingMessenger = new SteppingMessenger(this);
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 SteppingAction::~SteppingAction()
60 {
61 delete fSteppingMessenger;
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 void SteppingAction::UserSteppingAction(const G4Step* step)
66 {
67 static G4ParticleDefinition* opticalphoton = G4OpticalPhoton::OpticalPhotonDefinition();
68
69 G4AnalysisManager* analysisMan = G4AnalysisManager::Instance();
70 Run* run = static_cast<Run*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
71
72 G4Track* track = step->GetTrack();
73 G4StepPoint* endPoint = step->GetPostStepPoint();
74 G4StepPoint* startPoint = step->GetPreStepPoint();
75
76 const G4DynamicParticle* theParticle = track->GetDynamicParticle();
77 const G4ParticleDefinition* particleDef = theParticle->GetParticleDefinition();
78
79 auto trackInfo = (TrackInformation*)(track->GetUserInformation());
80
81 if (particleDef == opticalphoton) {
82 const G4VProcess* pds = endPoint->GetProcessDefinedStep();
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); // absorption energy
98 // loop over secondaries, create statistics
99 // const std::vector<const G4Track*>* secondaries =
100 auto secondaries = step->GetSecondaryInCurrentStep();
101 for (auto sec : *secondaries) {
102 en = sec->GetKineticEnergy();
103 run->AddWLSEmission();
104 run->AddWLSEmissionEnergy(en);
105 analysisMan->FillH1(5, en / eV); // emission energy
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); // absorption energy
115 // loop over secondaries, create statistics
116 // const std::vector<const G4Track*>* secondaries =
117 auto secondaries = step->GetSecondaryInCurrentStep();
118 for (auto sec : *secondaries) {
119 en = sec->GetKineticEnergy();
120 run->AddWLS2Emission();
121 run->AddWLS2EmissionEnergy(en);
122 analysisMan->FillH1(8, en / eV); // emission energy
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() == fGeomBoundary) {
130 G4ThreeVector p0 = startPoint->GetMomentumDirection();
131 G4ThreeVector p1 = endPoint->GetMomentumDirection();
132
133 G4OpBoundaryProcessStatus theStatus = Undefined;
134
135 G4ProcessManager* OpManager = opticalphoton->GetProcessManager();
136 G4ProcessVector* postStepDoItVector = OpManager->GetPostStepProcessVector(typeDoIt);
137 G4int n_proc = postStepDoItVector->entries();
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 photons here
144 if (track->GetTrackStatus() != fStopAndKill) {
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 = (*postStepDoItVector)[i];
162
163 auto opProc = dynamic_cast<G4OpBoundaryProcess*>(currentProcess);
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 / deg);
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 / deg);
179 break;
180 case FresnelReflection:
181 run->AddFresnelReflection();
182 analysisMan->FillH1(21, angle / deg);
183 analysisMan->FillH1(23, angle / deg);
184 break;
185 case TotalInternalReflection:
186 run->AddTotalInternalReflection();
187 analysisMan->FillH1(22, angle / deg);
188 analysisMan->FillH1(23, angle / deg);
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 / deg);
199 break;
200 case BackScattering:
201 run->AddBackScattering();
202 break;
203 case Absorption:
204 analysisMan->FillH1(24, angle / deg);
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 PolishedLumirrorAirReflection:
223 run->AddPolishedLumirrorAirReflection();
224 break;
225 case PolishedLumirrorGlueReflection:
226 run->AddPolishedLumirrorGlueReflection();
227 break;
228 case PolishedAirReflection:
229 run->AddPolishedAirReflection();
230 break;
231 case PolishedTeflonAirReflection:
232 run->AddPolishedTeflonAirReflection();
233 break;
234 case PolishedTiOAirReflection:
235 run->AddPolishedTiOAirReflection();
236 break;
237 case PolishedTyvekAirReflection:
238 run->AddPolishedTyvekAirReflection();
239 break;
240 case PolishedVM2000AirReflection:
241 run->AddPolishedVM2000AirReflection();
242 break;
243 case PolishedVM2000GlueReflection:
244 run->AddPolishedVM2000AirReflection();
245 break;
246 case EtchedLumirrorAirReflection:
247 run->AddEtchedLumirrorAirReflection();
248 break;
249 case EtchedLumirrorGlueReflection:
250 run->AddEtchedLumirrorGlueReflection();
251 break;
252 case EtchedAirReflection:
253 run->AddEtchedAirReflection();
254 break;
255 case EtchedTeflonAirReflection:
256 run->AddEtchedTeflonAirReflection();
257 break;
258 case EtchedTiOAirReflection:
259 run->AddEtchedTiOAirReflection();
260 break;
261 case EtchedTyvekAirReflection:
262 run->AddEtchedTyvekAirReflection();
263 break;
264 case EtchedVM2000AirReflection:
265 run->AddEtchedVM2000AirReflection();
266 break;
267 case EtchedVM2000GlueReflection:
268 run->AddEtchedVM2000AirReflection();
269 break;
270 case GroundLumirrorAirReflection:
271 run->AddGroundLumirrorAirReflection();
272 break;
273 case GroundLumirrorGlueReflection:
274 run->AddGroundLumirrorGlueReflection();
275 break;
276 case GroundAirReflection:
277 run->AddGroundAirReflection();
278 break;
279 case GroundTeflonAirReflection:
280 run->AddGroundTeflonAirReflection();
281 break;
282 case GroundTiOAirReflection:
283 run->AddGroundTiOAirReflection();
284 break;
285 case GroundTyvekAirReflection:
286 run->AddGroundTyvekAirReflection();
287 break;
288 case GroundVM2000AirReflection:
289 run->AddGroundVM2000AirReflection();
290 break;
291 case GroundVM2000GlueReflection:
292 run->AddGroundVM2000AirReflection();
293 break;
294 case Dichroic:
295 run->AddDichroic();
296 break;
297 case CoatedDielectricReflection:
298 run->AddCoatedDielectricReflection();
299 break;
300 case CoatedDielectricRefraction:
301 run->AddCoatedDielectricRefraction();
302 break;
303 case CoatedDielectricFrustratedTransmission:
304 run->AddCoatedDielectricFrustratedTransmission();
305 break;
306 default:
307 G4cout << "theStatus: " << theStatus << " was none of the above." << G4endl;
308 break;
309 }
310 }
311 }
312 }
313 // when studying boundary scattering, it can be useful to only
314 // visualize the photon before and after the first surface. If
315 // selected, kill the photon when reaching the second surface
316 // (note that there are 2 steps at the boundary, so the counter
317 // equals 0 and 1 on the first surface)
318 if (fKillOnSecondSurface) {
319 if (trackInfo->GetReflectionNumber() >= 2) {
320 track->SetTrackStatus(fStopAndKill);
321 }
322 }
323 trackInfo->IncrementReflectionNumber();
324 }
325
326 // This block serves to test that G4OpBoundaryProcess sets the group
327 // velocity correctly. It is not necessary to include in user code.
328 // Only steps where pre- and post- are the same material, to avoid
329 // incorrect checks (so, in practice, set e.g. OpRayleigh low enough
330 // for particles to step in the interior of each volume.
331 if (endPoint->GetMaterial() == startPoint->GetMaterial()) {
332 G4double trackVelocity = track->GetVelocity();
333 G4double materialVelocity = CLHEP::c_light;
334 G4MaterialPropertyVector* velVector =
335 endPoint->GetMaterial()->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
336 if (velVector) {
337 materialVelocity = velVector->Value(theParticle->GetTotalMomentum(), fIdxVelocity);
338 }
339
340 if (std::abs(trackVelocity - materialVelocity) > 1e-9 * CLHEP::c_light) {
341 G4ExceptionDescription ed;
342 ed << "Optical photon group velocity: " << trackVelocity / (cm / ns)
343 << " cm/ns is not what is expected from " << G4endl << "the material properties, "
344 << materialVelocity / (cm / ns) << " cm/ns";
345 G4Exception("OpNovice2 SteppingAction", "OpNovice2_1", FatalException, ed);
346 }
347 }
348 // end of group velocity test
349 }
350
351 else { // particle != opticalphoton
352 // print how many Cerenkov and scint photons produced this step
353 // this demonstrates use of GetNumPhotons()
354 auto proc_man = track->GetDynamicParticle()->GetParticleDefinition()->GetProcessManager();
355 G4ProcessVector* proc_vec = proc_man->GetPostStepProcessVector(typeDoIt);
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]->GetProcessName();
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_vec)[i];
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 << " Cerenkov and " << n_scint
374 << " scintillation photons were produced." << G4endl;
375 }
376 }
377
378 // loop over secondaries, create statistics
379 const std::vector<const G4Track*>* secondaries = step->GetSecondaryInCurrentStep();
380
381 for (auto sec : *secondaries) {
382 if (sec->GetDynamicParticle()->GetParticleDefinition() == opticalphoton) {
383 G4String creator_process = sec->GetCreatorProcess()->GetProcessName();
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 == "Scintillation") {
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........oooOO0OOooo........oooOO0OOooo......
407