Geant4 11.4.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=nullptr, const G4String &nam="muPairProd")
 ~G4MuPairProductionModel () override=default
void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel) override
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double) override
virtual G4double ComputeDMicroscopicCrossSection (G4double tkin, G4double Z, G4double pairEnergy)
void SetLowestKineticEnergy (G4double e)
void SetParticle (const G4ParticleDefinition *)
G4MuPairProductionModeloperator= (const G4MuPairProductionModel &right)=delete
 G4MuPairProductionModel (const G4MuPairProductionModel &)=delete
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 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 G4Material *, const G4ParticleDefinition *, const G4double kinEnergy, const G4double cutEnergy, const G4double &length, G4double &eloss)
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 *)
virtual void FillNumberOfSecondaries (G4int &numberOfTriplets, G4int &numberOfRecoil)
virtual void ModelDescription (std::ostream &outFile) const
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)
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4ElementSelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4ElementGetCurrentElement (const G4Material *mat=nullptr) const
G4int SelectRandomAtomNumber (const G4Material *) const
const G4IsotopeGetCurrentIsotope (const G4Element *elm=nullptr) const
G4int SelectIsotopeNumber (const G4Element *) const
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
G4ElementDataGetElementData ()
G4PhysicsTableGetCrossSectionTable ()
G4VEmFluctuationModelGetModelOfFluctuations ()
G4VEmAngularDistributionGetAngularDistribution ()
G4VEmModelGetTripletModel ()
void SetTripletModel (G4VEmModel *)
void SetAngularDistribution (G4VEmAngularDistribution *)
G4double HighEnergyLimit () const
G4double LowEnergyLimit () const
G4double HighEnergyActivationLimit () const
G4double LowEnergyActivationLimit () const
G4double PolarAngleLimit () const
G4double SecondaryThreshold () 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) const
void SetPolarAngleLimit (G4double)
void SetSecondaryThreshold (G4double)
void SetDeexcitationFlag (G4bool val)
void SetForceBuildTable (G4bool val)
void SetFluctuationFlag (G4bool val)
G4bool IsMaster () const
void SetUseBaseMaterials (G4bool val)
G4bool UseBaseMaterials () const
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
const G4StringGetName () const
void SetCurrentCouple (const G4MaterialCutsCouple *)
G4bool IsLocked () const
void SetLocked (G4bool)
void SetLPMFlag (G4bool)
void SetMasterThread (G4bool)
G4VEmModeloperator= (const G4VEmModel &right)=delete
 G4VEmModel (const G4VEmModel &)=delete

Protected Member Functions

G4double ComputMuPairLoss (G4double Z, G4double tkin, G4double cut, G4double tmax)
G4double ComputeMicroscopicCrossSection (G4double tkin, G4double Z, G4double cut)
G4double FindScaledEnergy (G4int Z, G4double rand, G4double logTkin, G4double yymin, G4double yymax)
G4double MaxSecondaryEnergyForElement (G4double kineticEnergy, G4double Z)
void MakeSamplingTables ()
void StoreTables () const
G4bool RetrieveTables ()
virtual void DataCorrupted (G4int Z, G4double logTkin) const
Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
G4ParticleChangeForGammaGetParticleChangeForGamma ()
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
const G4MaterialCutsCoupleCurrentCouple () const
void SetCurrentElement (const G4Element *)

Protected Attributes

G4ParticleChangeForLossfParticleChange = nullptr
const G4ParticleDefinitionparticle = nullptr
G4NistManagernist = nullptr
G4double factorForCross
G4double sqrte
G4double particleMass = 0.0
G4double z13 = 0.0
G4double z23 = 0.0
G4double lnZ = 0.0
G4double minPairEnergy
G4double lowestKinEnergy
G4double emin
G4double emax
G4double ymin = -5.0
G4double dy = 0.005
G4int currentZ = 0
G4int nYBinPerDecade = 4
std::size_t nbiny = 1000
std::size_t nbine = 0
G4bool fTableToFile = false
Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData = nullptr
G4VParticleChangepParticleChange = nullptr
G4PhysicsTablexSectionTable = nullptr
const G4MaterialpBaseMaterial = nullptr
const std::vector< G4double > * theDensityFactor = nullptr
const std::vector< G4int > * theDensityIdx = nullptr
G4double inveplus
G4double pFactor = 1.0
std::size_t currentCoupleIndex = 0
std::size_t basedCoupleIndex = 0
G4bool lossFlucFlag = true

Static Protected Attributes

static const G4int NZDATPAIR = 5
static const G4int NINTPAIR = 8
static const G4int ZDATPAIR [NZDATPAIR] = {1, 4, 13, 29, 92}
static const G4double xgi [NINTPAIR]
static const G4double wgi [NINTPAIR]

