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

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

geant4.9.4 beta rc0
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.