| 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/TestEm15/src/RunActi 26 /// \file electromagnetic/TestEm15/src/RunAction.cc
27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class
28 // 28 //
29 // << 29 // $Id: RunAction.cc 82283 2014-06-13 14:49:40Z 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
35 #include "DetectorConstruction.hh" 36 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh" <<
37 #include "PrimaryGeneratorAction.hh" 37 #include "PrimaryGeneratorAction.hh"
>> 38 #include "HistoManager.hh"
38 39
39 #include "G4EmCalculator.hh" <<
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"
44 #include "Randomize.hh" << 43 #include "G4EmCalculator.hh"
45 44
>> 45 #include "Randomize.hh"
>> 46 #include "G4SystemOfUnits.hh"
46 #include <iomanip> 47 #include <iomanip>
47 48
48 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 50
50 RunAction::RunAction(DetectorConstruction* det 51 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
51 : G4UserRunAction(), fDetector(det), fPrimar << 52 : G4UserRunAction(),fDetector(det), fPrimary(prim), fProcCounter(0),
>> 53 fHistoManager(0)
52 { 54 {
53 fHistoManager = new HistoManager(); << 55 fHistoManager = new HistoManager();
54 } 56 }
55 57
56 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 59
58 RunAction::~RunAction() 60 RunAction::~RunAction()
59 { 61 {
60 delete fHistoManager; << 62 delete fHistoManager;
61 } 63 }
62 64
63 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 66
65 void RunAction::BeginOfRunAction(const G4Run*) 67 void RunAction::BeginOfRunAction(const G4Run*)
66 { << 68 {
67 // save Rndm status 69 // save Rndm status
68 ////G4RunManager::GetRunManager()->SetRandom 70 ////G4RunManager::GetRunManager()->SetRandomNumberStore(true);
69 CLHEP::HepRandom::showEngineStatus(); 71 CLHEP::HepRandom::showEngineStatus();
70 72
71 fProcCounter = new ProcessesCount; 73 fProcCounter = new ProcessesCount;
72 fTotalCount = 0; 74 fTotalCount = 0;
73 << 75
74 fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0. 76 fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0.;
75 fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0. 77 fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0.;
76 fTetPrj = fTetPrj2 = 0.; 78 fTetPrj = fTetPrj2 = 0.;
77 fPhiCor = fPhiCor2 = 0.; 79 fPhiCor = fPhiCor2 = 0.;
78 << 80
79 // histograms << 81 //histograms
80 // 82 //
81 G4AnalysisManager* analysisManager = G4Analy 83 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
82 if (analysisManager->IsActive()) { << 84 if ( analysisManager->IsActive() ) {
83 analysisManager->OpenFile(); 85 analysisManager->OpenFile();
84 } << 86 }
85 } 87 }
86 88
87 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 90
89 void RunAction::CountProcesses(G4String procNa 91 void RunAction::CountProcesses(G4String procName)
90 { 92 {
91 // does the process already encounted ? << 93 //does the process already encounted ?
92 size_t nbProc = fProcCounter->size(); << 94 size_t nbProc = fProcCounter->size();
93 size_t i = 0; << 95 size_t i = 0;
94 while ((i < nbProc) && ((*fProcCounter)[i]-> << 96 while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
95 i++; << 97 if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName));
96 if (i == nbProc) fProcCounter->push_back(new <<
97 98
98 (*fProcCounter)[i]->Count(); << 99 (*fProcCounter)[i]->Count();
99 } 100 }
100 101
101 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 103
103 void RunAction::EndOfRunAction(const G4Run* aR 104 void RunAction::EndOfRunAction(const G4Run* aRun)
104 { 105 {
105 G4int NbOfEvents = aRun->GetNumberOfEvent(); 106 G4int NbOfEvents = aRun->GetNumberOfEvent();
106 if (NbOfEvents == 0) return; 107 if (NbOfEvents == 0) return;
107 << 108
108 G4int prec = G4cout.precision(5); << 109 G4int prec = G4cout.precision(5);
109 << 110
110 G4Material* material = fDetector->GetMateria 111 G4Material* material = fDetector->GetMaterial();
111 G4double density = material->GetDensity(); 112 G4double density = material->GetDensity();
112 << 113
113 G4ParticleDefinition* particle = fPrimary->G << 114 G4ParticleDefinition* particle =
114 G4String Particle = particle->GetParticleNam << 115 fPrimary->GetParticleGun()->GetParticleDefinition();
>> 116 G4String Particle = particle->GetParticleName();
115 G4double energy = fPrimary->GetParticleGun() 117 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
116 G4cout << "\n The run consists of " << NbOfE << 118 G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of "
117 << G4BestUnit(energy, "Energy") << " << 119 << G4BestUnit(energy,"Energy") << " through "
118 << G4BestUnit(fDetector->GetBoxSize() << 120 << G4BestUnit(fDetector->GetBoxSize(),"Length") << " of "
119 << " (density: " << G4BestUnit(densit << 121 << material->GetName() << " (density: "
120 << 122 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
121 // frequency of processes << 123
>> 124 //frequency of processes
122 G4cout << "\n Process calls frequency --->"; 125 G4cout << "\n Process calls frequency --->";
123 for (size_t i = 0; i < fProcCounter->size(); << 126 for (size_t i=0; i< fProcCounter->size();i++) {
124 G4String procName = (*fProcCounter)[i]->Ge << 127 G4String procName = (*fProcCounter)[i]->GetName();
125 G4int count = (*fProcCounter)[i]->GetCount << 128 G4int count = (*fProcCounter)[i]->GetCounter();
126 G4cout << "\t" << procName << " = " << cou << 129 G4cout << "\t" << procName << " = " << count;
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 } 130 }
193 else << 131
194 G4cout << G4endl; << 132 if (fTotalCount == 0) return;
195 << 133
196 // restore default format << 134 //compute path length and related quantities
197 G4cout.precision(prec); << 135 //
>> 136 G4double MeanTPL = fTruePL /fTotalCount;
>> 137 G4double MeanTPL2 = fTruePL2/fTotalCount;
>> 138 G4double rmsTPL = std::sqrt(std::fabs(MeanTPL2 - MeanTPL*MeanTPL));
>> 139
>> 140 G4double MeanGPL = fGeomPL /fTotalCount;
>> 141 G4double MeanGPL2 = fGeomPL2/fTotalCount;
>> 142 G4double rmsGPL = std::sqrt(std::fabs(MeanGPL2 - MeanGPL*MeanGPL));
>> 143
>> 144 G4double MeanLaD = fLDispl /fTotalCount;
>> 145 G4double MeanLaD2 = fLDispl2/fTotalCount;
>> 146 G4double rmsLaD = std::sqrt(std::fabs(MeanLaD2 - MeanLaD*MeanLaD));
>> 147
>> 148 G4double MeanPsi = fPsiSpa /(fTotalCount);
>> 149 G4double MeanPsi2 = fPsiSpa2/(fTotalCount);
>> 150 G4double rmsPsi = std::sqrt(std::fabs(MeanPsi2 - MeanPsi*MeanPsi));
>> 151
>> 152 G4double MeanTeta = fTetPrj /(2*fTotalCount);
>> 153 G4double MeanTeta2 = fTetPrj2/(2*fTotalCount);
>> 154 G4double rmsTeta = std::sqrt(std::fabs(MeanTeta2 - MeanTeta*MeanTeta));
>> 155
>> 156 G4double MeanCorrel = fPhiCor /(fTotalCount);
>> 157 G4double MeanCorrel2 = fPhiCor2/(fTotalCount);
>> 158 G4double rmsCorrel = std::sqrt(std::fabs(MeanCorrel2-MeanCorrel*MeanCorrel));
>> 159
>> 160 G4cout << "\n\n truePathLength :\t" << G4BestUnit(MeanTPL,"Length")
>> 161 << " +- " << G4BestUnit( rmsTPL,"Length")
>> 162 << "\n geomPathLength :\t" << G4BestUnit(MeanGPL,"Length")
>> 163 << " +- " << G4BestUnit( rmsGPL,"Length")
>> 164 << "\n lateralDisplac :\t" << G4BestUnit(MeanLaD,"Length")
>> 165 << " +- " << G4BestUnit( rmsLaD,"Length")
>> 166 << "\n Psi :\t" << MeanPsi/mrad << " mrad"
>> 167 << " +- " << rmsPsi /mrad << " mrad"
>> 168 << " (" << MeanPsi/deg << " deg"
>> 169 << " +- " << rmsPsi /deg << " deg)"
>> 170 << G4endl;
>> 171
>> 172 G4cout << "\n Theta_plane :\t" << rmsTeta/mrad << " mrad"
>> 173 << " (" << rmsTeta/deg << " deg)"
>> 174 << "\n phi correlation:\t" << MeanCorrel
>> 175 << " +- " << rmsCorrel
>> 176 << " (std::cos(phi_pos - phi_dir))"
>> 177 << G4endl;
>> 178
198 179
199 // delete and remove all contents in fProcCo << 180 //cross check from G4EmCalculator
200 while (fProcCounter->size() > 0) { << 181 //
201 OneProcessCount* aProcCount = fProcCounter << 182 G4cout << "\n Verification from G4EmCalculator. \n";
>> 183
>> 184 G4EmCalculator emCal;
>> 185
>> 186 //get transport mean free path (for multiple scattering)
>> 187 G4double MSmfp = emCal.GetMeanFreePath(energy,particle,"msc",material);
>> 188
>> 189 //get range from restricted dedx
>> 190 G4double range = emCal.GetRangeFromRestricteDEDX(energy,particle,material);
>> 191
>> 192 //effective facRange
>> 193 G4double efFacrange = MeanTPL/std::max(MSmfp, range);
>> 194 if (MeanTPL/range >= 0.99) efFacrange = 1.;
>> 195
>> 196 G4cout << "\n transport mean free path :\t" << G4BestUnit(MSmfp,"Length")
>> 197 << "\n range from restrict dE/dx:\t" << G4BestUnit(range,"Length")
>> 198 << "\n ---> effective facRange :\t" << efFacrange
>> 199 << G4endl;
>> 200
>> 201 G4cout << "\n compute theta0 from Highland :\t"
>> 202 << ComputeMscHighland(MeanTPL)/mrad << " mrad"
>> 203 << " (" << ComputeMscHighland(MeanTPL)/deg << " deg)"
>> 204 << G4endl;
>> 205
>> 206 //restore default format
>> 207 G4cout.precision(prec);
>> 208
>> 209 // delete and remove all contents in fProcCounter
>> 210 while (fProcCounter->size()>0){
>> 211 OneProcessCount* aProcCount=fProcCounter->back();
202 fProcCounter->pop_back(); 212 fProcCounter->pop_back();
203 delete aProcCount; 213 delete aProcCount;
204 } 214 }
205 delete fProcCounter; 215 delete fProcCounter;
206 << 216
207 // save histograms << 217 //save histograms
208 G4AnalysisManager* analysisManager = G4Analy << 218 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
209 if (analysisManager->IsActive()) { << 219 if ( analysisManager->IsActive() ) {
210 analysisManager->Write(); 220 analysisManager->Write();
211 analysisManager->CloseFile(); 221 analysisManager->CloseFile();
212 } << 222 }
213 223
214 // show Rndm status 224 // show Rndm status
215 CLHEP::HepRandom::showEngineStatus(); 225 CLHEP::HepRandom::showEngineStatus();
216 } 226 }
217 227
218 //....oooOO0OOooo........oooOO0OOooo........oo 228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219 229
220 G4double RunAction::ComputeMscHighland(G4doubl 230 G4double RunAction::ComputeMscHighland(G4double pathLength)
221 { 231 {
222 // compute the width of the Gaussian central << 232 //compute the width of the Gaussian central part of the MultipleScattering
223 // projected angular distribution. << 233 //projected angular distribution.
224 // Eur. Phys. Jour. C15 (2000) page 166, for << 234 //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
225 << 235
226 G4double t = pathLength / (fDetector->GetMat << 236 G4double t = pathLength/(fDetector->GetMaterial()->GetRadlen());
227 if (t < DBL_MIN) return 0.; << 237 if (t < DBL_MIN) return 0.;
228 << 238
229 G4ParticleGun* particle = fPrimary->GetParti << 239 G4ParticleGun* particle = fPrimary->GetParticleGun();
230 G4double T = particle->GetParticleEnergy(); << 240 G4double T = particle->GetParticleEnergy();
231 G4double M = particle->GetParticleDefinition << 241 G4double M = particle->GetParticleDefinition()->GetPDGMass();
232 G4double z = std::abs(particle->GetParticleD << 242 G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge()/eplus);
233 << 243
234 G4double bpc = T * (T + 2 * M) / (T + M); << 244 G4double bpc = T*(T+2*M)/(T+M);
235 G4double teta0 = 13.6 * MeV * z * std::sqrt( << 245 G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
236 return teta0; << 246 return teta0;
237 } 247 }
238 248
239 //....oooOO0OOooo........oooOO0OOooo........oo 249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240 250