Detailed Description

Definition at line 72 of file G4MuPairProductionModel.hh.

Constructor & Destructor Documentation

◆ G4MuPairProductionModel() [1/2]

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

Definition at line 118 of file G4MuPairProductionModel.cc.

120 : G4VEmModel(nam),
121 factorForCross(CLHEP::fine_structure_const*CLHEP::fine_structure_const*
122 CLHEP::classic_electr_radius*CLHEP::classic_electr_radius*
123 4./(3.*CLHEP::pi)),
124 sqrte(std::sqrt(G4Exp(1.))),
125 minPairEnergy(4.*CLHEP::electron_mass_c2),
126 lowestKinEnergy(0.85*CLHEP::GeV)
127{
129
130 theElectron = G4Electron::Electron();
131 thePositron = G4Positron::Positron();
132
133 if(nullptr != p) {
134 SetParticle(p);
135 lowestKinEnergy = std::max(lowestKinEnergy, p->GetPDGMass()*8.0);
136 }
138 emax = emin*10000.;
139 SetAngularDistribution(new G4ModifiedMephi());
140}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
static G4Electron * Electron()
Definition G4Electron.cc:91
void SetParticle(const G4ParticleDefinition *)
static G4NistManager * Instance()
static G4Positron * Positron()
Definition G4Positron.cc:90
G4VEmModel(const G4String &nam)
Definition G4VEmModel.cc:67
void SetAngularDistribution(G4VEmAngularDistribution *)

Referenced by G4hPairProductionModel::G4hPairProductionModel(), G4MuonToMuonPairProductionModel::G4MuonToMuonPairProductionModel(), G4MuPairProductionModel(), and operator=().

◆ ~G4MuPairProductionModel()

G4MuPairProductionModel::~G4MuPairProductionModel ( )
overridedefault

◆ G4MuPairProductionModel() [2/2]

G4MuPairProductionModel::G4MuPairProductionModel ( const G4MuPairProductionModel & )
delete

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

Reimplemented from G4VEmModel.

Definition at line 456 of file G4MuPairProductionModel.cc.

462{
463 G4double cross = 0.0;
464 if (kineticEnergy <= lowestKinEnergy) { return cross; }
465
466 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy, Z);
467 G4double tmax = std::min(maxEnergy, maxPairEnergy);
468 G4double cut = std::max(cutEnergy, minPairEnergy);
469 if (cut >= tmax) { return cross; }
470
471 cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
472 if(tmax < kineticEnergy) {
473 cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
474 }
475 return cross;
476}
double G4double
Definition G4Types.hh:83
G4double ComputeMicroscopicCrossSection(G4double tkin, G4double Z, G4double cut)
G4double MaxSecondaryEnergyForElement(G4double kineticEnergy, G4double Z)

◆ ComputeDEDXPerVolume()

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

Reimplemented from G4VEmModel.

Definition at line 228 of file G4MuPairProductionModel.cc.

233{
234 G4double dedx = 0.0;
235 if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
236 { return dedx; }
237
238 const G4ElementVector* theElementVector = material->GetElementVector();
239 const G4double* theAtomicNumDensityVector =
240 material->GetAtomicNumDensityVector();
241
242 // loop for elements in the material
243 for (std::size_t i=0; i<material->GetNumberOfElements(); ++i) {
244 G4double Z = (*theElementVector)[i]->GetZ();
245 G4double tmax = MaxSecondaryEnergyForElement(kineticEnergy, Z);
246 G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
247 dedx += loss*theAtomicNumDensityVector[i];
248 }
249 dedx = std::max(dedx, 0.0);
250 return dedx;
251}
std::vector< const G4Element * > G4ElementVector
const G4ElementVector * GetElementVector() const
const G4double * GetAtomicNumDensityVector() const
std::size_t GetNumberOfElements() const
G4double ComputMuPairLoss(G4double Z, G4double tkin, G4double cut, G4double tmax)

◆ ComputeDMicroscopicCrossSection()

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

Reimplemented in G4MuonToMuonPairProductionModel.

Definition at line 320 of file G4MuPairProductionModel.cc.

