Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Histo.cc
Go to the documentation of this file.
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 //
28 //
29 // $Id$
30 //
31 //---------------------------------------------------------------------------
32 //
33 // ClassName: Histo
34 //
35 //
36 // Author: V.Ivanchenko 30/01/01
37 //
38 //----------------------------------------------------------------------------
39 //
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43 
44 #include "Histo.hh"
45 #include "G4Gamma.hh"
46 #include "G4Electron.hh"
47 #include "G4Positron.hh"
48 #include "G4RootAnalysisManager.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include <iomanip>
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 
55 Histo* Histo::fManager = 0;
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58 
59 Histo* Histo::GetPointer()
60 {
61  if(!fManager) {
62  static Histo manager;
63  fManager = &manager;
64  }
65  return fManager;
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 
70 Histo::Histo()
71 {
72  fVerbose = 1;
73  fHistName = "";
74  fNHisto = 10;
75  fMaxEnergy = 50.0*MeV;
76  fNBinsZ = 60;
77  fNBinsR = 80;
78  fNBinsE = 200;
79  fScoreBin = 0;
80 
81  fAbsorberZ = 300.*mm;
82  fAbsorberR = 200.*mm;
83  fScoreZ = 100.*mm;
84 
85  fGamma = G4Gamma::Gamma();
86  fElectron = G4Electron::Electron();
87  fPositron = G4Positron::Positron();
88 
89  fCheckVolume = fGasVolume = fPhantom = fTarget1 = fTarget2 = 0;
90 
91  fAnalysisManager = 0;
92 
93  fStepZ = fStepR = fStepE = fNormZ = fSumR = 0.0;
94  fNevt = fNelec = fNposit = fNgam = fNstep = fNgamPh =
95  fNgamTar = fNeTar = fNePh = fNstepTarget = 0;
96 
97  fHisto.resize(fNHisto, 0);
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
102 Histo::~Histo()
103 {}
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106 
107 void Histo::BeginOfHisto()
108 {
109  G4cout << "### Histo start initialisation nHisto= " << fNHisto << G4endl;
110 
111  fNevt = fNelec = fNposit= fNgam = fNstep = fNgamPh = fNgamTar =
112  fNeTar = fNePh = fNstepTarget = 0;
113  fSumR = 0.0;
114  if(fNBinsR>1000) { SetNumberDivR(40); }
115 
116  fStepZ = fAbsorberZ/(G4double)fNBinsZ;
117  fStepR = fAbsorberR/(G4double)fNBinsR;
118  fStepE = fMaxEnergy/(G4double)fNBinsE;
119  fScoreBin = (G4int)(fScoreZ/fStepZ + 0.5);
120 
121  G4cout << " "<< fNBinsR << " bins R stepR= " << fStepR/mm << " mm "
122  << G4endl;
123  G4cout << " "<< fNBinsZ << " bins Z stepZ= " << fStepZ/mm << " mm "
124  << G4endl;
125  G4cout << " "<< fNBinsE << " bins E stepE= " << fStepE/MeV << " MeV "
126  << G4endl;
127  G4cout << " "<< fScoreBin << "-th bin in Z is used for R distribution"
128  << G4endl;
129 
130  fVolumeR.resize(fNBinsR);
131  fEdep.resize(fNBinsR, 0.0);
132  fGammaE.resize(fNBinsE, 0.0);
133 
134  G4double r1 = 0.0;
135  G4double r2 = fStepR;
136  for(G4int i=0; i<fNBinsR; ++i) {
137  fVolumeR[i] = cm*cm/(pi*(r2*r2 - r1*r1));
138  r1 = r2;
139  r2 += fStepR;
140  }
141 
142  if(fHistName != "") {
143  if(!fAnalysisManager) {
144  fAnalysisManager = G4RootAnalysisManager::Instance();
145  BookHisto();
146  }
147  if(fVerbose > 0) {
148  G4cout << "Histo: Histograms are booked and run has been started"
149  << G4endl;
150  }
151  } else if(fVerbose > 0) {
152  G4cout << "Histo: Histograms are not booked because file name is not set"
153  << G4endl;
154  }
155 }
156 
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158 
159 void Histo::EndOfHisto()
160 {
161 
162  G4cout << "Histo: End of run actions are started" << G4endl;
163 
164  // average
165 
166  G4cout<<"========================================================"<<G4endl;
167  G4double x = (G4double)fNevt;
168  if(fNevt > 0) { x = 1.0/x; }
169  G4double xe = x*(G4double)fNelec;
170  G4double xg = x*(G4double)fNgam;
171  G4double xp = x*(G4double)fNposit;
172  G4double xs = x*(G4double)fNstep;
173  G4double xph= x*(G4double)fNgamPh;
174  G4double xes= x*(G4double)fNstepTarget;
175  G4double xgt= x*(G4double)fNgamTar;
176  G4double xet= x*(G4double)fNeTar;
177  G4double xphe= x*(G4double)fNePh;
178 
179  G4cout << "Number of events "
180  << std::setprecision(8) << fNevt <<G4endl;
181  G4cout << std::setprecision(4) << "Average number of e- "
182  << xe << G4endl;
183  G4cout << std::setprecision(4) << "Average number of gamma "
184  << xg << G4endl;
185  G4cout << std::setprecision(4) << "Average number of e+ "
186  << xp << G4endl;
187  G4cout << std::setprecision(4) << "Average number of steps in the phantom "
188  << xs << G4endl;
189  G4cout << std::setprecision(4) << "Average number of e- steps in the target "
190  << xes << G4endl;
191  G4cout << std::setprecision(4) << "Average number of g produced in the target "
192  << xgt << G4endl;
193  G4cout << std::setprecision(4) << "Average number of e- produced in the target "
194  << xet << G4endl;
195  G4cout << std::setprecision(4) << "Average number of g produced in the phantom "
196  << xph << G4endl;
197  G4cout << std::setprecision(4) << "Average number of e- produced in the phantom "
198  << xphe << G4endl;
199  G4cout << std::setprecision(4) << "Total fGamma fluence in front of the phantom "
200  << x*fSumR/MeV
201  << " MeV " << G4endl;
202  G4cout<<"========================================================"<<G4endl;
203  G4cout<<G4endl;
204  G4cout<<G4endl;
205 
206  G4double sum = 0.0;
207  for(G4int i=0; i<fNBinsR; i++) {
208  fEdep[i] *= x;
209  sum += fEdep[i];
210  }
211 
212  if(fAnalysisManager) {
213 
214  // normalise histograms
215  for(G4int i=0; i<fNHisto; i++) {
216  fAnalysisManager->GetH1(fHisto[i])->scale(x);
217  }
218  G4double nr = fEdep[0];
219  if(nr > 0.0) { nr = 1./nr; }
220  fAnalysisManager->GetH1(fHisto[0])->scale(nr);
221 
222  nr = sum*fStepR;
223  if(nr > 0.0) { nr = 1./nr; }
224  fAnalysisManager->GetH1(fHisto[1])->scale(nr);
225 
226  fAnalysisManager->GetH1(fHisto[3])
227  ->scale(1000.0*cm3/(pi*fAbsorberR*fAbsorberR*fStepZ));
228  fAnalysisManager->GetH1(fHisto[4])
229  ->scale(1000.0*cm3*fVolumeR[0]/fStepZ);
230 
231  // Write histogram file
232  if(!fAnalysisManager->Write()) {
233  G4cout << "Histo::Save: FATAL ERROR writing ROOT file" << G4endl;
234  exit(1);
235  }
236  G4cout << "### Histo::Save: Histograms are saved" << G4endl;
237  if(fAnalysisManager->CloseFile() && fVerbose > 0) {
238  G4cout << " File is closed" << G4endl;
239  }
240  delete fAnalysisManager;
241  fAnalysisManager = 0;
242  }
243 }
244 
245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
246 
247 void Histo::BookHisto()
248 {
249  G4String nam = fHistName + ".root";
250 
251  if(!fAnalysisManager->OpenFile(nam)) {
252  G4cout << "Histo::bookHisto: ERROR open file <" << nam << ">" << G4endl;
253  exit(0);
254  }
255  G4cout << "### Histo::Save: Opended file <" << nam << "> for "
256  << fNHisto << " histograms " << G4endl;
257 
258  // Creating an 1-dimensional histograms in the root directory of the tree
259  fHisto[0] = fAnalysisManager->CreateH1("10",
260  "Energy deposit at radius (mm) normalised on 1st channel",fNBinsR,0.,fAbsorberR/mm);
261 
262  fHisto[1] = fAnalysisManager->CreateH1("11",
263  "Energy deposit at radius (mm) normalised to integral",fNBinsR,0.,fAbsorberR/mm);
264 
265  fHisto[2] = fAnalysisManager->CreateH1("12",
266  "Energy deposit (MeV/kg/electron) at radius (mm)",fNBinsR,0.,fAbsorberR/mm);
267 
268  fHisto[3] = fAnalysisManager->CreateH1("13",
269  "Energy profile (MeV/kg/electron) over Z (mm)",fNBinsZ,0.,fAbsorberZ/mm);
270 
271  fHisto[4] = fAnalysisManager->CreateH1("14",
272  "Energy profile (MeV/kg/electron) over Z (mm) at Central Voxel",fNBinsZ,0.,fAbsorberZ/mm);
273 
274  fHisto[5] = fAnalysisManager->CreateH1("15",
275  "Energy (MeV) of fGamma produced in the target",fNBinsE,0.,fMaxEnergy/MeV);
276 
277  fHisto[6] = fAnalysisManager->CreateH1("16",
278  "Energy (MeV) of fGamma before phantom",fNBinsE,0.,fMaxEnergy/MeV);
279 
280  fHisto[7] = fAnalysisManager->CreateH1("17",
281  "Energy (MeV) of electrons produced in phantom",fNBinsE,0.,fMaxEnergy/MeV);
282 
283  fHisto[8] = fAnalysisManager->CreateH1("18",
284  "Energy (MeV) of electrons produced in target",fNBinsE,0.,fMaxEnergy/MeV);
285 
286  fHisto[9] = fAnalysisManager->CreateH1("19",
287  "Gamma Energy Fluence (MeV/cm2) at radius(mm) in front of phantom",fNBinsR,0.,fAbsorberR/mm);
288 
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
292 
293 void Histo::AddTargetPhoton(const G4DynamicParticle* p)
294 {
295  ++fNgamTar;
296  if(fAnalysisManager) {
297  fAnalysisManager->FillH1(fHisto[5],p->GetKineticEnergy()/MeV,1.0);
298  }
299 }
300 
301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
302 
303 void Histo::AddPhantomPhoton(const G4DynamicParticle*)
304 {
305  ++fNgamPh;
306 }
307 
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
309 
310 void Histo::AddTargetElectron(const G4DynamicParticle* p)
311 {
312  ++fNeTar;
313  if(fAnalysisManager) {
314  fAnalysisManager->FillH1(fHisto[8],p->GetKineticEnergy()/MeV,1.0);
315  }
316 }
317 
318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
319 
320 void Histo::AddPhantomElectron(const G4DynamicParticle* p)
321 {
322  ++fNePh;
323  if(fAnalysisManager) {
324  fAnalysisManager->FillH1(fHisto[7],p->GetKineticEnergy()/MeV,1.0);
325  }
326 }
327 
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329 
330 void Histo::ScoreNewTrack(const G4Track* aTrack)
331 {
332  //Save primary parameters
333  const G4ParticleDefinition* particle = aTrack->GetParticleDefinition();
334  G4int pid = aTrack->GetParentID();
335  G4VPhysicalVolume* pv = aTrack->GetVolume();
336  const G4DynamicParticle* dp = aTrack->GetDynamicParticle();
337 
338  //primary particle
339  if(0 == pid) {
340 
341  ++fNevt;
342  if(1 < fVerbose) {
343  G4ThreeVector pos = aTrack->GetVertexPosition();
344  G4ThreeVector dir = aTrack->GetMomentumDirection();
345  G4cout << "TrackingAction: Primary "
346  << particle->GetParticleName()
347  << " Ekin(MeV)= "
348  << aTrack->GetKineticEnergy()/MeV
349  << "; pos= " << pos << "; dir= " << dir << G4endl;
350  }
351 
352  // secondary electron
353  } else if (0 < pid && particle == fElectron) {
354  if(1 < fVerbose) {
355  G4cout << "TrackingAction: Secondary electron " << G4endl;
356  }
357  AddElectron(dp);
358  if(pv == fPhantom) { AddPhantomElectron(dp); }
359  else if(pv == fTarget1 || pv == fTarget2) { AddTargetElectron(dp); }
360 
361  // secondary positron
362  } else if (0 < pid && particle == fPositron) {
363  if(1 < fVerbose) {
364  G4cout << "TrackingAction: Secondary positron " << G4endl;
365  }
366  AddPositron(dp);
367 
368  // secondary gamma
369  } else if (0 < pid && particle == fGamma) {
370  if(1 < fVerbose) {
371  G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
372  << " E= " << aTrack->GetKineticEnergy() << G4endl;
373  }
374  AddPhoton(dp);
375  if(pv == fPhantom) { AddPhantomPhoton(dp); }
376  else if(pv == fTarget1 || pv == fTarget2) { AddTargetPhoton(dp); }
377 
378  }
379 }
380 
381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
382 
383 void Histo::AddPhantomGamma(G4double e, G4double r)
384 {
385  e /= MeV;
386  fSumR += e;
387  G4int bin = (G4int)(e/fStepE);
388  if(bin >= fNBinsE) { bin = fNBinsE-1; }
389  fGammaE[bin] += e;
390  G4int bin1 = (G4int)(r/fStepR);
391  if(bin1 >= fNBinsR) { bin1 = fNBinsR-1; }
392  if(fAnalysisManager) {
393  fAnalysisManager->FillH1(fHisto[6],e,1.0);
394  fAnalysisManager->FillH1(fHisto[9],r,e*fVolumeR[bin1]);
395  }
396 }
397 
398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
399 
400 void Histo::AddPhantomStep(G4double edep, G4double r1, G4double z1,
401  G4double r2, G4double z2,
402  G4double r0, G4double z0)
403 {
404  ++fNstep;
405  G4int nzbin = (G4int)(z0/fStepZ);
406  if(fVerbose > 1) {
407  G4cout << "Histo: edep(MeV)= " << edep/MeV << " at binz= " << nzbin
408  << " r1= " << r1 << " z1= " << z1
409  << " r2= " << r2 << " z2= " << z2
410  << " r0= " << r0 << " z0= " << z0
411  << G4endl;
412  }
413  if(nzbin == fScoreBin) {
414  G4int bin = (G4int)(r0/fStepR);
415  if(bin >= fNBinsR) { bin = fNBinsR-1; }
416  double w = edep*fVolumeR[bin];
417  fEdep[bin] += w;
418  if(fAnalysisManager) {
419  fAnalysisManager->FillH1(fHisto[0],r0,w);
420  fAnalysisManager->FillH1(fHisto[1],r0,w);
421  fAnalysisManager->FillH1(fHisto[2],r0,w);
422  }
423  }
424  G4int bin1 = (G4int)(z1/fStepZ);
425  if(bin1 >= fNBinsZ) { bin1 = fNBinsZ-1; }
426  G4int bin2 = (G4int)(z2/fStepZ);
427  if(bin2 >= fNBinsZ) { bin2 = fNBinsZ-1; }
428  if(bin1 == bin2) {
429  if(fAnalysisManager) {
430  fAnalysisManager->FillH1(fHisto[3],z0,edep);
431  if(r1 < fStepR) {
432  G4double w = edep;
433  if(r2 > fStepR) { w *= (fStepR - r1)/(r2 - r1); }
434  fAnalysisManager->FillH1(fHisto[4],z0,w);
435  }
436  }
437  } else {
438  G4int bin;
439 
440  if(bin2 < bin1) {
441  bin = bin2;
442  G4double z = z2;
443  bin2 = bin1;
444  z2 = z1;
445  bin1 = bin;
446  z1 = z;
447  }
448  G4double zz1 = z1;
449  G4double zz2 = (bin1+1)*fStepZ;
450  G4double rr1 = r1;
451  G4double dz = z2 - z1;
452  G4double dr = r2 - r1;
453  G4double rr2 = r1 + dr*(zz2-zz1)/dz;
454  for(bin=bin1; bin<=bin2; bin++) {
455  if(fAnalysisManager) {
456  G4double de = edep*(zz2 - zz1)/dz;
457  G4double zf = (zz1+zz2)*0.5;
458  { fAnalysisManager->FillH1(fHisto[3],zf,de); }
459  if(rr1 < fStepR) {
460  G4double w = de;
461  if(rr2 > fStepR) w *= (fStepR - rr1)/(rr2 - rr1);
462  { fAnalysisManager->FillH1(fHisto[4],zf,w); }
463  }
464  }
465  zz1 = zz2;
466  zz2 = std::min(z2, zz1+fStepZ);
467  rr1 = rr2;
468  rr2 = rr1 + dr*(zz2 - zz1)/dz;
469  }
470  }
471 }
472 
473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
474 
Generated on Sat May 25 2013 14:32:21 for Geant4 by   doxygen 1.8.4