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

Last change on this file since 1275 was 1230, checked in by garnier, 16 years ago

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