Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PairProductionRelModel Class Reference

#include <G4PairProductionRelModel.hh>

Inheritance diagram for G4PairProductionRelModel:
Collaboration diagram for G4PairProductionRelModel:

Public Member Functions

 G4PairProductionRelModel (const G4ParticleDefinition *p=0, const G4String &nam="BetheHeitlerLPM")
 
virtual ~G4PairProductionRelModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel) override
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cut=0., G4double emax=DBL_MAX) override
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double) override
 
void SetCurrentElement (G4double Z)
 
void SetLPMconstant (G4double val)
 
G4double LPMconstant () const
 
void SetLPMflag (G4bool)
 
G4bool LPMflag () const
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector
< G4EmElementSelector * > * 
GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
virtual G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy)
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Member Functions

G4double Phi1 (G4double delta) const
 
G4double Phi2 (G4double delta) const
 
G4double ScreenFunction1 (G4double ScreenVariable)
 
G4double ScreenFunction2 (G4double ScreenVariable)
 
G4double DeltaMax () const
 
G4double DeltaMin (G4double) const
 
void CalcLPMFunctions (G4double k, G4double eplus)
 
G4double ComputeXSectionPerAtom (G4double totalEnergy, G4double Z)
 
G4double ComputeDXSectionPerAtom (G4double eplusEnergy, G4double totalEnergy, G4double Z)
 
G4double ComputeRelDXSectionPerAtom (G4double eplusEnergy, G4double totalEnergy, G4double Z)
 
G4PairProductionRelModeloperator= (const G4PairProductionRelModel &right)=delete
 
 G4PairProductionRelModel (const G4PairProductionRelModel &)=delete
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 

Protected Attributes

G4NistManagernist
 
G4ParticleDefinitiontheGamma
 
G4ParticleDefinitiontheElectron
 
G4ParticleDefinitionthePositron
 
G4ParticleChangeForGammafParticleChange
 
G4double fLPMconstant
 
G4bool fLPMflag
 
G4double z13
 
G4double z23
 
G4double lnZ
 
G4double Fel
 
G4double Finel
 
G4double fCoulomb
 
G4double currentZ
 
G4double lpmEnergy
 
G4double xiLPM
 
G4double phiLPM
 
G4double gLPM
 
G4bool use_completescreening
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 

Static Protected Attributes

static const G4double xgi [8]
 
static const G4double wgi [8]
 
static const G4double Fel_light [5] = {0., 5.31 , 4.79 , 4.74 , 4.71}
 
static const G4double Finel_light [5] = {0., 6.144 , 5.621 , 5.805 , 5.924}
 
static const G4double facFel = G4Log(184.15)
 
static const G4double facFinel = G4Log(1194.)
 
static const G4double preS1 = 1./(184.15*184.15)
 
static const G4double logTwo = G4Log(2.)
 
static const G4double xsfactor
 
static const G4double Egsmall =2.*MeV
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 61 of file G4PairProductionRelModel.hh.

Constructor & Destructor Documentation

G4PairProductionRelModel::G4PairProductionRelModel ( const G4ParticleDefinition p = 0,
const G4String nam = "BetheHeitlerLPM" 
)
explicit

Definition at line 90 of file G4PairProductionRelModel.cc.

92  : G4VEmModel(nam),
94  fLPMflag(true),
95  lpmEnergy(0.),
97 {
98  fParticleChange = nullptr;
102 
104 
105  currentZ = z13 = z23 = lnZ = Fel = Finel = fCoulomb = phiLPM = gLPM = xiLPM = 0;
106 
107 }
G4ParticleDefinition * thePositron
static G4NistManager * Instance()
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
float electron_mass_c2
Definition: hepunit.py:274
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4ParticleChangeForGamma * fParticleChange
static G4Electron * Electron()
Definition: G4Electron.cc:94
static constexpr double pi
Definition: G4SIunits.hh:75
G4ParticleDefinition * theElectron

Here is the call graph for this function:

G4PairProductionRelModel::~G4PairProductionRelModel ( )
virtual

Definition at line 111 of file G4PairProductionRelModel.cc.

112 {}
G4PairProductionRelModel::G4PairProductionRelModel ( const G4PairProductionRelModel )
protecteddelete

Member Function Documentation

