source: trunk/examples/extended/electromagnetic/TestEm0/src/RunAction.cc@ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 13.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: RunAction.cc,v 1.15 2010/05/10 13:45:49 maire Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-01 $
28//
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31
32#include "RunAction.hh"
33#include "DetectorConstruction.hh"
34#include "PrimaryGeneratorAction.hh"
35
36#include "G4Run.hh"
37#include "G4ProcessManager.hh"
38#include "G4UnitsTable.hh"
39#include "G4EmCalculator.hh"
40#include "G4Electron.hh"
41
42#include <vector>
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45
46RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
47:detector(det), primary(kin)
48{ }
49
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52RunAction::~RunAction()
53{ }
54
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57void RunAction::BeginOfRunAction(const G4Run*)
58{
59 //set precision for printing
60 G4int prec = G4cout.precision(6);
61
62 //instanciate EmCalculator
63 G4EmCalculator emCal;
64 // emCal.SetVerbose(2);
65
66 // get particle
67 G4ParticleDefinition* particle = primary->GetParticleGun()
68 ->GetParticleDefinition();
69 G4String partName = particle->GetParticleName();
70 G4double charge = particle->GetPDGCharge();
71 G4double energy = primary->GetParticleGun()->GetParticleEnergy();
72
73 // get material
74 G4Material* material = detector->GetMaterial();
75 G4String matName = material->GetName();
76 G4double density = material->GetDensity();
77 G4double radl = material->GetRadlen();
78
79 G4cout << "\n " << partName << " ("
80 << G4BestUnit(energy,"Energy") << ") in "
81 << material->GetName() << " (density: "
82 << G4BestUnit(density,"Volumic Mass") << "; radiation length: "
83 << G4BestUnit(radl, "Length") << ")" << G4endl;
84
85 // get cuts
86 GetCuts();
87 if (charge != 0.) {
88 G4cout << "\n Range cuts : \t gamma "
89 << std::setw(8) << G4BestUnit(rangeCut[0],"Length")
90 << "\t e- " << std::setw(8) << G4BestUnit(rangeCut[1],"Length");
91 G4cout << "\n Energy cuts : \t gamma "
92 << std::setw(8) << G4BestUnit(energyCut[0],"Energy")
93 << "\t e- " << std::setw(8) << G4BestUnit(energyCut[1],"Energy")
94 << G4endl;
95 }
96
97 // max energy transfert
98 if (charge != 0.) {
99 G4double Mass_c2 = particle->GetPDGMass();
100 G4double moverM = electron_mass_c2/Mass_c2;
101 G4double gamM1 = energy/Mass_c2, gam = gamM1 + 1., gamP1 = gam + 1.;
102 G4double Tmax =
103 (2*electron_mass_c2*gamM1*gamP1)/(1.+2*gam*moverM+moverM*moverM);
104 G4double range = emCal.GetCSDARange(Tmax,G4Electron::Electron(),material);
105
106 G4cout << "\n Max_energy _transferable : " << G4BestUnit(Tmax,"Energy")
107 << " (" << G4BestUnit(range,"Length") << ")" << G4endl;
108 }
109
110 // get processList and extract EM processes (but not MultipleScattering)
111 G4ProcessVector* plist = particle->GetProcessManager()->GetProcessList();
112 G4String procName;
113 G4double cut;
114 std::vector<G4String> emName;
115 std::vector<G4double> enerCut;
116 size_t length = plist->size();
117 for (size_t j=0; j<length; j++) {
118 procName = (*plist)[j]->GetProcessName();
119 cut = energyCut[1];
120 if ((procName == "eBrem")||(procName == "muBrems")) cut = energyCut[0];
121 if (((*plist)[j]->GetProcessType() == fElectromagnetic) &&
122 (procName != "msc")) {
123 emName.push_back(procName);
124 enerCut.push_back(cut);
125 }
126 }
127
128 // print list of processes
129 G4cout << "\n processes : ";
130 for (size_t j=0; j<emName.size();j++)
131 G4cout << "\t" << std::setw(13) << emName[j] << "\t";
132 G4cout << "\t" << std::setw(13) <<"total";
133
134 //compute cross section per atom (only for single material)
135 if (material->GetNumberOfElements() == 1) {
136 G4double Z = material->GetZ();
137 G4double A = material->GetA();
138
139 std::vector<G4double> sigma0;
140 G4double sig, sigtot = 0.;
141
142 for (size_t j=0; j<emName.size();j++) {
143 sig = emCal.ComputeCrossSectionPerAtom
144 (energy,particle,emName[j],Z,A,enerCut[j]);
145 sigtot += sig;
146 sigma0.push_back(sig);
147 }
148 sigma0.push_back(sigtot);
149
150 G4cout << "\n \n cross section per atom : ";
151 for (size_t j=0; j<sigma0.size();j++) {
152 G4cout << "\t" << std::setw(13) << G4BestUnit(sigma0[j], "Surface");
153 }
154 G4cout << G4endl;
155 }
156
157 //get cross section per volume
158 std::vector<G4double> sigma0;
159 std::vector<G4double> sigma1;
160 std::vector<G4double> sigma2;
161 G4double Sig, SigtotComp = 0., Sigtot = 0.;
162
163 for (size_t j=0; j<emName.size();j++) {
164 Sig = emCal.ComputeCrossSectionPerVolume
165 (energy,particle,emName[j],material,enerCut[j]);
166 SigtotComp += Sig;
167 sigma0.push_back(Sig);
168 Sig = emCal.GetCrossSectionPerVolume(energy,particle,emName[j],material);
169 Sigtot += Sig;
170 sigma1.push_back(Sig);
171 sigma2.push_back(Sig/density);
172 }
173 sigma0.push_back(SigtotComp);
174 sigma1.push_back(Sigtot);
175 sigma2.push_back(Sigtot/density);
176
177 //print cross sections
178 G4cout << "\n \n compCrossSectionPerVolume : ";
179 for (size_t j=0; j<sigma0.size();j++) {
180 G4cout << "\t" << std::setw(13) << sigma0[j]*cm << " cm^-1";
181 }
182 G4cout << "\n cross section per volume : ";
183 for (size_t j=0; j<sigma1.size();j++) {
184 G4cout << "\t" << std::setw(13) << sigma1[j]*cm << " cm^-1";
185 }
186
187 G4cout << "\n cross section per mass : ";
188 for (size_t j=0; j<sigma2.size();j++) {
189 G4cout << "\t" << std::setw(13)
190 << G4BestUnit(sigma2[j], "Surface/Mass");
191 }
192
193 //print mean free path
194
195 G4double lambda;
196
197 G4cout << "\n \n mean free path : ";
198 for (size_t j=0; j<sigma1.size();j++) {
199 lambda = DBL_MAX;
200 if (sigma1[j] > 0.) lambda = 1/sigma1[j];
201 G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Length");
202 }
203
204 //mean free path (g/cm2)
205 G4cout << "\n (g/cm2) : ";
206 for (size_t j=0; j<sigma2.size();j++) {
207 lambda = DBL_MAX;
208 if (sigma2[j] > 0.) lambda = 1/sigma2[j];
209 G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Mass/Surface");
210 }
211 G4cout << G4endl;
212
213 if (charge == 0.) {
214 G4cout.precision(prec);
215 G4cout << "\n-------------------------------------------------------------\n"
216 << G4endl;
217 return;
218 }
219
220 //get stopping power
221 std::vector<G4double> dedx1;
222 std::vector<G4double> dedx2;
223 G4double dedx, dedxtot = 0.;
224
225 for (size_t j=0; j<emName.size();j++) {
226 dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,enerCut[j]);
227 dedx1.push_back(dedx);
228 dedx2.push_back(dedx/density);
229 }
230 dedxtot = emCal.GetDEDX(energy,particle,material);
231 dedx1.push_back(dedxtot);
232 dedx2.push_back(dedxtot/density);
233
234 //print stopping power
235 G4cout << "\n \n restricted dE/dx : ";
236 for (size_t j=0; j<sigma1.size();j++) {
237 G4cout << "\t" << std::setw(13) << G4BestUnit(dedx1[j],"Energy/Length");
238 }
239
240 G4cout << "\n (MeV/g/cm2) : ";
241 for (size_t j=0; j<sigma2.size();j++) {
242 G4cout << "\t" << std::setw(13) << G4BestUnit(dedx2[j],"Energy*Surface/Mass");
243 }
244
245 //get range from restricted dedx
246 G4double range1 = emCal.GetRangeFromRestricteDEDX(energy,particle,material);
247 G4double range2 = range1*density;
248
249 //get range from full dedx
250 G4double Range1 = emCal.GetCSDARange(energy,particle,material);
251 G4double Range2 = Range1*density;
252
253 //print range
254 G4cout << "\n \n range from restrict dE/dx: "
255 << "\t" << std::setw(8) << G4BestUnit(range1,"Length")
256 << " (" << std::setw(8) << G4BestUnit(range2,"Mass/Surface") << ")";
257
258 G4cout << "\n range from full dE/dx : "
259 << "\t" << std::setw(8) << G4BestUnit(Range1,"Length")
260 << " (" << std::setw(8) << G4BestUnit(Range2,"Mass/Surface") << ")";
261
262 //get transport mean free path (for multiple scattering)
263 G4double MSmfp1 = emCal.GetMeanFreePath(energy,particle,"msc",material);
264 G4double MSmfp2 = MSmfp1*density;
265
266 //print transport mean free path
267 G4cout << "\n \n transport mean free path : "
268 << "\t" << std::setw(8) << G4BestUnit(MSmfp1,"Length")
269 << " (" << std::setw(8) << G4BestUnit(MSmfp2,"Mass/Surface") << ")";
270
271 if (particle == G4Electron::Electron()) CriticalEnergy();
272
273 G4cout << "\n-------------------------------------------------------------\n";
274 G4cout << G4endl;
275
276 // reset default precision
277 G4cout.precision(prec);
278}
279
280//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
281
282void RunAction::EndOfRunAction(const G4Run* )
283{ }
284
285//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
286
287#include "G4ProductionCutsTable.hh"
288
289//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
290
291void RunAction::GetCuts()
292{
293 G4ProductionCutsTable* theCoupleTable =
294 G4ProductionCutsTable::GetProductionCutsTable();
295
296 size_t numOfCouples = theCoupleTable->GetTableSize();
297 const G4MaterialCutsCouple* couple = 0;
298 G4int index = 0;
299 for (size_t i=0; i<numOfCouples; i++) {
300 couple = theCoupleTable->GetMaterialCutsCouple(i);
301 if (couple->GetMaterial() == detector->GetMaterial()) {index = i; break;}
302 }
303
304 rangeCut[0] =
305 (*(theCoupleTable->GetRangeCutsVector(idxG4GammaCut)))[index];
306 rangeCut[1] =
307 (*(theCoupleTable->GetRangeCutsVector(idxG4ElectronCut)))[index];
308 rangeCut[2] =
309 (*(theCoupleTable->GetRangeCutsVector(idxG4PositronCut)))[index];
310
311 energyCut[0] =
312 (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
313 energyCut[1] =
314 (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
315 energyCut[2] =
316 (*(theCoupleTable->GetEnergyCutsVector(idxG4PositronCut)))[index];
317
318}
319
320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
321
322void RunAction::CriticalEnergy()
323{
324 // compute e- critical energy (Rossi definition) and Moliere radius.
325 // Review of Particle Physics - Eur. Phys. J. C3 (1998) page 147
326 //
327 G4EmCalculator emCal;
328
329 const G4Material* material = detector->GetMaterial();
330 const G4double radl = material->GetRadlen();
331 G4double ekin = 5*MeV;
332 G4double deioni;
333 G4double err = 1., errmax = 0.001;
334 G4int iter = 0 , itermax = 10;
335 while (err > errmax && iter < itermax) {
336 iter++;
337 deioni = radl*
338 emCal.ComputeDEDX(ekin,G4Electron::Electron(),"eIoni",material);
339 err = std::abs(deioni - ekin)/ekin;
340 ekin = deioni;
341 }
342 G4cout << "\n \n critical energy (Rossi) : "
343 << "\t" << std::setw(8) << G4BestUnit(ekin,"Energy");
344
345 //Pdg formula (only for single material)
346 G4double pdga[2] = { 610*MeV, 710*MeV };
347 G4double pdgb[2] = { 1.24, 0.92 };
348 G4double EcPdg = 0.;
349
350 if (material->GetNumberOfElements() == 1) {
351 G4int istat = 0;
352 if (material->GetState() == kStateGas) istat = 1;
353 G4double Zeff = material->GetZ() + pdgb[istat];
354 EcPdg = pdga[istat]/Zeff;
355 G4cout << "\t\t\t (from Pdg formula : "
356 << std::setw(8) << G4BestUnit(EcPdg,"Energy") << ")";
357 }
358
359 const G4double Es = 21.2052*MeV;
360 G4double rMolier1 = Es/ekin, rMolier2 = rMolier1*radl;
361 G4cout << "\n Moliere radius : "
362 << "\t" << std::setw(8) << rMolier1 << " X0 "
363 << "= " << std::setw(8) << G4BestUnit(rMolier2,"Length");
364
365 if (material->GetNumberOfElements() == 1) {
366 G4double rMPdg = radl*Es/EcPdg;
367 G4cout << "\t (from Pdg formula : "
368 << std::setw(8) << G4BestUnit(rMPdg,"Length") << ")";
369 }
370}
371
372//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.