G4MuBetheBlochModel.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 header file
00031 //
00032 //
00033 // File name:     G4MuBetheBlochModel
00034 //
00035 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
00036 //
00037 // Creation date: 09.08.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 // 10-02-04 Calculation of radiative corrections using R.Kokoulin model (V.I)
00046 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
00047 // 12-04-05 Add usage of G4EmCorrections (V.Ivanchenko)
00048 // 13-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
00049 //
00050 
00051 //
00052 // -------------------------------------------------------------------
00053 //
00054 
00055 
00056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00058 
00059 #include "G4MuBetheBlochModel.hh"
00060 #include "G4PhysicalConstants.hh"
00061 #include "G4SystemOfUnits.hh"
00062 #include "Randomize.hh"
00063 #include "G4Electron.hh"
00064 #include "G4LossTableManager.hh"
00065 #include "G4EmCorrections.hh"
00066 #include "G4ParticleChangeForLoss.hh"
00067 
00068 G4double G4MuBetheBlochModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083, 0.5917,
00069                                       0.7628, 0.8983, 0.9801 };
00070                                       
00071 G4double G4MuBetheBlochModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813, 0.1813,
00072                                       0.1569, 0.1112, 0.0506 };
00073 
00074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00075 
00076 using namespace std;
00077 
00078 G4MuBetheBlochModel::G4MuBetheBlochModel(const G4ParticleDefinition* p,
00079                                          const G4String& nam)
00080   : G4VEmModel(nam),
00081   particle(0),
00082   limitKinEnergy(100.*keV),
00083   logLimitKinEnergy(log(limitKinEnergy)),
00084   twoln10(2.0*log(10.0)),
00085   bg2lim(0.0169),
00086   taulim(8.4146e-3),
00087   alphaprime(fine_structure_const/twopi)
00088 {
00089   theElectron = G4Electron::Electron();
00090   corr = G4LossTableManager::Instance()->EmCorrections();
00091   fParticleChange = 0;
00092 
00093   // initial initialisation of memeber should be overwritten
00094   // by SetParticle
00095   mass = massSquare = ratio = 1.0;
00096 
00097   if(p) { SetParticle(p); }
00098 }
00099 
00100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00101 
00102 G4MuBetheBlochModel::~G4MuBetheBlochModel()
00103 {}
00104 
00105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00106 
00107 G4double G4MuBetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
00108                                            const G4MaterialCutsCouple* couple)
00109 {
00110   return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
00111 }
00112 
00113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00114 
00115 G4double G4MuBetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
00116                                                  G4double kinEnergy) 
00117 {
00118   G4double tau  = kinEnergy/mass;
00119   G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
00120                   (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
00121   return tmax;
00122 }
00123 
00124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00125 
00126 void G4MuBetheBlochModel::Initialise(const G4ParticleDefinition* p,
00127                                      const G4DataVector&)
00128 {
00129   if(p) { SetParticle(p); }
00130   if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
00131 }
00132 
00133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00134 
00135 G4double G4MuBetheBlochModel::ComputeCrossSectionPerElectron(
00136                                            const G4ParticleDefinition* p,
00137                                                  G4double kineticEnergy,
00138                                                  G4double cutEnergy,
00139                                                  G4double maxKinEnergy)
00140 {
00141   G4double cross = 0.0;
00142   G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00143   G4double maxEnergy = min(tmax,maxKinEnergy);
00144   if(cutEnergy < maxEnergy) {
00145 
00146     G4double totEnergy = kineticEnergy + mass;
00147     G4double energy2   = totEnergy*totEnergy;
00148     G4double beta2     = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00149 
00150     cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax
00151           + 0.5*(maxEnergy - cutEnergy)/energy2;
00152 
00153     // radiative corrections of R. Kokoulin
00154     if (maxEnergy > limitKinEnergy) {
00155 
00156       G4double logtmax = log(maxEnergy);
00157       G4double logtmin = log(max(cutEnergy,limitKinEnergy));
00158       G4double logstep = logtmax - logtmin;
00159       G4double dcross  = 0.0;
00160 
00161       for (G4int ll=0; ll<8; ll++)
00162       {
00163         G4double ep = exp(logtmin + xgi[ll]*logstep);
00164         G4double a1 = log(1.0 + 2.0*ep/electron_mass_c2);
00165         G4double a3 = log(4.0*totEnergy*(totEnergy - ep)/massSquare);
00166         dcross += wgi[ll]*(1.0/ep - beta2/tmax + 0.5*ep/energy2)*a1*(a3 - a1);
00167       }
00168 
00169       cross += dcross*logstep*alphaprime;
00170     }
00171 
00172     cross *= twopi_mc2_rcl2/beta2;
00173 
00174   }
00175 
00176   //  G4cout << "tmin= " << cutEnergy << " tmax= " << tmax
00177   //         << " cross= " << cross << G4endl;
00178   
00179   return cross;
00180 }
00181 
00182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00183 
00184 G4double G4MuBetheBlochModel::ComputeCrossSectionPerAtom(
00185                                            const G4ParticleDefinition* p,
00186                                                  G4double kineticEnergy,
00187                                                  G4double Z, G4double,
00188                                                  G4double cutEnergy,
00189                                                  G4double maxEnergy)
00190 {
00191   G4double cross = Z*ComputeCrossSectionPerElectron
00192                                          (p,kineticEnergy,cutEnergy,maxEnergy);
00193   return cross;
00194 }
00195 
00196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00197 
00198 G4double G4MuBetheBlochModel::CrossSectionPerVolume(
00199                                            const G4Material* material,
00200                                            const G4ParticleDefinition* p,
00201                                                  G4double kineticEnergy,
00202                                                  G4double cutEnergy,
00203                                                  G4double maxEnergy)
00204 {
00205   G4double eDensity = material->GetElectronDensity();
00206   G4double cross = eDensity*ComputeCrossSectionPerElectron
00207                                          (p,kineticEnergy,cutEnergy,maxEnergy);
00208   return cross;
00209 }
00210 
00211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00212 
00213 G4double G4MuBetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
00214                                                   const G4ParticleDefinition* p,
00215                                                   G4double kineticEnergy,
00216                                                   G4double cut)
00217 {
00218   G4double tmax  = MaxSecondaryEnergy(p, kineticEnergy);
00219   G4double tau   = kineticEnergy/mass;
00220   G4double cutEnergy = min(cut,tmax);
00221   G4double gam   = tau + 1.0;
00222   G4double bg2   = tau * (tau+2.0);
00223   G4double beta2 = bg2/(gam*gam);
00224 
00225   G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
00226   G4double eexc2 = eexc*eexc;
00227   //G4double cden  = material->GetIonisation()->GetCdensity();
00228   //G4double mden  = material->GetIonisation()->GetMdensity();
00229   //G4double aden  = material->GetIonisation()->GetAdensity();
00230   //G4double x0den = material->GetIonisation()->GetX0density();
00231   //G4double x1den = material->GetIonisation()->GetX1density();
00232 
00233   G4double eDensity = material->GetElectronDensity();
00234 
00235   G4double dedx = log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
00236                  -(1.0 + cutEnergy/tmax)*beta2;
00237 
00238   G4double totEnergy = kineticEnergy + mass;
00239   G4double del = 0.5*cutEnergy/totEnergy;
00240   dedx += del*del;
00241 
00242   // density correction
00243   G4double x = log(bg2)/twoln10;
00244   //if ( x >= x0den ) {
00245   //  dedx -= twoln10*x - cden ;
00246   //  if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ;
00247   //}
00248   dedx -= material->GetIonisation()->DensityCorrection(x);
00249 
00250   // shell correction
00251   dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
00252 
00253   // now compute the total ionization loss
00254 
00255   if (dedx < 0.0) dedx = 0.0 ;
00256 
00257   // radiative corrections of R. Kokoulin
00258   if (cutEnergy > limitKinEnergy) {
00259 
00260     G4double logtmax = log(cutEnergy);
00261     G4double logstep = logtmax - logLimitKinEnergy;
00262     G4double dloss = 0.0;
00263     G4double ftot2= 0.5/(totEnergy*totEnergy);
00264 
00265     for (G4int ll=0; ll<8; ll++)
00266     {
00267       G4double ep = exp(logLimitKinEnergy + xgi[ll]*logstep);
00268       G4double a1 = log(1.0 + 2.0*ep/electron_mass_c2);
00269       G4double a3 = log(4.0*totEnergy*(totEnergy - ep)/massSquare);
00270       dloss += wgi[ll]*(1.0 - beta2*ep/tmax + ep*ep*ftot2)*a1*(a3 - a1);
00271     }
00272     dedx += dloss*logstep*alphaprime;
00273   }
00274 
00275   dedx *= twopi_mc2_rcl2*eDensity/beta2;
00276 
00277   //High order corrections
00278   dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
00279 
00280   return dedx;
00281 }
00282 
00283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00284 
00285 void G4MuBetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
00286                                             const G4MaterialCutsCouple*,
00287                                             const G4DynamicParticle* dp,
00288                                             G4double minKinEnergy,
00289                                             G4double maxEnergy)
00290 {
00291   G4double tmax = MaxSecondaryKinEnergy(dp);
00292   G4double maxKinEnergy = min(maxEnergy,tmax);
00293   if(minKinEnergy >= maxKinEnergy) { return; }
00294 
00295   G4double kineticEnergy = dp->GetKineticEnergy();
00296   G4double totEnergy     = kineticEnergy + mass;
00297   G4double etot2         = totEnergy*totEnergy;
00298   G4double beta2         = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
00299  
00300   G4double grej  = 1.;
00301   if(tmax > limitKinEnergy) {
00302     G4double a0    = log(2.*totEnergy/mass);
00303     grej  += alphaprime*a0*a0;
00304   }
00305 
00306   G4double deltaKinEnergy, f;
00307 
00308   // sampling follows ...
00309   do {
00310     G4double q = G4UniformRand();
00311     deltaKinEnergy = minKinEnergy*maxKinEnergy
00312                     /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
00313 
00314 
00315     f = 1.0 - beta2*deltaKinEnergy/tmax 
00316             + 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
00317 
00318     if(deltaKinEnergy > limitKinEnergy) {
00319       G4double a1 = log(1.0 + 2.0*deltaKinEnergy/electron_mass_c2);
00320       G4double a3 = log(4.0*totEnergy*(totEnergy - deltaKinEnergy)/massSquare);
00321       f *= (1. + alphaprime*a1*(a3 - a1));
00322     }
00323 
00324     if(f > grej) {
00325         G4cout << "G4MuBetheBlochModel::SampleSecondary Warning! "
00326                << "Majorant " << grej << " < "
00327                << f << " for edelta= " << deltaKinEnergy
00328                << " tmin= " << minKinEnergy << " max= " << maxKinEnergy
00329                << G4endl;
00330     }
00331 
00332 
00333   } while( grej*G4UniformRand() > f );
00334 
00335   G4double deltaMomentum =
00336            sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
00337   G4double totalMomentum = totEnergy*sqrt(beta2);
00338   G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
00339                                    (deltaMomentum * totalMomentum);
00340 
00341   G4double sint = sqrt(1.0 - cost*cost);
00342 
00343   G4double phi = twopi * G4UniformRand() ;
00344 
00345   G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
00346   G4ThreeVector direction = dp->GetMomentumDirection();
00347   deltaDirection.rotateUz(direction);
00348 
00349   // primary change
00350   kineticEnergy -= deltaKinEnergy;
00351   G4ThreeVector dir = totalMomentum*direction - deltaMomentum*deltaDirection;
00352   direction = dir.unit();
00353   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00354   fParticleChange->SetProposedMomentumDirection(direction);
00355 
00356   // create G4DynamicParticle object for delta ray
00357   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
00358                                                  deltaDirection,deltaKinEnergy);
00359   vdp->push_back(delta);
00360 }
00361 
00362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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