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