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 ]

Diff markup

Differences between /examples/extended/electromagnetic/TestEm0/src/RunAction.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm0/src/RunAction.cc (Version 6.2.p2)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 /// \file electromagnetic/TestEm0/src/RunActio    
 27 /// \brief Implementation of the RunAction cla    
 28 //                                                
 29 //                                                
 30 //....oooOO0OOooo........oooOO0OOooo........oo    
 31 //....oooOO0OOooo........oooOO0OOooo........oo    
 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........oo    
 51                                                   
 52 RunAction::RunAction(DetectorConstruction* det    
 53   : fDetector(det), fPrimary(kin)                 
 54 {}                                                
 55 //....oooOO0OOooo........oooOO0OOooo........oo    
 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->G    
 68   G4String partName = particle->GetParticleNam    
 69   G4double charge = particle->GetPDGCharge();     
 70   G4double energy = fPrimary->GetParticleGun()    
 71                                                   
 72   // get material                                 
 73   const G4Material* material = fDetector->GetM    
 74   G4String matName = material->GetName();         
 75   G4double density = material->GetDensity();      
 76   G4double radl = material->GetRadlen();          
 77                                                   
 78   G4cout << "\n " << partName << " (" << G4Bes    
 79          << material->GetName() << " (density:    
 80          << ";   radiation length: " << G4Best    
 81                                                   
 82   // get cuts                                     
 83   GetCuts();                                      
 84   if (charge != 0.) {                             
 85     G4cout << "\n  Range cuts: \t gamma " << s    
 86            << "\t e- " << std::setw(12) << G4B    
 87     G4cout << "\n Energy cuts: \t gamma " << s    
 88            << "\t e- " << std::setw(12) << G4B    
 89   }                                               
 90                                                   
 91   // max energy transfert                         
 92   if (charge != 0.) {                             
 93     G4double Mass_c2 = particle->GetPDGMass();    
 94     G4double moverM = electron_mass_c2 / Mass_    
 95     G4double gamM1 = energy / Mass_c2, gam = g    
 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 * g    
102     }                                             
103     G4double range = emCal.GetCSDARange(Tmax,     
104                                                   
105     G4cout << "\n  Max_energy _transferable  :    
106            << G4BestUnit(range, "Length") << "    
107   }                                               
108                                                   
109   // get processList and extract EM processes     
110   G4ProcessVector* plist = particle->GetProces    
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 ==     
120     if (((*plist)[j]->GetProcessType() == fEle    
121       emName.push_back(procName);                 
122       enerCut.push_back(cut);                     
123     }                                             
124   }                                               
125                                                   
126   // write html documentation, if requested       
127   char* htmlDocName = std::getenv("G4PhysListN    
128   char* htmlDocDir = std::getenv("G4PhysListDo    
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[    
137   }                                               
138   G4cout << "\t" << std::setw(14) << "total";     
139                                                   
140   // compute cross section per atom (only for     
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(e    
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) << G4Best    
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(e    
170     SigtotComp += Sig;                            
171     sigma0.push_back(Sig);                        
172     Sig = emCal.GetCrossSectionPerVolume(energ    
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    
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    
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) << G4BestUn    
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) << G4BestUn    
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) << G4BestUn    
213   }                                               
214   G4cout << G4endl;                               
215                                                   
216   if (charge == 0.) {                             
217     G4cout.precision(prec);                       
218     G4cout << "\n-----------------------------    
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,    
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) << G4BestUn    
241   }                                               
242                                                   
243   G4cout << "\n      (MeV/g/cm2)          : ";    
244   for (size_t j = 0; j <= nproc; ++j) {           
245     G4cout << "\t" << std::setw(9) << G4BestUn    
246   }                                               
247   dedxtot = 0.;                                   
248                                                   
249   for (size_t j = 0; j < nproc; ++j) {            
250     dedx = emCal.ComputeDEDX(energy, particle,    
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) << G4BestUn    
262   }                                               
263                                                   
264   G4cout << "\n      (MeV/g/cm2)          : ";    
265   for (size_t j = 0; j <= nproc; ++j) {           
266     G4cout << "\t" << std::setw(9) << G4BestUn    
267   }                                               
268                                                   
269   // get range from restricted dedx               
270   G4double range1 = emCal.GetRangeFromRestrict    
271   G4double range2 = range1 * density;             
272                                                   
273   // print range                                  
274   G4cout << "\n \n  range from restrict dE/dx:    
275          << "\t" << std::setw(9) << G4BestUnit    
276          << G4BestUnit(range2, "Mass/Surface")    
277                                                   
278   // get range from full dedx                     
279   G4double EmaxTable = G4EmParameters::Instanc    
280   if (energy < EmaxTable) {                       
281     G4double Range1 = emCal.GetCSDARange(energ    
282     G4double Range2 = Range1 * density;           
283                                                   
284     G4cout << "\n  range from full dE/dx    :     
285            << "\t" << std::setw(9) << G4BestUn    
286            << G4BestUnit(Range2, "Mass/Surface    
287   }                                               
288                                                   
289   // get transport mean free path (for multipl    
290   G4double MSmfp1 = emCal.GetMeanFreePath(ener    
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    
296          << G4BestUnit(MSmfp2, "Mass/Surface")    
297                                                   
298   if (particle == G4Electron::Electron()) Crit    
299                                                   
300   G4cout << "\n-------------------------------    
301   G4cout << G4endl;                               
302                                                   
303   // reset default precision                      
304   G4cout.precision(prec);                         
305 }                                                 
306                                                   
307 //....oooOO0OOooo........oooOO0OOooo........oo    
308                                                   
309 void RunAction::EndOfRunAction(const G4Run*) {    
310                                                   
311 //....oooOO0OOooo........oooOO0OOooo........oo    
312                                                   
313 #include "G4ProductionCutsTable.hh"               
314                                                   
315 //....oooOO0OOooo........oooOO0OOooo........oo    
316                                                   
317 void RunAction::GetCuts()                         
318 {                                                 
319   G4ProductionCutsTable* theCoupleTable = G4Pr    
320                                                   
321   size_t numOfCouples = theCoupleTable->GetTab    
322   const G4MaterialCutsCouple* couple = 0;         
323   G4int index = 0;                                
324   for (size_t i = 0; i < numOfCouples; i++) {     
325     couple = theCoupleTable->GetMaterialCutsCo    
326     if (couple->GetMaterial() == fDetector->Ge    
327       index = i;                                  
328       break;                                      
329     }                                             
330   }                                               
331                                                   
332   fRangeCut[0] = (*(theCoupleTable->GetRangeCu    
333   fRangeCut[1] = (*(theCoupleTable->GetRangeCu    
334   fRangeCut[2] = (*(theCoupleTable->GetRangeCu    
335                                                   
336   fEnergyCut[0] = (*(theCoupleTable->GetEnergy    
337   fEnergyCut[1] = (*(theCoupleTable->GetEnergy    
338   fEnergyCut[2] = (*(theCoupleTable->GetEnergy    
339 }                                                 
340                                                   
341 //....oooOO0OOooo........oooOO0OOooo........oo    
342                                                   
343 void RunAction::CriticalEnergy()                  
344 {                                                 
345   // compute e- critical energy (Rossi definit    
346   // Review of Particle Physics - Eur. Phys. J    
347   //                                              
348   G4EmCalculator emCal;                           
349                                                   
350   const G4Material* material = fDetector->GetM    
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, G4    
359     err = std::abs(deioni - ekin) / ekin;         
360     ekin = deioni;                                
361   }                                               
362   G4cout << "\n \n  critical energy (Rossi)  :    
363          << "\t" << std::setw(8) << G4BestUnit    
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) ist    
373     G4double Zeff = material->GetZ() + pdgb[is    
374     EcPdg = pdga[istat] / Zeff;                   
375     G4cout << "\t\t\t (from Pdg formula : " <<    
376   }                                               
377                                                   
378   const G4double Es = 21.2052 * MeV;              
379   G4double rMolier1 = Es / ekin, rMolier2 = rM    
380   G4cout << "\n  Moliere radius           : "     
381          << "\t" << std::setw(8) << rMolier1 <    
382          << "= " << std::setw(8) << G4BestUnit    
383                                                   
384   if (material->GetNumberOfElements() == 1) {     
385     G4double rMPdg = radl * Es / EcPdg;           
386     G4cout << "\t (from Pdg formula : " << std    
387   }                                               
388 }                                                 
389                                                   
390 //....oooOO0OOooo........oooOO0OOooo........oo    
391