ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GammaConversionToMuons.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GammaConversionToMuons.cc
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 // ------------ G4GammaConversionToMuons physics process ------
28 // by H.Burkhardt, S. Kelner and R. Kokoulin, April 2002
29 //
30 //
31 // 07-08-02: missprint in OR condition in DoIt : f1<0 || f1>f1_max ..etc ...
32 // 25-10-04: migrade to new interfaces of ParticleChange (vi)
33 // ---------------------------------------------------------------------------
34 
35 #include "G4GammaConversionToMuons.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4UnitsTable.hh"
39 #include "G4MuonPlus.hh"
40 #include "G4MuonMinus.hh"
41 #include "G4EmProcessSubType.hh"
42 #include "G4EmParameters.hh"
43 #include "G4LossTableManager.hh"
44 #include "G4BetheHeitler5DModel.hh"
45 #include "G4Gamma.hh"
46 #include "G4Electron.hh"
47 #include "G4Positron.hh"
48 #include "G4NistManager.hh"
49 #include "G4Log.hh"
50 #include "G4Exp.hh"
51 #include "G4ProductionCutsTable.hh"
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
54 
55 using namespace std;
56 
57 static const G4double sqrte=sqrt(exp(1.));
58 static const G4double PowSat=-0.88;
59 
60 G4GammaConversionToMuons::G4GammaConversionToMuons(const G4String& processName,
61  G4ProcessType type)
62  : G4VDiscreteProcess (processName, type),
63  Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass()),
64  Rc(elm_coupling/Mmuon),
65  LimitEnergy (5.*Mmuon),
66  LowestEnergyLimit (2.*Mmuon),
67  HighestEnergyLimit(1e12*GeV), // ok to 1e12GeV, then LPM suppression
68  Energy5DLimit(0.0),
69  CrossSecFactor(1.),
70  f5Dmodel(nullptr),
71  theGamma(G4Gamma::Gamma()),
72  theMuonPlus(G4MuonPlus::MuonPlus()),
73  theMuonMinus(G4MuonMinus::MuonMinus())
74 {
75  SetProcessSubType(fGammaConversionToMuMu);
76  MeanFreePath = DBL_MAX;
77  fManager = G4LossTableManager::Instance();
78  fManager->Register(this);
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
82 
83 G4GammaConversionToMuons::~G4GammaConversionToMuons()
84 {
85  fManager->DeRegister(this);
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
89 
90 G4bool G4GammaConversionToMuons::IsApplicable(const G4ParticleDefinition& part)
91 {
92  return (&part == theGamma);
93 }
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 
97 void G4GammaConversionToMuons::BuildPhysicsTable(const G4ParticleDefinition& p)
98 // Build cross section and mean free path tables
99 { //here no tables, just calling PrintInfoDefinition
100  Energy5DLimit = G4EmParameters::Instance()->MaxEnergyFor5DMuPair();
101  if(Energy5DLimit > 0.0 && !f5Dmodel) {
102  f5Dmodel = new G4BetheHeitler5DModel();
103  f5Dmodel->SetLeptonPair(theMuonPlus, theMuonMinus);
104  const size_t numElems = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
105  const G4DataVector cuts(numElems);
106  f5Dmodel->Initialise(&p, cuts);
107  }
108  PrintInfoDefinition();
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112 
113 G4double G4GammaConversionToMuons::GetMeanFreePath(const G4Track& aTrack,
114  G4double, G4ForceCondition*)
115 
116 // returns the photon mean free path in GEANT4 internal units
117 // (MeanFreePath is a private member of the class)
118 
119 {
120  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
121  G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
122  const G4Material* aMaterial = aTrack.GetMaterial();
123 
124  MeanFreePath = (GammaEnergy <= LowestEnergyLimit)
125  ? DBL_MAX : ComputeMeanFreePath(GammaEnergy,aMaterial);
126 
127  return MeanFreePath;
128 }
129 
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 
132 G4double
133 G4GammaConversionToMuons::ComputeMeanFreePath(G4double GammaEnergy,
134  const G4Material* aMaterial)
135 
136 // computes and returns the photon mean free path in GEANT4 internal units
137 {
138  if(GammaEnergy <= LowestEnergyLimit) { return DBL_MAX; }
139  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
140  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
141 
142  G4double SIGMA = 0.0;
143  G4double fact = 1.0;
144  G4double e = GammaEnergy;
145  // low energy approximation as in Bethe-Heitler model
146  if(e < LimitEnergy) {
147  G4double y = (e - LowestEnergyLimit)/(LimitEnergy - LowestEnergyLimit);
148  fact = y*y;
149  e = LimitEnergy;
150  }
151 
152  for ( size_t i=0 ; i < aMaterial->GetNumberOfElements(); ++i)
153  {
154  SIGMA += NbOfAtomsPerVolume[i] * fact *
155  ComputeCrossSectionPerAtom(e, (*theElementVector)[i]->GetZasInt());
156  }
157  return (SIGMA > 0.0) ? 1./SIGMA : DBL_MAX;
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161 
162 G4double G4GammaConversionToMuons::GetCrossSectionPerAtom(
163  const G4DynamicParticle* aDynamicGamma,
164  const G4Element* anElement)
165 
166 // gives the total cross section per atom in GEANT4 internal units
167 {
168  return ComputeCrossSectionPerAtom(aDynamicGamma->GetKineticEnergy(),
169  anElement->GetZasInt());
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
173 
174 G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom(
175  G4double Egam, G4int Z)
176 
177 // Calculates the microscopic cross section in GEANT4 internal units.
178 // Total cross section parametrisation from H.Burkhardt
179 // It gives a good description at any energy (from 0 to 10**21 eV)
180 {
181  if(Egam <= LowestEnergyLimit) { return 0.0; }
182 
183  G4NistManager* nist = G4NistManager::Instance();
184 
185  G4double PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr,
186  Wsatur,sigfac;
187 
188  if(Z==1) // special case of Hydrogen
189  { B=202.4;
190  Dn=1.49;
191  }
192  else
193  { B=183.;
194  Dn=1.54*nist->GetA27(Z);
195  }
196  Zthird=1./nist->GetZ13(Z); // Z**(-1/3)
197  Winfty=B*Zthird*Mmuon/(Dn*electron_mass_c2);
198  WMedAppr=1./(4.*Dn*sqrte*Mmuon);
199  Wsatur=Winfty/WMedAppr;
200  sigfac=4.*fine_structure_const*Z*Z*Rc*Rc;
201  PowThres=1.479+0.00799*Dn;
202  Ecor=-18.+4347./(B*Zthird);
203 
204  G4double CorFuc=1.+.04*G4Log(1.+Ecor/Egam);
205  //G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+
206  // pow(Egam,PowSat),1./PowSat); // threshold and saturation
207  G4double Eg=G4Exp(G4Log(1.-4.*Mmuon/Egam)*PowThres)*
208  G4Exp(G4Log( G4Exp(G4Log(Wsatur)*PowSat)+G4Exp(G4Log(Egam)*PowSat))/PowSat);
209  G4double CrossSection=7./9.*sigfac*G4Log(1.+WMedAppr*CorFuc*Eg);
210  CrossSection*=CrossSecFactor; // increase the CrossSection by (by default 1)
211  return CrossSection;
212 }
213 
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
215 
216 void G4GammaConversionToMuons::SetCrossSecFactor(G4double fac)
217 // Set the factor to artificially increase the cross section
218 {
219  CrossSecFactor=fac;
220  G4cout << "The cross section for GammaConversionToMuons is artificially "
221  << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
222 }
223 
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
225 
226 G4VParticleChange* G4GammaConversionToMuons::PostStepDoIt(
227  const G4Track& aTrack,
228  const G4Step& aStep)
229 //
230 // generation of gamma->mu+mu-
231 //
232 {
233  aParticleChange.Initialize(aTrack);
234  const G4Material* aMaterial = aTrack.GetMaterial();
235 
236  // current Gamma energy and direction, return if energy too low
237  const G4DynamicParticle *aDynamicGamma = aTrack.GetDynamicParticle();
238  G4double Egam = aDynamicGamma->GetKineticEnergy();
239  if (Egam <= LowestEnergyLimit) {
240  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
241  }
242  //
243  // Kill the incident photon
244  //
245  aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ;
246  aParticleChange.ProposeEnergy( 0. ) ;
247  aParticleChange.ProposeTrackStatus( fStopAndKill ) ;
248 
249  if (Egam <= Energy5DLimit) {
250  std::vector<G4DynamicParticle*> fvect;
251  f5Dmodel->SampleSecondaries(&fvect, aTrack.GetMaterialCutsCouple(),
252  aTrack.GetDynamicParticle(), 0.0, DBL_MAX);
253  aParticleChange.SetNumberOfSecondaries(fvect.size());
254  for(auto dp : fvect) { aParticleChange.AddSecondary(dp); }
255  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
256  }
257 
258  G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
259 
260  // select randomly one element constituting the material
261  const G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial);
262  G4int Z = anElement->GetZasInt();
263  G4NistManager* nist = G4NistManager::Instance();
264 
265  G4double B,Dn;
266  G4double A027 = nist->GetA27(Z);
267 
268  if(Z==1) // special case of Hydrogen
269  { B=202.4;
270  Dn=1.49;
271  }
272  else
273  { B=183.;
274  Dn=1.54*A027;
275  }
276  G4double Zthird=1./nist->GetZ13(Z); // Z**(-1/3)
277  G4double Winfty=B*Zthird*Mmuon/(Dn*electron_mass_c2);
278 
279  G4double C1Num=0.138*A027;
280  G4double C1Num2=C1Num*C1Num;
281  G4double C2Term2=electron_mass_c2/(183.*Zthird*Mmuon);
282 
283  G4double GammaMuonInv=Mmuon/Egam;
284 
285  // generate xPlus according to the differential cross section by rejection
286  G4double xmin=(Egam < LimitEnergy) ? GammaMuonInv : .5-sqrt(.25-GammaMuonInv);
287  G4double xmax=1.-xmin;
288 
289  G4double Ds2=(Dn*sqrte-2.);
290  G4double sBZ=sqrte*B*Zthird/electron_mass_c2;
291  G4double LogWmaxInv=1./G4Log(Winfty*(1.+2.*Ds2*GammaMuonInv)
292  /(1.+2.*sBZ*Mmuon*GammaMuonInv));
293  G4double xPlus,xMinus,xPM,result,W;
294  G4int nn = 0;
295  const G4int nmax = 1000;
296  do {
297  xPlus=xmin+G4UniformRand()*(xmax-xmin);
298  xMinus=1.-xPlus;
299  xPM=xPlus*xMinus;
300  G4double del=Mmuon*Mmuon/(2.*Egam*xPM);
301  W=Winfty*(1.+Ds2*del/Mmuon)/(1.+sBZ*del);
302  G4double xxp=1.-4./3.*xPM; // the main xPlus dependence
303  result=(xxp > 0.) ? xxp*G4Log(W)*LogWmaxInv : 0.0;
304  if(result>1.) {
305  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
306  << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
307  }
308  ++nn;
309  if(nn >= nmax) { break; }
310  }
311  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
312  while (G4UniformRand() > result);
313 
314  // now generate the angular variables via the auxilary variables t,psi,rho
315  G4double t;
316  G4double psi;
317  G4double rho;
318 
319  G4double a3 = (GammaMuonInv/(2.*xPM));
320  G4double a33 = a3*a3;
321  G4double f1;
322  G4double b1 = 1./(4.*C1Num2);
323  G4double b3 = b1*b1*b1;
324  G4double a21 = a33 + b1;
325 
326  G4double f1_max=-(1.-xPM)*(2.*b1+(a21+a33)*G4Log(a33/a21))/(2*b3);
327 
328  G4double thetaPlus,thetaMinus,phiHalf; // final angular variables
329  nn = 0;
330  // t, psi, rho generation start (while angle < pi)
331  do {
332  //generate t by the rejection method
333  do {
334  ++nn;
335  t=G4UniformRand();
336  G4double a34=a33/(t*t);
337  G4double a22 = a34 + b1;
338  if(std::abs(b1)<0.0001*a34)
339  // special case of a34=a22 because of logarithm accuracy
340  {
341  f1=(1.-2.*xPM+4.*xPM*t*(1.-t))/(12.*a34*a34*a34*a34);
342  }
343  else
344  {
345  f1=-(1.-2.*xPM+4.*xPM*t*(1.-t))*(2.*b1+(a22+a34)*G4Log(a34/a22))/(2*b3);
346  }
347  if(f1<0.0 || f1> f1_max) // should never happend
348  {
349  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
350  << "outside allowed range f1=" << f1
351  << " is set to zero, a34 = "<< a34 << " a22 = "<<a22<<"."
352  << G4endl;
353  f1 = 0.0;
354  }
355  if(nn > nmax) { break; }
356  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
357  } while ( G4UniformRand()*f1_max > f1);
358  // generate psi by the rejection method
359  G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
360  // long version
361  G4double f2;
362  do {
363  ++nn;
364  psi=twopi*G4UniformRand();
365  f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+cos(2.*psi));
366  if(f2<0 || f2> f2_max) // should never happend
367  {
368  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
369  << "outside allowed range f2=" << f2 << " is set to zero"
370  << G4endl;
371  f2 = 0.0;
372  }
373  if(nn >= nmax) { break; }
374  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
375  } while ( G4UniformRand()*f2_max > f2);
376 
377  // generate rho by direct transformation
378  G4double C2Term1=GammaMuonInv/(2.*xPM*t);
379  G4double C22 = C2Term1*C2Term1+C2Term2*C2Term2;
380  G4double C2=4.*C22*C22/sqrt(xPM);
381  G4double rhomax=(1./t-1.)*1.9/A027;
382  G4double beta=G4Log( (C2+rhomax*rhomax*rhomax*rhomax)/C2 );
383  rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4UniformRand())-1. ))*0.25);
384 
385  //now get from t and psi the kinematical variables
386  G4double u=sqrt(1./t-1.);
387  G4double xiHalf=0.5*rho*cos(psi);
388  phiHalf=0.5*rho/u*sin(psi);
389 
390  thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus;
391  thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus;
392 
393  // protection against infinite loop
394  if(nn > nmax) {
395  if(std::abs(thetaPlus)>pi) { thetaPlus = 0.0; }
396  if(std::abs(thetaMinus)>pi) { thetaMinus = 0.0; }
397  }
398 
399  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
400  } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi);
401 
402  // now construct the vectors
403  // azimuthal symmetry, take phi0 at random between 0 and 2 pi
404  G4double phi0=twopi*G4UniformRand();
405  G4double EPlus=xPlus*Egam;
406  G4double EMinus=xMinus*Egam;
407 
408  // mu+ mu- directions for gamma in z-direction
409  G4ThreeVector MuPlusDirection ( sin(thetaPlus) *cos(phi0+phiHalf),
410  sin(thetaPlus) *sin(phi0+phiHalf), cos(thetaPlus) );
411  G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf),
412  -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) );
413  // rotate to actual gamma direction
414  MuPlusDirection.rotateUz(GammaDirection);
415  MuMinusDirection.rotateUz(GammaDirection);
416  aParticleChange.SetNumberOfSecondaries(2);
417  // create G4DynamicParticle object for the particle1
418  G4DynamicParticle* aParticle1 =
419  new G4DynamicParticle(theMuonPlus,MuPlusDirection,EPlus-Mmuon);
420  aParticleChange.AddSecondary(aParticle1);
421  // create G4DynamicParticle object for the particle2
422  G4DynamicParticle* aParticle2 =
423  new G4DynamicParticle(theMuonMinus,MuMinusDirection,EMinus-Mmuon);
424  aParticleChange.AddSecondary(aParticle2);
425  // Reset NbOfInteractionLengthLeft and return aParticleChange
426  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
427 }
428 
429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
430 
431 const G4Element* G4GammaConversionToMuons::SelectRandomAtom(
432  const G4DynamicParticle* aDynamicGamma,
433  const G4Material* aMaterial)
434 {
435  // select randomly 1 element within the material, invoked by PostStepDoIt
436 
437  const G4int NumberOfElements = aMaterial->GetNumberOfElements();
438  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
439  const G4Element* elm = (*theElementVector)[0];
440 
441  if (NumberOfElements > 1) {
442  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
443 
444  G4double PartialSumSigma = 0.;
445  G4double rval = G4UniformRand()/MeanFreePath;
446 
447  for (G4int i=0; i<NumberOfElements; ++i)
448  {
449  elm = (*theElementVector)[i];
450  PartialSumSigma += NbOfAtomsPerVolume[i]
451  *GetCrossSectionPerAtom(aDynamicGamma, elm);
452  if (rval <= PartialSumSigma) { break; }
453  }
454  }
455  return elm;
456 }
457 
458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
459 
460 void G4GammaConversionToMuons::PrintInfoDefinition()
461 {
462  G4String comments ="gamma->mu+mu- Bethe Heitler process, SubType= ";
463  G4cout << G4endl << GetProcessName() << ": " << comments
464  << GetProcessSubType() << G4endl;
465  G4cout << " good cross section parametrization from "
466  << G4BestUnit(LowestEnergyLimit,"Energy")
467  << " to " << HighestEnergyLimit/GeV << " GeV for all Z." << G4endl;
468 }
469 
470 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....