Geant4  10.03.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 //
26 // $Id: Histo.cc 82277 2014-06-13 14:40:54Z gcosmo $
27 //
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
182  << std::setprecision(4) << "Average number of e- "
183  << xe << G4endl;
184  G4cout
185  << std::setprecision(4) << "Average number of gamma "
186  << xg << G4endl;
187  G4cout
188  << std::setprecision(4) << "Average number of e+ "
189  << xp << G4endl;
190  G4cout
191  << std::setprecision(4) << "Average number of steps in the phantom "
192  << xs << G4endl;
193  G4cout
194  << std::setprecision(4) << "Average number of e- steps in the target "
195  << xes << G4endl;
196  G4cout
197  << std::setprecision(4) << "Average number of g produced in the target "
198  << xgt << G4endl;
199  G4cout
200  << std::setprecision(4) << "Average number of e- produced in the target "
201  << xet << G4endl;
202  G4cout
203  << std::setprecision(4) << "Average number of g produced in the phantom "
204  << xph << G4endl;
205  G4cout
206  << std::setprecision(4) << "Average number of e- produced in the phantom "
207  << xphe << G4endl;
208  G4cout
209  << std::setprecision(4) << "Total fGamma fluence in front of the phantom "
210  << x*fSumR/MeV << " MeV " << G4endl;
211  G4cout<<"========================================================"<<G4endl;
212  G4cout<<G4endl;
213  G4cout<<G4endl;
214 
215  G4double sum = 0.0;
216  for(G4int i=0; i<fNBinsR; i++) {
217  fEdep[i] *= x;
218  sum += fEdep[i];
219  }
220 
221  if(fAnalysisManager) {
222 
223  // normalise histograms
224  for(G4int i=0; i<fNHisto; i++) {
225  fAnalysisManager->GetH1(fHisto[i])->scale(x);
226  }
227  G4double nr = fEdep[0];
228  if(nr > 0.0) { nr = 1./nr; }
229  fAnalysisManager->GetH1(fHisto[0])->scale(nr);
230 
231  nr = sum*fStepR;
232  if(nr > 0.0) { nr = 1./nr; }
233  fAnalysisManager->GetH1(fHisto[1])->scale(nr);
234 
235  fAnalysisManager->GetH1(fHisto[3])
236  ->scale(1000.0*cm3/(pi*fAbsorberR*fAbsorberR*fStepZ));
237  fAnalysisManager->GetH1(fHisto[4])
238  ->scale(1000.0*cm3*fVolumeR[0]/fStepZ);
239 
240  // Write histogram file
241  if(!fAnalysisManager->Write()) {
242  G4Exception ("Histo::Save()", "hist01", FatalException,
243  "Cannot write ROOT file.");
244  }
245  G4cout << "### Histo::Save: Histograms are saved" << G4endl;
246  if(fAnalysisManager->CloseFile() && fVerbose > 0) {
247  G4cout << " File is closed" << G4endl;
248  }
249  delete fAnalysisManager;
250  fAnalysisManager = 0;
251  }
252 }
253 
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255 
256 void Histo::BookHisto()
257 {
258  G4String nam = fHistName + ".root";
259 
260  if(!fAnalysisManager->OpenFile(nam)) {
261  G4cout << "Histo::bookHisto: ERROR open file <" << nam << ">" << G4endl;
262  G4Exception ("Histo::BookHisto()", "hist01", FatalException,
263  "Cannot open ROOT file.");
264  }
265  G4cout << "### Histo::Save: Opended file <" << nam << "> for "
266  << fNHisto << " histograms " << G4endl;
267 
268  // Creating an 1-dimensional histograms in the root directory of the tree
269  fHisto[0] = fAnalysisManager->CreateH1("10",
270  "Energy deposit at radius (mm) normalised on 1st channel",
271  fNBinsR, 0., fAbsorberR/mm);
272 
273  fHisto[1] = fAnalysisManager->CreateH1("11",
274  "Energy deposit at radius (mm) normalised to integral",
275  fNBinsR, 0., fAbsorberR/mm);
276 
277  fHisto[2] = fAnalysisManager->CreateH1("12",
278  "Energy deposit (MeV/kg/electron) at radius (mm)",
279  fNBinsR, 0., fAbsorberR/mm);
280 
281  fHisto[3] = fAnalysisManager->CreateH1("13",
282  "Energy profile (MeV/kg/electron) over Z (mm)",fNBinsZ,0.,fAbsorberZ/mm);
283 
284  fHisto[4] = fAnalysisManager->CreateH1("14",
285  "Energy profile (MeV/kg/electron) over Z (mm) at Central Voxel",
286  fNBinsZ, 0., fAbsorberZ/mm);
287 
288  fHisto[5] = fAnalysisManager->CreateH1("15",
289  "Energy (MeV) of fGamma produced in the target",
290  fNBinsE, 0., fMaxEnergy/MeV);
291 
292  fHisto[6] = fAnalysisManager->CreateH1("16",
293  "Energy (MeV) of fGamma before phantom",fNBinsE,0.,fMaxEnergy/MeV);
294 
295  fHisto[7] = fAnalysisManager->CreateH1("17",
296  "Energy (MeV) of electrons produced in phantom",fNBinsE,0.,fMaxEnergy/MeV);
297 
298  fHisto[8] = fAnalysisManager->CreateH1("18",
299  "Energy (MeV) of electrons produced in target",fNBinsE,0.,fMaxEnergy/MeV);
300 
301  fHisto[9] = fAnalysisManager->CreateH1("19",
302  "Gamma Energy Fluence (MeV/cm2) at radius(mm) in front of phantom",
303  fNBinsR, 0., fAbsorberR/mm);
304 
305 }
306 
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 
309 void Histo::AddTargetPhoton(const G4DynamicParticle* p)
310 {
311  ++fNgamTar;
312  if(fAnalysisManager) {
313  fAnalysisManager->FillH1(fHisto[5],p->GetKineticEnergy()/MeV,1.0);
314  }
315 }
316 
317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
318 
319 void Histo::AddPhantomPhoton(const G4DynamicParticle*)
320 {
321  ++fNgamPh;
322 }
323 
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
325 
326 void Histo::AddTargetElectron(const G4DynamicParticle* p)
327 {
328  ++fNeTar;
329  if(fAnalysisManager) {
330  fAnalysisManager->FillH1(fHisto[8],p->GetKineticEnergy()/MeV,1.0);
331  }
332 }
333 
334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
335 
336 void Histo::AddPhantomElectron(const G4DynamicParticle* p)
337 {
338  ++fNePh;
339  if(fAnalysisManager) {
340  fAnalysisManager->FillH1(fHisto[7],p->GetKineticEnergy()/MeV,1.0);
341  }
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
345 
346 void Histo::ScoreNewTrack(const G4Track* aTrack)
347 {
348  //Save primary parameters
349  const G4ParticleDefinition* particle = aTrack->GetParticleDefinition();
350  G4int pid = aTrack->GetParentID();
351  G4VPhysicalVolume* pv = aTrack->GetVolume();
352  const G4DynamicParticle* dp = aTrack->GetDynamicParticle();
353 
354  //primary particle
355  if(0 == pid) {
356 
357  ++fNevt;
358  if(1 < fVerbose) {
359  G4ThreeVector pos = aTrack->GetVertexPosition();
360  G4ThreeVector dir = aTrack->GetMomentumDirection();
361  G4cout << "TrackingAction: Primary "
362  << particle->GetParticleName()
363  << " Ekin(MeV)= "
364  << aTrack->GetKineticEnergy()/MeV
365  << "; pos= " << pos << "; dir= " << dir << G4endl;
366  }
367 
368  // secondary electron
369  } else if (0 < pid && particle == fElectron) {
370  if(1 < fVerbose) {
371  G4cout << "TrackingAction: Secondary electron " << G4endl;
372  }
373  AddElectron(dp);
374  if(pv == fPhantom) { AddPhantomElectron(dp); }
375  else if(pv == fTarget1 || pv == fTarget2) { AddTargetElectron(dp); }
376 
377  // secondary positron
378  } else if (0 < pid && particle == fPositron) {
379  if(1 < fVerbose) {
380  G4cout << "TrackingAction: Secondary positron " << G4endl;
381  }
382  AddPositron(dp);
383 
384  // secondary gamma
385  } else if (0 < pid && particle == fGamma) {
386  if(1 < fVerbose) {
387  G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
388  << " E= " << aTrack->GetKineticEnergy() << G4endl;
389  }
390  AddPhoton(dp);
391  if(pv == fPhantom) { AddPhantomPhoton(dp); }
392  else if(pv == fTarget1 || pv == fTarget2) { AddTargetPhoton(dp); }
393 
394  }
395 }
396 
397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
398 
399 void Histo::AddPhantomGamma(G4double e, G4double r)
400 {
401  e /= MeV;
402  fSumR += e;
403  G4int bin = (G4int)(e/fStepE);
404  if(bin >= fNBinsE) { bin = fNBinsE-1; }
405  fGammaE[bin] += e;
406  G4int bin1 = (G4int)(r/fStepR);
407  if(bin1 >= fNBinsR) { bin1 = fNBinsR-1; }
408  if(fAnalysisManager) {
409  fAnalysisManager->FillH1(fHisto[6],e,1.0);
410  fAnalysisManager->FillH1(fHisto[9],r,e*fVolumeR[bin1]);
411  }
412 }
413 
414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
415 
416 void Histo::AddPhantomStep(G4double edep, G4double r1, G4double z1,
417  G4double r2, G4double z2,
418  G4double r0, G4double z0)
419 {
420  ++fNstep;
421  G4int nzbin = (G4int)(z0/fStepZ);
422  if(fVerbose > 1) {
423  G4cout << "Histo: edep(MeV)= " << edep/MeV << " at binz= " << nzbin
424  << " r1= " << r1 << " z1= " << z1
425  << " r2= " << r2 << " z2= " << z2
426  << " r0= " << r0 << " z0= " << z0
427  << G4endl;
428  }
429  if(nzbin == fScoreBin) {
430  G4int bin = (G4int)(r0/fStepR);
431  if(bin >= fNBinsR) { bin = fNBinsR-1; }
432  double w = edep*fVolumeR[bin];
433  fEdep[bin] += w;
434  if(fAnalysisManager) {
435  fAnalysisManager->FillH1(fHisto[0],r0,w);
436  fAnalysisManager->FillH1(fHisto[1],r0,w);
437  fAnalysisManager->FillH1(fHisto[2],r0,w);
438  }
439  }
440  G4int bin1 = (G4int)(z1/fStepZ);
441  if(bin1 >= fNBinsZ) { bin1 = fNBinsZ-1; }
442  G4int bin2 = (G4int)(z2/fStepZ);
443  if(bin2 >= fNBinsZ) { bin2 = fNBinsZ-1; }
444  if(bin1 == bin2) {
445  if(fAnalysisManager) {
446  fAnalysisManager->FillH1(fHisto[3],z0,edep);
447  if(r1 < fStepR) {
448  G4double w = edep;
449  if(r2 > fStepR) { w *= (fStepR - r1)/(r2 - r1); }
450  fAnalysisManager->FillH1(fHisto[4],z0,w);
451  }
452  }
453  } else {
454  G4int bin;
455 
456  if(bin2 < bin1) {
457  bin = bin2;
458  G4double z = z2;
459  bin2 = bin1;
460  z2 = z1;
461  bin1 = bin;
462  z1 = z;
463  }
464  G4double zz1 = z1;
465  G4double zz2 = (bin1+1)*fStepZ;
466  G4double rr1 = r1;
467  G4double dz = z2 - z1;
468  G4double dr = r2 - r1;
469  G4double rr2 = r1 + dr*(zz2-zz1)/dz;
470  for(bin=bin1; bin<=bin2; bin++) {
471  if(fAnalysisManager) {
472  G4double de = edep*(zz2 - zz1)/dz;
473  G4double zf = (zz1+zz2)*0.5;
474  { fAnalysisManager->FillH1(fHisto[3],zf,de); }
475  if(rr1 < fStepR) {
476  G4double w = de;
477  if(rr2 > fStepR) w *= (fStepR - rr1)/(rr2 - rr1);
478  { fAnalysisManager->FillH1(fHisto[4],zf,w); }
479  }
480  }
481  zz1 = zz2;
482  zz2 = std::min(z2, zz1+fStepZ);
483  rr1 = rr2;
484  rr2 = rr1 + dr*(zz2 - zz1)/dz;
485  }
486  }
487 }
488 
489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
490 
Histo::SetNumberDivR
void SetNumberDivR(G4int val)
Definition: Histo.hh:99
G4Track::GetParentID
G4int GetParentID() const
Histo::AddPhantomGamma
void AddPhantomGamma(G4double e, G4double r)
Definition: Histo.cc:399
Histo::GetPointer
static Histo * GetPointer()
Definition: Histo.cc:59
mm
static constexpr double mm
Definition: G4SIunits.hh:115
plottest35.bin
tuple bin
Definition: plottest35.py:22
Histo::~Histo
~Histo()
Definition: Histo.cc:70
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4VAnalysisManager::CreateH1
G4int CreateH1(const G4String &name, const G4String &title, G4int nbins, G4double xmin, G4double xmax, const G4String &unitName="none", const G4String &fcnName="none", const G4String &binSchemeName="linear")
Definition: G4VAnalysisManager.cc:292
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
Histo
Definition: Histo.hh:56
Histo::Histo
Histo()
Definition: Histo.cc:51
p
const char * p
Definition: xmltok.h:285
test.x
tuple x
Definition: test.py:50
G4int
int G4int
Definition: G4Types.hh:78
G4VAnalysisManager::OpenFile
G4bool OpenFile(const G4String &fileName="")
Definition: G4VAnalysisManager.cc:185
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
Histo::ScoreNewTrack
void ScoreNewTrack(const G4Track *aTrack)
Definition: Histo.cc:346
G4Track::GetKineticEnergy
G4double GetKineticEnergy() const
G4cout
G4GLOB_DLL std::ostream G4cout
cm
static constexpr double cm
Definition: G4SIunits.hh:119
G4Positron.hh
G4Track::GetParticleDefinition
const G4ParticleDefinition * GetParticleDefinition() const
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4RootAnalysisManager.hh
cm3
static constexpr double cm3
Definition: G4SIunits.hh:121
G4Track::GetVertexPosition
const G4ThreeVector & GetVertexPosition() const
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
Histo::BeginOfHisto
void BeginOfHisto()
Definition: Histo.cc:107
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Electron.hh
G4Track::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4VAnalysisManager::FillH1
G4bool FillH1(G4int id, G4double value, G4double weight=1.0)
FatalException
Definition: G4ExceptionSeverity.hh:60
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4RootAnalysisManager::Instance
static G4RootAnalysisManager * Instance()
Definition: G4RootAnalysisManager.cc:63
G4Track::GetVolume
G4VPhysicalVolume * GetVolume() const
test.z
tuple z
Definition: test.py:28
Histo::EndOfHisto
void EndOfHisto()
Definition: Histo.cc:159
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
pi
static constexpr double pi
Definition: G4SIunits.hh:75
Histo::AddPhantomStep
void AddPhantomStep(G4double e, G4double r1, G4double z1, G4double r2, G4double z2, G4double r0, G4double z0)
Definition: Histo.cc:416
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4ToolsAnalysisManager::GetH1
tools::histo::h1d * GetH1(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66
G4Gamma.hh