| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file electromagnetic/TestEm7/src/RunActio << 26 // $Id: RunAction.cc,v 1.20 2007/06/12 14:01:13 vnivanch Exp $
27 /// \brief Implementation of the RunAction cla << 27 // GEANT4 tag $Name: geant4-09-00 $
28 // << 28 //
29 // <<
30 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 31
33 #include "RunAction.hh" 32 #include "RunAction.hh"
34 <<
35 #include "DetectorConstruction.hh" 33 #include "DetectorConstruction.hh"
36 #include "PhysicsList.hh" 34 #include "PhysicsList.hh"
37 #include "PrimaryGeneratorAction.hh" <<
38 #include "StepMax.hh" 35 #include "StepMax.hh"
>> 36 #include "PrimaryGeneratorAction.hh"
39 37
40 #include "G4Run.hh" 38 #include "G4Run.hh"
41 #include "G4RunManager.hh" 39 #include "G4RunManager.hh"
42 #include "G4SystemOfUnits.hh" <<
43 #include "G4UnitsTable.hh" 40 #include "G4UnitsTable.hh"
44 #include "G4ios.hh" 41 #include "G4ios.hh"
>> 42
45 #include "Randomize.hh" 43 #include "Randomize.hh"
46 44
>> 45 #ifdef G4ANALYSIS_USE
>> 46 #include "AIDA/AIDA.h"
>> 47 #endif
>> 48
47 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 50
49 RunAction::RunAction(DetectorConstruction* det << 51 RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
50 : G4UserRunAction(), << 52 PrimaryGeneratorAction* kin)
51 fAnalysisManager(0), << 53 :detector(det), physics(phys), kinematic(kin), af(0), tree(0)
52 fDetector(det), << 54 {
53 fPhysics(phys), << 55 tallyEdep = new G4double[MaxTally];
54 fKinematic(kin), << 56 binLength = offsetX = 0.;
55 fTallyEdep(new G4double[kMaxTally]), << 57 histo[0] = 0;
56 fProjRange(0.), << 58
57 fProjRange2(0.), << 59 #ifdef G4ANALYSIS_USE
58 fEdeptot(0.), << 60 // Creating the analysis factory
59 fEniel(0.), << 61 af = AIDA_createAnalysisFactory();
60 fNbPrimarySteps(0), << 62 if(!af) {
61 fRange(0) << 63 G4cout << "RunAction::RunAction() :"
62 { << 64 << " problem creating the AIDA analysis factory."
63 // Book predefined histograms << 65 << G4endl;
64 BookHisto(); << 66 }
>> 67 #endif
65 } 68 }
66 69
67 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 71
69 RunAction::~RunAction() 72 RunAction::~RunAction()
70 { 73 {
71 delete[] fTallyEdep; << 74 delete [] tallyEdep;
>> 75
>> 76 #ifdef G4ANALYSIS_USE
>> 77 delete af;
>> 78 #endif
72 } 79 }
73 80
74 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 82
76 void RunAction::BeginOfRunAction(const G4Run* << 83 void RunAction::bookHisto()
77 { 84 {
78 G4cout << "### Run " << aRun->GetRunID() << << 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 }
79 121
80 if (!fAnalysisManager) { << 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81 BookHisto(); <<
82 } <<
83 123
84 CLHEP::HepRandom::showEngineStatus(); << 124 void 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 }
85 135
86 // initialize projected range, tallies, Ebea << 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 // <<
88 fNbPrimarySteps = 0; <<
89 fRange = 0; <<
90 fProjRange = fProjRange2 = 0.; <<
91 fEdeptot = fEniel = 0.; <<
92 for (G4int j = 0; j < kMaxTally; ++j) { <<
93 fTallyEdep[j] = 0.; <<
94 } <<
95 fKinematic->ResetEbeamCumul(); <<
96 137
97 if (fAnalysisManager->IsActive()) { << 138 void RunAction::FillHisto(G4int ih, G4double x, G4double weight)
98 fAnalysisManager->OpenFile(); << 139 {
>> 140 #ifdef G4ANALYSIS_USE
>> 141 if(histo[ih]) histo[ih]->fill(x, weight);
>> 142 #endif
>> 143 }
99 144
100 // histogram "1" is defined by the length << 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 // zoomed histograms are defined by UI com << 146
102 G4double length = fDetector->GetAbsorSizeX << 147 void RunAction::BeginOfRunAction(const G4Run* aRun)
103 G4double stepMax = fPhysics->GetStepMaxPro << 148 {
104 G4int nbBins = 100; << 149 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
105 if (stepMax < DBL_MAX) { << 150
106 G4int nb = (G4int)(0.5 + length / stepMa << 151 // save Rndm status
107 nbBins = std::min(std::max(nbBins, nb), << 152 G4RunManager::GetRunManager()->SetRandomNumberStore(true);
108 } << 153 CLHEP::HepRandom::showEngineStatus();
109 fAnalysisManager->SetH1(1, nbBins, 0., len << 154
110 } << 155 //initialize projected range, tallies, Ebeam, and book histograms
>> 156 //
>> 157 nPrimarySteps = 0;
>> 158 projRange = projRange2 = 0.;
>> 159 for (G4int j=0; j<MaxTally; j++) tallyEdep[j] = 0.;
>> 160 kinematic->ResetEbeamCumul();
>> 161 bookHisto();
111 } 162 }
112 163
113 //....oooOO0OOooo........oooOO0OOooo........oo 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 165
115 void RunAction::EndOfRunAction(const G4Run* aR 166 void RunAction::EndOfRunAction(const G4Run* aRun)
116 { 167 {
117 G4int nbofEvents = aRun->GetNumberOfEvent(); << 168 G4int NbofEvents = aRun->GetNumberOfEvent();
118 if (nbofEvents == 0) return; << 169 if (NbofEvents == 0) return;
119 170
120 // run conditions << 171 //run conditions
121 // << 172 //
122 const G4Material* material = fDetector->GetA << 173 G4Material* material = detector->GetAbsorMaterial();
123 G4double density = material->GetDensity(); 174 G4double density = material->GetDensity();
124 << 175
125 G4String particle = fKinematic->GetParticleG << 176 G4String particle = kinematic->GetParticleGun()->GetParticleDefinition()
126 G4double energy = fKinematic->GetParticleGun << 177 ->GetParticleName();
127 G4cout << "\n The run consists of " << nbofE << 178 G4double energy = kinematic->GetParticleGun()->GetParticleEnergy();
128 << G4BestUnit(energy, "Energy") << " << 179 G4cout << "\n The run consists of " << NbofEvents << " "<< particle << " of "
129 << G4BestUnit(fDetector->GetAbsorSize << 180 << G4BestUnit(energy,"Energy") << " through "
130 << " (density: " << G4BestUnit(densit << 181 << G4BestUnit(detector->GetAbsorSizeX(),"Length") << " of "
131 << 182 << material->GetName() << " (density: "
132 // compute projected range and straggling << 183 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
133 // << 184
134 if (fRange > 0) { << 185 //compute projected range and straggling
135 fProjRange /= fRange; << 186 //
136 fProjRange2 /= fRange; << 187 projRange /= NbofEvents; projRange2 /= NbofEvents;
137 } << 188 G4double rms = projRange2 - projRange*projRange;
138 G4double rms = fProjRange2 - fProjRange * fP << 189 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
139 if (rms > 0.) << 190
140 rms = std::sqrt(rms); << 191 G4double nstep = G4double(nPrimarySteps)/G4double(NbofEvents);
141 else << 192
142 rms = 0.; << 193 G4cout.precision(6);
143 << 194 G4cout << "\n projected Range= "<< G4BestUnit(projRange,"Length")
144 G4double nstep = G4double(fNbPrimarySteps) / << 195 << " rms= " << G4BestUnit( rms,"Length")
145 << 196 << G4endl;
146 G4cout.precision(6); << 197 G4cout << " mean number of primary steps = "<< nstep << G4endl;
147 G4cout << "\n Projected Range= " << G4BestUn << 198
148 << " rms= " << G4BestUnit(rms, "Len << 199 //print dose in tallies
149 G4cout << " Mean number of primary steps = " <<
150 <<
151 // compute energy deposition and niel <<
152 // <<
153 fEdeptot /= nbofEvents; <<
154 G4cout << " Total energy deposit= " << G4Bes <<
155 fEniel /= nbofEvents; <<
156 G4cout << " niel energy deposit = " << G4Bes <<
157 <<
158 // print dose in tallies <<
159 // 200 //
160 G4int tallyNumber = fDetector->GetTallyNumbe << 201 G4int tallyNumber = detector->GetTallyNumber();
161 if (tallyNumber > 0) { 202 if (tallyNumber > 0) {
162 G4double Ebeam = fKinematic->GetEbeamCumul << 203 G4double tallyMass = detector->GetTallyMass();
>> 204 G4double Ebeam = kinematic->GetEbeamCumul();
163 G4cout << "\n----------------------------- 205 G4cout << "\n---------------------------------------------------------\n";
164 G4cout << " Cumulated Doses : \tEdep 206 G4cout << " Cumulated Doses : \tEdep \tEdep/Ebeam \tDose" << G4endl;
165 for (G4int j = 0; j < tallyNumber; ++j) { << 207 for (G4int j=0; j<tallyNumber; j++) {
166 G4double Edep = fTallyEdep[j], ratio = 1 << 208 G4double Edep = tallyEdep[j], ratio = 100*Edep/Ebeam;
167 G4double tallyMass = fDetector->GetTally << 209 G4double Dose = Edep/tallyMass;
168 G4double Dose = Edep / tallyMass; << 210 G4cout << " tally " << j << ": \t \t"
169 G4cout << " tally " << j << ": \t \t" << << 211 << G4BestUnit(Edep,"Energy") << "\t"
170 << " % \t" << G4BestUnit(Dose, "D << 212 << ratio << " % \t"
>> 213 << G4BestUnit(Dose,"Dose") << G4endl;
171 } 214 }
172 G4cout << "\n----------------------------- 215 G4cout << "\n---------------------------------------------------------\n";
173 G4cout << G4endl; << 216 G4cout << G4endl;
174 } <<
175 <<
176 if (fAnalysisManager->IsActive()) { <<
177 // normalize histograms <<
178 // <<
179 for (G4int j = 1; j < 3; ++j) { <<
180 G4double binWidth = fAnalysisManager->Ge <<
181 G4double fac = (mm / MeV) / (nbofEvents <<
182 fAnalysisManager->ScaleH1(j, fac); <<
183 } <<
184 fAnalysisManager->ScaleH1(3, 1. / nbofEven <<
185 <<
186 // save histograms <<
187 fAnalysisManager->Write(); <<
188 fAnalysisManager->CloseFile(); <<
189 } 217 }
190 218
>> 219 #ifdef G4ANALYSIS_USE
>> 220 // normalize histogram
>> 221 G4double fac = (mm/MeV)/(NbofEvents * binLength);
>> 222 histo[0]->scale(fac);
>> 223 #endif
>> 224
>> 225
>> 226 // save and clean histo
>> 227 cleanHisto();
>> 228
191 // show Rndm status 229 // show Rndm status
192 // <<
193 CLHEP::HepRandom::showEngineStatus(); 230 CLHEP::HepRandom::showEngineStatus();
194 } <<
195 <<
196 //....oooOO0OOooo........oooOO0OOooo........oo <<
197 <<
198 void RunAction::BookHisto() <<
199 { <<
200 // Create or get analysis manager <<
201 // The choice of analysis technology is done <<
202 // in HistoManager.hh <<
203 fAnalysisManager = G4AnalysisManager::Instan <<
204 fAnalysisManager->SetDefaultFileType("root") <<
205 fAnalysisManager->SetFileName("testem7"); <<
206 fAnalysisManager->SetVerboseLevel(1); <<
207 fAnalysisManager->SetActivation(true); // e <<
208 <<
209 // Define histograms start values <<
210 const G4int kMaxHisto = 4; <<
211 const G4String id[] = {"h0", "h1", "h2", "h3 <<
212 const G4String title[] = { <<
213 "dummy", // 0 <<
214 "Edep (MeV/mm) along absorber ", // 1 <<
215 "Edep (MeV/mm) along absorber zoomed", // <<
216 "projectile range" // 3 <<
217 }; <<
218 <<
219 // Default values (to be reset via /analysis <<
220 G4int nbins = 100; <<
221 G4double vmin = 0.; <<
222 G4double vmax = 100.; <<
223 <<
224 // Create all histograms as inactivated <<
225 // as we have not yet set nbins, vmin, vmax <<
226 for (G4int k = 0; k < kMaxHisto; ++k) { <<
227 G4int ih = fAnalysisManager->CreateH1(id[k <<
228 G4bool activ = false; <<
229 if (k == 1) activ = true; <<
230 fAnalysisManager->SetH1Activation(ih, acti <<
231 } <<
232 } 231 }
233 232
234 //....oooOO0OOooo........oooOO0OOooo........oo 233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 234