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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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