G4WentzelOKandVIxSection.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: G4WentzelOKandVIxSection.cc 69141 2013-04-19 09:36:47Z vnivanch $
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:   G4WentzelOKandVIxSection
00034 //
00035 // Author:      V.Ivanchenko 
00036 //
00037 // Creation date: 09.04.2008 from G4MuMscModel
00038 //
00039 // Modifications:
00040 //
00041 //
00042 
00043 // -------------------------------------------------------------------
00044 //
00045 
00046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00048 
00049 #include "G4WentzelOKandVIxSection.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "Randomize.hh"
00053 #include "G4Electron.hh"
00054 #include "G4Positron.hh"
00055 #include "G4Proton.hh"
00056 #include "G4LossTableManager.hh"
00057 
00058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00059 
00060 G4double G4WentzelOKandVIxSection::ScreenRSquareElec[] = {0.0};
00061 G4double G4WentzelOKandVIxSection::ScreenRSquare[]     = {0.0};
00062 G4double G4WentzelOKandVIxSection::FormFactor[]        = {0.0};
00063 
00064 using namespace std;
00065 
00066 G4WentzelOKandVIxSection::G4WentzelOKandVIxSection() :
00067   numlimit(0.1),
00068   nwarnings(0),
00069   nwarnlimit(50),
00070   alpha2(fine_structure_const*fine_structure_const)
00071 {
00072   fNistManager = G4NistManager::Instance();
00073   fG4pow = G4Pow::GetInstance();
00074   theElectron = G4Electron::Electron();
00075   thePositron = G4Positron::Positron();
00076   theProton   = G4Proton::Proton();
00077   lowEnergyLimit = 1.0*eV;
00078   G4double p0 = electron_mass_c2*classic_electr_radius;
00079   coeff = twopi*p0*p0;
00080   particle = 0;
00081 
00082   // Thomas-Fermi screening radii
00083   // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
00084 
00085   if(0.0 == ScreenRSquare[0]) {
00086     G4double a0 = electron_mass_c2/0.88534; 
00087     G4double constn = 6.937e-6/(MeV*MeV);
00088 
00089     ScreenRSquare[0] = alpha2*a0*a0;
00090     ScreenRSquareElec[0] = ScreenRSquare[0];
00091     for(G4int j=1; j<100; ++j) {
00092       G4double x = a0*fG4pow->Z13(j);
00093       if(1 == j) { ScreenRSquare[j] = 0.5*alpha2*a0*a0; }
00094       else {
00095         ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x;
00096         ScreenRSquareElec[j] = 0.5*alpha2*x*x;
00097       }
00098       x = fNistManager->GetA27(j);
00099       FormFactor[j] = constn*x*x;
00100     } 
00101   }
00102   currentMaterial = 0;
00103   elecXSRatio = factB = factD = formfactA = screenZ = 0.0;
00104   cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
00105 
00106   factB1= 0.5*CLHEP::pi*fine_structure_const;
00107 
00108   tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
00109   ecut = etag = DBL_MAX;
00110   targetZ = 0;
00111   cosThetaMax = 1.0;
00112   targetMass = proton_mass_c2;
00113   particle = 0;
00114 }
00115 
00116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00117 
00118 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection()
00119 {}
00120 
00121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00122 
00123 void G4WentzelOKandVIxSection::Initialise(const G4ParticleDefinition* p, 
00124                                           G4double CosThetaLim)
00125 {
00126   SetupParticle(p);
00127   tkin = mom2 = momCM2 = 0.0;
00128   ecut = etag = DBL_MAX;
00129   targetZ = 0;
00130   cosThetaMax = CosThetaLim;
00131   G4double a = G4LossTableManager::Instance()->FactorForAngleLimit()
00132     *CLHEP::hbarc/CLHEP::fermi;
00133   factorA2 = 0.5*a*a;
00134   currentMaterial = 0;
00135   
00136   //G4cout << "G4WentzelOKandVIxSection::Initialise  mass= " << mass
00137   //     << "  " << p->GetParticleName() << G4endl; 
00138   
00139 }
00140 
00141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00142 
00143 void G4WentzelOKandVIxSection::SetupParticle(const G4ParticleDefinition* p)
00144 {
00145   particle = p;
00146   mass = particle->GetPDGMass();
00147   spin = particle->GetPDGSpin();
00148   if(0.0 != spin) { spin = 0.5; }
00149   G4double q = std::fabs(particle->GetPDGCharge()/eplus);
00150   chargeSquare = q*q;
00151   charge3 = chargeSquare*q;
00152   tkin = 0.0;
00153   currentMaterial = 0;
00154   targetZ = 0;
00155 }
00156 
00157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00158   
00159 G4double
00160 G4WentzelOKandVIxSection::SetupTarget(G4int Z, G4double cut)
00161 {
00162   G4double cosTetMaxNuc2 = cosTetMaxNuc;
00163   if(Z != targetZ || tkin != etag) {
00164     etag    = tkin; 
00165     targetZ = Z;
00166     if(targetZ > 99) { targetZ = 99; }
00167     SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2);
00168     //G4double tmass2 = targetMass*targetMass;
00169     //G4double etot = tkin + mass;
00170     //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass;
00171     //momCM2 = mom2*tmass2/invmass2;
00172     //gam0pcmp = (etot + targetMass)*targetMass/invmass2;
00173     //pcmp2    = tmass2/invmass2;
00174 
00175     kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
00176 
00177     if(1 == Z) {
00178       screenZ = ScreenRSquare[targetZ]/mom2;
00179     } else if(mass > MeV) {
00180       screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
00181         ScreenRSquare[targetZ]/mom2;
00182     } else {
00183       G4double tau = tkin/mass;
00184       screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
00185           *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
00186         ScreenRSquareElec[targetZ]/mom2;
00187     }
00188     if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
00189       cosTetMaxNuc2 = 0.0;
00190     }
00191     formfactA = FormFactor[targetZ]*mom2;
00192 
00193     cosTetMaxElec = 1.0;
00194     ComputeMaxElectronScattering(cut); 
00195   }
00196   return cosTetMaxNuc2;
00197 } 
00198 
00199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00200 
00201 G4double 
00202 G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax)
00203 {
00204   G4double xsec = 0.0;
00205   if(cosTMax >= 1.0) { return xsec; }
00206  
00207   G4double xSection = 0.0;
00208   G4double x = 0; 
00209   G4double y = 0;
00210   G4double x1= 0;
00211   G4double x2= 0;
00212   G4double xlog = 0.0;
00213 
00214   G4double costm = std::max(cosTMax,cosTetMaxElec); 
00215   G4double fb = screenZ*factB;
00216 
00217   // scattering off electrons
00218   if(costm < 1.0) {
00219     x = (1.0 - costm)/screenZ;
00220     if(x < numlimit) { 
00221       x2 = 0.5*x*x;
00222       y  = x2*(1.0 - 1.3333333*x + 3*x2);
00223       if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
00224     } else { 
00225       x1= x/(1 + x);
00226       xlog = log(1.0 + x);  
00227       y = xlog - x1; 
00228       if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
00229     }
00230 
00231     if(y < 0.0) {
00232       ++nwarnings;
00233       if(nwarnings < nwarnlimit) {
00234         G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
00235                << G4endl;
00236         G4cout << "y= " << y 
00237                << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
00238                << " Z= " << targetZ << "  " 
00239                << particle->GetParticleName() << G4endl;
00240         G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ 
00241                << " x= " << x << G4endl;
00242       }
00243       y = 0.0;
00244     }
00245     xSection = y;
00246   }
00247   /* 
00248   G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 
00249          << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
00250          << " cut(MeV)= " << ecut/MeV  
00251          << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 
00252          << " zmaxN= " << (1.0 - cosThetaMax)/screenZ 
00253          << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
00254   */
00255   // scattering off nucleus
00256   if(cosTMax < 1.0) {
00257     x = (1.0 - cosTMax)/screenZ;
00258     if(x < numlimit) { 
00259       x2 = 0.5*x*x;
00260       y  = x2*(1.0 - 1.3333333*x + 3*x2); 
00261       if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
00262     } else { 
00263       x1= x/(1 + x);
00264       xlog = log(1.0 + x);  
00265       y = xlog - x1; 
00266       if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
00267     }
00268 
00269     if(y < 0.0) {
00270       ++nwarnings;
00271       if(nwarnings < nwarnlimit) {
00272         G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
00273                << G4endl;
00274         G4cout << "y= " << y 
00275                << " e(MeV)= " << tkin << " Z= " << targetZ << "  " 
00276                << particle->GetParticleName() << G4endl;
00277         G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ 
00278                << " x= " << " x1= " << x1 <<G4endl;
00279       }
00280       y = 0.0;
00281     }
00282     xSection += y*targetZ; 
00283   }
00284   xSection *= kinFactor;
00285   /*
00286   G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 
00287          << " screenZ= " << screenZ << " formF= " << formfactA 
00288          << " for " << particle->GetParticleName() 
00289          << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
00290          << " x= " << x 
00291          << G4endl;
00292   */
00293   return xSection; 
00294 }
00295 
00296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00297 
00298 G4ThreeVector
00299 G4WentzelOKandVIxSection::SampleSingleScattering(G4double cosTMin,
00300                                                  G4double cosTMax,
00301                                                  G4double elecRatio)
00302 {
00303   G4ThreeVector v(0.0,0.0,1.0);
00304  
00305   G4double formf = formfactA;
00306   G4double cost1 = cosTMin;
00307   G4double cost2 = cosTMax;
00308   if(elecRatio > 0.0) {
00309     if(G4UniformRand() <= elecRatio) {
00310       formf = 0.0;
00311       cost1 = std::max(cost1,cosTetMaxElec);
00312       cost2 = std::max(cost2,cosTetMaxElec);
00313     }
00314   }
00315   if(cost1 < cost2) { return v; }
00316 
00317   G4double w1 = 1. - cost1 + screenZ;
00318   G4double w2 = 1. - cost2 + screenZ;
00319   G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
00320 
00321   //G4double fm =  1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass);
00322   G4double fm =  1.0 + formf*z1;
00323   //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm );
00324   G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm );
00325   // "false" scattering
00326   if( G4UniformRand() > grej ) { return v; }
00327     // } 
00328   G4double cost = 1.0 - z1;
00329 
00330   if(cost > 1.0)       { cost = 1.0; }
00331   else if(cost < -1.0) { cost =-1.0; }
00332   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
00333   //G4cout << "sint= " << sint << G4endl;
00334   G4double phi  = twopi*G4UniformRand();
00335   G4double vx1 = sint*cos(phi);
00336   G4double vy1 = sint*sin(phi);
00337 
00338   // only direction is changed
00339   v.set(vx1,vy1,cost);
00340   return v;
00341 }
00342 
00343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00344 
00345 void 
00346 G4WentzelOKandVIxSection::ComputeMaxElectronScattering(G4double cutEnergy)
00347 {
00348   if(mass > MeV) {
00349     G4double ratio = electron_mass_c2/mass;
00350     G4double tau = tkin/mass;
00351     G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
00352       (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
00353     //tmax = std::min(tmax, targetZ*targetZ*10*eV); 
00354     cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
00355   } else {
00356 
00357     G4double tmax = tkin;
00358     if(particle == theElectron) { tmax *= 0.5; }
00359     //tmax = std::min(tmax, targetZ*targetZ*10*eV); 
00360     G4double t = std::min(cutEnergy, tmax);
00361     G4double mom21 = t*(t + 2.0*electron_mass_c2);
00362     G4double t1 = tkin - t;
00363     //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= " 
00364     //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
00365     if(t1 > 0.0) {
00366       G4double mom22 = t1*(t1 + 2.0*mass);
00367       G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
00368       if(ctm <  1.0) { cosTetMaxElec = ctm; }
00369       if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; }
00370     }
00371   }
00372 }
00373 
00374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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