327{
328 static const G4double bbbtf= 183. ;
329 static const G4double bbbh = 202.4 ;
330 static const G4double g1tf = 1.95e-5 ;
331 static const G4double g2tf = 5.3e-5 ;
332 static const G4double g1h = 4.4e-5 ;
333 static const G4double g2h = 4.8e-5 ;
334
335 if (pairEnergy <= minPairEnergy)
336 return 0.0;
337
338 G4double totalEnergy = tkin + particleMass;
339 G4double residEnergy = totalEnergy - pairEnergy;
340
341 if (residEnergy <= 0.75*sqrte*z13*particleMass)
342 return 0.0;
343
344 G4double a0 = 1.0 / (totalEnergy * residEnergy);
345 G4double alf = 4.0 * electron_mass_c2 / pairEnergy;
346 G4double rt = std::sqrt(1.0 - alf);
347 G4double delta = 6.0 * particleMass * particleMass * a0;
348 G4double tmnexp = alf/(1.0 + rt) + delta*rt;
349
350 if(tmnexp >= 1.0) { return 0.0; }
351
352 G4double tmn = G4Log(tmnexp);
353
354 G4double massratio = particleMass/CLHEP::electron_mass_c2;
355 G4double massratio2 = massratio*massratio;
356 G4double inv_massratio2 = 1.0 / massratio2;
357
358 // zeta calculation
359 G4double bbb,g1,g2;
360 if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
361 else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
362
363 G4double zeta = 0.0;
364 G4double z1exp = totalEnergy / (particleMass + g1*z23*totalEnergy);
365
366 // 35.221047195922 is the root of zeta1(x) = 0.073 * log(x) - 0.26, so the
367 // condition below is the same as zeta1 > 0.0, but without calling log(x)
368 if (z1exp > 35.221047195922)
369 {
370 G4double z2exp = totalEnergy / (particleMass + g2*z13*totalEnergy);
371 zeta = (0.073 * G4Log(z1exp) - 0.26) / (0.058 * G4Log(z2exp) - 0.14);
372 }
373
374 G4double z2 = Z*(Z+zeta);
375 G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
376 G4double beta = 0.5*pairEnergy*pairEnergy*a0;
377 G4double xi0 = 0.5*massratio2*beta;
378
379 // Gaussian integration in ln(1-ro) ( with 8 points)
380 G4double rho[NINTPAIR];
381 G4double rho2[NINTPAIR];
382 G4double xi[NINTPAIR];
383 G4double xi1[NINTPAIR];
384 G4double xii[NINTPAIR];
385
386 for (G4int i = 0; i < NINTPAIR; ++i)
387 {
388 rho[i] = G4Exp(tmn*xgi[i]) - 1.0; // rho = -asymmetry
389 rho2[i] = rho[i] * rho[i];
390 xi[i] = xi0*(1.0-rho2[i]);
391 xi1[i] = 1.0 + xi[i];
392 xii[i] = 1.0 / xi[i];
393 }
394
395 G4double ye1[NINTPAIR];
396 G4double ym1[NINTPAIR];
397
398 G4double b40 = 4.0 * beta;
399 G4double b62 = 6.0 * beta + 2.0;
400
401 for (G4int i = 0; i < NINTPAIR; ++i)
402 {
403 G4double yeu = (b40 + 5.0) + (b40 - 1.0) * rho2[i];
404 G4double yed = b62*G4Log(3.0 + xii[i]) + (2.0 * beta - 1.0)*rho2[i] - b40;
405
406 G4double ymu = b62 * (1.0 + rho2[i]) + 6.0;
407 G4double ymd = (b40 + 3.0)*(1.0 + rho2[i])*G4Log(3.0 + xi[i])
408 + 2.0 - 3.0 * rho2[i];
409
410 ye1[i] = 1.0 + yeu / yed;
411 ym1[i] = 1.0 + ymu / ymd;
412 }
413
414 G4double be[NINTPAIR];
415 G4double bm[NINTPAIR];
416
417 for(G4int i = 0; i < NINTPAIR; ++i) {
418 if(xi[i] <= 1000.0) {
419 be[i] = ((2.0 + rho2[i])*(1.0 + beta) +
420 xi[i]*(3.0 + rho2[i]))*G4Log(1.0 + xii[i]) +
421 (1.0 - rho2[i] - beta)/xi1[i] - (3.0 + rho2[i]);
422 } else {
423 be[i] = 0.5*(3.0 - rho2[i] + 2.0*beta*(1.0 + rho2[i]))*xii[i];
424 }
425
426 if(xi[i] >= 0.001) {
427 G4double a10 = (1.0 + 2.0 * beta) * (1.0 - rho2[i]);
428 bm[i] = ((1.0 + rho2[i])*(1.0 + 1.5 * beta) - a10*xii[i])*G4Log(xi1[i]) +
429 xi[i] * (1.0 - rho2[i] - beta)/xi1[i] + a10;
430 } else {
431 bm[i] = 0.5*(5.0 - rho2[i] + beta * (3.0 + rho2[i]))*xi[i];
432 }
433 }
434
435 G4double sum = 0.0;
436
437 for (G4int i = 0; i < NINTPAIR; ++i) {
438 G4double screen = screen0*xi1[i]/(1.0 - rho2[i]);
439 G4double ale = G4Log(bbb/z13*std::sqrt(xi1[i]*ye1[i])/(1. + screen*ye1[i]));
440 G4double cre = 0.5*G4Log(1. + 2.25*z23*xi1[i]*ye1[i]*inv_massratio2);
441
442 G4double fe = (ale-cre)*be[i];
443 fe = std::max(fe, 0.0);
444
445 G4double alm_crm = G4Log(bbb*massratio/(1.5*z23*(1. + screen*ym1[i])));
446 G4double fm = std::max(alm_crm*bm[i], 0.0)*inv_massratio2;
447
448 sum += wgi[i]*(1.0 + rho[i])*(fe + fm);
449 }
450
451 return -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
452}
const G4double a0
G4fissionEvent * fe
G4double G4Log(G4double x)
Definition G4Log.hh:169
int G4int
Definition G4Types.hh:85
static const G4double xgi[NINTPAIR]
static const G4double wgi[NINTPAIR]

