| 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 {
>> 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
47 fBoundaryProcs.assign(43, 0); 73 fBoundaryProcs.assign(43, 0);
48 } 74 }
49 75
50 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 void Run::SetPrimary(G4ParticleDefinition* par << 77 Run::~Run() {}
52 G4double polarization) << 78
>> 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 80 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy,
>> 81 G4bool polarized, G4double polarization)
53 { 82 {
54 fParticle = particle; << 83 fParticle = particle;
55 fEkin = energy; << 84 fEkin = energy;
56 fPolarized = polarized; << 85 fPolarized = polarized;
57 fPolarization = polarization; 86 fPolarization = polarization;
58 } 87 }
59 88
60 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 void Run::Merge(const G4Run* run) 90 void Run::Merge(const G4Run* run)
62 { 91 {
63 const Run* localRun = static_cast<const Run* 92 const Run* localRun = static_cast<const Run*>(run);
64 93
65 // pass information about primary particle 94 // pass information about primary particle
66 fParticle = localRun->fParticle; << 95 fParticle = localRun->fParticle;
67 fEkin = localRun->fEkin; << 96 fEkin = localRun->fEkin;
68 fPolarized = localRun->fPolarized; << 97 fPolarized = localRun->fPolarized;
69 fPolarization = localRun->fPolarization; 98 fPolarization = localRun->fPolarization;
70 99
71 fCerenkovEnergy += localRun->fCerenkovEnergy 100 fCerenkovEnergy += localRun->fCerenkovEnergy;
72 fScintEnergy += localRun->fScintEnergy; 101 fScintEnergy += localRun->fScintEnergy;
73 fWLSAbsorptionEnergy += localRun->fWLSAbsorp 102 fWLSAbsorptionEnergy += localRun->fWLSAbsorptionEnergy;
74 fWLSEmissionEnergy += localRun->fWLSEmission 103 fWLSEmissionEnergy += localRun->fWLSEmissionEnergy;
75 fWLS2AbsorptionEnergy += localRun->fWLS2Abso 104 fWLS2AbsorptionEnergy += localRun->fWLS2AbsorptionEnergy;
76 fWLS2EmissionEnergy += localRun->fWLS2Emissi 105 fWLS2EmissionEnergy += localRun->fWLS2EmissionEnergy;
77 106
78 fCerenkovCount += localRun->fCerenkovCount; 107 fCerenkovCount += localRun->fCerenkovCount;
79 fScintCount += localRun->fScintCount; 108 fScintCount += localRun->fScintCount;
80 fWLSAbsorptionCount += localRun->fWLSAbsorpt 109 fWLSAbsorptionCount += localRun->fWLSAbsorptionCount;
81 fWLSEmissionCount += localRun->fWLSEmissionC 110 fWLSEmissionCount += localRun->fWLSEmissionCount;
82 fWLS2AbsorptionCount += localRun->fWLS2Absor 111 fWLS2AbsorptionCount += localRun->fWLS2AbsorptionCount;
83 fWLS2EmissionCount += localRun->fWLS2Emissio 112 fWLS2EmissionCount += localRun->fWLS2EmissionCount;
84 fRayleighCount += localRun->fRayleighCount; 113 fRayleighCount += localRun->fRayleighCount;
85 114
86 fTotalSurface += localRun->fTotalSurface; 115 fTotalSurface += localRun->fTotalSurface;
87 116
88 fOpAbsorption += localRun->fOpAbsorption; 117 fOpAbsorption += localRun->fOpAbsorption;
89 fOpAbsorptionPrior += localRun->fOpAbsorptio 118 fOpAbsorptionPrior += localRun->fOpAbsorptionPrior;
90 119
91 for (size_t i = 0; i < fBoundaryProcs.size() << 120 for(size_t i = 0; i < fBoundaryProcs.size(); ++i)
>> 121 {
92 fBoundaryProcs[i] += localRun->fBoundaryPr 122 fBoundaryProcs[i] += localRun->fBoundaryProcs[i];
93 } 123 }
94 124
95 G4Run::Merge(run); 125 G4Run::Merge(run);
96 } 126 }
97 127
98 //....oooOO0OOooo........oooOO0OOooo........oo 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 void Run::EndOfRun() 129 void Run::EndOfRun()
100 { 130 {
101 if (numberOfEvent == 0) return; << 131 if(numberOfEvent == 0)
102 auto TotNbofEvents = (G4double)numberOfEvent << 132 return;
>> 133 G4double TotNbofEvents = (G4double) numberOfEvent;
103 134
104 G4AnalysisManager* analysisMan = G4AnalysisM 135 G4AnalysisManager* analysisMan = G4AnalysisManager::Instance();
105 G4int id = analysisMan->GetH1Id("Cerenkov sp << 136 G4int id = analysisMan->GetH1Id("Cerenkov spectrum");
106 analysisMan->SetH1XAxisTitle(id, "Energy [eV 137 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
107 analysisMan->SetH1YAxisTitle(id, "Number of 138 analysisMan->SetH1YAxisTitle(id, "Number of photons");
108 139
109 id = analysisMan->GetH1Id("Scintillation spe 140 id = analysisMan->GetH1Id("Scintillation spectrum");
110 analysisMan->SetH1XAxisTitle(id, "Energy [eV 141 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
111 analysisMan->SetH1YAxisTitle(id, "Number of 142 analysisMan->SetH1YAxisTitle(id, "Number of photons");
112 143
113 id = analysisMan->GetH1Id("Scintillation tim 144 id = analysisMan->GetH1Id("Scintillation time");
114 analysisMan->SetH1XAxisTitle(id, "Creation t 145 analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
115 analysisMan->SetH1YAxisTitle(id, "Number of 146 analysisMan->SetH1YAxisTitle(id, "Number of photons");
116 147
117 id = analysisMan->GetH1Id("WLS abs"); 148 id = analysisMan->GetH1Id("WLS abs");
118 analysisMan->SetH1XAxisTitle(id, "Energy [eV 149 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
119 analysisMan->SetH1YAxisTitle(id, "Number of 150 analysisMan->SetH1YAxisTitle(id, "Number of photons");
120 151
121 id = analysisMan->GetH1Id("WLS em"); 152 id = analysisMan->GetH1Id("WLS em");
122 analysisMan->SetH1XAxisTitle(id, "Energy [eV 153 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
123 analysisMan->SetH1YAxisTitle(id, "Number of 154 analysisMan->SetH1YAxisTitle(id, "Number of photons");
124 155
125 id = analysisMan->GetH1Id("WLS time"); 156 id = analysisMan->GetH1Id("WLS time");
126 analysisMan->SetH1XAxisTitle(id, "Creation t 157 analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
127 analysisMan->SetH1YAxisTitle(id, "Number of 158 analysisMan->SetH1YAxisTitle(id, "Number of photons");
128 159
129 id = analysisMan->GetH1Id("WLS2 abs"); 160 id = analysisMan->GetH1Id("WLS2 abs");
130 analysisMan->SetH1XAxisTitle(id, "Energy [eV 161 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
131 analysisMan->SetH1YAxisTitle(id, "Number of 162 analysisMan->SetH1YAxisTitle(id, "Number of photons");
132 163
133 id = analysisMan->GetH1Id("WLS2 em"); 164 id = analysisMan->GetH1Id("WLS2 em");
134 analysisMan->SetH1XAxisTitle(id, "Energy [eV 165 analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
135 analysisMan->SetH1YAxisTitle(id, "Number of 166 analysisMan->SetH1YAxisTitle(id, "Number of photons");
136 167
137 id = analysisMan->GetH1Id("WLS2 time"); 168 id = analysisMan->GetH1Id("WLS2 time");
138 analysisMan->SetH1XAxisTitle(id, "Creation t 169 analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
139 analysisMan->SetH1YAxisTitle(id, "Number of 170 analysisMan->SetH1YAxisTitle(id, "Number of photons");
140 171
141 id = analysisMan->GetH1Id("bdry status"); 172 id = analysisMan->GetH1Id("bdry status");
142 analysisMan->SetH1XAxisTitle(id, "Status cod 173 analysisMan->SetH1XAxisTitle(id, "Status code");
143 analysisMan->SetH1YAxisTitle(id, "Number of 174 analysisMan->SetH1YAxisTitle(id, "Number of photons");
144 175
145 id = analysisMan->GetH1Id("x_backward"); 176 id = analysisMan->GetH1Id("x_backward");
146 analysisMan->SetH1XAxisTitle(id, "Direction 177 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
147 analysisMan->SetH1YAxisTitle(id, "Number of 178 analysisMan->SetH1YAxisTitle(id, "Number of photons");
148 179
149 id = analysisMan->GetH1Id("y_backward"); 180 id = analysisMan->GetH1Id("y_backward");
150 analysisMan->SetH1XAxisTitle(id, "Direction 181 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
151 analysisMan->SetH1YAxisTitle(id, "Number of 182 analysisMan->SetH1YAxisTitle(id, "Number of photons");
152 183
153 id = analysisMan->GetH1Id("z_backward"); 184 id = analysisMan->GetH1Id("z_backward");
154 analysisMan->SetH1XAxisTitle(id, "Direction 185 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
155 analysisMan->SetH1YAxisTitle(id, "Number of 186 analysisMan->SetH1YAxisTitle(id, "Number of photons");
156 187
157 id = analysisMan->GetH1Id("x_forward"); 188 id = analysisMan->GetH1Id("x_forward");
158 analysisMan->SetH1XAxisTitle(id, "Direction 189 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
159 analysisMan->SetH1YAxisTitle(id, "Number of 190 analysisMan->SetH1YAxisTitle(id, "Number of photons");
160 191
161 id = analysisMan->GetH1Id("y_forward"); 192 id = analysisMan->GetH1Id("y_forward");
162 analysisMan->SetH1XAxisTitle(id, "Direction 193 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
163 analysisMan->SetH1YAxisTitle(id, "Number of 194 analysisMan->SetH1YAxisTitle(id, "Number of photons");
164 195
165 id = analysisMan->GetH1Id("z_forward"); 196 id = analysisMan->GetH1Id("z_forward");
166 analysisMan->SetH1XAxisTitle(id, "Direction 197 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
167 analysisMan->SetH1YAxisTitle(id, "Number of 198 analysisMan->SetH1YAxisTitle(id, "Number of photons");
168 199
169 id = analysisMan->GetH1Id("x_fresnel"); 200 id = analysisMan->GetH1Id("x_fresnel");
170 analysisMan->SetH1XAxisTitle(id, "Direction 201 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
171 analysisMan->SetH1YAxisTitle(id, "Number of 202 analysisMan->SetH1YAxisTitle(id, "Number of photons");
172 203
173 id = analysisMan->GetH1Id("y_fresnel"); 204 id = analysisMan->GetH1Id("y_fresnel");
174 analysisMan->SetH1XAxisTitle(id, "Direction 205 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
175 analysisMan->SetH1YAxisTitle(id, "Number of 206 analysisMan->SetH1YAxisTitle(id, "Number of photons");
176 207
177 id = analysisMan->GetH1Id("z_fresnel"); 208 id = analysisMan->GetH1Id("z_fresnel");
178 analysisMan->SetH1XAxisTitle(id, "Direction 209 analysisMan->SetH1XAxisTitle(id, "Direction cosine");
179 analysisMan->SetH1YAxisTitle(id, "Number of 210 analysisMan->SetH1YAxisTitle(id, "Number of photons");
180 211
181 id = analysisMan->GetH1Id("Fresnel reflectio 212 id = analysisMan->GetH1Id("Fresnel reflection");
182 analysisMan->SetH1XAxisTitle(id, "Angle [deg 213 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
183 analysisMan->SetH1YAxisTitle(id, "Fraction o 214 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
184 215
185 id = analysisMan->GetH1Id("Fresnel refractio 216 id = analysisMan->GetH1Id("Fresnel refraction");
186 analysisMan->SetH1XAxisTitle(id, "Angle [deg 217 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
187 analysisMan->SetH1YAxisTitle(id, "Fraction o 218 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
188 219
189 id = analysisMan->GetH1Id("Total internal re 220 id = analysisMan->GetH1Id("Total internal reflection");
190 analysisMan->SetH1XAxisTitle(id, "Angle [deg 221 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
191 analysisMan->SetH1YAxisTitle(id, "Fraction o 222 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
192 223
193 id = analysisMan->GetH1Id("Fresnel reflectio 224 id = analysisMan->GetH1Id("Fresnel reflection plus TIR");
194 analysisMan->SetH1XAxisTitle(id, "Angle [deg 225 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
195 analysisMan->SetH1YAxisTitle(id, "Fraction o 226 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
196 227
197 id = analysisMan->GetH1Id("Absorption"); 228 id = analysisMan->GetH1Id("Absorption");
198 analysisMan->SetH1XAxisTitle(id, "Angle [deg 229 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
199 analysisMan->SetH1YAxisTitle(id, "Fraction o 230 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
200 231
201 id = analysisMan->GetH1Id("Transmitted"); 232 id = analysisMan->GetH1Id("Transmitted");
202 analysisMan->SetH1XAxisTitle(id, "Angle [deg 233 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
203 analysisMan->SetH1YAxisTitle(id, "Fraction o 234 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
204 235
205 id = analysisMan->GetH1Id("Spike reflection" 236 id = analysisMan->GetH1Id("Spike reflection");
206 analysisMan->SetH1XAxisTitle(id, "Angle [deg 237 analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
207 analysisMan->SetH1YAxisTitle(id, "Fraction o 238 analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
208 239
209 const auto det = << 240 const DetectorConstruction* det =
210 (const DetectorConstruction*)(G4RunManager << 241 (const DetectorConstruction*) (G4RunManager::GetRunManager()
>> 242 ->GetUserDetectorConstruction());
211 243
212 std::ios::fmtflags mode = G4cout.flags(); 244 std::ios::fmtflags mode = G4cout.flags();
213 G4int prec = G4cout.precision(2); << 245 G4int prec = G4cout.precision(2);
214 246
215 G4cout << "\n Run Summary\n"; 247 G4cout << "\n Run Summary\n";
216 G4cout << "--------------------------------- 248 G4cout << "---------------------------------\n";
217 G4cout << "Primary particle was: " << fParti << 249 G4cout << "Primary particle was: " << fParticle->GetParticleName()
218 << G4BestUnit(fEkin, "Energy") << "." << 250 << " with energy " << G4BestUnit(fEkin, "Energy") << "." << G4endl;
219 G4cout << "Number of events: " << numberOfEv 251 G4cout << "Number of events: " << numberOfEvent << G4endl;
220 252
221 G4cout << "Material of world: " << det->GetW << 253 G4cout << "Material of world: " << det->GetWorldMaterial()->GetName()
222 G4cout << "Material of tank: " << det->GetT << 254 << G4endl;
>> 255 G4cout << "Material of tank: " << det->GetTankMaterial()->GetName() << G4endl
>> 256 << G4endl;
223 257
224 if (fParticle->GetParticleName() != "optical << 258 if(fParticle->GetParticleName() != "opticalphoton")
>> 259 {
225 G4cout << "Average energy of Cerenkov phot 260 G4cout << "Average energy of Cerenkov photons created per event: "
226 << (fCerenkovEnergy / eV) / TotNbof 261 << (fCerenkovEnergy / eV) / TotNbofEvents << " eV." << G4endl;
227 G4cout << "Average number of Cerenkov phot 262 G4cout << "Average number of Cerenkov photons created per event: "
228 << fCerenkovCount / TotNbofEvents < 263 << fCerenkovCount / TotNbofEvents << G4endl;
229 if (fCerenkovCount > 0) { << 264 if(fCerenkovCount > 0)
230 G4cout << " Average energy per photon: " << 265 {
231 << G4endl; << 266 G4cout << " Average energy per photon: "
>> 267 << (fCerenkovEnergy / eV) / fCerenkovCount << " eV." << G4endl;
232 } 268 }
233 G4cout << "Average energy of scintillation 269 G4cout << "Average energy of scintillation photons created per event: "
234 << (fScintEnergy / eV) / TotNbofEve 270 << (fScintEnergy / eV) / TotNbofEvents << " eV." << G4endl;
235 G4cout << "Average number of scintillation 271 G4cout << "Average number of scintillation photons created per event: "
236 << fScintCount / TotNbofEvents << G 272 << fScintCount / TotNbofEvents << G4endl;
237 if (fScintCount > 0) { << 273 if(fScintCount > 0)
238 G4cout << " Average energy per photon: " << 274 {
239 << G4endl; << 275 G4cout << " Average energy per photon: "
>> 276 << (fScintEnergy / eV) / fScintCount << " eV." << G4endl;
240 } 277 }
241 } 278 }
242 279
243 G4cout << "Average number of photons absorbe 280 G4cout << "Average number of photons absorbed by WLS per event: "
244 << fWLSAbsorptionCount / G4double(Tot 281 << fWLSAbsorptionCount / G4double(TotNbofEvents) << " " << G4endl;
245 if (fWLSAbsorptionCount > 0) { << 282 if(fWLSAbsorptionCount > 0)
246 G4cout << " Average energy per photon: " < << 283 {
247 << " eV." << G4endl; << 284 G4cout << " Average energy per photon: "
>> 285 << (fWLSAbsorptionEnergy / eV) / fWLSAbsorptionCount << " eV."
>> 286 << G4endl;
248 } 287 }
249 G4cout << "Average number of photons created 288 G4cout << "Average number of photons created by WLS per event: "
250 << fWLSEmissionCount / TotNbofEvents 289 << fWLSEmissionCount / TotNbofEvents << G4endl;
251 if (fWLSEmissionCount > 0) { << 290 if(fWLSEmissionCount > 0)
252 G4cout << " Average energy per photon: " < << 291 {
253 << " eV." << G4endl; << 292 G4cout << " Average energy per photon: "
>> 293 << (fWLSEmissionEnergy / eV) / fWLSEmissionCount << " eV." << G4endl;
254 } 294 }
255 G4cout << "Average energy of WLS photons cre 295 G4cout << "Average energy of WLS photons created per event: "
256 << (fWLSEmissionEnergy / eV) / TotNbo 296 << (fWLSEmissionEnergy / eV) / TotNbofEvents << " eV." << G4endl;
257 297
258 G4cout << "Average number of photons absorbe 298 G4cout << "Average number of photons absorbed by WLS2 per event: "
259 << fWLS2AbsorptionCount / G4double(To 299 << fWLS2AbsorptionCount / G4double(TotNbofEvents) << " " << G4endl;
260 if (fWLS2AbsorptionCount > 0) { << 300 if(fWLS2AbsorptionCount > 0)
261 G4cout << " Average energy per photon: " < << 301 {
262 << " eV." << G4endl; << 302 G4cout << " Average energy per photon: "
>> 303 << (fWLS2AbsorptionEnergy / eV) / fWLS2AbsorptionCount << " eV."
>> 304 << G4endl;
263 } 305 }
264 G4cout << "Average number of photons created 306 G4cout << "Average number of photons created by WLS2 per event: "
265 << fWLS2EmissionCount / TotNbofEvents 307 << fWLS2EmissionCount / TotNbofEvents << G4endl;
266 if (fWLS2EmissionCount > 0) { << 308 if(fWLS2EmissionCount > 0)
267 G4cout << " Average energy per photon: " < << 309 {
268 << " eV." << G4endl; << 310 G4cout << " Average energy per photon: "
>> 311 << (fWLS2EmissionEnergy / eV) / fWLS2EmissionCount << " eV."
>> 312 << G4endl;
269 } 313 }
270 G4cout << "Average energy of WLS2 photons cr 314 G4cout << "Average energy of WLS2 photons created per event: "
271 << (fWLS2EmissionEnergy / eV) / TotNb 315 << (fWLS2EmissionEnergy / eV) / TotNbofEvents << " eV." << G4endl;
272 316
273 G4cout << "Average number of OpRayleigh per << 317 G4cout << "Average number of OpRayleigh per event: "
>> 318 << fRayleighCount / TotNbofEvents << G4endl;
>> 319 G4cout << "Average number of OpAbsorption per event: "
>> 320 << fOpAbsorption / TotNbofEvents << G4endl;
>> 321 G4cout << "\nSurface events (on +X surface, maximum one per photon) this run:"
>> 322 << G4endl;
>> 323 G4cout << "# of primary particles: " << std::setw(8) << TotNbofEvents
>> 324 << G4endl;
>> 325 G4cout << "OpAbsorption before surface: " << std::setw(8)
>> 326 << fOpAbsorptionPrior << G4endl;
>> 327 G4cout << "Total # of surface events: " << std::setw(8) << fTotalSurface
274 << G4endl; 328 << G4endl;
275 G4cout << "Average number of OpAbsorption pe << 329 if(fParticle->GetParticleName() == "opticalphoton")
276 G4cout << "\nSurface events (on +X surface, << 330 {
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: " 331 G4cout << "Unaccounted for: " << std::setw(8)
282 << fTotalSurface + fOpAbsorptionPri 332 << fTotalSurface + fOpAbsorptionPrior - TotNbofEvents << G4endl;
283 } 333 }
284 G4cout << "\nSurface events by process:" << 334 G4cout << "\nSurface events by process:" << G4endl;
285 if (fBoundaryProcs[Transmission] > 0) { << 335 if(fBoundaryProcs[Transmission] > 0)
286 G4cout << " Transmission: " << 336 {
287 << G4endl; << 337 G4cout << " Transmission: " << std::setw(8)
>> 338 << fBoundaryProcs[Transmission] << G4endl;
288 } 339 }
289 if (fBoundaryProcs[FresnelRefraction] > 0) { << 340 if(fBoundaryProcs[FresnelRefraction] > 0)
290 G4cout << " Fresnel refraction: " << 341 {
291 << G4endl; << 342 G4cout << " Fresnel refraction: " << std::setw(8)
>> 343 << fBoundaryProcs[FresnelRefraction] << G4endl;
292 } 344 }
293 if (fBoundaryProcs[FresnelReflection] > 0) { << 345 if(fBoundaryProcs[FresnelReflection] > 0)
294 G4cout << " Fresnel reflection: " << 346 {
295 << G4endl; << 347 G4cout << " Fresnel reflection: " << std::setw(8)
>> 348 << fBoundaryProcs[FresnelReflection] << G4endl;
296 } 349 }
297 if (fBoundaryProcs[TotalInternalReflection] << 350 if(fBoundaryProcs[TotalInternalReflection] > 0)
>> 351 {
298 G4cout << " Total internal reflection: " 352 G4cout << " Total internal reflection: " << std::setw(8)
299 << fBoundaryProcs[TotalInternalRefl 353 << fBoundaryProcs[TotalInternalReflection] << G4endl;
300 } 354 }
301 if (fBoundaryProcs[LambertianReflection] > 0 << 355 if(fBoundaryProcs[LambertianReflection] > 0)
>> 356 {
302 G4cout << " Lambertian reflection: " 357 G4cout << " Lambertian reflection: " << std::setw(8)
303 << fBoundaryProcs[LambertianReflect 358 << fBoundaryProcs[LambertianReflection] << G4endl;
304 } 359 }
305 if (fBoundaryProcs[LobeReflection] > 0) { << 360 if(fBoundaryProcs[LobeReflection] > 0)
306 G4cout << " Lobe reflection: " << 361 {
307 << G4endl; << 362 G4cout << " Lobe reflection: " << std::setw(8)
>> 363 << fBoundaryProcs[LobeReflection] << G4endl;
308 } 364 }
309 if (fBoundaryProcs[SpikeReflection] > 0) { << 365 if(fBoundaryProcs[SpikeReflection] > 0)
310 G4cout << " Spike reflection: " << 366 {
311 << G4endl; << 367 G4cout << " Spike reflection: " << std::setw(8)
>> 368 << fBoundaryProcs[SpikeReflection] << G4endl;
312 } 369 }
313 if (fBoundaryProcs[BackScattering] > 0) { << 370 if(fBoundaryProcs[BackScattering] > 0)
314 G4cout << " Backscattering: " << 371 {
315 << G4endl; << 372 G4cout << " Backscattering: " << std::setw(8)
>> 373 << fBoundaryProcs[BackScattering] << G4endl;
316 } 374 }
317 if (fBoundaryProcs[Absorption] > 0) { << 375 if(fBoundaryProcs[Absorption] > 0)
318 G4cout << " Absorption: " << 376 {
319 << G4endl; << 377 G4cout << " Absorption: " << std::setw(8)
>> 378 << fBoundaryProcs[Absorption] << G4endl;
320 } 379 }
321 if (fBoundaryProcs[Detection] > 0) { << 380 if(fBoundaryProcs[Detection] > 0)
322 G4cout << " Detection: " << 381 {
323 << G4endl; << 382 G4cout << " Detection: " << std::setw(8)
>> 383 << fBoundaryProcs[Detection] << G4endl;
324 } 384 }
325 if (fBoundaryProcs[NotAtBoundary] > 0) { << 385 if(fBoundaryProcs[NotAtBoundary] > 0)
326 G4cout << " Not at boundary: " << 386 {
327 << G4endl; << 387 G4cout << " Not at boundary: " << std::setw(8)
>> 388 << fBoundaryProcs[NotAtBoundary] << G4endl;
328 } 389 }
329 if (fBoundaryProcs[SameMaterial] > 0) { << 390 if(fBoundaryProcs[SameMaterial] > 0)
330 G4cout << " Same material: " << 391 {
331 << G4endl; << 392 G4cout << " Same material: " << std::setw(8)
>> 393 << fBoundaryProcs[SameMaterial] << G4endl;
332 } 394 }
333 if (fBoundaryProcs[StepTooSmall] > 0) { << 395 if(fBoundaryProcs[StepTooSmall] > 0)
334 G4cout << " Step too small: " << 396 {
335 << G4endl; << 397 G4cout << " Step too small: " << std::setw(8)
>> 398 << fBoundaryProcs[StepTooSmall] << G4endl;
336 } 399 }
337 if (fBoundaryProcs[NoRINDEX] > 0) { << 400 if(fBoundaryProcs[NoRINDEX] > 0)
338 G4cout << " No RINDEX: " << 401 {
>> 402 G4cout << " No RINDEX: " << std::setw(8)
>> 403 << fBoundaryProcs[NoRINDEX] << G4endl;
339 } 404 }
340 // LBNL polished 405 // LBNL polished
341 if (fBoundaryProcs[PolishedLumirrorAirReflec << 406 if(fBoundaryProcs[PolishedLumirrorAirReflection] > 0)
>> 407 {
342 G4cout << " Polished Lumirror Air reflect 408 G4cout << " Polished Lumirror Air reflection: " << std::setw(8)
343 << fBoundaryProcs[PolishedLumirrorA 409 << fBoundaryProcs[PolishedLumirrorAirReflection] << G4endl;
344 } 410 }
345 if (fBoundaryProcs[PolishedLumirrorGlueRefle << 411 if(fBoundaryProcs[PolishedLumirrorGlueReflection] > 0)
>> 412 {
346 G4cout << " Polished Lumirror Glue reflec 413 G4cout << " Polished Lumirror Glue reflection: " << std::setw(8)
347 << fBoundaryProcs[PolishedLumirrorG 414 << fBoundaryProcs[PolishedLumirrorGlueReflection] << G4endl;
348 } 415 }
349 if (fBoundaryProcs[PolishedAirReflection] > << 416 if(fBoundaryProcs[PolishedAirReflection] > 0)
350 G4cout << " Polished Air reflection: " << << 417 {
351 << G4endl; << 418 G4cout << " Polished Air reflection: " << std::setw(8)
>> 419 << fBoundaryProcs[PolishedAirReflection] << G4endl;
352 } 420 }
353 if (fBoundaryProcs[PolishedTeflonAirReflecti << 421 if(fBoundaryProcs[PolishedTeflonAirReflection] > 0)
>> 422 {
354 G4cout << " Polished Teflon Air reflectio 423 G4cout << " Polished Teflon Air reflection: " << std::setw(8)
355 << fBoundaryProcs[PolishedTeflonAir 424 << fBoundaryProcs[PolishedTeflonAirReflection] << G4endl;
356 } 425 }
357 if (fBoundaryProcs[PolishedTiOAirReflection] << 426 if(fBoundaryProcs[PolishedTiOAirReflection] > 0)
>> 427 {
358 G4cout << " Polished TiO Air reflection: 428 G4cout << " Polished TiO Air reflection: " << std::setw(8)
359 << fBoundaryProcs[PolishedTiOAirRef 429 << fBoundaryProcs[PolishedTiOAirReflection] << G4endl;
360 } 430 }
361 if (fBoundaryProcs[PolishedTyvekAirReflectio << 431 if(fBoundaryProcs[PolishedTyvekAirReflection] > 0)
>> 432 {
362 G4cout << " Polished Tyvek Air reflection 433 G4cout << " Polished Tyvek Air reflection: " << std::setw(8)
363 << fBoundaryProcs[PolishedTyvekAirR 434 << fBoundaryProcs[PolishedTyvekAirReflection] << G4endl;
364 } 435 }
365 if (fBoundaryProcs[PolishedVM2000AirReflecti << 436 if(fBoundaryProcs[PolishedVM2000AirReflection] > 0)
>> 437 {
366 G4cout << " Polished VM2000 Air reflectio 438 G4cout << " Polished VM2000 Air reflection: " << std::setw(8)
367 << fBoundaryProcs[PolishedVM2000Air 439 << fBoundaryProcs[PolishedVM2000AirReflection] << G4endl;
368 } 440 }
369 if (fBoundaryProcs[PolishedVM2000GlueReflect << 441 if(fBoundaryProcs[PolishedVM2000GlueReflection] > 0)
>> 442 {
370 G4cout << " Polished VM2000 Glue reflecti 443 G4cout << " Polished VM2000 Glue reflection: " << std::setw(8)
371 << fBoundaryProcs[PolishedVM2000Glu 444 << fBoundaryProcs[PolishedVM2000GlueReflection] << G4endl;
372 } 445 }
373 // LBNL etched 446 // LBNL etched
374 if (fBoundaryProcs[EtchedLumirrorAirReflecti << 447 if(fBoundaryProcs[EtchedLumirrorAirReflection] > 0)
>> 448 {
375 G4cout << " Etched Lumirror Air reflectio 449 G4cout << " Etched Lumirror Air reflection: " << std::setw(8)
376 << fBoundaryProcs[EtchedLumirrorAir 450 << fBoundaryProcs[EtchedLumirrorAirReflection] << G4endl;
377 } 451 }
378 if (fBoundaryProcs[EtchedLumirrorGlueReflect << 452 if(fBoundaryProcs[EtchedLumirrorGlueReflection] > 0)
>> 453 {
379 G4cout << " Etched Lumirror Glue reflecti 454 G4cout << " Etched Lumirror Glue reflection: " << std::setw(8)
380 << fBoundaryProcs[EtchedLumirrorGlu 455 << fBoundaryProcs[EtchedLumirrorGlueReflection] << G4endl;
381 } 456 }
382 if (fBoundaryProcs[EtchedAirReflection] > 0) << 457 if(fBoundaryProcs[EtchedAirReflection] > 0)
383 G4cout << " Etched Air reflection: " << s << 458 {
384 << G4endl; << 459 G4cout << " Etched Air reflection: " << std::setw(8)
>> 460 << fBoundaryProcs[EtchedAirReflection] << G4endl;
385 } 461 }
386 if (fBoundaryProcs[EtchedTeflonAirReflection << 462 if(fBoundaryProcs[EtchedTeflonAirReflection] > 0)
>> 463 {
387 G4cout << " Etched Teflon Air reflection: 464 G4cout << " Etched Teflon Air reflection: " << std::setw(8)
388 << fBoundaryProcs[EtchedTeflonAirRe 465 << fBoundaryProcs[EtchedTeflonAirReflection] << G4endl;
389 } 466 }
390 if (fBoundaryProcs[EtchedTiOAirReflection] > << 467 if(fBoundaryProcs[EtchedTiOAirReflection] > 0)
>> 468 {
391 G4cout << " Etched TiO Air reflection: " 469 G4cout << " Etched TiO Air reflection: " << std::setw(8)
392 << fBoundaryProcs[EtchedTiOAirRefle 470 << fBoundaryProcs[EtchedTiOAirReflection] << G4endl;
393 } 471 }
394 if (fBoundaryProcs[EtchedTyvekAirReflection] << 472 if(fBoundaryProcs[EtchedTyvekAirReflection] > 0)
>> 473 {
395 G4cout << " Etched Tyvek Air reflection: 474 G4cout << " Etched Tyvek Air reflection: " << std::setw(8)
396 << fBoundaryProcs[EtchedTyvekAirRef 475 << fBoundaryProcs[EtchedTyvekAirReflection] << G4endl;
397 } 476 }
398 if (fBoundaryProcs[EtchedVM2000AirReflection << 477 if(fBoundaryProcs[EtchedVM2000AirReflection] > 0)
>> 478 {
399 G4cout << " Etched VM2000 Air reflection: 479 G4cout << " Etched VM2000 Air reflection: " << std::setw(8)
400 << fBoundaryProcs[EtchedVM2000AirRe 480 << fBoundaryProcs[EtchedVM2000AirReflection] << G4endl;
401 } 481 }
402 if (fBoundaryProcs[EtchedVM2000GlueReflectio << 482 if(fBoundaryProcs[EtchedVM2000GlueReflection] > 0)
>> 483 {
403 G4cout << " Etched VM2000 Glue reflection 484 G4cout << " Etched VM2000 Glue reflection: " << std::setw(8)
404 << fBoundaryProcs[EtchedVM2000GlueR 485 << fBoundaryProcs[EtchedVM2000GlueReflection] << G4endl;
405 } 486 }
406 // LBNL ground 487 // LBNL ground
407 if (fBoundaryProcs[GroundLumirrorAirReflecti << 488 if(fBoundaryProcs[GroundLumirrorAirReflection] > 0)
>> 489 {
408 G4cout << " Ground Lumirror Air reflectio 490 G4cout << " Ground Lumirror Air reflection: " << std::setw(8)
409 << fBoundaryProcs[GroundLumirrorAir 491 << fBoundaryProcs[GroundLumirrorAirReflection] << G4endl;
410 } 492 }
411 if (fBoundaryProcs[GroundLumirrorGlueReflect << 493 if(fBoundaryProcs[GroundLumirrorGlueReflection] > 0)
>> 494 {
412 G4cout << " Ground Lumirror Glue reflecti 495 G4cout << " Ground Lumirror Glue reflection: " << std::setw(8)
413 << fBoundaryProcs[GroundLumirrorGlu 496 << fBoundaryProcs[GroundLumirrorGlueReflection] << G4endl;
414 } 497 }
415 if (fBoundaryProcs[GroundAirReflection] > 0) << 498 if(fBoundaryProcs[GroundAirReflection] > 0)
416 G4cout << " Ground Air reflection: " << s << 499 {
417 << G4endl; << 500 G4cout << " Ground Air reflection: " << std::setw(8)
>> 501 << fBoundaryProcs[GroundAirReflection] << G4endl;
418 } 502 }
419 if (fBoundaryProcs[GroundTeflonAirReflection << 503 if(fBoundaryProcs[GroundTeflonAirReflection] > 0)
>> 504 {
420 G4cout << " Ground Teflon Air reflection: 505 G4cout << " Ground Teflon Air reflection: " << std::setw(8)
421 << fBoundaryProcs[GroundTeflonAirRe 506 << fBoundaryProcs[GroundTeflonAirReflection] << G4endl;
422 } 507 }
423 if (fBoundaryProcs[GroundTiOAirReflection] > << 508 if(fBoundaryProcs[GroundTiOAirReflection] > 0)
>> 509 {
424 G4cout << " Ground TiO Air reflection: " 510 G4cout << " Ground TiO Air reflection: " << std::setw(8)
425 << fBoundaryProcs[GroundTiOAirRefle 511 << fBoundaryProcs[GroundTiOAirReflection] << G4endl;
426 } 512 }
427 if (fBoundaryProcs[GroundTyvekAirReflection] << 513 if(fBoundaryProcs[GroundTyvekAirReflection] > 0)
>> 514 {
428 G4cout << " Ground Tyvek Air reflection: 515 G4cout << " Ground Tyvek Air reflection: " << std::setw(8)
429 << fBoundaryProcs[GroundTyvekAirRef 516 << fBoundaryProcs[GroundTyvekAirReflection] << G4endl;
430 } 517 }
431 if (fBoundaryProcs[GroundVM2000AirReflection << 518 if(fBoundaryProcs[GroundVM2000AirReflection] > 0)
>> 519 {
432 G4cout << " Ground VM2000 Air reflection: 520 G4cout << " Ground VM2000 Air reflection: " << std::setw(8)
433 << fBoundaryProcs[GroundVM2000AirRe 521 << fBoundaryProcs[GroundVM2000AirReflection] << G4endl;
434 } 522 }
435 if (fBoundaryProcs[GroundVM2000GlueReflectio << 523 if(fBoundaryProcs[GroundVM2000GlueReflection] > 0)
>> 524 {
436 G4cout << " Ground VM2000 Glue reflection 525 G4cout << " Ground VM2000 Glue reflection: " << std::setw(8)
437 << fBoundaryProcs[GroundVM2000GlueR 526 << fBoundaryProcs[GroundVM2000GlueReflection] << G4endl;
438 } 527 }
439 if (fBoundaryProcs[CoatedDielectricRefractio << 528 if(fBoundaryProcs[CoatedDielectricRefraction] > 0)
>> 529 {
440 G4cout << " CoatedDielectricRefraction: " 530 G4cout << " CoatedDielectricRefraction: " << std::setw(8)
441 << fBoundaryProcs[CoatedDielectricR 531 << fBoundaryProcs[CoatedDielectricRefraction] << G4endl;
442 } 532 }
443 if (fBoundaryProcs[CoatedDielectricReflectio << 533 if(fBoundaryProcs[CoatedDielectricReflection] > 0)
>> 534 {
444 G4cout << " CoatedDielectricReflection: " 535 G4cout << " CoatedDielectricReflection: " << std::setw(8)
445 << fBoundaryProcs[CoatedDielectricR 536 << fBoundaryProcs[CoatedDielectricReflection] << G4endl;
446 } 537 }
447 if (fBoundaryProcs[CoatedDielectricFrustrate << 538 if(fBoundaryProcs[CoatedDielectricFrustratedTransmission] > 0)
>> 539 {
448 G4cout << " CoatedDielectricFrustratedTra 540 G4cout << " CoatedDielectricFrustratedTransmission: " << std::setw(8)
449 << fBoundaryProcs[CoatedDielectricF 541 << fBoundaryProcs[CoatedDielectricFrustratedTransmission] << G4endl;
450 } 542 }
451 543
452 G4int sum = std::accumulate(fBoundaryProcs.b 544 G4int sum = std::accumulate(fBoundaryProcs.begin(), fBoundaryProcs.end(), 0);
453 G4cout << " Sum: " << 545 G4cout << " Sum: " << std::setw(8) << sum << G4endl;
454 G4cout << " Unaccounted for: " << << 546 G4cout << " Unaccounted for: " << std::setw(8)
>> 547 << fTotalSurface - sum << G4endl;
455 548
456 G4cout << "--------------------------------- 549 G4cout << "---------------------------------\n";
457 G4cout.setf(mode, std::ios::floatfield); 550 G4cout.setf(mode, std::ios::floatfield);
458 G4cout.precision(prec); 551 G4cout.precision(prec);
459 552
460 G4int histo_id_refract = analysisMan->GetH1I 553 G4int histo_id_refract = analysisMan->GetH1Id("Fresnel refraction");
461 G4int histo_id_reflect = analysisMan->GetH1I 554 G4int histo_id_reflect = analysisMan->GetH1Id("Fresnel reflection plus TIR");
462 G4int histo_id_spike = analysisMan->GetH1Id( << 555 G4int histo_id_spike = analysisMan->GetH1Id("Spike reflection");
463 G4int histo_id_absorption = analysisMan->Get 556 G4int histo_id_absorption = analysisMan->GetH1Id("Absorption");
464 557
465 if (analysisMan->GetH1Activation(histo_id_re << 558 if(analysisMan->GetH1Activation(histo_id_refract) &&
466 && analysisMan->GetH1Activation(histo_id << 559 analysisMan->GetH1Activation(histo_id_reflect))
467 { 560 {
468 G4double rindex1 = << 561 G4double rindex1 = det->GetTankMaterial()
469 det->GetTankMaterial()->GetMaterialPrope << 562 ->GetMaterialPropertiesTable()
470 G4double rindex2 = << 563 ->GetProperty(kRINDEX)
471 det->GetWorldMaterial()->GetMaterialProp << 564 ->Value(fEkin);
>> 565 G4double rindex2 = det->GetWorldMaterial()
>> 566 ->GetMaterialPropertiesTable()
>> 567 ->GetProperty(kRINDEX)
>> 568 ->Value(fEkin);
472 569
473 auto histo_refract = analysisMan->GetH1(hi 570 auto histo_refract = analysisMan->GetH1(histo_id_refract);
474 auto histo_reflect = analysisMan->GetH1(hi 571 auto histo_reflect = analysisMan->GetH1(histo_id_reflect);
475 // std::vector<G4double> refract; 572 // std::vector<G4double> refract;
476 std::vector<G4double> reflect; 573 std::vector<G4double> reflect;
477 // std::vector<G4double> tir; 574 // std::vector<G4double> tir;
478 std::vector<G4double> tot; 575 std::vector<G4double> tot;
479 for (size_t i = 0; i < histo_refract->axis << 576 for(size_t i = 0; i < histo_refract->axis().bins(); ++i)
>> 577 {
480 // refract.push_back(histo_refract->bin_ 578 // refract.push_back(histo_refract->bin_height(i));
481 reflect.push_back(histo_reflect->bin_hei 579 reflect.push_back(histo_reflect->bin_height(i));
482 // tir.push_back(histo_TIR->bin_height(i 580 // tir.push_back(histo_TIR->bin_height(i));
483 tot.push_back(histo_refract->bin_height( << 581 tot.push_back(histo_refract->bin_height(i) +
>> 582 histo_reflect->bin_height(i));
484 } 583 }
485 584
486 // find Brewster angle: Rp = 0 585 // find Brewster angle: Rp = 0
487 // need enough statistics for this method 586 // need enough statistics for this method to work
488 G4double min_angle = -1.; 587 G4double min_angle = -1.;
489 G4double min_val = DBL_MAX; << 588 G4double min_val = DBL_MAX;
490 G4double bin_width = 0.; 589 G4double bin_width = 0.;
491 for (size_t i = 0; i < reflect.size(); ++i << 590 for(size_t i = 0; i < reflect.size(); ++i)
492 if (reflect[i] < min_val) { << 591 {
493 min_val = reflect[i]; << 592 if(reflect[i] < min_val)
>> 593 {
>> 594 min_val = reflect[i];
494 min_angle = histo_reflect->axis().bin_ 595 min_angle = histo_reflect->axis().bin_lower_edge(i);
495 bin_width = << 596 bin_width = histo_reflect->axis().bin_upper_edge(i) -
496 histo_reflect->axis().bin_upper_edge << 597 histo_reflect->axis().bin_lower_edge(i);
497 min_angle += bin_width / 2.; 598 min_angle += bin_width / 2.;
498 } 599 }
499 } 600 }
500 G4cout << "Polarization of primary optical << 601 G4cout << "Polarization of primary optical photons: "
501 << G4endl; << 602 << fPolarization / deg << " deg." << G4endl;
502 if (fPolarized && fPolarization == 0.0) { << 603 if(fPolarized && fPolarization == 0.0)
503 G4cout << "Reflectance shows a minimum a << 604 {
>> 605 G4cout << "Reflectance shows a minimum at: " << min_angle << " +/- "
>> 606 << bin_width / 2;
504 G4cout << " deg. Expected Brewster angle 607 G4cout << " deg. Expected Brewster angle: "
505 << (360. / CLHEP::twopi) * std::a << 608 << (360. / CLHEP::twopi) * std::atan(rindex2 / rindex1)
>> 609 << " deg. " << G4endl;
506 } 610 }
507 611
508 // find angle of total internal reflection 612 // find angle of total internal reflection: T -> 0
509 // last bin for T > 0 613 // last bin for T > 0
510 min_angle = -1.; 614 min_angle = -1.;
511 min_val = DBL_MAX; << 615 min_val = DBL_MAX;
512 for (size_t i = 0; i < histo_refract->axis << 616 for(size_t i = 0; i < histo_refract->axis().bins() - 1; ++i)
513 if (histo_refract->bin_height(i) > 0. && << 617 {
>> 618 if(histo_refract->bin_height(i) > 0. &&
>> 619 histo_refract->bin_height(i + 1) == 0.)
>> 620 {
514 min_angle = histo_refract->axis().bin_ 621 min_angle = histo_refract->axis().bin_lower_edge(i);
515 bin_width = << 622 bin_width = histo_reflect->axis().bin_upper_edge(i) -
516 histo_reflect->axis().bin_upper_edge << 623 histo_reflect->axis().bin_lower_edge(i);
517 min_angle += bin_width / 2.; 624 min_angle += bin_width / 2.;
518 break; 625 break;
519 } 626 }
520 } 627 }
521 if (fPolarized) { << 628 if(fPolarized)
522 G4cout << "Fresnel transmission goes to << 629 {
523 << " deg." << 630 G4cout << "Fresnel transmission goes to 0 at: " << min_angle << " +/- "
524 << " Expected: " << (360. / CLHEP << 631 << bin_width / 2. << " deg."
525 << G4endl; << 632 << " Expected: "
>> 633 << (360. / CLHEP::twopi) * std::asin(rindex2 / rindex1)
>> 634 << " deg." << G4endl;
526 } 635 }
527 636
528 // Normalize the transmission/reflection h 637 // Normalize the transmission/reflection histos so that max is 1.
529 // Only if x values are the same 638 // Only if x values are the same
530 if ((analysisMan->GetH1Nbins(histo_id_refr << 639 if((analysisMan->GetH1Nbins(histo_id_refract) ==
531 && (analysisMan->GetH1Xmin(histo_id_re << 640 analysisMan->GetH1Nbins(histo_id_reflect)) &&
532 && (analysisMan->GetH1Xmax(histo_id_re << 641 (analysisMan->GetH1Xmin(histo_id_refract) ==
>> 642 analysisMan->GetH1Xmin(histo_id_reflect)) &&
>> 643 (analysisMan->GetH1Xmax(histo_id_refract) ==
>> 644 analysisMan->GetH1Xmax(histo_id_reflect)))
533 { 645 {
534 unsigned int ent; 646 unsigned int ent;
535 G4double sw; 647 G4double sw;
536 G4double sw2; 648 G4double sw2;
537 G4double sx2; 649 G4double sx2;
538 G4double sx2w; 650 G4double sx2w;
539 for (size_t bin = 0; bin < histo_refract << 651 for(size_t bin = 0; bin < histo_refract->axis().bins(); ++bin)
>> 652 {
540 // "bin+1" below because bin 0 is unde 653 // "bin+1" below because bin 0 is underflow bin
541 // NB. We are ignoring underflow/overf 654 // NB. We are ignoring underflow/overflow bins
542 histo_refract->get_bin_content(bin + 1 655 histo_refract->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
543 if (tot[bin] > 0) { << 656 if(tot[bin] > 0)
>> 657 {
544 sw /= tot[bin]; 658 sw /= tot[bin];
545 // bin error is sqrt(sw2) 659 // bin error is sqrt(sw2)
546 sw2 /= (tot[bin] * tot[bin]); 660 sw2 /= (tot[bin] * tot[bin]);
547 sx2 /= (tot[bin] * tot[bin]); 661 sx2 /= (tot[bin] * tot[bin]);
548 sx2w /= (tot[bin] * tot[bin]); 662 sx2w /= (tot[bin] * tot[bin]);
549 histo_refract->set_bin_content(bin + 663 histo_refract->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
550 } 664 }
551 665
552 histo_reflect->get_bin_content(bin + 1 666 histo_reflect->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
553 if (tot[bin] > 0) { << 667 if(tot[bin] > 0)
>> 668 {
554 sw /= tot[bin]; 669 sw /= tot[bin];
555 // bin error is sqrt(sw2) 670 // bin error is sqrt(sw2)
556 sw2 /= (tot[bin] * tot[bin]); 671 sw2 /= (tot[bin] * tot[bin]);
557 sx2 /= (tot[bin] * tot[bin]); 672 sx2 /= (tot[bin] * tot[bin]);
558 sx2w /= (tot[bin] * tot[bin]); 673 sx2w /= (tot[bin] * tot[bin]);
559 histo_reflect->set_bin_content(bin + 674 histo_reflect->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
560 } 675 }
561 676
562 G4int histo_id_fresnelrefl = analysisM << 677 G4int histo_id_fresnelrefl =
>> 678 analysisMan->GetH1Id("Fresnel reflection");
563 auto histo_fresnelreflect = analysisMa 679 auto histo_fresnelreflect = analysisMan->GetH1(histo_id_fresnelrefl);
564 histo_fresnelreflect->get_bin_content( << 680 histo_fresnelreflect->get_bin_content(bin + 1, ent, sw, sw2, sx2,
565 if (tot[bin] > 0) { << 681 sx2w);
>> 682 if(tot[bin] > 0)
>> 683 {
566 sw /= tot[bin]; 684 sw /= tot[bin];
567 // bin error is sqrt(sw2) 685 // bin error is sqrt(sw2)
568 sw2 /= (tot[bin] * tot[bin]); 686 sw2 /= (tot[bin] * tot[bin]);
569 sx2 /= (tot[bin] * tot[bin]); 687 sx2 /= (tot[bin] * tot[bin]);
570 sx2w /= (tot[bin] * tot[bin]); 688 sx2w /= (tot[bin] * tot[bin]);
571 histo_fresnelreflect->set_bin_conten << 689 histo_fresnelreflect->set_bin_content(bin + 1, ent, sw, sw2, sx2,
>> 690 sx2w);
572 } 691 }
573 692
574 G4int histo_id_TIR = analysisMan->GetH << 693 G4int histo_id_TIR =
>> 694 analysisMan->GetH1Id("Total internal reflection");
575 auto histo_TIR = analysisMan->GetH1(hi 695 auto histo_TIR = analysisMan->GetH1(histo_id_TIR);
576 if (analysisMan->GetH1Activation(histo << 696 if(analysisMan->GetH1Activation(histo_id_TIR))
>> 697 {
577 histo_TIR->get_bin_content(bin + 1, 698 histo_TIR->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
578 if (tot[bin] > 0) { << 699 if(tot[bin] > 0)
>> 700 {
579 sw /= tot[bin]; 701 sw /= tot[bin];
580 // bin error is sqrt(sw2) 702 // bin error is sqrt(sw2)
581 sw2 /= (tot[bin] * tot[bin]); 703 sw2 /= (tot[bin] * tot[bin]);
582 sx2 /= (tot[bin] * tot[bin]); 704 sx2 /= (tot[bin] * tot[bin]);
583 sx2w /= (tot[bin] * tot[bin]); 705 sx2w /= (tot[bin] * tot[bin]);
584 histo_TIR->set_bin_content(bin + 1 706 histo_TIR->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
585 } 707 }
586 } 708 }
587 } 709 }
588 } 710 }
589 else { << 711 else
>> 712 {
590 G4cout << "Not going to normalize refrac 713 G4cout << "Not going to normalize refraction and reflection "
591 << "histograms because bins are n 714 << "histograms because bins are not the same." << G4endl;
592 } 715 }
593 } 716 }
594 717
595 // complex index of refraction; have spike r 718 // complex index of refraction; have spike reflection and absorption
596 // Only works for polished surfaces. Ground 719 // Only works for polished surfaces. Ground surfaces neglected.
597 else if (analysisMan->GetH1Activation(histo_ << 720 else if(analysisMan->GetH1Activation(histo_id_absorption) &&
598 && analysisMan->GetH1Activation(his << 721 analysisMan->GetH1Activation(histo_id_spike))
599 { 722 {
600 auto histo_spike = analysisMan->GetH1(hist << 723 auto histo_spike = analysisMan->GetH1(histo_id_spike);
601 auto histo_absorption = analysisMan->GetH1 724 auto histo_absorption = analysisMan->GetH1(histo_id_absorption);
602 725
603 std::vector<G4double> tot; 726 std::vector<G4double> tot;
604 for (size_t i = 0; i < histo_absorption->a << 727 for(size_t i = 0; i < histo_absorption->axis().bins(); ++i)
605 tot.push_back(histo_absorption->bin_heig << 728 {
>> 729 tot.push_back(histo_absorption->bin_height(i) +
>> 730 histo_spike->bin_height(i));
606 } 731 }
607 732
608 if ((analysisMan->GetH1Nbins(histo_id_abso << 733 if((analysisMan->GetH1Nbins(histo_id_absorption) ==
609 && (analysisMan->GetH1Xmin(histo_id_ab << 734 analysisMan->GetH1Nbins(histo_id_spike)) &&
610 && (analysisMan->GetH1Xmax(histo_id_ab << 735 (analysisMan->GetH1Xmin(histo_id_absorption) ==
>> 736 analysisMan->GetH1Xmin(histo_id_spike)) &&
>> 737 (analysisMan->GetH1Xmax(histo_id_absorption) ==
>> 738 analysisMan->GetH1Xmax(histo_id_spike)))
611 { 739 {
612 unsigned int ent; 740 unsigned int ent;
613 G4double sw; 741 G4double sw;
614 G4double sw2; 742 G4double sw2;
615 G4double sx2; 743 G4double sx2;
616 G4double sx2w; 744 G4double sx2w;
617 for (size_t bin = 0; bin < histo_absorpt << 745 for(size_t bin = 0; bin < histo_absorption->axis().bins(); ++bin)
>> 746 {
618 // "bin+1" below because bin 0 is unde 747 // "bin+1" below because bin 0 is underflow bin
619 // NB. We are ignoring underflow/overf 748 // NB. We are ignoring underflow/overflow bins
620 histo_absorption->get_bin_content(bin 749 histo_absorption->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
621 if (tot[bin] > 0) { << 750 if(tot[bin] > 0)
>> 751 {
622 sw /= tot[bin]; 752 sw /= tot[bin];
623 // bin error is sqrt(sw2) 753 // bin error is sqrt(sw2)
624 sw2 /= (tot[bin] * tot[bin]); 754 sw2 /= (tot[bin] * tot[bin]);
625 sx2 /= (tot[bin] * tot[bin]); 755 sx2 /= (tot[bin] * tot[bin]);
626 sx2w /= (tot[bin] * tot[bin]); 756 sx2w /= (tot[bin] * tot[bin]);
627 histo_absorption->set_bin_content(bi 757 histo_absorption->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
628 } 758 }
629 759
630 histo_spike->get_bin_content(bin + 1, 760 histo_spike->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
631 if (tot[bin] > 0) { << 761 if(tot[bin] > 0)
>> 762 {
632 sw /= tot[bin]; 763 sw /= tot[bin];
633 // bin error is sqrt(sw2) 764 // bin error is sqrt(sw2)
634 sw2 /= (tot[bin] * tot[bin]); 765 sw2 /= (tot[bin] * tot[bin]);
635 sx2 /= (tot[bin] * tot[bin]); 766 sx2 /= (tot[bin] * tot[bin]);
636 sx2w /= (tot[bin] * tot[bin]); 767 sx2w /= (tot[bin] * tot[bin]);
637 histo_spike->set_bin_content(bin + 1 768 histo_spike->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
638 } 769 }
639 } 770 }
640 } 771 }
641 else { << 772 else
>> 773 {
642 G4cout << "Not going to normalize spike 774 G4cout << "Not going to normalize spike reflection and absorption "
643 << "histograms because bins are n 775 << "histograms because bins are not the same." << G4endl;
644 } 776 }
645 } 777 }
646 } 778 }
647 779