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

Last change on this file since 1337 was 1337, checked in by garnier, 15 years ago

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