Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4BetheHeitlerModel.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 // $Id: G4BetheHeitlerModel.cc 110939 2018-06-27 12:02:21Z gunter $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4BetheHeitlerModel
34 //
35 // Author: Vladimir Ivanchenko on base of Michel Maire code
36 //
37 // Creation date: 15.03.2005
38 //
39 // Modifications:
40 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivantchenko)
41 // 24-06-05 Increase number of bins to 200 (V.Ivantchenko)
42 // 16-11-05 replace shootBit() by G4UniformRand() mma
43 // 04-12-05 SetProposedKineticEnergy(0.) for the killed photon (mma)
44 // 20-02-07 SelectRandomElement is called for any initial gamma energy
45 // in order to have selected element for polarized model (VI)
46 // 25-10-10 Removed unused table, added element selector (VI)
47 // 28-05-18 New version with improved screening function approximation, improved
48 // efficiency, documentation and cleanup. Corrected call to selecting
49 // target atom in the final state sampling. (M. Novak)
50 //
51 // Class Description:
52 //
53 // -------------------------------------------------------------------
54 //
55 
56 #include "G4BetheHeitlerModel.hh"
57 #include "G4PhysicalConstants.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "G4Electron.hh"
60 #include "G4Positron.hh"
61 #include "G4Gamma.hh"
62 #include "Randomize.hh"
64 #include "G4Pow.hh"
65 #include "G4Exp.hh"
66 #include "G4ModifiedTsai.hh"
67 
68 
70 std::vector<G4BetheHeitlerModel::ElementData*> G4BetheHeitlerModel::gElementData;
71 
72 
73 
75  const G4String& nam)
76 : G4VEmModel(nam),
77  fG4Calc(G4Pow::GetInstance()), fTheGamma(G4Gamma::Gamma()),
78  fTheElectron(G4Electron::Electron()), fThePositron(G4Positron::Positron()),
79  fParticleChange(nullptr)
80 {
82 }
83 
84 
85 
87 {
88  if (IsMaster()) {
89  // clear ElementData container
90  for (size_t iz = 0; iz < gElementData.size(); ++iz) {
91  if (gElementData[iz]) delete gElementData[iz];
92  }
93  gElementData.clear();
94  }
95 }
96 
97 
98 
100  const G4DataVector& cuts)
101 {
102  if (IsMaster()) {
104  }
106  if (IsMaster()) {
107  InitialiseElementSelectors(p, cuts);
108  }
109 }
110 
111 
112 
114  G4VEmModel* masterModel)
115 {
117 }
118 
119 
120 
121 // Calculates the microscopic cross section in GEANT4 internal units.
122 // A parametrized formula from L. Urban is used to estimate
123 // the total cross section.
124 // It gives a good description of the data from 1.5 MeV to 100 GeV.
125 // below 1.5 MeV: sigma=sigma(1.5MeV)*(GammaEnergy-2electronmass)
126 // *(GammaEnergy-2electronmass)
127 G4double
129  G4double gammaEnergy, G4double Z,
131 {
132  G4double xSection = 0.0 ;
133  // short versions
134  static const G4double kMC2 = CLHEP::electron_mass_c2;
135  // zero cross section below the kinematical limit: Eg<2mc^2
136  if (Z < 0.9 || gammaEnergy <= 2.0*kMC2) { return xSection; }
137  //
138  static const G4double gammaEnergyLimit = 1.5*CLHEP::MeV;
139  // set coefficients a, b c
140  static const G4double a0 = 8.7842e+2*CLHEP::microbarn;
141  static const G4double a1 = -1.9625e+3*CLHEP::microbarn;
142  static const G4double a2 = 1.2949e+3*CLHEP::microbarn;
143  static const G4double a3 = -2.0028e+2*CLHEP::microbarn;
144  static const G4double a4 = 1.2575e+1*CLHEP::microbarn;
145  static const G4double a5 = -2.8333e-1*CLHEP::microbarn;
146 
147  static const G4double b0 = -1.0342e+1*CLHEP::microbarn;
148  static const G4double b1 = 1.7692e+1*CLHEP::microbarn;
149  static const G4double b2 = -8.2381 *CLHEP::microbarn;
150  static const G4double b3 = 1.3063 *CLHEP::microbarn;
151  static const G4double b4 = -9.0815e-2*CLHEP::microbarn;
152  static const G4double b5 = 2.3586e-3*CLHEP::microbarn;
153 
154  static const G4double c0 = -4.5263e+2*CLHEP::microbarn;
155  static const G4double c1 = 1.1161e+3*CLHEP::microbarn;
156  static const G4double c2 = -8.6749e+2*CLHEP::microbarn;
157  static const G4double c3 = 2.1773e+2*CLHEP::microbarn;
158  static const G4double c4 = -2.0467e+1*CLHEP::microbarn;
159  static const G4double c5 = 6.5372e-1*CLHEP::microbarn;
160  // check low energy limit of the approximation (1.5 MeV)
161  G4double gammaEnergyOrg = gammaEnergy;
162  if (gammaEnergy < gammaEnergyLimit) { gammaEnergy = gammaEnergyLimit; }
163  // compute gamma energy variables
164  const G4double x = G4Log(gammaEnergy/kMC2);
165  const G4double x2 = x *x;
166  const G4double x3 = x2*x;
167  const G4double x4 = x3*x;
168  const G4double x5 = x4*x;
169  //
170  const G4double F1 = a0 + a1*x + a2*x2 + a3*x3 + a4*x4 + a5*x5;
171  const G4double F2 = b0 + b1*x + b2*x2 + b3*x3 + b4*x4 + b5*x5;
172  const G4double F3 = c0 + c1*x + c2*x2 + c3*x3 + c4*x4 + c5*x5;
173  // compute the approximated cross section
174  xSection = (Z + 1.)*(F1*Z + F2*Z*Z + F3);
175  // check if we are below the limit of the approximation and apply correction
176  if (gammaEnergyOrg < gammaEnergyLimit) {
177  const G4double dum = (gammaEnergyOrg-2.*kMC2)/(gammaEnergyLimit-2.*kMC2);
178  xSection *= dum*dum;
179  }
180  // make sure that the cross section is never negative
181  xSection = std::max(xSection, 0.);
182  return xSection;
183 }
184 
185 
186 
187 // The secondaries e+e- energies are sampled using the Bethe - Heitler
188 // cross sections with Coulomb correction.
189 // A modified version of the random number techniques of Butcher & Messel
190 // is used (Nuc Phys 20(1960),15).
191 //
192 // GEANT4 internal units.
193 //
194 // Note 1 : Effects due to the breakdown of the Born approximation at
195 // low energy are ignored.
196 // Note 2 : The differential cross section implicitly takes account of
197 // pair creation in both nuclear and atomic electron fields.
198 // However triplet prodution is not generated.
199 void G4BetheHeitlerModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
200  const G4MaterialCutsCouple* couple,
201  const G4DynamicParticle* aDynamicGamma,
203 {
204  // set some constant values
205  const G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
206  const G4double eps0 = CLHEP::electron_mass_c2/gammaEnergy;
207  //
208  // check kinematical limit: gamma energy(Eg) must be at least 2 e- rest mass
209  if (eps0 > 0.5) { return; }
210  //
211  // select target element of the material (probs. are based on partial x-secs)
212  const G4Element* anElement = SelectRandomAtom(couple,fTheGamma,gammaEnergy);
213  //
214  // get the random engine
215  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
216  //
217  // 'eps' is the total energy transferred to one of the e-/e+ pair in initial
218  // gamma energy units Eg. Since the corresponding DCS is symmetric on eps=0.5,
219  // the kinematical limits for eps0=mc^2/Eg <= eps <= 0.5
220  // 1. 'eps' is sampled uniformly on the [eps0, 0.5] inteval if Eg<Egsmall
221  // 2. otherwise, on the [eps_min, 0.5] interval according to the DCS (case 2.)
222  G4double eps;
223  // case 1.
224  static const G4double Egsmall = 2.*CLHEP::MeV;
225  if (gammaEnergy < Egsmall) {
226  eps = eps0 + (0.5-eps0)*rndmEngine->flat();
227  } else {
228  // case 2.
229  // get the Coulomb factor for the target element (Z) and gamma energy (Eg)
230  // F(Z) = 8*ln(Z)/3 if Eg <= 50 [MeV] => no Coulomb correction
231  // F(Z) = 8*ln(Z)/3 + 8*fc(Z) if Eg > 50 [MeV] => fc(Z) is the Coulomb cor.
232  //
233  // The screening variable 'delta(eps)' = 136*Z^{-1/3}*eps0/[eps(1-eps)]
234  // Due to the Coulomb correction, the DCS can go below zero even at
235  // kinematicaly allowed eps > eps0 values. In order to exclude this eps
236  // range with negative DCS, the minimum eps value will be set to eps_min =
237  // max[eps0, epsp] with epsp is the solution of SF(delta(epsp)) - F(Z)/2 = 0
238  // with SF being the screening function (SF1=SF2 at high value of delta).
239  // The solution is epsp = 0.5 - 0.5*sqrt[ 1 - 4*136*Z^{-1/3}eps0/deltap]
240  // with deltap = Exp[(42.038-F(Z))/8.29]-0.958. So the limits are:
241  // - when eps=eps_max = 0.5 => delta_min = 136*Z^{-1/3}*eps0/4
242  // - epsp = 0.5 - 0.5*sqrt[ 1 - delta_min/deltap]
243  // - and eps_min = max[eps0, epsp]
244  static const G4double midEnergy = 50.*CLHEP::MeV;
245  const G4int iZet = std::min(gMaxZet, anElement->GetZasInt());
246  const G4double deltaFactor = 136.*eps0/anElement->GetIonisation()->GetZ3();
247  G4double deltaMax = gElementData[iZet]->fDeltaMaxLow;
248  G4double FZ = 8.*anElement->GetIonisation()->GetlogZ3();
249  if (gammaEnergy > midEnergy) {
250  FZ += 8.*(anElement->GetfCoulomb());
251  deltaMax = gElementData[iZet]->fDeltaMaxHigh;
252  }
253  const G4double deltaMin = 4.*deltaFactor;
254  //
255  // compute the limits of eps
256  const G4double epsp = 0.5 - 0.5*std::sqrt(1. - deltaMin/deltaMax) ;
257  const G4double epsMin = std::max(eps0,epsp);
258  const G4double epsRange = 0.5 - epsMin;
259  //
260  // sample the energy rate (eps) of the created electron (or positron)
261  G4double F10, F20;
262  ScreenFunction12(deltaMin, F10, F20);
263  F10 -= FZ;
264  F20 -= FZ;
265  const G4double NormF1 = std::max(F10 * epsRange * epsRange, 0.);
266  const G4double NormF2 = std::max(1.5 * F20 , 0.);
267  const G4double NormCond = NormF1/(NormF1 + NormF2);
268  // we will need 3 uniform random number for each trial of sampling
269  G4double rndmv[3];
270  G4double greject = 0.;
271  do {
272  rndmEngine->flatArray(3, rndmv);
273  if (NormCond > rndmv[0]) {
274  eps = 0.5 - epsRange * fG4Calc->A13(rndmv[1]);
275  const G4double delta = deltaFactor/(eps*(1.-eps));
276  greject = (ScreenFunction1(delta)-FZ)/F10;
277  } else {
278  eps = epsMin + epsRange*rndmv[1];
279  const G4double delta = deltaFactor/(eps*(1.-eps));
280  greject = (ScreenFunction2(delta)-FZ)/F20;
281  }
282  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
283  } while (greject < rndmv[2]);
284  } // end of eps sampling
285  //
286  // select charges randomly
287  G4double eTotEnergy, pTotEnergy;
288  if (rndmEngine->flat() > 0.5) {
289  eTotEnergy = (1.-eps)*gammaEnergy;
290  pTotEnergy = eps*gammaEnergy;
291  } else {
292  pTotEnergy = (1.-eps)*gammaEnergy;
293  eTotEnergy = eps*gammaEnergy;
294  }
295  //
296  // sample pair kinematics
297  const G4double eKinEnergy = std::max(0.,eTotEnergy - CLHEP::electron_mass_c2);
298  const G4double pKinEnergy = std::max(0.,pTotEnergy - CLHEP::electron_mass_c2);
299  //
300  G4ThreeVector eDirection, pDirection;
301  //
303  eKinEnergy, pKinEnergy,
304  eDirection, pDirection);
305  // create G4DynamicParticle object for the particle1
306  G4DynamicParticle* aParticle1= new G4DynamicParticle(
307  fTheElectron,eDirection,eKinEnergy);
308  // create G4DynamicParticle object for the particle2
309  G4DynamicParticle* aParticle2= new G4DynamicParticle(
310  fThePositron,pDirection,pKinEnergy);
311  // Fill output vector
312  fvect->push_back(aParticle1);
313  fvect->push_back(aParticle2);
314  // kill incident photon
317 }
318 
319 
320 
321 // should be called only by the master and at initialisation
323 {
324  G4int size = gElementData.size();
325  if (size < gMaxZet+1) {
326  gElementData.resize(gMaxZet+1, nullptr);
327  }
328  // create for all elements that are in the detector
329  const G4ElementTable* elemTable = G4Element::GetElementTable();
330  size_t numElems = (*elemTable).size();
331  for (size_t ie = 0; ie < numElems; ++ie) {
332  const G4Element* elem = (*elemTable)[ie];
333  const G4int iz = std::min(gMaxZet, elem->GetZasInt());
334  if (!gElementData[iz]) { // create it if doesn't exist yet
335  G4double FZLow = 8.*elem->GetIonisation()->GetlogZ3();
336  G4double FZHigh = FZLow + 8.*elem->GetfCoulomb();
337  ElementData* elD = new ElementData();
338  elD->fDeltaMaxLow = G4Exp((42.038 - FZLow )/8.29) - 0.958;
339  elD->fDeltaMaxHigh = G4Exp((42.038 - FZHigh)/8.29) - 0.958;
340  gElementData[iz] = elD;
341  }
342  }
343 }
344 
Float_t x
Definition: compare.C:6
#define F10
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:792
virtual void SamplePairDirections(const G4DynamicParticle *dp, G4double elecKinEnergy, G4double posiKinEnergy, G4ThreeVector &dirElectron, G4ThreeVector &dirPositron, G4int Z=0, const G4Material *mat=nullptr)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static constexpr double microbarn
Definition: SystemOfUnits.h:87
G4double ScreenFunction2(const G4double delta)
G4double A13(G4double A) const
Definition: G4Pow.cc:138
G4double ScreenFunction1(const G4double delta)
void ScreenFunction12(const G4double delta, G4double &f1, G4double &f2)
G4ParticleDefinition * fTheGamma
const char * p
Definition: xmltok.h:285
static std::vector< ElementData * > gElementData
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual void flatArray(const int size, double *vect)=0
G4BetheHeitlerModel(const G4ParticleDefinition *p=0, const G4String &nam="BetheHeitler")
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:552
static const G4int gMaxZet
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:784
G4double G4Log(G4double x)
Definition: G4Log.hh:230
Float_t Z
double G4double
Definition: G4Types.hh:76
virtual double flat()=0
TCanvas * c2
Definition: plot_hist.C:75
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cut=0., G4double emax=DBL_MAX) override
static constexpr double electron_mass_c2
G4double GetZ3() const
std::vector< G4Element * > G4ElementTable
static constexpr double MeV
#define F20
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:582
G4int GetZasInt() const
Definition: G4Element.hh:132
std::map< G4String, G4AttDef > * GetInstance(const G4String &storeKey, G4bool &isNew)
const G4double a0
G4ParticleChangeForGamma * fParticleChange
int G4int
Definition: G4Types.hh:78
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
Definition: G4Pow.hh:56
G4IonisParamElm * GetIonisation() const
Definition: G4Element.hh:199
G4double GetKineticEnergy() const
G4ParticleDefinition * fThePositron
G4ParticleDefinition * fTheElectron
G4bool IsMaster() const
Definition: G4VEmModel.hh:700
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:575
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:131
static const G4double eps
G4double GetfCoulomb() const
Definition: G4Element.hh:191
G4double GetlogZ3() const
Float_t x2[n_points_geant4]
Definition: compare.C:26
void ProposeTrackStatus(G4TrackStatus status)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398