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

Last change on this file since 1230 was 1230, checked in by garnier, 14 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.