Geant4  10.00.p02
HadrontherapyAnalysisManager.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 // This is the *BASIC* version of Hadrontherapy, a Geant4-based application
27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
28 //
29 // Visit the Hadrontherapy web site (http://www.lns.infn.it/link/Hadrontherapy) to request
30 // the *COMPLETE* version of this program, together with its documentation;
31 // Hadrontherapy (both basic and full version) are supported by the Italian INFN
32 // Institute in the framework of the MC-INFN Group
33 //
34 
36 #include "HadrontherapyMatrix.hh"
38 #include "G4SystemOfUnits.hh"
39 #include <time.h>
40 
42 
44 #ifdef G4ANALYSIS_USE_ROOT
45 :
46 analysisFileName("DoseDistribution.root"),theTFile(0), histo1(0), histo2(0), histo3(0),
47 histo4(0), histo5(0), histo6(0), histo7(0), histo8(0), histo9(0), histo10(0), histo11(0), histo12(0), histo13(0), histo14(0), histo15(0), histo16(0),
48 kinFragNtuple(0),
49 kineticEnergyPrimaryNtuple(0),
50 doseFragNtuple(0),
51 fluenceFragNtuple(0),
52 letFragNtuple(0),
53 theROOTNtuple(0),
54 theROOTIonTuple(0),
55 fragmentNtuple(0),
56 metaData(0),
57 eventCounter(0)
58 #endif
59 {
60  fMess = new HadrontherapyAnalysisFileMessenger(this);
61 }
63 
65 {
66  delete fMess;
67 #ifdef G4ANALYSIS_USE_ROOT
68  Clear();
69 #endif
70 }
71 
73 {
75  return instance;
76 }
77 #ifdef G4ANALYSIS_USE_ROOT
79 {
80  if (theTFile)
81  {
82  delete metaData;
83  metaData = 0;
84 
85  delete fragmentNtuple;
86  fragmentNtuple = 0;
87 
88  delete theROOTIonTuple;
89  theROOTIonTuple = 0;
90 
91  delete theROOTNtuple;
92  theROOTNtuple = 0;
93 
94 
95  delete histo16;
96  histo16 = 0;
97 
98  delete histo15;
99  histo15 = 0;
100 
101  delete histo14;
102  histo14 = 0;
103 
104  delete histo13;
105  histo13 = 0;
106 
107  delete histo12;
108  histo12 = 0;
109 
110  delete histo11;
111  histo11 = 0;
112 
113  delete histo10;
114  histo10 = 0;
115 
116  delete histo9;
117  histo9 = 0;
118 
119  delete histo8;
120  histo8 = 0;
121 
122  delete histo7;
123  histo7 = 0;
124 
125  delete histo6;
126  histo6 = 0;
127 
128  delete histo5;
129  histo5 = 0;
130 
131  delete histo4;
132  histo4 = 0;
133 
134  delete histo3;
135  histo3 = 0;
136 
137  delete histo2;
138  histo2 = 0;
139 
140  delete histo1;
141  histo1 = 0;
142  }
143 }
145 
146 void HadrontherapyAnalysisManager::SetAnalysisFileName(G4String aFileName)
147 {
148  this->analysisFileName = aFileName;
149 }
150 
152 G4bool HadrontherapyAnalysisManager::IsTheTFile()
153 {
154  return (theTFile) ? true:false;
155 }
156 void HadrontherapyAnalysisManager::book()
157 {
158  delete theTFile; // this is similar to theTFile->Close() => delete all associated variables created via new, moreover it delete itself.
159 
160  theTFile = new TFile(analysisFileName, "RECREATE");
161 
162  // Create the histograms with the energy deposit along the X axis
163  histo1 = createHistogram1D("braggPeak","slice, energy", 400, 0., 80); //<different waterthicknesses are accoutned for in ROOT-analysis stage
164  histo2 = createHistogram1D("h20","Secondary protons - slice, energy", 400, 0., 400.);
165  histo3 = createHistogram1D("h30","Secondary neutrons - slice, energy", 400, 0., 400.);
166  histo4 = createHistogram1D("h40","Secondary alpha - slice, energy", 400, 0., 400.);
167  histo5 = createHistogram1D("h50","Secondary gamma - slice, energy", 400, 0., 400.);
168  histo6 = createHistogram1D("h60","Secondary electron - slice, energy", 400, 0., 400.);
169  histo7 = createHistogram1D("h70","Secondary triton - slice, energy", 400, 0., 400.);
170  histo8 = createHistogram1D("h80","Secondary deuteron - slice, energy", 400, 0., 400.);
171  histo9 = createHistogram1D("h90","Secondary pion - slice, energy", 400, 0., 400.);
172  histo10 = createHistogram1D("h100","Energy distribution of secondary electrons", 70, 0., 70.);
173  histo11 = createHistogram1D("h110","Energy distribution of secondary photons", 70, 0., 70.);
174  histo12 = createHistogram1D("h120","Energy distribution of secondary deuterons", 70, 0., 70.);
175  histo13 = createHistogram1D("h130","Energy distribution of secondary tritons", 70, 0., 70.);
176  histo14 = createHistogram1D("h140","Energy distribution of secondary alpha particles", 70, 0., 70.);
177  histo15 = createHistogram1D("heliumEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
178  70, 0., 500.);
179  histo16 = createHistogram1D("hydrogenEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
180  70, 0., 500.);
181 
182  kinFragNtuple = new TNtuple("kinFragNtuple",
183  "Kinetic energy by voxel & fragment",
184  "i:j:k:A:Z:kineticEnergy");
185  kineticEnergyPrimaryNtuple= new TNtuple("kineticEnergyPrimaryNtuple",
186  "Kinetic energy by voxel of primary",
187  "i:j:k:kineticEnergy");
188  doseFragNtuple = new TNtuple("doseFragNtuple",
189  "Energy deposit by voxel & fragment",
190  "i:j:k:A:Z:energy");
191 
192  fluenceFragNtuple = new TNtuple("fluenceFragNtuple",
193  "Fluence by voxel & fragment",
194  "i:j:k:A:Z:fluence");
195 
196  letFragNtuple = new TNtuple("letFragNtuple",
197  "Let by voxel & fragment",
198  "i:j:k:A:Z:letT:letD");
199 
200  theROOTNtuple = new TNtuple("theROOTNtuple",
201  "Energy deposit by slice",
202  "i:j:k:energy");
203 
204  theROOTIonTuple = new TNtuple("theROOTIonTuple",
205  "Generic ion information",
206  "a:z:occupancy:energy");
207 
208  fragmentNtuple = new TNtuple("fragmentNtuple",
209  "Fragments",
210  "A:Z:energy:posX:posY:posZ");
211 
212  metaData = new TNtuple("metaData",
213  "Metadata",
214  "events:detectorDistance:waterThickness:beamEnergy:energyError:phantomCenterDistance");
215 }
216 
218 void HadrontherapyAnalysisManager::FillEnergyDeposit(G4int i,
219  G4int j,
220  G4int k,
222 {
223  if (theROOTNtuple)
224  {
225  theROOTNtuple->Fill(i, j, k, energy);
226  }
227 }
228 
230 void HadrontherapyAnalysisManager::BraggPeak(G4int slice, G4double energy)
231 {
232  histo1->SetBinContent(slice, energy); //This uses setbincontent instead of fill to get labels correct
233 }
234 
236 void HadrontherapyAnalysisManager::SecondaryProtonEnergyDeposit(G4int slice, G4double energy)
237 {
238  histo2->Fill(slice, energy);
239 }
240 
242 void HadrontherapyAnalysisManager::SecondaryNeutronEnergyDeposit(G4int slice, G4double energy)
243 {
244  histo3->Fill(slice, energy);
245 }
246 
248 void HadrontherapyAnalysisManager::SecondaryAlphaEnergyDeposit(G4int slice, G4double energy)
249 {
250  histo4->Fill(slice, energy);
251 }
252 
254 void HadrontherapyAnalysisManager::SecondaryGammaEnergyDeposit(G4int slice, G4double energy)
255 {
256  histo5->Fill(slice, energy);
257 }
258 
260 void HadrontherapyAnalysisManager::SecondaryElectronEnergyDeposit(G4int slice, G4double energy)
261 {
262  histo6->Fill(slice, energy);
263 }
264 
266 void HadrontherapyAnalysisManager::SecondaryTritonEnergyDeposit(G4int slice, G4double energy)
267 {
268  histo7->Fill(slice, energy);
269 }
270 
272 void HadrontherapyAnalysisManager::SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy)
273 {
274  histo8->Fill(slice, energy);
275 }
276 
278 void HadrontherapyAnalysisManager::SecondaryPionEnergyDeposit(G4int slice, G4double energy)
279 {
280  histo9->Fill(slice, energy);
281 }
282 
284 void HadrontherapyAnalysisManager::electronEnergyDistribution(G4double energy)
285 {
286  histo10->Fill(energy);
287 }
288 
290 void HadrontherapyAnalysisManager::gammaEnergyDistribution(G4double energy)
291 {
292  histo11->Fill(energy);
293 }
294 
296 void HadrontherapyAnalysisManager::deuteronEnergyDistribution(G4double energy)
297 {
298  histo12->Fill(energy);
299 }
300 
302 void HadrontherapyAnalysisManager::tritonEnergyDistribution(G4double energy)
303 {
304  histo13->Fill(energy);
305 }
306 
308 void HadrontherapyAnalysisManager::alphaEnergyDistribution(G4double energy)
309 {
310  histo14->Fill(energy);
311 }
313 void HadrontherapyAnalysisManager::heliumEnergy(G4double secondaryParticleKineticEnergy)
314 {
315  histo15->Fill(secondaryParticleKineticEnergy);
316 }
317 
319 void HadrontherapyAnalysisManager::hydrogenEnergy(G4double secondaryParticleKineticEnergy)
320 {
321  histo16->Fill(secondaryParticleKineticEnergy);
322 }
323 
325  // FillKineticFragmentTuple create an ntuple where the voxel indexs, the atomic number and mass and the kinetic
326  // energy of all the particles interacting with the phantom, are stored
327 void HadrontherapyAnalysisManager::FillKineticFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double kinEnergy)
328 {
329  kinFragNtuple -> Fill(i, j, k, A, Z, kinEnergy);
330 }
331 
333  // FillKineticEnergyPrimaryNTuple creates a ntuple where the voxel indexs and the kinetic
334  // energies of ONLY primary particles interacting with the phantom, are stored
335 void HadrontherapyAnalysisManager::FillKineticEnergyPrimaryNTuple(G4int i, G4int j, G4int k, G4double kinEnergy)
336 {
337  kineticEnergyPrimaryNtuple -> Fill(i, j, k, kinEnergy);
338 }
339 
341  // This function is called only if ROOT is activated.
342  // It is called by the HadrontherapyMatric.cc class file and it is used to create two ntuples containing
343  // the total energy deposited and the fluence values, in each voxel and per any particle (primary
344  // and secondary particles beam)
345 void HadrontherapyAnalysisManager::FillVoxelFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double energy, G4double fluence)
346 {
347  // Fill the ntuple containing the voxel, mass and atomic number and the energy deposited
348  doseFragNtuple -> Fill( i, j, k, A, Z, energy );
349 
350  // Fill the ntuple containing the voxel, mass and atomic number and the fluence
351  if (i==1 && Z==1) {
352  fluenceFragNtuple -> Fill( i, j, k, A, Z, fluence );
353 
354  }
355 }
356 
357 void HadrontherapyAnalysisManager::FillLetFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double letD)
358 {
359  letFragNtuple -> Fill( i, j, k, A, Z, letD/MeV);
360 }
362 void HadrontherapyAnalysisManager::FillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ)
363 {
364  fragmentNtuple->Fill(A, Z, energy, posX, posY, posZ);
365 }
366 
368 void HadrontherapyAnalysisManager::genericIonInformation(G4int a,
369  G4double z,
370  G4int electronOccupancy,
371  G4double energy)
372 {
373  if (theROOTIonTuple) {
374  theROOTIonTuple->Fill(a, z, electronOccupancy, energy);
375  }
376 }
377 
379 void HadrontherapyAnalysisManager::startNewEvent()
380 {
381  eventCounter++;
382 }
384 void HadrontherapyAnalysisManager::setGeometryMetaData(G4double endDetectorPosition, G4double waterThickness, G4double phantomCenter)
385 {
386  this->detectorDistance = endDetectorPosition;
387  this->phantomDepth = waterThickness;
388  this->phantomCenterDistance = phantomCenter;
389 }
390 void HadrontherapyAnalysisManager::setBeamMetaData(G4double meanKineticEnergy,G4double sigmaEnergy)
391 {
392  this->beamEnergy = meanKineticEnergy;
393  this->energyError = sigmaEnergy;
394 }
396 // Flush data & close the file
397 void HadrontherapyAnalysisManager::flush()
398 {
399  if (theTFile)
400  {
401  theTFile -> Write();
402  theTFile -> Close();
403  }
404  theTFile = 0;
405  eventCounter = 0;
406 }
407 
408 #endif
static HadrontherapyAnalysisManager * GetInstance()
Get the pointer to the analysis manager.
static const double MeV
Definition: G4SIunits.hh:193
G4double z
Definition: TRTMaterials.hh:39
void Clear(Node *)
static HadrontherapyAnalysisManager * instance
G4double a
Definition: TRTMaterials.hh:39
A messenger object of this class is created by the AnalysisManager.
int G4int
Definition: G4Types.hh:78
bool G4bool
Definition: G4Types.hh:79
HadrontherapyAnalysisFileMessenger * fMess
static const G4double A[nN]
G4double energy(const ThreeVector &p, const G4double m)
HadrontherapyAnalysisManager()
Analysis manager is a singleton object (there is only one instance).
A class for connecting the simulation to an analysis package.
double G4double
Definition: G4Types.hh:76