source: trunk/examples/extended/electromagnetic/TestEm9/src/HistoManager.cc@ 1326

Last change on this file since 1326 was 1230, checked in by garnier, 16 years ago

update to geant4.9.3
File size: 13.8 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//---------------------------------------------------------------------------
27//
28// ClassName: HistoManager
29//
30//
31// Author: V.Ivanchenko 30/01/01
32//
33//----------------------------------------------------------------------------
34//
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38
39#include "HistoManager.hh"
40#include "G4UnitsTable.hh"
41#include "Histo.hh"
42#include "EmAcceptance.hh"
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45
46HistoManager* HistoManager::fManager = 0;
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50HistoManager* HistoManager::GetPointer()
51{
52 if(!fManager) {
53 fManager = new HistoManager();
54 }
55 return fManager;
56}
57
58//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59
60HistoManager::HistoManager()
61{
62 verbose = 1;
63 nEvt1 = -1;
64 nEvt2 = -1;
65 nmax = 3;
66 histo = new Histo();
67 bookHisto();
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72HistoManager::~HistoManager()
73{
74 delete histo;
75}
76
77//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78
79void HistoManager::bookHisto()
80{
81 maxEnergy = 50.0*MeV;
82 beamEnergy= 1.*GeV;
83 maxEnergyAbs = 10.*MeV;
84 thKinE = 1.*MeV;
85 nBinsE = 100;
86 nBinsEA= 40;
87 nBinsED= 100;
88 nTuple = false;
89 nHisto = 13;
90
91 // initialise acceptance
92 for(G4int i=0; i<nmax; i++) {
93 edeptrue[i] = 1.0;
94 rmstrue[i] = 1.0;
95 limittrue[i]= DBL_MAX;
96 }
97
98 histo->add1D("10",
99 "Evis/E0 in central crystal",nBinsED,0.0,1,1.0);
100
101 histo->add1D("11",
102 "Evis/E0 in 3x3",nBinsED,0.0,1.0,1.0);
103
104 histo->add1D("12",
105 "Evis/E0 in 5x5",nBinsED,0.0,1.0,1.0);
106
107 histo->add1D("13",
108 "Energy (MeV) of delta-electrons",nBinsE,0.0,maxEnergy,MeV);
109
110 histo->add1D("14",
111 "Energy (MeV) of gammas",nBinsE,0.0,maxEnergy,MeV);
112
113 histo->add1D("15",
114 "Energy (MeV) in abs1",nBinsEA,0.0,maxEnergyAbs,MeV);
115
116 histo->add1D("16",
117 "Energy (MeV) in abs2",nBinsEA,0.0,maxEnergyAbs,MeV);
118
119 histo->add1D("17",
120 "Energy (MeV) in abs3",nBinsEA,0.0,maxEnergyAbs,MeV);
121
122 histo->add1D("18",
123 "Energy (MeV) in abs4",nBinsEA,0.0,maxEnergyAbs,MeV);
124
125 histo->add1D("19",
126 "Number of vertex hits",20,-0.5,19.5,1.0);
127
128 histo->add1D("20",
129 "E1/E9 Ratio",nBinsED,0.0,1,1.0);
130
131 histo->add1D("21",
132 "E1/25 Ratio",nBinsED,0.0,1.0,1.0);
133
134 histo->add1D("22",
135 "E9/E25 Ratio",nBinsED,0.0,1.0,1.0);
136
137 if(nTuple) {
138 histo->addTuple( "100", "Dose deposite","float r, z, e" );
139 }
140}
141
142//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143
144void HistoManager::BeginOfRun()
145{
146 // initilise scoring
147 n_evt = 0;
148 n_elec = 0;
149 n_posit= 0;
150 n_gam = 0;
151 n_step = 0;
152 n_lowe = 0;
153
154 for(G4int i=0; i<6; i++) {
155 stat[i] = 0;
156 edep[i] = 0.0;
157 erms[i] = 0.0;
158 if(i < 3) {
159 edeptr[i] = 0.0;
160 ermstr[i] = 0.0;
161 }
162 }
163
164 histo->book();
165
166 if(verbose > 0) {
167 G4cout << "HistoManager: Histograms are booked and run has been started"
168 << G4endl;
169 }
170}
171
172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173
174void HistoManager::EndOfRun()
175{
176
177 G4cout << "HistoManager: End of run actions are started" << G4endl;
178 G4String nam[6] = {"1x1", "3x3", "5x5", "E1/E9 ", "E1/E25", "E9/E25"};
179
180 // average
181
182 G4cout<<"================================================================="<<G4endl;
183 G4double x = (G4double)n_evt;
184 if(n_evt > 0) x = 1.0/x;
185 G4int j;
186 for(j=0; j<nmax; j++) {
187
188 // total mean
189 edep[j] *= x;
190 G4double y = erms[j]*x - edep[j]*edep[j];
191 if(y < 0.0) y = 0.0;
192 erms[j] = std::sqrt(y);
193
194 // trancated mean
195 G4double xx = G4double(stat[j]);
196 if(xx > 0.0) xx = 1.0/xx;
197 edeptr[j] *= xx;
198 y = ermstr[j]*xx - edeptr[j]*edeptr[j];
199 if(y < 0.0) y = 0.0;
200 ermstr[j] = std::sqrt(y);
201 }
202 G4double xe = x*(G4double)n_elec;
203 G4double xg = x*(G4double)n_gam;
204 G4double xp = x*(G4double)n_posit;
205 G4double xs = x*n_step;
206
207 G4double f = 100.*std::sqrt(beamEnergy/GeV);
208
209 G4cout << "Number of events " << n_evt <<G4endl;
210 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
211 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
212 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
213 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
214
215 for(j=0; j<3; j++) {
216 G4double e = edeptr[j];
217 G4double s = ermstr[j];
218 G4double xx= G4double(stat[j]);
219 if(xx > 0.0) xx = 1.0/xx;
220 G4double r = s*std::sqrt(xx);
221 G4cout << std::setprecision(4) << "Edep " << nam[j] << " = " << e
222 << " +- " << r;
223 if(e > 0.0) G4cout << " res= " << f*s/e << " %";
224 G4cout << G4endl;
225 }
226 if(limittrue[0] != DBL_MAX || limittrue[1] != DBL_MAX || limittrue[2] != DBL_MAX) {
227 G4cout<<"=========== Mean values without trancating ====================="<<G4endl;
228 for(j=0; j<nmax; j++) {
229 G4double e = edep[j];
230 G4double s = erms[j];
231 G4double r = s*std::sqrt(x);
232 G4cout << std::setprecision(4) << "Edep " << nam[j] << " = " << e
233 << " +- " << r;
234 if(e > 0.0) G4cout << " res= " << f*s/e << " %";
235 G4cout << G4endl;
236 }
237 }
238 G4cout<<"=========== Ratios without trancating ==========================="<<G4endl;
239 for(j=3; j<6; j++) {
240 G4double e = edep[j];
241 G4double xx= G4double(stat[j]);
242 if(xx > 0.0) xx = 1.0/xx;
243 e *= xx;
244 G4double y = erms[j]*xx - e*e;
245 G4double r = 0.0;
246 if(y > 0.0) r = std::sqrt(y*xx);
247 G4cout << " " << nam[j] << " = " << e
248 << " +- " << r;
249 G4cout << G4endl;
250 }
251 G4cout << std::setprecision(4) << "Beam Energy " << beamEnergy/GeV
252 << " GeV" << G4endl;
253 if(n_lowe > 0) G4cout << "Number of events E/E0<0.8 " << n_lowe << G4endl;
254 G4cout<<"=================================================================="<<G4endl;
255 G4cout<<G4endl;
256
257 // normalise histograms
258 for(G4int i=0; i<nHisto; i++) {
259 histo->scale(i,x);
260 }
261
262 histo->save();
263
264 // Acceptance
265 EmAcceptance acc;
266 G4bool isStarted = false;
267 for (j=0; j<nmax; j++) {
268
269 G4double ltrue = limittrue[j];
270 if (ltrue < DBL_MAX) {
271 if (!isStarted) {
272 acc.BeginOfAcceptance("Crystal Calorimeter",n_evt);
273 isStarted = true;
274 }
275 G4double etrue = edeptrue[j];
276 G4double rtrue = rmstrue[j];
277 acc.EmAcceptanceGauss("Edep"+nam[j],n_evt,edeptr[j],etrue,rtrue,ltrue);
278 acc.EmAcceptanceGauss("Erms"+nam[j],n_evt,ermstr[j],rtrue,rtrue,2.0*ltrue);
279 }
280 }
281 if(isStarted) acc.EndOfAcceptance();
282}
283
284//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
285
286void HistoManager::BeginOfEvent()
287{
288 n_evt++;
289
290 Eabs1 = 0.0;
291 Eabs2 = 0.0;
292 Eabs3 = 0.0;
293 Eabs4 = 0.0;
294 Evertex.clear();
295 Nvertex.clear();
296 for (G4int i=0; i<25; i++) {
297 E[i] = 0.0;
298 }
299}
300
301//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
302
303void HistoManager::EndOfEvent()
304{
305 G4double e9 = 0.0;
306 G4double e25= 0.0;
307 for (G4int i=0; i<25; i++) {
308 E[i] /= beamEnergy;
309 e25 += E[i];
310 if( ( 6<=i && 8>=i) || (11<=i && 13>=i) || (16<=i && 18>=i)) e9 += E[i];
311 }
312
313 if(e25 < 0.8) {
314 n_lowe++;
315 G4cout << "### in the event# " << n_evt << " E25= " << e25 << G4endl;
316 }
317
318 // compute ratios
319 G4double e0 = E[12];
320 G4double e19 = 0.0;
321 G4double e125 = 0.0;
322 G4double e925 = 0.0;
323 if(e9 > 0.0) {
324 e19 = e0/e9;
325 e125 = e0/e25;
326 e925 = e9/e25;
327 edep[3] += e19;
328 erms[3] += e19*e19;
329 edep[4] += e125;
330 erms[4] += e125*e125;
331 edep[5] += e925;
332 erms[5] += e925*e925;
333 stat[3] += 1;
334 stat[4] += 1;
335 stat[5] += 1;
336 }
337
338 // fill histo
339 histo->fill(0,e0,1.0);
340 histo->fill(1,e9,1.0);
341 histo->fill(2,e25,1.0);
342 histo->fill(5,Eabs1,1.0);
343 histo->fill(6,Eabs2,1.0);
344 histo->fill(7,Eabs3,1.0);
345 histo->fill(8,Eabs4,1.0);
346 histo->fill(9,G4double(Nvertex.size()),1.0);
347 histo->fill(10,e19,1.0);
348 histo->fill(11,e125,1.0);
349 histo->fill(12,e925,1.0);
350
351 // compute sums
352 edep[0] += e0;
353 erms[0] += e0*e0;
354 edep[1] += e9;
355 erms[1] += e9*e9;
356 edep[2] += e25;
357 erms[2] += e25*e25;
358
359 // trancated mean
360 if(limittrue[0] == DBL_MAX || std::abs(e0-edeptrue[0])<rmstrue[0]*limittrue[0]) {
361 stat[0] += 1;
362 edeptr[0] += e0;
363 ermstr[0] += e0*e0;
364 }
365 if(limittrue[1] == DBL_MAX || std::abs(e9-edeptrue[1])<rmstrue[1]*limittrue[1]) {
366 stat[1] += 1;
367 edeptr[1] += e9;
368 ermstr[1] += e9*e9;
369 }
370 if(limittrue[2] == DBL_MAX || std::abs(e25-edeptrue[2])<rmstrue[2]*limittrue[2]) {
371 stat[2] += 1;
372 edeptr[2] += e25;
373 ermstr[2] += e25*e25;
374 }
375 if(nTuple) histo->addRow();
376
377}
378
379//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
380
381void HistoManager::ScoreNewTrack(const G4Track* aTrack)
382{
383 //Save primary parameters
384 ResetTrackLength();
385 const G4ParticleDefinition* particle = aTrack->GetDefinition();
386 const G4DynamicParticle* dynParticle = aTrack->GetDynamicParticle();
387 G4String name = particle->GetParticleName();
388 G4int pid = aTrack->GetParentID();
389 G4double kinE = dynParticle->GetKineticEnergy();
390 G4ThreeVector pos = aTrack->GetVertexPosition();
391
392 if(0 == pid) {
393
394 beamEnergy = kinE;
395 histo->fillTuple("TKIN", kinE/MeV);
396
397 G4double mass = 0.0;
398 if(particle) {
399 mass = particle->GetPDGMass();
400 histo->fillTuple("MASS", mass/MeV);
401 histo->fillTuple("CHAR",(particle->GetPDGCharge())/eplus);
402 G4double beta = 1.;
403 if(mass > 0.) {
404 G4double gamma = kinE/mass + 1.;
405 beta = std::sqrt(1. - 1./(gamma*gamma));
406 }
407 }
408
409 G4ThreeVector dir = dynParticle->GetMomentumDirection();
410 if(1 < verbose) {
411 G4cout << "TrackingAction: Primary kinE(MeV)= " << kinE/MeV
412 << "; m(MeV)= " << mass/MeV
413 << "; pos= " << pos << "; dir= " << dir << G4endl;
414 }
415
416 // delta-electron
417 } else if (0 < pid && "e-" == name) {
418 if(1 < verbose) {
419 G4cout << "TrackingAction: Secondary electron " << G4endl;
420 }
421 AddDeltaElectron(dynParticle);
422
423 } else if (0 < pid && "e+" == name) {
424 if(1 < verbose) {
425 G4cout << "TrackingAction: Secondary positron " << G4endl;
426 }
427 AddPositron(dynParticle);
428
429 } else if (0 < pid && "gamma" == name) {
430 if(1 < verbose) {
431 G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
432 << " E= " << kinE << G4endl;
433 }
434 AddPhoton(dynParticle);
435
436 }
437}
438
439//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
440
441void HistoManager::AddEnergy(G4double edep, G4int volIndex, G4int copyNo)
442{
443 if(1 < verbose) {
444 G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV
445 << "; volIdx= " << volIndex
446 << "; copyNo= " << copyNo
447 << G4endl;
448 }
449 if(0 == volIndex) {
450 E[copyNo] += edep;
451 } else if (1 == volIndex) {
452 Eabs1 += edep;
453 } else if (2 == volIndex) {
454 Eabs2 += edep;
455 } else if (3 == volIndex) {
456 Eabs3 += edep;
457 } else if (4 == volIndex) {
458 Eabs4 += edep;
459 } else if (5 == volIndex) {
460 G4int n = Nvertex.size();
461 G4bool newPad = true;
462 if (n > 0) {
463 for(G4int i=0; i<n; i++) {
464 if (copyNo == Nvertex[i]) {
465 newPad = false;
466 Evertex[i] += edep;
467 break;
468 }
469 }
470 }
471 if(newPad) {
472 Nvertex.push_back(copyNo);
473 Evertex.push_back(edep);
474 }
475 }
476}
477
478//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
479
480void HistoManager::AddDeltaElectron(const G4DynamicParticle* elec)
481{
482 G4double e = elec->GetKineticEnergy()/MeV;
483 if(e > 0.0) n_elec++;
484 histo->fill(3,e,1.0);
485}
486
487//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
488
489void HistoManager::AddPhoton(const G4DynamicParticle* ph)
490{
491 G4double e = ph->GetKineticEnergy()/MeV;
492 if(e > 0.0) n_gam++;
493 histo->fill(4,e,1.0);
494}
495
496//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
497
498void HistoManager::SetEdepAndRMS(G4int i, G4ThreeVector val)
499{
500 if(i<nmax && i>=0) {
501 if(val[0] > 0.0) edeptrue[i] = val[0];
502 if(val[1] > 0.0) rmstrue[i] = val[1];
503 if(val[2] > 0.0) limittrue[i] = val[2];
504 }
505}
506
507//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
508
Note: See TracBrowser for help on using the repository browser.