Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes
G4MuPairProductionModel Class Reference

#include <G4MuPairProductionModel.hh>

Inheritance diagram for G4MuPairProductionModel:
G4VEmModel G4hPairProductionModel

Public Member Functions

 G4MuPairProductionModel (const G4ParticleDefinition *p=0, const G4String &nam="muPairProd")
 
virtual ~G4MuPairProductionModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double)
 
void SetLowestKineticEnergy (G4double e)
 
void SetParticle (const G4ParticleDefinition *)
 
- 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 CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, 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 MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector
< G4EmElementSelector * > * 		GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
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 G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
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 SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4Element * GetCurrentElement () const
 

Protected Member Functions

G4double ComputMuPairLoss (G4double Z, G4double tkin, G4double cut, G4double tmax)
 
G4double ComputeMicroscopicCrossSection (G4double tkin, G4double Z, G4double cut)
 
virtual G4double ComputeDMicroscopicCrossSection (G4double tkin, G4double Z, G4double pairEnergy)
 
G4double MaxSecondaryEnergyForElement (G4double kineticEnergy, G4double Z)
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 

Protected Attributes

const G4ParticleDefinition * particle
 
G4NistManager * nist
 
G4double factorForCross
 
G4double sqrte
 
G4double particleMass
 
G4double z13
 
G4double z23
 
G4double lnZ
 
G4int currentZ
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 

Static Protected Attributes

static const G4double xgi [8]
 
static const G4double wgi [8]
 

Detailed Description

Definition at line 73 of file G4MuPairProductionModel.hh.

Constructor & Destructor Documentation

G4MuPairProductionModel::G4MuPairProductionModel ( const G4ParticleDefinition *  p = 0,
const G4String &  nam = "muPairProd" 
)

Definition at line 104 of file G4MuPairProductionModel.cc.

References G4Electron::Electron(), G4ParticleDefinition::GetPDGMass(), G4NistManager::Instance(), lnZ, nist, particleMass, G4Positron::Positron(), SetParticle(), python.hepunit::TeV, z13, and z23.

106  : G4VEmModel(nam),
107  particle(0),
108  factorForCross(4.*fine_structure_const*fine_structure_const
109  *classic_electr_radius*classic_electr_radius/(3.*pi)),
110  sqrte(sqrt(G4Exp(1.))),
111  currentZ(0),
112  fParticleChange(0),
113  minPairEnergy(4.*electron_mass_c2),
114  lowestKinEnergy(1.0*GeV),
115  nzdat(5),
116  nYBinPerDecade(4),
117  nbiny(1000),
118  nbine(0),
119  ymin(-5.),
120  dy(0.005)
121 {
122  nist = G4NistManager::Instance();
123 
124  theElectron = G4Electron::Electron();
125  thePositron = G4Positron::Positron();
126 
127  particleMass = lnZ = z13 = z23 = 0;
128 
129  // setup lowest limit dependent on particle mass
130  if(p) {
131  SetParticle(p);
132  G4double limit = p->GetPDGMass()*8;
133  if(limit > lowestKinEnergy) { lowestKinEnergy = limit; }
134  }
135  emin = lowestKinEnergy;
136  emax = 10*TeV;
137 }
G4MuPairProductionModel::particle
const G4ParticleDefinition * particle
Definition: G4MuPairProductionModel.hh:144
G4MuPairProductionModel::SetParticle
void SetParticle(const G4ParticleDefinition *)
Definition: G4MuPairProductionModel.hh:191
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:65
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4double
double G4double
Definition: G4Types.hh:76
G4MuPairProductionModel::~G4MuPairProductionModel ( )
virtual

Definition at line 141 of file G4MuPairProductionModel.cc.

142 {}

Member Function Documentation

G4double G4MuPairProductionModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  ,
G4double  kineticEnergy,
G4double  Z,
G4double  A,
G4double  cutEnergy,
G4double  maxEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 415 of file G4MuPairProductionModel.cc.

References ComputeMicroscopicCrossSection(), G4INCL::Math::max(), MaxSecondaryEnergyForElement(), and G4INCL::Math::min().

