G4BraggModel.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:   G4BraggModel
00034 //
00035 // Author:        Vladimir Ivanchenko
00036 //
00037 // Creation date: 03.01.2002
00038 //
00039 // Modifications: 
00040 //
00041 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
00042 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
00043 // 27-01-03 Make models region aware (V.Ivanchenko)
00044 // 13-02-03 Add name (V.Ivanchenko)
00045 // 04-06-03 Fix compilation warnings (V.Ivanchenko)
00046 // 12-09-04 Add lowestKinEnergy and change order of if in DEDX method (VI)
00047 // 11-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
00048 // 16-06-05 Fix problem of chemical formula (V.Ivantchenko)
00049 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
00050 // 25-04-06 Add stopping data from PSTAR (V.Ivanchenko)
00051 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio, 
00052 //          CorrectionsAlongStep needed for ions(V.Ivanchenko)
00053 
00054 // Class Description:
00055 //
00056 // Implementation of energy loss and delta-electron production by
00057 // slow charged heavy particles
00058 
00059 // -------------------------------------------------------------------
00060 //
00061 
00062 
00063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00065 
00066 #include "G4BraggModel.hh"
00067 #include "G4PhysicalConstants.hh"
00068 #include "G4SystemOfUnits.hh"
00069 #include "Randomize.hh"
00070 #include "G4Electron.hh"
00071 #include "G4ParticleChangeForLoss.hh"
00072 #include "G4LossTableManager.hh"
00073 #include "G4EmCorrections.hh"
00074 
00075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00076 
00077 using namespace std;
00078 
00079 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
00080   : G4VEmModel(nam),
00081     particle(0),
00082     currentMaterial(0),
00083     protonMassAMU(1.007276),
00084     iMolecula(-1),
00085     iPSTAR(-1),
00086     isIon(false),
00087     isInitialised(false)
00088 {
00089   fParticleChange = 0;
00090   SetHighEnergyLimit(2.0*MeV);
00091 
00092   lowestKinEnergy  = 1.0*keV;
00093   theZieglerFactor = eV*cm2*1.0e-15;
00094   theElectron = G4Electron::Electron();
00095   expStopPower125 = 0.0;
00096 
00097   corr = G4LossTableManager::Instance()->EmCorrections();
00098   if(p) { SetParticle(p); }
00099   else  { SetParticle(theElectron); }
00100 }
00101 
00102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00103 
00104 G4BraggModel::~G4BraggModel()
00105 {}
00106 
00107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00108 
00109 void G4BraggModel::Initialise(const G4ParticleDefinition* p,
00110                               const G4DataVector&)
00111 {
00112   if(p != particle) { SetParticle(p); }
00113 
00114   // always false before the run
00115   SetDeexcitationFlag(false);
00116 
00117   if(!isInitialised) {
00118     isInitialised = true;
00119 
00120     G4String pname = particle->GetParticleName();
00121     if(particle->GetParticleType() == "nucleus" && 
00122        pname != "deuteron" && pname != "triton" &&
00123        pname != "alpha+"   && pname != "helium" &&
00124        pname != "hydrogen") { isIon = true; }
00125 
00126     fParticleChange = GetParticleChangeForLoss();
00127   }
00128 }
00129 
00130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00131 
00132 G4double G4BraggModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
00133                                             const G4Material* mat,
00134                                             G4double kineticEnergy)
00135 {
00136   // this method is called only for ions
00137   G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
00138   GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
00139   return q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
00140 }
00141 
00142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00143 
00144 G4double G4BraggModel::GetParticleCharge(const G4ParticleDefinition* p,
00145                                          const G4Material* mat,
00146                                          G4double kineticEnergy)
00147 {
00148   // this method is called only for ions, so no check if it is an ion 
00149   return corr->GetParticleCharge(p,mat,kineticEnergy);
00150 }
00151 
00152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00153 
00154 G4double G4BraggModel::ComputeCrossSectionPerElectron(
00155                                            const G4ParticleDefinition* p,
00156                                                  G4double kineticEnergy,
00157                                                  G4double cutEnergy,
00158                                                  G4double maxKinEnergy)
00159 {
00160   G4double cross     = 0.0;
00161   G4double tmax      = MaxSecondaryEnergy(p, kineticEnergy);
00162   G4double maxEnergy = std::min(tmax,maxKinEnergy);
00163   if(cutEnergy < maxEnergy) {
00164 
00165     G4double energy  = kineticEnergy + mass;
00166     G4double energy2 = energy*energy;
00167     G4double beta2   = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00168     cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
00169 
00170     cross *= twopi_mc2_rcl2*chargeSquare/beta2;
00171   }
00172  //   G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy 
00173  //          << " tmax= " << tmax << " cross= " << cross << G4endl;
00174  
00175   return cross;
00176 }
00177 
00178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00179 
00180 G4double G4BraggModel::ComputeCrossSectionPerAtom(
00181                                            const G4ParticleDefinition* p,
00182                                                  G4double kineticEnergy,
00183                                                  G4double Z, G4double,
00184                                                  G4double cutEnergy,
00185                                                  G4double maxEnergy)
00186 {
00187   G4double cross = Z*ComputeCrossSectionPerElectron
00188                                          (p,kineticEnergy,cutEnergy,maxEnergy);
00189   return cross;
00190 }
00191 
00192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00193 
00194 G4double G4BraggModel::CrossSectionPerVolume(
00195                                            const G4Material* material,
00196                                            const G4ParticleDefinition* p,
00197                                                  G4double kineticEnergy,
00198                                                  G4double cutEnergy,
00199                                                  G4double maxEnergy)
00200 {
00201   G4double eDensity = material->GetElectronDensity();
00202   G4double cross = eDensity*ComputeCrossSectionPerElectron
00203                                          (p,kineticEnergy,cutEnergy,maxEnergy);
00204   return cross;
00205 }
00206 
00207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00208 
00209 G4double G4BraggModel::ComputeDEDXPerVolume(const G4Material* material,
00210                                             const G4ParticleDefinition* p,
00211                                             G4double kineticEnergy,
00212                                             G4double cutEnergy)
00213 {
00214   G4double tmax  = MaxSecondaryEnergy(p, kineticEnergy);
00215   G4double tkin  = kineticEnergy/massRate;
00216   G4double dedx  = 0.0;
00217 
00218   if(tkin < lowestKinEnergy) {
00219     dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
00220   } else {
00221     dedx = DEDX(material, tkin); 
00222   }
00223 
00224   if (cutEnergy < tmax) {
00225 
00226     G4double tau   = kineticEnergy/mass;
00227     G4double gam   = tau + 1.0;
00228     G4double bg2   = tau * (tau+2.0);
00229     G4double beta2 = bg2/(gam*gam);
00230     G4double x     = cutEnergy/tmax;
00231 
00232     dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
00233           * (material->GetElectronDensity())/beta2;
00234   }
00235 
00236   // now compute the total ionization loss
00237 
00238   if (dedx < 0.0) { dedx = 0.0; }
00239 
00240   dedx *= chargeSquare;
00241 
00242   //G4cout << "E(MeV)= " << tkin/MeV << " dedx= " << dedx 
00243   //     << "  " << material->GetName() << G4endl;
00244 
00245   return dedx;
00246 }
00247 
00248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00249 
00250 void G4BraggModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
00251                                      const G4MaterialCutsCouple*,
00252                                      const G4DynamicParticle* dp,
00253                                      G4double xmin,
00254                                      G4double maxEnergy)
00255 {
00256   G4double tmax = MaxSecondaryKinEnergy(dp);
00257   G4double xmax = std::min(tmax, maxEnergy);
00258   if(xmin >= xmax) { return; }
00259 
00260   G4double kineticEnergy = dp->GetKineticEnergy();
00261   G4double energy  = kineticEnergy + mass;
00262   G4double energy2 = energy*energy;
00263   G4double beta2   = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00264   G4double grej    = 1.0;
00265   G4double deltaKinEnergy, f;
00266 
00267   G4ThreeVector direction = dp->GetMomentumDirection();
00268 
00269   // sampling follows ...
00270   do {
00271     G4double q = G4UniformRand();
00272     deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
00273 
00274     f = 1.0 - beta2*deltaKinEnergy/tmax;
00275 
00276     if(f > grej) {
00277         G4cout << "G4BraggModel::SampleSecondary Warning! "
00278                << "Majorant " << grej << " < "
00279                << f << " for e= " << deltaKinEnergy
00280                << G4endl;
00281     }
00282 
00283   } while( grej*G4UniformRand() >= f );
00284 
00285   G4double deltaMomentum =
00286            sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
00287   G4double totMomentum = energy*sqrt(beta2);
00288   G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
00289                                    (deltaMomentum * totMomentum);
00290   if(cost > 1.0) cost = 1.0;
00291   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
00292 
00293   G4double phi = twopi * G4UniformRand() ;
00294 
00295   G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
00296   deltaDirection.rotateUz(direction);
00297 
00298   // Change kinematics of primary particle
00299   kineticEnergy       -= deltaKinEnergy;
00300   G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
00301   finalP               = finalP.unit();
00302   
00303   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00304   fParticleChange->SetProposedMomentumDirection(finalP);
00305 
00306   // create G4DynamicParticle object for delta ray
00307   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,
00308                                                    deltaKinEnergy);
00309 
00310   vdp->push_back(delta);
00311 }
00312 
00313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00314 
00315 G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
00316                                           G4double kinEnergy)
00317 {
00318   if(pd != particle) { SetParticle(pd); }
00319   G4double tau  = kinEnergy/mass;
00320   G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
00321                   (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
00322   return tmax;
00323 }
00324 
00325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00326 
00327 G4bool G4BraggModel::HasMaterial(const G4Material* material)
00328 {
00329   G4String chFormula = material->GetChemicalFormula();
00330   if("" == chFormula) { return false; }
00331 
00332   // ICRU Report N49, 1993. Power's model for He.
00333   const size_t numberOfMolecula = 11;
00334   const G4String molName[numberOfMolecula] = {
00335     "Al_2O_3",                 "CO_2",                      "CH_4",
00336     "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene",  "(C_8H_8)_N",
00337     "C_3H_8",                  "SiO_2",                     "H_2O",
00338     "H_2O-Gas",                "Graphite" } ;
00339   const G4int idxPSTAR[numberOfMolecula] = {
00340     6,  16,  36,
00341     52, 55,  56,
00342     59, 62,  71,
00343     72, 13};
00344 
00345   // Special treatment for water in gas state
00346   const G4State theState = material->GetState() ;
00347 
00348   if( theState == kStateGas && "H_2O" == chFormula) {
00349     chFormula = G4String("H_2O-Gas");
00350   }
00351 
00352   // Search for the material in the table
00353   for (size_t i=0; i<numberOfMolecula; ++i) {
00354     if (chFormula == molName[i]) {
00355       iMolecula = -1;
00356       iPSTAR = idxPSTAR[i];  
00357       return true;
00358     }
00359   }
00360   return false;
00361 }
00362 
00363 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00364 
00365 G4double G4BraggModel::StoppingPower(const G4Material* material,
00366                                            G4double kineticEnergy) 
00367 {
00368   G4double ionloss = 0.0 ;
00369 
00370   if (iMolecula >= 0) {
00371   
00372     // The data and the fit from: 
00373     // ICRU Report N49, 1993. Ziegler's model for protons.
00374     // Proton kinetic energy for parametrisation (keV/amu)  
00375 
00376     G4double T = kineticEnergy/(keV*protonMassAMU) ; 
00377 
00378      static G4double a[11][5] = {
00379    {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2},
00380    {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3}, 
00381    {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3}, 
00382    {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3}, 
00383    {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2}, 
00384    {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1}, 
00385    {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2}, 
00386    {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2},
00387    {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3}, 
00388    {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3},
00389    {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} };
00390 
00391      static G4double atomicWeight[11] = {
00392     101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
00393     104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};       
00394 
00395     if ( T < 10.0 ) {
00396       ionloss = a[iMolecula][0] * sqrt(T) ;
00397     
00398     } else if ( T < 10000.0 ) {
00399       G4double slow  = a[iMolecula][1] * pow(T, 0.45) ;
00400       G4double shigh = log( 1.0 + a[iMolecula][3]/T  
00401                      + a[iMolecula][4]*T ) * a[iMolecula][2]/T ;
00402       ionloss = slow*shigh / (slow + shigh) ;     
00403     } 
00404 
00405     if ( ionloss < 0.0) ionloss = 0.0 ;
00406     if ( 10 == iMolecula ) { 
00407       if (T < 100.0) {
00408         ionloss *= (1.0+0.023+0.0066*log10(T));  
00409       }
00410       else if (T < 700.0) {   
00411         ionloss *=(1.0+0.089-0.0248*log10(T-99.));
00412       } 
00413       else if (T < 10000.0) {    
00414         ionloss *=(1.0+0.089-0.0248*log10(700.-99.));
00415       }
00416     }
00417     ionloss /= atomicWeight[iMolecula];
00418 
00419   // pure material (normally not the case for this function)
00420   } else if(1 == (material->GetNumberOfElements())) {
00421     G4double z = material->GetZ() ;
00422     ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;  
00423   }
00424   
00425   return ionloss;
00426 }
00427 
00428 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00429 
00430 G4double G4BraggModel::ElectronicStoppingPower(G4double z,
00431                                                G4double kineticEnergy) const
00432 {
00433   G4double ionloss ;
00434   G4int i = G4int(z)-1 ;  // index of atom
00435   if(i < 0)  i = 0 ;
00436   if(i > 91) i = 91 ;
00437   
00438   // The data and the fit from: 
00439   // ICRU Report 49, 1993. Ziegler's type of parametrisations.
00440   // Proton kinetic energy for parametrisation (keV/amu)  
00441 
00442   G4double T = kineticEnergy/(keV*protonMassAMU) ; 
00443   
00444   static G4double a[92][5] = {
00445    {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1},
00446    {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2},
00447    {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2},
00448    {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2},
00449    {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2},
00450    {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2},
00451    {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2},
00452    {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2},
00453    {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2},
00454    {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2},
00455        // Z= 11-20
00456    {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2},
00457    {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2},
00458    {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2},
00459    {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2},
00460    {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2},
00461    {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2},
00462    {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2},
00463    {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2},
00464    {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2},
00465    {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2},
00466        // Z= 21-30
00467    {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3},
00468    {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3},
00469    {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3},
00470    {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3},
00471    {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3},
00472    {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3},
00473    {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3},
00474    {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3},
00475    {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2},
00476    {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3},
00477        // Z= 31-40
00478    {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3},
00479    {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3},
00480    {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3},
00481    {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3},
00482    {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3},
00483    {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3},
00484    {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3},
00485    {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3},
00486    {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3},
00487    {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3},
00488        // Z= 41-50
00489    {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3},
00490    {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3},
00491    {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3},
00492    {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3},
00493    {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3},
00494    {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3},
00495    // {5.623,    6.354,    7160.0,   337.6,    0.013940}, // Ag Ziegler77
00496    {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2}, // Ag ICRU49
00497    {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3},
00498    {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3},
00499    {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3},
00500        // Z= 51-60
00501    {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3},
00502    {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3},
00503    {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3},
00504    {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3},
00505    {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3},
00506    {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3},
00507    {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3},
00508    {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3},
00509    {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3},
00510    {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3},
00511        // Z= 61-70
00512    {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3},
00513    {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3},
00514    {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3},
00515    {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3},
00516    {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3},
00517    {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3},
00518    {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3},
00519    {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3},
00520    {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3},
00521    {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3},
00522        // Z= 71-80
00523    {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3},
00524    {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3},
00525    {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3},
00526    {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3},
00527    {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3},
00528    {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3},
00529    {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3},
00530    {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3},
00531    //  {4.856,    5.460,    18320.0,  438.5,    0.002542}, //Ziegler77
00532    {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2}, //ICRU49
00533    {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3},
00534        // Z= 81-90
00535    {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3},
00536    {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3},
00537    {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3},
00538    {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3},
00539    {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3},
00540    {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3},
00541    {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3},
00542    {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3},
00543    {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3},
00544    {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3},
00545        // Z= 91-92
00546    {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3},
00547    {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3}
00548   };
00549 
00550   G4double fac = 1.0 ;
00551 
00552     // Carbon specific case for E < 40 keV
00553   if ( T < 40.0 && 5 == i) {
00554     fac = sqrt(T/40.0) ;
00555     T = 40.0 ;  
00556 
00557     // Free electron gas model
00558   } else if ( T < 10.0 ) { 
00559     fac = sqrt(T*0.1) ;
00560     T =10.0 ;
00561   }
00562 
00563   // Main parametrisation
00564   G4double slow  = a[i][1] * pow(T, 0.45) ;
00565   G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
00566   ionloss = slow*shigh*fac / (slow + shigh) ;     
00567   
00568   if ( ionloss < 0.0) { ionloss = 0.0; }
00569   
00570   return ionloss;
00571 }
00572 
00573 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00574 
00575 G4double G4BraggModel::DEDX(const G4Material* material, G4double kineticEnergy) 
00576 {
00577   G4double eloss = 0.0;
00578 
00579   // check DB
00580   if(material != currentMaterial) {
00581     currentMaterial = material;
00582     iPSTAR    = -1;
00583     iMolecula = -1;
00584     if( !HasMaterial(material) ) { iPSTAR = pstar.GetIndex(material); }
00585 
00586     //G4cout << "%%% " <<material->GetName() << "  iMolecula= " 
00587     //     << iMolecula << "  iPSTAR= " << iPSTAR << G4endl; 
00588 
00589   }
00590 
00591   const G4int numberOfElements = material->GetNumberOfElements();
00592   const G4double* theAtomicNumDensityVector =
00593                                  material->GetAtomicNumDensityVector();
00594   
00595   if( iPSTAR >= 0 ) {
00596     return pstar.GetElectronicDEDX(iPSTAR, kineticEnergy)*material->GetDensity();
00597 
00598   } else if(iMolecula >= 0) {
00599 
00600     eloss = StoppingPower(material, kineticEnergy)*
00601                           material->GetDensity()/amu;
00602 
00603   // Pure material ICRU49 paralmeterisation
00604   } else if(1 == numberOfElements) {
00605 
00606     G4double z = material->GetZ();
00607     eloss = ElectronicStoppingPower(z, kineticEnergy)
00608                                * (material->GetTotNbOfAtomsPerVolume());
00609 
00610 
00611   // Experimental data exist only for kinetic energy 125 keV
00612   } else if( MolecIsInZiegler1988(material) ) { 
00613 
00614     // Loop over elements - calculation based on Bragg's rule 
00615     G4double eloss125 = 0.0 ;
00616     const G4ElementVector* theElementVector =
00617                            material->GetElementVector();
00618   
00619     //  Loop for the elements in the material
00620     for (G4int i=0; i<numberOfElements; i++) {
00621       const G4Element* element = (*theElementVector)[i] ;
00622       G4double z = element->GetZ() ;
00623       eloss    += ElectronicStoppingPower(z,kineticEnergy)
00624                                     * theAtomicNumDensityVector[i] ;
00625       eloss125 += ElectronicStoppingPower(z,125.0*keV)
00626                                     * theAtomicNumDensityVector[i] ;
00627     }      
00628 
00629     // Chemical factor is taken into account
00630     eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
00631  
00632   // Brugg's rule calculation
00633   } else {
00634     const G4ElementVector* theElementVector =
00635                            material->GetElementVector() ;
00636   
00637     //  loop for the elements in the material
00638     for (G4int i=0; i<numberOfElements; i++)
00639     {
00640       const G4Element* element = (*theElementVector)[i] ;
00641       eloss   += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
00642                                    * theAtomicNumDensityVector[i];
00643     }      
00644   }
00645   return eloss*theZieglerFactor;
00646 }
00647 
00648 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00649 
00650 G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material) 
00651 {
00652   // The list of molecules from
00653   // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
00654   // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
00655   
00656   G4String myFormula = G4String(" ") ;
00657   const G4String chFormula = material->GetChemicalFormula() ;
00658   if (myFormula == chFormula ) { return false; }
00659   
00660   //  There are no evidence for difference of stopping power depended on
00661   //  phase of the compound except for water. The stopping power of the 
00662   //  water in gas phase can be predicted using Bragg's rule.
00663   //  
00664   //  No chemical factor for water-gas 
00665    
00666   myFormula = G4String("H_2O") ;
00667   const G4State theState = material->GetState() ;
00668   if( theState == kStateGas && myFormula == chFormula) return false ;
00669     
00670   const size_t numberOfMolecula = 53 ;
00671 
00672   // The coffecient from Table.4 of Ziegler & Manoyan
00673   const G4double HeEff = 2.8735 ;
00674   
00675   static G4String nameOfMol[numberOfMolecula] = {
00676     "H_2O",      "C_2H_4O",    "C_3H_6O",  "C_2H_2",             "C_H_3OH",
00677     "C_2H_5OH",  "C_3H_7OH",   "C_3H_4",   "NH_3",               "C_14H_10",
00678     "C_6H_6",    "C_4H_10",    "C_4H_6",   "C_4H_8O",            "CCl_4",
00679     "CF_4",      "C_6H_8",     "C_6H_12",  "C_6H_10O",           "C_6H_10",
00680     "C_8H_16",   "C_5H_10",    "C_5H_8",   "C_3H_6-Cyclopropane","C_2H_4F_2",
00681     "C_2H_2F_2", "C_4H_8O_2",  "C_2H_6",   "C_2F_6",             "C_2H_6O",
00682     "C_3H_6O",   "C_4H_10O",   "C_2H_4",   "C_2H_4O",            "C_2H_4S",
00683     "SH_2",      "CH_4",       "CCLF_3",   "CCl_2F_2",           "CHCl_2F",
00684     "(CH_3)_2S", "N_2O",       "C_5H_10O", "C_8H_6",             "(CH_2)_N",
00685     "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8",   "C_3H_6-Propylene",   "C_3H_6O",
00686     "C_3H_6S",   "C_4H_4S",    "C_7H_8"
00687   } ;
00688 
00689   static G4double expStopping[numberOfMolecula] = {
00690      66.1,  190.4, 258.7,  42.2, 141.5,
00691     210.9,  279.6, 198.8,  31.0, 267.5,
00692     122.8,  311.4, 260.3, 328.9, 391.3,
00693     206.6,  374.0, 422.0, 432.0, 398.0,
00694     554.0,  353.0, 326.0,  74.6, 220.5,
00695     197.4,  362.0, 170.0, 330.5, 211.3,
00696     262.3,  349.6,  51.3, 187.0, 236.9,
00697     121.9,   35.8, 247.0, 292.6, 268.0,
00698     262.3,   49.0, 398.9, 444.0,  22.91,
00699      68.0,  155.0,  84.0,  74.2, 254.7,
00700     306.8,  324.4, 420.0
00701   } ;
00702 
00703   static G4double expCharge[numberOfMolecula] = {
00704     HeEff, HeEff, HeEff,   1.0, HeEff,
00705     HeEff, HeEff, HeEff,   1.0,   1.0,
00706       1.0, HeEff, HeEff, HeEff, HeEff,
00707     HeEff, HeEff, HeEff, HeEff, HeEff,
00708     HeEff, HeEff, HeEff,   1.0, HeEff,
00709     HeEff, HeEff, HeEff, HeEff, HeEff,
00710     HeEff, HeEff,   1.0, HeEff, HeEff,
00711     HeEff,   1.0, HeEff, HeEff, HeEff,
00712     HeEff,   1.0, HeEff, HeEff,   1.0,
00713       1.0,   1.0,   1.0,   1.0, HeEff,
00714     HeEff, HeEff, HeEff
00715   } ;
00716 
00717   static G4double numberOfAtomsPerMolecula[numberOfMolecula] = {
00718     3.0,  7.0, 10.0,  4.0,  6.0,
00719     9.0, 12.0,  7.0,  4.0, 24.0,
00720     12.0, 14.0, 10.0, 13.0,  5.0,
00721     5.0, 14.0, 18.0, 17.0, 17.0,
00722     24.0, 15.0, 13.0,  9.0,  8.0,
00723     6.0, 14.0,  8.0,  8.0,  9.0,
00724     10.0, 15.0,  6.0,  7.0,  7.0,
00725     3.0,  5.0,  5.0,  5.0,  5.0,
00726     9.0,  3.0, 16.0, 14.0,  3.0,
00727     9.0, 16.0, 11.0,  9.0, 10.0,
00728     10.0,  9.0, 15.0
00729   } ;
00730 
00731   // Search for the compaund in the table
00732   for (size_t i=0; i<numberOfMolecula; i++)
00733     {
00734       if(chFormula == nameOfMol[i]) {
00735         G4double exp125 = expStopping[i] *
00736                           (material->GetTotNbOfAtomsPerVolume()) /
00737                           (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
00738         SetExpStopPower125(exp125);
00739         return true;
00740       }
00741     }
00742   
00743   return false;
00744 }
00745 
00746 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00747 
00748 G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy, 
00749                                       G4double eloss125) const
00750 {
00751   // Approximation of Chemical Factor according to
00752   // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
00753   // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
00754   
00755   G4double gamma    = 1.0 + kineticEnergy/proton_mass_c2 ;    
00756   G4double gamma25  = 1.0 + 25.0*keV /proton_mass_c2 ;
00757   G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
00758   G4double beta     = sqrt(1.0 - 1.0/(gamma*gamma)) ;
00759   G4double beta25   = sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
00760   G4double beta125  = sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
00761   
00762   G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
00763                    (1.0 + exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
00764                    (1.0 + exp( 1.48 * ( beta/beta25    - 7.0 ) ) ) ;
00765 
00766   return factor ;
00767 }
00768 
00769 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00770 

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