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