Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm0/src/RunAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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