#include <G4GammaConversionToMuons.hh>
Inheritance diagram for G4GammaConversionToMuons:
Public Member Functions | |
G4GammaConversionToMuons (const G4String &processName="GammaToMuPair", G4ProcessType type=fElectromagnetic) | |
~G4GammaConversionToMuons () | |
G4bool | IsApplicable (const G4ParticleDefinition &) |
void | BuildPhysicsTable (const G4ParticleDefinition &) |
void | PrintInfoDefinition () |
void | SetCrossSecFactor (G4double fac) |
G4double | GetCrossSecFactor () |
G4double | GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) |
G4double | GetCrossSectionPerAtom (const G4DynamicParticle *aDynamicGamma, G4Element *anElement) |
G4VParticleChange * | PostStepDoIt (const G4Track &aTrack, const G4Step &aStep) |
virtual G4double | ComputeCrossSectionPerAtom (G4double GammaEnergy, G4double AtomicZ, G4double AtomicA) |
G4double | ComputeMeanFreePath (G4double GammaEnergy, G4Material *aMaterial) |
Definition at line 61 of file G4GammaConversionToMuons.hh.
G4GammaConversionToMuons::G4GammaConversionToMuons | ( | const G4String & | processName = "GammaToMuPair" , |
|
G4ProcessType | type = fElectromagnetic | |||
) |
Definition at line 48 of file G4GammaConversionToMuons.cc.
References DBL_MAX, and G4VProcess::SetProcessSubType().
00049 :G4VDiscreteProcess (processName, type), 00050 LowestEnergyLimit (4*G4MuonPlus::MuonPlus()->GetPDGMass()), // 4*Mmuon 00051 HighestEnergyLimit(1e21*eV), // ok to 1e21eV=1e12GeV, then LPM suppression 00052 CrossSecFactor(1.) 00053 { 00054 SetProcessSubType(15); 00055 MeanFreePath = DBL_MAX; 00056 }
G4GammaConversionToMuons::~G4GammaConversionToMuons | ( | ) |
void G4GammaConversionToMuons::BuildPhysicsTable | ( | const G4ParticleDefinition & | ) | [virtual] |
Reimplemented from G4VProcess.
Definition at line 75 of file G4GammaConversionToMuons.cc.
References PrintInfoDefinition().
00077 { //here no tables, just calling PrintInfoDefinition 00078 PrintInfoDefinition(); 00079 }
G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom | ( | G4double | GammaEnergy, | |
G4double | AtomicZ, | |||
G4double | AtomicA | |||
) | [virtual] |
Definition at line 142 of file G4GammaConversionToMuons.cc.
References G4ParticleDefinition::GetPDGMass(), and G4MuonPlus::MuonPlus().
Referenced by ComputeMeanFreePath(), and GetCrossSectionPerAtom().
00148 { static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass(); 00149 static const G4double Mele=electron_mass_c2; 00150 static const G4double Rc=elm_coupling/Mmuon; // classical particle radius 00151 static const G4double sqrte=sqrt(exp(1.)); 00152 static const G4double PowSat=-0.88; 00153 00154 static G4double CrossSection = 0.0 ; 00155 00156 if ( A < 1. ) return 0; 00157 if ( Egam < 4*Mmuon ) return 0 ; // below threshold return 0 00158 00159 static G4double EgamLast=0,Zlast=0,PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr, 00160 Wsatur,sigfac; 00161 00162 if(Zlast==Z && Egam==EgamLast) return CrossSection; // already calculated 00163 EgamLast=Egam; 00164 00165 if(Zlast!=Z) // new element 00166 { Zlast=Z; 00167 if(Z==1) // special case of Hydrogen 00168 { B=202.4; 00169 Dn=1.49; 00170 } 00171 else 00172 { B=183.; 00173 Dn=1.54*pow(A,0.27); 00174 } 00175 Zthird=pow(Z,-1./3.); // Z**(-1/3) 00176 Winfty=B*Zthird*Mmuon/(Dn*Mele); 00177 WMedAppr=1./(4.*Dn*sqrte*Mmuon); 00178 Wsatur=Winfty/WMedAppr; 00179 sigfac=4.*fine_structure_const*Z*Z*Rc*Rc; 00180 PowThres=1.479+0.00799*Dn; 00181 Ecor=-18.+4347./(B*Zthird); 00182 } 00183 G4double CorFuc=1.+.04*log(1.+Ecor/Egam); 00184 G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+ 00185 pow(Egam,PowSat),1./PowSat); // threshold and saturation 00186 CrossSection=7./9.*sigfac*log(1.+WMedAppr*CorFuc*Eg); 00187 CrossSection*=CrossSecFactor; // increase the CrossSection by (by default 1) 00188 return CrossSection; 00189 }
G4double G4GammaConversionToMuons::ComputeMeanFreePath | ( | G4double | GammaEnergy, | |
G4Material * | aMaterial | |||
) |
Definition at line 104 of file G4GammaConversionToMuons.cc.
References ComputeCrossSectionPerAtom(), DBL_MAX, DBL_MIN, G4Material::GetElementVector(), G4Material::GetNumberOfElements(), and G4Material::GetVecNbOfAtomsPerVolume().
Referenced by GetMeanFreePath().
00108 { 00109 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); 00110 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume(); 00111 00112 G4double SIGMA = 0 ; 00113 00114 for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ ) 00115 { 00116 G4double AtomicZ = (*theElementVector)[i]->GetZ(); 00117 G4double AtomicA = (*theElementVector)[i]->GetA()/(g/mole); 00118 SIGMA += NbOfAtomsPerVolume[i] * 00119 ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA); 00120 } 00121 return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX; 00122 }
G4double G4GammaConversionToMuons::GetCrossSecFactor | ( | ) | [inline] |
G4double G4GammaConversionToMuons::GetCrossSectionPerAtom | ( | const G4DynamicParticle * | aDynamicGamma, | |
G4Element * | anElement | |||
) |
Definition at line 126 of file G4GammaConversionToMuons.cc.
References ComputeCrossSectionPerAtom(), G4Element::GetA(), G4DynamicParticle::GetKineticEnergy(), and G4Element::GetZ().
00131 { 00132 G4double GammaEnergy = aDynamicGamma->GetKineticEnergy(); 00133 G4double AtomicZ = anElement->GetZ(); 00134 G4double AtomicA = anElement->GetA()/(g/mole); 00135 G4double crossSection = 00136 ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA); 00137 return crossSection; 00138 }
G4double G4GammaConversionToMuons::GetMeanFreePath | ( | const G4Track & | aTrack, | |
G4double | previousStepSize, | |||
G4ForceCondition * | condition | |||
) | [virtual] |
Implements G4VDiscreteProcess.
Definition at line 83 of file G4GammaConversionToMuons.cc.
References ComputeMeanFreePath(), DBL_MAX, G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), and G4Track::GetMaterial().
00089 { 00090 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle(); 00091 G4double GammaEnergy = aDynamicGamma->GetKineticEnergy(); 00092 G4Material* aMaterial = aTrack.GetMaterial(); 00093 00094 if (GammaEnergy < LowestEnergyLimit) 00095 MeanFreePath = DBL_MAX; 00096 else 00097 MeanFreePath = ComputeMeanFreePath(GammaEnergy,aMaterial); 00098 00099 return MeanFreePath; 00100 }
G4bool G4GammaConversionToMuons::IsApplicable | ( | const G4ParticleDefinition & | ) | [virtual] |
Reimplemented from G4VProcess.
Definition at line 67 of file G4GammaConversionToMuons.cc.
References G4Gamma::Gamma().
00069 { 00070 return ( &particle == G4Gamma::Gamma() ); 00071 }
G4VParticleChange * G4GammaConversionToMuons::PostStepDoIt | ( | const G4Track & | aTrack, | |
const G4Step & | aStep | |||
) | [virtual] |
Reimplemented from G4VDiscreteProcess.
Definition at line 202 of file G4GammaConversionToMuons.cc.
References G4ParticleChange::AddSecondary(), G4VProcess::aParticleChange, C1, fStopAndKill, G4cout, G4endl, G4UniformRand, G4Element::GetA(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4Element::GetZ(), G4ParticleChange::Initialize(), G4MuonMinus::MuonMinus(), G4MuonPlus::MuonPlus(), G4INCL::Math::pi, G4VDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), and G4VParticleChange::SetNumberOfSecondaries().
00208 { 00209 aParticleChange.Initialize(aTrack); 00210 G4Material* aMaterial = aTrack.GetMaterial(); 00211 00212 static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass(); 00213 static const G4double Mele=electron_mass_c2; 00214 static const G4double sqrte=sqrt(exp(1.)); 00215 00216 // current Gamma energy and direction, return if energy too low 00217 const G4DynamicParticle *aDynamicGamma = aTrack.GetDynamicParticle(); 00218 G4double Egam = aDynamicGamma->GetKineticEnergy(); 00219 if (Egam < 4*Mmuon) return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); 00220 G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection(); 00221 00222 // select randomly one element constituting the material 00223 const G4Element& anElement = *SelectRandomAtom(aDynamicGamma, aMaterial); 00224 G4double Z = anElement.GetZ(); 00225 G4double A = anElement.GetA()/(g/mole); 00226 00227 static G4double Zlast=0,B,Dn,Zthird,Winfty,A027,C1Num2,C2Term2; 00228 if(Zlast!=Z) // the element has changed 00229 { Zlast=Z; 00230 if(Z==1) // special case of Hydrogen 00231 { B=202.4; 00232 Dn=1.49; 00233 } 00234 else 00235 { B=183.; 00236 Dn=1.54*pow(A,0.27); 00237 } 00238 Zthird=pow(Z,-1./3.); // Z**(-1/3) 00239 Winfty=B*Zthird*Mmuon/(Dn*Mele); 00240 A027=pow(A,0.27); 00241 G4double C1Num=0.35*A027; 00242 C1Num2=C1Num*C1Num; 00243 C2Term2=Mele/(183.*Zthird*Mmuon); 00244 } 00245 00246 G4double GammaMuonInv=Mmuon/Egam; 00247 G4double sqrtx=sqrt(.25-GammaMuonInv); 00248 G4double xmax=.5+sqrtx; 00249 G4double xmin=.5-sqrtx; 00250 00251 // generate xPlus according to the differential cross section by rejection 00252 G4double Ds2=(Dn*sqrte-2.); 00253 G4double sBZ=sqrte*B*Zthird/Mele; 00254 G4double LogWmaxInv=1./log(Winfty*(1.+2.*Ds2*GammaMuonInv) 00255 /(1.+2.*sBZ*Mmuon*GammaMuonInv)); 00256 G4double xPlus,xMinus,xPM,result,W; 00257 do 00258 { xPlus=xmin+G4UniformRand()*(xmax-xmin); 00259 xMinus=1.-xPlus; 00260 xPM=xPlus*xMinus; 00261 G4double del=Mmuon*Mmuon/(2.*Egam*xPM); 00262 W=Winfty*(1.+Ds2*del/Mmuon)/(1.+sBZ*del); 00263 if(W<1.) W=1.; // to avoid negative cross section at xmin 00264 G4double xxp=1.-4./3.*xPM; // the main xPlus dependence 00265 result=xxp*log(W)*LogWmaxInv; 00266 if(result>1.) { 00267 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" 00268 << " in dSigxPlusGen, result=" << result << " > 1" << G4endl; 00269 } 00270 } 00271 while (G4UniformRand() > result); 00272 00273 // now generate the angular variables via the auxilary variables t,psi,rho 00274 G4double t; 00275 G4double psi; 00276 G4double rho; 00277 00278 G4double thetaPlus,thetaMinus,phiHalf; // final angular variables 00279 00280 do // t, psi, rho generation start (while angle < pi) 00281 { 00282 //generate t by the rejection method 00283 G4double C1=C1Num2* GammaMuonInv/xPM; 00284 G4double f1_max=(1.-xPM) / (1.+C1); 00285 G4double f1; // the probability density 00286 do 00287 { t=G4UniformRand(); 00288 f1=(1.-2.*xPM+4.*xPM*t*(1.-t)) / (1.+C1/(t*t)); 00289 if(f1<0 || f1> f1_max) // should never happend 00290 { 00291 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" 00292 << "outside allowed range f1=" << f1 << " is set to zero" 00293 << G4endl; 00294 f1 = 0.0; 00295 } 00296 } 00297 while ( G4UniformRand()*f1_max > f1); 00298 // generate psi by the rejection method 00299 G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t)); 00300 00301 // long version 00302 G4double f2; 00303 do 00304 { psi=2.*pi*G4UniformRand(); 00305 f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+cos(2.*psi)); 00306 if(f2<0 || f2> f2_max) // should never happend 00307 { 00308 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" 00309 << "outside allowed range f2=" << f2 << " is set to zero" 00310 << G4endl; 00311 f2 = 0.0; 00312 } 00313 } 00314 while ( G4UniformRand()*f2_max > f2); 00315 00316 // generate rho by direct transformation 00317 G4double C2Term1=GammaMuonInv/(2.*xPM*t); 00318 G4double C2=4./sqrt(xPM)*pow(C2Term1*C2Term1+C2Term2*C2Term2,2.); 00319 G4double rhomax=1.9/A027*(1./t-1.); 00320 G4double beta=log( (C2+pow(rhomax,4.))/C2 ); 00321 rho=pow(C2 *( exp(beta*G4UniformRand())-1. ) ,0.25); 00322 00323 //now get from t and psi the kinematical variables 00324 G4double u=sqrt(1./t-1.); 00325 G4double xiHalf=0.5*rho*cos(psi); 00326 phiHalf=0.5*rho/u*sin(psi); 00327 00328 thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus; 00329 thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus; 00330 00331 } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi); 00332 00333 // now construct the vectors 00334 // azimuthal symmetry, take phi0 at random between 0 and 2 pi 00335 G4double phi0=2.*pi*G4UniformRand(); 00336 G4double EPlus=xPlus*Egam; 00337 G4double EMinus=xMinus*Egam; 00338 00339 // mu+ mu- directions for gamma in z-direction 00340 G4ThreeVector MuPlusDirection ( sin(thetaPlus) *cos(phi0+phiHalf), 00341 sin(thetaPlus) *sin(phi0+phiHalf), cos(thetaPlus) ); 00342 G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf), 00343 -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) ); 00344 // rotate to actual gamma direction 00345 MuPlusDirection.rotateUz(GammaDirection); 00346 MuMinusDirection.rotateUz(GammaDirection); 00347 aParticleChange.SetNumberOfSecondaries(2); 00348 // create G4DynamicParticle object for the particle1 00349 G4DynamicParticle* aParticle1= new G4DynamicParticle( 00350 G4MuonPlus::MuonPlus(),MuPlusDirection,EPlus-Mmuon); 00351 aParticleChange.AddSecondary(aParticle1); 00352 // create G4DynamicParticle object for the particle2 00353 G4DynamicParticle* aParticle2= new G4DynamicParticle( 00354 G4MuonMinus::MuonMinus(),MuMinusDirection,EMinus-Mmuon); 00355 aParticleChange.AddSecondary(aParticle2); 00356 // 00357 // Kill the incident photon 00358 // 00359 aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ; 00360 aParticleChange.ProposeEnergy( 0. ) ; 00361 aParticleChange.ProposeTrackStatus( fStopAndKill ) ; 00362 // Reset NbOfInteractionLengthLeft and return aParticleChange 00363 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); 00364 }
void G4GammaConversionToMuons::PrintInfoDefinition | ( | ) |
Definition at line 396 of file G4GammaConversionToMuons.cc.
References G4BestUnit, G4cout, G4endl, G4VProcess::GetProcessName(), and G4VProcess::GetProcessSubType().
Referenced by BuildPhysicsTable().
00397 { 00398 G4String comments ="gamma->mu+mu- Bethe Heitler process, SubType= "; 00399 G4cout << G4endl << GetProcessName() << ": " << comments 00400 << GetProcessSubType() << G4endl; 00401 G4cout << " good cross section parametrization from " 00402 << G4BestUnit(LowestEnergyLimit,"Energy") 00403 << " to " << HighestEnergyLimit/GeV << " GeV for all Z." << G4endl; 00404 }
void G4GammaConversionToMuons::SetCrossSecFactor | ( | G4double | fac | ) |