| 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 /// \file optical/OpNovice2/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
45 Run::Run() : G4Run()
46 {
47 fBoundaryProcs.assign(43, 0);
48 }
49
50 //....oooOO0OOooo........oooOO0OOooo........oo
51 void Run::SetPrimary(G4ParticleDefinition* par
52 G4double polarization)
53 {
54 fParticle = particle;
55 fEkin = energy;
56 fPolarized = polarized;
57 fPolarization = polarization;
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61 void Run::Merge(const G4Run* run)
62 {
63 const Run* localRun = static_cast<const 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->fWLSAbsorp
74 fWLSEmissionEnergy += localRun->fWLSEmission
75 fWLS2AbsorptionEnergy += localRun->fWLS2Abso
76 fWLS2EmissionEnergy += localRun->fWLS2Emissi
77
78 fCerenkovCount += localRun->fCerenkovCount;
79 fScintCount += localRun->fScintCount;
80 fWLSAbsorptionCount += localRun->fWLSAbsorpt
81 fWLSEmissionCount += localRun->fWLSEmissionC
82 fWLS2AbsorptionCount += localRun->fWLS2Absor
83 fWLS2EmissionCount += localRun->fWLS2Emissio
84 fRayleighCount += localRun->fRayleighCount;
85
86 fTotalSurface += localRun->fTotalSurface;
87
88 fOpAbsorption += localRun->fOpAbsorption;
89 fOpAbsorptionPrior += localRun->fOpAbsorptio
90
91 for (size_t i = 0; i < fBoundaryProcs.size()
92 fBoundaryProcs[i] += localRun->fBoundaryPr
93 }
94
95 G4Run::Merge(run);
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oo
99 void Run::EndOfRun()
100 {
101 if (numberOfEvent == 0) return;
102 auto TotNbofEvents = (G4double)numberOfEvent
103
104 G4AnalysisManager* analysisMan = G4AnalysisM
105 G4int id = analysisMan->GetH1Id("Cerenkov sp
106 analysisMan->SetH1XAxisTitle(id, "Energy [eV
107 analysisMan->SetH1YAxisTitle(id, "Number of
108
109 id = analysisMan->GetH1Id("Scintillation spe
110 analysisMan->SetH1XAxisTitle(id, "Energy [eV
111 analysisMan->SetH1YAxisTitle(id, "Number of
112
113 id = analysisMan->GetH1Id("Scintillation tim
114 analysisMan->SetH1XAxisTitle(id, "Creation t
115 analysisMan->SetH1YAxisTitle(id, "Number of
116
117 id = analysisMan->GetH1Id("WLS abs");
118 analysisMan->SetH1XAxisTitle(id, "Energy [eV
119 analysisMan->SetH1YAxisTitle(id, "Number of
120
121 id = analysisMan->GetH1Id("WLS em");
122 analysisMan->SetH1XAxisTitle(id, "Energy [eV
123 analysisMan->SetH1YAxisTitle(id, "Number of
124
125 id = analysisMan->GetH1Id("WLS time");
126 analysisMan->SetH1XAxisTitle(id, "Creation t
127 analysisMan->SetH1YAxisTitle(id, "Number of
128
129 id = analysisMan->GetH1Id("WLS2 abs");
130 analysisMan->SetH1XAxisTitle(id, "Energy [eV
131 analysisMan->SetH1YAxisTitle(id, "Number of
132
133 id = analysisMan->GetH1Id("WLS2 em");
134 analysisMan->SetH1XAxisTitle(id, "Energy [eV
135 analysisMan->SetH1YAxisTitle(id, "Number of
136
137 id = analysisMan->GetH1Id("WLS2 time");
138 analysisMan->SetH1XAxisTitle(id, "Creation t
139 analysisMan->SetH1YAxisTitle(id, "Number of
140
141 id = analysisMan->GetH1Id("bdry status");
142 analysisMan->SetH1XAxisTitle(id, "Status cod
143 analysisMan->SetH1YAxisTitle(id, "Number of
144
145 id = analysisMan->GetH1Id("x_backward");
146 analysisMan->SetH1XAxisTitle(id, "Direction
147 analysisMan->SetH1YAxisTitle(id, "Number of
148
149 id = analysisMan->GetH1Id("y_backward");
150 analysisMan->SetH1XAxisTitle(id, "Direction
151 analysisMan->SetH1YAxisTitle(id, "Number of
152
153 id = analysisMan->GetH1Id("z_backward");
154 analysisMan->SetH1XAxisTitle(id, "Direction
155 analysisMan->SetH1YAxisTitle(id, "Number of
156
157 id = analysisMan->GetH1Id("x_forward");
158 analysisMan->SetH1XAxisTitle(id, "Direction
159 analysisMan->SetH1YAxisTitle(id, "Number of
160
161 id = analysisMan->GetH1Id("y_forward");
162 analysisMan->SetH1XAxisTitle(id, "Direction
163 analysisMan->SetH1YAxisTitle(id, "Number of
164
165 id = analysisMan->GetH1Id("z_forward");
166 analysisMan->SetH1XAxisTitle(id, "Direction
167 analysisMan->SetH1YAxisTitle(id, "Number of
168
169 id = analysisMan->GetH1Id("x_fresnel");
170 analysisMan->SetH1XAxisTitle(id, "Direction
171 analysisMan->SetH1YAxisTitle(id, "Number of
172
173 id = analysisMan->GetH1Id("y_fresnel");
174 analysisMan->SetH1XAxisTitle(id, "Direction
175 analysisMan->SetH1YAxisTitle(id, "Number of
176
177 id = analysisMan->GetH1Id("z_fresnel");
178 analysisMan->SetH1XAxisTitle(id, "Direction
179 analysisMan->SetH1YAxisTitle(id, "Number of
180
181 id = analysisMan->GetH1Id("Fresnel reflectio
182 analysisMan->SetH1XAxisTitle(id, "Angle [deg
183 analysisMan->SetH1YAxisTitle(id, "Fraction o
184
185 id = analysisMan->GetH1Id("Fresnel refractio
186 analysisMan->SetH1XAxisTitle(id, "Angle [deg
187 analysisMan->SetH1YAxisTitle(id, "Fraction o
188
189 id = analysisMan->GetH1Id("Total internal re
190 analysisMan->SetH1XAxisTitle(id, "Angle [deg
191 analysisMan->SetH1YAxisTitle(id, "Fraction o
192
193 id = analysisMan->GetH1Id("Fresnel reflectio
194 analysisMan->SetH1XAxisTitle(id, "Angle [deg
195 analysisMan->SetH1YAxisTitle(id, "Fraction o
196
197 id = analysisMan->GetH1Id("Absorption");
198 analysisMan->SetH1XAxisTitle(id, "Angle [deg
199 analysisMan->SetH1YAxisTitle(id, "Fraction o
200
201 id = analysisMan->GetH1Id("Transmitted");
202 analysisMan->SetH1XAxisTitle(id, "Angle [deg
203 analysisMan->SetH1YAxisTitle(id, "Fraction o
204
205 id = analysisMan->GetH1Id("Spike reflection"
206 analysisMan->SetH1XAxisTitle(id, "Angle [deg
207 analysisMan->SetH1YAxisTitle(id, "Fraction o
208
209 const auto det =
210 (const DetectorConstruction*)(G4RunManager
211
212 std::ios::fmtflags mode = G4cout.flags();
213 G4int prec = G4cout.precision(2);
214
215 G4cout << "\n Run Summary\n";
216 G4cout << "---------------------------------
217 G4cout << "Primary particle was: " << fParti
218 << G4BestUnit(fEkin, "Energy") << "."
219 G4cout << "Number of events: " << numberOfEv
220
221 G4cout << "Material of world: " << det->GetW
222 G4cout << "Material of tank: " << det->GetT
223
224 if (fParticle->GetParticleName() != "optical
225 G4cout << "Average energy of Cerenkov phot
226 << (fCerenkovEnergy / eV) / TotNbof
227 G4cout << "Average number of Cerenkov phot
228 << fCerenkovCount / TotNbofEvents <
229 if (fCerenkovCount > 0) {
230 G4cout << " Average energy per photon: "
231 << G4endl;
232 }
233 G4cout << "Average energy of scintillation
234 << (fScintEnergy / eV) / TotNbofEve
235 G4cout << "Average number of scintillation
236 << fScintCount / TotNbofEvents << G
237 if (fScintCount > 0) {
238 G4cout << " Average energy per photon: "
239 << G4endl;
240 }
241 }
242
243 G4cout << "Average number of photons absorbe
244 << fWLSAbsorptionCount / G4double(Tot
245 if (fWLSAbsorptionCount > 0) {
246 G4cout << " Average energy per photon: " <
247 << " eV." << G4endl;
248 }
249 G4cout << "Average number of photons created
250 << fWLSEmissionCount / TotNbofEvents
251 if (fWLSEmissionCount > 0) {
252 G4cout << " Average energy per photon: " <
253 << " eV." << G4endl;
254 }
255 G4cout << "Average energy of WLS photons cre
256 << (fWLSEmissionEnergy / eV) / TotNbo
257
258 G4cout << "Average number of photons absorbe
259 << fWLS2AbsorptionCount / G4double(To
260 if (fWLS2AbsorptionCount > 0) {
261 G4cout << " Average energy per photon: " <
262 << " eV." << G4endl;
263 }
264 G4cout << "Average number of photons created
265 << fWLS2EmissionCount / TotNbofEvents
266 if (fWLS2EmissionCount > 0) {
267 G4cout << " Average energy per photon: " <
268 << " eV." << G4endl;
269 }
270 G4cout << "Average energy of WLS2 photons cr
271 << (fWLS2EmissionEnergy / eV) / TotNb
272
273 G4cout << "Average number of OpRayleigh per
274 << G4endl;
275 G4cout << "Average number of OpAbsorption pe
276 G4cout << "\nSurface events (on +X surface,
277 G4cout << "# of primary particles: " <<
278 G4cout << "OpAbsorption before surface: " <<
279 G4cout << "Total # of surface events: " <<
280 if (fParticle->GetParticleName() == "optical
281 G4cout << "Unaccounted for: "
282 << fTotalSurface + fOpAbsorptionPri
283 }
284 G4cout << "\nSurface events by process:" <<
285 if (fBoundaryProcs[Transmission] > 0) {
286 G4cout << " Transmission: "
287 << G4endl;
288 }
289 if (fBoundaryProcs[FresnelRefraction] > 0) {
290 G4cout << " Fresnel refraction: "
291 << G4endl;
292 }
293 if (fBoundaryProcs[FresnelReflection] > 0) {
294 G4cout << " Fresnel reflection: "
295 << G4endl;
296 }
297 if (fBoundaryProcs[TotalInternalReflection]
298 G4cout << " Total internal reflection: "
299 << fBoundaryProcs[TotalInternalRefl
300 }
301 if (fBoundaryProcs[LambertianReflection] > 0
302 G4cout << " Lambertian reflection: "
303 << fBoundaryProcs[LambertianReflect
304 }
305 if (fBoundaryProcs[LobeReflection] > 0) {
306 G4cout << " Lobe reflection: "
307 << G4endl;
308 }
309 if (fBoundaryProcs[SpikeReflection] > 0) {
310 G4cout << " Spike reflection: "
311 << G4endl;
312 }
313 if (fBoundaryProcs[BackScattering] > 0) {
314 G4cout << " Backscattering: "
315 << G4endl;
316 }
317 if (fBoundaryProcs[Absorption] > 0) {
318 G4cout << " Absorption: "
319 << G4endl;
320 }
321 if (fBoundaryProcs[Detection] > 0) {
322 G4cout << " Detection: "
323 << G4endl;
324 }
325 if (fBoundaryProcs[NotAtBoundary] > 0) {
326 G4cout << " Not at boundary: "
327 << G4endl;
328 }
329 if (fBoundaryProcs[SameMaterial] > 0) {
330 G4cout << " Same material: "
331 << G4endl;
332 }
333 if (fBoundaryProcs[StepTooSmall] > 0) {
334 G4cout << " Step too small: "
335 << G4endl;
336 }
337 if (fBoundaryProcs[NoRINDEX] > 0) {
338 G4cout << " No RINDEX: "
339 }
340 // LBNL polished
341 if (fBoundaryProcs[PolishedLumirrorAirReflec
342 G4cout << " Polished Lumirror Air reflect
343 << fBoundaryProcs[PolishedLumirrorA
344 }
345 if (fBoundaryProcs[PolishedLumirrorGlueRefle
346 G4cout << " Polished Lumirror Glue reflec
347 << fBoundaryProcs[PolishedLumirrorG
348 }
349 if (fBoundaryProcs[PolishedAirReflection] >
350 G4cout << " Polished Air reflection: " <<
351 << G4endl;
352 }
353 if (fBoundaryProcs[PolishedTeflonAirReflecti
354 G4cout << " Polished Teflon Air reflectio
355 << fBoundaryProcs[PolishedTeflonAir
356 }
357 if (fBoundaryProcs[PolishedTiOAirReflection]
358 G4cout << " Polished TiO Air reflection:
359 << fBoundaryProcs[PolishedTiOAirRef
360 }
361 if (fBoundaryProcs[PolishedTyvekAirReflectio
362 G4cout << " Polished Tyvek Air reflection
363 << fBoundaryProcs[PolishedTyvekAirR
364 }
365 if (fBoundaryProcs[PolishedVM2000AirReflecti
366 G4cout << " Polished VM2000 Air reflectio
367 << fBoundaryProcs[PolishedVM2000Air
368 }
369 if (fBoundaryProcs[PolishedVM2000GlueReflect
370 G4cout << " Polished VM2000 Glue reflecti
371 << fBoundaryProcs[PolishedVM2000Glu
372 }
373 // LBNL etched
374 if (fBoundaryProcs[EtchedLumirrorAirReflecti
375 G4cout << " Etched Lumirror Air reflectio
376 << fBoundaryProcs[EtchedLumirrorAir
377 }
378 if (fBoundaryProcs[EtchedLumirrorGlueReflect
379 G4cout << " Etched Lumirror Glue reflecti
380 << fBoundaryProcs[EtchedLumirrorGlu
381 }
382 if (fBoundaryProcs[EtchedAirReflection] > 0)
383 G4cout << " Etched Air reflection: " << s
384 << G4endl;
385 }
386 if (fBoundaryProcs[EtchedTeflonAirReflection
387 G4cout << " Etched Teflon Air reflection:
388 << fBoundaryProcs[EtchedTeflonAirRe
389 }
390 if (fBoundaryProcs[EtchedTiOAirReflection] >
391 G4cout << " Etched TiO Air reflection: "
392 << fBoundaryProcs[EtchedTiOAirRefle
393 }
394 if (fBoundaryProcs[EtchedTyvekAirReflection]
395 G4cout << " Etched Tyvek Air reflection:
396 << fBoundaryProcs[EtchedTyvekAirRef
397 }
398 if (fBoundaryProcs[EtchedVM2000AirReflection
399 G4cout << " Etched VM2000 Air reflection:
400 << fBoundaryProcs[EtchedVM2000AirRe
401 }
402 if (fBoundaryProcs[EtchedVM2000GlueReflectio
403 G4cout << " Etched VM2000 Glue reflection
404 << fBoundaryProcs[EtchedVM2000GlueR
405 }
406 // LBNL ground
407 if (fBoundaryProcs[GroundLumirrorAirReflecti
408 G4cout << " Ground Lumirror Air reflectio
409 << fBoundaryProcs[GroundLumirrorAir
410 }
411 if (fBoundaryProcs[GroundLumirrorGlueReflect
412 G4cout << " Ground Lumirror Glue reflecti
413 << fBoundaryProcs[GroundLumirrorGlu
414 }
415 if (fBoundaryProcs[GroundAirReflection] > 0)
416 G4cout << " Ground Air reflection: " << s
417 << G4endl;
418 }
419 if (fBoundaryProcs[GroundTeflonAirReflection
420 G4cout << " Ground Teflon Air reflection:
421 << fBoundaryProcs[GroundTeflonAirRe
422 }
423 if (fBoundaryProcs[GroundTiOAirReflection] >
424 G4cout << " Ground TiO Air reflection: "
425 << fBoundaryProcs[GroundTiOAirRefle
426 }
427 if (fBoundaryProcs[GroundTyvekAirReflection]
428 G4cout << " Ground Tyvek Air reflection:
429 << fBoundaryProcs[GroundTyvekAirRef
430 }
431 if (fBoundaryProcs[GroundVM2000AirReflection
432 G4cout << " Ground VM2000 Air reflection:
433 << fBoundaryProcs[GroundVM2000AirRe
434 }
435 if (fBoundaryProcs[GroundVM2000GlueReflectio
436 G4cout << " Ground VM2000 Glue reflection
437 << fBoundaryProcs[GroundVM2000GlueR
438 }
439 if (fBoundaryProcs[CoatedDielectricRefractio
440 G4cout << " CoatedDielectricRefraction: "
441 << fBoundaryProcs[CoatedDielectricR
442 }
443 if (fBoundaryProcs[CoatedDielectricReflectio
444 G4cout << " CoatedDielectricReflection: "
445 << fBoundaryProcs[CoatedDielectricR
446 }
447 if (fBoundaryProcs[CoatedDielectricFrustrate
448 G4cout << " CoatedDielectricFrustratedTra
449 << fBoundaryProcs[CoatedDielectricF
450 }
451
452 G4int sum = std::accumulate(fBoundaryProcs.b
453 G4cout << " Sum: " <<
454 G4cout << " Unaccounted for: " <<
455
456 G4cout << "---------------------------------
457 G4cout.setf(mode, std::ios::floatfield);
458 G4cout.precision(prec);
459
460 G4int histo_id_refract = analysisMan->GetH1I
461 G4int histo_id_reflect = analysisMan->GetH1I
462 G4int histo_id_spike = analysisMan->GetH1Id(
463 G4int histo_id_absorption = analysisMan->Get
464
465 if (analysisMan->GetH1Activation(histo_id_re
466 && analysisMan->GetH1Activation(histo_id
467 {
468 G4double rindex1 =
469 det->GetTankMaterial()->GetMaterialPrope
470 G4double rindex2 =
471 det->GetWorldMaterial()->GetMaterialProp
472
473 auto histo_refract = analysisMan->GetH1(hi
474 auto histo_reflect = analysisMan->GetH1(hi
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
480 // refract.push_back(histo_refract->bin_
481 reflect.push_back(histo_reflect->bin_hei
482 // tir.push_back(histo_TIR->bin_height(i
483 tot.push_back(histo_refract->bin_height(
484 }
485
486 // find Brewster angle: Rp = 0
487 // need enough statistics for this method
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_
495 bin_width =
496 histo_reflect->axis().bin_upper_edge
497 min_angle += bin_width / 2.;
498 }
499 }
500 G4cout << "Polarization of primary optical
501 << G4endl;
502 if (fPolarized && fPolarization == 0.0) {
503 G4cout << "Reflectance shows a minimum a
504 G4cout << " deg. Expected Brewster angle
505 << (360. / CLHEP::twopi) * std::a
506 }
507
508 // find angle of total internal reflection
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
513 if (histo_refract->bin_height(i) > 0. &&
514 min_angle = histo_refract->axis().bin_
515 bin_width =
516 histo_reflect->axis().bin_upper_edge
517 min_angle += bin_width / 2.;
518 break;
519 }
520 }
521 if (fPolarized) {
522 G4cout << "Fresnel transmission goes to
523 << " deg."
524 << " Expected: " << (360. / CLHEP
525 << G4endl;
526 }
527
528 // Normalize the transmission/reflection h
529 // Only if x values are the same
530 if ((analysisMan->GetH1Nbins(histo_id_refr
531 && (analysisMan->GetH1Xmin(histo_id_re
532 && (analysisMan->GetH1Xmax(histo_id_re
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
540 // "bin+1" below because bin 0 is unde
541 // NB. We are ignoring underflow/overf
542 histo_refract->get_bin_content(bin + 1
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 +
550 }
551
552 histo_reflect->get_bin_content(bin + 1
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 +
560 }
561
562 G4int histo_id_fresnelrefl = analysisM
563 auto histo_fresnelreflect = analysisMa
564 histo_fresnelreflect->get_bin_content(
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_conten
572 }
573
574 G4int histo_id_TIR = analysisMan->GetH
575 auto histo_TIR = analysisMan->GetH1(hi
576 if (analysisMan->GetH1Activation(histo
577 histo_TIR->get_bin_content(bin + 1,
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
585 }
586 }
587 }
588 }
589 else {
590 G4cout << "Not going to normalize refrac
591 << "histograms because bins are n
592 }
593 }
594
595 // complex index of refraction; have spike r
596 // Only works for polished surfaces. Ground
597 else if (analysisMan->GetH1Activation(histo_
598 && analysisMan->GetH1Activation(his
599 {
600 auto histo_spike = analysisMan->GetH1(hist
601 auto histo_absorption = analysisMan->GetH1
602
603 std::vector<G4double> tot;
604 for (size_t i = 0; i < histo_absorption->a
605 tot.push_back(histo_absorption->bin_heig
606 }
607
608 if ((analysisMan->GetH1Nbins(histo_id_abso
609 && (analysisMan->GetH1Xmin(histo_id_ab
610 && (analysisMan->GetH1Xmax(histo_id_ab
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_absorpt
618 // "bin+1" below because bin 0 is unde
619 // NB. We are ignoring underflow/overf
620 histo_absorption->get_bin_content(bin
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(bi
628 }
629
630 histo_spike->get_bin_content(bin + 1,
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
638 }
639 }
640 }
641 else {
642 G4cout << "Not going to normalize spike
643 << "histograms because bins are n
644 }
645 }
646 }
647