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

Last change on this file was 1337, checked in by garnier, 14 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.