| 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 electromagnetic/TestEm9/src/HistoMan 26 /// \file electromagnetic/TestEm9/src/HistoManager.cc
27 /// \brief Implementation of the HistoManager 27 /// \brief Implementation of the HistoManager class
28 // 28 //
29 // 29 //
30 //-------------------------------------------- 30 //---------------------------------------------------------------------------
31 // 31 //
32 // ClassName: HistoManager 32 // ClassName: HistoManager
33 // 33 //
34 // 34 //
35 // Author: V.Ivanchenko 30/01/01 35 // Author: V.Ivanchenko 30/01/01
36 // 36 //
37 //-------------------------------------------- 37 //----------------------------------------------------------------------------
38 // 38 //
39 39
40 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 42
43 #include "HistoManager.hh" 43 #include "HistoManager.hh"
44 << 44 #include "G4MaterialCutsCouple.hh"
45 #include "EmAcceptance.hh" <<
46 #include "Histo.hh" <<
47 <<
48 #include "G4Electron.hh" <<
49 #include "G4EmProcessSubType.hh" 45 #include "G4EmProcessSubType.hh"
>> 46 #include "G4VProcess.hh"
>> 47 #include "G4VEmProcess.hh"
>> 48 #include "G4VEnergyLossProcess.hh"
>> 49 #include "G4UnitsTable.hh"
>> 50 #include "Histo.hh"
>> 51 #include "EmAcceptance.hh"
50 #include "G4Gamma.hh" 52 #include "G4Gamma.hh"
51 #include "G4GammaGeneralProcess.hh" << 53 #include "G4Electron.hh"
52 #include "G4MaterialCutsCouple.hh" <<
53 #include "G4Positron.hh" 54 #include "G4Positron.hh"
54 #include "G4SystemOfUnits.hh" 55 #include "G4SystemOfUnits.hh"
55 #include "G4UnitsTable.hh" << 56 #include "G4GammaGeneralProcess.hh"
56 #include "G4VEmProcess.hh" <<
57 #include "G4VEnergyLossProcess.hh" <<
58 #include "G4VProcess.hh" <<
59 57
60 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 59
62 HistoManager* HistoManager::fManager = nullptr 60 HistoManager* HistoManager::fManager = nullptr;
63 61
64 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 63
66 HistoManager* HistoManager::GetPointer() 64 HistoManager* HistoManager::GetPointer()
67 { 65 {
68 if (nullptr == fManager) { << 66 if(nullptr == fManager) {
69 fManager = new HistoManager(); 67 fManager = new HistoManager();
70 } 68 }
71 return fManager; 69 return fManager;
72 } 70 }
73 71
74 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 73
76 HistoManager::HistoManager() 74 HistoManager::HistoManager()
77 : fGamma(G4Gamma::Gamma()), << 75 : fGamma(G4Gamma::Gamma()),
78 fElectron(G4Electron::Electron()), << 76 fElectron(G4Electron::Electron()),
79 fPositron(G4Positron::Positron()), << 77 fPositron(G4Positron::Positron()),
80 fHisto(new Histo()) << 78 fHisto(new Histo())
81 { 79 {
82 fVerbose = 1; 80 fVerbose = 1;
83 fEvt1 = -1; << 81 fEvt1 = -1;
84 fEvt2 = -1; << 82 fEvt2 = -1;
85 fNmax = 3; << 83 fNmax = 3;
86 fMaxEnergy = 50.0 * MeV; << 84 fMaxEnergy = 50.0*MeV;
87 fBeamEnergy = 1. * GeV; << 85 fBeamEnergy= 1.*GeV;
88 fMaxEnergyAbs = 10. * MeV; << 86 fMaxEnergyAbs = 10.*MeV;
89 fBinsE = 100; 87 fBinsE = 100;
90 fBinsEA = 40; << 88 fBinsEA= 40;
91 fBinsED = 100; << 89 fBinsED= 100;
92 fNHisto = 13; 90 fNHisto = 13;
93 91
94 BookHisto(); 92 BookHisto();
95 } 93 }
96 94
97 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 96
99 HistoManager::~HistoManager() 97 HistoManager::~HistoManager()
100 { 98 {
101 delete fHisto; 99 delete fHisto;
102 } 100 }
103 101
104 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 103
106 void HistoManager::BookHisto() 104 void HistoManager::BookHisto()
107 { 105 {
108 fHisto->Add1D("10", "Evis/E0 in central crys << 106 fHisto->Add1D("10","Evis/E0 in central crystal",fBinsED,0.0,1,1.0);
109 fHisto->Add1D("11", "Evis/E0 in 3x3", fBinsE << 107 fHisto->Add1D("11","Evis/E0 in 3x3",fBinsED,0.0,1.0,1.0);
110 fHisto->Add1D("12", "Evis/E0 in 5x5", fBinsE << 108 fHisto->Add1D("12","Evis/E0 in 5x5",fBinsED,0.0,1.0,1.0);
111 fHisto->Add1D("13", "Energy (MeV) of delta-e << 109 fHisto->Add1D("13","Energy (MeV) of delta-electrons",
112 fHisto->Add1D("14", "Energy (MeV) of gammas" << 110 fBinsE,0.0,fMaxEnergy,MeV);
113 fHisto->Add1D("15", "Energy (MeV) in abs1", << 111 fHisto->Add1D("14","Energy (MeV) of gammas",fBinsE,0.0,fMaxEnergy,MeV);
114 fHisto->Add1D("16", "Energy (MeV) in abs2", << 112 fHisto->Add1D("15","Energy (MeV) in abs1",fBinsEA,0.0,fMaxEnergyAbs,MeV);
115 fHisto->Add1D("17", "Energy (MeV) in abs3", << 113 fHisto->Add1D("16","Energy (MeV) in abs2",fBinsEA,0.0,fMaxEnergyAbs,MeV);
116 fHisto->Add1D("18", "Energy (MeV) in abs4", << 114 fHisto->Add1D("17","Energy (MeV) in abs3",fBinsEA,0.0,fMaxEnergyAbs,MeV);
117 fHisto->Add1D("19", "Number of vertex hits", << 115 fHisto->Add1D("18","Energy (MeV) in abs4",fBinsEA,0.0,fMaxEnergyAbs,MeV);
118 fHisto->Add1D("20", "E1/E9 Ratio", fBinsED, << 116 fHisto->Add1D("19","Number of vertex hits",20,-0.5,19.5,1.0);
119 fHisto->Add1D("21", "E1/E25 Ratio", fBinsED, << 117 fHisto->Add1D("20","E1/E9 Ratio",fBinsED,0.0,1,1.0);
120 fHisto->Add1D("22", "E9/E25 Ratio", fBinsED, << 118 fHisto->Add1D("21","E1/E25 Ratio",fBinsED,0.0,1.0,1.0);
>> 119 fHisto->Add1D("22","E9/E25 Ratio",fBinsED,0.0,1.0,1.0);
121 } 120 }
122 121
123 //....oooOO0OOooo........oooOO0OOooo........oo 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 123
125 void HistoManager::BeginOfRun() 124 void HistoManager::BeginOfRun()
126 { 125 {
127 // initilise scoring 126 // initilise scoring
128 fEvt = 0; << 127 fEvt = 0;
129 fElec = 0; 128 fElec = 0;
130 fPosit = 0; << 129 fPosit= 0;
131 fGam = 0; << 130 fGam = 0;
132 fStep = 0; 131 fStep = 0;
133 fLowe = 0; 132 fLowe = 0;
134 133
135 for (G4int i = 0; i < 6; i++) { << 134 for(G4int i=0; i<6; i++) {
136 fStat[i] = 0; 135 fStat[i] = 0;
137 fEdep[i] = 0.0; 136 fEdep[i] = 0.0;
138 fErms[i] = 0.0; 137 fErms[i] = 0.0;
139 if (i < 3) { << 138 if(i < 3) {
140 fEdeptr[i] = 0.0; 139 fEdeptr[i] = 0.0;
141 fErmstr[i] = 0.0; 140 fErmstr[i] = 0.0;
142 } 141 }
143 } 142 }
144 143
145 // initialise counters 144 // initialise counters
146 fBrem.resize(93, 0.0); << 145 fBrem.resize(93,0.0);
147 fPhot.resize(93, 0.0); << 146 fPhot.resize(93,0.0);
148 fComp.resize(93, 0.0); << 147 fComp.resize(93,0.0);
149 fConv.resize(93, 0.0); << 148 fConv.resize(93,0.0);
150 149
151 // initialise acceptance - by default is not 150 // initialise acceptance - by default is not applied
152 for (G4int i = 0; i < fNmax; i++) { << 151 for(G4int i=0; i<fNmax; i++) {
153 fEdeptrue[i] = 1.0; 152 fEdeptrue[i] = 1.0;
154 fRmstrue[i] = 1.0; << 153 fRmstrue[i] = 1.0;
155 fLimittrue[i] = 10.; << 154 fLimittrue[i]= 10.;
156 } 155 }
157 156
158 if (fHisto->IsActive()) { << 157 if(fHisto->IsActive()) {
159 for (G4int i = 0; i < fNHisto; ++i) { << 158 for(G4int i=0; i<fNHisto; ++i) {fHisto->Activate(i, true); }
160 fHisto->Activate(i, true); <<
161 } <<
162 fHisto->Book(); 159 fHisto->Book();
163 160
164 if (fVerbose > 0) { << 161 if(fVerbose > 0) {
165 G4cout << "HistoManager: Histograms are << 162 G4cout << "HistoManager: Histograms are booked and run has been started"
>> 163 << G4endl;
166 } 164 }
167 } 165 }
168 } 166 }
169 167
170 //....oooOO0OOooo........oooOO0OOooo........oo 168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 169
172 void HistoManager::EndOfRun(G4int runID) 170 void HistoManager::EndOfRun(G4int runID)
173 { 171 {
174 G4cout << "HistoManager: End of run actions << 172
>> 173 G4cout << "HistoManager: End of run actions are started RunID# "
>> 174 << runID << G4endl;
175 G4String nam[6] = {"1x1", "3x3", "5x5", "E1/ 175 G4String nam[6] = {"1x1", "3x3", "5x5", "E1/E9 ", "E1/E25", "E9/E25"};
176 176
177 // average 177 // average
178 178
179 G4cout << "================================= << 179 G4cout<<"================================================================="
>> 180 <<G4endl;
180 G4double x = (G4double)fEvt; 181 G4double x = (G4double)fEvt;
181 if (fEvt > 0) x = 1.0 / x; << 182 if(fEvt > 0) x = 1.0/x;
182 G4int j; 183 G4int j;
183 for (j = 0; j < fNmax; j++) { << 184 for(j=0; j<fNmax; j++) {
>> 185
184 // total mean 186 // total mean
185 fEdep[j] *= x; 187 fEdep[j] *= x;
186 G4double y = fErms[j] * x - fEdep[j] * fEd << 188 G4double y = fErms[j]*x - fEdep[j]*fEdep[j];
187 if (y < 0.0) y = 0.0; << 189 if(y < 0.0) y = 0.0;
188 fErms[j] = std::sqrt(y); 190 fErms[j] = std::sqrt(y);
189 191
190 // trancated mean 192 // trancated mean
191 G4double xx = G4double(fStat[j]); 193 G4double xx = G4double(fStat[j]);
192 if (xx > 0.0) xx = 1.0 / xx; << 194 if(xx > 0.0) xx = 1.0/xx;
193 fEdeptr[j] *= xx; 195 fEdeptr[j] *= xx;
194 y = fErmstr[j] * xx - fEdeptr[j] * fEdeptr << 196 y = fErmstr[j]*xx - fEdeptr[j]*fEdeptr[j];
195 if (y < 0.0) y = 0.0; << 197 if(y < 0.0) y = 0.0;
196 fErmstr[j] = std::sqrt(y); 198 fErmstr[j] = std::sqrt(y);
197 } 199 }
198 G4double xe = x * (G4double)fElec; << 200 G4double xe = x*(G4double)fElec;
199 G4double xg = x * (G4double)fGam; << 201 G4double xg = x*(G4double)fGam;
200 G4double xp = x * (G4double)fPosit; << 202 G4double xp = x*(G4double)fPosit;
201 G4double xs = x * fStep; << 203 G4double xs = x*fStep;
202 204
203 G4double f = 100. * std::sqrt(fBeamEnergy / << 205 G4double f = 100.*std::sqrt(fBeamEnergy/GeV);
204 206
205 G4cout << "Number of events " << << 207 G4cout << "Number of events " << fEvt <<G4endl;
206 G4cout << std::setprecision(4) << "Average n 208 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
207 G4cout << std::setprecision(4) << "Average n 209 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
208 G4cout << std::setprecision(4) << "Average n 210 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
209 G4cout << std::setprecision(4) << "Average n 211 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
210 << 212
211 for (j = 0; j < 3; ++j) { << 213 for(j=0; j<3; j++) {
212 G4double ex = fEdeptr[j]; 214 G4double ex = fEdeptr[j];
213 G4double sx = fErmstr[j]; 215 G4double sx = fErmstr[j];
214 G4double xx = G4double(fStat[j]); << 216 G4double xx= G4double(fStat[j]);
215 if (xx > 0.0) xx = 1.0 / xx; << 217 if(xx > 0.0) xx = 1.0/xx;
216 G4double r = sx * std::sqrt(xx); << 218 G4double r = sx*std::sqrt(xx);
217 G4cout << std::setprecision(4) << "Edep " << 219 G4cout << std::setprecision(4) << "Edep " << nam[j]
218 << r; << 220 << " = " << ex
219 if (ex > 0.1) G4cout << " res= " << f * << 221 << " +- " << r;
>> 222 if(ex > 0.1) G4cout << " res= " << f*sx/ex << " % " << fStat[j];
220 G4cout << G4endl; 223 G4cout << G4endl;
221 } 224 }
222 if (fLimittrue[0] < 10. || fLimittrue[1] < 1 << 225 if(fLimittrue[0] < 10. || fLimittrue[1] < 10. || fLimittrue[2] < 10.) {
223 G4cout << "=========== Mean values withou << 226 G4cout
224 for (j = 0; j < fNmax; j++) { << 227 <<"=========== Mean values without trancating ====================="<<G4endl;
>> 228 for(j=0; j<fNmax; j++) {
225 G4double ex = fEdep[j]; 229 G4double ex = fEdep[j];
226 G4double sx = fErms[j]; 230 G4double sx = fErms[j];
227 G4double rx = sx * std::sqrt(x); << 231 G4double rx = sx*std::sqrt(x);
228 G4cout << std::setprecision(4) << "Edep << 232 G4cout << std::setprecision(4) << "Edep " << nam[j]
229 << rx; << 233 << " = " << ex
230 if (ex > 0.0) G4cout << " res= " << f << 234 << " +- " << rx;
>> 235 if(ex > 0.0) G4cout << " res= " << f*sx/ex << " %";
231 G4cout << G4endl; 236 G4cout << G4endl;
232 } 237 }
233 } 238 }
234 G4cout << "=========== Ratios without tranc << 239 G4cout
235 for (j = 3; j < 6; ++j) { << 240 <<"=========== Ratios without trancating ==========================="<<G4endl;
>> 241 for(j=3; j<6; j++) {
236 G4double e = fEdep[j]; 242 G4double e = fEdep[j];
237 G4double xx = G4double(fStat[j]); << 243 G4double xx= G4double(fStat[j]);
238 if (xx > 0.0) xx = 1.0 / xx; << 244 if(xx > 0.0) xx = 1.0/xx;
239 e *= xx; 245 e *= xx;
240 G4double y = fErms[j] * xx - e * e; << 246 G4double y = fErms[j]*xx - e*e;
241 G4double r = 0.0; 247 G4double r = 0.0;
242 if (y > 0.0) r = std::sqrt(y * xx); << 248 if(y > 0.0) r = std::sqrt(y*xx);
243 G4cout << " " << nam[j] << " = << 249 G4cout << " " << nam[j] << " = " << e
>> 250 << " +- " << r;
244 G4cout << G4endl; 251 G4cout << G4endl;
245 } 252 }
246 G4cout << std::setprecision(4) << "Beam Ener << 253 G4cout << std::setprecision(4) << "Beam Energy " << fBeamEnergy/GeV
247 << G4endl; << 254 << " GeV" << G4endl;
248 if (fLowe > 0) G4cout << "Number of events E << 255 if(fLowe > 0) G4cout << "Number of events E/E0<0.8 " << fLowe << G4endl;
249 G4cout << "================================= << 256 G4cout
250 G4cout << G4endl; << 257 <<"=================================================================="<<G4endl;
>> 258 G4cout<<G4endl;
251 259
252 // normalise histograms 260 // normalise histograms
253 if (fHisto->IsActive()) { << 261 if(fHisto->IsActive()) {
254 for (G4int i = 0; i < fNHisto; ++i) { << 262 for(G4int i=0; i<fNHisto; i++) {
255 fHisto->ScaleH1(i, x); << 263 fHisto->ScaleH1(i,x);
256 } 264 }
257 fHisto->Save(); 265 fHisto->Save();
258 } 266 }
259 if (0 < runID) { << 267 if(0 < runID) { return; }
260 return; <<
261 } <<
262 268
263 // Acceptance only for the first run 269 // Acceptance only for the first run
264 EmAcceptance acc; 270 EmAcceptance acc;
265 G4bool isStarted = false; 271 G4bool isStarted = false;
266 for (j = 0; j < fNmax; j++) { << 272 for (j=0; j<fNmax; j++) {
>> 273
267 G4double ltrue = fLimittrue[j]; 274 G4double ltrue = fLimittrue[j];
268 if (ltrue < DBL_MAX) { 275 if (ltrue < DBL_MAX) {
269 if (!isStarted) { 276 if (!isStarted) {
270 acc.BeginOfAcceptance("Crystal Calorim << 277 acc.BeginOfAcceptance("Crystal Calorimeter",fEvt);
271 isStarted = true; 278 isStarted = true;
272 } 279 }
273 G4double etrue = fEdeptrue[j]; 280 G4double etrue = fEdeptrue[j];
274 G4double rtrue = fRmstrue[j]; 281 G4double rtrue = fRmstrue[j];
275 acc.EmAcceptanceGauss("Edep" + nam[j], f << 282 acc.EmAcceptanceGauss("Edep"+nam[j],fEvt,fEdeptr[j],etrue,rtrue,ltrue);
276 acc.EmAcceptanceGauss("Erms" + nam[j], f << 283 acc.EmAcceptanceGauss("Erms"+nam[j],fEvt,fErmstr[j],rtrue,rtrue,2.0*ltrue);
277 } 284 }
278 } 285 }
279 if (isStarted) acc.EndOfAcceptance(); << 286 if(isStarted) acc.EndOfAcceptance();
280 287
281 // atom frequency 288 // atom frequency
282 G4cout << " Z bremsstrahlung photoeffect 289 G4cout << " Z bremsstrahlung photoeffect compton conversion" << G4endl;
283 for (j = 1; j < 93; ++j) { << 290 for(j=1; j<93; ++j) {
284 G4int n1 = G4int(fBrem[j] * x); << 291 G4int n1 = G4int(fBrem[j]*x);
285 G4int n2 = G4int(fPhot[j] * x); << 292 G4int n2 = G4int(fPhot[j]*x);
286 G4int n3 = G4int(fComp[j] * x); << 293 G4int n3 = G4int(fComp[j]*x);
287 G4int n4 = G4int(fConv[j] * x); << 294 G4int n4 = G4int(fConv[j]*x);
288 if (n1 + n2 + n3 + n4 > 0) { << 295 if(n1 + n2 + n3 + n4 > 0) {
289 G4cout << std::setw(4) << j << std::setw << 296 G4cout << std::setw(4) << j << std::setw(12) << n1 << std::setw(12) << n2
290 << n3 << std::setw(12) << n4 << G << 297 << std::setw(12) << n3 << std::setw(12) << n4 << G4endl;
291 } 298 }
292 } 299 }
293 } 300 }
294 301
295 //....oooOO0OOooo........oooOO0OOooo........oo 302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296 303
297 void HistoManager::BeginOfEvent() 304 void HistoManager::BeginOfEvent()
298 { 305 {
299 ++fEvt; 306 ++fEvt;
300 307
301 fEabs1 = 0.0; << 308 fEabs1 = 0.0;
302 fEabs2 = 0.0; << 309 fEabs2 = 0.0;
303 fEabs3 = 0.0; << 310 fEabs3 = 0.0;
304 fEabs4 = 0.0; << 311 fEabs4 = 0.0;
305 fEvertex.clear(); 312 fEvertex.clear();
306 fNvertex.clear(); 313 fNvertex.clear();
307 for (G4int i = 0; i < 25; i++) { << 314 for (G4int i=0; i<25; i++) {
308 fE[i] = 0.0; 315 fE[i] = 0.0;
309 } 316 }
310 } 317 }
311 318
312 //....oooOO0OOooo........oooOO0OOooo........oo 319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
313 320
314 void HistoManager::EndOfEvent() 321 void HistoManager::EndOfEvent()
315 { 322 {
316 G4double e9 = 0.0; 323 G4double e9 = 0.0;
317 G4double e25 = 0.0; << 324 G4double e25= 0.0;
318 for (G4int i = 0; i < 25; i++) { << 325 for (G4int i=0; i<25; i++) {
319 fE[i] /= fBeamEnergy; 326 fE[i] /= fBeamEnergy;
320 e25 += fE[i]; 327 e25 += fE[i];
321 if ((6 <= i && 8 >= i) || (11 <= i && 13 > << 328 if( ( 6<=i && 8>=i) || (11<=i && 13>=i) || (16<=i && 18>=i)) e9 += fE[i];
322 } 329 }
323 330
324 if (1 < fVerbose && e25 < 0.8) { << 331 if(1 < fVerbose && e25 < 0.8) {
325 ++fLowe; 332 ++fLowe;
326 G4cout << "### in the event# " << fEvt << 333 G4cout << "### in the event# " << fEvt << " E25= " << e25 << G4endl;
327 } 334 }
328 335
329 // compute ratios 336 // compute ratios
330 G4double e0 = fE[12]; 337 G4double e0 = fE[12];
331 G4double e19 = 0.0; << 338 G4double e19 = 0.0;
332 G4double e125 = 0.0; 339 G4double e125 = 0.0;
333 G4double e925 = 0.0; 340 G4double e925 = 0.0;
334 if (e9 > 0.0) { << 341 if(e9 > 0.0) {
335 e19 = e0 / e9; << 342 e19 = e0/e9;
336 e125 = e0 / e25; << 343 e125 = e0/e25;
337 e925 = e9 / e25; << 344 e925 = e9/e25;
338 fEdep[3] += e19; 345 fEdep[3] += e19;
339 fErms[3] += e19 * e19; << 346 fErms[3] += e19*e19;
340 fEdep[4] += e125; 347 fEdep[4] += e125;
341 fErms[4] += e125 * e125; << 348 fErms[4] += e125*e125;
342 fEdep[5] += e925; 349 fEdep[5] += e925;
343 fErms[5] += e925 * e925; << 350 fErms[5] += e925*e925;
344 fStat[3] += 1; 351 fStat[3] += 1;
345 fStat[4] += 1; 352 fStat[4] += 1;
346 fStat[5] += 1; 353 fStat[5] += 1;
347 } 354 }
348 355
349 // Fill histo 356 // Fill histo
350 fHisto->Fill(0, e0, 1.0); << 357 fHisto->Fill(0,e0,1.0);
351 fHisto->Fill(1, e9, 1.0); << 358 fHisto->Fill(1,e9,1.0);
352 fHisto->Fill(2, e25, 1.0); << 359 fHisto->Fill(2,e25,1.0);
353 fHisto->Fill(5, fEabs1, 1.0); << 360 fHisto->Fill(5,fEabs1,1.0);
354 fHisto->Fill(6, fEabs2, 1.0); << 361 fHisto->Fill(6,fEabs2,1.0);
355 fHisto->Fill(7, fEabs3, 1.0); << 362 fHisto->Fill(7,fEabs3,1.0);
356 fHisto->Fill(8, fEabs4, 1.0); << 363 fHisto->Fill(8,fEabs4,1.0);
357 fHisto->Fill(9, G4double(fNvertex.size()), 1 << 364 fHisto->Fill(9,G4double(fNvertex.size()),1.0);
358 fHisto->Fill(10, e19, 1.0); << 365 fHisto->Fill(10,e19,1.0);
359 fHisto->Fill(11, e125, 1.0); << 366 fHisto->Fill(11,e125,1.0);
360 fHisto->Fill(12, e925, 1.0); << 367 fHisto->Fill(12,e925,1.0);
361 368
362 // compute sums 369 // compute sums
363 fEdep[0] += e0; 370 fEdep[0] += e0;
364 fErms[0] += e0 * e0; << 371 fErms[0] += e0*e0;
365 fEdep[1] += e9; 372 fEdep[1] += e9;
366 fErms[1] += e9 * e9; << 373 fErms[1] += e9*e9;
367 fEdep[2] += e25; 374 fEdep[2] += e25;
368 fErms[2] += e25 * e25; << 375 fErms[2] += e25*e25;
369 376
370 // trancated mean 377 // trancated mean
371 if (std::abs(e0 - fEdeptrue[0]) < fRmstrue[0 << 378 if(std::abs(e0-fEdeptrue[0])<fRmstrue[0]*fLimittrue[0]) {
372 fStat[0] += 1; 379 fStat[0] += 1;
373 fEdeptr[0] += e0; 380 fEdeptr[0] += e0;
374 fErmstr[0] += e0 * e0; << 381 fErmstr[0] += e0*e0;
375 } 382 }
376 if (std::abs(e9 - fEdeptrue[1]) < fRmstrue[1 << 383 if(std::abs(e9-fEdeptrue[1])<fRmstrue[1]*fLimittrue[1]) {
377 fStat[1] += 1; 384 fStat[1] += 1;
378 fEdeptr[1] += e9; 385 fEdeptr[1] += e9;
379 fErmstr[1] += e9 * e9; << 386 fErmstr[1] += e9*e9;
380 } 387 }
381 if (std::abs(e25 - fEdeptrue[2]) < fRmstrue[ << 388 if(std::abs(e25-fEdeptrue[2])<fRmstrue[2]*fLimittrue[2]) {
382 fStat[2] += 1; 389 fStat[2] += 1;
383 fEdeptr[2] += e25; 390 fEdeptr[2] += e25;
384 fErmstr[2] += e25 * e25; << 391 fErmstr[2] += e25*e25;
385 } 392 }
386 } 393 }
387 394
388 //....oooOO0OOooo........oooOO0OOooo........oo 395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389 396
390 void HistoManager::ScoreNewTrack(const G4Track 397 void HistoManager::ScoreNewTrack(const G4Track* aTrack)
391 { 398 {
392 // Save primary parameters << 399 //Save primary parameters
393 ResetTrackLength(); 400 ResetTrackLength();
394 const G4ParticleDefinition* particle = aTrac 401 const G4ParticleDefinition* particle = aTrack->GetDefinition();
395 const G4DynamicParticle* dynParticle = aTrac 402 const G4DynamicParticle* dynParticle = aTrack->GetDynamicParticle();
396 403
397 G4int pid = aTrack->GetParentID(); 404 G4int pid = aTrack->GetParentID();
398 G4double kinE = dynParticle->GetKineticEnerg 405 G4double kinE = dynParticle->GetKineticEnergy();
399 G4ThreeVector pos = aTrack->GetVertexPositio 406 G4ThreeVector pos = aTrack->GetVertexPosition();
400 407
401 // primary 408 // primary
402 if (0 == pid) { << 409 if(0 == pid) {
>> 410
403 G4double mass = 0.0; 411 G4double mass = 0.0;
404 if (particle) { << 412 if(particle) {
405 mass = particle->GetPDGMass(); 413 mass = particle->GetPDGMass();
406 } 414 }
407 415
408 G4ThreeVector dir = dynParticle->GetMoment 416 G4ThreeVector dir = dynParticle->GetMomentumDirection();
409 if (1 < fVerbose) { << 417 if(1 < fVerbose) {
410 G4cout << "TrackingAction: Primary kinE( << 418 G4cout << "TrackingAction: Primary kinE(MeV)= " << kinE/MeV
411 << "; pos= " << pos << "; dir= " << 419 << "; m(MeV)= " << mass/MeV
>> 420 << "; pos= " << pos << "; dir= " << dir << G4endl;
412 } 421 }
413 422
414 // secondary 423 // secondary
415 } << 424 } else {
416 else { <<
417 const G4VProcess* proc = aTrack->GetCreato 425 const G4VProcess* proc = aTrack->GetCreatorProcess();
418 G4int type = proc->GetProcessSubType(); 426 G4int type = proc->GetProcessSubType();
>> 427
>> 428 if(type == fGammaGeneralProcess) {
>> 429 type = static_cast<const G4GammaGeneralProcess*>(proc)->GetSubProcessSubType();
>> 430 proc = static_cast<const G4GammaGeneralProcess*>(proc)->GetSelectedProcess();
>> 431 }
419 432
420 if (type == fBremsstrahlung) { << 433 if(type == fBremsstrahlung) {
421 auto elm = static_cast<const G4VEnergyLo << 434 const G4Element* elm =
422 if (nullptr != elm) { << 435 static_cast<const G4VEnergyLossProcess*>(proc)->GetCurrentElement();
423 G4int Z = elm->GetZasInt(); << 436 if(elm) {
424 if (Z > 0 && Z < 93) { << 437 G4int Z = G4lrint(elm->GetZ());
425 fBrem[Z] += 1.0; << 438 if(Z > 0 && Z < 93) { fBrem[Z] += 1.0; }
426 } <<
427 } 439 }
428 } << 440 } else if(type == fPhotoElectricEffect) {
429 else if (type == fPhotoElectricEffect) { << 441 const G4Element* elm =
430 auto elm = static_cast<const G4VEmProces << 442 static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
431 if (nullptr != elm) { << 443 if(elm) {
432 G4int Z = elm->GetZasInt(); << 444 G4int Z = G4lrint(elm->GetZ());
433 if (Z > 0 && Z < 93) { << 445 if(Z > 0 && Z < 93) { fPhot[Z] += 1.0; }
434 fPhot[Z] += 1.0; <<
435 } <<
436 } 446 }
437 } << 447 } else if(type == fGammaConversion) {
438 else if (type == fGammaConversion) { << 448 const G4Element* elm =
439 auto elm = static_cast<const G4VEmProces << 449 static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
440 if (nullptr != elm) { << 450 if(elm) {
441 G4int Z = elm->GetZasInt(); << 451 G4int Z = G4lrint(elm->GetZ());
442 if (Z > 0 && Z < 93) { << 452 if(Z > 0 && Z < 93) { fConv[Z] += 1.0; }
443 fConv[Z] += 1.0; <<
444 } <<
445 } 453 }
446 } << 454 } else if(type == fComptonScattering) {
447 else if (type == fComptonScattering) { << 455 const G4Element* elm =
448 auto elm = static_cast<const G4VEmProces << 456 static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
449 if (nullptr != elm) { << 457 if(elm) {
450 G4int Z = elm->GetZasInt(); << 458 G4int Z = G4lrint(elm->GetZ());
451 if (Z > 0 && Z < 93) { << 459 if(Z > 0 && Z < 93) { fComp[Z] += 1.0; }
452 fComp[Z] += 1.0; <<
453 } <<
454 } 460 }
455 } 461 }
456 462
457 // delta-electron 463 // delta-electron
458 if (particle == fElectron) { 464 if (particle == fElectron) {
459 if (1 < fVerbose) { << 465 if(1 < fVerbose) {
460 G4cout << "TrackingAction: Secondary e 466 G4cout << "TrackingAction: Secondary electron " << G4endl;
461 } 467 }
462 AddDeltaElectron(dynParticle); 468 AddDeltaElectron(dynParticle);
463 } << 469
464 else if (particle == fPositron) { << 470 } else if (particle == fPositron) {
465 if (1 < fVerbose) { << 471 if(1 < fVerbose) {
466 G4cout << "TrackingAction: Secondary p 472 G4cout << "TrackingAction: Secondary positron " << G4endl;
467 } 473 }
468 AddPositron(dynParticle); 474 AddPositron(dynParticle);
469 } << 475
470 else if (particle == fGamma) { << 476 } else if (particle == fGamma) {
471 if (1 < fVerbose) { << 477 if(1 < fVerbose) {
472 G4cout << "TrackingAction: Secondary g << 478 G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
>> 479 << " E= " << kinE << G4endl;
473 } 480 }
474 AddPhoton(dynParticle); 481 AddPhoton(dynParticle);
475 } 482 }
476 } 483 }
477 } 484 }
478 485
479 //....oooOO0OOooo........oooOO0OOooo........oo 486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
480 487
481 void HistoManager::AddEnergy(G4double edep, G4 488 void HistoManager::AddEnergy(G4double edep, G4int volIndex, G4int copyNo)
482 { 489 {
483 if (1 < fVerbose) { << 490 if(1 < fVerbose) {
484 G4cout << "HistoManager::AddEnergy: e(keV) << 491 G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV
485 << "; copyNo= " << copyNo << G4endl << 492 << "; volIdx= " << volIndex
>> 493 << "; copyNo= " << copyNo
>> 494 << G4endl;
486 } 495 }
487 if (0 == volIndex) { << 496 if(0 == volIndex) {
488 fE[copyNo] += edep; 497 fE[copyNo] += edep;
489 } << 498 } else if (1 == volIndex) {
490 else if (1 == volIndex) { <<
491 fEabs1 += edep; 499 fEabs1 += edep;
492 } << 500 } else if (2 == volIndex) {
493 else if (2 == volIndex) { <<
494 fEabs2 += edep; 501 fEabs2 += edep;
495 } << 502 } else if (3 == volIndex) {
496 else if (3 == volIndex) { <<
497 fEabs3 += edep; 503 fEabs3 += edep;
498 } << 504 } else if (4 == volIndex) {
499 else if (4 == volIndex) { <<
500 fEabs4 += edep; 505 fEabs4 += edep;
501 } << 506 } else if (5 == volIndex) {
502 else if (5 == volIndex) { <<
503 G4int n = fNvertex.size(); 507 G4int n = fNvertex.size();
504 G4bool newPad = true; 508 G4bool newPad = true;
505 if (n > 0) { 509 if (n > 0) {
506 for (G4int i = 0; i < n; i++) { << 510 for(G4int i=0; i<n; i++) {
507 if (copyNo == fNvertex[i]) { 511 if (copyNo == fNvertex[i]) {
508 newPad = false; 512 newPad = false;
509 fEvertex[i] += edep; 513 fEvertex[i] += edep;
510 break; 514 break;
511 } 515 }
512 } 516 }
513 } 517 }
514 if (newPad) { << 518 if(newPad) {
515 fNvertex.push_back(copyNo); 519 fNvertex.push_back(copyNo);
516 fEvertex.push_back(edep); 520 fEvertex.push_back(edep);
517 } 521 }
518 } 522 }
519 } 523 }
520 524
521 //....oooOO0OOooo........oooOO0OOooo........oo 525 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
522 526
523 void HistoManager::AddDeltaElectron(const G4Dy 527 void HistoManager::AddDeltaElectron(const G4DynamicParticle* elec)
524 { 528 {
525 G4double e = elec->GetKineticEnergy() / MeV; << 529 G4double e = elec->GetKineticEnergy()/MeV;
526 if (e > 0.0) { << 530 if(e > 0.0) fElec++;
527 ++fElec; << 531 fHisto->Fill(3,e,1.0);
528 fHisto->Fill(3, e, 1.0); <<
529 } <<
530 } 532 }
531 533
532 //....oooOO0OOooo........oooOO0OOooo........oo 534 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
533 535
534 void HistoManager::AddPhoton(const G4DynamicPa 536 void HistoManager::AddPhoton(const G4DynamicParticle* ph)
535 { 537 {
536 G4double e = ph->GetKineticEnergy() / MeV; << 538 G4double e = ph->GetKineticEnergy()/MeV;
537 if (e > 0.0) { << 539 if(e > 0.0) fGam++;
538 ++fGam; << 540 fHisto->Fill(4,e,1.0);
539 fHisto->Fill(4, e, 1.0); <<
540 } <<
541 } 541 }
542 542
543 //....oooOO0OOooo........oooOO0OOooo........oo 543 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
544 544
545 void HistoManager::SetEdepAndRMS(G4int i, cons << 545 void HistoManager::SetEdepAndRMS(G4int i, G4ThreeVector val)
546 { 546 {
547 if (i < fNmax && i >= 0) { << 547 if(i<fNmax && i>=0) {
548 if (val[0] > 0.0) fEdeptrue[i] = val[0]; << 548 if(val[0] > 0.0) fEdeptrue[i] = val[0];
549 if (val[1] > 0.0) fRmstrue[i] = val[1]; << 549 if(val[1] > 0.0) fRmstrue[i] = val[1];
550 if (val[2] > 0.0) fLimittrue[i] = val[2]; << 550 if(val[2] > 0.0) fLimittrue[i] = val[2];
551 } 551 }
552 } 552 }
553 553
554 //....oooOO0OOooo........oooOO0OOooo........oo 554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 555
555 556