source: trunk/examples/extended/electromagnetic/TestEm7/src/RunAction.cc@ 1036

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

update
File size: 8.1 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: RunAction.cc,v 1.21 2008/01/14 12:11:39 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31
32#include "RunAction.hh"
33#include "DetectorConstruction.hh"
34#include "PhysicsList.hh"
35#include "StepMax.hh"
36#include "PrimaryGeneratorAction.hh"
37
38#include "G4Run.hh"
39#include "G4RunManager.hh"
40#include "G4UnitsTable.hh"
41#include "G4ios.hh"
42
43#include "Randomize.hh"
44
45#ifdef G4ANALYSIS_USE
46#include "AIDA/AIDA.h"
47#endif
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
52 PrimaryGeneratorAction* kin)
53:detector(det), physics(phys), kinematic(kin), af(0), tree(0)
54{
55 tallyEdep = new G4double[MaxTally];
56 binLength = offsetX = 0.;
57 histo[0] = 0;
58
59#ifdef G4ANALYSIS_USE
60 // Creating the analysis factory
61 af = AIDA_createAnalysisFactory();
62 if(!af) {
63 G4cout << "RunAction::RunAction() :"
64 << " problem creating the AIDA analysis factory."
65 << G4endl;
66 }
67#endif
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71
72RunAction::~RunAction()
73{
74 delete [] tallyEdep;
75
76#ifdef G4ANALYSIS_USE
77 delete af;
78#endif
79}
80
81//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83void RunAction::bookHisto()
84{
85 length = detector->GetAbsorSizeX();
86 G4double stepMax = physics->GetStepMaxProcess()->GetMaxStep();
87 const G4int nbmin = 100;
88 G4int nbBins = (int)(0.5 + length/stepMax);
89 if (nbBins < nbmin) nbBins = nbmin;
90 binLength = length/nbBins;
91 offsetX = 0.5*length;
92
93#ifdef G4ANALYSIS_USE
94 if (!af) return;
95
96 // Create a tree mapped to an hbook file.
97 G4bool readOnly = false;
98 G4bool createNew = true;
99 G4String options = "--noErrors uncompress";
100 AIDA::ITreeFactory* tf = af->createTreeFactory();
101 tree = tf->create("testem7.hbook","hbook", readOnly, createNew, options);
102 //tree = tf->create("testem7.root", "root",readOnly, createNew, options);
103 //tree = tf->create("testem7.XML" , "XML" ,readOnly, createNew, options);
104 delete tf;
105 if (!tree) {
106 G4cout << "RunAction::bookHisto()" << G4endl;
107 return;
108 }
109
110 // Create a histogram factory, whose histograms will be handled by the tree
111 AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree);
112
113 // Create histogram
114 histo[0] = hf->createHistogram1D("1","Edep (MeV/mm) along absorber (mm)",
115 nbBins, 0, length/mm);
116
117 delete hf;
118 G4cout << "\n----> Histogram Tree opened" << G4endl;
119#endif
120}
121
122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123
124void RunAction::cleanHisto()
125{
126#ifdef G4ANALYSIS_USE
127 tree->commit(); // Writing the histograms to the file
128 tree->close(); // and closing the tree (and the file)
129 delete tree;
130 tree = 0;
131
132 G4cout << "\n----> Histogram Tree saved" << G4endl;
133#endif
134}
135
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137
138void RunAction::FillHisto(G4int ih, G4double x, G4double weight)
139{
140#ifdef G4ANALYSIS_USE
141 if(histo[ih]) histo[ih]->fill(x, weight);
142#endif
143}
144
145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146
147void RunAction::BeginOfRunAction(const G4Run* aRun)
148{
149 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
150
151 // save Rndm status
152 G4RunManager::GetRunManager()->SetRandomNumberStore(true);
153 CLHEP::HepRandom::showEngineStatus();
154
155 //initialize projected range, tallies, Ebeam, and book histograms
156 //
157 nPrimarySteps = 0;
158 projRange = projRange2 = 0.;
159 edeptot = eniel = 0.;
160 for (G4int j=0; j<MaxTally; j++) tallyEdep[j] = 0.;
161 kinematic->ResetEbeamCumul();
162 bookHisto();
163}
164
165//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166
167void RunAction::EndOfRunAction(const G4Run* aRun)
168{
169 G4int NbofEvents = aRun->GetNumberOfEvent();
170 if (NbofEvents == 0) return;
171
172 //run conditions
173 //
174 G4Material* material = detector->GetAbsorMaterial();
175 G4double density = material->GetDensity();
176
177 G4String particle = kinematic->GetParticleGun()->GetParticleDefinition()
178 ->GetParticleName();
179 G4double energy = kinematic->GetParticleGun()->GetParticleEnergy();
180 G4cout << "\n The run consists of " << NbofEvents << " "<< particle << " of "
181 << G4BestUnit(energy,"Energy") << " through "
182 << G4BestUnit(detector->GetAbsorSizeX(),"Length") << " of "
183 << material->GetName() << " (density: "
184 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
185
186 //compute projected range and straggling
187 //
188 projRange /= NbofEvents; projRange2 /= NbofEvents;
189 G4double rms = projRange2 - projRange*projRange;
190 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
191
192 G4double nstep = G4double(nPrimarySteps)/G4double(NbofEvents);
193
194 G4cout.precision(6);
195 G4cout << "\n Projected Range= "<< G4BestUnit(projRange,"Length")
196 << " rms= " << G4BestUnit( rms,"Length")
197 << G4endl;
198 G4cout << " Mean number of primary steps = "<< nstep << G4endl;
199
200 //compute energy deposition and NIEL
201 //
202 edeptot /= NbofEvents;
203 G4cout << " Total energy deposit= "<< G4BestUnit(edeptot,"Energy")
204 << G4endl;
205 eniel /= NbofEvents;
206 G4cout << " NIEL energy deposit = "<< G4BestUnit(eniel,"Energy")
207 << G4endl;
208
209 //print dose in tallies
210 //
211 G4int tallyNumber = detector->GetTallyNumber();
212 if (tallyNumber > 0) {
213 G4double tallyMass = detector->GetTallyMass();
214 G4double Ebeam = kinematic->GetEbeamCumul();
215 G4cout << "\n---------------------------------------------------------\n";
216 G4cout << " Cumulated Doses : \tEdep \tEdep/Ebeam \tDose" << G4endl;
217 for (G4int j=0; j<tallyNumber; j++) {
218 G4double Edep = tallyEdep[j], ratio = 100*Edep/Ebeam;
219 G4double Dose = Edep/tallyMass;
220 G4cout << " tally " << j << ": \t \t"
221 << G4BestUnit(Edep,"Energy") << "\t"
222 << ratio << " % \t"
223 << G4BestUnit(Dose,"Dose") << G4endl;
224 }
225 G4cout << "\n---------------------------------------------------------\n";
226 G4cout << G4endl;
227 }
228
229#ifdef G4ANALYSIS_USE
230 // normalize histogram
231 G4double fac = (mm/MeV)/(NbofEvents * binLength);
232 histo[0]->scale(fac);
233#endif
234
235
236 // save and clean histo
237 cleanHisto();
238
239 // show Rndm status
240 CLHEP::HepRandom::showEngineStatus();
241}
242
243//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.