| 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 /// \file electromagnetic/TestEm7/src/RunAction.cc
27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class
28 // 28 //
29 // << 29 // $Id: RunAction.cc 101250 2016-11-10 08:54:02Z gcosmo $
>> 30 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 33
33 #include "RunAction.hh" 34 #include "RunAction.hh"
34 <<
35 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh"
36 #include "PhysicsList.hh" 36 #include "PhysicsList.hh"
37 #include "PrimaryGeneratorAction.hh" <<
38 #include "StepMax.hh" 37 #include "StepMax.hh"
>> 38 #include "PrimaryGeneratorAction.hh"
39 39
40 #include "G4Run.hh" 40 #include "G4Run.hh"
41 #include "G4RunManager.hh" 41 #include "G4RunManager.hh"
42 #include "G4SystemOfUnits.hh" <<
43 #include "G4UnitsTable.hh" 42 #include "G4UnitsTable.hh"
>> 43 #include "G4SystemOfUnits.hh"
44 #include "G4ios.hh" 44 #include "G4ios.hh"
>> 45
45 #include "Randomize.hh" 46 #include "Randomize.hh"
46 47
47 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 49
49 RunAction::RunAction(DetectorConstruction* det << 50 RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
50 : G4UserRunAction(), << 51 PrimaryGeneratorAction* kin)
51 fAnalysisManager(0), << 52 : G4UserRunAction(),
52 fDetector(det), << 53 fAnalysisManager(0), fDetector(det), fPhysics(phys), fKinematic(kin),
53 fPhysics(phys), << 54 fTallyEdep(new G4double[kMaxTally]), fProjRange(0.), fProjRange2(0.),
54 fKinematic(kin), << 55 fEdeptot(0.), fEniel(0.), fNbPrimarySteps(0), fRange(0)
55 fTallyEdep(new G4double[kMaxTally]), << 56 {
56 fProjRange(0.), <<
57 fProjRange2(0.), <<
58 fEdeptot(0.), <<
59 fEniel(0.), <<
60 fNbPrimarySteps(0), <<
61 fRange(0) <<
62 { <<
63 // Book predefined histograms 57 // Book predefined histograms
64 BookHisto(); << 58 BookHisto();
65 } 59 }
66 60
67 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 62
69 RunAction::~RunAction() 63 RunAction::~RunAction()
70 { 64 {
71 delete[] fTallyEdep; << 65 delete [] fTallyEdep;
72 } 66 }
73 67
74 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 69
76 void RunAction::BeginOfRunAction(const G4Run* 70 void RunAction::BeginOfRunAction(const G4Run* aRun)
77 { << 71 {
78 G4cout << "### Run " << aRun->GetRunID() << 72 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
79 73
80 if (!fAnalysisManager) { << 74 if(!fAnalysisManager) { BookHisto(); }
81 BookHisto(); << 75
82 } <<
83 <<
84 CLHEP::HepRandom::showEngineStatus(); 76 CLHEP::HepRandom::showEngineStatus();
85 << 77
86 // initialize projected range, tallies, Ebea << 78 //initialize projected range, tallies, Ebeam, and book histograms
87 // 79 //
88 fNbPrimarySteps = 0; 80 fNbPrimarySteps = 0;
89 fRange = 0; 81 fRange = 0;
90 fProjRange = fProjRange2 = 0.; 82 fProjRange = fProjRange2 = 0.;
91 fEdeptot = fEniel = 0.; 83 fEdeptot = fEniel = 0.;
92 for (G4int j = 0; j < kMaxTally; ++j) { << 84 for (G4int j=0; j<kMaxTally; ++j) { fTallyEdep[j] = 0.; }
93 fTallyEdep[j] = 0.; <<
94 } <<
95 fKinematic->ResetEbeamCumul(); 85 fKinematic->ResetEbeamCumul();
96 << 86
97 if (fAnalysisManager->IsActive()) { 87 if (fAnalysisManager->IsActive()) {
98 fAnalysisManager->OpenFile(); << 88 fAnalysisManager->OpenFile();
99 89
100 // histogram "1" is defined by the length 90 // histogram "1" is defined by the length of the target
101 // zoomed histograms are defined by UI com << 91 // zoomed histograms are defined by UI command
102 G4double length = fDetector->GetAbsorSizeX << 92 G4double length = fDetector->GetAbsorSizeX();
103 G4double stepMax = fPhysics->GetStepMaxPro 93 G4double stepMax = fPhysics->GetStepMaxProcess()->GetMaxStep();
104 G4int nbBins = 100; << 94 G4int nbmin = 100;
105 if (stepMax < DBL_MAX) { << 95 G4int nbBins = (G4int)(0.5 + length/stepMax);
106 G4int nb = (G4int)(0.5 + length / stepMa << 96 if (nbBins < nbmin) nbBins = nbmin;
107 nbBins = std::min(std::max(nbBins, nb), <<
108 } <<
109 fAnalysisManager->SetH1(1, nbBins, 0., len 97 fAnalysisManager->SetH1(1, nbBins, 0., length, "mm");
110 } 98 }
111 } 99 }
112 100
113 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 102
115 void RunAction::EndOfRunAction(const G4Run* aR 103 void RunAction::EndOfRunAction(const G4Run* aRun)
116 { 104 {
117 G4int nbofEvents = aRun->GetNumberOfEvent(); 105 G4int nbofEvents = aRun->GetNumberOfEvent();
118 if (nbofEvents == 0) return; 106 if (nbofEvents == 0) return;
119 107
120 // run conditions << 108 //run conditions
121 // << 109 //
122 const G4Material* material = fDetector->GetA 110 const G4Material* material = fDetector->GetAbsorMaterial();
123 G4double density = material->GetDensity(); 111 G4double density = material->GetDensity();
124 << 112
125 G4String particle = fKinematic->GetParticleG << 113 G4String particle = fKinematic->GetParticleGun()->GetParticleDefinition()
>> 114 ->GetParticleName();
126 G4double energy = fKinematic->GetParticleGun 115 G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
127 G4cout << "\n The run consists of " << nbofE << 116 G4cout << "\n The run consists of " << nbofEvents << " "<< particle << " of "
128 << G4BestUnit(energy, "Energy") << " << 117 << G4BestUnit(energy,"Energy") << " through "
129 << G4BestUnit(fDetector->GetAbsorSize << 118 << G4BestUnit(fDetector->GetAbsorSizeX(),"Length") << " of "
130 << " (density: " << G4BestUnit(densit << 119 << material->GetName() << " (density: "
131 << 120 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
132 // compute projected range and straggling << 121
>> 122 //compute projected range and straggling
133 // 123 //
134 if (fRange > 0) { << 124 if(fRange > 0) {
135 fProjRange /= fRange; << 125 fProjRange /= fRange;
136 fProjRange2 /= fRange; 126 fProjRange2 /= fRange;
137 } 127 }
138 G4double rms = fProjRange2 - fProjRange * fP << 128 G4double rms = fProjRange2 - fProjRange*fProjRange;
139 if (rms > 0.) << 129 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
140 rms = std::sqrt(rms); <<
141 else <<
142 rms = 0.; <<
143 <<
144 G4double nstep = G4double(fNbPrimarySteps) / <<
145 <<
146 G4cout.precision(6); <<
147 G4cout << "\n Projected Range= " << G4BestUn <<
148 << " rms= " << G4BestUnit(rms, "Len <<
149 G4cout << " Mean number of primary steps = " <<
150 130
151 // compute energy deposition and niel << 131 G4double nstep = G4double(fNbPrimarySteps)/G4double(nbofEvents);
152 // <<
153 fEdeptot /= nbofEvents; <<
154 G4cout << " Total energy deposit= " << G4Bes <<
155 fEniel /= nbofEvents; <<
156 G4cout << " niel energy deposit = " << G4Bes <<
157 132
158 // print dose in tallies << 133 G4cout.precision(6);
>> 134 G4cout << "\n Projected Range= "<< G4BestUnit(fProjRange,"Length")
>> 135 << " rms= " << G4BestUnit( rms,"Length")
>> 136 << G4endl;
>> 137 G4cout << " Mean number of primary steps = "<< nstep << G4endl;
>> 138
>> 139 //compute energy deposition and niel
>> 140 //
>> 141 fEdeptot /= nbofEvents;
>> 142 G4cout << " Total energy deposit= "<< G4BestUnit(fEdeptot,"Energy")
>> 143 << G4endl;
>> 144 fEniel /= nbofEvents;
>> 145 G4cout << " niel energy deposit = "<< G4BestUnit(fEniel,"Energy")
>> 146 << G4endl;
>> 147
>> 148 //print dose in tallies
159 // 149 //
160 G4int tallyNumber = fDetector->GetTallyNumbe 150 G4int tallyNumber = fDetector->GetTallyNumber();
161 if (tallyNumber > 0) { 151 if (tallyNumber > 0) {
162 G4double Ebeam = fKinematic->GetEbeamCumul 152 G4double Ebeam = fKinematic->GetEbeamCumul();
163 G4cout << "\n----------------------------- 153 G4cout << "\n---------------------------------------------------------\n";
164 G4cout << " Cumulated Doses : \tEdep 154 G4cout << " Cumulated Doses : \tEdep \tEdep/Ebeam \tDose" << G4endl;
165 for (G4int j = 0; j < tallyNumber; ++j) { << 155 for (G4int j=0; j < tallyNumber; ++j) {
166 G4double Edep = fTallyEdep[j], ratio = 1 << 156 G4double Edep = fTallyEdep[j], ratio = 100*Edep/Ebeam;
167 G4double tallyMass = fDetector->GetTally << 157 G4double tallyMass = fDetector->GetTallyMass(j);
168 G4double Dose = Edep / tallyMass; << 158 G4double Dose = Edep/tallyMass;
169 G4cout << " tally " << j << ": \t \t" << << 159 G4cout << " tally " << j << ": \t \t"
170 << " % \t" << G4BestUnit(Dose, "D << 160 << G4BestUnit(Edep,"Energy") << "\t"
>> 161 << ratio << " % \t"
>> 162 << G4BestUnit(Dose,"Dose") << G4endl;
171 } 163 }
172 G4cout << "\n----------------------------- 164 G4cout << "\n---------------------------------------------------------\n";
173 G4cout << G4endl; << 165 G4cout << G4endl;
174 } 166 }
175 167
176 if (fAnalysisManager->IsActive()) { << 168 if (fAnalysisManager->IsActive() ) {
177 // normalize histograms 169 // normalize histograms
178 // 170 //
179 for (G4int j = 1; j < 3; ++j) { << 171 for (G4int j=1; j<3; ++j) {
180 G4double binWidth = fAnalysisManager->Ge 172 G4double binWidth = fAnalysisManager->GetH1Width(j);
181 G4double fac = (mm / MeV) / (nbofEvents << 173 G4double fac = (mm/MeV)/(nbofEvents * binWidth);
182 fAnalysisManager->ScaleH1(j, fac); 174 fAnalysisManager->ScaleH1(j, fac);
183 } 175 }
184 fAnalysisManager->ScaleH1(3, 1. / nbofEven << 176 fAnalysisManager->ScaleH1(3, 1./nbofEvents);
185 << 177
186 // save histograms 178 // save histograms
187 fAnalysisManager->Write(); 179 fAnalysisManager->Write();
188 fAnalysisManager->CloseFile(); 180 fAnalysisManager->CloseFile();
>> 181 delete fAnalysisManager;
>> 182 fAnalysisManager = 0;
189 } 183 }
190 << 184
191 // show Rndm status 185 // show Rndm status
192 // 186 //
193 CLHEP::HepRandom::showEngineStatus(); 187 CLHEP::HepRandom::showEngineStatus();
194 } 188 }
195 189
196 //....oooOO0OOooo........oooOO0OOooo........oo 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
197 191
198 void RunAction::BookHisto() 192 void RunAction::BookHisto()
199 { 193 {
200 // Create or get analysis manager 194 // Create or get analysis manager
201 // The choice of analysis technology is done 195 // The choice of analysis technology is done via selection of a namespace
202 // in HistoManager.hh 196 // in HistoManager.hh
203 fAnalysisManager = G4AnalysisManager::Instan 197 fAnalysisManager = G4AnalysisManager::Instance();
204 fAnalysisManager->SetDefaultFileType("root") <<
205 fAnalysisManager->SetFileName("testem7"); 198 fAnalysisManager->SetFileName("testem7");
206 fAnalysisManager->SetVerboseLevel(1); 199 fAnalysisManager->SetVerboseLevel(1);
207 fAnalysisManager->SetActivation(true); // e 200 fAnalysisManager->SetActivation(true); // enable inactivation of histograms
208 201
209 // Define histograms start values 202 // Define histograms start values
210 const G4int kMaxHisto = 4; 203 const G4int kMaxHisto = 4;
211 const G4String id[] = {"h0", "h1", "h2", "h3 << 204 const G4String id[] = { "h0", "h1", "h2", "h3" };
212 const G4String title[] = { << 205 const G4String title[] =
213 "dummy", // 0 << 206 { "dummy", //0
214 "Edep (MeV/mm) along absorber ", // 1 << 207 "Edep (MeV/mm) along absorber ", //1
215 "Edep (MeV/mm) along absorber zoomed", // << 208 "Edep (MeV/mm) along absorber zoomed", //2
216 "projectile range" // 3 << 209 "projectile range" //3
217 }; << 210 };
218 211
219 // Default values (to be reset via /analysis << 212 // Default values (to be reset via /analysis/h1/set command)
220 G4int nbins = 100; 213 G4int nbins = 100;
221 G4double vmin = 0.; 214 G4double vmin = 0.;
222 G4double vmax = 100.; 215 G4double vmax = 100.;
223 216
224 // Create all histograms as inactivated << 217 // Create all histograms as inactivated
225 // as we have not yet set nbins, vmin, vmax 218 // as we have not yet set nbins, vmin, vmax
226 for (G4int k = 0; k < kMaxHisto; ++k) { << 219 for (G4int k=0; k<kMaxHisto; ++k) {
227 G4int ih = fAnalysisManager->CreateH1(id[k 220 G4int ih = fAnalysisManager->CreateH1(id[k], title[k], nbins, vmin, vmax);
228 G4bool activ = false; 221 G4bool activ = false;
229 if (k == 1) activ = true; 222 if (k == 1) activ = true;
230 fAnalysisManager->SetH1Activation(ih, acti 223 fAnalysisManager->SetH1Activation(ih, activ);
231 } 224 }
232 } 225 }
233 226
234 //....oooOO0OOooo........oooOO0OOooo........oo 227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 228