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