Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
medical/GammaTherapy/src/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 68699 2013-04-05 08:42:28Z gcosmo $
27 //
28 /// \file medical/GammaTherapy/src/Histo.cc
29 /// \brief Implementation of the Histo class
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  G4cout << "Histo::Save: FATAL ERROR writing ROOT file" << G4endl;
243  exit(1);
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  exit(0);
263  }
264  G4cout << "### Histo::Save: Opended file <" << nam << "> for "
265  << fNHisto << " histograms " << G4endl;
266 
267  // Creating an 1-dimensional histograms in the root directory of the tree
268  fHisto[0] = fAnalysisManager->CreateH1("10",
269  "Energy deposit at radius (mm) normalised on 1st channel",
270  fNBinsR, 0., fAbsorberR/mm);
271 
272  fHisto[1] = fAnalysisManager->CreateH1("11",
273  "Energy deposit at radius (mm) normalised to integral",
274  fNBinsR, 0., fAbsorberR/mm);
275 
276  fHisto[2] = fAnalysisManager->CreateH1("12",
277  "Energy deposit (MeV/kg/electron) at radius (mm)",
278  fNBinsR, 0., fAbsorberR/mm);
279 
280  fHisto[3] = fAnalysisManager->CreateH1("13",
281  "Energy profile (MeV/kg/electron) over Z (mm)",fNBinsZ,0.,fAbsorberZ/mm);
282 
283  fHisto[4] = fAnalysisManager->CreateH1("14",
284  "Energy profile (MeV/kg/electron) over Z (mm) at Central Voxel",
285  fNBinsZ, 0., fAbsorberZ/mm);
286 
287  fHisto[5] = fAnalysisManager->CreateH1("15",
288  "Energy (MeV) of fGamma produced in the target",
289  fNBinsE, 0., fMaxEnergy/MeV);
290 
291  fHisto[6] = fAnalysisManager->CreateH1("16",
292  "Energy (MeV) of fGamma before phantom",fNBinsE,0.,fMaxEnergy/MeV);
293 
294  fHisto[7] = fAnalysisManager->CreateH1("17",
295  "Energy (MeV) of electrons produced in phantom",fNBinsE,0.,fMaxEnergy/MeV);
296 
297  fHisto[8] = fAnalysisManager->CreateH1("18",
298  "Energy (MeV) of electrons produced in target",fNBinsE,0.,fMaxEnergy/MeV);
299 
300  fHisto[9] = fAnalysisManager->CreateH1("19",
301  "Gamma Energy Fluence (MeV/cm2) at radius(mm) in front of phantom",
302  fNBinsR, 0., fAbsorberR/mm);
303 
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307 
308 void Histo::AddTargetPhoton(const G4DynamicParticle* p)
309 {
310  ++fNgamTar;
311  if(fAnalysisManager) {
312  fAnalysisManager->FillH1(fHisto[5],p->GetKineticEnergy()/MeV,1.0);
313  }
314 }
315 
316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
317 
318 void Histo::AddPhantomPhoton(const G4DynamicParticle*)
319 {
320  ++fNgamPh;
321 }
322 
323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324 
325 void Histo::AddTargetElectron(const G4DynamicParticle* p)
326 {
327  ++fNeTar;
328  if(fAnalysisManager) {
329  fAnalysisManager->FillH1(fHisto[8],p->GetKineticEnergy()/MeV,1.0);
330  }
331 }
332 
333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
334 
335 void Histo::AddPhantomElectron(const G4DynamicParticle* p)
336 {
337  ++fNePh;
338  if(fAnalysisManager) {
339  fAnalysisManager->FillH1(fHisto[7],p->GetKineticEnergy()/MeV,1.0);
340  }
341 }
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
344 
345 void Histo::ScoreNewTrack(const G4Track* aTrack)
346 {
347  //Save primary parameters
348  const G4ParticleDefinition* particle = aTrack->GetParticleDefinition();
349  G4int pid = aTrack->GetParentID();
350  G4VPhysicalVolume* pv = aTrack->GetVolume();
351  const G4DynamicParticle* dp = aTrack->GetDynamicParticle();
352 
353  //primary particle
354  if(0 == pid) {
355 
356  ++fNevt;
357  if(1 < fVerbose) {
358  G4ThreeVector pos = aTrack->GetVertexPosition();
359  G4ThreeVector dir = aTrack->GetMomentumDirection();
360  G4cout << "TrackingAction: Primary "
361  << particle->GetParticleName()
362  << " Ekin(MeV)= "
363  << aTrack->GetKineticEnergy()/MeV
364  << "; pos= " << pos << "; dir= " << dir << G4endl;
365  }
366 
367  // secondary electron
368  } else if (0 < pid && particle == fElectron) {
369  if(1 < fVerbose) {
370  G4cout << "TrackingAction: Secondary electron " << G4endl;
371  }
372  AddElectron(dp);
373  if(pv == fPhantom) { AddPhantomElectron(dp); }
374  else if(pv == fTarget1 || pv == fTarget2) { AddTargetElectron(dp); }
375 
376  // secondary positron
377  } else if (0 < pid && particle == fPositron) {
378  if(1 < fVerbose) {
379  G4cout << "TrackingAction: Secondary positron " << G4endl;
380  }
381  AddPositron(dp);
382 
383  // secondary gamma
384  } else if (0 < pid && particle == fGamma) {
385  if(1 < fVerbose) {
386  G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
387  << " E= " << aTrack->GetKineticEnergy() << G4endl;
388  }
389  AddPhoton(dp);
390  if(pv == fPhantom) { AddPhantomPhoton(dp); }
391  else if(pv == fTarget1 || pv == fTarget2) { AddTargetPhoton(dp); }
392 
393  }
394 }
395 
396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
397 
398 void Histo::AddPhantomGamma(G4double e, G4double r)
399 {
400  e /= MeV;
401  fSumR += e;
402  G4int bin = (G4int)(e/fStepE);
403  if(bin >= fNBinsE) { bin = fNBinsE-1; }
404  fGammaE[bin] += e;
405  G4int bin1 = (G4int)(r/fStepR);
406  if(bin1 >= fNBinsR) { bin1 = fNBinsR-1; }
407  if(fAnalysisManager) {
408  fAnalysisManager->FillH1(fHisto[6],e,1.0);
409  fAnalysisManager->FillH1(fHisto[9],r,e*fVolumeR[bin1]);
410  }
411 }
412 
413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
414 
415 void Histo::AddPhantomStep(G4double edep, G4double r1, G4double z1,
416  G4double r2, G4double z2,
417  G4double r0, G4double z0)
418 {
419  ++fNstep;
420  G4int nzbin = (G4int)(z0/fStepZ);
421  if(fVerbose > 1) {
422  G4cout << "Histo: edep(MeV)= " << edep/MeV << " at binz= " << nzbin
423  << " r1= " << r1 << " z1= " << z1
424  << " r2= " << r2 << " z2= " << z2
425  << " r0= " << r0 << " z0= " << z0
426  << G4endl;
427  }
428  if(nzbin == fScoreBin) {
429  G4int bin = (G4int)(r0/fStepR);
430  if(bin >= fNBinsR) { bin = fNBinsR-1; }
431  double w = edep*fVolumeR[bin];
432  fEdep[bin] += w;
433  if(fAnalysisManager) {
434  fAnalysisManager->FillH1(fHisto[0],r0,w);
435  fAnalysisManager->FillH1(fHisto[1],r0,w);
436  fAnalysisManager->FillH1(fHisto[2],r0,w);
437  }
438  }
439  G4int bin1 = (G4int)(z1/fStepZ);
440  if(bin1 >= fNBinsZ) { bin1 = fNBinsZ-1; }
441  G4int bin2 = (G4int)(z2/fStepZ);
442  if(bin2 >= fNBinsZ) { bin2 = fNBinsZ-1; }
443  if(bin1 == bin2) {
444  if(fAnalysisManager) {
445  fAnalysisManager->FillH1(fHisto[3],z0,edep);
446  if(r1 < fStepR) {
447  G4double w = edep;
448  if(r2 > fStepR) { w *= (fStepR - r1)/(r2 - r1); }
449  fAnalysisManager->FillH1(fHisto[4],z0,w);
450  }
451  }
452  } else {
453  G4int bin;
454 
455  if(bin2 < bin1) {
456  bin = bin2;
457  G4double z = z2;
458  bin2 = bin1;
459  z2 = z1;
460  bin1 = bin;
461  z1 = z;
462  }
463  G4double zz1 = z1;
464  G4double zz2 = (bin1+1)*fStepZ;
465  G4double rr1 = r1;
466  G4double dz = z2 - z1;
467  G4double dr = r2 - r1;
468  G4double rr2 = r1 + dr*(zz2-zz1)/dz;
469  for(bin=bin1; bin<=bin2; bin++) {
470  if(fAnalysisManager) {
471  G4double de = edep*(zz2 - zz1)/dz;
472  G4double zf = (zz1+zz2)*0.5;
473  { fAnalysisManager->FillH1(fHisto[3],zf,de); }
474  if(rr1 < fStepR) {
475  G4double w = de;
476  if(rr2 > fStepR) w *= (fStepR - rr1)/(rr2 - rr1);
477  { fAnalysisManager->FillH1(fHisto[4],zf,w); }
478  }
479  }
480  zz1 = zz2;
481  zz2 = std::min(z2, zz1+fStepZ);
482  rr1 = rr2;
483  rr2 = rr1 + dr*(zz2 - zz1)/dz;
484  }
485  }
486 }
487 
488 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
489 
Histo::SetNumberDivR
void SetNumberDivR(G4int val)
Definition: medical/GammaTherapy/include/Histo.hh:99
G4Track::GetParentID
G4int GetParentID() const
Histo::AddPhantomGamma
void AddPhantomGamma(G4double e, G4double r)
Definition: medical/GammaTherapy/src/Histo.cc:398
Histo::GetPointer
static Histo * GetPointer()
Definition: medical/GammaTherapy/src/Histo.cc:59
plottest35.bin
tuple bin
Definition: plottest35.py:22
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:215
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
z
G4double z
Definition: TRTMaterials.hh:39
p
const char * p
Definition: xmltok.h:285
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
Histo::ScoreNewTrack
void ScoreNewTrack(const G4Track *aTrack)
Definition: medical/GammaTherapy/src/Histo.cc:345
G4Track::GetKineticEnergy
G4double GetKineticEnergy() const
G4cout
G4GLOB_DLL std::ostream G4cout
G4Positron.hh
G4Track::GetParticleDefinition
const G4ParticleDefinition * GetParticleDefinition() const
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4RootAnalysisManager.hh
G4Track::GetVertexPosition
const G4ThreeVector & GetVertexPosition() const
G4PhysicalConstants.hh
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)
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4RootAnalysisManager::Instance
static G4RootAnalysisManager * Instance()
Definition: G4RootAnalysisManager.cc:58
G4Track::GetVolume
G4VPhysicalVolume * GetVolume() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
Histo::AddPhantomStep
void AddPhantomStep(G4double e, G4double r1, G4double z1, G4double r2, G4double z2, G4double r0, G4double z0)
Definition: medical/GammaTherapy/src/Histo.cc:415
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4RootAnalysisManager::GetH1
tools::histo::h1d * GetH1(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
G4Gamma.hh
Generated on Wed Apr 30 2014 15:55:15 for Geant4.10 by   doxygen 1.8.7