Referenced by ComputeMicroscopicCrossSection(), ComputMuPairLoss(), and MakeSamplingTables().

◆ ComputeMicroscopicCrossSection()

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

Definition at line 288 of file G4MuPairProductionModel.cc.

292{
293 G4double cross = 0.;
295 G4double cut = std::max(cutEnergy, minPairEnergy);
296 if (tmax <= cut) { return cross; }
297
298 G4double aaa = G4Log(cut);
299 G4double bbb = G4Log(tmax);
300 G4int kkk = std::min(std::max(G4lrint((bbb-aaa)/ak1 + ak2), 8), 1);
301
302 G4double hhh = (bbb-aaa)/(kkk);
303 G4double x = aaa;
304
305 for (G4int l=0; l<kkk; ++l) {
306 for (G4int i=0; i<NINTPAIR; ++i) {
307 G4double ep = G4Exp(x + xgi[i]*hhh);
308 cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
309 }
310 x += hhh;
311 }
312
313 cross *= hhh;
314 cross = std::max(cross, 0.0);
315 return cross;
316}
virtual G4double ComputeDMicroscopicCrossSection(G4double tkin, G4double Z, G4double pairEnergy)
int G4lrint(double ad)
Definition templates.hh:134

Referenced by ComputeCrossSectionPerAtom().

◆ ComputMuPairLoss()

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

Definition at line 255 of file G4MuPairProductionModel.cc.

259{
260 G4double loss = 0.0;
261
262 G4double cut = std::min(cutEnergy, tmax);
263 if(cut <= minPairEnergy) { return loss; }
264
265 // calculate the rectricted loss
266 // numerical integration in log(PairEnergy)
268 G4double bbb = G4Log(cut);
269
270 G4int kkk = std::min(std::max(G4lrint((bbb-aaa)/ak1 + ak2), 8), 1);
271 G4double hhh = (bbb-aaa)/kkk;
272 G4double x = aaa;
273
274 for (G4int l=0 ; l<kkk; ++l) {
275 for (G4int ll=0; ll<NINTPAIR; ++ll) {
276 G4double ep = G4Exp(x+xgi[ll]*hhh);
277 loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
278 }
279 x += hhh;
280 }
281 loss *= hhh;
282 loss = std::max(loss, 0.0);
283 return loss;
284}

Referenced by ComputeDEDXPerVolume().

◆ DataCorrupted()

void G4MuPairProductionModel::DataCorrupted ( G4int Z,
G4double logTkin ) const
protectedvirtual

Reimplemented in G4MuonToMuonPairProductionModel.

Definition at line 701 of file G4MuPairProductionModel.cc.

702{
704 ed << "G4ElementData is not properly initialized Z= " << Z
705 << " Ekin(MeV)= " << G4Exp(logTkin)
706 << " IsMasterThread= " << IsMaster()
707 << " Model " << GetName();
708 G4Exception("G4MuPairProductionModel::()", "em0033", FatalException, ed, "");
709}
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
G4bool IsMaster() const
const G4String & GetName() const

Referenced by FindScaledEnergy(), and StoreTables().

◆ FindScaledEnergy()

G4double G4MuPairProductionModel::FindScaledEnergy ( G4int Z,
G4double rand,
G4double logTkin,
G4double yymin,
G4double yymax )
protected

Definition at line 681 of file G4MuPairProductionModel.cc.

