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