| 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 /// \file optical/OpNovice2/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33 #include "Run.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37
38 #include "G4OpBoundaryProcess.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4UnitsTable.hh"
41
42 #include <numeric>
43
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 Run::Run() : G4Run()
46 {
47 fBoundaryProcs.assign(43, 0);
48 }
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy, G4bool polarized,
52 G4double polarization)
53 {
54 fParticle = particle;
55 fEkin = energy;
56 fPolarized = polarized;
57 fPolarization = polarization;
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 void Run::Merge(const G4Run* run)
62 {
63 const Run* localRun = static_cast<const Run*>(run);
64
65 // pass information about primary particle
66 fParticle = localRun->fParticle;
67 fEkin = localRun->fEkin;
68 fPolarized = localRun->fPolarized;
69 fPolarization = localRun->fPolarization;
70
71 fCerenkovEnergy += localRun->fCerenkovEnergy;
72 fScintEnergy += localRun->fScintEnergy;
73 fWLSAbsorptionEnergy += localRun->fWLSAbsorptionEnergy;
74 fWLSEmissionEnergy += localRun->fWLSEmissionEnergy;
75 fWLS2AbsorptionEnergy += localRun->fWLS2AbsorptionEnergy;
76 fWLS2EmissionEnergy += localRun->fWLS2EmissionEnergy;
77
78 fCerenkovCount += localRun->fCerenkovCount;
79 fScintCount += localRun->fScintCount;
80 fWLSAbsorptionCount += localRun->fWLSAbsorptionCount;
81 fWLSEmissionCount += localRun->fWLSEmissionCount;
82 fWLS2AbsorptionCount += localRun->fWLS2AbsorptionCount;
83 fWLS2EmissionCount += localRun->fWLS2EmissionCount;
84 fRayleighCount += localRun->fRayleighCount;
85
86 fTotalSurface += localRun->fTotalSurface;
87
88 fOpAbsorption += localRun->fOpAbsorption;
89 fOpAbsorptionPrior += localRun->fOpAbsorptionPrior;
90
91 for (size_t i = 0; i < fBoundaryProcs.size(); ++i) {
92 fBoundaryProcs[i] += localRun->fBoundaryProcs[i];
93 }
94
95 G4Run::Merge(run);
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 void Run::EndOfRun()
100 {
101 if (numberOfEvent == 0) return;
102 auto TotNbofEvents = (G4double)numberOfEvent;
103
104 G4AnalysisManager* analysisMan = G4AnalysisManager::Instance();
105 G4int id = analysisMan->GetH1Id("Cerenkov spectrum");
106 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
107 analysisMan->SetH1YAxisTitle(id, "Number of photons");
108
109 id = analysisMan->GetH1Id("Scintillation spectrum");
110 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
111 analysisMan->SetH1YAxisTitle(id, "Number of photons");
112
113 id = analysisMan->GetH1Id("Scintillation time");
114 analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
115 analysisMan->SetH1YAxisTitle(id, "Number of photons");
116
117 id = analysisMan->GetH1Id("WLS abs");
118 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
119 analysisMan->SetH1YAxisTitle(id, "Number of photons");
120
121 id = analysisMan->GetH1Id("WLS em");
122 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
123 analysisMan->SetH1YAxisTitle(id, "Number of photons");
124
125 id = analysisMan->GetH1Id("WLS time");
126 analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
127 analysisMan->SetH1YAxisTitle(id, "Number of photons");
128
129 id = analysisMan->GetH1Id("WLS2 abs");
130 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
131 analysisMan->SetH1YAxisTitle(id, "Number of photons");
132
133 id = analysisMan->GetH1Id("WLS2 em");
134 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
135 analysisMan->SetH1YAxisTitle(id, "Number of photons");
136
137 id = analysisMan->GetH1Id("WLS2 time");
138 analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
139 analysisMan->SetH1YAxisTitle(id, "Number of photons");
140
141 id = analysisMan->GetH1Id("bdry status");
142 analysisMan->SetH1XAxisTitle(id, "Status code");
143 analysisMan->SetH1YAxisTitle(id, "Number of photons");
144
145 id = analysisMan->GetH1Id("x_backward");
146 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
147 analysisMan->SetH1YAxisTitle(id, "Number of photons");
148
149 id = analysisMan->GetH1Id("y_backward");
150 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
151 analysisMan->SetH1YAxisTitle(id, "Number of photons");
152
153 id = analysisMan->GetH1Id("z_backward");
154 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
155 analysisMan->SetH1YAxisTitle(id, "Number of photons");
156
157 id = analysisMan->GetH1Id("x_forward");
158 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
159 analysisMan->SetH1YAxisTitle(id, "Number of photons");
160
161 id = analysisMan->GetH1Id("y_forward");
162 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
163 analysisMan->SetH1YAxisTitle(id, "Number of photons");
164
165 id = analysisMan->GetH1Id("z_forward");
166 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
167 analysisMan->SetH1YAxisTitle(id, "Number of photons");
168
169 id = analysisMan->GetH1Id("x_fresnel");
170 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
171 analysisMan->SetH1YAxisTitle(id, "Number of photons");
172
173 id = analysisMan->GetH1Id("y_fresnel");
174 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
175 analysisMan->SetH1YAxisTitle(id, "Number of photons");
176
177 id = analysisMan->GetH1Id("z_fresnel");
178 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
179 analysisMan->SetH1YAxisTitle(id, "Number of photons");
180
181 id = analysisMan->GetH1Id("Fresnel reflection");
182 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
183 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
184
185 id = analysisMan->GetH1Id("Fresnel refraction");
186 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
187 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
188
189 id = analysisMan->GetH1Id("Total internal reflection");
190 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
191 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
192
193 id = analysisMan->GetH1Id("Fresnel reflection plus TIR");
194 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
195 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
196
197 id = analysisMan->GetH1Id("Absorption");
198 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
199 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
200
201 id = analysisMan->GetH1Id("Transmitted");
202 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
203 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
204
205 id = analysisMan->GetH1Id("Spike reflection");
206 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
207 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
208
209 const auto det =
210 (const DetectorConstruction*)(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
211
212 std::ios::fmtflags mode = G4cout.flags();
213 G4int prec = G4cout.precision(2);
214
215 G4cout << "\n Run Summary\n";
216 G4cout << "---------------------------------\n";
217 G4cout << "Primary particle was: " << fParticle->GetParticleName() << " with energy "
218 << G4BestUnit(fEkin, "Energy") << "." << G4endl;
219 G4cout << "Number of events: " << numberOfEvent << G4endl;
220
221 G4cout << "Material of world: " << det->GetWorldMaterial()->GetName() << G4endl;
222 G4cout << "Material of tank: " << det->GetTankMaterial()->GetName() << G4endl << G4endl;
223
224 if (fParticle->GetParticleName() != "opticalphoton") {
225 G4cout << "Average energy of Cerenkov photons created per event: "
226 << (fCerenkovEnergy / eV) / TotNbofEvents << " eV." << G4endl;
227 G4cout << "Average number of Cerenkov photons created per event: "
228 << fCerenkovCount / TotNbofEvents << G4endl;
229 if (fCerenkovCount > 0) {
230 G4cout << " Average energy per photon: " << (fCerenkovEnergy / eV) / fCerenkovCount << " eV."
231 << G4endl;
232 }
233 G4cout << "Average energy of scintillation photons created per event: "
234 << (fScintEnergy / eV) / TotNbofEvents << " eV." << G4endl;
235 G4cout << "Average number of scintillation photons created per event: "
236 << fScintCount / TotNbofEvents << G4endl;
237 if (fScintCount > 0) {
238 G4cout << " Average energy per photon: " << (fScintEnergy / eV) / fScintCount << " eV."
239 << G4endl;
240 }
241 }
242
243 G4cout << "Average number of photons absorbed by WLS per event: "
244 << fWLSAbsorptionCount / G4double(TotNbofEvents) << " " << G4endl;
245 if (fWLSAbsorptionCount > 0) {
246 G4cout << " Average energy per photon: " << (fWLSAbsorptionEnergy / eV) / fWLSAbsorptionCount
247 << " eV." << G4endl;
248 }
249 G4cout << "Average number of photons created by WLS per event: "
250 << fWLSEmissionCount / TotNbofEvents << G4endl;
251 if (fWLSEmissionCount > 0) {
252 G4cout << " Average energy per photon: " << (fWLSEmissionEnergy / eV) / fWLSEmissionCount
253 << " eV." << G4endl;
254 }
255 G4cout << "Average energy of WLS photons created per event: "
256 << (fWLSEmissionEnergy / eV) / TotNbofEvents << " eV." << G4endl;
257
258 G4cout << "Average number of photons absorbed by WLS2 per event: "
259 << fWLS2AbsorptionCount / G4double(TotNbofEvents) << " " << G4endl;
260 if (fWLS2AbsorptionCount > 0) {
261 G4cout << " Average energy per photon: " << (fWLS2AbsorptionEnergy / eV) / fWLS2AbsorptionCount
262 << " eV." << G4endl;
263 }
264 G4cout << "Average number of photons created by WLS2 per event: "
265 << fWLS2EmissionCount / TotNbofEvents << G4endl;
266 if (fWLS2EmissionCount > 0) {
267 G4cout << " Average energy per photon: " << (fWLS2EmissionEnergy / eV) / fWLS2EmissionCount
268 << " eV." << G4endl;
269 }
270 G4cout << "Average energy of WLS2 photons created per event: "
271 << (fWLS2EmissionEnergy / eV) / TotNbofEvents << " eV." << G4endl;
272
273 G4cout << "Average number of OpRayleigh per event: " << fRayleighCount / TotNbofEvents
274 << G4endl;
275 G4cout << "Average number of OpAbsorption per event: " << fOpAbsorption / TotNbofEvents << G4endl;
276 G4cout << "\nSurface events (on +X surface, maximum one per photon) this run:" << G4endl;
277 G4cout << "# of primary particles: " << std::setw(8) << TotNbofEvents << G4endl;
278 G4cout << "OpAbsorption before surface: " << std::setw(8) << fOpAbsorptionPrior << G4endl;
279 G4cout << "Total # of surface events: " << std::setw(8) << fTotalSurface << G4endl;
280 if (fParticle->GetParticleName() == "opticalphoton") {
281 G4cout << "Unaccounted for: " << std::setw(8)
282 << fTotalSurface + fOpAbsorptionPrior - TotNbofEvents << G4endl;
283 }
284 G4cout << "\nSurface events by process:" << G4endl;
285 if (fBoundaryProcs[Transmission] > 0) {
286 G4cout << " Transmission: " << std::setw(8) << fBoundaryProcs[Transmission]
287 << G4endl;
288 }
289 if (fBoundaryProcs[FresnelRefraction] > 0) {
290 G4cout << " Fresnel refraction: " << std::setw(8) << fBoundaryProcs[FresnelRefraction]
291 << G4endl;
292 }
293 if (fBoundaryProcs[FresnelReflection] > 0) {
294 G4cout << " Fresnel reflection: " << std::setw(8) << fBoundaryProcs[FresnelReflection]
295 << G4endl;
296 }
297 if (fBoundaryProcs[TotalInternalReflection] > 0) {
298 G4cout << " Total internal reflection: " << std::setw(8)
299 << fBoundaryProcs[TotalInternalReflection] << G4endl;
300 }
301 if (fBoundaryProcs[LambertianReflection] > 0) {
302 G4cout << " Lambertian reflection: " << std::setw(8)
303 << fBoundaryProcs[LambertianReflection] << G4endl;
304 }
305 if (fBoundaryProcs[LobeReflection] > 0) {
306 G4cout << " Lobe reflection: " << std::setw(8) << fBoundaryProcs[LobeReflection]
307 << G4endl;
308 }
309 if (fBoundaryProcs[SpikeReflection] > 0) {
310 G4cout << " Spike reflection: " << std::setw(8) << fBoundaryProcs[SpikeReflection]
311 << G4endl;
312 }
313 if (fBoundaryProcs[BackScattering] > 0) {
314 G4cout << " Backscattering: " << std::setw(8) << fBoundaryProcs[BackScattering]
315 << G4endl;
316 }
317 if (fBoundaryProcs[Absorption] > 0) {
318 G4cout << " Absorption: " << std::setw(8) << fBoundaryProcs[Absorption]
319 << G4endl;
320 }
321 if (fBoundaryProcs[Detection] > 0) {
322 G4cout << " Detection: " << std::setw(8) << fBoundaryProcs[Detection]
323 << G4endl;
324 }
325 if (fBoundaryProcs[NotAtBoundary] > 0) {
326 G4cout << " Not at boundary: " << std::setw(8) << fBoundaryProcs[NotAtBoundary]
327 << G4endl;
328 }
329 if (fBoundaryProcs[SameMaterial] > 0) {
330 G4cout << " Same material: " << std::setw(8) << fBoundaryProcs[SameMaterial]
331 << G4endl;
332 }
333 if (fBoundaryProcs[StepTooSmall] > 0) {
334 G4cout << " Step too small: " << std::setw(8) << fBoundaryProcs[StepTooSmall]
335 << G4endl;
336 }
337 if (fBoundaryProcs[NoRINDEX] > 0) {
338 G4cout << " No RINDEX: " << std::setw(8) << fBoundaryProcs[NoRINDEX] << G4endl;
339 }
340 // LBNL polished
341 if (fBoundaryProcs[PolishedLumirrorAirReflection] > 0) {
342 G4cout << " Polished Lumirror Air reflection: " << std::setw(8)
343 << fBoundaryProcs[PolishedLumirrorAirReflection] << G4endl;
344 }
345 if (fBoundaryProcs[PolishedLumirrorGlueReflection] > 0) {
346 G4cout << " Polished Lumirror Glue reflection: " << std::setw(8)
347 << fBoundaryProcs[PolishedLumirrorGlueReflection] << G4endl;
348 }
349 if (fBoundaryProcs[PolishedAirReflection] > 0) {
350 G4cout << " Polished Air reflection: " << std::setw(8) << fBoundaryProcs[PolishedAirReflection]
351 << G4endl;
352 }
353 if (fBoundaryProcs[PolishedTeflonAirReflection] > 0) {
354 G4cout << " Polished Teflon Air reflection: " << std::setw(8)
355 << fBoundaryProcs[PolishedTeflonAirReflection] << G4endl;
356 }
357 if (fBoundaryProcs[PolishedTiOAirReflection] > 0) {
358 G4cout << " Polished TiO Air reflection: " << std::setw(8)
359 << fBoundaryProcs[PolishedTiOAirReflection] << G4endl;
360 }
361 if (fBoundaryProcs[PolishedTyvekAirReflection] > 0) {
362 G4cout << " Polished Tyvek Air reflection: " << std::setw(8)
363 << fBoundaryProcs[PolishedTyvekAirReflection] << G4endl;
364 }
365 if (fBoundaryProcs[PolishedVM2000AirReflection] > 0) {
366 G4cout << " Polished VM2000 Air reflection: " << std::setw(8)
367 << fBoundaryProcs[PolishedVM2000AirReflection] << G4endl;
368 }
369 if (fBoundaryProcs[PolishedVM2000GlueReflection] > 0) {
370 G4cout << " Polished VM2000 Glue reflection: " << std::setw(8)
371 << fBoundaryProcs[PolishedVM2000GlueReflection] << G4endl;
372 }
373 // LBNL etched
374 if (fBoundaryProcs[EtchedLumirrorAirReflection] > 0) {
375 G4cout << " Etched Lumirror Air reflection: " << std::setw(8)
376 << fBoundaryProcs[EtchedLumirrorAirReflection] << G4endl;
377 }
378 if (fBoundaryProcs[EtchedLumirrorGlueReflection] > 0) {
379 G4cout << " Etched Lumirror Glue reflection: " << std::setw(8)
380 << fBoundaryProcs[EtchedLumirrorGlueReflection] << G4endl;
381 }
382 if (fBoundaryProcs[EtchedAirReflection] > 0) {
383 G4cout << " Etched Air reflection: " << std::setw(8) << fBoundaryProcs[EtchedAirReflection]
384 << G4endl;
385 }
386 if (fBoundaryProcs[EtchedTeflonAirReflection] > 0) {
387 G4cout << " Etched Teflon Air reflection: " << std::setw(8)
388 << fBoundaryProcs[EtchedTeflonAirReflection] << G4endl;
389 }
390 if (fBoundaryProcs[EtchedTiOAirReflection] > 0) {
391 G4cout << " Etched TiO Air reflection: " << std::setw(8)
392 << fBoundaryProcs[EtchedTiOAirReflection] << G4endl;
393 }
394 if (fBoundaryProcs[EtchedTyvekAirReflection] > 0) {
395 G4cout << " Etched Tyvek Air reflection: " << std::setw(8)
396 << fBoundaryProcs[EtchedTyvekAirReflection] << G4endl;
397 }
398 if (fBoundaryProcs[EtchedVM2000AirReflection] > 0) {
399 G4cout << " Etched VM2000 Air reflection: " << std::setw(8)
400 << fBoundaryProcs[EtchedVM2000AirReflection] << G4endl;
401 }
402 if (fBoundaryProcs[EtchedVM2000GlueReflection] > 0) {
403 G4cout << " Etched VM2000 Glue reflection: " << std::setw(8)
404 << fBoundaryProcs[EtchedVM2000GlueReflection] << G4endl;
405 }
406 // LBNL ground
407 if (fBoundaryProcs[GroundLumirrorAirReflection] > 0) {
408 G4cout << " Ground Lumirror Air reflection: " << std::setw(8)
409 << fBoundaryProcs[GroundLumirrorAirReflection] << G4endl;
410 }
411 if (fBoundaryProcs[GroundLumirrorGlueReflection] > 0) {
412 G4cout << " Ground Lumirror Glue reflection: " << std::setw(8)
413 << fBoundaryProcs[GroundLumirrorGlueReflection] << G4endl;
414 }
415 if (fBoundaryProcs[GroundAirReflection] > 0) {
416 G4cout << " Ground Air reflection: " << std::setw(8) << fBoundaryProcs[GroundAirReflection]
417 << G4endl;
418 }
419 if (fBoundaryProcs[GroundTeflonAirReflection] > 0) {
420 G4cout << " Ground Teflon Air reflection: " << std::setw(8)
421 << fBoundaryProcs[GroundTeflonAirReflection] << G4endl;
422 }
423 if (fBoundaryProcs[GroundTiOAirReflection] > 0) {
424 G4cout << " Ground TiO Air reflection: " << std::setw(8)
425 << fBoundaryProcs[GroundTiOAirReflection] << G4endl;
426 }
427 if (fBoundaryProcs[GroundTyvekAirReflection] > 0) {
428 G4cout << " Ground Tyvek Air reflection: " << std::setw(8)
429 << fBoundaryProcs[GroundTyvekAirReflection] << G4endl;
430 }
431 if (fBoundaryProcs[GroundVM2000AirReflection] > 0) {
432 G4cout << " Ground VM2000 Air reflection: " << std::setw(8)
433 << fBoundaryProcs[GroundVM2000AirReflection] << G4endl;
434 }
435 if (fBoundaryProcs[GroundVM2000GlueReflection] > 0) {
436 G4cout << " Ground VM2000 Glue reflection: " << std::setw(8)
437 << fBoundaryProcs[GroundVM2000GlueReflection] << G4endl;
438 }
439 if (fBoundaryProcs[CoatedDielectricRefraction] > 0) {
440 G4cout << " CoatedDielectricRefraction: " << std::setw(8)
441 << fBoundaryProcs[CoatedDielectricRefraction] << G4endl;
442 }
443 if (fBoundaryProcs[CoatedDielectricReflection] > 0) {
444 G4cout << " CoatedDielectricReflection: " << std::setw(8)
445 << fBoundaryProcs[CoatedDielectricReflection] << G4endl;
446 }
447 if (fBoundaryProcs[CoatedDielectricFrustratedTransmission] > 0) {
448 G4cout << " CoatedDielectricFrustratedTransmission: " << std::setw(8)
449 << fBoundaryProcs[CoatedDielectricFrustratedTransmission] << G4endl;
450 }
451
452 G4int sum = std::accumulate(fBoundaryProcs.begin(), fBoundaryProcs.end(), 0);
453 G4cout << " Sum: " << std::setw(8) << sum << G4endl;
454 G4cout << " Unaccounted for: " << std::setw(8) << fTotalSurface - sum << G4endl;
455
456 G4cout << "---------------------------------\n";
457 G4cout.setf(mode, std::ios::floatfield);
458 G4cout.precision(prec);
459
460 G4int histo_id_refract = analysisMan->GetH1Id("Fresnel refraction");
461 G4int histo_id_reflect = analysisMan->GetH1Id("Fresnel reflection plus TIR");
462 G4int histo_id_spike = analysisMan->GetH1Id("Spike reflection");
463 G4int histo_id_absorption = analysisMan->GetH1Id("Absorption");
464
465 if (analysisMan->GetH1Activation(histo_id_refract)
466 && analysisMan->GetH1Activation(histo_id_reflect))
467 {
468 G4double rindex1 =
469 det->GetTankMaterial()->GetMaterialPropertiesTable()->GetProperty(kRINDEX)->Value(fEkin);
470 G4double rindex2 =
471 det->GetWorldMaterial()->GetMaterialPropertiesTable()->GetProperty(kRINDEX)->Value(fEkin);
472
473 auto histo_refract = analysisMan->GetH1(histo_id_refract);
474 auto histo_reflect = analysisMan->GetH1(histo_id_reflect);
475 // std::vector<G4double> refract;
476 std::vector<G4double> reflect;
477 // std::vector<G4double> tir;
478 std::vector<G4double> tot;
479 for (size_t i = 0; i < histo_refract->axis().bins(); ++i) {
480 // refract.push_back(histo_refract->bin_height(i));
481 reflect.push_back(histo_reflect->bin_height(i));
482 // tir.push_back(histo_TIR->bin_height(i));
483 tot.push_back(histo_refract->bin_height(i) + histo_reflect->bin_height(i));
484 }
485
486 // find Brewster angle: Rp = 0
487 // need enough statistics for this method to work
488 G4double min_angle = -1.;
489 G4double min_val = DBL_MAX;
490 G4double bin_width = 0.;
491 for (size_t i = 0; i < reflect.size(); ++i) {
492 if (reflect[i] < min_val) {
493 min_val = reflect[i];
494 min_angle = histo_reflect->axis().bin_lower_edge(i);
495 bin_width =
496 histo_reflect->axis().bin_upper_edge(i) - histo_reflect->axis().bin_lower_edge(i);
497 min_angle += bin_width / 2.;
498 }
499 }
500 G4cout << "Polarization of primary optical photons: " << fPolarization / deg << " deg."
501 << G4endl;
502 if (fPolarized && fPolarization == 0.0) {
503 G4cout << "Reflectance shows a minimum at: " << min_angle << " +/- " << bin_width / 2;
504 G4cout << " deg. Expected Brewster angle: "
505 << (360. / CLHEP::twopi) * std::atan(rindex2 / rindex1) << " deg. " << G4endl;
506 }
507
508 // find angle of total internal reflection: T -> 0
509 // last bin for T > 0
510 min_angle = -1.;
511 min_val = DBL_MAX;
512 for (size_t i = 0; i < histo_refract->axis().bins() - 1; ++i) {
513 if (histo_refract->bin_height(i) > 0. && histo_refract->bin_height(i + 1) == 0.) {
514 min_angle = histo_refract->axis().bin_lower_edge(i);
515 bin_width =
516 histo_reflect->axis().bin_upper_edge(i) - histo_reflect->axis().bin_lower_edge(i);
517 min_angle += bin_width / 2.;
518 break;
519 }
520 }
521 if (fPolarized) {
522 G4cout << "Fresnel transmission goes to 0 at: " << min_angle << " +/- " << bin_width / 2.
523 << " deg."
524 << " Expected: " << (360. / CLHEP::twopi) * std::asin(rindex2 / rindex1) << " deg."
525 << G4endl;
526 }
527
528 // Normalize the transmission/reflection histos so that max is 1.
529 // Only if x values are the same
530 if ((analysisMan->GetH1Nbins(histo_id_refract) == analysisMan->GetH1Nbins(histo_id_reflect))
531 && (analysisMan->GetH1Xmin(histo_id_refract) == analysisMan->GetH1Xmin(histo_id_reflect))
532 && (analysisMan->GetH1Xmax(histo_id_refract) == analysisMan->GetH1Xmax(histo_id_reflect)))
533 {
534 unsigned int ent;
535 G4double sw;
536 G4double sw2;
537 G4double sx2;
538 G4double sx2w;
539 for (size_t bin = 0; bin < histo_refract->axis().bins(); ++bin) {
540 // "bin+1" below because bin 0 is underflow bin
541 // NB. We are ignoring underflow/overflow bins
542 histo_refract->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
543 if (tot[bin] > 0) {
544 sw /= tot[bin];
545 // bin error is sqrt(sw2)
546 sw2 /= (tot[bin] * tot[bin]);
547 sx2 /= (tot[bin] * tot[bin]);
548 sx2w /= (tot[bin] * tot[bin]);
549 histo_refract->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
550 }
551
552 histo_reflect->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
553 if (tot[bin] > 0) {
554 sw /= tot[bin];
555 // bin error is sqrt(sw2)
556 sw2 /= (tot[bin] * tot[bin]);
557 sx2 /= (tot[bin] * tot[bin]);
558 sx2w /= (tot[bin] * tot[bin]);
559 histo_reflect->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
560 }
561
562 G4int histo_id_fresnelrefl = analysisMan->GetH1Id("Fresnel reflection");
563 auto histo_fresnelreflect = analysisMan->GetH1(histo_id_fresnelrefl);
564 histo_fresnelreflect->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
565 if (tot[bin] > 0) {
566 sw /= tot[bin];
567 // bin error is sqrt(sw2)
568 sw2 /= (tot[bin] * tot[bin]);
569 sx2 /= (tot[bin] * tot[bin]);
570 sx2w /= (tot[bin] * tot[bin]);
571 histo_fresnelreflect->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
572 }
573
574 G4int histo_id_TIR = analysisMan->GetH1Id("Total internal reflection");
575 auto histo_TIR = analysisMan->GetH1(histo_id_TIR);
576 if (analysisMan->GetH1Activation(histo_id_TIR)) {
577 histo_TIR->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
578 if (tot[bin] > 0) {
579 sw /= tot[bin];
580 // bin error is sqrt(sw2)
581 sw2 /= (tot[bin] * tot[bin]);
582 sx2 /= (tot[bin] * tot[bin]);
583 sx2w /= (tot[bin] * tot[bin]);
584 histo_TIR->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
585 }
586 }
587 }
588 }
589 else {
590 G4cout << "Not going to normalize refraction and reflection "
591 << "histograms because bins are not the same." << G4endl;
592 }
593 }
594
595 // complex index of refraction; have spike reflection and absorption
596 // Only works for polished surfaces. Ground surfaces neglected.
597 else if (analysisMan->GetH1Activation(histo_id_absorption)
598 && analysisMan->GetH1Activation(histo_id_spike))
599 {
600 auto histo_spike = analysisMan->GetH1(histo_id_spike);
601 auto histo_absorption = analysisMan->GetH1(histo_id_absorption);
602
603 std::vector<G4double> tot;
604 for (size_t i = 0; i < histo_absorption->axis().bins(); ++i) {
605 tot.push_back(histo_absorption->bin_height(i) + histo_spike->bin_height(i));
606 }
607
608 if ((analysisMan->GetH1Nbins(histo_id_absorption) == analysisMan->GetH1Nbins(histo_id_spike))
609 && (analysisMan->GetH1Xmin(histo_id_absorption) == analysisMan->GetH1Xmin(histo_id_spike))
610 && (analysisMan->GetH1Xmax(histo_id_absorption) == analysisMan->GetH1Xmax(histo_id_spike)))
611 {
612 unsigned int ent;
613 G4double sw;
614 G4double sw2;
615 G4double sx2;
616 G4double sx2w;
617 for (size_t bin = 0; bin < histo_absorption->axis().bins(); ++bin) {
618 // "bin+1" below because bin 0 is underflow bin
619 // NB. We are ignoring underflow/overflow bins
620 histo_absorption->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
621 if (tot[bin] > 0) {
622 sw /= tot[bin];
623 // bin error is sqrt(sw2)
624 sw2 /= (tot[bin] * tot[bin]);
625 sx2 /= (tot[bin] * tot[bin]);
626 sx2w /= (tot[bin] * tot[bin]);
627 histo_absorption->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
628 }
629
630 histo_spike->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
631 if (tot[bin] > 0) {
632 sw /= tot[bin];
633 // bin error is sqrt(sw2)
634 sw2 /= (tot[bin] * tot[bin]);
635 sx2 /= (tot[bin] * tot[bin]);
636 sx2w /= (tot[bin] * tot[bin]);
637 histo_spike->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
638 }
639 }
640 }
641 else {
642 G4cout << "Not going to normalize spike reflection and absorption "
643 << "histograms because bins are not the same." << G4endl;
644 }
645 }
646 }
647