421 {
422  G4double cross = 0.0;
423  if (kineticEnergy <= lowestKinEnergy) { return cross; }
424 
425  G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy, Z);
426  G4double tmax = std::min(maxEnergy, maxPairEnergy);
427  G4double cut = std::max(cutEnergy, minPairEnergy);
428  if (cut >= tmax) { return cross; }
429 
430  cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
431  if(tmax < kineticEnergy) {
432  cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
433  }
434  return cross;
435 }
G4MuPairProductionModel::MaxSecondaryEnergyForElement
G4double MaxSecondaryEnergyForElement(G4double kineticEnergy, G4double Z)
Definition: G4MuPairProductionModel.hh:202
G4MuPairProductionModel::ComputeMicroscopicCrossSection
G4double ComputeMicroscopicCrossSection(G4double tkin, G4double Z, G4double cut)
Definition: G4MuPairProductionModel.cc:264
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4double
double G4double
Definition: G4Types.hh:76
G4double G4MuPairProductionModel::ComputeDEDXPerVolume ( const G4Material *  material,
const G4ParticleDefinition *  ,
G4double  kineticEnergy,
G4double  cutEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 196 of file G4MuPairProductionModel.cc.

References ComputMuPairLoss(), G4Material::GetAtomicNumDensityVector(), G4Material::GetElementVector(), G4Material::GetNumberOfElements(), and MaxSecondaryEnergyForElement().

201 {
202  G4double dedx = 0.0;
203  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
204  { return dedx; }
205 
206  const G4ElementVector* theElementVector = material->GetElementVector();
207  const G4double* theAtomicNumDensityVector =
208  material->GetAtomicNumDensityVector();
209 
210  // loop for elements in the material
211  for (size_t i=0; i<material->GetNumberOfElements(); ++i) {
212  G4double Z = (*theElementVector)[i]->GetZ();
213  G4double tmax = MaxSecondaryEnergyForElement(kineticEnergy, Z);
214  G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
215  dedx += loss*theAtomicNumDensityVector[i];
216  }
217  if (dedx < 0.) { dedx = 0.; }
218  return dedx;
219 }
G4MuPairProductionModel::MaxSecondaryEnergyForElement
G4double MaxSecondaryEnergyForElement(G4double kineticEnergy, G4double Z)
Definition: G4MuPairProductionModel.hh:202
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
G4MuPairProductionModel::ComputMuPairLoss
G4double ComputMuPairLoss(G4double Z, G4double tkin, G4double cut, G4double tmax)
Definition: G4MuPairProductionModel.cc:223
G4Material::GetAtomicNumDensityVector
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:214
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
G4double
double G4double
Definition: G4Types.hh:76
G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection ( G4double  tkin,
G4double  Z,
G4double  pairEnergy 
)
protectedvirtual

Reimplemented in G4hPairProductionModel.

Definition at line 302 of file G4MuPairProductionModel.cc.

References python.hepunit::electron_mass_c2, factorForCross, fe, fm, G4Exp(), G4Log(), particleMass, sqrte, wgi, xgi, G4InuclParticleNames::xi0, z13, and z23.

Referenced by ComputeMicroscopicCrossSection(), and ComputMuPairLoss().

309 {
310  static const G4double bbbtf= 183. ;
311  static const G4double bbbh = 202.4 ;
312  static const G4double g1tf = 1.95e-5 ;
313  static const G4double g2tf = 5.3e-5 ;
314  static const G4double g1h = 4.4e-5 ;
315  static const G4double g2h = 4.8e-5 ;
316 
317  G4double totalEnergy = tkin + particleMass;
318  G4double residEnergy = totalEnergy - pairEnergy;
319  G4double massratio = particleMass/electron_mass_c2 ;
320  G4double massratio2 = massratio*massratio ;
321  G4double cross = 0.;
322 
323  G4double c3 = 0.75*sqrte*particleMass;
324  if (residEnergy <= c3*z13) { return cross; }
325 
326  G4double c7 = 4.*CLHEP::electron_mass_c2;
327  G4double c8 = 6.*particleMass*particleMass;
328  G4double alf = c7/pairEnergy;
329  G4double a3 = 1. - alf;
330  if (a3 <= 0.) { return cross; }
331 
332  // zeta calculation
333  G4double bbb,g1,g2;
334  if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
335  else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
336 
337  G4double zeta = 0;
338  G4double zeta1 =
339  0.073*G4Log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26;
340  if ( zeta1 > 0.)
341  {
342  G4double zeta2 =
343  0.058*G4Log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14;
344  zeta = zeta1/zeta2 ;
345  }
346 
347  G4double z2 = Z*(Z+zeta);
348  G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
349  G4double a0 = totalEnergy*residEnergy;
350  G4double a1 = pairEnergy*pairEnergy/a0;
351  G4double bet = 0.5*a1;
352  G4double xi0 = 0.25*massratio2*a1;
353  G4double del = c8/a0;
354 
355  G4double rta3 = sqrt(a3);
356  G4double tmnexp = alf/(1. + rta3) + del*rta3;
357  if(tmnexp >= 1.0) { return cross; }
358 
359  G4double tmn = G4Log(tmnexp);
360  G4double sum = 0.;
361 
362  // Gaussian integration in ln(1-ro) ( with 8 points)
363  for (G4int i=0; i<8; ++i)
364  {
365  G4double a4 = G4Exp(tmn*xgi[i]); // a4 = (1.-asymmetry)
366  G4double a5 = a4*(2.-a4) ;
367  G4double a6 = 1.-a5 ;
368  G4double a7 = 1.+a6 ;
369  G4double a9 = 3.+a6 ;
370  G4double xi = xi0*a5 ;
371  G4double xii = 1./xi ;
372  G4double xi1 = 1.+xi ;
373  G4double screen = screen0*xi1/a5 ;
374  G4double yeu = 5.-a6+4.*bet*a7 ;
375  G4double yed = 2.*(1.+3.*bet)*G4Log(3.+xii)-a6-a1*(2.-a6) ;
376  G4double ye1 = 1.+yeu/yed ;
377  G4double ale = G4Log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
378  G4double cre = 0.5*G4Log(1.+2.25*z23*xi1*ye1/massratio2) ;
379  G4double be;
380 
381  if (xi <= 1.e3) {
382  be = ((2.+a6)*(1.+bet)+xi*a9)*G4Log(1.+xii)+(a5-bet)/xi1-a9;
383  } else {
384  be = (3.-a6+a1*a7)/(2.*xi);
385  }
386  G4double fe = (ale-cre)*be;
387  if ( fe < 0.) fe = 0. ;
388 
389  G4double ymu = 4.+a6 +3.*bet*a7 ;
390  G4double ymd = a7*(1.5+a1)*G4Log(3.+xi)+1.-1.5*a6 ;
391  G4double ym1 = 1.+ymu/ymd ;
392  G4double alm_crm = G4Log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
393  G4double a10,bm;
394  if ( xi >= 1.e-3)
395  {
396  a10 = (1.+a1)*a5 ;
397  bm = (a7*(1.+1.5*bet)-a10*xii)*G4Log(xi1)+xi*(a5-bet)/xi1+a10;
398  } else {
399  bm = (5.-a6+bet*a9)*(xi/2.);
400  }
401 
402  G4double fm = alm_crm*bm;
403  if ( fm < 0.) { fm = 0.; }
404 
405  sum += wgi[i]*a4*(fe+fm/massratio2);
406  }
407 
408  cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
409  if(cross < 0.0) { cross = 0.0; }
410  return cross;
411 }
G4MuPairProductionModel::xgi
static const G4double xgi[8]
Definition: G4MuPairProductionModel.hh:155
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
fm
#define fm
G4MuPairProductionModel::wgi
static const G4double wgi[8]
Definition: G4MuPairProductionModel.hh:155
G4double
double G4double
Definition: G4Types.hh:76
fe
G4fissionEvent * fe
Definition: G4LLNLFission.cc:63
G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection ( G4double  tkin,
G4double  Z,
G4double  cut 
)
protected

Definition at line 264 of file G4MuPairProductionModel.cc.

References ComputeDMicroscopicCrossSection(), G4Exp(), G4Log(), G4INCL::Math::max(), MaxSecondaryEnergyForElement(), wgi, test::x, and xgi.

Referenced by ComputeCrossSectionPerAtom().

268 {
269  G4double cross = 0.;
270  G4double tmax = MaxSecondaryEnergyForElement(tkin, Z);
271  G4double cut = std::max(cutEnergy, minPairEnergy);
272  if (tmax <= cut) { return cross; }
273 
274  G4double ak1=6.9 ;
275  G4double ak2=1.0 ;
276  G4double aaa = G4Log(cut);
277  G4double bbb = G4Log(tmax);
278  G4int kkk = (G4int)((bbb-aaa)/ak1 + ak2);
279  if(kkk > 8) { kkk = 8; }
280  else if (kkk < 1) { kkk = 1; }
281 
282  G4double hhh = (bbb-aaa)/G4double(kkk);
283  G4double x = aaa;
284 
285  for(G4int l=0; l<kkk; ++l)
286  {
287  for(G4int i=0; i<8; ++i)
288  {
289  G4double ep = G4Exp(x + xgi[i]*hhh);
290  cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
291  }
292  x += hhh;
293  }
294 
295  cross *= hhh;
296  if(cross < 0.0) { cross = 0.0; }
297  return cross;
298 }
G4MuPairProductionModel::MaxSecondaryEnergyForElement
G4double MaxSecondaryEnergyForElement(G4double kineticEnergy, G4double Z)
Definition: G4MuPairProductionModel.hh:202
G4MuPairProductionModel::xgi
static const G4double xgi[8]
Definition: G4MuPairProductionModel.hh:155
G4int
int G4int
Definition: G4Types.hh:78
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4MuPairProductionModel::ComputeDMicroscopicCrossSection
virtual G4double ComputeDMicroscopicCrossSection(G4double tkin, G4double Z, G4double pairEnergy)
Definition: G4MuPairProductionModel.cc:302
G4MuPairProductionModel::wgi
static const G4double wgi[8]
Definition: G4MuPairProductionModel.hh:155
G4double
double G4double
Definition: G4Types.hh:76
G4double G4MuPairProductionModel::ComputMuPairLoss ( G4double  Z,
G4double  tkin,
G4double  cut,
G4double  tmax 
)
protected

Definition at line 223 of file G4MuPairProductionModel.cc.

References ComputeDMicroscopicCrossSection(), G4Exp(), G4Log(), G4INCL::Math::min(), wgi, test::x, and xgi.

Referenced by ComputeDEDXPerVolume().

227 {
228  G4double loss = 0.0;
229 
230  G4double cut = std::min(cutEnergy,tmax);
231  if(cut <= minPairEnergy) { return loss; }
232 
233  // calculate the rectricted loss
234  // numerical integration in log(PairEnergy)
235  G4double ak1=6.9;
236  G4double ak2=1.0;
237  G4double aaa = G4Log(minPairEnergy);
238  G4double bbb = G4Log(cut);
239 
240  G4int kkk = (G4int)((bbb-aaa)/ak1+ak2);
241  if (kkk > 8) { kkk = 8; }
242  else if (kkk < 1) { kkk = 1; }
243 
244  G4double hhh = (bbb-aaa)/(G4double)kkk;
245  G4double x = aaa;
246 
247  for (G4int l=0 ; l<kkk; l++)
248  {
249 
250  for (G4int ll=0; ll<8; ll++)
251  {
252  G4double ep = G4Exp(x+xgi[ll]*hhh);
253  loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
254  }
255  x += hhh;
256  }
257  loss *= hhh;
258  if (loss < 0.) loss = 0.;
259  return loss;
260 }
G4MuPairProductionModel::xgi
static const G4double xgi[8]
Definition: G4MuPairProductionModel.hh:155
G4int
int G4int
Definition: G4Types.hh:78
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4MuPairProductionModel::ComputeDMicroscopicCrossSection
virtual G4double ComputeDMicroscopicCrossSection(G4double tkin, G4double Z, G4double pairEnergy)
Definition: G4MuPairProductionModel.cc:302
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4MuPairProductionModel::wgi
static const G4double wgi[8]
Definition: G4MuPairProductionModel.hh:155
G4double
double G4double
Definition: G4Types.hh:76
void G4MuPairProductionModel::Initialise ( const G4ParticleDefinition *  p,
const G4DataVector &  cuts 
)
virtual

Implements G4VEmModel.

Definition at line 155 of file G4MuPairProductionModel.cc.

References G4VEmModel::fElementData, G4Log(), G4VEmModel::GetParticleChangeForLoss(), G4VEmModel::HighEnergyLimit(), G4VEmModel::InitialiseElementSelectors(), G4VEmModel::IsMaster(), G4VEmModel::LowEnergyLimit(), G4INCL::Math::max(), particle, and SetParticle().

157 {
158  SetParticle(p);
159 
160  // define scale of internal table for each thread only once
161  if(0 == nbine) {
162  emax = HighEnergyLimit();
163  emin = std::max(lowestKinEnergy, LowEnergyLimit());
164  nbine = size_t(nYBinPerDecade*std::log10(emax/emin));
165  if(nbine < 3) { nbine = 3; }
166 
167  ymin = G4Log(minPairEnergy/emin);
168  dy = -ymin/G4double(nbiny);
169  }
170 
171  if(IsMaster() && p == particle) {
172 
173  if(!fElementData) {
174  fElementData = new G4ElementData();
175  MakeSamplingTables();
176  }
177  InitialiseElementSelectors(p, cuts);
178  }
179 
180  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
181 }
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
G4VEmModel::GetParticleChangeForLoss
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:107
G4VEmModel::InitialiseElementSelectors
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:135
G4MuPairProductionModel::particle
const G4ParticleDefinition * particle
Definition: G4MuPairProductionModel.hh:144
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
G4VEmModel::IsMaster
G4bool IsMaster() const
Definition: G4VEmModel.hh:676
G4MuPairProductionModel::SetParticle
void SetParticle(const G4ParticleDefinition *)
Definition: G4MuPairProductionModel.hh:191
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4double
double G4double
Definition: G4Types.hh:76
G4VEmModel::fElementData
G4ElementData * fElementData
Definition: G4VEmModel.hh:397
void G4MuPairProductionModel::InitialiseLocal ( const G4ParticleDefinition *  p,
G4VEmModel *  masterModel 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 185 of file G4MuPairProductionModel.cc.

References G4VEmModel::fElementData, G4VEmModel::GetElementData(), G4VEmModel::GetElementSelectors(), particle, and G4VEmModel::SetElementSelectors().

187 {
188  if(p == particle) {
189  SetElementSelectors(masterModel->GetElementSelectors());
190  fElementData = masterModel->GetElementData();
191  }
192 }
G4MuPairProductionModel::particle
const G4ParticleDefinition * particle
Definition: G4MuPairProductionModel.hh:144
G4VEmModel::GetElementData
G4ElementData * GetElementData()
Definition: G4VEmModel.hh:777
G4VEmModel::GetElementSelectors
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:760
G4VEmModel::SetElementSelectors
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:768
G4VEmModel::fElementData
G4ElementData * fElementData
Definition: G4VEmModel.hh:397
G4double G4MuPairProductionModel::MaxSecondaryEnergyForElement ( G4double  kineticEnergy,
G4double  Z 
)
inlineprotected

Definition at line 202 of file G4MuPairProductionModel.hh.

References currentZ, G4lrint(), G4NistManager::GetLOGZ(), G4NistManager::GetZ13(), lnZ, nist, particleMass, sqrte, z13, and z23.

Referenced by ComputeCrossSectionPerAtom(), ComputeDEDXPerVolume(), ComputeMicroscopicCrossSection(), and SampleSecondaries().

204 {
205  G4int Z = G4lrint(ZZ);
206  if(Z != currentZ) {
207  currentZ = Z;
208  z13 = nist->GetZ13(Z);
209  z23 = z13*z13;
210  lnZ = nist->GetLOGZ(Z);
211  }
212  return kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
213 }
G4NistManager::GetZ13
G4double GetZ13(G4double Z)
Definition: G4NistManager.hh:528
G4int
int G4int
Definition: G4Types.hh:78
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
G4NistManager::GetLOGZ
G4double GetLOGZ(G4int Z)
Definition: G4NistManager.hh:551
G4double G4MuPairProductionModel::MinPrimaryEnergy ( const G4Material *  ,
const G4ParticleDefinition *  ,
G4double  cut 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 146 of file G4MuPairProductionModel.cc.

References G4INCL::Math::max().

149 {
150  return std::max(lowestKinEnergy,cut);
151 }
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
void G4MuPairProductionModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  vdp,
const G4MaterialCutsCouple *  couple,
const G4DynamicParticle *  aDynamicParticle,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 508 of file G4MuPairProductionModel.cc.

References currentZ, python.hepunit::electron_mass_c2, G4Exp(), G4Log(), G4UniformRand, G4InuclParticleNames::gam, G4DynamicParticle::GetKineticEnergy(), G4NistManager::GetLOGZ(), G4DynamicParticle::GetMomentum(), G4DynamicParticle::GetMomentumDirection(), G4Element::GetZ(), iz, lnZ, G4INCL::Math::max(), MaxSecondaryEnergyForElement(), python.hepunit::MeV, G4INCL::Math::min(), nist, particle, particleMass, CLHEP::Hep3Vector::rotateUz(), G4VEmModel::SelectRandomAtom(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), G4ParticleChangeForLoss::SetProposedMomentumDirection(), python.hepunit::twopi, CLHEP::Hep3Vector::unit(), and test::x.

514 {
515  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
516  //G4cout << "------- G4MuPairProductionModel::SampleSecondaries E(MeV)= "
517  // << kineticEnergy << " "
518  // << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
519  G4double totalEnergy = kineticEnergy + particleMass;
520  G4double totalMomentum =
521  sqrt(kineticEnergy*(kineticEnergy + 2.0*particleMass));
522 
523  G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
524 
525  // select randomly one element constituing the material
526  const G4Element* anElement = SelectRandomAtom(couple,particle,kineticEnergy);
527 
528  // define interval of energy transfer
529  G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy,
530  anElement->GetZ());
531  G4double maxEnergy = std::min(tmax, maxPairEnergy);
532  G4double minEnergy = std::max(tmin, minPairEnergy);
533 
534  if(minEnergy >= maxEnergy) { return; }
535  //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
536  // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
537  // << " ymin= " << ymin << " dy= " << dy << G4endl;
538 
539  G4double coeff = G4Log(minPairEnergy/kineticEnergy)/ymin;
540 
541  // compute limits
542  G4double yymin = G4Log(minEnergy/kineticEnergy)/coeff;
543  G4double yymax = G4Log(maxEnergy/kineticEnergy)/coeff;
544 
545  //G4cout << "yymin= " << yymin << " yymax= " << yymax << G4endl;
546 
547  // units should not be used, bacause table was built without
548  G4double logTkin = G4Log(kineticEnergy/MeV);
549 
550  // sample e-e+ energy, pair energy first
551 
552  // select sample table via Z
553  G4int iz1(0), iz2(0);
554  for(G4int iz=0; iz<nzdat; ++iz) {
555  if(currentZ == zdat[iz]) {
556  iz1 = iz2 = currentZ;
557  break;
558  } else if(currentZ < zdat[iz]) {
559  iz2 = zdat[iz];
560  if(iz > 0) { iz1 = zdat[iz-1]; }
561  else { iz1 = iz2; }
562  break;
563  }
564  }
565  if(0 == iz1) { iz1 = iz2 = zdat[nzdat-1]; }
566 
567  G4double PairEnergy = 0.0;
568  G4int count = 0;
569  //G4cout << "start loop Z1= " << iz1 << " Z2= " << iz2 << G4endl;
570  do {
571  ++count;
572  // sampling using only one random number
573  G4double rand = G4UniformRand();
574 
575  G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
576  if(iz1 != iz2) {
577  G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
578  G4double lz1= nist->GetLOGZ(iz1);
579  G4double lz2= nist->GetLOGZ(iz2);
580  //G4cout << count << ". x= " << x << " x2= " << x2
581  // << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
582  x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
583  }
584  //G4cout << "x= " << x << " coeff= " << coeff << G4endl;
585  PairEnergy = kineticEnergy*G4Exp(x*coeff);
586 
587  } while((PairEnergy < minEnergy || PairEnergy > maxEnergy) && 10 > count);
588 
589  //G4cout << "## PairEnergy(GeV)= " << PairEnergy/GeV
590  // << " Etot(GeV)= " << totalEnergy/GeV << G4endl;
591 
592  // sample r=(E+-E-)/PairEnergy ( uniformly .....)
593  G4double rmax =
594  (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
595  *sqrt(1.-minPairEnergy/PairEnergy);
596  G4double r = rmax * (-1.+2.*G4UniformRand()) ;
597 
598  // compute energies from PairEnergy,r
599  G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
600  G4double PositronEnergy = PairEnergy - ElectronEnergy;
601 
602  // The angle of the emitted virtual photon is sampled
603  // according to the muon bremsstrahlung model
604 
605  G4double gam = totalEnergy/particleMass;
606  G4double gmax = gam*std::min(1.0, totalEnergy/PairEnergy - 1.0);
607  G4double gmax2= gmax*gmax;
608  G4double x = G4UniformRand()*gmax2/(1.0 + gmax2);
609 
610  G4double theta = sqrt(x/(1.0 - x))/gam;
611  G4double sint = sin(theta);
612  G4double phi = twopi * G4UniformRand() ;
613  G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
614 
615  G4ThreeVector gDirection(dirx, diry, dirz);
616  gDirection.rotateUz(partDirection);
617 
618  // the angles of e- and e+ assumed to be the same as virtual gamma
619 
620  // create G4DynamicParticle object for the particle1
621  G4DynamicParticle* aParticle1 =
622  new G4DynamicParticle(theElectron, gDirection,
623  ElectronEnergy - electron_mass_c2);
624 
625  // create G4DynamicParticle object for the particle2
626  G4DynamicParticle* aParticle2 =
627  new G4DynamicParticle(thePositron, gDirection,
628  PositronEnergy - electron_mass_c2);
629 
630  // primary change
631  kineticEnergy -= (ElectronEnergy + PositronEnergy);
632  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
633 
634  partDirection *= totalMomentum;
635  partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
636  partDirection = partDirection.unit();
637  fParticleChange->SetProposedMomentumDirection(partDirection);
638 
639  // add secondary
640  vdp->push_back(aParticle1);
641  vdp->push_back(aParticle2);
642  //G4cout << "-- G4MuPairProductionModel::SampleSecondaries done" << G4endl;
643 }
G4MuPairProductionModel::MaxSecondaryEnergyForElement
G4double MaxSecondaryEnergyForElement(G4double kineticEnergy, G4double Z)
Definition: G4MuPairProductionModel.hh:202
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4MuPairProductionModel::particle
const G4ParticleDefinition * particle
Definition: G4MuPairProductionModel.hh:144
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
iz
G4double iz
Definition: TRTMaterials.hh:39
G4ParticleChangeForLoss::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForLoss.hh:150
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4ParticleChangeForLoss::SetProposedMomentumDirection
void SetProposedMomentumDirection(const G4ThreeVector &dir)
Definition: G4ParticleChangeForLoss.hh:194
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4NistManager::GetLOGZ
G4double GetLOGZ(G4int Z)
Definition: G4NistManager.hh:551
G4double
double G4double
Definition: G4Types.hh:76
G4VEmModel::SelectRandomAtom
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:510
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
void G4MuPairProductionModel::SetLowestKineticEnergy ( G4double  e)
inline

Definition at line 183 of file G4MuPairProductionModel.hh.

184 {
185  lowestKinEnergy = e;
186 }
void G4MuPairProductionModel::SetParticle ( const G4ParticleDefinition *  p)
inline

Definition at line 191 of file G4MuPairProductionModel.hh.

References G4ParticleDefinition::GetPDGMass(), particle, and particleMass.

Referenced by G4MuPairProductionModel(), and Initialise().

192 {
193  if(!particle) {
194  particle = p;
195  particleMass = particle->GetPDGMass();
196  }
197 }
G4MuPairProductionModel::particle
const G4ParticleDefinition * particle
Definition: G4MuPairProductionModel.hh:144
p
const char * p
Definition: xmltok.h:285
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161

Field Documentation

G4int G4MuPairProductionModel::currentZ
protected

Definition at line 153 of file G4MuPairProductionModel.hh.

Referenced by MaxSecondaryEnergyForElement(), and SampleSecondaries().

G4double G4MuPairProductionModel::factorForCross
protected

Definition at line 147 of file G4MuPairProductionModel.hh.

Referenced by G4hPairProductionModel::ComputeDMicroscopicCrossSection(), and ComputeDMicroscopicCrossSection().

G4double G4MuPairProductionModel::lnZ
protected

Definition at line 152 of file G4MuPairProductionModel.hh.

Referenced by G4MuPairProductionModel(), MaxSecondaryEnergyForElement(), and SampleSecondaries().

G4NistManager* G4MuPairProductionModel::nist
protected

Definition at line 145 of file G4MuPairProductionModel.hh.

Referenced by G4MuPairProductionModel(), MaxSecondaryEnergyForElement(), and SampleSecondaries().

const G4ParticleDefinition* G4MuPairProductionModel::particle
protected

Definition at line 144 of file G4MuPairProductionModel.hh.

Referenced by Initialise(), InitialiseLocal(), SampleSecondaries(), and SetParticle().

G4double G4MuPairProductionModel::particleMass
protected

Definition at line 149 of file G4MuPairProductionModel.hh.

Referenced by G4hPairProductionModel::ComputeDMicroscopicCrossSection(), ComputeDMicroscopicCrossSection(), G4MuPairProductionModel(), MaxSecondaryEnergyForElement(), SampleSecondaries(), and SetParticle().

G4double G4MuPairProductionModel::sqrte
protected

Definition at line 148 of file G4MuPairProductionModel.hh.

Referenced by G4hPairProductionModel::ComputeDMicroscopicCrossSection(), ComputeDMicroscopicCrossSection(), and MaxSecondaryEnergyForElement().

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

Definition at line 155 of file G4MuPairProductionModel.hh.

Referenced by G4hPairProductionModel::ComputeDMicroscopicCrossSection(), ComputeDMicroscopicCrossSection(), ComputeMicroscopicCrossSection(), and ComputMuPairLoss().

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

Definition at line 155 of file G4MuPairProductionModel.hh.

Referenced by G4hPairProductionModel::ComputeDMicroscopicCrossSection(), ComputeDMicroscopicCrossSection(), ComputeMicroscopicCrossSection(), and ComputMuPairLoss().

G4double G4MuPairProductionModel::z13
protected

Definition at line 150 of file G4MuPairProductionModel.hh.

Referenced by G4hPairProductionModel::ComputeDMicroscopicCrossSection(), ComputeDMicroscopicCrossSection(), G4MuPairProductionModel(), and MaxSecondaryEnergyForElement().

G4double G4MuPairProductionModel::z23
protected

Definition at line 151 of file G4MuPairProductionModel.hh.

Referenced by G4hPairProductionModel::ComputeDMicroscopicCrossSection(), ComputeDMicroscopicCrossSection(), G4MuPairProductionModel(), and MaxSecondaryEnergyForElement().

The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:45 for Geant4.10 by   doxygen 1.8.7