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 10.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 /// \file electromagnetic/TestEm0/src/RunActio     26 /// \file electromagnetic/TestEm0/src/RunAction.cc
 27 /// \brief Implementation of the RunAction cla     27 /// \brief Implementation of the RunAction class
 28 //                                                 28 //
 29 //                                             <<  29 // $Id: RunAction.cc 68243 2013-03-19 18:03:11Z vnivanch $
                                                   >>  30 // 
 30 //....oooOO0OOooo........oooOO0OOooo........oo     31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oo     32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32                                                    33 
 33 #include "RunAction.hh"                            34 #include "RunAction.hh"
 34                                                << 
 35 #include "DetectorConstruction.hh"                 35 #include "DetectorConstruction.hh"
 36 #include "PrimaryGeneratorAction.hh"               36 #include "PrimaryGeneratorAction.hh"
 37                                                    37 
 38 #include "G4Electron.hh"                       <<  38 #include "G4Run.hh"
                                                   >>  39 #include "G4ProcessManager.hh"
                                                   >>  40 #include "G4UnitsTable.hh"
 39 #include "G4EmCalculator.hh"                       41 #include "G4EmCalculator.hh"
 40 #include "G4LossTableManager.hh"               <<  42 #include "G4Electron.hh"
 41 #include "G4PhysicalConstants.hh"                  43 #include "G4PhysicalConstants.hh"
 42 #include "G4Positron.hh"                       << 
 43 #include "G4ProcessManager.hh"                 << 
 44 #include "G4Run.hh"                            << 
 45 #include "G4SystemOfUnits.hh"                      44 #include "G4SystemOfUnits.hh"
 46 #include "G4UnitsTable.hh"                     << 
 47                                                    45 
 48 #include <vector>                                  46 #include <vector>
 49                                                    47 
 50 //....oooOO0OOooo........oooOO0OOooo........oo     48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 51                                                    49 
 52 RunAction::RunAction(DetectorConstruction* det     50 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
 53   : fDetector(det), fPrimary(kin)              <<  51 :G4UserRunAction(),fDetector(det), fPrimary(kin)
 54 {}                                             <<  52 { }
                                                   >>  53 
                                                   >>  54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >>  55 
                                                   >>  56 RunAction::~RunAction()
                                                   >>  57 { }
                                                   >>  58 
 55 //....oooOO0OOooo........oooOO0OOooo........oo     59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 56                                                    60 
 57 void RunAction::BeginOfRunAction(const G4Run*)     61 void RunAction::BeginOfRunAction(const G4Run*)
 58 {                                                  62 {
 59   // set precision for printing                <<  63   //set precision for printing
 60   G4int prec = G4cout.precision(6);                64   G4int prec = G4cout.precision(6);
 61                                                <<  65   
 62   // instanciate EmCalculator                  <<  66   //instanciate EmCalculator
 63   G4EmCalculator emCal;                            67   G4EmCalculator emCal;
 64   //  emCal.SetVerbose(2);                         68   //  emCal.SetVerbose(2);
 65                                                <<  69      
 66   // get particle                              <<  70   // get particle 
 67   G4ParticleDefinition* particle = fPrimary->G <<  71   G4ParticleDefinition* particle = fPrimary->GetParticleGun()
                                                   >>  72                                           ->GetParticleDefinition();
 68   G4String partName = particle->GetParticleNam     73   G4String partName = particle->GetParticleName();
 69   G4double charge = particle->GetPDGCharge();  <<  74   G4double charge   = particle->GetPDGCharge();    
 70   G4double energy = fPrimary->GetParticleGun() <<  75   G4double energy   = fPrimary->GetParticleGun()->GetParticleEnergy();
 71                                                <<  76  
 72   // get material                                  77   // get material
 73   const G4Material* material = fDetector->GetM <<  78   G4Material* material = fDetector->GetMaterial();
 74   G4String matName = material->GetName();      <<  79   G4String matName     = material->GetName();
 75   G4double density = material->GetDensity();   <<  80   G4double density     = material->GetDensity();
 76   G4double radl = material->GetRadlen();       <<  81   G4double radl        = material->GetRadlen();  
 77                                                <<  82 
 78   G4cout << "\n " << partName << " (" << G4Bes <<  83   G4cout << "\n " << partName << " ("
 79          << material->GetName() << " (density: <<  84          << G4BestUnit(energy,"Energy") << ") in " 
 80          << ";   radiation length: " << G4Best <<  85          << material->GetName() << " (density: " 
                                                   >>  86          << G4BestUnit(density,"Volumic Mass") << ";   radiation length: "
                                                   >>  87          << G4BestUnit(radl,   "Length")       << ")" << G4endl;
 81                                                    88 
 82   // get cuts                                  <<  89   // get cuts         
 83   GetCuts();                                       90   GetCuts();
 84   if (charge != 0.) {                              91   if (charge != 0.) {
 85     G4cout << "\n  Range cuts: \t gamma " << s <<  92    G4cout << "\n  Range cuts : \t gamma "  
 86            << "\t e- " << std::setw(12) << G4B <<  93                       << std::setw(8) << G4BestUnit(fRangeCut[0],"Length")
 87     G4cout << "\n Energy cuts: \t gamma " << s <<  94           << "\t e- " << std::setw(8) << G4BestUnit(fRangeCut[1],"Length");
 88            << "\t e- " << std::setw(12) << G4B <<  95    G4cout << "\n Energy cuts : \t gamma " 
 89   }                                            <<  96                       << std::setw(8) << G4BestUnit(fEnergyCut[0],"Energy")
 90                                                <<  97           << "\t e- " << std::setw(8) << G4BestUnit(fEnergyCut[1],"Energy")
                                                   >>  98           << G4endl;
                                                   >>  99    }
                                                   >> 100    
 91   // max energy transfert                         101   // max energy transfert
 92   if (charge != 0.) {                             102   if (charge != 0.) {
 93     G4double Mass_c2 = particle->GetPDGMass(); << 103   G4double Mass_c2 = particle->GetPDGMass();
 94     G4double moverM = electron_mass_c2 / Mass_ << 104   G4double moverM = electron_mass_c2/Mass_c2;
 95     G4double gamM1 = energy / Mass_c2, gam = g << 105   G4double gamM1 = energy/Mass_c2, gam = gamM1 + 1., gamP1 = gam + 1.;
 96     G4double Tmax = energy;                    << 106   G4double Tmax = 
 97     if (particle == G4Electron::Electron()) {  << 107             (2*electron_mass_c2*gamM1*gamP1)/(1.+2*gam*moverM+moverM*moverM);
 98       Tmax *= 0.5;                             << 108   G4double range = emCal.GetCSDARange(Tmax,G4Electron::Electron(),material);
 99     }                                          << 109   
100     else if (particle != G4Positron::Positron( << 110   G4cout << "\n  Max_energy _transferable  : " << G4BestUnit(Tmax,"Energy")
101       Tmax = (2 * electron_mass_c2 * gamM1 * g << 111          << " (" << G4BestUnit(range,"Length") << ")" << G4endl;               
102     }                                          << 
103     G4double range = emCal.GetCSDARange(Tmax,  << 
104                                                << 
105     G4cout << "\n  Max_energy _transferable  : << 
106            << G4BestUnit(range, "Length") << " << 
107   }                                               112   }
108                                                << 113                
109   // get processList and extract EM processes     114   // get processList and extract EM processes (but not MultipleScattering)
110   G4ProcessVector* plist = particle->GetProces    115   G4ProcessVector* plist = particle->GetProcessManager()->GetProcessList();
111   G4String procName;                              116   G4String procName;
112   G4double cut;                                   117   G4double cut;
113   std::vector<G4String> emName;                   118   std::vector<G4String> emName;
114   std::vector<G4double> enerCut;                  119   std::vector<G4double> enerCut;
115   size_t length = plist->size();                  120   size_t length = plist->size();
116   for (size_t j = 0; j < length; j++) {        << 121   for (size_t j=0; j<length; j++) {
117     procName = (*plist)[j]->GetProcessName();  << 122      procName = (*plist)[j]->GetProcessName();
118     cut = fEnergyCut[1];                       << 123      cut = fEnergyCut[1];
119     if ((procName == "eBrem") || (procName ==  << 124      if ((procName == "eBrem")||(procName == "muBrems")) cut = fEnergyCut[0];
120     if (((*plist)[j]->GetProcessType() == fEle << 125      if (((*plist)[j]->GetProcessType() == fElectromagnetic) &&
121       emName.push_back(procName);              << 126          (procName != "msc")) {
122       enerCut.push_back(cut);                  << 127        emName.push_back(procName);
123     }                                          << 128        enerCut.push_back(cut);
                                                   >> 129      }  
124   }                                               130   }
125                                                << 131   
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                      132   // print list of processes
134   G4cout << "\n  processes :                ";    133   G4cout << "\n  processes :                ";
135   for (size_t j = 0; j < emName.size(); ++j) { << 134   for (size_t j=0; j<emName.size();j++)
136     G4cout << "\t" << std::setw(14) << emName[ << 135     G4cout << "\t" << std::setw(13) << emName[j] << "\t";
137   }                                            << 136   G4cout << "\t" << std::setw(13) <<"total";
138   G4cout << "\t" << std::setw(14) << "total";  << 137   
139                                                << 138   //compute cross section per atom (only for single material)
140   // compute cross section per atom (only for  << 
141   if (material->GetNumberOfElements() == 1) {     139   if (material->GetNumberOfElements() == 1) {
142     G4double Z = material->GetZ();                140     G4double Z = material->GetZ();
143     G4double A = material->GetA();                141     G4double A = material->GetA();
144                                                << 142      
145     std::vector<G4double> sigma0;                 143     std::vector<G4double> sigma0;
146     G4double sig, sigtot = 0.;                    144     G4double sig, sigtot = 0.;
147                                                   145 
148     for (size_t j = 0; j < emName.size(); j++) << 146     for (size_t j=0; j<emName.size();j++) {
149       sig = emCal.ComputeCrossSectionPerAtom(e << 147       sig = emCal.ComputeCrossSectionPerAtom
150       sigtot += sig;                           << 148                       (energy,particle,emName[j],Z,A,enerCut[j]);
151       sigma0.push_back(sig);                   << 149       sigtot += sig;                              
                                                   >> 150       sigma0.push_back(sig);                
152     }                                             151     }
153     sigma0.push_back(sigtot);                     152     sigma0.push_back(sigtot);
154                                                   153 
155     G4cout << "\n \n  cross section per atom   << 154     G4cout << "\n \n  cross section per atom    : ";
156     for (size_t j = 0; j < sigma0.size(); ++j) << 155     for (size_t j=0; j<sigma0.size();j++) {             
157       G4cout << "\t" << std::setw(9) << G4Best << 156       G4cout << "\t" << std::setw(13) << G4BestUnit(sigma0[j], "Surface");
158     }                                             157     }
159     G4cout << G4endl;                             158     G4cout << G4endl;
160   }                                               159   }
161                                                << 160     
162   // get cross section per volume              << 161   //get cross section per volume 
163   std::vector<G4double> sigma0;                   162   std::vector<G4double> sigma0;
164   std::vector<G4double> sigma1;                   163   std::vector<G4double> sigma1;
165   std::vector<G4double> sigma2;                   164   std::vector<G4double> sigma2;
166   G4double Sig, SigtotComp = 0., Sigtot = 0.;     165   G4double Sig, SigtotComp = 0., Sigtot = 0.;
167                                                   166 
168   for (size_t j = 0; j < emName.size(); ++j) { << 167   for (size_t j=0; j<emName.size();j++) {
169     Sig = emCal.ComputeCrossSectionPerVolume(e << 168     Sig = emCal.ComputeCrossSectionPerVolume
170     SigtotComp += Sig;                         << 169       (energy,particle,emName[j],material,enerCut[j]);  
                                                   >> 170     SigtotComp += Sig;    
171     sigma0.push_back(Sig);                        171     sigma0.push_back(Sig);
172     Sig = emCal.GetCrossSectionPerVolume(energ << 172     Sig = emCal.GetCrossSectionPerVolume(energy,particle,emName[j],material);      
173     Sigtot += Sig;                             << 173     Sigtot += Sig;    
174     sigma1.push_back(Sig);                        174     sigma1.push_back(Sig);
175     sigma2.push_back(Sig / density);           << 175     sigma2.push_back(Sig/density);                        
176   }                                               176   }
177   sigma0.push_back(SigtotComp);                   177   sigma0.push_back(SigtotComp);
178   sigma1.push_back(Sigtot);                       178   sigma1.push_back(Sigtot);
179   sigma2.push_back(Sigtot / density);          << 179   sigma2.push_back(Sigtot/density);          
180                                                << 180     
181   // print cross sections                      << 181   //print cross sections
182   G4cout << "\n  compCrossSectionPerVolume: "; << 182   G4cout << "\n \n  compCrossSectionPerVolume : ";
183   for (size_t j = 0; j < sigma0.size(); ++j) { << 183   for (size_t j=0; j<sigma0.size();j++) {             
184     G4cout << "\t" << std::setw(9) << sigma0[j << 184     G4cout << "\t" << std::setw(13) << sigma0[j]*cm << " cm^-1";
185   }                                               185   }
186   G4cout << "\n  cross section per volume : ";    186   G4cout << "\n  cross section per volume : ";
187   for (size_t j = 0; j < sigma1.size(); ++j) { << 187   for (size_t j=0; j<sigma1.size();j++) {             
188     G4cout << "\t" << std::setw(9) << sigma1[j << 188     G4cout << "\t" << std::setw(13) << sigma1[j]*cm << " cm^-1";
189   }                                               189   }
190                                                << 190   
191   G4cout << "\n  cross section per mass   : ";    191   G4cout << "\n  cross section per mass   : ";
192   for (size_t j = 0; j < sigma2.size(); ++j) { << 192   for (size_t j=0; j<sigma2.size();j++) {
193     G4cout << "\t" << std::setw(9) << G4BestUn << 193     G4cout << "\t" << std::setw(13) 
194   }                                            << 194            << G4BestUnit(sigma2[j], "Surface/Mass");
195                                                << 195   }
196   // print mean free path                      << 196    
197                                                << 197   //print mean free path
                                                   >> 198   
198   G4double lambda;                                199   G4double lambda;
199                                                << 200   
200   G4cout << "\n \n  mean free path           :    201   G4cout << "\n \n  mean free path           : ";
201   for (size_t j = 0; j < sigma1.size(); ++j) { << 202   for (size_t j=0; j<sigma1.size();j++) {
202     lambda = DBL_MAX;                          << 203     lambda = DBL_MAX; 
203     if (sigma1[j] > 0.) lambda = 1 / sigma1[j] << 204     if (sigma1[j] > 0.) lambda = 1/sigma1[j];
204     G4cout << "\t" << std::setw(9) << G4BestUn << 205     G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Length");
205   }                                            << 206   }
206                                                << 207   
207   // mean free path (g/cm2)                    << 208   //mean free path (g/cm2)
208   G4cout << "\n        (g/cm2)            : "; << 209   G4cout << "\n        (g/cm2)            : ";  
209   for (size_t j = 0; j < sigma2.size(); ++j) { << 210   for (size_t j=0; j<sigma2.size();j++) {
210     lambda = DBL_MAX;                          << 211     lambda =  DBL_MAX;
211     if (sigma2[j] > 0.) lambda = 1 / sigma2[j] << 212     if (sigma2[j] > 0.) lambda = 1/sigma2[j];                       
212     G4cout << "\t" << std::setw(9) << G4BestUn << 213     G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Mass/Surface");    
213   }                                               214   }
214   G4cout << G4endl;                               215   G4cout << G4endl;
215                                                << 216   
216   if (charge == 0.) {                             217   if (charge == 0.) {
217     G4cout.precision(prec);                       218     G4cout.precision(prec);
218     G4cout << "\n----------------------------- << 219     G4cout << "\n-------------------------------------------------------------\n"
                                                   >> 220            << G4endl;
219     return;                                       221     return;
220   }                                               222   }
221                                                << 223   
222   // get stopping power                        << 224   //get stopping power 
223   std::vector<G4double> dedx1;                    225   std::vector<G4double> dedx1;
224   std::vector<G4double> dedx2;                 << 226   std::vector<G4double> dedx2;  
225   G4double dedx, dedxtot = 0.;                    227   G4double dedx, dedxtot = 0.;
226   size_t nproc = emName.size();                   228   size_t nproc = emName.size();
227                                                   229 
228   for (size_t j = 0; j < nproc; ++j) {         << 230   for (size_t j=0; j<nproc; j++) {
229     dedx = emCal.ComputeDEDX(energy, particle, << 231     dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,enerCut[j]);
230     dedxtot += dedx;                              232     dedxtot += dedx;
231     dedx1.push_back(dedx);                        233     dedx1.push_back(dedx);
232     dedx2.push_back(dedx / density);           << 234     dedx2.push_back(dedx/density);                        
233   }                                               235   }
234   dedx1.push_back(dedxtot);                       236   dedx1.push_back(dedxtot);
235   dedx2.push_back(dedxtot / density);          << 237   dedx2.push_back(dedxtot/density);          
236                                                << 238     
237   // print stopping power                      << 239   //print stopping power
238   G4cout << "\n \n  restricted dE/dx         :    240   G4cout << "\n \n  restricted dE/dx         : ";
239   for (size_t j = 0; j <= nproc; ++j) {        << 241   for (size_t j=0; j<=nproc; j++) {             
240     G4cout << "\t" << std::setw(9) << G4BestUn << 242     G4cout << "\t" << std::setw(13) 
                                                   >> 243      << G4BestUnit(dedx1[j],"Energy/Length");
241   }                                               244   }
242                                                << 245   
243   G4cout << "\n      (MeV/g/cm2)          : ";    246   G4cout << "\n      (MeV/g/cm2)          : ";
244   for (size_t j = 0; j <= nproc; ++j) {        << 247   for (size_t j=0; j<=nproc; j++) {
245     G4cout << "\t" << std::setw(9) << G4BestUn << 248     G4cout << "\t" << std::setw(13) 
                                                   >> 249      << G4BestUnit(dedx2[j],"Energy*Surface/Mass");
246   }                                               250   }
247   dedxtot = 0.;                                   251   dedxtot = 0.;
248                                                   252 
249   for (size_t j = 0; j < nproc; ++j) {         << 253   for (size_t j=0; j<nproc; j++) {
250     dedx = emCal.ComputeDEDX(energy, particle, << 254     dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,energy);
251     dedxtot += dedx;                              255     dedxtot += dedx;
252     dedx1[j] = dedx;                              256     dedx1[j] = dedx;
253     dedx2[j] = dedx / density;                 << 257     dedx2[j] = dedx/density;                        
254   }                                               258   }
255   dedx1[nproc] = dedxtot;                         259   dedx1[nproc] = dedxtot;
256   dedx2[nproc] = dedxtot / density;            << 260   dedx2[nproc] = dedxtot/density;          
257                                                << 261     
258   // print stopping power                      << 262   //print stopping power
259   G4cout << "\n \n  unrestricted dE/dx       :    263   G4cout << "\n \n  unrestricted dE/dx       : ";
260   for (size_t j = 0; j <= nproc; ++j) {        << 264   for (size_t j=0; j<=nproc; j++) {             
261     G4cout << "\t" << std::setw(9) << G4BestUn << 265     G4cout << "\t" << std::setw(13) << G4BestUnit(dedx1[j],"Energy/Length");
262   }                                               266   }
263                                                << 267   
264   G4cout << "\n      (MeV/g/cm2)          : ";    268   G4cout << "\n      (MeV/g/cm2)          : ";
265   for (size_t j = 0; j <= nproc; ++j) {        << 269   for (size_t j=0; j<=nproc; j++) {
266     G4cout << "\t" << std::setw(9) << G4BestUn << 270     G4cout << "\t" << std::setw(13) 
267   }                                            << 271      << G4BestUnit(dedx2[j],"Energy*Surface/Mass");
268                                                << 272   }
269   // get range from restricted dedx            << 273   
270   G4double range1 = emCal.GetRangeFromRestrict << 274   //get range from restricted dedx
271   G4double range2 = range1 * density;          << 275   G4double range1 = emCal.GetRangeFromRestricteDEDX(energy,particle,material);
272                                                << 276   G4double range2 = range1*density;
273   // print range                               << 277 
274   G4cout << "\n \n  range from restrict dE/dx: << 278   //print range
275          << "\t" << std::setw(9) << G4BestUnit << 279   G4cout << "\n \n  range from restrict dE/dx: " 
276          << G4BestUnit(range2, "Mass/Surface") << 280          << "\t" << std::setw(8) << G4BestUnit(range1,"Length")
277                                                << 281          << " (" << std::setw(8) << G4BestUnit(range2,"Mass/Surface") << ")";
278   // get range from full dedx                  << 282   
279   G4double EmaxTable = G4EmParameters::Instanc << 283    //get range from full dedx
280   if (energy < EmaxTable) {                    << 284   if(energy < GeV) {
281     G4double Range1 = emCal.GetCSDARange(energ << 285     G4double Range1 = emCal.GetCSDARange(energy,particle,material);
282     G4double Range2 = Range1 * density;        << 286     G4double Range2 = Range1*density;
283                                                << 287      
284     G4cout << "\n  range from full dE/dx    :  << 288     G4cout << "\n  range from full dE/dx    : " 
285            << "\t" << std::setw(9) << G4BestUn << 289      << "\t" << std::setw(8) << G4BestUnit(Range1,"Length")
286            << G4BestUnit(Range2, "Mass/Surface << 290      << " (" << std::setw(8) << G4BestUnit(Range2,"Mass/Surface") << ")";
287   }                                            << 291   }         
288                                                << 292 
289   // get transport mean free path (for multipl << 293   //get transport mean free path (for multiple scattering)
290   G4double MSmfp1 = emCal.GetMeanFreePath(ener << 294   G4double MSmfp1 = emCal.GetMeanFreePath(energy,particle,"msc",material);
291   G4double MSmfp2 = MSmfp1 * density;          << 295   G4double MSmfp2 = MSmfp1*density;
292                                                << 296   
293   // print transport mean free path            << 297   //print transport mean free path
294   G4cout << "\n \n  transport mean free path : << 298   G4cout << "\n \n  transport mean free path : " 
295          << "\t" << std::setw(9) << G4BestUnit << 299          << "\t" << std::setw(8) << G4BestUnit(MSmfp1,"Length")
296          << G4BestUnit(MSmfp2, "Mass/Surface") << 300          << " (" << std::setw(8) << G4BestUnit(MSmfp2,"Mass/Surface") << ")";
297                                                   301 
298   if (particle == G4Electron::Electron()) Crit    302   if (particle == G4Electron::Electron()) CriticalEnergy();
299                                                << 303            
300   G4cout << "\n-------------------------------    304   G4cout << "\n-------------------------------------------------------------\n";
301   G4cout << G4endl;                               305   G4cout << G4endl;
302                                                << 306        
303   // reset default precision                   << 307  // reset default precision
304   G4cout.precision(prec);                      << 308  G4cout.precision(prec);    
305 }                                                 309 }
306                                                   310 
307 //....oooOO0OOooo........oooOO0OOooo........oo    311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308                                                   312 
309 void RunAction::EndOfRunAction(const G4Run*) { << 313 void RunAction::EndOfRunAction(const G4Run* )
                                                   >> 314 { }
310                                                   315 
311 //....oooOO0OOooo........oooOO0OOooo........oo    316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312                                                   317 
313 #include "G4ProductionCutsTable.hh"               318 #include "G4ProductionCutsTable.hh"
314                                                   319 
315 //....oooOO0OOooo........oooOO0OOooo........oo    320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
316                                                   321 
317 void RunAction::GetCuts()                         322 void RunAction::GetCuts()
318 {                                              << 323 {  
319   G4ProductionCutsTable* theCoupleTable = G4Pr << 324   G4ProductionCutsTable* theCoupleTable =
320                                                << 325         G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 326         
321   size_t numOfCouples = theCoupleTable->GetTab    327   size_t numOfCouples = theCoupleTable->GetTableSize();
322   const G4MaterialCutsCouple* couple = 0;         328   const G4MaterialCutsCouple* couple = 0;
323   G4int index = 0;                                329   G4int index = 0;
324   for (size_t i = 0; i < numOfCouples; i++) {  << 330   for (size_t i=0; i<numOfCouples; i++) {
325     couple = theCoupleTable->GetMaterialCutsCo << 331      couple = theCoupleTable->GetMaterialCutsCouple(i);
326     if (couple->GetMaterial() == fDetector->Ge << 332      if (couple->GetMaterial() == fDetector->GetMaterial()) {index = i; break;}
327       index = i;                               << 333   }
328       break;                                   << 334   
329     }                                          << 335   fRangeCut[0] =
330   }                                            << 336          (*(theCoupleTable->GetRangeCutsVector(idxG4GammaCut)))[index];
                                                   >> 337   fRangeCut[1] =      
                                                   >> 338          (*(theCoupleTable->GetRangeCutsVector(idxG4ElectronCut)))[index];
                                                   >> 339   fRangeCut[2] =      
                                                   >> 340          (*(theCoupleTable->GetRangeCutsVector(idxG4PositronCut)))[index]; 
                                                   >> 341 
                                                   >> 342   fEnergyCut[0] =
                                                   >> 343          (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
                                                   >> 344   fEnergyCut[1] =      
                                                   >> 345          (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
                                                   >> 346   fEnergyCut[2] =      
                                                   >> 347          (*(theCoupleTable->GetEnergyCutsVector(idxG4PositronCut)))[index];
331                                                   348 
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 }                                                 349 }
340                                                   350 
341 //....oooOO0OOooo........oooOO0OOooo........oo    351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
342                                                   352 
343 void RunAction::CriticalEnergy()                  353 void RunAction::CriticalEnergy()
344 {                                                 354 {
345   // compute e- critical energy (Rossi definit    355   // compute e- critical energy (Rossi definition) and Moliere radius.
346   // Review of Particle Physics - Eur. Phys. J    356   // Review of Particle Physics - Eur. Phys. J. C3 (1998) page 147
347   //                                              357   //
348   G4EmCalculator emCal;                           358   G4EmCalculator emCal;
349                                                << 359     
350   const G4Material* material = fDetector->GetM    360   const G4Material* material = fDetector->GetMaterial();
351   const G4double radl = material->GetRadlen();    361   const G4double radl = material->GetRadlen();
352   G4double ekin = 5 * MeV;                     << 362   G4double ekin = 5*MeV;
353   G4double deioni;                                363   G4double deioni;
354   G4double err = 1., errmax = 0.001;           << 364   G4double err  = 1., errmax = 0.001;
355   G4int iter = 0, itermax = 10;                << 365   G4int    iter = 0 , itermax = 10;  
356   while (err > errmax && iter < itermax) {        366   while (err > errmax && iter < itermax) {
357     iter++;                                    << 367     iter++;          
358     deioni = radl * emCal.ComputeDEDX(ekin, G4 << 368     deioni  = radl*
359     err = std::abs(deioni - ekin) / ekin;      << 369               emCal.ComputeDEDX(ekin,G4Electron::Electron(),"eIoni",material);
                                                   >> 370     err = std::abs(deioni - ekin)/ekin;
360     ekin = deioni;                                371     ekin = deioni;
361   }                                               372   }
362   G4cout << "\n \n  critical energy (Rossi)  : << 373   G4cout << "\n \n  critical energy (Rossi)  : " 
363          << "\t" << std::setw(8) << G4BestUnit << 374          << "\t" << std::setw(8) << G4BestUnit(ekin,"Energy");
364                                                << 375          
365   // Pdg formula (only for single material)    << 376   //Pdg formula (only for single material)
366   G4double pdga[2] = {610 * MeV, 710 * MeV};   << 377   G4double pdga[2] = { 610*MeV, 710*MeV };
367   G4double pdgb[2] = {1.24, 0.92};             << 378   G4double pdgb[2] = { 1.24, 0.92 };
368   G4double EcPdg = 0.;                            379   G4double EcPdg = 0.;
369                                                << 380   
370   if (material->GetNumberOfElements() == 1) {     381   if (material->GetNumberOfElements() == 1) {
371     G4int istat = 0;                              382     G4int istat = 0;
372     if (material->GetState() == kStateGas) ist << 383     if (material->GetState() == kStateGas) istat = 1;  
373     G4double Zeff = material->GetZ() + pdgb[is    384     G4double Zeff = material->GetZ() + pdgb[istat];
374     EcPdg = pdga[istat] / Zeff;                << 385     EcPdg = pdga[istat]/Zeff;
375     G4cout << "\t\t\t (from Pdg formula : " << << 386     G4cout << "\t\t\t (from Pdg formula : " 
376   }                                            << 387            << std::setw(8) << G4BestUnit(EcPdg,"Energy") << ")";    
377                                                << 388   }
378   const G4double Es = 21.2052 * MeV;           << 389      
379   G4double rMolier1 = Es / ekin, rMolier2 = rM << 390  const G4double Es = 21.2052*MeV;
380   G4cout << "\n  Moliere radius           : "  << 391  G4double rMolier1 = Es/ekin, rMolier2 = rMolier1*radl;
381          << "\t" << std::setw(8) << rMolier1 < << 392  G4cout << "\n  Moliere radius           : "
382          << "= " << std::setw(8) << G4BestUnit << 393         << "\t" << std::setw(8) << rMolier1 << " X0 "   
383                                                << 394         << "= " << std::setw(8) << G4BestUnit(rMolier2,"Length");
384   if (material->GetNumberOfElements() == 1) {  << 395         
385     G4double rMPdg = radl * Es / EcPdg;        << 396  if (material->GetNumberOfElements() == 1) {
386     G4cout << "\t (from Pdg formula : " << std << 397     G4double rMPdg = radl*Es/EcPdg;
387   }                                            << 398     G4cout << "\t (from Pdg formula : " 
                                                   >> 399            << std::setw(8) << G4BestUnit(rMPdg,"Length") << ")";    
                                                   >> 400   }         
388 }                                                 401 }
389                                                   402 
390 //....oooOO0OOooo........oooOO0OOooo........oo    403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
391                                                   404