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

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

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