Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4GammaConversionToMuons.cc 106961 2017-10-31 08:36:29Z 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 
71 {}
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
74 
76  const G4ParticleDefinition& particle)
77 {
78  return ( &particle == G4Gamma::Gamma() );
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 
84 // Build cross section and mean free path tables
85 { //here no tables, just calling PrintInfoDefinition
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
93 
94 // returns the photon mean free path in GEANT4 internal units
95 // (MeanFreePath is a private member of the class)
96 
97 {
98  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
99  G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
100  const G4Material* aMaterial = aTrack.GetMaterial();
101 
102  MeanFreePath = (GammaEnergy <= LowestEnergyLimit)
103  ? DBL_MAX : ComputeMeanFreePath(GammaEnergy,aMaterial);
104 
105  return MeanFreePath;
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109 
110 G4double
112  const G4Material* aMaterial)
113 
114 // computes and returns the photon mean free path in GEANT4 internal units
115 {
116  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
117  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
118 
119  G4double SIGMA = 0.0;
120 
121  for ( size_t i=0 ; i < aMaterial->GetNumberOfElements(); ++i)
122  {
123  SIGMA += NbOfAtomsPerVolume[i] *
124  ComputeCrossSectionPerAtom(GammaEnergy,
125  (*theElementVector)[i]->GetZasInt());
126  }
127  return (SIGMA > 0.0) ? 1./SIGMA : DBL_MAX;
128 }
129 
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 
133  const G4DynamicParticle* aDynamicGamma,
134  const G4Element* anElement)
135 
136 // gives the total cross section per atom in GEANT4 internal units
137 {
138  return ComputeCrossSectionPerAtom(aDynamicGamma->GetKineticEnergy(),
139  anElement->GetZasInt());
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
143 
145  G4double Egam, G4int Z)
146 
147 // Calculates the microscopic cross section in GEANT4 internal units.
148 // Total cross section parametrisation from H.Burkhardt
149 // It gives a good description at any energy (from 0 to 10**21 eV)
150 {
151  if(Egam <= LowestEnergyLimit) return 0.0; // below threshold return 0
152 
153  G4double CrossSection = 0.0;
155 
156  G4double PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr,
157  Wsatur,sigfac;
158 
159  if(Z==1) // special case of Hydrogen
160  { B=202.4;
161  Dn=1.49;
162  }
163  else
164  { B=183.;
165  Dn=1.54*nist->GetA27(Z);
166  }
167  Zthird=1./nist->GetZ13(Z); // Z**(-1/3)
168  Winfty=B*Zthird*Mmuon/(Dn*electron_mass_c2);
169  WMedAppr=1./(4.*Dn*sqrte*Mmuon);
170  Wsatur=Winfty/WMedAppr;
171  sigfac=4.*fine_structure_const*Z*Z*Rc*Rc;
172  PowThres=1.479+0.00799*Dn;
173  Ecor=-18.+4347./(B*Zthird);
174 
175  G4double CorFuc=1.+.04*G4Log(1.+Ecor/Egam);
176  //G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+
177  // pow(Egam,PowSat),1./PowSat); // threshold and saturation
178  G4double Eg=G4Exp(G4Log(1.-4.*Mmuon/Egam)*PowThres)*
179  G4Exp(G4Log( G4Exp(G4Log(Wsatur)*PowSat)+G4Exp(G4Log(Egam)*PowSat))/PowSat);
180  CrossSection=7./9.*sigfac*G4Log(1.+WMedAppr*CorFuc*Eg);
181  CrossSection*=CrossSecFactor; // increase the CrossSection by (by default 1)
182  return CrossSection;
183 }
184 
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
186 
188 // Set the factor to artificially increase the cross section
189 {
191  G4cout << "The cross section for GammaConversionToMuons is artificially "
192  << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
193 }
194 
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
196 
198  const G4Track& aTrack,
199  const G4Step& aStep)
200 //
201 // generation of gamma->mu+mu-
202 //
203 {
204  aParticleChange.Initialize(aTrack);
205  const G4Material* aMaterial = aTrack.GetMaterial();
206 
207  // current Gamma energy and direction, return if energy too low
208  const G4DynamicParticle *aDynamicGamma = aTrack.GetDynamicParticle();
209  G4double Egam = aDynamicGamma->GetKineticEnergy();
210  if (Egam <= LowestEnergyLimit) {
211  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
212  }
213  G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
214 
215  // select randomly one element constituting the material
216  const G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial);
217  G4int Z = anElement->GetZasInt();
219 
220  G4double B,Dn;
221  G4double A027 = nist->GetA27(Z);
222 
223  if(Z==1) // special case of Hydrogen
224  { B=202.4;
225  Dn=1.49;
226  }
227  else
228  { B=183.;
229  Dn=1.54*A027;
230  }
231  G4double Zthird=1./nist->GetZ13(Z); // Z**(-1/3)
232  G4double Winfty=B*Zthird*Mmuon/(Dn*electron_mass_c2);
233 
234  G4double C1Num=0.138*A027;
235  G4double C1Num2=C1Num*C1Num;
236  G4double C2Term2=electron_mass_c2/(183.*Zthird*Mmuon);
237 
238  G4double GammaMuonInv=Mmuon/Egam;
239  G4double sqrtx=sqrt(.25-GammaMuonInv);
240  G4double xmax=.5+sqrtx;
241  G4double xmin=.5-sqrtx;
242 
243  // generate xPlus according to the differential cross section by rejection
244  G4double Ds2=(Dn*sqrte-2.);
245  G4double sBZ=sqrte*B*Zthird/electron_mass_c2;
246  G4double LogWmaxInv=1./G4Log(Winfty*(1.+2.*Ds2*GammaMuonInv)
247  /(1.+2.*sBZ*Mmuon*GammaMuonInv));
248  G4double xPlus,xMinus,xPM,result,W;
249  G4int nn = 0;
250  const G4int nmax = 1000;
251  do
252  { xPlus=xmin+G4UniformRand()*(xmax-xmin);
253  xMinus=1.-xPlus;
254  xPM=xPlus*xMinus;
255  G4double del=Mmuon*Mmuon/(2.*Egam*xPM);
256  W=Winfty*(1.+Ds2*del/Mmuon)/(1.+sBZ*del);
257  if(W<=1. || nn > nmax) { break; } // to avoid negative cross section at xmin
258  G4double xxp=1.-4./3.*xPM; // the main xPlus dependence
259  result=xxp*G4Log(W)*LogWmaxInv;
260  if(result>1.) {
261  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
262  << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
263  }
264  ++nn;
265  if(nn >= nmax) { break; }
266  }
267  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
268  while (G4UniformRand() > result);
269 
270  // now generate the angular variables via the auxilary variables t,psi,rho
271  G4double t;
272  G4double psi;
273  G4double rho;
274 
275  G4double a3 = (GammaMuonInv/(2.*xPM));
276  G4double a33 = a3*a3;
277  G4double f1;
278  G4double b1 = 1./(4.*C1Num2);
279  G4double b3 = b1*b1*b1;
280  G4double a21 = a33 + b1;
281 
282  G4double f1_max=-(1.-xPM)*(2.*b1+(a21+a33)*G4Log(a33/a21))/(2*b3);
283 
284  G4double thetaPlus,thetaMinus,phiHalf; // final angular variables
285  nn = 0;
286  // t, psi, rho generation start (while angle < pi)
287  do {
288  //generate t by the rejection method
289  do {
290  ++nn;
291  t=G4UniformRand();
292  G4double a34=a33/(t*t);
293  G4double a22 = a34 + b1;
294  if(std::abs(b1)<0.0001*a34)
295  // special case of a34=a22 because of logarithm accuracy
296  {
297  f1=(1.-2.*xPM+4.*xPM*t*(1.-t))/(12.*a34*a34*a34*a34);
298  }
299  else
300  {
301  f1=-(1.-2.*xPM+4.*xPM*t*(1.-t))*(2.*b1+(a22+a34)*G4Log(a34/a22))/(2*b3);
302  }
303  if(f1<0.0 || f1> f1_max) // should never happend
304  {
305  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
306  << "outside allowed range f1=" << f1
307  << " is set to zero, a34 = "<< a34 << " a22 = "<<a22<<"."
308  << G4endl;
309  f1 = 0.0;
310  }
311  if(nn > nmax) { break; }
312  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
313  } while ( G4UniformRand()*f1_max > f1);
314  // generate psi by the rejection method
315  G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
316  // long version
317  G4double f2;
318  do {
319  ++nn;
320  psi=twopi*G4UniformRand();
321  f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+cos(2.*psi));
322  if(f2<0 || f2> f2_max) // should never happend
323  {
324  G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
325  << "outside allowed range f2=" << f2 << " is set to zero"
326  << G4endl;
327  f2 = 0.0;
328  }
329  if(nn >= nmax) { break; }
330  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
331  } while ( G4UniformRand()*f2_max > f2);
332 
333  // generate rho by direct transformation
334  G4double C2Term1=GammaMuonInv/(2.*xPM*t);
335  G4double C22 = C2Term1*C2Term1+C2Term2*C2Term2;
336  G4double C2=4.*C22*C22/sqrt(xPM);
337  G4double rhomax=(1./t-1.)*1.9/A027;
338  G4double beta=G4Log( (C2+rhomax*rhomax*rhomax*rhomax)/C2 );
339  rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4UniformRand())-1. ))*0.25);
340 
341  //now get from t and psi the kinematical variables
342  G4double u=sqrt(1./t-1.);
343  G4double xiHalf=0.5*rho*cos(psi);
344  phiHalf=0.5*rho/u*sin(psi);
345 
346  thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus;
347  thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus;
348 
349  // protection against infinite loop
350  if(nn > nmax) {
351  if(std::abs(thetaPlus)>pi) { thetaPlus = 0.0; }
352  if(std::abs(thetaMinus)>pi) { thetaMinus = 0.0; }
353  }
354 
355  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
356  } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi);
357 
358  // now construct the vectors
359  // azimuthal symmetry, take phi0 at random between 0 and 2 pi
360  G4double phi0=twopi*G4UniformRand();
361  G4double EPlus=xPlus*Egam;
362  G4double EMinus=xMinus*Egam;
363 
364  // mu+ mu- directions for gamma in z-direction
365  G4ThreeVector MuPlusDirection ( sin(thetaPlus) *cos(phi0+phiHalf),
366  sin(thetaPlus) *sin(phi0+phiHalf), cos(thetaPlus) );
367  G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf),
368  -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) );
369  // rotate to actual gamma direction
370  MuPlusDirection.rotateUz(GammaDirection);
371  MuMinusDirection.rotateUz(GammaDirection);
373  // create G4DynamicParticle object for the particle1
374  G4DynamicParticle* aParticle1= new G4DynamicParticle(
375  G4MuonPlus::MuonPlus(),MuPlusDirection,EPlus-Mmuon);
376  aParticleChange.AddSecondary(aParticle1);
377  // create G4DynamicParticle object for the particle2
378  G4DynamicParticle* aParticle2= new G4DynamicParticle(
379  G4MuonMinus::MuonMinus(),MuMinusDirection,EMinus-Mmuon);
380  aParticleChange.AddSecondary(aParticle2);
381  //
382  // Kill the incident photon
383  //
387  // Reset NbOfInteractionLengthLeft and return aParticleChange
388  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
389 }
390 
391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
392 
394  const G4DynamicParticle* aDynamicGamma,
395  const G4Material* aMaterial)
396 {
397  // select randomly 1 element within the material, invoked by PostStepDoIt
398 
399  const G4int NumberOfElements = aMaterial->GetNumberOfElements();
400  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
401  const G4Element* elm = (*theElementVector)[0];
402 
403  if (NumberOfElements > 1) {
404  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
405 
406  G4double PartialSumSigma = 0.;
408 
409  for (G4int i=0; i<NumberOfElements; ++i)
410  {
411  elm = (*theElementVector)[i];
412  PartialSumSigma += NbOfAtomsPerVolume[i]
413  *GetCrossSectionPerAtom(aDynamicGamma, elm);
414  if (rval <= PartialSumSigma) { break; }
415  }
416  }
417  return elm;
418 }
419 
420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
421 
423 {
424  G4String comments ="gamma->mu+mu- Bethe Heitler process, SubType= ";
425  G4cout << G4endl << GetProcessName() << ": " << comments
426  << GetProcessSubType() << G4endl;
427  G4cout << " good cross section parametrization from "
428  << G4BestUnit(LowestEnergyLimit,"Energy")
429  << " to " << HighestEnergyLimit/GeV << " GeV for all Z." << G4endl;
430 }
431 
432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void BuildPhysicsTable(const G4ParticleDefinition &) override
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double ComputeCrossSectionPerAtom(G4double GammaEnergy, G4int Z)
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const G4Element * SelectRandomAtom(const G4DynamicParticle *aDynamicGamma, const G4Material *aMaterial)
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:100
Float_t f2
Double_t beta
G4double G4Log(G4double x)
Definition: G4Log.hh:230
void AddSecondary(G4Track *aSecondary)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4bool IsApplicable(const G4ParticleDefinition &) override
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4int nmax
Float_t Z
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4double GetCrossSectionPerAtom(const G4DynamicParticle *aDynamicGamma, const G4Element *anElement)
virtual void Initialize(const G4Track &)
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) override
G4ProcessType
static constexpr double electron_mass_c2
static constexpr double elm_coupling
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
const G4String & GetProcessName() const
Definition: G4VProcess.hh:411
static constexpr double eV
Definition: G4SIunits.hh:215
Float_t f1
Definition: G4Step.hh:76
void ProposeEnergy(G4double finalEnergy)
G4Material * GetMaterial() const
static const G4double sqrte
G4int GetZasInt() const
Definition: G4Element.hh:132
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
G4double G4ParticleHPJENDLHEData::G4double result
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
std::vector< G4Element * > G4ElementVector
int G4int
Definition: G4Types.hh:78
G4double GetA27(G4int Z) const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
static const G4double fac
static const G4double PowSat
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:99
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:435
G4ForceCondition
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4double GetZ13(G4double Z) const
static constexpr double pi
Definition: G4SIunits.hh:75
static constexpr double fine_structure_const
#define DBL_MAX
Definition: templates.hh:83
double B(double temperature)
void ProposeTrackStatus(G4TrackStatus status)
G4GammaConversionToMuons(const G4String &processName="GammaToMuPair", G4ProcessType type=fElectromagnetic)
static constexpr double GeV
Definition: G4SIunits.hh:217
const G4DynamicParticle * GetDynamicParticle() const
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
const double C2
static G4NistManager * Instance()
G4int GetProcessSubType() const
Definition: G4VProcess.hh:429
void SetNumberOfSecondaries(G4int totSecondaries)
G4double ComputeMeanFreePath(G4double GammaEnergy, const G4Material *aMaterial)