684{
685 G4double res = yymin;
686 G4Physics2DVector* pv = fElementData->GetElement2DData(iz);
687 if (nullptr != pv) {
688 G4double pmin = pv->Value(yymin, logTkin);
689 G4double pmax = pv->Value(yymax, logTkin);
690 G4double p0 = pv->Value(0.0, logTkin);
691 if(p0 <= 0.0) { DataCorrupted(ZDATPAIR[iz], logTkin); }
692 else { res = pv->FindLinearX((pmin + rand*(pmax - pmin))/p0, logTkin); }
693 } else {
694 DataCorrupted(ZDATPAIR[iz], logTkin);
695 }
696 return res;
697}
virtual void DataCorrupted(G4int Z, G4double logTkin) const
static const G4int ZDATPAIR[NZDATPAIR]
G4double Value(G4double x, G4double y, std::size_t &lastidx, std::size_t &lastidy) const
G4double FindLinearX(G4double rand, G4double y, std::size_t &lastidy) const
G4ElementData * fElementData

Referenced by G4MuonToMuonPairProductionModel::SampleSecondaries(), and SampleSecondaries().

◆ Initialise()

void G4MuPairProductionModel::Initialise ( const G4ParticleDefinition * p,
const G4DataVector & cuts )
overridevirtual

Implements G4VEmModel.

Definition at line 153 of file G4MuPairProductionModel.cc.

155{
156 SetParticle(p);
157
158 if (nullptr == fParticleChange) {
160
161 // define scale of internal table for each thread only once
162 if (0 == nbine) {
163 emin = std::max(lowestKinEnergy, LowEnergyLimit());
164 emax = std::max(HighEnergyLimit(), emin*2);
165 nbine = std::size_t(nYBinPerDecade*std::log10(emax/emin));
166 if(nbine < 3) { nbine = 3; }
167
169 dy = -ymin/G4double(nbiny);
170 }
171 if (p == particle) {
172 G4int pdg = std::abs(p->GetPDGEncoding());
173 if (pdg == 2212) {
174 dataName = "pEEPairProd";
175 } else if (pdg == 321) {
176 dataName = "kaonEEPairProd";
177 } else if (pdg == 211) {
178 dataName = "pionEEPairProd";
179 } else if (pdg == 11) {
180 dataName = "eEEPairProd";
181 } else if (pdg == 13) {
182 if (GetName() == "muToMuonPairProd") {
183 dataName = "muMuMuPairProd";
184 } else {
185 dataName = "muEEPairProd";
186 }
187 }
188 }
189 }
190
191 // for low-energy application this process should not work
192 if(lowestKinEnergy >= HighEnergyLimit()) { return; }
193
194 if (p == particle) {
196 fElementData = reg->GetElementDataByName(dataName);
197 if (nullptr == fElementData) {
198 G4AutoLock l(&theMuPairMutex);
199 fElementData = reg->NewElementData(dataName, NZDATPAIR);
201 if (useDataFile) {
202 useDataFile = RetrieveTables();
203 }
204 else {
206 }
207 if (fTableToFile) { StoreTables(); }
208 l.unlock();
209 }
210 if (IsMaster()) {
212 }
213 }
214}
G4TemplateAutoLock< G4Mutex > G4AutoLock
bool G4bool
Definition G4Types.hh:86
static G4ElementDataRegistry * Instance()
static G4EmParameters * Instance()
G4bool RetrieveMuDataFromFile() const
const G4ParticleDefinition * particle
G4ParticleChangeForLoss * fParticleChange
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
G4ParticleChangeForLoss * GetParticleChangeForLoss()

◆ InitialiseLocal()

void G4MuPairProductionModel::InitialiseLocal ( const G4ParticleDefinition * p,
G4VEmModel * masterModel )
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 218 of file G4MuPairProductionModel.cc.

220{
223 }
224}
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
std::vector< G4EmElementSelector * > * GetElementSelectors()

◆ MakeSamplingTables()

void G4MuPairProductionModel::MakeSamplingTables ( )
protected

Definition at line 480 of file G4MuPairProductionModel.cc.

481{
483
484 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
485
486 G4double Z = ZDATPAIR[iz];
487 G4Physics2DVector* pv = fElementData->New2DVector(iz, G4int(nbiny)+1, G4int(nbine)+1);
488 G4double kinEnergy = emin;
489
490 for (std::size_t it=0; it<=nbine; ++it) {
491
492 pv->PutY(it, G4Log(kinEnergy/CLHEP::MeV));
493 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy, Z);
494 /*
495 G4cout << "it= " << it << " E= " << kinEnergy
496 << " " << particle->GetParticleName()
497 << " maxE= " << maxPairEnergy << " minE= " << minPairEnergy
498 << " ymin= " << ymin << G4endl;
499 */
500 G4double coef = G4Log(minPairEnergy/kinEnergy)/ymin;
501 G4double ymax = G4Log(maxPairEnergy/kinEnergy)/coef;
502 G4double fac = (ymax - ymin)/dy;
503 std::size_t imax = (std::size_t)fac;
504 fac -= (G4double)imax;
505
506 G4double xSec = 0.0;
507 G4double x = ymin;
508 /*
509 G4cout << "Z= " << currentZ << " z13= " << z13
510 << " mE= " << maxPairEnergy << " ymin= " << ymin
511 << " dy= " << dy << " c= " << coef << G4endl;
512 */
513 // start from zero
514 pv->PutValue(0, it, 0.0);
515 if(0 == it) { pv->PutX(nbiny, 0.0); }
516
517 for (std::size_t i=0; i<nbiny; ++i) {
518
519 if(0 == it) { pv->PutX(i, x); }
520
521 if(i < imax) {
522 G4double ep = kinEnergy*G4Exp(coef*(x + dy*0.5));
523
524 // not multiplied by interval, because table
525 // will be used only for sampling
526 //G4cout << "i= " << i << " x= " << x << "E= " << kinEnergy
527 // << " Egamma= " << ep << G4endl;
528 xSec += ep*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
529
530 // last bin before the kinematic limit
531 } else if(i == imax) {
532 G4double ep = kinEnergy*G4Exp(coef*(x + fac*dy*0.5));
533 xSec += ep*fac*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
534 }
535 pv->PutValue(i + 1, it, xSec);
536 x += dy;
537 }
538 kinEnergy *= factore;
539
540 // to avoid precision lost
541 if(it+1 == nbine) { kinEnergy = emax; }
542 }
543 }
544}
void PutY(std::size_t idy, G4double value)
void PutValue(std::size_t idx, std::size_t idy, G4double value)
void PutX(std::size_t idx, G4double value)

Referenced by Initialise().

◆ MaxSecondaryEnergyForElement()

G4double G4MuPairProductionModel::MaxSecondaryEnergyForElement ( G4double kineticEnergy,
G4double Z )
inlineprotected

Definition at line 204 of file G4MuPairProductionModel.hh.

206{
207 G4int Z = G4lrint(ZZ);
208 if(Z != currentZ) {
209 currentZ = Z;
210 z13 = nist->GetZ13(Z);
211 z23 = z13*z13;
212 lnZ = nist->GetLOGZ(Z);
213 }
214 return kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
215}

Referenced by ComputeCrossSectionPerAtom(), ComputeDEDXPerVolume(), ComputeMicroscopicCrossSection(), MakeSamplingTables(), G4MuonToMuonPairProductionModel::SampleSecondaries(), and SampleSecondaries().

◆ MinPrimaryEnergy()

G4double G4MuPairProductionModel::MinPrimaryEnergy ( const G4Material * ,
const G4ParticleDefinition * ,
G4double cut )
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 144 of file G4MuPairProductionModel.cc.

147{
148 return std::max(lowestKinEnergy, cut);
149}

◆ operator=()

G4MuPairProductionModel & G4MuPairProductionModel::operator= ( const G4MuPairProductionModel & right)
delete

◆ RetrieveTables()

G4bool G4MuPairProductionModel::RetrieveTables ( )
protected

Definition at line 731 of file G4MuPairProductionModel.cc.

732{
733 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
734 G4double Z = ZDATPAIR[iz];
735 G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
736 std::ostringstream ss;
737 ss << G4EmParameters::Instance()->GetDirLEDATA() << "/mupair/"
738 << particle->GetParticleName() << Z << ".dat";
739 std::ifstream infile(ss.str(), std::ios::in);
740 if(!pv->Retrieve(infile)) {
741 delete pv;
742 return false;
743 }
744 fElementData->InitialiseForElement(iz, pv);
745 }
746 return true;
747}
const G4String & GetDirLEDATA() const
G4bool Retrieve(std::ifstream &fIn)

Referenced by Initialise().

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 548 of file G4MuPairProductionModel.cc.

