00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include "G4GammaConversionToMuons.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4UnitsTable.hh"
00041 #include "G4MuonPlus.hh"
00042 #include "G4MuonMinus.hh"
00043
00044
00045
00046 using namespace std;
00047
00048 G4GammaConversionToMuons::G4GammaConversionToMuons(const G4String& processName,
00049 G4ProcessType type):G4VDiscreteProcess (processName, type),
00050 LowestEnergyLimit (4*G4MuonPlus::MuonPlus()->GetPDGMass()),
00051 HighestEnergyLimit(1e21*eV),
00052 CrossSecFactor(1.)
00053 {
00054 SetProcessSubType(15);
00055 MeanFreePath = DBL_MAX;
00056 }
00057
00058
00059
00060
00061
00062 G4GammaConversionToMuons::~G4GammaConversionToMuons()
00063 { }
00064
00065
00066
00067 G4bool G4GammaConversionToMuons::IsApplicable(
00068 const G4ParticleDefinition& particle)
00069 {
00070 return ( &particle == G4Gamma::Gamma() );
00071 }
00072
00073
00074
00075 void G4GammaConversionToMuons::BuildPhysicsTable(const G4ParticleDefinition&)
00076
00077 {
00078 PrintInfoDefinition();
00079 }
00080
00081
00082
00083 G4double G4GammaConversionToMuons::GetMeanFreePath(const G4Track& aTrack,
00084 G4double, G4ForceCondition*)
00085
00086
00087
00088
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 }
00101
00102
00103
00104 G4double G4GammaConversionToMuons::ComputeMeanFreePath(G4double GammaEnergy,
00105 G4Material* aMaterial)
00106
00107
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 }
00123
00124
00125
00126 G4double G4GammaConversionToMuons::GetCrossSectionPerAtom(
00127 const G4DynamicParticle* aDynamicGamma,
00128 G4Element* anElement)
00129
00130
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 }
00139
00140
00141
00142 G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom(
00143 G4double Egam, G4double Z, G4double A)
00144
00145
00146
00147
00148 { static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
00149 static const G4double Mele=electron_mass_c2;
00150 static const G4double Rc=elm_coupling/Mmuon;
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 ;
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;
00163 EgamLast=Egam;
00164
00165 if(Zlast!=Z)
00166 { Zlast=Z;
00167 if(Z==1)
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.);
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);
00186 CrossSection=7./9.*sigfac*log(1.+WMedAppr*CorFuc*Eg);
00187 CrossSection*=CrossSecFactor;
00188 return CrossSection;
00189 }
00190
00191
00192
00193 void G4GammaConversionToMuons::SetCrossSecFactor(G4double fac)
00194
00195 { CrossSecFactor=fac;
00196 G4cout << "The cross section for GammaConversionToMuons is artificially "
00197 << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
00198 }
00199
00200
00201
00202 G4VParticleChange* G4GammaConversionToMuons::PostStepDoIt(
00203 const G4Track& aTrack,
00204 const G4Step& aStep)
00205
00206
00207
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
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
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)
00229 { Zlast=Z;
00230 if(Z==1)
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.);
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
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.;
00264 G4double xxp=1.-4./3.*xPM;
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
00274 G4double t;
00275 G4double psi;
00276 G4double rho;
00277
00278 G4double thetaPlus,thetaMinus,phiHalf;
00279
00280 do
00281 {
00282
00283 G4double C1=C1Num2* GammaMuonInv/xPM;
00284 G4double f1_max=(1.-xPM) / (1.+C1);
00285 G4double f1;
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)
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
00299 G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
00300
00301
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)
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
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
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
00334
00335 G4double phi0=2.*pi*G4UniformRand();
00336 G4double EPlus=xPlus*Egam;
00337 G4double EMinus=xMinus*Egam;
00338
00339
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
00345 MuPlusDirection.rotateUz(GammaDirection);
00346 MuMinusDirection.rotateUz(GammaDirection);
00347 aParticleChange.SetNumberOfSecondaries(2);
00348
00349 G4DynamicParticle* aParticle1= new G4DynamicParticle(
00350 G4MuonPlus::MuonPlus(),MuPlusDirection,EPlus-Mmuon);
00351 aParticleChange.AddSecondary(aParticle1);
00352
00353 G4DynamicParticle* aParticle2= new G4DynamicParticle(
00354 G4MuonMinus::MuonMinus(),MuMinusDirection,EMinus-Mmuon);
00355 aParticleChange.AddSecondary(aParticle2);
00356
00357
00358
00359 aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ;
00360 aParticleChange.ProposeEnergy( 0. ) ;
00361 aParticleChange.ProposeTrackStatus( fStopAndKill ) ;
00362
00363 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
00364 }
00365
00366
00367
00368 G4Element* G4GammaConversionToMuons::SelectRandomAtom(
00369 const G4DynamicParticle* aDynamicGamma,
00370 G4Material* aMaterial)
00371 {
00372
00373
00374 const G4int NumberOfElements = aMaterial->GetNumberOfElements();
00375 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00376 if (NumberOfElements == 1) return (*theElementVector)[0];
00377
00378 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
00379
00380 G4double PartialSumSigma = 0. ;
00381 G4double rval = G4UniformRand()/MeanFreePath;
00382
00383
00384 for ( G4int i=0 ; i < NumberOfElements ; i++ )
00385 { PartialSumSigma += NbOfAtomsPerVolume[i] *
00386 GetCrossSectionPerAtom(aDynamicGamma, (*theElementVector)[i]);
00387 if (rval <= PartialSumSigma) return ((*theElementVector)[i]);
00388 }
00389 G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName()
00390 << "' has no elements, NULL pointer returned." << G4endl;
00391 return NULL;
00392 }
00393
00394
00395
00396 void G4GammaConversionToMuons::PrintInfoDefinition()
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 }
00405
00406