G4ICRU49NuclearStoppingModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:   G4ICRU49NuclearStoppingModel
00034 //
00035 // Author:      V.Ivanchenko 
00036 //
00037 // Creation date: 09.04.2008 from G4MuMscModel
00038 //
00039 // Modifications:
00040 //
00041 //
00042 // Class Description:
00043 //
00044 // Implementation of the model of ICRU'49 nuclear stopping
00045 
00046 // -------------------------------------------------------------------
00047 //
00048 
00049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00050 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00051 
00052 #include "G4ICRU49NuclearStoppingModel.hh"
00053 #include "G4PhysicalConstants.hh"
00054 #include "G4SystemOfUnits.hh"
00055 #include "Randomize.hh"
00056 #include "G4LossTableManager.hh"
00057 #include "G4ParticleChangeForLoss.hh"
00058 #include "G4ElementVector.hh"
00059 #include "G4ProductionCutsTable.hh"
00060 #include "G4Step.hh"
00061 #include "G4Pow.hh"
00062 
00063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00064 
00065 G4double G4ICRU49NuclearStoppingModel::ad[] = {0.0};
00066 G4double G4ICRU49NuclearStoppingModel::ed[] = {0.0};
00067 
00068 using namespace std;
00069 
00070 G4ICRU49NuclearStoppingModel::G4ICRU49NuclearStoppingModel(const G4String& nam) 
00071   : G4VEmModel(nam),lossFlucFlag(false)
00072 {
00073   theZieglerFactor = eV*cm2*1.0e-15;
00074   g4pow = G4Pow::GetInstance();
00075   if(ad[0] == 0.0) { InitialiseNuclearStopping(); }
00076 }
00077 
00078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00079 
00080 G4ICRU49NuclearStoppingModel::~G4ICRU49NuclearStoppingModel()
00081 {}
00082 
00083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00084 
00085 void G4ICRU49NuclearStoppingModel::Initialise(const G4ParticleDefinition*, 
00086                                               const G4DataVector&)
00087 {}
00088 
00089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00090 
00091 void 
00092 G4ICRU49NuclearStoppingModel::SampleSecondaries(
00093                          std::vector<G4DynamicParticle*>*,
00094                          const G4MaterialCutsCouple*,
00095                          const G4DynamicParticle*,
00096                          G4double, G4double)
00097 {}
00098 
00099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00100 
00101 G4double 
00102 G4ICRU49NuclearStoppingModel::ComputeDEDXPerVolume(
00103                          const G4Material* mat,
00104                          const G4ParticleDefinition* p,
00105                          G4double kinEnergy,
00106                          G4double)
00107 {
00108   G4double nloss = 0.0;
00109   if(kinEnergy <= 0.0) { return nloss; }
00110 
00111   // projectile
00112   G4double mass1 = p->GetPDGMass();
00113   G4double z1 = std::fabs(p->GetPDGCharge()/eplus);
00114 
00115   if(kinEnergy*proton_mass_c2/mass1 > z1*z1*MeV) { return nloss; }
00116 
00117   // Projectile nucleus
00118   mass1 /= amu_c2;
00119 
00120   //  loop for the elements in the material
00121   G4int numberOfElements = mat->GetNumberOfElements();
00122   const G4ElementVector* theElementVector = mat->GetElementVector();
00123   const G4double* atomDensity  = mat->GetAtomicNumDensityVector();
00124  
00125   for (G4int iel=0; iel<numberOfElements; iel++) {
00126     const G4Element* element = (*theElementVector)[iel] ;
00127     G4double z2 = element->GetZ();
00128     G4double mass2 = element->GetA()*mole/g ;
00129     nloss += (NuclearStoppingPower(kinEnergy, z1, z2, mass1, mass2))
00130            * atomDensity[iel] ;
00131   }
00132   nloss *= theZieglerFactor;
00133   return nloss;
00134 }
00135 
00136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00137 
00138 G4double 
00139 G4ICRU49NuclearStoppingModel::NuclearStoppingPower(G4double kineticEnergy,
00140                                                    G4double z1, G4double z2,
00141                                                    G4double mass1, G4double mass2)
00142 {
00143   G4double energy = kineticEnergy/keV ;  // energy in keV
00144   G4double nloss = 0.0;
00145   G4double z12 = z1*z2;
00146   G4int iz1 = G4int(z1);
00147   G4int iz2 = G4int(z2);
00148   
00149   G4double rm;
00150   if(iz1 > 1) { rm = (mass1 + mass2)*(g4pow->Z23(iz1) + g4pow->Z23(iz2)); }
00151   else        { rm = (mass1 + mass2)*g4pow->Z13(iz2); }
00152 
00153   G4double er = 32.536 * mass2 * energy / ( z12 * rm ) ;  // reduced energy
00154 
00155   if (er >= ed[0]) { nloss = ad[0]; }
00156   else {
00157     // the table is inverse in energy
00158     for (G4int i=102; i>=0; --i)
00159     {
00160       if (er <= ed[i]) {
00161         nloss = (ad[i] - ad[i+1])*(er - ed[i+1])/(ed[i] - ed[i+1]) + ad[i+1];
00162         break;
00163       }
00164     }
00165   }
00166 
00167   // Stragling
00168   if(lossFlucFlag) {
00169     // G4double sig = 4.0 * mass1 * mass2 / ((mass1 + mass2)*(mass1 + mass2)*
00170     // (4.0 + 0.197*std::pow(er,-1.6991)+6.584*std::pow(er,-1.0494))) ;
00171     G4double sig = 4.0 * mass1 * mass2 / ((mass1 + mass2)*(mass1 + mass2)*
00172                                     (4.0 + 0.197/(er*er) + 6.584/er));
00173 
00174     nloss *= G4RandGauss::shoot(1.0,sig);
00175     lossFlucFlag = false;
00176   }
00177    
00178   nloss *= 8.462 * z12 * mass1 / rm; // Return to [ev/(10^15 atoms/cm^2]
00179 
00180   if ( nloss < 0.0) { nloss = 0.0; }
00181 
00182   return nloss;
00183 }
00184 
00185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00186 
00187 void G4ICRU49NuclearStoppingModel::InitialiseNuclearStopping()
00188 {
00189   const G4double nuca[104][2] = {
00190   { 1.0E+8, 5.831E-8},
00191   { 8.0E+7, 7.288E-8},
00192   { 6.0E+7, 9.719E-8},
00193   { 5.0E+7, 1.166E-7},
00194   { 4.0E+7, 1.457E-7},
00195   { 3.0E+7, 1.942E-7},
00196   { 2.0E+7, 2.916E-7},
00197   { 1.5E+7, 3.887E-7},
00198 
00199   { 1.0E+7, 5.833E-7},
00200   { 8.0E+6, 7.287E-7},
00201   { 6.0E+6, 9.712E-7},
00202   { 5.0E+6, 1.166E-6},
00203   { 4.0E+6, 1.457E-6},
00204   { 3.0E+6, 1.941E-6},
00205   { 2.0E+6, 2.911E-6},
00206   { 1.5E+6, 3.878E-6},
00207 
00208   { 1.0E+6, 5.810E-6},
00209   { 8.0E+5, 7.262E-6},
00210   { 6.0E+5, 9.663E-6},
00211   { 5.0E+5, 1.157E-5},
00212   { 4.0E+5, 1.442E-5},
00213   { 3.0E+5, 1.913E-5},
00214   { 2.0E+5, 2.845E-5},
00215   { 1.5E+5, 3.762E-5},
00216 
00217   { 1.0E+5, 5.554E-5},
00218   { 8.0E+4, 6.866E-5},
00219   { 6.0E+4, 9.020E-5},
00220   { 5.0E+4, 1.070E-4},
00221   { 4.0E+4, 1.319E-4},
00222   { 3.0E+4, 1.722E-4},
00223   { 2.0E+4, 2.499E-4},
00224   { 1.5E+4, 3.248E-4},
00225 
00226   { 1.0E+4, 4.688E-4},
00227   { 8.0E+3, 5.729E-4},
00228   { 6.0E+3, 7.411E-4},
00229   { 5.0E+3, 8.718E-4},
00230   { 4.0E+3, 1.063E-3},
00231   { 3.0E+3, 1.370E-3},
00232   { 2.0E+3, 1.955E-3},
00233   { 1.5E+3, 2.511E-3},
00234 
00235   { 1.0E+3, 3.563E-3},
00236   { 8.0E+2, 4.314E-3},
00237   { 6.0E+2, 5.511E-3},
00238   { 5.0E+2, 6.430E-3},
00239   { 4.0E+2, 7.756E-3},
00240   { 3.0E+2, 9.855E-3},
00241   { 2.0E+2, 1.375E-2},
00242   { 1.5E+2, 1.736E-2},
00243 
00244   { 1.0E+2, 2.395E-2},
00245   { 8.0E+1, 2.850E-2},
00246   { 6.0E+1, 3.552E-2},
00247   { 5.0E+1, 4.073E-2},
00248   { 4.0E+1, 4.802E-2},
00249   { 3.0E+1, 5.904E-2},
00250   { 1.5E+1, 9.426E-2},
00251 
00252   { 1.0E+1, 1.210E-1},
00253   { 8.0E+0, 1.377E-1},
00254   { 6.0E+0, 1.611E-1},
00255   { 5.0E+0, 1.768E-1},
00256   { 4.0E+0, 1.968E-1},
00257   { 3.0E+0, 2.235E-1},
00258   { 2.0E+0, 2.613E-1},
00259   { 1.5E+0, 2.871E-1},
00260 
00261   { 1.0E+0, 3.199E-1},
00262   { 8.0E-1, 3.354E-1},
00263   { 6.0E-1, 3.523E-1},
00264   { 5.0E-1, 3.609E-1},
00265   { 4.0E-1, 3.693E-1},
00266   { 3.0E-1, 3.766E-1},
00267   { 2.0E-1, 3.803E-1},
00268   { 1.5E-1, 3.788E-1},
00269 
00270   { 1.0E-1, 3.711E-1},
00271   { 8.0E-2, 3.644E-1},
00272   { 6.0E-2, 3.530E-1},
00273   { 5.0E-2, 3.444E-1},
00274   { 4.0E-2, 3.323E-1},
00275   { 3.0E-2, 3.144E-1},
00276   { 2.0E-2, 2.854E-1},
00277   { 1.5E-2, 2.629E-1},
00278 
00279   { 1.0E-2, 2.298E-1},
00280   { 8.0E-3, 2.115E-1},
00281   { 6.0E-3, 1.883E-1},
00282   { 5.0E-3, 1.741E-1},
00283   { 4.0E-3, 1.574E-1},
00284   { 3.0E-3, 1.372E-1},
00285   { 2.0E-3, 1.116E-1},
00286   { 1.5E-3, 9.559E-2},
00287 
00288   { 1.0E-3, 7.601E-2},
00289   { 8.0E-4, 6.668E-2},
00290   { 6.0E-4, 5.605E-2},
00291   { 5.0E-4, 5.008E-2},
00292   { 4.0E-4, 4.352E-2},
00293   { 3.0E-4, 3.617E-2},
00294   { 2.0E-4, 2.768E-2},
00295   { 1.5E-4, 2.279E-2},
00296 
00297   { 1.0E-4, 1.723E-2},
00298   { 8.0E-5, 1.473E-2},
00299   { 6.0E-5, 1.200E-2},
00300   { 5.0E-5, 1.052E-2},
00301   { 4.0E-5, 8.950E-3},
00302   { 3.0E-5, 7.246E-3},
00303   { 2.0E-5, 5.358E-3},
00304   { 1.5E-5, 4.313E-3},
00305   { 0.0, 3.166E-3}
00306   };
00307 
00308   for(G4int i=0; i<104; ++i) {
00309     ed[i] = nuca[i][0];
00310     ad[i] = nuca[i][1];
00311   }
00312 }
00313 
00314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

Generated on Mon May 27 17:48:33 2013 for Geant4 by  doxygen 1.4.7