void G4PairProductionRelModel::CalcLPMFunctions ( G4double  k,
G4double  eplus 
)
protected

Definition at line 238 of file G4PairProductionRelModel.cc.

239 {
240  // *** calculate lpm variable s & sprime ***
241  // Klein eqs. (78) & (79)
242  G4double sprime = sqrt(0.125*k*lpmEnergy/(eplusEnergy*(k-eplusEnergy)));
243 
244  G4double s1 = preS1*z23;
245  G4double logS1 = 2./3.*lnZ-2.*facFel;
246  G4double logTS1 = logTwo+logS1;
247 
248  xiLPM = 2.;
249 
250  if (sprime>1)
251  xiLPM = 1.;
252  else if (sprime>sqrt(2.)*s1) {
253  G4double h = G4Log(sprime)/logTS1;
254  xiLPM = 1+h-0.08*(1-h)*(1-sqr(1-h))/logTS1;
255  }
256 
257  G4double s0 = sprime/sqrt(xiLPM);
258  // G4cout<<"k="<<k<<" y="<<eplusEnergy/k<<G4endl;
259  // G4cout<<"s0="<<s0<<G4endl;
260 
261  // *** calculate supression functions phi and G ***
262  // Klein eqs. (77)
263  G4double s2=s0*s0;
264  G4double s3=s0*s2;
265  G4double s4=s2*s2;
266 
267  if (s0<0.1) {
268  // high suppression limit
269  phiLPM = 6.*s0 - 18.84955592153876*s2 + 39.47841760435743*s3
270  - 57.69873135166053*s4;
271  gLPM = 37.69911184307752*s2 - 236.8705056261446*s3 + 807.7822389*s4;
272  }
273  else if (s0<1.9516) {
274  // intermediate suppression
275  // using eq.77 approxim. valid s0<2.
276  phiLPM = 1.-G4Exp(-6.*s0*(1.+(3.-pi)*s0)
277  +s3/(0.623+0.795*s0+0.658*s2));
278  if (s0<0.415827397755) {
279  // using eq.77 approxim. valid 0.07<s<2
280  G4double psiLPM = 1-G4Exp(-4*s0-8*s2/(1+3.936*s0+4.97*s2-0.05*s3+7.50*s4));
281  gLPM = 3*psiLPM-2*phiLPM;
282  }
283  else {
284  // using alternative parametrisiation
285  G4double pre = -0.16072300849123999 + s0*3.7550300067531581 + s2*-1.7981383069010097
286  + s3*0.67282686077812381 + s4*-0.1207722909879257;
287  gLPM = tanh(pre);
288  }
289  }
290  else {
291  // low suppression limit valid s>2.
292  phiLPM = 1. - 0.0119048/s4;
293  gLPM = 1. - 0.0230655/s4;
294  }
295 
296  // *** make sure suppression is smaller than 1 ***
297  // *** caused by Migdal approximation in xi ***
298  if (xiLPM*phiLPM>1. || s0>0.57) xiLPM=1./phiLPM;
299 }
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static constexpr double pi
Definition: G4SIunits.hh:75
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PairProductionRelModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition ,
G4double  kinEnergy,
G4double  Z,
G4double  A = 0.,
G4double  cut = 0.,
G4double  emax = DBL_MAX 
)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 305 of file G4PairProductionRelModel.cc.

307 {
308  G4double crossSection = 0.0 ;
309  if ( gammaEnergy <= 2.0*electron_mass_c2 ) return crossSection;
310 
312  // choose calculator according to parameters and switches
313  // in the moment only one calculator:
314  crossSection=ComputeXSectionPerAtom(gammaEnergy,Z);
315 
316  G4double xi = Finel/(Fel - fCoulomb); // inelastic contribution
317  crossSection *= xsfactor*Z*(Z+xi);
318 
319  return crossSection;
320 }
G4double ComputeXSectionPerAtom(G4double totalEnergy, G4double Z)
float electron_mass_c2
Definition: hepunit.py:274
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4PairProductionRelModel::ComputeDXSectionPerAtom ( G4double  eplusEnergy,
G4double  totalEnergy,
G4double  Z 
)
protected

Definition at line 181 of file G4PairProductionRelModel.cc.

