Geant4_10
G4GammaConversionToMuons.cc
Go to the documentation of this file.
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 66872 2013-01-15 01:25:57Z japost $
28 //
29 // ------------ G4GammaConversionToMuons physics process ------
30 // by H.Burkhardt, S. Kelner and R. Kokoulin, April 2002
31 //
32 //
33 // 07-08-02: missprint in OR condition in DoIt : f1<0 || f1>f1_max ..etc ...
34 // 25-10-04: migrade to new interfaces of ParticleChange (vi)
35 // ---------------------------------------------------------------------------
36 
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4UnitsTable.hh"
41 #include "G4MuonPlus.hh"
42 #include "G4MuonMinus.hh"
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
45 
46 using namespace std;
47 
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 {
55  MeanFreePath = DBL_MAX;
56 }
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
59 
60 // destructor
61 
63 { }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
66 
68  const G4ParticleDefinition& particle)
69 {
70  return ( &particle == G4Gamma::Gamma() );
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 
76 // Build cross section and mean free path tables
77 { //here no tables, just calling PrintInfoDefinition
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 
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 
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 
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 
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 G4ThreadLocal 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 G4ThreadLocal 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 
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 
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 G4ThreadLocal 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);
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  //
362  // Reset NbOfInteractionLengthLeft and return aParticleChange
363  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
364 }
365 
366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
367 
368 G4Element* 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 
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....
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:99
G4bool IsApplicable(const G4ParticleDefinition &)
std::vector< G4Element * > G4ElementVector
G4double GetKineticEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
const G4String & GetName() const
Definition: G4Material.hh:176
G4double GetZ() const
Definition: G4Element.hh:131
tuple elm_coupling
Definition: hepunit.py:286
G4double G4NeutronHPJENDLHEData::G4double result
G4double ComputeMeanFreePath(G4double GammaEnergy, G4Material *aMaterial)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
#define G4ThreadLocal
Definition: tls.hh:52
G4double GetA() const
Definition: G4Element.hh:138
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
void BuildPhysicsTable(const G4ParticleDefinition &)
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:204
Float_t f1
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
Float_t f2
Float_t Z
Definition: plot.C:39
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
virtual G4double ComputeCrossSectionPerAtom(G4double GammaEnergy, G4double AtomicZ, G4double AtomicA)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4double GetCrossSectionPerAtom(const G4DynamicParticle *aDynamicGamma, G4Element *anElement)
float electron_mass_c2
Definition: hepunit.py:274
Definition: G4Step.hh:76
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4Material * GetMaterial() const
#define C1
virtual void Initialize(const G4Track &)
G4double GetPDGMass() const
void SetNumberOfSecondaries(G4int totSecondaries)
void ProposeEnergy(G4double finalEnergy)
#define DBL_MIN
Definition: templates.hh:75
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
void AddSecondary(G4Track *aSecondary)
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:100
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4GammaConversionToMuons(const G4String &processName="GammaToMuPair", G4ProcessType type=fElectromagnetic)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ForceCondition
#define DBL_MAX
Definition: templates.hh:83
G4int GetProcessSubType() const
Definition: G4VProcess.hh:426
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4ProcessType