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

Last change on this file since 812 was 807, checked in by garnier, 16 years ago

update
File size: 12.4 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 = 10;
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    "Energy deposit (MeV) in central crystal",nBinsED,0.0,beamEnergy,MeV);
100
101  histo->add1D("11",
102    "Energy deposit (MeV) in 3x3",nBinsED,0.0,beamEnergy,MeV);
103
104  histo->add1D("12",
105    "Energy deposit (MeV) in 5x5",nBinsED,0.0,beamEnergy,MeV);
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  if(nTuple) {
129    histo->addTuple( "100", "Dose deposite","float r, z, e" );
130  }
131}
132
133//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134
135void HistoManager::BeginOfRun()
136{
137  // initilise scoring
138  n_evt  = 0;
139  n_elec = 0;
140  n_posit= 0;
141  n_gam  = 0;
142  n_step = 0;
143
144  for(G4int i=0; i<nmax; i++) {
145    stat[i] = 0;
146    edep[i] = 0.0;
147    erms[i] = 0.0;
148    edeptr[i] = 0.0;
149    ermstr[i] = 0.0;
150  }
151
152  histo->book();
153
154  if(verbose > 0) {
155    G4cout << "HistoManager: Histograms are booked and run has been started"
156           << G4endl;
157  }
158}
159
160//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161
162void HistoManager::EndOfRun()
163{
164
165  G4cout << "HistoManager: End of run actions are started" << G4endl;
166  G4String nam[3] = {"1x1", "3x3", "5x5"};
167
168  // average
169
170  G4cout<<"================================================================="<<G4endl;
171  G4double x = (G4double)n_evt;
172  if(n_evt > 0) x = 1.0/x;
173  G4int j;
174  for(j=0; j<nmax; j++) {
175
176    // total mean
177    edep[j] *= x/beamEnergy;
178    G4double y = erms[j]*x/(beamEnergy*beamEnergy) - edep[j]*edep[j];
179    if(y < 0.0) y = 0.0;
180    erms[j] = std::sqrt(y);
181
182    // trancated mean
183    G4double xx = G4double(stat[j]);
184    if(xx > 0.0) xx = 1.0/xx;
185    edeptr[j] *= xx/beamEnergy;
186    y = ermstr[j]*xx/(beamEnergy*beamEnergy) - edeptr[j]*edeptr[j];
187    if(y < 0.0) y = 0.0;
188    ermstr[j] = std::sqrt(y);
189  }
190  G4double xe = x*(G4double)n_elec;
191  G4double xg = x*(G4double)n_gam;
192  G4double xp = x*(G4double)n_posit;
193  G4double xs = x*(G4double)n_step;
194
195  G4double f = 100.*std::sqrt(beamEnergy/GeV);
196
197  G4cout                         << "Number of events             " << n_evt <<G4endl;
198  G4cout << std::setprecision(4) << "Average number of e-         " << xe << G4endl;
199  G4cout << std::setprecision(4) << "Average number of gamma      " << xg << G4endl;
200  G4cout << std::setprecision(4) << "Average number of e+         " << xp << G4endl;
201  G4cout << std::setprecision(4) << "Average number of steps      " << xs << G4endl;
202 
203  for(j=0; j<nmax; j++) {
204    G4double e = edeptr[j];
205    G4double s = ermstr[j];
206    G4double r = s*std::sqrt(x);
207    G4cout << std::setprecision(4) << "Edep " << nam[j] << " =                   " << e
208           << " +- " << r;
209    if(e > 0.0) G4cout << "  res=  " << f*s/e << " %";
210    G4cout << G4endl;
211  }
212  if(limittrue[0] != DBL_MAX || limittrue[1] != DBL_MAX || limittrue[2] != DBL_MAX) {
213    G4cout<<"===========  Mean values without trancating ====================="<<G4endl;
214    for(j=0; j<nmax; j++) {
215      G4double e = edep[j];
216      G4double s = erms[j];
217      G4double r = s*std::sqrt(x);
218      G4cout << std::setprecision(4) << "Edep " << nam[j] << " =                   " << e
219             << " +- " << r;
220      if(e > 0.0) G4cout << "  res=  " << f*s/e << " %";
221      G4cout << G4endl;
222    }
223  }
224  G4cout<<"=================================================================="<<G4endl;
225  G4cout<<G4endl;
226
227  // normalise histograms
228  for(G4int i=0; i<nHisto; i++) {
229    histo->scale(i,x);
230  }
231
232  histo->save();
233
234  // Acceptance
235  EmAcceptance acc;
236  G4bool isStarted = false;
237  for (j=0; j<nmax; j++) {
238
239    G4double ltrue = limittrue[j];
240    if (ltrue < DBL_MAX) {
241      if (!isStarted) {
242        acc.BeginOfAcceptance("Crystal Calorimeter",n_evt);
243        isStarted = true;
244      }
245      G4double etrue = edeptrue[j];
246      G4double rtrue = rmstrue[j];
247      acc.EmAcceptanceGauss("Edep"+nam[j],n_evt,edeptr[j],etrue,rtrue,ltrue);
248      acc.EmAcceptanceGauss("Erms"+nam[j],n_evt,ermstr[j],rtrue,rtrue,2.0*ltrue);
249    }
250  }
251  if(isStarted) acc.EndOfAcceptance();
252}
253
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255
256void HistoManager::BeginOfEvent()
257{
258  n_evt++;
259
260  Eabs1  = 0.0;
261  Eabs2  = 0.0;
262  Eabs3  = 0.0;
263  Eabs4  = 0.0;
264  Evertex.clear();
265  Nvertex.clear();
266  for (int i=0; i<25; i++) {
267    E[i] = 0.0;
268  }
269}
270
271//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
272
273void HistoManager::EndOfEvent()
274{
275  G4double e9 = 0.0;
276  G4double e25= 0.0;
277  for (int i=0; i<25; i++) {
278    e25 += E[i];
279    if( ( 6<=i &&  8>=i) || (11<=i && 13>=i) || (16<=i && 18>=i)) e9 += E[i];
280  }
281  histo->fill(0,E[12],1.0);
282  histo->fill(1,e9,1.0);
283  histo->fill(2,e25,1.0);
284  histo->fill(5,Eabs1,1.0);
285  histo->fill(6,Eabs2,1.0);
286  histo->fill(7,Eabs3,1.0);
287  histo->fill(8,Eabs4,1.0);
288  float nn = (double)(Nvertex.size());
289  histo->fill(9,nn,1.0);
290
291  G4double e0 = E[12];
292
293  edep[0] += e0;
294  erms[0] += e0*e0;
295  edep[1] += e9;
296  erms[1] += e9*e9;
297  edep[2] += e25;
298  erms[2] += e25*e25;
299
300  // trancated mean
301  if(limittrue[0] == DBL_MAX || std::abs(e0/beamEnergy-edeptrue[0])<rmstrue[0]*limittrue[0]) {
302    stat[0] += 1;
303    edeptr[0] += e0;
304    ermstr[0] += e0*e0;
305  }
306  if(limittrue[1] == DBL_MAX || std::abs(e9/beamEnergy-edeptrue[1])<rmstrue[1]*limittrue[1]) {
307    stat[1] += 1;
308    edeptr[1] += e9;
309    ermstr[1] += e9*e9;
310  }
311  if(limittrue[2] == DBL_MAX || std::abs(e25/beamEnergy-edeptrue[2])<rmstrue[2]*limittrue[2]) {
312    stat[2] += 1;
313    edeptr[2] += e25;
314    ermstr[2] += e25*e25;
315  }
316  if(nTuple) histo->addRow();
317
318}
319
320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
321
322void HistoManager::ScoreNewTrack(const G4Track* aTrack)
323{
324  //Save primary parameters
325  ResetTrackLength();
326  const G4ParticleDefinition* particle = aTrack->GetDefinition();
327  const G4DynamicParticle* dynParticle = aTrack->GetDynamicParticle();
328  G4String name = particle->GetParticleName();
329  G4int pid = aTrack->GetParentID();
330  G4double kinE = dynParticle->GetKineticEnergy();
331  G4ThreeVector pos = aTrack->GetVertexPosition();
332
333  if(0 == pid) {
334
335    histo->fillTuple("TKIN", kinE/MeV);
336
337    G4double mass = 0.0;
338    if(particle) {
339      mass = particle->GetPDGMass();
340      histo->fillTuple("MASS", mass/MeV);
341      histo->fillTuple("CHAR",(particle->GetPDGCharge())/eplus);
342      G4double beta = 1.;
343        if(mass > 0.) {
344          G4double gamma = kinE/mass + 1.;
345          beta = std::sqrt(1. - 1./(gamma*gamma));
346        }
347    }
348
349    G4ThreeVector dir = dynParticle->GetMomentumDirection();
350    if(1 < verbose) {
351      G4cout << "TrackingAction: Primary kinE(MeV)= " << kinE/MeV
352           << "; m(MeV)= " << mass/MeV
353           << "; pos= " << pos << ";  dir= " << dir << G4endl;
354    }
355
356    // delta-electron
357  } else if (0 < pid && "e-" == name) {
358    if(1 < verbose) {
359      G4cout << "TrackingAction: Secondary electron " << G4endl;
360    }
361    AddDeltaElectron(dynParticle);
362
363  } else if (0 < pid && "e+" == name) {
364    if(1 < verbose) {
365      G4cout << "TrackingAction: Secondary positron " << G4endl;
366    }
367    AddPositron(dynParticle);
368
369  } else if (0 < pid && "gamma" == name) {
370    if(1 < verbose) {
371      G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
372             << " E= " << kinE << G4endl;
373    }
374    AddPhoton(dynParticle);
375
376  }
377}
378
379//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
380
381void HistoManager::AddEnergy(G4double edep, G4int volIndex, G4int copyNo)
382{
383  if(1 < verbose) {
384    G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV
385           << "; volIdx= " << volIndex
386           << "; copyNo= " << copyNo
387           << G4endl;
388  }
389  if(0 == volIndex) {
390    E[copyNo] += edep;
391  } else if (1 == volIndex) {
392    Eabs1 += edep;
393  } else if (2 == volIndex) {
394    Eabs2 += edep;
395  } else if (3 == volIndex) {
396    Eabs3 += edep;
397  } else if (4 == volIndex) {
398    Eabs4 += edep;
399  } else if (5 == volIndex) {
400    G4int n = Nvertex.size();
401    G4bool newPad = true;
402    if (n > 0) {
403      for(G4int i=0; i<n; i++) {
404        if (copyNo == Nvertex[i]) {
405          newPad = false;
406          Evertex[i] += edep;
407          break;
408        }
409      }
410    }
411    if(newPad) {
412      Nvertex.push_back(copyNo);
413      Evertex.push_back(edep);
414    }
415  }
416}
417
418//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
419
420void HistoManager::AddDeltaElectron(const G4DynamicParticle* elec)
421{
422  G4double e = elec->GetKineticEnergy()/MeV;
423  if(e > 0.0) n_elec++;
424  histo->fill(3,e,1.0);
425}
426
427//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
428
429void HistoManager::AddPhoton(const G4DynamicParticle* ph)
430{
431  G4double e = ph->GetKineticEnergy()/MeV;
432  if(e > 0.0) n_gam++;
433  histo->fill(4,e,1.0);
434}
435
436//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
437
438void HistoManager::SetEdepAndRMS(G4int i, G4ThreeVector val)
439{
440  if(i<nmax && i>=0) {
441    if(val[0] > 0.0) edeptrue[i] = val[0];
442    if(val[1] > 0.0) rmstrue[i] = val[1];
443    if(val[2] > 0.0) limittrue[i] = val[2];
444  }
445}
446
447//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
448
Note: See TracBrowser for help on using the repository browser.