00001 00017 #include <cassert> 00018 00019 #include "BetheBlochModel.h" 00020 #include "Constants.h" 00021 #include "Utils.h" 00022 00023 ClassImp(BetheBlochModel) 00024 00025 //_________________________________________________________________________________________ 00026 BetheBlochModel::BetheBlochModel() : 00027 ProcessModel() 00028 { 00029 00030 } 00031 //_________________________________________________________________________________________ 00032 BetheBlochModel::BetheBlochModel(const Material & material) : 00033 ProcessModel(material, Process::eIonization) 00034 { 00035 00036 } 00037 //_________________________________________________________________________________________ 00038 BetheBlochModel::~BetheBlochModel() 00039 { 00040 00041 } 00042 //_________________________________________________________________________________________ 00043 ValidityRange_t BetheBlochModel::ValidityRange(void) const 00044 { 00045 ValidityRange_t validity_range; 00046 00047 validity_range.Emin = 0; 00048 validity_range.Emax = 10000000.; // MeV 00049 00050 return validity_range; 00051 } 00052 //_________________________________________________________________________________________ 00053 double BetheBlochModel::dE_dx(double E) const 00054 { 00055 // Calculates ionization dE/dx for muons via Bethe-Bloch formula in MeV / (gr/cm^2) 00056 // Input : E, muon energy in MeV 00057 00058 double a_2 = MuELoss::a_2; 00059 double pi = MuELoss::pi; 00060 double Na = MuELoss::Na; // in 10^-23 00061 double lamda_2 = MuELoss::Le_2; // in 10^22 cm^2 00062 double Z_A = fMaterial.Z() / fMaterial.A(); // in mol/gr 00063 double me = MuELoss::Me; // in MeV 00064 double beta = Utils::BetaFactor(MuELoss::Mm, E); 00065 double beta_2 = beta*beta; 00066 double gamma = Utils::GammaFactor(MuELoss::Mm, E); 00067 double gamma_2 = gamma*gamma; 00068 double I = 1e-6 * fMaterial.IonPotential(); // in MeV (ionization potential) 00069 double I_2 = I*I; // in MeV^2 00070 double Em = MaximumEnergyTransfer(E); // in MeV^2 00071 double Em_2 = Em*Em; // in MeV^2 00072 double E_2 = E*E; // in MeV^2 00073 double d = DensityCorrectionFactor(E); 00074 00075 double de_dx = a_2 * 2*pi * Na * lamda_2 * Z_A * (me / beta_2) * 00076 ( log( 2*me*beta_2*gamma_2*Em/I_2 ) - 2*beta_2 + 0.25*(Em_2/E_2) - d ); 00077 00078 // fix units... Compton wavelength ^2 is multiplied by 10^22 cm^2, 00079 // Avogadro's number is in 10^-23. 00080 // need to multiply by 10 00081 00082 return 10*de_dx; // in MeV/(gr/cm^2) 00083 } 00084 //_________________________________________________________________________________________ 00085 double BetheBlochModel::MaximumEnergyTransfer(double E) const 00086 { 00087 // Calculates the maximum energy transfer to the electron (in MeV) 00088 // Input : muon energy (MeV) 00089 00090 double me = MuELoss::Me; 00091 double me_2 = MuELoss::Me_2; 00092 double mm = MuELoss::Mm; 00093 double mm_2 = MuELoss::Mm_2; 00094 double p = Utils::E2P(mm, E); 00095 double p_2 = p*p; 00096 00097 return 2 * me * p_2 / ( me_2 + mm_2 + 2*me*E ); 00098 } 00099 //___________________________________________________________________________________________ 00100 double BetheBlochModel::DensityCorrectionFactor(double E) const 00101 { 00102 // Calculates the density correction factor delta 00103 // Input : muon energy (MeV) 00104 00105 // Looking at R.M.Sternheimer (BNL), Density Effect for the BetheBlochModel Loss in Various 00106 // Materials, Phys. Rev. 103, 511 (1956): 00107 // For Fe: X1 = 3, X0 = 0.10. a = 0.127, m = 3.29, C = -4.62 00108 00109 // Also Table 1, at YR. 00110 00111 const double X0 = fMaterial.X0(); 00112 const double X1 = fMaterial.X1(); 00113 const double a = fMaterial.a(); 00114 const double m = fMaterial.m(); 00115 const double C = fMaterial.C(); 00116 const double beta = Utils::BetaFactor(MuELoss::Mm, E); 00117 const double gamma = Utils::GammaFactor(MuELoss::Mm, E); 00118 const double X = log10( beta*gamma ); 00119 00120 if(X0 < X && X < X1) return ( 4.6052 * X + a * pow(X1-X,m) + C ); // X e (X0, X1) 00121 if(X > X1) return ( 4.6052 * X + C); // X e (X1, oo) 00122 00123 return 0; // presumably valid in case X < X0 00124 } 00125 //___________________________________________________________________________________________ 00126
1.3.9.1