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 9.2.p3)


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