source: trunk/examples/advanced/hadrontherapy/src/HadrontherapyAnalysisManager.cc@ 1344

Last change on this file since 1344 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 13.9 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: HadrontherapyAnalisysManager.cc;
27// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
28
29#include "HadrontherapyAnalysisManager.hh"
30#include "HadrontherapyMatrix.hh"
31#include "HadrontherapyAnalysisFileMessenger.hh"
32#include <time.h>
33
34HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::instance = 0;
35
36HadrontherapyAnalysisManager::HadrontherapyAnalysisManager():
37#ifdef G4ANALYSIS_USE_ROOT
38analysisFileName("DoseDistribution.root"),theTFile(0), histo1(0), histo2(0), histo3(0),
39histo4(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),
40kinFragNtuple(0),
41kineticEnergyPrimaryNtuple(0),
42doseFragNtuple(0),
43fluenceFragNtuple(0),
44letFragNtuple(0),
45theROOTNtuple(0),
46theROOTIonTuple(0),
47fragmentNtuple(0),
48metaData(0),
49eventCounter(0)
50{
51 fMess = new HadrontherapyAnalysisFileMessenger(this);
52}
53#endif
54/////////////////////////////////////////////////////////////////////////////
55
56HadrontherapyAnalysisManager::~HadrontherapyAnalysisManager()
57{
58 delete(fMess);
59#ifdef G4ANALYSIS_USE_ROOT
60 Clear();
61#endif
62}
63
64HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::GetInstance()
65{
66 if (instance == 0) instance = new HadrontherapyAnalysisManager;
67 return instance;
68}
69#ifdef G4ANALYSIS_USE_ROOT
70void HadrontherapyAnalysisManager::Clear()
71{
72 // XXX delete ALL variables!!
73 // but this seems to be automatically done by
74 // theTFile->Close() command!
75 // so we get a *double free or corruption*
76
77 delete metaData;
78 metaData = 0;
79
80 delete fragmentNtuple;
81 fragmentNtuple = 0;
82
83 delete theROOTIonTuple;
84 theROOTIonTuple = 0;
85
86 delete theROOTNtuple;
87 theROOTNtuple = 0;
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/////////////////////////////////////////////////////////////////////////////
138
139void HadrontherapyAnalysisManager::SetAnalysisFileName(G4String aFileName)
140{
141 this->analysisFileName = aFileName;
142}
143
144 /////////////////////////////////////////////////////////////////////////////
145void HadrontherapyAnalysisManager::book()
146{
147 delete theTFile; // this is similar to theTFile->Close() => delete all associated variables created via new, moreover it delete itself.
148
149 theTFile = new TFile(analysisFileName, "RECREATE");
150
151 // Create the histograms with the energy deposit along the X axis
152 histo1 = createHistogram1D("braggPeak","slice, energy", 400, 0., 80); //<different waterthicknesses are accoutned for in ROOT-analysis stage
153 histo2 = createHistogram1D("h20","Secondary protons - slice, energy", 400, 0., 400.);
154 histo3 = createHistogram1D("h30","Secondary neutrons - slice, energy", 400, 0., 400.);
155 histo4 = createHistogram1D("h40","Secondary alpha - slice, energy", 400, 0., 400.);
156 histo5 = createHistogram1D("h50","Secondary gamma - slice, energy", 400, 0., 400.);
157 histo6 = createHistogram1D("h60","Secondary electron - slice, energy", 400, 0., 400.);
158 histo7 = createHistogram1D("h70","Secondary triton - slice, energy", 400, 0., 400.);
159 histo8 = createHistogram1D("h80","Secondary deuteron - slice, energy", 400, 0., 400.);
160 histo9 = createHistogram1D("h90","Secondary pion - slice, energy", 400, 0., 400.);
161 histo10 = createHistogram1D("h100","Energy distribution of secondary electrons", 70, 0., 70.);
162 histo11 = createHistogram1D("h110","Energy distribution of secondary photons", 70, 0., 70.);
163 histo12 = createHistogram1D("h120","Energy distribution of secondary deuterons", 70, 0., 70.);
164 histo13 = createHistogram1D("h130","Energy distribution of secondary tritons", 70, 0., 70.);
165 histo14 = createHistogram1D("h140","Energy distribution of secondary alpha particles", 70, 0., 70.);
166 histo15 = createHistogram1D("heliumEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
167 70, 0., 500.);
168 histo16 = createHistogram1D("hydrogenEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
169 70, 0., 500.);
170
171 kinFragNtuple = new TNtuple("kinFragNtuple",
172 "Kinetic energy by voxel & fragment",
173 "i:j:k:A:Z:kineticEnergy");
174 kineticEnergyPrimaryNtuple= new TNtuple("kineticEnergyPrimaryNtuple",
175 "Kinetic energy by voxel of primary",
176 "i:j:k:kineticEnergy");
177 doseFragNtuple = new TNtuple("doseFragNtuple",
178 "Energy deposit by voxel & fragment",
179 "i:j:k:A:Z:energy");
180
181 fluenceFragNtuple = new TNtuple("fluenceFragNtuple",
182 "Fluence by voxel & fragment",
183 "i:j:k:A:Z:fluence");
184
185 letFragNtuple = new TNtuple("letFragNtuple",
186 "Let by voxel & fragment",
187 "i:j:k:A:Z:letT:letD");
188
189 theROOTNtuple = new TNtuple("theROOTNtuple",
190 "Energy deposit by slice",
191 "i:j:k:energy");
192
193 theROOTIonTuple = new TNtuple("theROOTIonTuple",
194 "Generic ion information",
195 "a:z:occupancy:energy");
196
197 fragmentNtuple = new TNtuple("fragmentNtuple",
198 "Fragments",
199 "A:Z:energy:posX:posY:posZ");
200
201 metaData = new TNtuple("metaData",
202 "Metadata",
203 "events:detectorDistance:waterThickness:beamEnergy:energyError:phantomCenterDistance");
204}
205
206 /////////////////////////////////////////////////////////////////////////////
207void HadrontherapyAnalysisManager::FillEnergyDeposit(G4int i,
208 G4int j,
209 G4int k,
210 G4double energy)
211{
212 if (theROOTNtuple) {
213 theROOTNtuple->Fill(i, j, k, energy);
214 }
215}
216
217 /////////////////////////////////////////////////////////////////////////////
218void HadrontherapyAnalysisManager::BraggPeak(G4int slice, G4double energy)
219{
220 histo1->SetBinContent(slice, energy); //This uses setbincontent instead of fill to get labels correct
221}
222
223 /////////////////////////////////////////////////////////////////////////////
224void HadrontherapyAnalysisManager::SecondaryProtonEnergyDeposit(G4int slice, G4double energy)
225{
226 histo2->Fill(slice, energy);
227}
228
229 /////////////////////////////////////////////////////////////////////////////
230void HadrontherapyAnalysisManager::SecondaryNeutronEnergyDeposit(G4int slice, G4double energy)
231{
232 histo3->Fill(slice, energy);
233}
234
235 /////////////////////////////////////////////////////////////////////////////
236void HadrontherapyAnalysisManager::SecondaryAlphaEnergyDeposit(G4int slice, G4double energy)
237{
238 histo4->Fill(slice, energy);
239}
240
241 /////////////////////////////////////////////////////////////////////////////
242void HadrontherapyAnalysisManager::SecondaryGammaEnergyDeposit(G4int slice, G4double energy)
243{
244 histo5->Fill(slice, energy);
245}
246
247 /////////////////////////////////////////////////////////////////////////////
248void HadrontherapyAnalysisManager::SecondaryElectronEnergyDeposit(G4int slice, G4double energy)
249{
250 histo6->Fill(slice, energy);
251}
252
253 /////////////////////////////////////////////////////////////////////////////
254void HadrontherapyAnalysisManager::SecondaryTritonEnergyDeposit(G4int slice, G4double energy)
255{
256 histo7->Fill(slice, energy);
257}
258
259 /////////////////////////////////////////////////////////////////////////////
260void HadrontherapyAnalysisManager::SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy)
261{
262 histo8->Fill(slice, energy);
263}
264
265 /////////////////////////////////////////////////////////////////////////////
266void HadrontherapyAnalysisManager::SecondaryPionEnergyDeposit(G4int slice, G4double energy)
267{
268 histo9->Fill(slice, energy);
269}
270
271 /////////////////////////////////////////////////////////////////////////////
272void HadrontherapyAnalysisManager::electronEnergyDistribution(G4double energy)
273{
274 histo10->Fill(energy);
275}
276
277 /////////////////////////////////////////////////////////////////////////////
278void HadrontherapyAnalysisManager::gammaEnergyDistribution(G4double energy)
279{
280 histo11->Fill(energy);
281}
282
283 /////////////////////////////////////////////////////////////////////////////
284void HadrontherapyAnalysisManager::deuteronEnergyDistribution(G4double energy)
285{
286 histo12->Fill(energy);
287}
288
289 /////////////////////////////////////////////////////////////////////////////
290void HadrontherapyAnalysisManager::tritonEnergyDistribution(G4double energy)
291{
292 histo13->Fill(energy);
293}
294
295 /////////////////////////////////////////////////////////////////////////////
296void HadrontherapyAnalysisManager::alphaEnergyDistribution(G4double energy)
297{
298 histo14->Fill(energy);
299}
300 /////////////////////////////////////////////////////////////////////////////
301void HadrontherapyAnalysisManager::heliumEnergy(G4double secondaryParticleKineticEnergy)
302{
303 histo15->Fill(secondaryParticleKineticEnergy);
304}
305
306 /////////////////////////////////////////////////////////////////////////////
307void HadrontherapyAnalysisManager::hydrogenEnergy(G4double secondaryParticleKineticEnergy)
308{
309 histo16->Fill(secondaryParticleKineticEnergy);
310}
311
312 /////////////////////////////////////////////////////////////////////////////
313 // FillKineticFragmentTuple create an ntuple where the voxel indexs, the atomic number and mass and the kinetic
314 // energy of all the particles interacting with the phantom, are stored
315void HadrontherapyAnalysisManager::FillKineticFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double kinEnergy)
316{
317 kinFragNtuple -> Fill(i, j, k, A, Z, kinEnergy);
318}
319
320 /////////////////////////////////////////////////////////////////////////////
321 // FillKineticEnergyPrimaryNTuple creates a ntuple where the voxel indexs and the kinetic
322 // energies of ONLY primary particles interacting with the phantom, are stored
323void HadrontherapyAnalysisManager::FillKineticEnergyPrimaryNTuple(G4int i, G4int j, G4int k, G4double kinEnergy)
324{
325 kineticEnergyPrimaryNtuple -> Fill(i, j, k, kinEnergy);
326}
327
328 /////////////////////////////////////////////////////////////////////////////
329 // This function is called only if ROOT is activated.
330 // It is called by the HadrontherapyMatric.cc class file and it is used to create two ntuples containing
331 // the total energy deposited and the fluence values, in each voxel and per any particle (primary beam
332 // and secondaries )
333void HadrontherapyAnalysisManager::FillVoxelFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double energy, G4double fluence)
334{
335 // Fill the ntuple containing the voxel, mass and atomic number and the energy deposited
336 doseFragNtuple -> Fill( i, j, k, A, Z, energy );
337
338 // Fill the ntuple containing the voxel, mass and atomic number and the fluence
339 if (i==1 && Z==1) {
340 fluenceFragNtuple -> Fill( i, j, k, A, Z, fluence );
341
342 }
343}
344
345void HadrontherapyAnalysisManager::FillLetFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double letT, G4double letD)
346{
347 letFragNtuple -> Fill( i, j, k, A, Z, letT, letD);
348
349}
350 /////////////////////////////////////////////////////////////////////////////
351void HadrontherapyAnalysisManager::FillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ)
352{
353 fragmentNtuple->Fill(A, Z, energy, posX, posY, posZ);
354}
355
356 /////////////////////////////////////////////////////////////////////////////
357void HadrontherapyAnalysisManager::genericIonInformation(G4int a,
358 G4double z,
359 G4int electronOccupancy,
360 G4double energy)
361{
362 if (theROOTIonTuple) {
363 theROOTIonTuple->Fill(a, z, electronOccupancy, energy);
364 }
365}
366
367 /////////////////////////////////////////////////////////////////////////////
368void HadrontherapyAnalysisManager::startNewEvent()
369{
370 eventCounter++;
371}
372 /////////////////////////////////////////////////////////////////////////////
373void HadrontherapyAnalysisManager::setGeometryMetaData(G4double endDetectorPosition, G4double waterThickness, G4double phantomCenter)
374{
375 this->detectorDistance = endDetectorPosition;
376 this->phantomDepth = waterThickness;
377 this->phantomCenterDistance = phantomCenter;
378}
379void HadrontherapyAnalysisManager::setBeamMetaData(G4double meanKineticEnergy,G4double sigmaEnergy)
380{
381 this->beamEnergy = meanKineticEnergy;
382 this->energyError = sigmaEnergy;
383}
384/////////////////////////////////////////////////////////////////////////////
385// Flush data & close the file
386void HadrontherapyAnalysisManager::flush()
387{
388 if (theTFile)
389 {
390 theTFile -> Write();
391 theTFile -> Close();
392 }
393 eventCounter = 0;
394}
395
396#endif
Note: See TracBrowser for help on using the repository browser.