184 {
185  // most simple case - complete screening:
186 
187  // dsig/dE+ = 4 * alpha * Z**2 * r0**2 / k
188  // * [ (y**2 + (1-y**2) + 2/3*y*(1-y) ) * ( log (183 * Z**-1/3) + 1/9 * y*(1-y) ]
189  // y = E+/k
190  G4double yp=eplusEnergy/totalEnergy;
191  G4double ym=1.-yp;
192 
193  G4double cross = 0.;
195  cross = (yp*yp + ym*ym + 2./3.*ym*yp)*(Fel - fCoulomb) + yp*ym/9.;
196  else {
197  G4double delta = 0.25*DeltaMin(totalEnergy)/(yp*ym);
198  cross = (yp*yp + ym*ym)*(0.25*Phi1(delta) - lnZ/3. - fCoulomb)
199  + 2./3.*ym*yp*(0.25*Phi2(delta) - lnZ/3. - fCoulomb);
200  }
201  return cross/totalEnergy;
202 
203 }
G4double Phi1(G4double delta) const
G4double Phi2(G4double delta) const
G4double DeltaMin(G4double) const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PairProductionRelModel::ComputeRelDXSectionPerAtom ( G4double  eplusEnergy,
G4double  totalEnergy,
G4double  Z 
)
protected

Definition at line 207 of file G4PairProductionRelModel.cc.

210 {
211  // most simple case - complete screening:
212 
213  // dsig/dE+ = 4 * alpha * Z**2 * r0**2 / k
214  // * [ (y**2 + (1-y**2) + 2/3*y*(1-y) ) * ( log (183 * Z**-1/3) + 1/9 * y*(1-y) ]
215  // y = E+/k
216  G4double yp=eplusEnergy/totalEnergy;
217  G4double ym=1.-yp;
218 
219  CalcLPMFunctions(totalEnergy,eplusEnergy); // gamma
220 
221  G4double cross = 0.;
223  cross = xiLPM*(2./3.*phiLPM*(yp*yp + ym*ym) + gLPM)*(Fel - fCoulomb);
224  else {
225  G4double delta = 0.25*DeltaMin(totalEnergy)/(yp*ym);
226  cross = (1./3.*gLPM + 2./3.*phiLPM)*(yp*yp + ym*ym)
227  *(0.25*Phi1(delta) - lnZ/3. - fCoulomb)
228  + 2./3.*gLPM*ym*yp*(0.25*Phi2(delta) - lnZ/3. - fCoulomb);
229  cross *= xiLPM;
230  }
231  return cross/totalEnergy;
232 
233 }
G4double Phi1(G4double delta) const
G4double Phi2(G4double delta) const
G4double DeltaMin(G4double) const
void CalcLPMFunctions(G4double k, G4double eplus)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PairProductionRelModel::ComputeXSectionPerAtom ( G4double  totalEnergy,
G4double  Z 
)
protected

Definition at line 137 of file G4PairProductionRelModel.cc.

138 {
139  G4double cross = 0.0;
140 
141  // number of intervals and integration step
142  G4double vcut = electron_mass_c2/totalEnergy ;
143 
144  // limits by the screening variable
145  G4double dmax = DeltaMax();
146  G4double dmin = min(DeltaMin(totalEnergy),dmax);
147  G4double vcut1 = 0.5 - 0.5*sqrt(1. - dmin/dmax) ;
148  vcut = max(vcut, vcut1);
149 
150 
151  G4double vmax = 0.5;
152  G4int n = 1; // needs optimisation
153 
154  G4double delta = (vmax - vcut)*totalEnergy/G4double(n);
155 
156  G4double e0 = vcut*totalEnergy;
157  G4double xs;
158 
159  // simple integration
160  for(G4int l=0; l<n; l++,e0 += delta) {
161  for(G4int i=0; i<8; i++) {
162 
163  G4double eg = (e0 + xgi[i]*delta);
164  if (fLPMflag && totalEnergy>100.*GeV)
165  xs = ComputeRelDXSectionPerAtom(eg,totalEnergy,Z);
166  else
167  xs = ComputeDXSectionPerAtom(eg,totalEnergy,Z);
168  cross += wgi[i]*xs;
169 
170  }
171  }
172 
173  cross *= delta*2.;
174 
175  return cross;
176 }
int G4int
Definition: G4Types.hh:78
static const G4double xgi[8]
static const G4double wgi[8]
G4double DeltaMin(G4double) const
float electron_mass_c2
Definition: hepunit.py:274
const G4int n
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
G4double ComputeRelDXSectionPerAtom(G4double eplusEnergy, G4double totalEnergy, G4double Z)
double G4double
Definition: G4Types.hh:76
G4double ComputeDXSectionPerAtom(G4double eplusEnergy, G4double totalEnergy, G4double Z)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PairProductionRelModel::DeltaMax ( ) const
inlineprotected

