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