| 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/TestEm0/src/RunActio 26 /// \file electromagnetic/TestEm0/src/RunAction.cc
27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class
28 // 28 //
29 // << 29 // $Id: RunAction.cc 99373 2016-09-20 07:13:41Z 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 "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh"
37 37
38 #include "G4Electron.hh" << 38 #include "G4Run.hh"
>> 39 #include "G4ProcessManager.hh"
>> 40 #include "G4UnitsTable.hh"
39 #include "G4EmCalculator.hh" 41 #include "G4EmCalculator.hh"
40 #include "G4LossTableManager.hh" << 42 #include "G4Electron.hh"
41 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh"
42 #include "G4Positron.hh" <<
43 #include "G4ProcessManager.hh" <<
44 #include "G4Run.hh" <<
45 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh"
46 #include "G4UnitsTable.hh" << 45 #include "G4Electron.hh"
>> 46 #include "G4Positron.hh"
47 47
48 #include <vector> 48 #include <vector>
49 49
50 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 51
52 RunAction::RunAction(DetectorConstruction* det 52 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
53 : fDetector(det), fPrimary(kin) << 53 :G4UserRunAction(),fDetector(det), fPrimary(kin)
54 {} << 54 { }
>> 55
>> 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 57
>> 58 RunAction::~RunAction()
>> 59 { }
>> 60
55 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 62
57 void RunAction::BeginOfRunAction(const G4Run*) 63 void RunAction::BeginOfRunAction(const G4Run*)
58 { 64 {
59 // set precision for printing << 65 //set precision for printing
60 G4int prec = G4cout.precision(6); 66 G4int prec = G4cout.precision(6);
61 << 67
62 // instanciate EmCalculator << 68 //instanciate EmCalculator
63 G4EmCalculator emCal; 69 G4EmCalculator emCal;
64 // emCal.SetVerbose(2); 70 // emCal.SetVerbose(2);
65 << 71
66 // get particle << 72 // get particle
67 G4ParticleDefinition* particle = fPrimary->G << 73 G4ParticleDefinition* particle = fPrimary->GetParticleGun()
>> 74 ->GetParticleDefinition();
68 G4String partName = particle->GetParticleNam 75 G4String partName = particle->GetParticleName();
69 G4double charge = particle->GetPDGCharge(); << 76 G4double charge = particle->GetPDGCharge();
70 G4double energy = fPrimary->GetParticleGun() << 77 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
71 << 78
72 // get material 79 // get material
73 const G4Material* material = fDetector->GetM 80 const G4Material* material = fDetector->GetMaterial();
74 G4String matName = material->GetName(); << 81 G4String matName = material->GetName();
75 G4double density = material->GetDensity(); << 82 G4double density = material->GetDensity();
76 G4double radl = material->GetRadlen(); << 83 G4double radl = material->GetRadlen();
77 << 84
78 G4cout << "\n " << partName << " (" << G4Bes << 85 G4cout << "\n " << partName << " ("
79 << material->GetName() << " (density: << 86 << G4BestUnit(energy,"Energy") << ") in "
80 << "; radiation length: " << G4Best << 87 << material->GetName() << " (density: "
>> 88 << G4BestUnit(density,"Volumic Mass") << "; radiation length: "
>> 89 << G4BestUnit(radl, "Length") << ")" << G4endl;
81 90
82 // get cuts << 91 // get cuts
83 GetCuts(); 92 GetCuts();
84 if (charge != 0.) { 93 if (charge != 0.) {
85 G4cout << "\n Range cuts: \t gamma " << s << 94 G4cout << "\n Range cuts : \t gamma "
86 << "\t e- " << std::setw(12) << G4B << 95 << std::setw(8) << G4BestUnit(fRangeCut[0],"Length")
87 G4cout << "\n Energy cuts: \t gamma " << s << 96 << "\t e- " << std::setw(8) << G4BestUnit(fRangeCut[1],"Length");
88 << "\t e- " << std::setw(12) << G4B << 97 G4cout << "\n Energy cuts : \t gamma "
89 } << 98 << std::setw(8) << G4BestUnit(fEnergyCut[0],"Energy")
90 << 99 << "\t e- " << std::setw(8) << G4BestUnit(fEnergyCut[1],"Energy")
>> 100 << G4endl;
>> 101 }
>> 102
91 // max energy transfert 103 // max energy transfert
92 if (charge != 0.) { 104 if (charge != 0.) {
93 G4double Mass_c2 = particle->GetPDGMass(); << 105 G4double Mass_c2 = particle->GetPDGMass();
94 G4double moverM = electron_mass_c2 / Mass_ << 106 G4double moverM = electron_mass_c2/Mass_c2;
95 G4double gamM1 = energy / Mass_c2, gam = g << 107 G4double gamM1 = energy/Mass_c2, gam = gamM1 + 1., gamP1 = gam + 1.;
96 G4double Tmax = energy; << 108 G4double Tmax = energy;
97 if (particle == G4Electron::Electron()) { << 109 if(particle == G4Electron::Electron()) {
98 Tmax *= 0.5; << 110 Tmax *= 0.5;
99 } << 111 } else if(particle != G4Positron::Positron()) {
100 else if (particle != G4Positron::Positron( << 112 Tmax = (2*electron_mass_c2*gamM1*gamP1)/(1.+2*gam*moverM+moverM*moverM);
101 Tmax = (2 * electron_mass_c2 * gamM1 * g << 113 }
102 } << 114 G4double range = emCal.GetCSDARange(Tmax,G4Electron::Electron(),material);
103 G4double range = emCal.GetCSDARange(Tmax, << 115
104 << 116 G4cout << "\n Max_energy _transferable : " << G4BestUnit(Tmax,"Energy")
105 G4cout << "\n Max_energy _transferable : << 117 << " (" << G4BestUnit(range,"Length") << ")" << G4endl;
106 << G4BestUnit(range, "Length") << " <<
107 } 118 }
108 << 119
109 // get processList and extract EM processes 120 // get processList and extract EM processes (but not MultipleScattering)
110 G4ProcessVector* plist = particle->GetProces 121 G4ProcessVector* plist = particle->GetProcessManager()->GetProcessList();
111 G4String procName; 122 G4String procName;
112 G4double cut; 123 G4double cut;
113 std::vector<G4String> emName; 124 std::vector<G4String> emName;
114 std::vector<G4double> enerCut; 125 std::vector<G4double> enerCut;
115 size_t length = plist->size(); 126 size_t length = plist->size();
116 for (size_t j = 0; j < length; j++) { << 127 for (size_t j=0; j<length; j++) {
117 procName = (*plist)[j]->GetProcessName(); << 128 procName = (*plist)[j]->GetProcessName();
118 cut = fEnergyCut[1]; << 129 cut = fEnergyCut[1];
119 if ((procName == "eBrem") || (procName == << 130 if ((procName == "eBrem")||(procName == "muBrems")) cut = fEnergyCut[0];
120 if (((*plist)[j]->GetProcessType() == fEle << 131 if (((*plist)[j]->GetProcessType() == fElectromagnetic) &&
121 emName.push_back(procName); << 132 (procName != "msc")) {
122 enerCut.push_back(cut); << 133 emName.push_back(procName);
123 } << 134 enerCut.push_back(cut);
124 } << 135 }
125 <<
126 // write html documentation, if requested <<
127 char* htmlDocName = std::getenv("G4PhysListN <<
128 char* htmlDocDir = std::getenv("G4PhysListDo <<
129 if (htmlDocName && htmlDocDir) { <<
130 G4LossTableManager::Instance()->DumpHtml() <<
131 } 136 }
132 << 137
133 // print list of processes 138 // print list of processes
134 G4cout << "\n processes : "; 139 G4cout << "\n processes : ";
135 for (size_t j = 0; j < emName.size(); ++j) { << 140 for (size_t j=0; j<emName.size();j++)
136 G4cout << "\t" << std::setw(14) << emName[ << 141 G4cout << "\t" << std::setw(13) << emName[j] << "\t";
137 } << 142 G4cout << "\t" << std::setw(13) <<"total";
138 G4cout << "\t" << std::setw(14) << "total"; << 143
139 << 144 //compute cross section per atom (only for single material)
140 // compute cross section per atom (only for <<
141 if (material->GetNumberOfElements() == 1) { 145 if (material->GetNumberOfElements() == 1) {
142 G4double Z = material->GetZ(); 146 G4double Z = material->GetZ();
143 G4double A = material->GetA(); 147 G4double A = material->GetA();
144 << 148
145 std::vector<G4double> sigma0; 149 std::vector<G4double> sigma0;
146 G4double sig, sigtot = 0.; 150 G4double sig, sigtot = 0.;
147 151
148 for (size_t j = 0; j < emName.size(); j++) << 152 for (size_t j=0; j<emName.size();j++) {
149 sig = emCal.ComputeCrossSectionPerAtom(e << 153 sig = emCal.ComputeCrossSectionPerAtom
150 sigtot += sig; << 154 (energy,particle,emName[j],Z,A,enerCut[j]);
151 sigma0.push_back(sig); << 155 sigtot += sig;
>> 156 sigma0.push_back(sig);
152 } 157 }
153 sigma0.push_back(sigtot); 158 sigma0.push_back(sigtot);
154 159
155 G4cout << "\n \n cross section per atom << 160 G4cout << "\n \n cross section per atom : ";
156 for (size_t j = 0; j < sigma0.size(); ++j) << 161 for (size_t j=0; j<sigma0.size();j++) {
157 G4cout << "\t" << std::setw(9) << G4Best << 162 G4cout << "\t" << std::setw(13) << G4BestUnit(sigma0[j], "Surface");
158 } 163 }
159 G4cout << G4endl; 164 G4cout << G4endl;
160 } 165 }
161 << 166
162 // get cross section per volume << 167 //get cross section per volume
163 std::vector<G4double> sigma0; 168 std::vector<G4double> sigma0;
164 std::vector<G4double> sigma1; 169 std::vector<G4double> sigma1;
165 std::vector<G4double> sigma2; 170 std::vector<G4double> sigma2;
166 G4double Sig, SigtotComp = 0., Sigtot = 0.; 171 G4double Sig, SigtotComp = 0., Sigtot = 0.;
167 172
168 for (size_t j = 0; j < emName.size(); ++j) { << 173 for (size_t j=0; j<emName.size();j++) {
169 Sig = emCal.ComputeCrossSectionPerVolume(e << 174 Sig = emCal.ComputeCrossSectionPerVolume
170 SigtotComp += Sig; << 175 (energy,particle,emName[j],material,enerCut[j]);
>> 176 SigtotComp += Sig;
171 sigma0.push_back(Sig); 177 sigma0.push_back(Sig);
172 Sig = emCal.GetCrossSectionPerVolume(energ << 178 Sig = emCal.GetCrossSectionPerVolume(energy,particle,emName[j],material);
173 Sigtot += Sig; << 179 Sigtot += Sig;
174 sigma1.push_back(Sig); 180 sigma1.push_back(Sig);
175 sigma2.push_back(Sig / density); << 181 sigma2.push_back(Sig/density);
176 } 182 }
177 sigma0.push_back(SigtotComp); 183 sigma0.push_back(SigtotComp);
178 sigma1.push_back(Sigtot); 184 sigma1.push_back(Sigtot);
179 sigma2.push_back(Sigtot / density); << 185 sigma2.push_back(Sigtot/density);
180 << 186
181 // print cross sections << 187 //print cross sections
182 G4cout << "\n compCrossSectionPerVolume: "; << 188 G4cout << "\n \n compCrossSectionPerVolume : ";
183 for (size_t j = 0; j < sigma0.size(); ++j) { << 189 for (size_t j=0; j<sigma0.size();j++) {
184 G4cout << "\t" << std::setw(9) << sigma0[j << 190 G4cout << "\t" << std::setw(13) << sigma0[j]*cm << " cm^-1";
185 } 191 }
186 G4cout << "\n cross section per volume : "; 192 G4cout << "\n cross section per volume : ";
187 for (size_t j = 0; j < sigma1.size(); ++j) { << 193 for (size_t j=0; j<sigma1.size();j++) {
188 G4cout << "\t" << std::setw(9) << sigma1[j << 194 G4cout << "\t" << std::setw(13) << sigma1[j]*cm << " cm^-1";
189 } 195 }
190 << 196
191 G4cout << "\n cross section per mass : "; 197 G4cout << "\n cross section per mass : ";
192 for (size_t j = 0; j < sigma2.size(); ++j) { << 198 for (size_t j=0; j<sigma2.size();j++) {
193 G4cout << "\t" << std::setw(9) << G4BestUn << 199 G4cout << "\t" << std::setw(13)
194 } << 200 << G4BestUnit(sigma2[j], "Surface/Mass");
195 << 201 }
196 // print mean free path << 202
197 << 203 //print mean free path
>> 204
198 G4double lambda; 205 G4double lambda;
199 << 206
200 G4cout << "\n \n mean free path : 207 G4cout << "\n \n mean free path : ";
201 for (size_t j = 0; j < sigma1.size(); ++j) { << 208 for (size_t j=0; j<sigma1.size();j++) {
202 lambda = DBL_MAX; << 209 lambda = DBL_MAX;
203 if (sigma1[j] > 0.) lambda = 1 / sigma1[j] << 210 if (sigma1[j] > 0.) lambda = 1/sigma1[j];
204 G4cout << "\t" << std::setw(9) << G4BestUn << 211 G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Length");
205 } << 212 }
206 << 213
207 // mean free path (g/cm2) << 214 //mean free path (g/cm2)
208 G4cout << "\n (g/cm2) : "; << 215 G4cout << "\n (g/cm2) : ";
209 for (size_t j = 0; j < sigma2.size(); ++j) { << 216 for (size_t j=0; j<sigma2.size();j++) {
210 lambda = DBL_MAX; << 217 lambda = DBL_MAX;
211 if (sigma2[j] > 0.) lambda = 1 / sigma2[j] << 218 if (sigma2[j] > 0.) lambda = 1/sigma2[j];
212 G4cout << "\t" << std::setw(9) << G4BestUn << 219 G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Mass/Surface");
213 } 220 }
214 G4cout << G4endl; 221 G4cout << G4endl;
215 << 222
216 if (charge == 0.) { 223 if (charge == 0.) {
217 G4cout.precision(prec); 224 G4cout.precision(prec);
218 G4cout << "\n----------------------------- << 225 G4cout << "\n-----------------------------------------------------------\n"
>> 226 << G4endl;
219 return; 227 return;
220 } 228 }
221 << 229
222 // get stopping power << 230 //get stopping power
223 std::vector<G4double> dedx1; 231 std::vector<G4double> dedx1;
224 std::vector<G4double> dedx2; << 232 std::vector<G4double> dedx2;
225 G4double dedx, dedxtot = 0.; 233 G4double dedx, dedxtot = 0.;
226 size_t nproc = emName.size(); 234 size_t nproc = emName.size();
227 235
228 for (size_t j = 0; j < nproc; ++j) { << 236 for (size_t j=0; j<nproc; j++) {
229 dedx = emCal.ComputeDEDX(energy, particle, << 237 dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,enerCut[j]);
230 dedxtot += dedx; 238 dedxtot += dedx;
231 dedx1.push_back(dedx); 239 dedx1.push_back(dedx);
232 dedx2.push_back(dedx / density); << 240 dedx2.push_back(dedx/density);
233 } 241 }
234 dedx1.push_back(dedxtot); 242 dedx1.push_back(dedxtot);
235 dedx2.push_back(dedxtot / density); << 243 dedx2.push_back(dedxtot/density);
236 << 244
237 // print stopping power << 245 //print stopping power
238 G4cout << "\n \n restricted dE/dx : 246 G4cout << "\n \n restricted dE/dx : ";
239 for (size_t j = 0; j <= nproc; ++j) { << 247 for (size_t j=0; j<=nproc; j++) {
240 G4cout << "\t" << std::setw(9) << G4BestUn << 248 G4cout << "\t" << std::setw(13)
>> 249 << G4BestUnit(dedx1[j],"Energy/Length");
241 } 250 }
242 << 251
243 G4cout << "\n (MeV/g/cm2) : "; 252 G4cout << "\n (MeV/g/cm2) : ";
244 for (size_t j = 0; j <= nproc; ++j) { << 253 for (size_t j=0; j<=nproc; j++) {
245 G4cout << "\t" << std::setw(9) << G4BestUn << 254 G4cout << "\t" << std::setw(13)
>> 255 << G4BestUnit(dedx2[j],"Energy*Surface/Mass");
246 } 256 }
247 dedxtot = 0.; 257 dedxtot = 0.;
248 258
249 for (size_t j = 0; j < nproc; ++j) { << 259 for (size_t j=0; j<nproc; j++) {
250 dedx = emCal.ComputeDEDX(energy, particle, << 260 dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,energy);
251 dedxtot += dedx; 261 dedxtot += dedx;
252 dedx1[j] = dedx; 262 dedx1[j] = dedx;
253 dedx2[j] = dedx / density; << 263 dedx2[j] = dedx/density;
254 } 264 }
255 dedx1[nproc] = dedxtot; 265 dedx1[nproc] = dedxtot;
256 dedx2[nproc] = dedxtot / density; << 266 dedx2[nproc] = dedxtot/density;
257 << 267
258 // print stopping power << 268 //print stopping power
259 G4cout << "\n \n unrestricted dE/dx : 269 G4cout << "\n \n unrestricted dE/dx : ";
260 for (size_t j = 0; j <= nproc; ++j) { << 270 for (size_t j=0; j<=nproc; j++) {
261 G4cout << "\t" << std::setw(9) << G4BestUn << 271 G4cout << "\t" << std::setw(13) << G4BestUnit(dedx1[j],"Energy/Length");
262 } 272 }
263 << 273
264 G4cout << "\n (MeV/g/cm2) : "; 274 G4cout << "\n (MeV/g/cm2) : ";
265 for (size_t j = 0; j <= nproc; ++j) { << 275 for (size_t j=0; j<=nproc; j++) {
266 G4cout << "\t" << std::setw(9) << G4BestUn << 276 G4cout << "\t" << std::setw(13)
267 } << 277 << G4BestUnit(dedx2[j],"Energy*Surface/Mass");
268 << 278 }
269 // get range from restricted dedx << 279
270 G4double range1 = emCal.GetRangeFromRestrict << 280 //get range from restricted dedx
271 G4double range2 = range1 * density; << 281 G4double range1 = emCal.GetRangeFromRestricteDEDX(energy,particle,material);
272 << 282 G4double range2 = range1*density;
273 // print range << 283
274 G4cout << "\n \n range from restrict dE/dx: << 284 //print range
275 << "\t" << std::setw(9) << G4BestUnit << 285 G4cout << "\n \n range from restrict dE/dx: "
276 << G4BestUnit(range2, "Mass/Surface") << 286 << "\t" << std::setw(8) << G4BestUnit(range1,"Length")
277 << 287 << " (" << std::setw(8) << G4BestUnit(range2,"Mass/Surface") << ")";
278 // get range from full dedx << 288
279 G4double EmaxTable = G4EmParameters::Instanc << 289 //get range from full dedx
280 if (energy < EmaxTable) { << 290 if(energy < GeV) {
281 G4double Range1 = emCal.GetCSDARange(energ << 291 G4double Range1 = emCal.GetCSDARange(energy,particle,material);
282 G4double Range2 = Range1 * density; << 292 G4double Range2 = Range1*density;
283 << 293
284 G4cout << "\n range from full dE/dx : << 294 G4cout << "\n range from full dE/dx : "
285 << "\t" << std::setw(9) << G4BestUn << 295 << "\t" << std::setw(8) << G4BestUnit(Range1,"Length")
286 << G4BestUnit(Range2, "Mass/Surface << 296 << " (" << std::setw(8) << G4BestUnit(Range2,"Mass/Surface") << ")";
287 } << 297 }
288 << 298
289 // get transport mean free path (for multipl << 299 //get transport mean free path (for multiple scattering)
290 G4double MSmfp1 = emCal.GetMeanFreePath(ener << 300 G4double MSmfp1 = emCal.GetMeanFreePath(energy,particle,"msc",material);
291 G4double MSmfp2 = MSmfp1 * density; << 301 G4double MSmfp2 = MSmfp1*density;
292 << 302
293 // print transport mean free path << 303 //print transport mean free path
294 G4cout << "\n \n transport mean free path : << 304 G4cout << "\n \n transport mean free path : "
295 << "\t" << std::setw(9) << G4BestUnit << 305 << "\t" << std::setw(8) << G4BestUnit(MSmfp1,"Length")
296 << G4BestUnit(MSmfp2, "Mass/Surface") << 306 << " (" << std::setw(8) << G4BestUnit(MSmfp2,"Mass/Surface") << ")";
297 307
298 if (particle == G4Electron::Electron()) Crit 308 if (particle == G4Electron::Electron()) CriticalEnergy();
299 << 309
300 G4cout << "\n------------------------------- 310 G4cout << "\n-------------------------------------------------------------\n";
301 G4cout << G4endl; 311 G4cout << G4endl;
302 << 312
303 // reset default precision << 313 // reset default precision
304 G4cout.precision(prec); << 314 G4cout.precision(prec);
305 } 315 }
306 316
307 //....oooOO0OOooo........oooOO0OOooo........oo 317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 318
309 void RunAction::EndOfRunAction(const G4Run*) { << 319 void RunAction::EndOfRunAction(const G4Run* )
>> 320 { }
310 321
311 //....oooOO0OOooo........oooOO0OOooo........oo 322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312 323
313 #include "G4ProductionCutsTable.hh" 324 #include "G4ProductionCutsTable.hh"
314 325
315 //....oooOO0OOooo........oooOO0OOooo........oo 326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
316 327
317 void RunAction::GetCuts() 328 void RunAction::GetCuts()
318 { << 329 {
319 G4ProductionCutsTable* theCoupleTable = G4Pr << 330 G4ProductionCutsTable* theCoupleTable =
320 << 331 G4ProductionCutsTable::GetProductionCutsTable();
>> 332
321 size_t numOfCouples = theCoupleTable->GetTab 333 size_t numOfCouples = theCoupleTable->GetTableSize();
322 const G4MaterialCutsCouple* couple = 0; 334 const G4MaterialCutsCouple* couple = 0;
323 G4int index = 0; 335 G4int index = 0;
324 for (size_t i = 0; i < numOfCouples; i++) { << 336 for (size_t i=0; i<numOfCouples; i++) {
325 couple = theCoupleTable->GetMaterialCutsCo << 337 couple = theCoupleTable->GetMaterialCutsCouple(i);
326 if (couple->GetMaterial() == fDetector->Ge << 338 if (couple->GetMaterial() == fDetector->GetMaterial()) {index = i; break;}
327 index = i; << 339 }
328 break; << 340
329 } << 341 fRangeCut[0] =
330 } << 342 (*(theCoupleTable->GetRangeCutsVector(idxG4GammaCut)))[index];
>> 343 fRangeCut[1] =
>> 344 (*(theCoupleTable->GetRangeCutsVector(idxG4ElectronCut)))[index];
>> 345 fRangeCut[2] =
>> 346 (*(theCoupleTable->GetRangeCutsVector(idxG4PositronCut)))[index];
>> 347
>> 348 fEnergyCut[0] =
>> 349 (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
>> 350 fEnergyCut[1] =
>> 351 (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
>> 352 fEnergyCut[2] =
>> 353 (*(theCoupleTable->GetEnergyCutsVector(idxG4PositronCut)))[index];
331 354
332 fRangeCut[0] = (*(theCoupleTable->GetRangeCu <<
333 fRangeCut[1] = (*(theCoupleTable->GetRangeCu <<
334 fRangeCut[2] = (*(theCoupleTable->GetRangeCu <<
335 <<
336 fEnergyCut[0] = (*(theCoupleTable->GetEnergy <<
337 fEnergyCut[1] = (*(theCoupleTable->GetEnergy <<
338 fEnergyCut[2] = (*(theCoupleTable->GetEnergy <<
339 } 355 }
340 356
341 //....oooOO0OOooo........oooOO0OOooo........oo 357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
342 358
343 void RunAction::CriticalEnergy() 359 void RunAction::CriticalEnergy()
344 { 360 {
345 // compute e- critical energy (Rossi definit 361 // compute e- critical energy (Rossi definition) and Moliere radius.
346 // Review of Particle Physics - Eur. Phys. J 362 // Review of Particle Physics - Eur. Phys. J. C3 (1998) page 147
347 // 363 //
348 G4EmCalculator emCal; 364 G4EmCalculator emCal;
349 << 365
350 const G4Material* material = fDetector->GetM 366 const G4Material* material = fDetector->GetMaterial();
351 const G4double radl = material->GetRadlen(); 367 const G4double radl = material->GetRadlen();
352 G4double ekin = 5 * MeV; << 368 G4double ekin = 5*MeV;
353 G4double deioni; 369 G4double deioni;
354 G4double err = 1., errmax = 0.001; << 370 G4double err = 1., errmax = 0.001;
355 G4int iter = 0, itermax = 10; << 371 G4int iter = 0 , itermax = 10;
356 while (err > errmax && iter < itermax) { 372 while (err > errmax && iter < itermax) {
357 iter++; << 373 iter++;
358 deioni = radl * emCal.ComputeDEDX(ekin, G4 << 374 deioni = radl*
359 err = std::abs(deioni - ekin) / ekin; << 375 emCal.ComputeDEDX(ekin,G4Electron::Electron(),"eIoni",material);
>> 376 err = std::abs(deioni - ekin)/ekin;
360 ekin = deioni; 377 ekin = deioni;
361 } 378 }
362 G4cout << "\n \n critical energy (Rossi) : << 379 G4cout << "\n \n critical energy (Rossi) : "
363 << "\t" << std::setw(8) << G4BestUnit << 380 << "\t" << std::setw(8) << G4BestUnit(ekin,"Energy");
364 << 381
365 // Pdg formula (only for single material) << 382 //Pdg formula (only for single material)
366 G4double pdga[2] = {610 * MeV, 710 * MeV}; << 383 G4double pdga[2] = { 610*MeV, 710*MeV };
367 G4double pdgb[2] = {1.24, 0.92}; << 384 G4double pdgb[2] = { 1.24, 0.92 };
368 G4double EcPdg = 0.; 385 G4double EcPdg = 0.;
369 << 386
370 if (material->GetNumberOfElements() == 1) { 387 if (material->GetNumberOfElements() == 1) {
371 G4int istat = 0; 388 G4int istat = 0;
372 if (material->GetState() == kStateGas) ist << 389 if (material->GetState() == kStateGas) istat = 1;
373 G4double Zeff = material->GetZ() + pdgb[is 390 G4double Zeff = material->GetZ() + pdgb[istat];
374 EcPdg = pdga[istat] / Zeff; << 391 EcPdg = pdga[istat]/Zeff;
375 G4cout << "\t\t\t (from Pdg formula : " << << 392 G4cout << "\t\t\t (from Pdg formula : "
376 } << 393 << std::setw(8) << G4BestUnit(EcPdg,"Energy") << ")";
377 << 394 }
378 const G4double Es = 21.2052 * MeV; << 395
379 G4double rMolier1 = Es / ekin, rMolier2 = rM << 396 const G4double Es = 21.2052*MeV;
380 G4cout << "\n Moliere radius : " << 397 G4double rMolier1 = Es/ekin, rMolier2 = rMolier1*radl;
381 << "\t" << std::setw(8) << rMolier1 < << 398 G4cout << "\n Moliere radius : "
382 << "= " << std::setw(8) << G4BestUnit << 399 << "\t" << std::setw(8) << rMolier1 << " X0 "
383 << 400 << "= " << std::setw(8) << G4BestUnit(rMolier2,"Length");
384 if (material->GetNumberOfElements() == 1) { << 401
385 G4double rMPdg = radl * Es / EcPdg; << 402 if (material->GetNumberOfElements() == 1) {
386 G4cout << "\t (from Pdg formula : " << std << 403 G4double rMPdg = radl*Es/EcPdg;
387 } << 404 G4cout << "\t (from Pdg formula : "
>> 405 << std::setw(8) << G4BestUnit(rMPdg,"Length") << ")";
>> 406 }
388 } 407 }
389 408
390 //....oooOO0OOooo........oooOO0OOooo........oo 409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
391 410