554{
555 G4double kinEnergy = aDynamicParticle->GetKineticEnergy();
556 //G4cout << "------- G4MuPairProductionModel::SampleSecondaries E(MeV)= "
557 // << kinEnergy << " "
558 // << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
559 G4double totalEnergy = kinEnergy + particleMass;
560 G4double totalMomentum =
561 std::sqrt(kinEnergy*(kinEnergy + 2.0*particleMass));
562
563 G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
564
565 // select randomly one element constituing the material
566 const G4Element* anElement = SelectRandomAtom(couple,particle,kinEnergy);
567
568 // define interval of energy transfer
569 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy,
570 anElement->GetZ());
571 G4double maxEnergy = std::min(tmax, maxPairEnergy);
572 G4double minEnergy = std::max(tmin, minPairEnergy);
573
574 if (minEnergy >= maxEnergy) { return; }
575 //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
576 // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
577 // << " ymin= " << ymin << " dy= " << dy << G4endl;
578
580
581 // compute limits
582 G4double yymin = G4Log(minEnergy/kinEnergy)/coeff;
583 G4double yymax = G4Log(maxEnergy/kinEnergy)/coeff;
584
585 //G4cout << "yymin= " << yymin << " yymax= " << yymax << G4endl;
586
587 // units should not be used, bacause table was built without
588 G4double logTkin = G4Log(kinEnergy/CLHEP::MeV);
589
590 // sample e-e+ energy, pair energy first
591
592 // select sample table via Z
593 G4int iz1(0), iz2(0);
594 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
595 if(currentZ == ZDATPAIR[iz]) {
596 iz1 = iz2 = iz;
597 break;
598 } else if(currentZ < ZDATPAIR[iz]) {
599 iz2 = iz;
600 if(iz > 0) { iz1 = iz-1; }
601 else { iz1 = iz2; }
602 break;
603 }
604 }
605 if (0 == iz1) { iz1 = iz2 = NZDATPAIR-1; }
606
607 G4double pairEnergy = 0.0;
608 G4int count = 0;
609 //G4cout << "start loop Z1= " << iz1 << " Z2= " << iz2 << G4endl;
610 do {
611 ++count;
612 // sampling using only one random number
613 G4double rand = G4UniformRand();
614
615 G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
616 if(iz1 != iz2) {
617 G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
618 G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1]);
619 G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2]);
620 //G4cout << count << ". x= " << x << " x2= " << x2
621 // << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
622 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
623 }
624 //G4cout << "x= " << x << " coeff= " << coeff << G4endl;
625 pairEnergy = kinEnergy*G4Exp(x*coeff);
626
627 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
628 } while((pairEnergy < minEnergy || pairEnergy > maxEnergy) && 10 > count);
629
630 //G4cout << "## pairEnergy(GeV)= " << pairEnergy/GeV
631 // << " Etot(GeV)= " << totalEnergy/GeV << G4endl;
632
633 // sample r=(E+-E-)/pairEnergy ( uniformly .....)
634 G4double rmax =
635 (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-pairEnergy)))
636 *std::sqrt(1.-minPairEnergy/pairEnergy);
637 G4double r = rmax * (-1.+2.*G4UniformRand()) ;
638
639 // compute energies from pairEnergy,r
640 G4double eEnergy = (1.-r)*pairEnergy*0.5;
641 G4double pEnergy = pairEnergy - eEnergy;
642
643 // Sample angles
644 G4ThreeVector eDirection, pDirection;
645 //
646 GetAngularDistribution()->SamplePairDirections(aDynamicParticle,
647 eEnergy, pEnergy,
648 eDirection, pDirection);
649 // create G4DynamicParticle object for e+e-
650 eEnergy = std::max(eEnergy - CLHEP::electron_mass_c2, 0.0);
651 pEnergy = std::max(pEnergy - CLHEP::electron_mass_c2, 0.0);
652 auto aParticle1 = new G4DynamicParticle(theElectron,eDirection,eEnergy);
653 auto aParticle2 = new G4DynamicParticle(thePositron,pDirection,pEnergy);
654 // Fill output vector
655 vdp->push_back(aParticle1);
656 vdp->push_back(aParticle2);
657
658 // primary change
659 kinEnergy -= pairEnergy;
660 partDirection *= totalMomentum;
661 partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
662 partDirection = partDirection.unit();
663
664 // if energy transfer is higher than threshold (very high by default)
665 // then stop tracking the primary particle and create a new secondary
666 if (pairEnergy > SecondaryThreshold()) {
667 fParticleChange->ProposeTrackStatus(fStopAndKill);
668 fParticleChange->SetProposedKineticEnergy(0.0);
669 auto newdp = new G4DynamicParticle(particle, partDirection, kinEnergy);
670 vdp->push_back(newdp);
671 } else { // continue tracking the primary e-/e+ otherwise
672 fParticleChange->SetProposedMomentumDirection(partDirection);
673 fParticleChange->SetProposedKineticEnergy(kinEnergy);
674 }
675 //G4cout << "-- G4MuPairProductionModel::SampleSecondaries done" << G4endl;
676}
CLHEP::Hep3Vector G4ThreeVector
@ fStopAndKill
#define G4UniformRand()
Definition Randomize.hh:52
Hep3Vector unit() const
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
G4double GetZ() const
Definition G4Element.hh:119
G4double FindScaledEnergy(G4int Z, G4double rand, G4double logTkin, G4double yymin, G4double yymax)
virtual void SamplePairDirections(const G4DynamicParticle *dp, G4double elecKinEnergy, G4double posiKinEnergy, G4ThreeVector &dirElectron, G4ThreeVector &dirPositron, G4int Z=0, const G4Material *mat=nullptr)
G4VEmAngularDistribution * GetAngularDistribution()
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double SecondaryThreshold() const

◆ SetLowestKineticEnergy()

void G4MuPairProductionModel::SetLowestKineticEnergy ( G4double e)
inline

Definition at line 185 of file G4MuPairProductionModel.hh.

186{
187 lowestKinEnergy = e;
188}

Referenced by G4ePairProduction::InitialiseEnergyLossProcess().

◆ SetParticle()

void G4MuPairProductionModel::SetParticle ( const G4ParticleDefinition * p)
inline

Definition at line 193 of file G4MuPairProductionModel.hh.

194{
195 if(nullptr == particle) {
196 particle = p;
197 particleMass = particle->GetPDGMass();
198 }
199}

Referenced by G4MuPairProductionModel(), and Initialise().

◆ StoreTables()

void G4MuPairProductionModel::StoreTables ( ) const
protected

Definition at line 713 of file G4MuPairProductionModel.cc.

714{
715 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
716 G4int Z = ZDATPAIR[iz];
717 G4Physics2DVector* pv = fElementData->GetElement2DData(Z);
718 if(nullptr == pv) {
719 DataCorrupted(Z, 1.0);
720 return;
721 }
722 std::ostringstream ss;
723 ss << "mupair/" << particle->GetParticleName() << Z << ".dat";
724 std::ofstream outfile(ss.str());
725 pv->Store(outfile);
726 }
727}
void Store(std::ofstream &fOut) const

Referenced by Initialise().

Member Data Documentation

◆ currentZ

G4int G4MuPairProductionModel::currentZ = 0
protected

◆ dy

G4double G4MuPairProductionModel::dy = 0.005
protected

Definition at line 160 of file G4MuPairProductionModel.hh.

Referenced by Initialise(), and MakeSamplingTables().

◆ emax

G4double G4MuPairProductionModel::emax
protected

◆ emin

G4double G4MuPairProductionModel::emin
protected

◆ factorForCross

G4double G4MuPairProductionModel::factorForCross
protected

◆ fParticleChange

G4ParticleChangeForLoss* G4MuPairProductionModel::fParticleChange = nullptr
protected

◆ fTableToFile

G4bool G4MuPairProductionModel::fTableToFile = false
protected

Definition at line 167 of file G4MuPairProductionModel.hh.

Referenced by Initialise().

◆ lnZ

G4double G4MuPairProductionModel::lnZ = 0.0
protected

◆ lowestKinEnergy

◆ minPairEnergy

◆ nbine

std::size_t G4MuPairProductionModel::nbine = 0
protected

Definition at line 165 of file G4MuPairProductionModel.hh.

Referenced by Initialise(), MakeSamplingTables(), and RetrieveTables().

◆ nbiny

std::size_t G4MuPairProductionModel::nbiny = 1000
protected

Definition at line 164 of file G4MuPairProductionModel.hh.

Referenced by Initialise(), MakeSamplingTables(), and RetrieveTables().

◆ NINTPAIR

◆ nist

◆ nYBinPerDecade

G4int G4MuPairProductionModel::nYBinPerDecade = 4
protected

Definition at line 163 of file G4MuPairProductionModel.hh.

Referenced by Initialise().

◆ NZDATPAIR

const G4int G4MuPairProductionModel::NZDATPAIR = 5
staticprotected

◆ particle

◆ particleMass

◆ sqrte

G4double G4MuPairProductionModel::sqrte
protected

◆ wgi

const G4double G4MuPairProductionModel::wgi
staticprotected
Initial value:
= {
0.0506142681451880, 0.1111905172266870, 0.1568533229389435, 0.1813418916891810,
0.1813418916891810, 0.1568533229389435, 0.1111905172266870, 0.0506142681451880
}

Definition at line 103 of file G4MuPairProductionModel.hh.

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

◆ xgi

const G4double G4MuPairProductionModel::xgi
staticprotected
Initial value:
= {
0.0198550717512320, 0.1016667612931865, 0.2372337950418355, 0.4082826787521750,
0.5917173212478250, 0.7627662049581645, 0.8983332387068135, 0.9801449282487680
}

Definition at line 98 of file G4MuPairProductionModel.hh.

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

◆ ymin

G4double G4MuPairProductionModel::ymin = -5.0
protected

◆ z13

G4double G4MuPairProductionModel::z13 = 0.0
protected

◆ z23

G4double G4MuPairProductionModel::z23 = 0.0
protected

◆ ZDATPAIR

const G4int G4MuPairProductionModel::ZDATPAIR = {1, 4, 13, 29, 92}
staticprotected

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