Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4MuPairProductionModel Class Reference

#include <G4MuPairProductionModel.hh>

+ Inheritance diagram for G4MuPairProductionModel:

Public Member Functions

 G4MuPairProductionModel (const G4ParticleDefinition *p=0, const G4String &nam="muPairProd")
 
virtual ~G4MuPairProductionModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
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 MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
void SetLowestKineticEnergy (G4double e)
 
void SetParticle (const G4ParticleDefinition *)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
 
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 ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., 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 *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
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)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *)
 
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
 
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 ForceBuildTable (G4bool val)
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () 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)
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
void SetCurrentElement (G4double Z)
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 

Protected Attributes

const G4ParticleDefinitionparticle
 
G4NistManagernist
 
G4double factorForCross
 
G4double sqrte
 
G4double particleMass
 
G4double currentZ
 
G4double z13
 
G4double z23
 
G4double lnZ
 
- Protected Attributes inherited from G4VEmModel
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 

Static Protected Attributes

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

Detailed Description

Definition at line 71 of file G4MuPairProductionModel.hh.

Constructor & Destructor Documentation

◆ G4MuPairProductionModel()

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

Definition at line 104 of file G4MuPairProductionModel.cc.

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(exp(1.))),
111 currentZ(0),
112 fParticleChange(0),
113 minPairEnergy(4.*electron_mass_c2),
114 lowestKinEnergy(GeV),
115 nzdat(5),
116 ntdat(8),
117 nbiny(1000),
118 nmaxElements(0),
119 ymin(-5.),
120 ymax(0.),
121 dy((ymax-ymin)/nbiny),
122 samplingTablesAreFilled(false)
123{
124 SetLowEnergyLimit(minPairEnergy);
126
127 theElectron = G4Electron::Electron();
128 thePositron = G4Positron::Positron();
129
130 particleMass = lnZ = z13 = z23 = 0;
131
132 for(size_t i=0; i<1001; ++i) { ya[i] = 0.0; }
133
134 if(p) { SetParticle(p); }
135}
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetParticle(const G4ParticleDefinition *)
const G4ParticleDefinition * particle
static G4NistManager * Instance()
static G4Positron * Positron()
Definition: G4Positron.cc:94
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:592

◆ ~G4MuPairProductionModel()

G4MuPairProductionModel::~G4MuPairProductionModel ( )
virtual

Definition at line 139 of file G4MuPairProductionModel.cc.

140{}

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

Reimplemented from G4VEmModel.

Definition at line 385 of file G4MuPairProductionModel.cc.

391{
392 G4double cross = 0.0;
393 if (kineticEnergy <= lowestKinEnergy) { return cross; }
394
396
397 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
398 G4double tmax = std::min(maxEnergy, maxPairEnergy);
399 G4double cut = std::max(cutEnergy, minPairEnergy);
400 if (cut >= tmax) return cross;
401
402 cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
403 if(tmax < kineticEnergy) {
404 cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
405 }
406 return cross;
407}
double G4double
Definition: G4Types.hh:64
G4double ComputeMicroscopicCrossSection(G4double tkin, G4double Z, G4double cut)
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kineticEnergy)

◆ ComputeDEDXPerVolume()

G4double G4MuPairProductionModel::ComputeDEDXPerVolume ( const G4Material material,
const G4ParticleDefinition ,
G4double  kineticEnergy,
G4double  cutEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 173 of file G4MuPairProductionModel.cc.

178{
179 G4double dedx = 0.0;
180 if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
181 { return dedx; }
182
183 const G4ElementVector* theElementVector = material->GetElementVector();
184 const G4double* theAtomicNumDensityVector =
185 material->GetAtomicNumDensityVector();
186
187 // loop for elements in the material
188 for (size_t i=0; i<material->GetNumberOfElements(); ++i) {
189 G4double Z = (*theElementVector)[i]->GetZ();
191 G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
192 G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
193 dedx += loss*theAtomicNumDensityVector[i];
194 }
195 if (dedx < 0.) { dedx = 0.; }
196 return dedx;
197}
std::vector< G4Element * > G4ElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:215
G4double ComputMuPairLoss(G4double Z, G4double tkin, G4double cut, G4double tmax)

◆ ComputeDMicroscopicCrossSection()

G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection ( G4double  tkin,
G4double  Z,
G4double  pairEnergy 
)
protectedvirtual

Reimplemented in G4hPairProductionModel.

Definition at line 274 of file G4MuPairProductionModel.cc.

281{
282 G4double bbbtf= 183. ;
283 G4double bbbh = 202.4 ;
284 G4double g1tf = 1.95e-5 ;
285 G4double g2tf = 5.3e-5 ;
286 G4double g1h = 4.4e-5 ;
287 G4double g2h = 4.8e-5 ;
288
289 G4double totalEnergy = tkin + particleMass;
290 G4double residEnergy = totalEnergy - pairEnergy;
291 G4double massratio = particleMass/electron_mass_c2 ;
292 G4double massratio2 = massratio*massratio ;
293 G4double cross = 0.;
294
296
297 G4double c3 = 0.75*sqrte*particleMass;
298 if (residEnergy <= c3*z13) { return cross; }
299
300 G4double c7 = 4.*electron_mass_c2;
302 G4double alf = c7/pairEnergy;
303 G4double a3 = 1. - alf;
304 if (a3 <= 0.) { return cross; }
305
306 // zeta calculation
307 G4double bbb,g1,g2;
308 if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
309 else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
310
311 G4double zeta = 0;
312 G4double zeta1 = 0.073*log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26;
313 if ( zeta1 > 0.)
314 {
315 G4double zeta2 = 0.058*log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14;
316 zeta = zeta1/zeta2 ;
317 }
318
319 G4double z2 = Z*(Z+zeta);
320 G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
321 G4double a0 = totalEnergy*residEnergy;
322 G4double a1 = pairEnergy*pairEnergy/a0;
323 G4double bet = 0.5*a1;
324 G4double xi0 = 0.25*massratio2*a1;
325 G4double del = c8/a0;
326
327 G4double rta3 = sqrt(a3);
328 G4double tmnexp = alf/(1. + rta3) + del*rta3;
329 if(tmnexp >= 1.0) return cross;
330
331 G4double tmn = log(tmnexp);
332 G4double sum = 0.;
333
334 // Gaussian integration in ln(1-ro) ( with 8 points)
335 for (G4int i=0; i<8; ++i)
336 {
337 G4double a4 = exp(tmn*xgi[i]); // a4 = (1.-asymmetry)
338 G4double a5 = a4*(2.-a4) ;
339 G4double a6 = 1.-a5 ;
340 G4double a7 = 1.+a6 ;
341 G4double a9 = 3.+a6 ;
342 G4double xi = xi0*a5 ;
343 G4double xii = 1./xi ;
344 G4double xi1 = 1.+xi ;
345 G4double screen = screen0*xi1/a5 ;
346 G4double yeu = 5.-a6+4.*bet*a7 ;
347 G4double yed = 2.*(1.+3.*bet)*log(3.+xii)-a6-a1*(2.-a6) ;
348 G4double ye1 = 1.+yeu/yed ;
349 G4double ale=log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
350 G4double cre = 0.5*log(1.+2.25*z23*xi1*ye1/massratio2) ;
351 G4double be;
352
353 if (xi <= 1.e3) be = ((2.+a6)*(1.+bet)+xi*a9)*log(1.+xii)+(a5-bet)/xi1-a9;
354 else be = (3.-a6+a1*a7)/(2.*xi);
355
356 G4double fe = (ale-cre)*be;
357 if ( fe < 0.) fe = 0. ;
358
359 G4double ymu = 4.+a6 +3.*bet*a7 ;
360 G4double ymd = a7*(1.5+a1)*log(3.+xi)+1.-1.5*a6 ;
361 G4double ym1 = 1.+ymu/ymd ;
362 G4double alm_crm = log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
363 G4double a10,bm;
364 if ( xi >= 1.e-3)
365 {
366 a10 = (1.+a1)*a5 ;
367 bm = (a7*(1.+1.5*bet)-a10*xii)*log(xi1)+xi*(a5-bet)/xi1+a10;
368 } else {
369 bm = (5.-a6+bet*a9)*(xi/2.);
370 }
371
372 G4double fm = alm_crm*bm;
373 if ( fm < 0.) fm = 0. ;
374
375 sum += wgi[i]*a4*(fe+fm/massratio2);
376 }
377
378 cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
379
380 return cross;
381}
G4fissionEvent * fe
int G4int
Definition: G4Types.hh:66

Referenced by ComputeMicroscopicCrossSection(), and ComputMuPairLoss().

◆ ComputeMicroscopicCrossSection()

G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection ( G4double  tkin,
G4double  Z,
G4double  cut 
)
protected

Definition at line 240 of file G4MuPairProductionModel.cc.

244{
245 G4double cross = 0.;
248 if (tmax <= cut) { return cross; }
249
250 G4double ak1=6.9 ;
251 G4double ak2=1.0 ;
252 G4double aaa = log(cut);
253 G4double bbb = log(tmax);
254 G4int kkk = (G4int)((bbb-aaa)/ak1 + ak2);
255 if(kkk > 8) { kkk = 8; }
256 G4double hhh = (bbb-aaa)/G4double(kkk);
257 G4double x = aaa;
258
259 for(G4int l=0; l<kkk; ++l)
260 {
261 for(G4int i=0; i<8; ++i)
262 {
263 G4double ep = exp(x + xgi[i]*hhh);
264 cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
265 }
266 x += hhh;
267 }
268
269 cross *= hhh;
270 if(cross < 0.0) { cross = 0.0; }
271 return cross;
272}
virtual G4double ComputeDMicroscopicCrossSection(G4double tkin, G4double Z, G4double pairEnergy)

Referenced by ComputeCrossSectionPerAtom().

◆ ComputMuPairLoss()

G4double G4MuPairProductionModel::ComputMuPairLoss ( G4double  Z,
G4double  tkin,
G4double  cut,
G4double  tmax 
)
protected

Definition at line 201 of file G4MuPairProductionModel.cc.

205{
207 G4double loss = 0.0;
208
209 G4double cut = std::min(cutEnergy,tmax);
210 if(cut <= minPairEnergy) { return loss; }
211
212 // calculate the rectricted loss
213 // numerical integration in log(PairEnergy)
214 G4double ak1=6.9;
215 G4double ak2=1.0;
216 G4double aaa = log(minPairEnergy);
217 G4double bbb = log(cut);
218 G4int kkk = (G4int)((bbb-aaa)/ak1+ak2);
219 if (kkk > 8) kkk = 8;
220 G4double hhh = (bbb-aaa)/(G4double)kkk;
221 G4double x = aaa;
222
223 for (G4int l=0 ; l<kkk; l++)
224 {
225
226 for (G4int ll=0; ll<8; ll++)
227 {
228 G4double ep = exp(x+xgi[ll]*hhh);
229 loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
230 }
231 x += hhh;
232 }
233 loss *= hhh;
234 if (loss < 0.) loss = 0.;
235 return loss;
236}

Referenced by ComputeDEDXPerVolume().

◆ Initialise()

void G4MuPairProductionModel::Initialise ( const G4ParticleDefinition p,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 161 of file G4MuPairProductionModel.cc.

163{
164 if (!samplingTablesAreFilled) {
165 if(p) { SetParticle(p); }
166 MakeSamplingTables();
167 }
168 if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
169}
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:95

◆ MaxSecondaryEnergy()

G4double G4MuPairProductionModel::MaxSecondaryEnergy ( const G4ParticleDefinition ,
G4double  kineticEnergy 
)
protectedvirtual

Reimplemented from G4VEmModel.

Definition at line 152 of file G4MuPairProductionModel.cc.

154{
155 G4double maxPairEnergy = kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
156 return maxPairEnergy;
157}

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

◆ MinEnergyCut()

G4double G4MuPairProductionModel::MinEnergyCut ( const G4ParticleDefinition ,
const G4MaterialCutsCouple  
)
virtual

Definition at line 144 of file G4MuPairProductionModel.cc.

146{
147 return minPairEnergy;
148}

◆ SampleSecondaries()

void G4MuPairProductionModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  vdp,
const G4MaterialCutsCouple couple,
const G4DynamicParticle aDynamicParticle,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 466 of file G4MuPairProductionModel.cc.

471{
472 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
473 G4double totalEnergy = kineticEnergy + particleMass;
474 G4double totalMomentum =
475 sqrt(kineticEnergy*(kineticEnergy + 2.0*particleMass));
476
477 G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
478
479 G4int it;
480 for(it=1; it<ntdat; ++it) { if(kineticEnergy <= tdat[it]) { break; } }
481 if(it == ntdat) { --it; }
482 G4double dt = log(kineticEnergy/tdat[it-1])/log(tdat[it]/tdat[it-1]);
483
484 // select randomly one element constituing the material
485 const G4Element* anElement =
486 SelectRandomAtom(kineticEnergy, dt, it, couple, tmin);
487 SetCurrentElement(anElement->GetZ());
488
489 // define interval of enegry transfer
490 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
491 G4double maxEnergy = std::min(tmax, maxPairEnergy);
492 G4double minEnergy = std::max(tmin, minPairEnergy);
493
494 if(minEnergy >= maxEnergy) { return; }
495 //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
496 // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
497 // << " ymin= " << ymin << " dy= " << dy << G4endl;
498
499 G4double logmaxmin = log(maxPairEnergy/minPairEnergy);
500
501 // select bins
502 G4int iymin = 0;
503 G4int iymax = nbiny-1;
504 if( minEnergy > minPairEnergy)
505 {
506 G4double xc = log(minEnergy/minPairEnergy)/logmaxmin;
507 iymin = (G4int)((log(xc) - ymin)/dy);
508 if(iymin >= nbiny) iymin = nbiny-1;
509 else if(iymin < 0) iymin = 0;
510 xc = log(maxEnergy/minPairEnergy)/logmaxmin;
511 iymax = (G4int)((log(xc) - ymin)/dy) + 1;
512 if(iymax >= nbiny) iymax = nbiny-1;
513 else if(iymax < 0) iymax = 0;
514 }
515
516 // sample e-e+ energy, pair energy first
517 G4int iz, iy;
518
519 for(iz=1; iz<nzdat; ++iz) { if(currentZ <= zdat[iz]) { break; } }
520 if(iz == nzdat) { --iz; }
521
522 G4double dz = log(currentZ/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
523
524 G4double pmin = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymin,currentZ);
525 G4double pmax = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymax,currentZ);
526
527 G4double p = pmin+G4UniformRand()*(pmax - pmin);
528
529 // interpolate sampling vector;
530 G4double p1 = pmin;
531 G4double p2 = pmin;
532 for(iy=iymin+1; iy<=iymax; ++iy) {
533 p1 = p2;
534 p2 = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, currentZ);
535 if(p <= p2) { break; }
536 }
537 // G4cout << "iy= " << iy << " iymin= " << iymin << " iymax= "
538 // << iymax << " Z= " << currentZ << G4endl;
539 G4double y = ya[iy-1] + dy*(p - p1)/(p2 - p1);
540
541 G4double PairEnergy = minPairEnergy*exp( exp(y)*logmaxmin );
542
543 if(PairEnergy < minEnergy) { PairEnergy = minEnergy; }
544 if(PairEnergy > maxEnergy) { PairEnergy = maxEnergy; }
545
546 // sample r=(E+-E-)/PairEnergy ( uniformly .....)
547 G4double rmax =
548 (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
549 *sqrt(1.-minPairEnergy/PairEnergy);
550 G4double r = rmax * (-1.+2.*G4UniformRand()) ;
551
552 // compute energies from PairEnergy,r
553 G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
554 G4double PositronEnergy = PairEnergy - ElectronEnergy;
555
556 // The angle of the emitted virtual photon is sampled
557 // according to the muon bremsstrahlung model
558
559 G4double gam = totalEnergy/particleMass;
560 G4double gmax = gam*std::min(1.0, totalEnergy/PairEnergy - 1.0);
561 G4double gmax2= gmax*gmax;
562 G4double x = G4UniformRand()*gmax2/(1.0 + gmax2);
563
564 G4double theta = sqrt(x/(1.0 - x))/gam;
565 G4double sint = sin(theta);
566 G4double phi = twopi * G4UniformRand() ;
567 G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
568
569 G4ThreeVector gDirection(dirx, diry, dirz);
570 gDirection.rotateUz(partDirection);
571
572 // the angles of e- and e+ assumed to be the same as virtual gamma
573
574 // create G4DynamicParticle object for the particle1
575 G4DynamicParticle* aParticle1 =
576 new G4DynamicParticle(theElectron, gDirection,
577 ElectronEnergy - electron_mass_c2);
578
579 // create G4DynamicParticle object for the particle2
580 G4DynamicParticle* aParticle2 =
581 new G4DynamicParticle(thePositron, gDirection,
582 PositronEnergy - electron_mass_c2);
583
584 // primary change
585 kineticEnergy -= (ElectronEnergy + PositronEnergy);
586 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
587
588 partDirection *= totalMomentum;
589 partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
590 partDirection = partDirection.unit();
591 fParticleChange->SetProposedMomentumDirection(partDirection);
592
593 // add secondary
594 vdp->push_back(aParticle1);
595 vdp->push_back(aParticle2);
596}
#define G4UniformRand()
Definition: Randomize.hh:53
Hep3Vector unit() const
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
G4double GetZ() const
Definition: G4Element.hh:131
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)

