Geant4  10.02
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 91869 2015-08-07 15:21:02Z gcosmo $
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 #include "G4EmProcessSubType.hh"
44 #include "G4NistManager.hh"
45 #include "G4Log.hh"
46 #include "G4Exp.hh"
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
49 
50 using namespace std;
51 
52 static const G4double sqrte=sqrt(exp(1.));
53 static const G4double PowSat=-0.88;
54 
56  G4ProcessType type)
57  : G4VDiscreteProcess (processName, type),
58  Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass()),
59  Rc(elm_coupling/Mmuon),
60  LowestEnergyLimit (4*Mmuon), // 4*Mmuon
61  HighestEnergyLimit(1e21*eV), // ok to 1e21eV=1e12GeV, then LPM suppression
62  CrossSecFactor(1.)
63 {
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
69 
70 // destructor
71 
73 {}
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
76 
78  const G4ParticleDefinition& particle)
79 {
80  return ( &particle == G4Gamma::Gamma() );
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 
86 // Build cross section and mean free path tables
87 { //here no tables, just calling PrintInfoDefinition
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 
95 
96 // returns the photon mean free path in GEANT4 internal units
97 // (MeanFreePath is a private member of the class)
98 
99 {
100  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
101  G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
102  G4Material* aMaterial = aTrack.GetMaterial();
103 
104  if (GammaEnergy <= LowestEnergyLimit)
106  else
107  MeanFreePath = ComputeMeanFreePath(GammaEnergy,aMaterial);
108 
109  return MeanFreePath;
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 
115  G4Material* aMaterial)
116 
117 // computes and returns the photon mean free path in GEANT4 internal units
118 {
119  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
120  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
121 
122  G4double SIGMA = 0 ;
123 
124  for ( size_t i=0 ; i < aMaterial->GetNumberOfElements(); ++i)
125  {
126  G4double AtomicZ = (*theElementVector)[i]->GetZ();
127  G4double AtomicA = (*theElementVector)[i]->GetA()/(g/mole);
128  SIGMA += NbOfAtomsPerVolume[i] *
129  ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA);
130  }
131  return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 
137  const G4DynamicParticle* aDynamicGamma,
138  G4Element* anElement)
139 
140 // gives the total cross section per atom in GEANT4 internal units
141 {
142  G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
143  G4double AtomicZ = anElement->GetZ();
144  G4double AtomicA = anElement->GetN();
145  G4double crossSection =
146  ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA);
147  return crossSection;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
151 
153  G4double Egam, G4double ZZ, G4double)
154 
155 // Calculates the microscopic cross section in GEANT4 internal units.
156 // Total cross section parametrisation from H.Burkhardt
157 // It gives a good description at any energy (from 0 to 10**21 eV)
158 {
159  if(Egam <= LowestEnergyLimit) return 0 ; // below threshold return 0
160 
161  G4int Z = G4lrint(ZZ);
162  G4double CrossSection = 0.0;
164 
165  G4double PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr,
166  Wsatur,sigfac;
167 
168  if(Z==1) // special case of Hydrogen
169  { B=202.4;
170  Dn=1.49;
171  }
172  else
173  { B=183.;
174  Dn=1.54*nist->GetA27(Z);
175  }
176  Zthird=1./nist->GetZ13(Z); // Z**(-1/3)
177  Winfty=B*Zthird*Mmuon/(Dn*electron_mass_c2);
178  WMedAppr=1./(4.*Dn*sqrte*Mmuon);
179  Wsatur=Winfty/WMedAppr;
180  sigfac=4.*fine_structure_const*Z*Z*Rc*Rc;
181  PowThres=1.479+0.00799*Dn;
182  Ecor=-18.+4347./(B*Zthird);
183 
184  G4double CorFuc=1.+.04*G4Log(1.+Ecor/Egam);
185  //G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+
186  // pow(Egam,PowSat),1./PowSat); // threshold and saturation
187  G4double Eg=G4Exp(G4Log(1.-4.*Mmuon/Egam)*PowThres)*
188  G4Exp(G4Log( G4Exp(G4Log(Wsatur)*PowSat)+G4Exp(G4Log(Egam)*PowSat))/PowSat);
189  CrossSection=7./9.*sigfac*G4Log(1.+WMedAppr*CorFuc*Eg);
190  CrossSection*=CrossSecFactor; // increase the CrossSection by (by default 1)
191  return CrossSection;
192 }
193 
194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
195 
197 // Set the factor to artificially increase the cross section
198 {
200  G4cout << "The cross section for GammaConversionToMuons is artificially "
201  << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
205 
207  const G4Track& aTrack,
208  const G4Step& aStep)
209 //
210 // generation of gamma->mu+mu-
211 //
212 {
213  aParticleChange.Initialize(aTrack);
214  G4Material* aMaterial = aTrack.GetMaterial();
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 <= LowestEnergyLimit) {
220  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
221  }
222  G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
223 
224  // select randomly one element constituting the material
225  const G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial);
226  G4int Z = G4lrint(anElement->GetZ());
228 
229  G4double B,Dn;
230  G4double A027 = nist->GetA27(Z);
231 
232  if(Z==1) // special case of Hydrogen
233  { B=202.4;
234  Dn=1.49;
235  }
236  else
237  { B=183.;
238  Dn=1.54*A027;
239  }
240  G4double Zthird=1./nist->GetZ13(Z); // Z**(-1/3)
241  G4double Winfty=B*Zthird*Mmuon/(Dn*electron_mass_c2);
242  G4double C1Num=0.35*A027;
243  G4double C1Num2=C1Num*C1Num;
244  G4double C2Term2=electron_mass_c2/(183.*Zthird*Mmuon);
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/electron_mass_c2;
254  G4double LogWmaxInv=1./G4Log(Winfty*(1.+2.*Ds2*GammaMuonInv)
255  /(1.+2.*sBZ*Mmuon*GammaMuonInv));
256  G4double xPlus,xMinus,xPM,result,W;
257  G4int nn = 0;
258  const G4int nmax = 1000;
259  do
260  { xPlus=xmin+G4UniformRand()*(xmax-xmin);
261  xMinus=1.-xPlus;
262  xPM=xPlus*xMinus;
263  G4double del=Mmuon*Mmuon/(2.*Egam*xPM);
264  W=Winfty*(1.+Ds2*del/Mmuon)/(1.+sBZ*del);
265  if(W<=1. || nn > nmax) { break; } // to avoid negative cross section at xmin
266  G4double xxp=1.-4./3.*xPM; // the main xPlus dependence
267  result=xxp*G4Log(W)*LogWmaxInv;
268  if(result>1.) {
269  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
270  << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
271  }
272  ++nn;
273  if(nn >= nmax) { break; }
274  }
275  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
276  while (G4UniformRand() > result);
277 
278  // now generate the angular variables via the auxilary variables t,psi,rho
279  G4double t;
280  G4double psi;
281  G4double rho;
282 
283  G4double thetaPlus,thetaMinus,phiHalf; // final angular variables
284  nn = 0;
285  do // t, psi, rho generation start (while angle < pi)
286  {
287  //generate t by the rejection method
288  G4double C1=C1Num2* GammaMuonInv/xPM;
289  G4double f1_max=(1.-xPM) / (1.+C1);
290  G4double f1; // the probability density
291  do
292  {
293  ++nn;
294  t=G4UniformRand();
295  f1=(1.-2.*xPM+4.*xPM*t*(1.-t)) / (1.+C1/(t*t));
296  if(f1<0 || f1> f1_max) // should never happend
297  {
298  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
299  << "outside allowed range f1=" << f1 << " is set to zero"
300  << G4endl;
301  f1 = 0.0;
302  }
303  if(nn > nmax) { break; }
304  }
305  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
306  while ( G4UniformRand()*f1_max > f1);
307  // generate psi by the rejection method
308  G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
309 
310  // long version
311  G4double f2;
312  do
313  {
314  ++nn;
315  psi=2.*pi*G4UniformRand();
316  f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+cos(2.*psi));
317  if(f2<0 || f2> f2_max) // should never happend
318  {
319  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
320  << "outside allowed range f2=" << f2 << " is set to zero"
321  << G4endl;
322  f2 = 0.0;
323  }
324  if(nn >= nmax) { break; }
325  }
326  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
327  while ( G4UniformRand()*f2_max > f2);
328 
329  // generate rho by direct transformation
330  G4double C2Term1=GammaMuonInv/(2.*xPM*t);
331  G4double C2=4./sqrt(xPM)*pow(C2Term1*C2Term1+C2Term2*C2Term2,2.);
332  G4double rhomax=1.9/A027*(1./t-1.);
333  G4double beta=G4Log( (C2+rhomax*rhomax*rhomax*rhomax)/C2 );
334  rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4UniformRand())-1. ))*0.25);
335 
336  //now get from t and psi the kinematical variables
337  G4double u=sqrt(1./t-1.);
338  G4double xiHalf=0.5*rho*cos(psi);
339  phiHalf=0.5*rho/u*sin(psi);
340 
341  thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus;
342  thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus;
343 
344  // protection against infinite loop
345  if(nn > nmax) {
346  if(std::abs(thetaPlus)>pi) { thetaPlus = 0.0; }
347  if(std::abs(thetaMinus)>pi) { thetaMinus = 0.0; }
348  }
349 
350  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
351  } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi);
352 
353  // now construct the vectors
354  // azimuthal symmetry, take phi0 at random between 0 and 2 pi
355  G4double phi0=2.*pi*G4UniformRand();
356  G4double EPlus=xPlus*Egam;
357  G4double EMinus=xMinus*Egam;
358 
359  // mu+ mu- directions for gamma in z-direction
360  G4ThreeVector MuPlusDirection ( sin(thetaPlus) *cos(phi0+phiHalf),
361  sin(thetaPlus) *sin(phi0+phiHalf), cos(thetaPlus) );
362  G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf),
363  -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) );
364  // rotate to actual gamma direction
365  MuPlusDirection.rotateUz(GammaDirection);
366  MuMinusDirection.rotateUz(GammaDirection);
368  // create G4DynamicParticle object for the particle1
369  G4DynamicParticle* aParticle1= new G4DynamicParticle(
370  G4MuonPlus::MuonPlus(),MuPlusDirection,EPlus-Mmuon);
371  aParticleChange.AddSecondary(aParticle1);
372  // create G4DynamicParticle object for the particle2
373  G4DynamicParticle* aParticle2= new G4DynamicParticle(
374  G4MuonMinus::MuonMinus(),MuMinusDirection,EMinus-Mmuon);
375  aParticleChange.AddSecondary(aParticle2);
376  //
377  // Kill the incident photon
378  //
382  // Reset NbOfInteractionLengthLeft and return aParticleChange
383  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
384 }
385 
386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
387 
389  const G4DynamicParticle* aDynamicGamma,
390  G4Material* aMaterial)
391 {
392  // select randomly 1 element within the material, invoked by PostStepDoIt
393 
394  const G4int NumberOfElements = aMaterial->GetNumberOfElements();
395  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
396  if (NumberOfElements == 1) return (*theElementVector)[0];
397 
398  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
399 
400  G4double PartialSumSigma = 0. ;
402 
403 
404  for ( G4int i=0 ; i < NumberOfElements ; ++i)
405  { PartialSumSigma += NbOfAtomsPerVolume[i] *
406  GetCrossSectionPerAtom(aDynamicGamma, (*theElementVector)[i]);
407  if (rval <= PartialSumSigma) return ((*theElementVector)[i]);
408  }
409  G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName()
410  << "' has no elements, NULL pointer returned." << G4endl;
411  return NULL;
412 }
413 
414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
415 
417 {
418  G4String comments ="gamma->mu+mu- Bethe Heitler process, SubType= ";
419  G4cout << G4endl << GetProcessName() << ": " << comments
420  << GetProcessSubType() << G4endl;
421  G4cout << " good cross section parametrization from "
422  << G4BestUnit(LowestEnergyLimit,"Energy")
423  << " to " << HighestEnergyLimit/GeV << " GeV for all Z." << G4endl;
424 }
425 
426 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
const double C2
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:99
G4bool IsApplicable(const G4ParticleDefinition &)
const double C1
std::vector< G4Element * > G4ElementVector
static const G4double f2
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
G4double GetN() const
Definition: G4Element.hh:134
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetZ() const
Definition: G4Element.hh:131
G4double GetZ13(G4double Z)
G4double ComputeMeanFreePath(G4double GammaEnergy, G4Material *aMaterial)
G4double GetA27(G4int Z)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
double B(double temperature)
static const G4double PowSat
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
void BuildPhysicsTable(const G4ParticleDefinition &)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
static const G4double sqrte
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
virtual G4double ComputeCrossSectionPerAtom(G4double GammaEnergy, G4double AtomicZ, G4double AtomicA)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
static const double GeV
Definition: G4SIunits.hh:214
G4double GetCrossSectionPerAtom(const G4DynamicParticle *aDynamicGamma, G4Element *anElement)
const G4int nmax
Definition: G4Step.hh:76
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static const G4double f1
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4Material * GetMaterial() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
virtual void Initialize(const G4Track &)
static const double eV
Definition: G4SIunits.hh:212
static const double pi
Definition: G4SIunits.hh:74
int G4lrint(double ad)
Definition: templates.hh:163
void SetNumberOfSecondaries(G4int totSecondaries)
static const double g
Definition: G4SIunits.hh:180
void ProposeEnergy(G4double finalEnergy)
#define DBL_MIN
Definition: templates.hh:75
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
void AddSecondary(G4Track *aSecondary)
static const double mole
Definition: G4SIunits.hh:283
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:100
static const G4double fac
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4Element * SelectRandomAtom(const G4DynamicParticle *aDynamicGamma, G4Material *aMaterial)
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
G4ThreeVector G4ParticleMomentum
#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