Definition at line 279 of file G4PairProductionRelModel.hh.

280 {
281  // k > 50 MeV
282  G4double FZ = 8.*(lnZ/3. + fCoulomb);
283  return std::exp( (42.24-FZ)/8.368 ) + 0.952;
284 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4PairProductionRelModel::DeltaMin ( G4double  k) const
inlineprotected

Definition at line 286 of file G4PairProductionRelModel.hh.

287 {
288  return 4.*136./z13*(CLHEP::electron_mass_c2/k);
289 }
static constexpr double electron_mass_c2

Here is the caller graph for this function:

void G4PairProductionRelModel::Initialise ( const G4ParticleDefinition p,
const G4DataVector cuts 
)
overridevirtual

Implements G4VEmModel.

Definition at line 116 of file G4PairProductionRelModel.cc.

118 {
120  if(IsMaster() && LowEnergyLimit() < HighEnergyLimit()) {
121  InitialiseElementSelectors(p, cuts);
122  }
123 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
G4ParticleChangeForGamma * fParticleChange
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132

Here is the call graph for this function:

void G4PairProductionRelModel::InitialiseLocal ( const G4ParticleDefinition ,
G4VEmModel masterModel 
)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 127 of file G4PairProductionRelModel.cc.

129 {
130  if(LowEnergyLimit() < HighEnergyLimit()) {
132  }
133 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:809
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:817

Here is the call graph for this function:

G4double G4PairProductionRelModel::LPMconstant ( ) const
inline

Definition at line 171 of file G4PairProductionRelModel.hh.

172 {
173  return fLPMconstant;
174 }
G4bool G4PairProductionRelModel::LPMflag ( ) const
inline

Definition at line 187 of file G4PairProductionRelModel.hh.

188 {
189  return fLPMflag;
190 }
G4PairProductionRelModel& G4PairProductionRelModel::operator= ( const G4PairProductionRelModel right)
protecteddelete
G4double G4PairProductionRelModel::Phi1 ( G4double  delta) const
inlineprotected

Definition at line 218 of file G4PairProductionRelModel.hh.

219 {
220  G4double screenVal;
221 
222  if (delta > 1.)
223  screenVal = 21.12 - 4.184*std::log(delta+0.952);
224  else
225  screenVal = 20.868 - delta*(3.242 - 0.625*delta);
226 
227  return screenVal;
228 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4PairProductionRelModel::Phi2 ( G4double  delta) const
inlineprotected

Definition at line 232 of file G4PairProductionRelModel.hh.

233 {
234  G4double screenVal;
235 
236  if (delta > 1.)
237  screenVal = 21.12 - 4.184*std::log(delta+0.952);
238  else
239  screenVal = 20.209 - delta*(1.930 + 0.086*delta);
240 
241  return screenVal;
242 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

void G4PairProductionRelModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  fvect,
const G4MaterialCutsCouple couple,
const G4DynamicParticle aDynamicGamma,
G4double  tmin,
G4double  maxEnergy 
)
overridevirtual

Implements G4VEmModel.

Definition at line 325 of file G4PairProductionRelModel.cc.

342 {
343  const G4Material* aMaterial = couple->GetMaterial();
344 
345  G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
346  G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
347 
348  G4double epsil ;
349  G4double epsil0 = electron_mass_c2/GammaEnergy ;
350  if(epsil0 > 1.0) { return; }
351 
352  // do it fast if GammaEnergy < 2. MeV
353  SetupForMaterial(theGamma, aMaterial, GammaEnergy);
354 
355  // select randomly one element constituing the material
356  const G4Element* anElement =
357  SelectRandomAtom(aMaterial, theGamma, GammaEnergy);
358 
359  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
360 
361  if (GammaEnergy < Egsmall) {
362 
363  epsil = epsil0 + (0.5-epsil0)*rndmEngine->flat();
364 
365  } else {
366  // now comes the case with GammaEnergy >= 2. MeV
367 
368  // Extract Coulomb factor for this Element
369  G4double FZ = 8.*(anElement->GetIonisation()->GetlogZ3());
370  if (GammaEnergy > 50.*MeV) { FZ += 8.*(anElement->GetfCoulomb()); }
371 
372  // limits of the screening variable
373  G4double screenfac = 136.*epsil0/(anElement->GetIonisation()->GetZ3());
374  G4double screenmax = G4Exp ((42.24 - FZ)/8.368) - 0.952 ;
375  G4double screenmin = min(4.*screenfac,screenmax);
376 
377  // limits of the energy sampling
378  G4double epsil1 = 0.5 - 0.5*sqrt(1. - screenmin/screenmax) ;
379  G4double epsilmin = max(epsil0,epsil1) , epsilrange = 0.5 - epsilmin;
380 
381  //
382  // sample the energy rate of the created electron (or positron)
383  //
384  //G4double epsil, screenvar, greject ;
385  G4double screenvar, greject ;
386 
387  G4double F10 = ScreenFunction1(screenmin) - FZ;
388  G4double F20 = ScreenFunction2(screenmin) - FZ;
389  G4double NormF1 = max(F10*epsilrange*epsilrange,0.);
390  G4double NormF2 = max(1.5*F20,0.);
391 
392  do {
393  if ( NormF1/(NormF1+NormF2) > rndmEngine->flat() ) {
394  epsil = 0.5 - epsilrange*nist->GetZ13(rndmEngine->flat());
395  screenvar = screenfac/(epsil*(1-epsil));
396  if (fLPMflag && GammaEnergy>100.*GeV) {
397  CalcLPMFunctions(GammaEnergy,GammaEnergy*epsil);
398  greject = xiLPM*((gLPM+2.*phiLPM)*Phi1(screenvar) -
399  gLPM*Phi2(screenvar) - phiLPM*FZ)/F10;
400  }
401  else {
402  greject = (ScreenFunction1(screenvar) - FZ)/F10;
403  }
404 
405  } else {
406  epsil = epsilmin + epsilrange*rndmEngine->flat();
407  screenvar = screenfac/(epsil*(1-epsil));
408  if (fLPMflag && GammaEnergy>100.*GeV) {
409  CalcLPMFunctions(GammaEnergy,GammaEnergy*epsil);
410  greject = xiLPM*((0.5*gLPM+phiLPM)*Phi1(screenvar) +
411  0.5*gLPM*Phi2(screenvar) - 0.5*(gLPM+phiLPM)*FZ)/F20;
412  }
413  else {
414  greject = (ScreenFunction2(screenvar) - FZ)/F20;
415  }
416  }
417 
418  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
419  } while( greject < rndmEngine->flat());
420 
421  } // end of epsil sampling
422 
423  //
424  // fixe charges randomly
425  //
426 
427  G4double ElectTotEnergy, PositTotEnergy;
428  if (rndmEngine->flat() > 0.5) {
429 
430  ElectTotEnergy = (1.-epsil)*GammaEnergy;
431  PositTotEnergy = epsil*GammaEnergy;
432 
433  } else {
434 
435  PositTotEnergy = (1.-epsil)*GammaEnergy;
436  ElectTotEnergy = epsil*GammaEnergy;
437  }
438 
439  //
440  // scattered electron (positron) angles. ( Z - axis along the parent photon)
441  //
442  // universal distribution suggested by L. Urban
443  // (Geant3 manual (1993) Phys211),
444  // derived from Tsai distribution (Rev Mod Phys 49,421(1977))
445 
446  G4double u = -G4Log(rndmEngine->flat()*rndmEngine->flat());
447 
448  if (9. > 36.*rndmEngine->flat()) { u *= 1.6; }
449  else { u *= 0.53333; }
450 
451  G4double TetEl = u*electron_mass_c2/ElectTotEnergy;
452  G4double TetPo = u*electron_mass_c2/PositTotEnergy;
453  G4double Phi = twopi * rndmEngine->flat();
454  G4double dxEl= sin(TetEl)*cos(Phi),dyEl= sin(TetEl)*sin(Phi),dzEl=cos(TetEl);
455  G4double dxPo=-sin(TetPo)*cos(Phi),dyPo=-sin(TetPo)*sin(Phi),dzPo=cos(TetPo);
456 
457  //
458  // kinematic of the created pair
459  //
460  // the electron and positron are assumed to have a symetric
461  // angular distribution with respect to the Z axis along the parent photon.
462 
463  G4double ElectKineEnergy = max(0.,ElectTotEnergy - electron_mass_c2);
464 
465  G4ThreeVector ElectDirection (dxEl, dyEl, dzEl);
466  ElectDirection.rotateUz(GammaDirection);
467 
468  // create G4DynamicParticle object for the particle1
469  G4DynamicParticle* aParticle1= new G4DynamicParticle(
470  theElectron,ElectDirection,ElectKineEnergy);
471 
472  // the e+ is always created (even with Ekine=0) for further annihilation.
473 
474  G4double PositKineEnergy = max(0.,PositTotEnergy - electron_mass_c2);
475 
476  G4ThreeVector PositDirection (dxPo, dyPo, dzPo);
477  PositDirection.rotateUz(GammaDirection);
478 
479  // create G4DynamicParticle object for the particle2
480  G4DynamicParticle* aParticle2= new G4DynamicParticle(
481  thePositron,PositDirection,PositKineEnergy);
482 
483  // Fill output vector
484  fvect->push_back(aParticle1);
485  fvect->push_back(aParticle2);
486 
487  // kill incident photon
490 }
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double) override
G4double Phi1(G4double delta) const
G4double GetKineticEnergy() const
G4double Phi2(G4double delta) const
G4double GetfCoulomb() const
Definition: G4Element.hh:191
virtual double flat()=0
#define F20
G4ParticleDefinition * thePositron
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double ScreenFunction1(G4double ScreenVariable)
G4double GetZ13(G4double Z) const
const G4ThreeVector & GetMomentumDirection() const
void CalcLPMFunctions(G4double k, G4double eplus)
float electron_mass_c2
Definition: hepunit.py:274
G4double ScreenFunction2(G4double ScreenVariable)
#define F10
double flat()
Definition: G4AblaRandom.cc:47
G4double GetlogZ3() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4ParticleChangeForGamma * fParticleChange
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4IonisParamElm * GetIonisation() const
Definition: G4Element.hh:199
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static constexpr double MeV
Definition: G4SIunits.hh:214
G4ParticleDefinition * theElectron
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4double GetZ3() const
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
const G4Material * GetMaterial() const

Here is the call graph for this function:

G4double G4PairProductionRelModel::ScreenFunction1 ( G4double  ScreenVariable)
inlineprotected

Definition at line 244 of file G4PairProductionRelModel.hh.

248 {
249  G4double screenVal;
250 
251  if (ScreenVariable > 1.)
252  screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
253  else
254  screenVal = 42.392 - ScreenVariable*(7.796 - 1.961*ScreenVariable);
255 
256  return screenVal;
257 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4PairProductionRelModel::ScreenFunction2 ( G4double  ScreenVariable)
inlineprotected

Definition at line 261 of file G4PairProductionRelModel.hh.

265 {
266  G4double screenVal;
267 
268  if (ScreenVariable > 1.)
269  screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
270  else
271  screenVal = 41.405 - ScreenVariable*(5.828 - 0.8945*ScreenVariable);
272 
273  return screenVal;
274 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

void G4PairProductionRelModel::SetCurrentElement ( G4double  Z)
inline

Definition at line 194 of file G4PairProductionRelModel.hh.

195 {
196  if(Z != currentZ) {
197  currentZ = Z;
198 
199  G4int iz = G4lrint(Z);
200  z13 = nist->GetZ13(iz);
201  z23 = z13*z13;
202  lnZ = nist->GetLOGZ(iz);
203 
204  if (iz <= 4) {
205  Fel = Fel_light[iz];
206  Finel = Finel_light[iz] ;
207  }
208  else {
209  Fel = facFel - lnZ/3. ;
210  Finel = facFinel - 2.*lnZ/3. ;
211  }
213  }
214 }
G4double GetfCoulomb() const
Definition: G4Element.hh:191
static const G4double Finel_light[5]
int G4int
Definition: G4Types.hh:78
G4double GetZ13(G4double Z) const
G4double GetLOGZ(G4int Z) const
int G4lrint(double ad)
Definition: templates.hh:163
static const G4double Fel_light[5]
const G4Element * GetCurrentElement() const
Definition: G4VEmModel.hh:469

Here is the call graph for this function:

Here is the caller graph for this function:

void G4PairProductionRelModel::SetLPMconstant ( G4double  val)
inline

Definition at line 163 of file G4PairProductionRelModel.hh.

164 {
165  fLPMconstant = val;
166 }
void G4PairProductionRelModel::SetLPMflag ( G4bool  val)
inline

Definition at line 179 of file G4PairProductionRelModel.hh.

180 {
181  fLPMflag = val;
182 }
void G4PairProductionRelModel::SetupForMaterial ( const G4ParticleDefinition ,
const G4Material mat,
G4double   
)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 495 of file G4PairProductionRelModel.cc.

497 {
499  // G4cout<<" lpmEnergy="<<lpmEnergy<<G4endl;
500 }
G4double GetRadlen() const
Definition: G4Material.hh:220

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

G4double G4PairProductionRelModel::currentZ
protected

Definition at line 140 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::Egsmall =2.*MeV
staticprotected

Definition at line 155 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::facFel = G4Log(184.15)
staticprotected

Definition at line 152 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::facFinel = G4Log(1194.)
staticprotected

Definition at line 153 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::fCoulomb
protected

Definition at line 139 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::Fel
protected

Definition at line 139 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::Fel_light = {0., 5.31 , 4.79 , 4.74 , 4.71}
staticprotected

Definition at line 150 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::Finel
protected

Definition at line 139 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::Finel_light = {0., 6.144 , 5.621 , 5.805 , 5.924}
staticprotected

Definition at line 151 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::fLPMconstant
protected

Definition at line 134 of file G4PairProductionRelModel.hh.

G4bool G4PairProductionRelModel::fLPMflag
protected

Definition at line 135 of file G4PairProductionRelModel.hh.

G4ParticleChangeForGamma* G4PairProductionRelModel::fParticleChange
protected

Definition at line 132 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::gLPM
protected

Definition at line 144 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::lnZ
protected

Definition at line 138 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::logTwo = G4Log(2.)
staticprotected

Definition at line 155 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::lpmEnergy
protected

Definition at line 143 of file G4PairProductionRelModel.hh.

G4NistManager* G4PairProductionRelModel::nist
protected

Definition at line 127 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::phiLPM
protected

Definition at line 144 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::preS1 = 1./(184.15*184.15)
staticprotected

Definition at line 155 of file G4PairProductionRelModel.hh.

G4ParticleDefinition* G4PairProductionRelModel::theElectron
protected

Definition at line 130 of file G4PairProductionRelModel.hh.

G4ParticleDefinition* G4PairProductionRelModel::theGamma
protected

Definition at line 129 of file G4PairProductionRelModel.hh.

G4ParticleDefinition* G4PairProductionRelModel::thePositron
protected

Definition at line 131 of file G4PairProductionRelModel.hh.

G4bool G4PairProductionRelModel::use_completescreening
protected

Definition at line 147 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::wgi
staticprotected
Initial value:
={ 0.0506, 0.1112, 0.1569, 0.1813,
0.1813, 0.1569, 0.1112, 0.0506 }

Definition at line 149 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::xgi
staticprotected
Initial value:
={ 0.0199, 0.1017, 0.2372, 0.4083,
0.5917, 0.7628, 0.8983, 0.9801 }

Definition at line 149 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::xiLPM
protected

Definition at line 144 of file G4PairProductionRelModel.hh.

const G4double G4PairProductionRelModel::xsfactor
staticprotected
G4double G4PairProductionRelModel::z13
protected

Definition at line 138 of file G4PairProductionRelModel.hh.

G4double G4PairProductionRelModel::z23
protected

Definition at line 138 of file G4PairProductionRelModel.hh.


The documentation for this class was generated from the following files: