| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file electromagnetic/TestEm15/src/RunActi
27 /// \brief Implementation of the RunAction cla
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
32
33 #include "RunAction.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh"
38
39 #include "G4EmCalculator.hh"
40 #include "G4Run.hh"
41 #include "G4RunManager.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UnitsTable.hh"
44 #include "Randomize.hh"
45
46 #include <iomanip>
47
48 //....oooOO0OOooo........oooOO0OOooo........oo
49
50 RunAction::RunAction(DetectorConstruction* det
51 : G4UserRunAction(), fDetector(det), fPrimar
52 {
53 fHistoManager = new HistoManager();
54 }
55
56 //....oooOO0OOooo........oooOO0OOooo........oo
57
58 RunAction::~RunAction()
59 {
60 delete fHistoManager;
61 }
62
63 //....oooOO0OOooo........oooOO0OOooo........oo
64
65 void RunAction::BeginOfRunAction(const G4Run*)
66 {
67 // save Rndm status
68 ////G4RunManager::GetRunManager()->SetRandom
69 CLHEP::HepRandom::showEngineStatus();
70
71 fProcCounter = new ProcessesCount;
72 fTotalCount = 0;
73
74 fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0.
75 fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0.
76 fTetPrj = fTetPrj2 = 0.;
77 fPhiCor = fPhiCor2 = 0.;
78
79 // histograms
80 //
81 G4AnalysisManager* analysisManager = G4Analy
82 if (analysisManager->IsActive()) {
83 analysisManager->OpenFile();
84 }
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oo
88
89 void RunAction::CountProcesses(G4String procNa
90 {
91 // does the process already encounted ?
92 size_t nbProc = fProcCounter->size();
93 size_t i = 0;
94 while ((i < nbProc) && ((*fProcCounter)[i]->
95 i++;
96 if (i == nbProc) fProcCounter->push_back(new
97
98 (*fProcCounter)[i]->Count();
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oo
102
103 void RunAction::EndOfRunAction(const G4Run* aR
104 {
105 G4int NbOfEvents = aRun->GetNumberOfEvent();
106 if (NbOfEvents == 0) return;
107
108 G4int prec = G4cout.precision(5);
109
110 G4Material* material = fDetector->GetMateria
111 G4double density = material->GetDensity();
112
113 G4ParticleDefinition* particle = fPrimary->G
114 G4String Particle = particle->GetParticleNam
115 G4double energy = fPrimary->GetParticleGun()
116 G4cout << "\n The run consists of " << NbOfE
117 << G4BestUnit(energy, "Energy") << "
118 << G4BestUnit(fDetector->GetBoxSize()
119 << " (density: " << G4BestUnit(densit
120
121 // frequency of processes
122 G4cout << "\n Process calls frequency --->";
123 for (size_t i = 0; i < fProcCounter->size();
124 G4String procName = (*fProcCounter)[i]->Ge
125 G4int count = (*fProcCounter)[i]->GetCount
126 G4cout << "\t" << procName << " = " << cou
127 }
128
129 if (fTotalCount > 0) {
130 // compute path length and related quantit
131 //
132 G4double MeanTPL = fTruePL / fTotalCount;
133 G4double MeanTPL2 = fTruePL2 / fTotalCount
134 G4double rmsTPL = std::sqrt(std::fabs(Mean
135
136 G4double MeanGPL = fGeomPL / fTotalCount;
137 G4double MeanGPL2 = fGeomPL2 / fTotalCount
138 G4double rmsGPL = std::sqrt(std::fabs(Mean
139
140 G4double MeanLaD = fLDispl / fTotalCount;
141 G4double MeanLaD2 = fLDispl2 / fTotalCount
142 G4double rmsLaD = std::sqrt(std::fabs(Mean
143
144 G4double MeanPsi = fPsiSpa / (fTotalCount)
145 G4double MeanPsi2 = fPsiSpa2 / (fTotalCoun
146 G4double rmsPsi = std::sqrt(std::fabs(Mean
147
148 G4double MeanTeta = fTetPrj / (2 * fTotalC
149 G4double MeanTeta2 = fTetPrj2 / (2 * fTota
150 G4double rmsTeta = std::sqrt(std::fabs(Mea
151
152 G4double MeanCorrel = fPhiCor / (fTotalCou
153 G4double MeanCorrel2 = fPhiCor2 / (fTotalC
154 G4double rmsCorrel = std::sqrt(std::fabs(M
155
156 G4cout << "\n\n truePathLength :\t" << G4B
157 << G4BestUnit(rmsTPL, "Length") <<
158 << G4BestUnit(MeanGPL, "Length") <<
159 << "\n lateralDisplac :\t" << G4Bes
160 << G4BestUnit(rmsLaD, "Length") <<
161 << " +- " << rmsPsi / mrad << " mra
162 << " (" << MeanPsi / deg << " deg"
163 << " +- " << rmsPsi / deg << " deg)
164
165 G4cout << "\n Theta_plane :\t" << rmsTe
166 << " (" << rmsTeta / deg << " deg)
167 << "\n phi correlation:\t" << MeanC
168 << " (std::cos(phi_pos - phi_dir))
169
170 // cross check from G4EmCalculator
171 //
172 G4cout << "\n Verification from G4EmCalcul
173
174 G4EmCalculator emCal;
175
176 // get transport mean free path (for multi
177 G4double MSmfp = emCal.GetMeanFreePath(ene
178
179 // get range from restricted dedx
180 G4double range = emCal.GetRangeFromRestric
181
182 // effective facRange
183 G4double efFacrange = MeanTPL / std::max(M
184 if (MeanTPL / range >= 0.99) efFacrange =
185
186 G4cout << "\n transport mean free path :\t
187 << "\n range from restrict dE/dx:\t
188 << "\n ---> effective facRange :\t
189
190 G4cout << "\n compute theta0 from Highland
191 << " (" << ComputeMscHighland(Mean
192 }
193 else
194 G4cout << G4endl;
195
196 // restore default format
197 G4cout.precision(prec);
198
199 // delete and remove all contents in fProcCo
200 while (fProcCounter->size() > 0) {
201 OneProcessCount* aProcCount = fProcCounter
202 fProcCounter->pop_back();
203 delete aProcCount;
204 }
205 delete fProcCounter;
206
207 // save histograms
208 G4AnalysisManager* analysisManager = G4Analy
209 if (analysisManager->IsActive()) {
210 analysisManager->Write();
211 analysisManager->CloseFile();
212 }
213
214 // show Rndm status
215 CLHEP::HepRandom::showEngineStatus();
216 }
217
218 //....oooOO0OOooo........oooOO0OOooo........oo
219
220 G4double RunAction::ComputeMscHighland(G4doubl
221 {
222 // compute the width of the Gaussian central
223 // projected angular distribution.
224 // Eur. Phys. Jour. C15 (2000) page 166, for
225
226 G4double t = pathLength / (fDetector->GetMat
227 if (t < DBL_MIN) return 0.;
228
229 G4ParticleGun* particle = fPrimary->GetParti
230 G4double T = particle->GetParticleEnergy();
231 G4double M = particle->GetParticleDefinition
232 G4double z = std::abs(particle->GetParticleD
233
234 G4double bpc = T * (T + 2 * M) / (T + M);
235 G4double teta0 = 13.6 * MeV * z * std::sqrt(
236 return teta0;
237 }
238
239 //....oooOO0OOooo........oooOO0OOooo........oo
240