source: trunk/examples/extended/medical/GammaTherapy/src/Histo.cc @ 1342

Last change on this file since 1342 was 1342, checked in by garnier, 14 years ago

update ti head
File size: 16.6 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// $Id: Histo.cc,v 1.10 2010/10/26 12:05:14 vnivanch Exp $
27// GEANT4 tag $Name: examples-V09-03-09 $
28//
29//---------------------------------------------------------------------------
30//
31// ClassName:   Histo
32//
33//
34// Author:      V.Ivanchenko 30/01/01
35//
36//----------------------------------------------------------------------------
37//
38
39//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41
42#include "Histo.hh"
43#include "G4Gamma.hh"
44#include "G4Electron.hh"
45#include "G4Positron.hh"
46#include "G4Neutron.hh"
47#include <iomanip>
48
49#ifdef G4ANALYSIS_USE
50#include "AIDA/AIDA.h"
51#endif
52
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55Histo* Histo::fManager = 0;
56
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58
59Histo* 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
70Histo::Histo()
71{
72  verbose   = 1;
73  histName  = G4String("histo");
74  histType  = G4String("root");
75  nHisto    = 10;
76  nHisto1   = 10;
77  maxEnergy = 50.0*MeV;
78  nTuple = false;
79  nBinsZ = 60;
80  nBinsR = 80;
81  nBinsE = 200;
82  absorberZ = 300.*mm;
83  absorberR = 200.*mm;
84  scoreZ    = 100.*mm;
85
86  gamma    = G4Gamma::Gamma();
87  electron = G4Electron::Electron();
88  positron = G4Positron::Positron();
89  neutron  = G4Neutron::Neutron();
90
91#ifdef G4ANALYSIS_USE
92  af   = 0;
93  tree = 0;
94  ntup = 0;
95#endif
96}
97
98//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99
100Histo::~Histo()
101{
102#ifdef G4ANALYSIS_USE
103  delete af;
104#endif
105}
106
107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108
109void Histo::BeginOfHisto()
110{
111  G4cout << "### Histo start initialisation nHisto= " << nHisto << G4endl;
112
113  n_evt  = 0;
114  n_elec = 0;
115  n_posit= 0;
116  n_gam  = 0;
117  n_step = 0;
118  n_gam_ph= 0;
119  n_gam_tar= 0;
120  n_e_tar = 0;
121  n_e_ph  = 0;
122  n_step_target = 0;
123  n_neutron = 0;
124  sumR = 0.0;
125  if(nBinsR>1000) SetNumberDivR(40);
126
127  stepZ = absorberZ/(G4double)nBinsZ;
128  stepR = absorberR/(G4double)nBinsR;
129  stepE = maxEnergy/(G4double)nBinsE;
130  nScoreBin = (G4int)(scoreZ/stepZ + 0.5);
131
132  G4cout << "   "<< nBinsR << " bins R   stepR= " << stepR/mm << " mm " << G4endl;
133  G4cout << "   "<< nBinsZ << " bins Z   stepZ= " << stepZ/mm << " mm " << G4endl;
134  G4cout << "   "<< nBinsE << " bins E   stepE= " << stepE/MeV << " MeV " << G4endl;
135  G4cout << "   "<< nScoreBin << "th bin in Z is used for R distribution" << G4endl;
136
137  G4int i;
138  G4double r1 = 0.0;
139  G4double r2 = stepR;
140  volumeR.clear();
141  for(i=0; i<nBinsR; i++) {
142    volumeR.push_back(cm*cm/(pi*(r2*r2 - r1*r1)));
143    r1 = r2;
144    r2 += stepR;
145  }
146  for(i=0; i<nBinsE; i++) {
147    gammaE.push_back(0.0);
148  }
149
150  bookHisto();
151
152  if(verbose > 0) {
153    G4cout << "Histo: Histograms are booked and run has been started"
154           << G4endl;
155  }
156}
157
158//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
159
160void Histo::EndOfHisto()
161{
162
163  G4cout << "Histo: End of run actions are started" << G4endl;
164
165  // average
166
167  G4cout<<"========================================================"<<G4endl;
168  G4double x = (G4double)n_evt;
169  if(n_evt > 0) x = 1.0/x;
170  G4double xe = x*(G4double)n_elec;
171  G4double xg = x*(G4double)n_gam;
172  G4double xp = x*(G4double)n_posit;
173  G4double xs = x*(G4double)n_step;
174  G4double xph= x*(G4double)n_gam_ph;
175  G4double xes= x*(G4double)n_step_target;
176  G4double xgt= x*(G4double)n_gam_tar;
177  G4double xet= x*(G4double)n_e_tar;
178  G4double xphe= x*(G4double)n_e_ph;
179  G4double xne= x*(G4double)n_neutron;
180  G4cout                    << "Number of events                                  " << n_evt <<G4endl;
181  G4cout << std::setprecision(4) << "Average number of e-                         " << xe << G4endl;
182  G4cout << std::setprecision(4) << "Average number of gamma                      " << xg << G4endl;
183  G4cout << std::setprecision(4) << "Average number of e+                         " << xp << G4endl;
184  G4cout << std::setprecision(4) << "Average number of neutrons                   " << xne << G4endl;
185  G4cout << std::setprecision(4) << "Average number of steps in absorber          " << xs << G4endl;
186  G4cout << std::setprecision(4) << "Average number of e- steps in target         " << xes << G4endl;
187  G4cout << std::setprecision(4) << "Average number of g  produced in the target  " << xgt << G4endl;
188  G4cout << std::setprecision(4) << "Average number of e- produced in the target  " << xet << G4endl;
189  G4cout << std::setprecision(4) << "Average number of g produced in the phantom  " << xph << G4endl;
190  G4cout << std::setprecision(4) << "Average number of e- produced in the phantom " << xphe << G4endl;
191  G4cout << std::setprecision(4) << "Total gamma fluence in front of phantom      " << x*sumR/MeV
192         << " MeV " << G4endl;
193  G4cout<<"========================================================"<<G4endl;
194  G4cout<<G4endl;
195  G4cout<<G4endl;
196
197#ifdef G4ANALYSIS_USE
198  if(tree) {
199    // normalise histograms
200    for(G4int i=0; i<nHisto1; i++) {
201      histo[i]->scale(x);
202    }
203    G4double nr = histo[0]->binHeight(0);
204    if(nr > 0.0) nr = 1./nr;
205    histo[0]->scale(nr);
206
207    nr = (histo[1]->sumAllBinHeights())*stepR;
208    if(nr > 0.0) nr = 1./nr;
209    histo[1]->scale(nr);
210
211    histo[3]->scale(1000.0*cm3/(pi*absorberR*absorberR*stepZ));
212    histo[4]->scale(1000.0*cm3*volumeR[0]/stepZ);
213
214    // Write histogram file
215    if(0 < nHisto) {
216      tree->commit();
217      G4cout << "Histograms and Ntuples are saved" << G4endl;
218    }
219    tree->close();
220    delete tree;
221    tree = 0;
222    G4cout << "Tree is closed" << G4endl;
223  }
224#endif
225}
226
227//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228
229void Histo::SaveEvent()
230{
231#ifdef G4ANALYSIS_USE
232  if(ntup) {
233    ntup->addRow();
234  }
235#endif
236}
237
238//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
239
240void Histo::SaveToTuple(const G4String& parname, G4double val)
241{
242  if(2 < verbose) G4cout << "Save to tuple " << parname << "   " << val << G4endl;
243#ifdef G4ANALYSIS_USE
244  if(ntup) ntup->fill( ntup->findColumn(parname), (float)val);
245#endif
246}
247
248//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249
250void Histo::SaveToTuple(const G4String& parname,G4double val, G4double)
251{
252  if(2 < verbose) G4cout << "Save to tuple " << parname << "   " << val << G4endl;
253#ifdef G4ANALYSIS_USE
254  if(ntup) ntup->fill( ntup->findColumn(parname), (float)val);
255#endif
256}
257
258//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
259
260void Histo::bookHisto()
261{
262
263#ifdef G4ANALYSIS_USE
264  G4cout << "### Histo books " << nHisto << " histograms " << G4endl;
265  // Creating the analysis factory
266  if(!af) af = AIDA_createAnalysisFactory();
267
268  // Creating the tree factory
269  AIDA::ITreeFactory* tf = af->createTreeFactory();
270
271  // Creating a tree mapped to a new hbook file.
272  G4String tt = histType;
273  G4String nn = histName + "." + histType;
274  if(histType == "xml" || histType == "XML" || histType == "aida" || 
275     histType == "AIDA") {
276    tt = "xml";
277    nn = histName + ".aida";
278  }
279
280  tree = tf->create(nn,tt,false,true, "");
281  if(tree) {
282    G4cout << "Tree store : " << tree->storeName() << G4endl;
283  } else {
284    G4cout << "Fail to open tree store " << nn << G4endl;
285    return;
286  }
287  delete tf;
288  histo.resize(nHisto1);
289
290  // Creating a histogram factory, whose histograms will be handled by the tree
291  AIDA::IHistogramFactory* hf = af->createHistogramFactory( *tree );
292
293  // Creating an 1-dimensional histograms in the root directory of the tree
294
295  histo[0] = hf->createHistogram1D("10",
296    "Energy deposit at radius (mm) normalised on 1st channel",nBinsR,0.,absorberR/mm);
297
298  histo[1] = hf->createHistogram1D("11",
299    "Energy deposit at radius (mm) normalised to integral",nBinsR,0.,absorberR/mm);
300
301  histo[2] = hf->createHistogram1D("12",
302    "Energy deposit (MeV/kg/electron) at radius (mm)",nBinsR,0.,absorberR/mm);
303
304  histo[3] = hf->createHistogram1D("13",
305    "Energy profile (MeV/kg/electron) over Z (mm)",nBinsZ,0.,absorberZ/mm);
306
307  histo[4] = hf->createHistogram1D("14",
308    "Energy profile (MeV/kg/electron) over Z (mm) at Central Voxel",nBinsZ,0.,absorberZ/mm);
309
310  histo[5] = hf->createHistogram1D("15",
311    "Energy (MeV) of gamma produced in the target",nBinsE,0.,maxEnergy/MeV);
312
313  histo[6] = hf->createHistogram1D("16",
314    "Energy (MeV) of gamma before phantom",nBinsE,0.,maxEnergy/MeV);
315
316  histo[7] = hf->createHistogram1D("17",
317    "Energy (MeV) of electrons produced in phantom",nBinsE,0.,maxEnergy/MeV);
318
319  histo[8] = hf->createHistogram1D("18",
320    "Energy (MeV) of electrons produced in target",nBinsE,0.,maxEnergy/MeV);
321
322  histo[9] = hf->createHistogram1D("19",
323    "Gamma Energy Fluence (MeV/cm2) at radius(mm) in front of phantom",nBinsR,0.,absorberR/mm);
324
325  // Creating a tuple factory, whose tuples will be handled by the tree
326  AIDA::ITupleFactory* tpf = af->createTupleFactory( *tree );
327
328  // If using Anaphe HBOOK implementation, there is a limitation on the
329  // length of the variable names in a ntuple
330  if(nTuple) {
331     ntup = tpf->create( "100", "Dose deposite","float r, z, e" );
332  }
333  delete hf;
334  delete tpf;
335#endif
336}
337
338//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339
340void Histo::AddDeltaElectron(const G4DynamicParticle* elec)
341{
342  G4double e = elec->GetKineticEnergy()/MeV;
343  if(e > 0.0) { ++n_elec; }
344}
345
346//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
347
348void Histo::AddPhoton(const G4DynamicParticle* ph)
349{
350  G4double e = ph->GetKineticEnergy()/MeV;
351  if(e > 0.0) { ++n_gam; }
352}
353
354//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355
356void Histo::AddTargetPhoton(const G4DynamicParticle* ph)
357{
358  G4double e = ph->GetKineticEnergy()/MeV;
359  if(e > 0.0) { ++n_gam_tar; }
360#ifdef G4ANALYSIS_USE
361  if(tree) histo[5]->fill(e,1.0);
362#endif
363}
364
365//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
366
367void Histo::AddPhantomPhoton(const G4DynamicParticle* ph)
368{
369  G4double e = ph->GetKineticEnergy()/MeV;
370  if(e > 0.0) { ++n_gam_ph; }
371}
372
373//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
374
375void Histo::AddTargetElectron(const G4DynamicParticle* ph)
376{
377  G4double e = ph->GetKineticEnergy()/MeV;
378  if(e > 0.0) { ++n_e_tar; }
379#ifdef G4ANALYSIS_USE
380  if(tree) histo[8]->fill(e,1.0);
381#endif
382}
383
384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
385
386void Histo::AddPhantomElectron(const G4DynamicParticle* ph)
387{
388  G4double e = ph->GetKineticEnergy()/MeV;
389  if(e > 0.0) { ++n_e_ph; }
390#ifdef G4ANALYSIS_USE
391  if(tree) histo[7]->fill(e,1.0);
392#endif
393}
394
395//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
396
397void Histo::ScoreNewTrack(const G4Track* aTrack)
398{
399  //Save primary parameters
400
401  ResetTrackLength();
402  const G4ParticleDefinition* particle = aTrack->GetParticleDefinition();
403  G4int pid = aTrack->GetParentID();
404  G4double kinE = aTrack->GetKineticEnergy();
405  G4ThreeVector pos = aTrack->GetVertexPosition();
406  G4VPhysicalVolume* pv = aTrack->GetVolume();
407  const G4DynamicParticle* dp = aTrack->GetDynamicParticle();
408
409  if(0 == pid) {
410
411    SaveToTuple("TKIN", kinE/MeV);
412
413    G4double mass = 0.0;
414    if(particle) {
415      mass = particle->GetPDGMass();
416      SaveToTuple("MASS", mass/MeV);
417      SaveToTuple("CHAR",(particle->GetPDGCharge())/eplus);
418    }
419
420    G4ThreeVector dir = aTrack->GetMomentumDirection();
421    if(1 < verbose) {
422      G4cout << "TrackingAction: Primary kinE(MeV)= " << kinE/MeV
423           << "; m(MeV)= " << mass/MeV
424           << "; pos= " << pos << ";  dir= " << dir << G4endl;
425    }
426
427    // delta-electron
428  } else if (0 < pid && particle == electron) {
429    if(1 < verbose) {
430      G4cout << "TrackingAction: Secondary electron " << G4endl;
431    }
432    AddDeltaElectron(dp);
433    if(pv == phantom)                       { AddPhantomElectron(dp); }
434    else if(pv == target1 || pv == target2) { AddTargetElectron(dp); }
435
436  } else if (0 < pid && particle == positron) {
437    if(1 < verbose) {
438      G4cout << "TrackingAction: Secondary positron " << G4endl;
439    }
440    AddPositron(dp);
441
442  } else if (0 < pid && particle == gamma) {
443    if(1 < verbose) {
444      G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
445             << " E= " << aTrack->GetKineticEnergy() << G4endl;
446    }
447    AddPhoton(dp);
448    if(pv == phantom)                       { AddPhantomPhoton(dp); }
449    else if(pv == target1 || pv == target2) { AddTargetPhoton(dp); }
450
451  } else if (0 < pid && particle == neutron) {
452    n_neutron++;
453    if(1 < verbose) {
454      G4cout << "TrackingAction: Secondary neutron; parentID= " << pid
455             << " E= " << aTrack->GetKineticEnergy() << G4endl;
456    }
457  }
458}
459
460//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
461
462void Histo::AddGamma(G4double e, G4double r)
463{
464  e /= MeV;
465  sumR += e;
466  G4int bin = (G4int)(e/stepE);
467  if(bin >= nBinsE) { bin = nBinsE-1; }
468  gammaE[bin] += e;
469  G4int bin1 = (G4int)(r/stepR);
470  if(bin1 >= nBinsR) { bin1 = nBinsR-1; }
471#ifdef G4ANALYSIS_USE
472  if(tree) {
473    histo[6]->fill(e,1.0);
474    histo[9]->fill(r,e*volumeR[bin1]);
475  }
476#endif
477
478}
479
480//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
481
482void Histo::AddStep(G4double edep, G4double r1, G4double z1, G4double r2, G4double z2,
483                                   G4double r0, G4double z0)
484{
485  n_step++;
486  G4int nzbin = (G4int)(z0/stepZ);
487  if(verbose > 1) {
488    G4cout << "Histo: edep(MeV)= " << edep/MeV << " at binz= " << nzbin
489           << " r1= " << r1 << " z1= " << z1
490           << " r2= " << r2 << " z2= " << z2
491           << " r0= " << r0 << " z0= " << z0
492           << G4endl;
493  }
494  if(nzbin == nScoreBin) {
495#ifdef G4ANALYSIS_USE
496    if(tree) {
497      G4int bin  = (G4int)(r0/stepR);
498      if(bin >= nBinsR) bin = nBinsR-1;
499      double w    = edep*volumeR[bin];
500      histo[0]->fill(r0,w);
501      histo[1]->fill(r0,w);
502      histo[2]->fill(r0,w);
503    }
504#endif
505  }
506  G4int bin1 = (G4int)(z1/stepZ);
507  if(bin1 >= nBinsZ) bin1 = nBinsZ-1;
508  G4int bin2 = (G4int)(z2/stepZ);
509  if(bin2 >= nBinsZ) bin2 = nBinsZ-1;
510  if(bin1 == bin2) {
511#ifdef G4ANALYSIS_USE
512    if(tree) {
513      histo[3]->fill(z0,edep);
514      if(r1 < stepR) {
515        G4double w = edep;
516        if(r2 > stepR) w *= (stepR - r1)/(r2 - r1);
517        histo[4]->fill(z0,w);
518      }
519    }
520#endif
521  } else {
522    G4int bin;
523
524    if(bin2 < bin1) {
525      bin = bin2;
526      G4double z = z2;
527      bin2 = bin1;
528      z2 = z1;
529      bin1 = bin;
530      z1 = z;
531    }
532    G4double zz1 = z1;
533    G4double zz2 = (bin1+1)*stepZ;
534    G4double rr1 = r1;
535    G4double dz  = z2 - z1;
536    G4double dr  = r2 - r1;
537    G4double rr2 = r1 + dr*(zz2-zz1)/dz;
538    for(bin=bin1; bin<=bin2; bin++) {
539#ifdef G4ANALYSIS_USE
540      if(tree) {
541        G4double de = edep*(zz2 - zz1)/dz;
542        G4double zf = (zz1+zz2)*0.5;
543        histo[3]->fill(zf,de);
544        if(rr1 < stepR) {
545          G4double w = de;
546          if(rr2 > stepR) w *= (stepR - rr1)/(rr2 - rr1);
547          histo[4]->fill(zf,w);
548        }
549      }
550#endif
551      zz1 = zz2;
552      zz2 = std::min(z2, zz1+stepZ);
553      rr1 = rr2;
554      rr2 = rr1 + dr*(zz2 - zz1)/dz;
555    }
556  }
557}
558
559//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
560
Note: See TracBrowser for help on using the repository browser.