source: trunk/source/processes/electromagnetic/highenergy/src/G4GammaConversionToMuons.cc@ 1201

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 15.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4GammaConversionToMuons.cc,v 1.7 2008/10/16 14:29:48 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// ------------ G4GammaConversionToMuons physics process ------
31// by H.Burkhardt, S. Kelner and R. Kokoulin, April 2002
32//
33//
34// 07-08-02: missprint in OR condition in DoIt : f1<0 || f1>f1_max ..etc ...
35// 25-10-04: migrade to new interfaces of ParticleChange (vi)
36// ---------------------------------------------------------------------------
37
38#include "G4GammaConversionToMuons.hh"
39#include "G4EnergyLossTables.hh"
40#include "G4UnitsTable.hh"
41#include "G4MuonPlus.hh"
42#include "G4MuonMinus.hh"
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
45
46using namespace std;
47
48G4GammaConversionToMuons::G4GammaConversionToMuons(const G4String& processName,
49 G4ProcessType type):G4VDiscreteProcess (processName, type),
50 LowestEnergyLimit (4*G4MuonPlus::MuonPlus()->GetPDGMass()), // 4*Mmuon
51 HighestEnergyLimit(1e21*eV), // ok to 1e21eV=1e12GeV, then LPM suppression
52 CrossSecFactor(1.)
53{
54 SetProcessSubType(15);
55}
56
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
58
59// destructor
60
61G4GammaConversionToMuons::~G4GammaConversionToMuons() // (empty) destructor
62{ }
63
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
65
66G4bool G4GammaConversionToMuons::IsApplicable(
67 const G4ParticleDefinition& particle)
68{
69 return ( &particle == G4Gamma::Gamma() );
70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74void G4GammaConversionToMuons::BuildPhysicsTable(const G4ParticleDefinition&)
75// Build cross section and mean free path tables
76{ //here no tables, just calling PrintInfoDefinition
77 PrintInfoDefinition();
78}
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82G4double G4GammaConversionToMuons::GetMeanFreePath(const G4Track& aTrack,
83 G4double, G4ForceCondition*)
84
85// returns the photon mean free path in GEANT4 internal units
86// (MeanFreePath is a private member of the class)
87
88{
89 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
90 G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
91 G4Material* aMaterial = aTrack.GetMaterial();
92
93 if (GammaEnergy < LowestEnergyLimit)
94 MeanFreePath = DBL_MAX;
95 else
96 MeanFreePath = ComputeMeanFreePath(GammaEnergy,aMaterial);
97
98 return MeanFreePath;
99}
100
101//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102
103G4double G4GammaConversionToMuons::ComputeMeanFreePath(G4double GammaEnergy,
104 G4Material* aMaterial)
105
106// computes and returns the photon mean free path in GEANT4 internal units
107{
108 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
109 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
110
111 G4double SIGMA = 0 ;
112
113 for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ )
114 {
115 G4double AtomicZ = (*theElementVector)[i]->GetZ();
116 G4double AtomicA = (*theElementVector)[i]->GetA()/(g/mole);
117 SIGMA += NbOfAtomsPerVolume[i] *
118 ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA);
119 }
120 return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
121}
122
123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124
125G4double G4GammaConversionToMuons::GetCrossSectionPerAtom(
126 const G4DynamicParticle* aDynamicGamma,
127 G4Element* anElement)
128
129// gives the total cross section per atom in GEANT4 internal units
130{
131 G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
132 G4double AtomicZ = anElement->GetZ();
133 G4double AtomicA = anElement->GetA()/(g/mole);
134 G4double crossSection =
135 ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA);
136 return crossSection;
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
140
141G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom(
142 G4double Egam, G4double Z, G4double A)
143
144// Calculates the microscopic cross section in GEANT4 internal units.
145// Total cross section parametrisation from H.Burkhardt
146// It gives a good description at any energy (from 0 to 10**21 eV)
147{ static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
148 static const G4double Mele=electron_mass_c2;
149 static const G4double Rc=elm_coupling/Mmuon; // classical particle radius
150 static const G4double sqrte=sqrt(exp(1.));
151 static const G4double PowSat=-0.88;
152
153 static G4double CrossSection = 0.0 ;
154
155 if ( A < 1. ) return 0;
156 if ( Egam < 4*Mmuon ) return 0 ; // below threshold return 0
157
158 static G4double EgamLast=0,Zlast=0,PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr,
159 Wsatur,sigfac;
160
161 if(Zlast==Z && Egam==EgamLast) return CrossSection; // already calculated
162 EgamLast=Egam;
163
164 if(Zlast!=Z) // new element
165 { Zlast=Z;
166 if(Z==1) // special case of Hydrogen
167 { B=202.4;
168 Dn=1.49;
169 }
170 else
171 { B=183.;
172 Dn=1.54*pow(A,0.27);
173 }
174 Zthird=pow(Z,-1./3.); // Z**(-1/3)
175 Winfty=B*Zthird*Mmuon/(Dn*Mele);
176 WMedAppr=1./(4.*Dn*sqrte*Mmuon);
177 Wsatur=Winfty/WMedAppr;
178 sigfac=4.*fine_structure_const*Z*Z*Rc*Rc;
179 PowThres=1.479+0.00799*Dn;
180 Ecor=-18.+4347./(B*Zthird);
181 }
182 G4double CorFuc=1.+.04*log(1.+Ecor/Egam);
183 G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+
184 pow(Egam,PowSat),1./PowSat); // threshold and saturation
185 CrossSection=7./9.*sigfac*log(1.+WMedAppr*CorFuc*Eg);
186 CrossSection*=CrossSecFactor; // increase the CrossSection by (by default 1)
187 return CrossSection;
188}
189
190//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
191
192void G4GammaConversionToMuons::SetCrossSecFactor(G4double fac)
193// Set the factor to artificially increase the cross section
194{ CrossSecFactor=fac;
195 G4cout << "The cross section for GammaConversionToMuons is artificially "
196 << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
197}
198
199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
200
201G4VParticleChange* G4GammaConversionToMuons::PostStepDoIt(
202 const G4Track& aTrack,
203 const G4Step& aStep)
204//
205// generation of gamma->mu+mu-
206//
207{
208 aParticleChange.Initialize(aTrack);
209 G4Material* aMaterial = aTrack.GetMaterial();
210
211 static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
212 static const G4double Mele=electron_mass_c2;
213 static const G4double sqrte=sqrt(exp(1.));
214
215 // current Gamma energy and direction, return if energy too low
216 const G4DynamicParticle *aDynamicGamma = aTrack.GetDynamicParticle();
217 G4double Egam = aDynamicGamma->GetKineticEnergy();
218 if (Egam < 4*Mmuon) return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
219 G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
220
221 // select randomly one element constituting the material
222 const G4Element& anElement = *SelectRandomAtom(aDynamicGamma, aMaterial);
223 G4double Z = anElement.GetZ();
224 G4double A = anElement.GetA()/(g/mole);
225
226 static G4double Zlast=0,B,Dn,Zthird,Winfty,A027,C1Num2,C2Term2;
227 if(Zlast!=Z) // the element has changed
228 { Zlast=Z;
229 if(Z==1) // special case of Hydrogen
230 { B=202.4;
231 Dn=1.49;
232 }
233 else
234 { B=183.;
235 Dn=1.54*pow(A,0.27);
236 }
237 Zthird=pow(Z,-1./3.); // Z**(-1/3)
238 Winfty=B*Zthird*Mmuon/(Dn*Mele);
239 A027=pow(A,0.27);
240 G4double C1Num=0.35*A027;
241 C1Num2=C1Num*C1Num;
242 C2Term2=Mele/(183.*Zthird*Mmuon);
243 }
244
245 G4double GammaMuonInv=Mmuon/Egam;
246 G4double sqrtx=sqrt(.25-GammaMuonInv);
247 G4double xmax=.5+sqrtx;
248 G4double xmin=.5-sqrtx;
249
250 // generate xPlus according to the differential cross section by rejection
251 G4double Ds2=(Dn*sqrte-2.);
252 G4double sBZ=sqrte*B*Zthird/Mele;
253 G4double LogWmaxInv=1./log(Winfty*(1.+2.*Ds2*GammaMuonInv)
254 /(1.+2.*sBZ*Mmuon*GammaMuonInv));
255 G4double xPlus,xMinus,xPM,result,W;
256 do
257 { xPlus=xmin+G4UniformRand()*(xmax-xmin);
258 xMinus=1.-xPlus;
259 xPM=xPlus*xMinus;
260 G4double del=Mmuon*Mmuon/(2.*Egam*xPM);
261 W=Winfty*(1.+Ds2*del/Mmuon)/(1.+sBZ*del);
262 if(W<1.) W=1.; // to avoid negative cross section at xmin
263 G4double xxp=1.-4./3.*xPM; // the main xPlus dependence
264 result=xxp*log(W)*LogWmaxInv;
265 if(result>1.) {
266 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
267 << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
268 }
269 }
270 while (G4UniformRand() > result);
271
272 // now generate the angular variables via the auxilary variables t,psi,rho
273 G4double t;
274 G4double psi;
275 G4double rho;
276
277 G4double thetaPlus,thetaMinus,phiHalf; // final angular variables
278
279 do // t, psi, rho generation start (while angle < pi)
280 {
281 //generate t by the rejection method
282 G4double C1=C1Num2* GammaMuonInv/xPM;
283 G4double f1_max=(1.-xPM) / (1.+C1);
284 G4double f1; // the probability density
285 do
286 { t=G4UniformRand();
287 f1=(1.-2.*xPM+4.*xPM*t*(1.-t)) / (1.+C1/(t*t));
288 if(f1<0 || f1> f1_max) // should never happend
289 {
290 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
291 << "outside allowed range f1=" << f1 << " is set to zero"
292 << G4endl;
293 f1 = 0.0;
294 }
295 }
296 while ( G4UniformRand()*f1_max > f1);
297 // generate psi by the rejection method
298 G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
299
300 // long version
301 G4double f2;
302 do
303 { psi=2.*pi*G4UniformRand();
304 f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+cos(2.*psi));
305 if(f2<0 || f2> f2_max) // should never happend
306 {
307 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
308 << "outside allowed range f2=" << f2 << " is set to zero"
309 << G4endl;
310 f2 = 0.0;
311 }
312 }
313 while ( G4UniformRand()*f2_max > f2);
314
315 // generate rho by direct transformation
316 G4double C2Term1=GammaMuonInv/(2.*xPM*t);
317 G4double C2=4./sqrt(xPM)*pow(C2Term1*C2Term1+C2Term2*C2Term2,2.);
318 G4double rhomax=1.9/A027*(1./t-1.);
319 G4double beta=log( (C2+pow(rhomax,4.))/C2 );
320 rho=pow(C2 *( exp(beta*G4UniformRand())-1. ) ,0.25);
321
322 //now get from t and psi the kinematical variables
323 G4double u=sqrt(1./t-1.);
324 G4double xiHalf=0.5*rho*cos(psi);
325 phiHalf=0.5*rho/u*sin(psi);
326
327 thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus;
328 thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus;
329
330 } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi);
331
332 // now construct the vectors
333 // azimuthal symmetry, take phi0 at random between 0 and 2 pi
334 G4double phi0=2.*pi*G4UniformRand();
335 G4double EPlus=xPlus*Egam;
336 G4double EMinus=xMinus*Egam;
337
338 // mu+ mu- directions for gamma in z-direction
339 G4ThreeVector MuPlusDirection ( sin(thetaPlus) *cos(phi0+phiHalf),
340 sin(thetaPlus) *sin(phi0+phiHalf), cos(thetaPlus) );
341 G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf),
342 -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) );
343 // rotate to actual gamma direction
344 MuPlusDirection.rotateUz(GammaDirection);
345 MuMinusDirection.rotateUz(GammaDirection);
346 aParticleChange.SetNumberOfSecondaries(2);
347 // create G4DynamicParticle object for the particle1
348 G4DynamicParticle* aParticle1= new G4DynamicParticle(
349 G4MuonPlus::MuonPlus(),MuPlusDirection,EPlus-Mmuon);
350 aParticleChange.AddSecondary(aParticle1);
351 // create G4DynamicParticle object for the particle2
352 G4DynamicParticle* aParticle2= new G4DynamicParticle(
353 G4MuonMinus::MuonMinus(),MuMinusDirection,EMinus-Mmuon);
354 aParticleChange.AddSecondary(aParticle2);
355 //
356 // Kill the incident photon
357 //
358 aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ;
359 aParticleChange.ProposeEnergy( 0. ) ;
360 aParticleChange.ProposeTrackStatus( fStopAndKill ) ;
361 // Reset NbOfInteractionLengthLeft and return aParticleChange
362 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
363}
364
365//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
366
367G4Element* G4GammaConversionToMuons::SelectRandomAtom(
368 const G4DynamicParticle* aDynamicGamma,
369 G4Material* aMaterial)
370{
371 // select randomly 1 element within the material, invoked by PostStepDoIt
372
373 const G4int NumberOfElements = aMaterial->GetNumberOfElements();
374 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
375 if (NumberOfElements == 1) return (*theElementVector)[0];
376
377 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
378
379 G4double PartialSumSigma = 0. ;
380 G4double rval = G4UniformRand()/MeanFreePath;
381
382
383 for ( G4int i=0 ; i < NumberOfElements ; i++ )
384 { PartialSumSigma += NbOfAtomsPerVolume[i] *
385 GetCrossSectionPerAtom(aDynamicGamma, (*theElementVector)[i]);
386 if (rval <= PartialSumSigma) return ((*theElementVector)[i]);
387 }
388 G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName()
389 << "' has no elements, NULL pointer returned." << G4endl;
390 return NULL;
391}
392
393//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
394
395void G4GammaConversionToMuons::PrintInfoDefinition()
396{
397 G4String comments ="gamma->mu+mu- Bethe Heitler process, SubType= ";
398 G4cout << G4endl << GetProcessName() << ": " << comments
399 << GetProcessSubType() << G4endl;
400 G4cout << " good cross section parametrization from "
401 << G4BestUnit(LowestEnergyLimit,"Energy")
402 << " to " << HighestEnergyLimit/GeV << " GeV for all Z." << G4endl;
403}
404
405//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.