◆ SetCurrentElement()

void G4MuPairProductionModel::SetCurrentElement ( G4double  Z)
inlineprotected

◆ SetLowestKineticEnergy()

void G4MuPairProductionModel::SetLowestKineticEnergy ( G4double  e)
inline

Definition at line 185 of file G4MuPairProductionModel.hh.

186{
187 lowestKinEnergy = e;
188}

◆ SetParticle()

void G4MuPairProductionModel::SetParticle ( const G4ParticleDefinition p)
inline

Definition at line 193 of file G4MuPairProductionModel.hh.

194{
195 if(!particle) {
196 particle = p;
198 }
199}

Referenced by G4MuPairProductionModel(), and Initialise().

Member Data Documentation

◆ currentZ

G4double G4MuPairProductionModel::currentZ
protected

Definition at line 151 of file G4MuPairProductionModel.hh.

Referenced by SampleSecondaries(), and SetCurrentElement().

◆ factorForCross

G4double G4MuPairProductionModel::factorForCross
protected

◆ lnZ

G4double G4MuPairProductionModel::lnZ
protected

Definition at line 154 of file G4MuPairProductionModel.hh.

Referenced by G4MuPairProductionModel(), and SetCurrentElement().

◆ nist

G4NistManager* G4MuPairProductionModel::nist
protected

Definition at line 146 of file G4MuPairProductionModel.hh.

Referenced by G4MuPairProductionModel(), and SetCurrentElement().

◆ particle

◆ particleMass

◆ sqrte

◆ wgi

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 156 of file G4MuPairProductionModel.hh.

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

◆ xgi

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 156 of file G4MuPairProductionModel.hh.

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

◆ z13

◆ z23


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