source: trunk/examples/advanced/radioprotection/src/RemSimAnalysisManager.cc@ 1219

Last change on this file since 1219 was 807, checked in by garnier, 17 years ago

update
File size: 19.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// * RemSimAnalysisManager.cc *
29// * *
30// *******************************
31//
32// $Id: RemSimAnalysisManager.cc,v 1.12 2006/06/29 16:23:31 gunter Exp $
33//
34// Author:Susanna Guatelli, guatelli@ge.infn.it
35//
36#ifdef G4ANALYSIS_USE
37#include <stdlib.h>
38#include <fstream>
39#include "RemSimAnalysisManager.hh"
40#include "G4ios.hh"
41#include <AIDA/AIDA.h>
42#include "G4RunManager.hh"
43#include <vector>
44
45RemSimAnalysisManager* RemSimAnalysisManager::instance = 0;
46
47RemSimAnalysisManager::RemSimAnalysisManager()
48 : aFact(0), treeFact(0), theTree(0), dataPointFactory(0),
49 histogramFactory(0), dataPoint(0), energyDeposit(0),
50 primary(0), secondaryDeposit(0),
51 primaryInitialE(0), primaryInitialEout(0), initialE(0),
52 initialEout(0), shape(0), energyShape(0) , histo_secondary_phantom(0),histo_secondary(0),
53 histo_proton(0), histo_neutron(0), histo_pion(0), histo_alpha(0),
54 histo_positron(0), histo_electron(0), histo_gamma(0), histo_muon(0), histo_other(0), histo_neutrino(),
55 histo_proton_reaching(0), histo_neutron_reaching(0), histo_pion_reaching(0), histo_alpha_reaching(0),
56 histo_positron_reaching(0), histo_electron_reaching(0), histo_gamma_reaching(0), histo_muon_reaching(0),
57 histo_other_reaching(0),
58 histo_proton_slice(0), histo_neutron_slice(0), histo_pion_slice(0), histo_alpha_slice(0),
59 histo_positron_slice(0), histo_electron_slice(0), histo_gamma_slice(0),
60 histo_muon_slice(0), histo_other_slice(0), histo_vehicle(0)
61{
62 aFact = AIDA_createAnalysisFactory();
63 fileFormat = "hbook";
64}
65
66RemSimAnalysisManager::~RemSimAnalysisManager()
67{
68 delete histo_vehicle; histo_vehicle = 0;
69 delete histo_proton_reaching; histo_proton_reaching = 0;
70 delete histo_neutron_reaching; histo_neutron_reaching = 0;
71 delete histo_pion_reaching; histo_pion_reaching = 0;
72 delete histo_alpha_reaching; histo_alpha_reaching = 0;
73 delete histo_positron_reaching; histo_positron_reaching = 0;
74 delete histo_electron_reaching; histo_electron_reaching = 0;
75 delete histo_gamma_reaching; histo_gamma_reaching = 0;
76 delete histo_muon_reaching; histo_muon_reaching = 0;
77 delete histo_other_reaching; histo_other_reaching = 0;
78 delete histo_proton_slice; histo_proton_slice = 0;
79 delete histo_neutron_slice; histo_neutron_slice = 0;
80 delete histo_pion_slice; histo_pion_slice = 0;
81 delete histo_alpha_slice; histo_alpha_slice = 0;
82 delete histo_positron_slice; histo_positron_slice = 0;
83 delete histo_electron_slice; histo_electron_slice = 0;
84 delete histo_gamma_slice; histo_gamma_slice = 0;
85 delete histo_muon_slice; histo_muon_slice = 0;
86 delete histo_other_slice; histo_other_slice = 0;
87 delete histo_proton; histo_proton = 0;
88 delete histo_neutron; histo_neutron = 0;
89 delete histo_pion; histo_pion = 0;
90 delete histo_alpha; histo_alpha =0;
91 delete histo_positron; histo_positron = 0;
92 delete histo_electron; histo_electron = 0;
93 delete histo_gamma; histo_gamma = 0;
94 delete histo_muon; histo_muon = 0;
95 delete histo_other; histo_other = 0;
96 delete histo_neutrino; histo_neutrino = 0;
97 delete histo_secondary; histo_secondary = 0;
98 delete histo_secondary_phantom; histo_secondary_phantom = 0;
99 delete energyShape; energyShape = 0;
100 delete shape; shape = 0;
101 delete initialEout; initialEout = 0;
102 delete initialE; initialE = 0;
103 delete primaryInitialEout; primaryInitialEout = 0;
104 delete primaryInitialE; primaryInitialE = 0;
105 delete secondaryDeposit; secondaryDeposit = 0;
106 delete primary; primary = 0;
107 delete energyDeposit; energyDeposit = 0;
108 delete dataPoint; dataPoint = 0;
109 delete histogramFactory; histogramFactory = 0;
110 delete dataPointFactory; dataPointFactory = 0;
111 delete theTree; theTree = 0;
112 delete treeFact; treeFact = 0;
113 delete aFact; aFact = 0;
114}
115
116RemSimAnalysisManager* RemSimAnalysisManager::getInstance()
117{
118 if (instance == 0) instance = new RemSimAnalysisManager;
119 return instance;
120}
121
122void RemSimAnalysisManager::book()
123{
124 treeFact = aFact -> createTreeFactory();
125
126
127 theTree = treeFact -> create("remsim.hbk", "hbook", false, true);
128 G4cout << "The format of the output file is hbook" << G4endl;
129
130 histogramFactory = aFact -> createHistogramFactory(*theTree);
131 energyDeposit = histogramFactory -> createHistogram1D("10",
132 "Energy Deposit",
133 30,// number of bins
134 0.,//xmin
135 30.);//xmax
136
137 /// A variable binning is created for some histograms
138 G4int size = 1901;
139
140 std::vector<double> vector_bin;
141
142 G4double x_value = 0.;
143 for (G4int ptn = 0; ptn < size; ptn++ )
144 {
145 G4int binNb = 1000;
146 if (ptn < binNb)
147 {
148 vector_bin.push_back(x_value);
149 x_value += 10;
150 }
151 else
152 { vector_bin.push_back(x_value); x_value += 100.;}
153 }
154
155 primary = histogramFactory -> createHistogram1D("20",
156 "Initial energy of primary particles",
157 vector_bin);
158
159 secondaryDeposit = histogramFactory -> createHistogram1D("30",
160 "EnergyDeposit given by secondaries",
161 30,0.,30.);
162
163 primaryInitialE = histogramFactory -> createHistogram1D("40",
164 "Initial energy of primaries reaching the phantom", vector_bin);
165
166 primaryInitialEout = histogramFactory -> createHistogram1D("50",
167 "Initial energy of primaries ougoing the phantom",
168 vector_bin);
169
170 initialE = histogramFactory -> createHistogram1D("60",
171 "Energy of primaries reaching the phantom",
172 vector_bin);
173
174 initialEout = histogramFactory -> createHistogram1D("70",
175 "Energy of primaries outgoing the phantom",
176 vector_bin);
177
178 shape = histogramFactory -> createHistogram2D("80",
179 "Shape",
180 320,-160.,160.,
181 320, -160.,160.);
182
183 energyShape = histogramFactory -> createHistogram2D("90",
184 "energyDepShape",
185 320, -160.,160.,
186 320, -160.,160.);
187
188 histo_secondary_phantom = histogramFactory -> createHistogram1D("100",
189 "secondary particles produced in the phantom",
190 10, 0., 10.);
191
192 histo_secondary = histogramFactory -> createHistogram1D("300",
193 "secondary particles reaching the phantom",
194 10, 0., 10.);
195
196 histo_vehicle = histogramFactory -> createHistogram1D("500",
197 "secondary particles in the vehicle",
198 10, 0., 10.);
199
200 histo_proton = histogramFactory -> createHistogram1D("110","Energy of secondary p produced in the phantom",
201 100, 0., 1000.); // up to 1 GeV, 10 MeV bin
202
203 histo_neutron= histogramFactory -> createHistogram1D("120","Energy of secondary n produced in the phantom",
204 100, 0., 1000.); // up to 1 GeV, 10 MeV bin
205
206 histo_pion = histogramFactory-> createHistogram1D("130","Energy of secondary pions produced in the phantom",
207 200, 0., 2000.);// up to 2 GeV, 10 MeV bin
208
209
210 histo_alpha = histogramFactory-> createHistogram1D("140","Energy of secondary alpha produced in the phantom",
211 100, 0.,100.); // up to 100 MeV, 1 MeV bin
212
213 histo_positron = histogramFactory-> createHistogram1D("150","Energy of secondary e+ produced in the phantom",
214 100, 0., 1000.); // up to 1 GeV, 10 MeV bin
215
216 histo_electron = histogramFactory-> createHistogram1D("160","Energy of secondary e- produced in the phantom",
217 1000, 0., 10.); // up to 10 MeV, 10. keV bin
218
219 histo_gamma = histogramFactory-> createHistogram1D("170","Energy of secondary gamma produced in the phantom",
220 100, 0., 10.); // up to 10 MeV, 1 MeV bin
221
222 histo_muon = histogramFactory-> createHistogram1D("180","Energy of secondary mu produced in the phantom",
223 100, 0., 1000.); // up to 1 GeV, 10 MeV bin
224
225 histo_other= histogramFactory -> createHistogram1D("190","Energy of secondary other particles produced in the phantom",
226 200,0., 200.); // up to 200 MeV, 1 MeV bin
227
228 histo_neutrino = histogramFactory -> createHistogram1D("400", "Energy of secondary neutrinos produced in the phantom", 1000, 0., 10000.);
229
230 histo_proton_slice = histogramFactory -> createHistogram1D("200",
231 " Phantom Slice where secondary protons are produced",
232 30, 0., 30.);
233
234 histo_neutron_slice = histogramFactory -> createHistogram1D("210",
235 "Phantom Slice where secondary neutrons are produced",
236 30, 0., 30.);
237
238 histo_pion_slice = histogramFactory -> createHistogram1D("220",
239 "Phantom Slice where secondary pions are produced",
240 30, 0., 30.);
241
242 histo_alpha_slice = histogramFactory -> createHistogram1D("230",
243 "Phantom Slice where secondary alpha are produced",
244 30, 0., 30.);
245
246 histo_positron_slice = histogramFactory -> createHistogram1D("240",
247 "Phantom Slice where secondary positrons are produced ",
248 30, 0.,30.);
249
250 histo_electron_slice = histogramFactory -> createHistogram1D("250",
251 "Phantom Slice where secondary electrons are produced",
252 30, 0., 30.);
253
254 histo_gamma_slice = histogramFactory -> createHistogram1D("260",
255 "Phantom Slice where secondary gamma are produced",
256 30, 0., 30.);
257
258 histo_muon_slice = histogramFactory -> createHistogram1D("270",
259 "Phantom Slice where secondary muons are produced",
260 30, 0., 30.);
261
262 histo_other_slice = histogramFactory -> createHistogram1D("280",
263 " Phantom Slice where secondary other particles are produced",
264 30,0.,30.);
265
266 histo_proton_reaching = histogramFactory -> createHistogram1D("310","Energy of secondary p reaching the phantom",
267 1000, 0., 10000.); // up to 10 GeV, 10 MeV bin
268
269 histo_neutron_reaching= histogramFactory -> createHistogram1D("320","Energy of secondary n reaching the phantom",
270 1000, 0., 10000.); // up to 10 GeV, 10 MeV bin
271
272 histo_pion_reaching = histogramFactory -> createHistogram1D("330","Energy of secondary pions reaching in the phantom",
273 1000, 0., 10000.); // up to 10 GeV, 10 MeV bin
274
275 histo_alpha_reaching = histogramFactory -> createHistogram1D("340","Energy of secondary alpha reaching the phantom",
276 100, 0., 100.); // up to 100 MeV, 1 MeV bin
277
278 histo_positron_reaching = histogramFactory -> createHistogram1D("350","Energy of secondary e+ reaching the phantom",
279 500, 0., 5000.); // up to 5 GeV, 10 MeV bin
280
281 histo_electron_reaching = histogramFactory -> createHistogram1D("360","Energy of secondary e- reaching the phantom",
282 200, 0., 2000.); //up to 2 GeV, 10 MeV bin
283
284 histo_gamma_reaching = histogramFactory -> createHistogram1D("370","Energy of secondary gamma reaching the phantom",
285 200, 0., 2000.);//up to 2 GeV, 10 MeV bin
286
287 histo_muon_reaching = histogramFactory -> createHistogram1D("380","Energy of secondary mu reaching the phantom",
288 200, 0., 2000.); // up to 2 GeV, 10 MeV bin
289
290 histo_other_reaching= histogramFactory -> createHistogram1D("390","Energy of secondary other particles reaching the phantom",
291 100, 0., 1000.); // up to 1 GeV, 10 MeV bin
292}
293
294void RemSimAnalysisManager::energyDepositStore(G4int layer, G4double eDep)
295{
296 energyDeposit -> fill(layer,eDep);
297}
298
299void RemSimAnalysisManager::primaryParticleEnergyDistribution(G4double energy)
300{
301 primary -> fill(energy);
302}
303
304void RemSimAnalysisManager::SecondaryEnergyDeposit(G4int i, G4double EnergyDep)
305{
306 secondaryDeposit -> fill(i,EnergyDep);
307}
308
309void RemSimAnalysisManager::PrimaryInitialEnergyIn(G4double EnergyDep)
310{
311 primaryInitialE -> fill(EnergyDep);
312}
313
314void RemSimAnalysisManager::PrimaryInitialEnergyOut(G4double EnergyDep)
315{
316 primaryInitialEout -> fill(EnergyDep);
317}
318
319void RemSimAnalysisManager::PrimaryEnergyIn(G4double EnergyDep)
320{
321 initialE -> fill(EnergyDep);
322}
323
324void RemSimAnalysisManager::PrimaryEnergyOut(G4double EnergyDep)
325{
326 initialEout -> fill(EnergyDep);
327}
328
329void RemSimAnalysisManager::particleShape(G4double x, G4double y)
330{
331 shape -> fill(x,y);
332}
333
334void RemSimAnalysisManager::energyDepShape(G4double x, G4double y, G4double energyDep)
335{
336 energyShape -> fill(x,y, energyDep);
337}
338
339void RemSimAnalysisManager:: SetFormat(G4String format)
340{
341 // fileFormat = format;
342
343 if (fileFormat != "hbook")
344 G4cout << format << "is not available at the moment !!!" << G4endl;
345}
346void RemSimAnalysisManager:: SecondaryInPhantom(G4int i)
347{
348 histo_secondary_phantom -> fill(i);
349
350 // if i = 0 secondary proton, i = 1 neutron, i= 2 pion, i= 3 alpha,
351 // i =4 other, i=5 electron, i = 6 gamma, i= 7 positrons,
352 // i = 8 muons, i= 9 neutrinos
353}
354
355void RemSimAnalysisManager::SecondaryProtonInPhantom(G4double energy)
356{
357 histo_proton -> fill(energy);
358 //G4cout<< "proton energy : "<< energy << G4endl;
359}
360
361void RemSimAnalysisManager::SecondaryNeutronInPhantom(G4double energy)
362{
363 histo_neutron -> fill(energy);
364 //G4cout<< "neutron energy : "<< energy << G4endl;
365}
366
367void RemSimAnalysisManager::SecondaryPionInPhantom(G4double energy)
368{
369 histo_pion -> fill(energy);
370 //G4cout<< "pion energy : "<< energy << G4endl;
371}
372
373void RemSimAnalysisManager::SecondaryAlphaInPhantom(G4double energy)
374{
375 histo_alpha -> fill(energy);
376 //G4cout<< "alpha energy : "<< energy << G4endl;
377}
378
379void RemSimAnalysisManager::SecondaryPositronInPhantom(G4double energy)
380{
381 histo_positron -> fill(energy);
382 //G4cout<< "positron energy : "<< energy << G4endl;
383}
384
385void RemSimAnalysisManager::SecondaryElectronInPhantom(G4double energy)
386{
387 histo_electron -> fill(energy);
388 //G4cout<< "electron energy : "<< energy << G4endl;
389}
390
391void RemSimAnalysisManager::SecondaryGammaInPhantom(G4double energy)
392{
393 histo_gamma -> fill(energy);
394 //G4cout<< "gamma energy : "<< energy << G4endl;
395}
396
397void RemSimAnalysisManager::SecondaryMuonInPhantom(G4double energy)
398{
399 histo_muon -> fill(energy);
400 //G4cout<< "muon energy : "<< energy << G4endl;
401}
402
403void RemSimAnalysisManager::SecondaryOtherInPhantom(G4double energy)
404{
405 histo_other -> fill(energy);
406 //G4cout<< "other energy : "<< energy << G4endl;
407}
408void RemSimAnalysisManager::SecondaryNeutrinoInPhantom (G4double energy)
409{
410 histo_neutrino -> fill(energy);
411 // G4cout<< "neutrino energy:" << energy << G4endl;
412}
413
414void RemSimAnalysisManager::SecondaryProtonInPhantomSlice(G4double slice)
415{
416 histo_proton_slice -> fill(slice);
417 //G4cout<< "proton slice : "<< slice << G4endl;
418}
419
420void RemSimAnalysisManager::SecondaryNeutronInPhantomSlice(G4double slice)
421{
422 histo_neutron_slice -> fill(slice);
423 //G4cout<< "neutron slice : "<< slice << G4endl;
424}
425
426void RemSimAnalysisManager::SecondaryPionInPhantomSlice(G4double slice)
427{
428 histo_pion_slice -> fill(slice);
429 //G4cout<< "pion slice : "<< slice << G4endl;
430}
431
432void RemSimAnalysisManager::SecondaryAlphaInPhantomSlice(G4double slice)
433{
434 histo_alpha_slice -> fill(slice);
435 //G4cout<< "alpha slice : "<< slice << G4endl;
436}
437
438void RemSimAnalysisManager::SecondaryPositronInPhantomSlice(G4double slice)
439{
440 histo_positron_slice -> fill(slice);
441 //G4cout<< "positron slice : "<< slice << G4endl;
442}
443
444void RemSimAnalysisManager::SecondaryElectronInPhantomSlice(G4double slice)
445{
446 histo_electron_slice -> fill(slice);
447 //G4cout<< "electron slice : "<< slice << G4endl;
448}
449
450void RemSimAnalysisManager::SecondaryGammaInPhantomSlice(G4double slice)
451{
452 histo_gamma_slice -> fill(slice);
453 //G4cout<< "gamma slice : "<< slice << G4endl;
454}
455
456void RemSimAnalysisManager::SecondaryMuonInPhantomSlice(G4double slice)
457{
458 histo_muon_slice -> fill(slice);
459 //G4cout<< "muon slice : "<< slice << G4endl;
460}
461
462void RemSimAnalysisManager::SecondaryOtherInPhantomSlice(G4double slice)
463{
464 histo_other_slice -> fill(slice);
465 //G4cout<< "other slice : "<< slice << G4endl;
466}
467
468void RemSimAnalysisManager::SecondaryReachingThePhantom(G4int i)
469{
470 histo_secondary -> fill(i);
471 // if i =0 secondary proton, i =1 neutron, i=2 pion, i=3 alpha,
472 // i =4 other, i=5 electron, i = 6 gamma, i=7 positrons,
473 // i=8 muons, i=9 neutrinos
474}
475
476void RemSimAnalysisManager::SecondaryProtonReachingThePhantom(G4double energy)
477{
478 histo_proton_reaching -> fill(energy);
479 //G4cout<< "proton energy reaching the phantom: "<< energy << G4endl;
480}
481
482void RemSimAnalysisManager::SecondaryNeutronReachingThePhantom(G4double energy)
483{
484 histo_neutron_reaching -> fill(energy);
485 //G4cout<< "neutron energy reaching the phantom: "<< energy << G4endl;
486}
487
488void RemSimAnalysisManager::SecondaryPionReachingThePhantom(G4double energy)
489{
490 histo_pion_reaching -> fill(energy);
491 //G4cout<< "pion energy reaching the phantom: "<< energy << G4endl;
492}
493
494void RemSimAnalysisManager::SecondaryAlphaReachingThePhantom(G4double energy)
495{
496 histo_alpha_reaching -> fill(energy);
497 //G4cout<< "alpha energy reaching the phantom: "<< energy << G4endl;
498}
499
500void RemSimAnalysisManager::SecondaryPositronReachingThePhantom(G4double energy)
501{
502 histo_positron_reaching -> fill(energy);
503 //G4cout<< "positron energy reaching the phantom: "<< energy << G4endl;
504}
505
506void RemSimAnalysisManager::SecondaryElectronReachingThePhantom(G4double energy)
507{
508 histo_electron_reaching -> fill(energy);
509 //G4cout<< "electron energy reaching the phantom: "<< energy << G4endl;
510}
511
512void RemSimAnalysisManager::SecondaryGammaReachingThePhantom(G4double energy)
513{
514 histo_gamma_reaching -> fill(energy);
515 //G4cout<< "gamma energy reaching the phantom: "<< energy << G4endl;
516}
517
518void RemSimAnalysisManager::SecondaryMuonReachingThePhantom(G4double energy)
519{
520 histo_muon_reaching -> fill(energy);
521 //G4cout<< "muon energy reaching the phantom: "<< energy << G4endl;
522}
523
524void RemSimAnalysisManager::SecondaryOtherReachingThePhantom(G4double energy)
525{
526 histo_other_reaching -> fill(energy);
527 //G4cout<< "other energy reaching the phantom: "<< energy << G4endl;
528}
529
530void RemSimAnalysisManager::SecondaryInVehicle(G4int i)
531{
532 histo_vehicle -> fill(i);
533}
534
535void RemSimAnalysisManager::finish()
536{
537 theTree -> commit();
538 theTree -> close();
539}
540#endif
541
542
543
544
545
546
547
548
549
550
551
Note: See TracBrowser for help on using the repository browser.