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

#include <G4RiGeMuPairProductionModel.hh>

Inheritance diagram for G4RiGeMuPairProductionModel:

Public Member Functions

 G4RiGeMuPairProductionModel (const G4ParticleDefinition *p=nullptr)
 ~G4RiGeMuPairProductionModel () 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
G4double ComputeDMicroscopicCrossSection (G4double tkin, G4double Z, G4double pairEnergy)
void SetLowestKineticEnergy (G4double e)
void SetParticle (const G4ParticleDefinition *)
G4RiGeMuPairProductionModeloperator= (const G4RiGeMuPairProductionModel &right)=delete
 G4RiGeMuPairProductionModel (const G4RiGeMuPairProductionModel &)=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
G4double randNumbs [9]
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 56 of file G4RiGeMuPairProductionModel.hh.

Constructor & Destructor Documentation

◆ G4RiGeMuPairProductionModel() [1/2]

G4RiGeMuPairProductionModel::G4RiGeMuPairProductionModel ( const G4ParticleDefinition * p = nullptr)
explicit

Definition at line 94 of file G4RiGeMuPairProductionModel.cc.

95 : G4VEmModel("muPairProdRiGe"),
96 factorForCross(CLHEP::fine_structure_const*CLHEP::fine_structure_const*
97 CLHEP::classic_electr_radius*CLHEP::classic_electr_radius*
98 4./(3.*CLHEP::pi)),
99 sqrte(std::sqrt(G4Exp(1.))),
100 minPairEnergy(4.*CLHEP::electron_mass_c2),
101 lowestKinEnergy(0.85*CLHEP::GeV)
102{
104
105 theElectron = G4Electron::Electron();
106 thePositron = G4Positron::Positron();
107
108 if (nullptr != p) {
109 SetParticle(p);
110 lowestKinEnergy = std::max(lowestKinEnergy, p->GetPDGMass()*8.0);
111 }
113 emax = emin*10000.;
114 fAngularGenerator = new G4RiGeAngularGenerator();
115 SetAngularDistribution(fAngularGenerator);
116 for (G4int i=0; i<9; ++i) { randNumbs[i] = 0.0; }
117}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
int G4int
Definition G4Types.hh:85
static G4Electron * Electron()
Definition G4Electron.cc:91
static G4NistManager * Instance()
static G4Positron * Positron()
Definition G4Positron.cc:90
void SetParticle(const G4ParticleDefinition *)
G4VEmModel(const G4String &nam)
Definition G4VEmModel.cc:67
void SetAngularDistribution(G4VEmAngularDistribution *)

Referenced by G4RiGeMuPairProductionModel(), and operator=().

◆ ~G4RiGeMuPairProductionModel()

G4RiGeMuPairProductionModel::~G4RiGeMuPairProductionModel ( )
overridedefault

◆ G4RiGeMuPairProductionModel() [2/2]

G4RiGeMuPairProductionModel::G4RiGeMuPairProductionModel ( const G4RiGeMuPairProductionModel & )
delete

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

Reimplemented from G4VEmModel.

Definition at line 434 of file G4RiGeMuPairProductionModel.cc.

439{
440 G4double cross = 0.0;
441 if (kineticEnergy <= lowestKinEnergy) { return cross; }
442
443 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy, Z);
444 G4double tmax = std::min(maxEnergy, maxPairEnergy);
445 G4double cut = std::max(cutEnergy, minPairEnergy);
446 if (cut >= tmax) { return cross; }
447
448 cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
449 if(tmax < kineticEnergy) {
450 cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
451 }
452 return cross;
453}
double G4double
Definition G4Types.hh:83
G4double MaxSecondaryEnergyForElement(G4double kineticEnergy, G4double Z)
G4double ComputeMicroscopicCrossSection(G4double tkin, G4double Z, G4double cut)

◆ ComputeDEDXPerVolume()

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

Reimplemented from G4VEmModel.

Definition at line 207 of file G4RiGeMuPairProductionModel.cc.

211{
212 G4double dedx = 0.0;
213 if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
214 { return dedx; }
215
216 const G4ElementVector* theElementVector = material->GetElementVector();
217 const G4double* theAtomicNumDensityVector =
218 material->GetAtomicNumDensityVector();
219
220 // loop for elements in the material
221 for (std::size_t i=0; i<material->GetNumberOfElements(); ++i) {
222 G4double Z = (*theElementVector)[i]->GetZ();
223 G4double tmax = MaxSecondaryEnergyForElement(kineticEnergy, Z);
224 G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
225 dedx += loss*theAtomicNumDensityVector[i];
226 }
227 dedx = std::max(dedx, 0.0);
228 return dedx;
229}
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 G4RiGeMuPairProductionModel::ComputeDMicroscopicCrossSection ( G4double tkin,
G4double Z,
G4double pairEnergy )

Definition at line 297 of file G4RiGeMuPairProductionModel.cc.

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

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

◆ ComputeMicroscopicCrossSection()

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

Definition at line 266 of file G4RiGeMuPairProductionModel.cc.

269{
270 G4double cross = 0.;
272 G4double cut = std::max(cutEnergy, minPairEnergy);
273 if (tmax <= cut) { return cross; }
274
275 G4double aaa = G4Log(cut);
276 G4double bbb = G4Log(tmax);
277 G4int kkk = std::min(std::max(G4lrint((bbb-aaa)/ak1 + ak2), 8), 1);
278
279 G4double hhh = (bbb-aaa)/(kkk);
280 G4double x = aaa;
281
282 for (G4int l=0; l<kkk; ++l) {
283 for (G4int i=0; i<NINTPAIR; ++i) {
284 G4double ep = G4Exp(x + xgi[i]*hhh);
285 cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
286 }
287 x += hhh;
288 }
289
290 cross *= hhh;
291 cross = std::max(cross, 0.0);
292 return cross;
293}
G4double ComputeDMicroscopicCrossSection(G4double tkin, G4double Z, G4double pairEnergy)
int G4lrint(double ad)
Definition templates.hh:134

Referenced by ComputeCrossSectionPerAtom().

◆ ComputMuPairLoss()

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

Definition at line 233 of file G4RiGeMuPairProductionModel.cc.

236{
237 G4double loss = 0.0;
238
239 G4double cut = std::min(cutEnergy, tmax);
240 if(cut <= minPairEnergy) { return loss; }
241
242 // calculate the rectricted loss
243 // numerical integration in log(PairEnergy)
245 G4double bbb = G4Log(cut);
246
247 G4int kkk = std::min(std::max(G4lrint((bbb-aaa)/ak1 + ak2), 8), 1);
248 G4double hhh = (bbb-aaa)/kkk;
249 G4double x = aaa;
250
251 for (G4int l=0 ; l<kkk; ++l) {
252 for (G4int ll=0; ll<NINTPAIR; ++ll) {
253 G4double ep = G4Exp(x+xgi[ll]*hhh);
254 loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
255 }
256 x += hhh;
257 }
258 loss *= hhh;
259 loss = std::max(loss, 0.0);
260 return loss;
261}

Referenced by ComputeDEDXPerVolume().

◆ DataCorrupted()

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

Definition at line 951 of file G4RiGeMuPairProductionModel.cc.

952{
954 ed << "G4ElementData is not properly initialized Z= " << Z
955 << " Ekin(MeV)= " << G4Exp(logTkin)
956 << " IsMasterThread= " << IsMaster()
957 << " Model " << GetName();
958 G4Exception("G4RiGeMuPairProductionModel::()", "em0033", FatalException, ed, "");
959}
@ 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 G4RiGeMuPairProductionModel::FindScaledEnergy ( G4int Z,
G4double rand,
G4double logTkin,
G4double yymin,
G4double yymax )
protected

Definition at line 931 of file G4RiGeMuPairProductionModel.cc.

934{
935 G4double res = yymin;
936 G4Physics2DVector* pv = fElementData->GetElement2DData(iz);
937 if (nullptr != pv) {
938 G4double pmin = pv->Value(yymin, logTkin);
939 G4double pmax = pv->Value(yymax, logTkin);
940 G4double p0 = pv->Value(0.0, logTkin);
941 if(p0 <= 0.0) { DataCorrupted(ZDATPAIR[iz], logTkin); }
942 else { res = pv->FindLinearX((pmin + rand*(pmax - pmin))/p0, logTkin); }
943 } else {
944 DataCorrupted(ZDATPAIR[iz], logTkin);
945 }
946 return res;
947}
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
static const G4int ZDATPAIR[NZDATPAIR]
virtual void DataCorrupted(G4int Z, G4double logTkin) const
G4ElementData * fElementData

Referenced by SampleSecondaries().

◆ Initialise()

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

Implements G4VEmModel.

Definition at line 131 of file G4RiGeMuPairProductionModel.cc.

133{
134 SetParticle(p);
135
136 if (nullptr == fParticleChange) {
138
139 // define scale of internal table for each thread only once
140 if (0 == nbine) {
141 emin = std::max(lowestKinEnergy, LowEnergyLimit());
142 emax = std::max(HighEnergyLimit(), emin*2);
143 nbine = std::size_t(nYBinPerDecade*std::log10(emax/emin));
144 if(nbine < 3) { nbine = 3; }
145
147 dy = -ymin/G4double(nbiny);
148 }
149 if (p == particle) {
150 G4int pdg = std::abs(p->GetPDGEncoding());
151 if (pdg == 2212) {
152 dataName = "pEEPairProd";
153 } else if (pdg == 321) {
154 dataName = "kaonEEPairProd";
155 } else if (pdg == 211) {
156 dataName = "pionEEPairProd";
157 } else if (pdg == 11) {
158 dataName = "eEEPairProd";
159 } else if (pdg == 13) {
160 if (GetName() == "muToMuonPairProd") {
161 dataName = "muMuMuPairProd";
162 } else {
163 dataName = "muEEPairProd";
164 }
165 }
166 }
167 }
168
169 // for low-energy application this process should not work
170 if(lowestKinEnergy >= HighEnergyLimit()) { return; }
171
172 if (p == particle) {
174 fElementData = data->GetElementDataByName(dataName);
175 if (nullptr == fElementData) {
176 G4AutoLock l(&theRiGeMuPairMutex);
177 fElementData = data->GetElementDataByName(dataName);
178 if (nullptr == fElementData) {
179 fElementData = new G4ElementData(NZDATPAIR);
180 fElementData->SetName(dataName);
181 }
183 if (useDataFile) { useDataFile = RetrieveTables(); }
184 if (!useDataFile) { MakeSamplingTables(); }
185 if (fTableToFile) { StoreTables(); }
186 l.unlock();
187 }
188 if (IsMaster()) {
190 }
191 }
192}
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 G4RiGeMuPairProductionModel::InitialiseLocal ( const G4ParticleDefinition * p,
G4VEmModel * masterModel )
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 196 of file G4RiGeMuPairProductionModel.cc.

198{
201 }
202}
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
std::vector< G4EmElementSelector * > * GetElementSelectors()

◆ MakeSamplingTables()

void G4RiGeMuPairProductionModel::MakeSamplingTables ( )
protected

Definition at line 457 of file G4RiGeMuPairProductionModel.cc.

458{
460
461 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
462
463 G4double Z = ZDATPAIR[iz];
464 G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
465 G4double kinEnergy = emin;
466
467 for (std::size_t it=0; it<=nbine; ++it) {
468
469 pv->PutY(it, G4Log(kinEnergy/CLHEP::MeV));
470 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy, Z);
471 /*
472 G4cout << "it= " << it << " E= " << kinEnergy
473 << " " << particle->GetParticleName()
474 << " maxE= " << maxPairEnergy << " minE= " << minPairEnergy
475 << " ymin= " << ymin << G4endl;
476 */
477 G4double coef = G4Log(minPairEnergy/kinEnergy)/ymin;
478 G4double ymax = G4Log(maxPairEnergy/kinEnergy)/coef;
479 G4double fac = (ymax - ymin)/dy;
480 std::size_t imax = (std::size_t)fac;
481 fac -= (G4double)imax;
482
483 G4double xSec = 0.0;
484 G4double x = ymin;
485 /*
486 G4cout << "Z= " << currentZ << " z13= " << z13
487 << " mE= " << maxPairEnergy << " ymin= " << ymin
488 << " dy= " << dy << " c= " << coef << G4endl;
489 */
490 // start from zero
491 pv->PutValue(0, it, 0.0);
492 if(0 == it) { pv->PutX(nbiny, 0.0); }
493
494 for (std::size_t i=0; i<nbiny; ++i) {
495
496 if(0 == it) { pv->PutX(i, x); }
497
498 if(i < imax) {
499 G4double ep = kinEnergy*G4Exp(coef*(x + dy*0.5));
500
501 // not multiplied by interval, because table
502 // will be used only for sampling
503 //G4cout << "i= " << i << " x= " << x << "E= " << kinEnergy
504 // << " Egamma= " << ep << G4endl;
505 xSec += ep*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
506
507 // last bin before the kinematic limit
508 } else if(i == imax) {
509 G4double ep = kinEnergy*G4Exp(coef*(x + fac*dy*0.5));
510 xSec += ep*fac*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
511 }
512 pv->PutValue(i + 1, it, xSec);
513 x += dy;
514 }
515 kinEnergy *= factore;
516
517 // to avoid precision lost
518 if(it+1 == nbine) { kinEnergy = emax; }
519 }
520 fElementData->InitialiseForElement(iz, pv);
521 }
522}
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 G4RiGeMuPairProductionModel::MaxSecondaryEnergyForElement ( G4double kineticEnergy,
G4double Z )
inlineprotected

Definition at line 191 of file G4RiGeMuPairProductionModel.hh.

193{
194 G4int Z = G4lrint(ZZ);
195 if(Z != currentZ) {
196 currentZ = Z;
197 z13 = nist->GetZ13(Z);
198 z23 = z13*z13;
199 lnZ = nist->GetLOGZ(Z);
200 }
201 return kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
202}

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

◆ MinPrimaryEnergy()

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

Reimplemented from G4VEmModel.

Definition at line 122 of file G4RiGeMuPairProductionModel.cc.

125{
126 return std::max(lowestKinEnergy, cut);
127}

◆ operator=()

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

◆ RetrieveTables()

G4bool G4RiGeMuPairProductionModel::RetrieveTables ( )
protected

Definition at line 981 of file G4RiGeMuPairProductionModel.cc.

982{
983 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
984 G4double Z = ZDATPAIR[iz];
985 G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
986 std::ostringstream ss;
987 ss << G4EmParameters::Instance()->GetDirLEDATA() << "/mupair/"
988 << particle->GetParticleName() << Z << ".dat";
989 std::ifstream infile(ss.str(), std::ios::in);
990 if(!pv->Retrieve(infile)) {
991 delete pv;
992 return false;
993 }
994 fElementData->InitialiseForElement(iz, pv);
995 }
996 return true;
997}
const G4String & GetDirLEDATA() const
G4bool Retrieve(std::ifstream &fIn)

Referenced by Initialise().

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 526 of file G4RiGeMuPairProductionModel.cc.

531{
532 G4double eMass = CLHEP::electron_mass_c2;
533 G4double eMass2 = eMass*eMass;
534
535 // Energy and momentum of the pramary particle
536 G4double kinEnergy = aDynamicParticle->GetKineticEnergy();
537 G4double totEnergy = aDynamicParticle->GetTotalEnergy();
538 G4double particleMomentum = aDynamicParticle->GetTotalMomentum();
539 G4ThreeVector particleMomentumVector = aDynamicParticle->GetMomentum();
540 G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
541
542 G4double minQ2 = 4.*eMass2;
543 G4double maxQ2 = (kinEnergy - particleMass)*(kinEnergy - particleMass);
544 G4double intervalQ2 = G4Log(maxQ2/minQ2);
545
546 // Square invariant of mass of the pair
547 G4double Q2 = minQ2*G4Exp(intervalQ2*randNumbs[4]);
548
549 G4double mingEnergy = std::sqrt(Q2);
550 G4double maxgEnergy = kinEnergy - particleMass;
551 G4double intervalgEnergy = maxgEnergy - mingEnergy;
552
553 // Energy of virtual gamma
554 G4double gEnergy = mingEnergy + intervalgEnergy*randNumbs[5];
555
556 // Momentum module of the virtual gamma
557 G4double gMomentum = std::sqrt(gEnergy*gEnergy - Q2);
558
559 // Energy and momentum module of the outgoing parent particle
560 G4double particleFinalEnergy = kinEnergy - gEnergy;
561 G4double particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
563
564 G4double mint3 = 0.;
565 G4double maxt3 = CLHEP::pi;
566 G4double Cmin = std::cos(maxt3);
567 G4double Cmax = std::cos(mint3);
568
569 //G4cout << "------- G4RiGeMuPairProductionModel::SampleSecondaries E(MeV)= "
570 // << kinEnergy << " "
571 // << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
572
573 // select randomly one element constituing the material
574 const G4Element* anElement = SelectRandomAtom(couple,particle,kinEnergy);
575
576 // define interval of energy transfer
577 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy,
578 anElement->GetZ());
579 G4double maxEnergy = std::min(tmax, maxPairEnergy);
580 G4double minEnergy = std::max(tmin, minPairEnergy);
581
582 if (minEnergy >= maxEnergy) { return; }
583
584 //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
585 // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
586 // << " ymin= " << ymin << " dy= " << dy << G4endl;
587
588 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
589
591
592 // compute limits
593 G4double yymin = G4Log(minEnergy/kinEnergy)/coeff;
594 G4double yymax = G4Log(maxEnergy/kinEnergy)/coeff;
595
596 //G4cout << "yymin= " << yymin << " yymax= " << yymax << G4endl;
597
598 // units should not be used, bacause table was built without
599 G4double logTkin = G4Log(kinEnergy/CLHEP::MeV);
600
601 // sample e-e+ energy, pair energy first
602
603 // select sample table via Z
604 G4int iz1(0), iz2(0);
605 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
606 if(currentZ == ZDATPAIR[iz]) {
607 iz1 = iz2 = iz;
608 break;
609 } else if(currentZ < ZDATPAIR[iz]) {
610 iz2 = iz;
611 if(iz > 0) { iz1 = iz-1; }
612 else { iz1 = iz2; }
613 break;
614 }
615 }
616 if (0 == iz1) { iz1 = iz2 = NZDATPAIR-1; }
617
618 G4double pairEnergy = 0.0;
619 G4int count = 0;
620 //G4cout << "start loop Z1= " << iz1 << " Z2= " << iz2 << G4endl;
621 do {
622 ++count;
623 // sampling using only one random number
624 G4double rand = rndmEngine->flat();
625
626 G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
627 if(iz1 != iz2) {
628 G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
629 G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1]);
630 G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2]);
631 //G4cout << count << ". x= " << x << " x2= " << x2
632 // << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
633 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
634 }
635 //G4cout << "x= " << x << " coeff= " << coeff << G4endl;
636 pairEnergy = kinEnergy*G4Exp(x*coeff);
637
638 // Loop checking, 30-Oct-2024, Vladimir Ivanchenko
639 } while((pairEnergy < minEnergy || pairEnergy > maxEnergy) && 50 > count);
640
641 //G4cout << "## pairEnergy(GeV)= " << pairEnergy/GeV
642 // << " Etot(GeV)= " << totalEnergy/GeV << G4endl;
643 rndmEngine->flatArray(9, randNumbs);
644 G4double phi3 = CLHEP::twopi*randNumbs[0];
645 fAngularGenerator->PhiRotation(partDirection, phi3);
646
647 G4LorentzVector muF;
648 G4ThreeVector eDirection, pDirection;
649 G4double eEnergy, pEnergy;
650
651 if (randNumbs[7] < W[0]) {
652 G4double A1 = -(Q2 - 2.*kinEnergy*gEnergy);
653 G4double B1 = -(2.*gMomentum*particleMomentum);
654 G4double tginterval = G4Log((A1 + B1)/(A1 - B1))/B1;
655
656 G4double costg = (-A1 + (A1 - B1)*G4Exp(B1*tginterval*randNumbs[1]))/B1;
657 G4double sintg = std::sqrt((1.0 - costg)*(1.0 + costg));
658 G4double phig = CLHEP::twopi*randNumbs[2];
659 G4double sinpg = std::sin(phig);
660 G4double cospg = std::cos(phig);
661
662 G4ThreeVector dirGamma;
663 dirGamma.set(sintg*cospg, sintg*sinpg, costg);
664 G4LorentzVector gFourMomentum(dirGamma*gMomentum, gEnergy);
665
666 G4double Ap = particleMomentum*particleMomentum +
667 particleFinalMomentum*particleFinalMomentum + gMomentum*gMomentum;
668 G4double A = Ap - 2.*particleMomentum*gMomentum*costg;
669 G4double B = 2.*particleMomentum*gMomentum*sintg*cospg;
670 G4double C = 2.*particleFinalMomentum*gMomentum*costg -
671 2.*particleMomentum*particleFinalMomentum;
672 G4double absB = std::abs(B);
673 G4double t1interval = (1./(A + C + absB*mint3) - 1./(A + C + absB*maxt3))/absB;
674 G4double t1 = (-(A + C) + 1./(1./(A + C + absB*mint3) - absB*t1interval*randNumbs[0]))/absB;
675 G4double sint1 = std::sin(t1);
676 G4double cost1 = std::cos(t1);
677
678 // Ingoing parent particle change
679 G4double Phi = CLHEP::twopi*randNumbs[3];
680 partDirection.set(sint1, 0., cost1);
681 fAngularGenerator->PhiRotation(partDirection, Phi);
682 kinEnergy = particleFinalEnergy;
683
684 G4double cost5 = -1. + 2.*randNumbs[6];
685 G4double phi5 = CLHEP::twopi*randNumbs[8];
686
687 G4LorentzVector eFourMomentumMQ = fAngularGenerator->eDP2(Q2, eMass2, eMass2, cost5, phi5);
688 G4LorentzVector pFourMomentumMQ = fAngularGenerator->pDP2(eMass2, eFourMomentumMQ);
689
690 G4LorentzVector eFourMomentum = eFourMomentumMQ.boost(gFourMomentum.boostVector());
691 G4LorentzVector pFourMomentum = pFourMomentumMQ.boost(gFourMomentum.boostVector());
692
693 eEnergy = eFourMomentum.t();
694 pEnergy = pFourMomentum.t();
695
696 } else if (randNumbs[7] >= W[0] && randNumbs[7] < W[1]) {
697 G4double A3 = Q2 + 2.*gEnergy*particleFinalEnergy;
698 G4double B3 = -2.*gMomentum*particleFinalMomentum;
699
700 G4double tQ3interval = G4Log((A3 + B3)/(A3 - B3))/B3;
701 G4double tQMG = (-A3 + (A3 - B3)*G4Exp(B3*tQ3interval*randNumbs[0]))/B3;
702 G4double phiQP = CLHEP::twopi*randNumbs[2];
703
704 G4double sintQ3 = std::sqrt(1. - tQMG*tQMG);
705 G4double cospQP = std::cos(phiQP);
706 G4double sinpQP = std::sin(phiQP);
707
708 G4double Ap = particleMomentum*particleMomentum +
709 particleFinalMomentum*particleFinalMomentum + gMomentum*gMomentum;
710 G4double A = Ap + 2.*particleFinalMomentum*gMomentum*tQMG;
711 G4double B = -2.*particleMomentum*gMomentum*sintQ3*cospQP;
712 G4double C = -2.*particleMomentum*gMomentum*tQMG - 2.*particleMomentum*particleFinalMomentum;
713
714 G4double absB = std::abs(B);
715 G4double t3interval = (1./(A + C + absB*mint3) - 1./(A + C + absB*maxt3))/absB;
716 G4double t3 = (-(A + C) + 1./(1./(A + C + absB*mint3) - absB*t3interval*randNumbs[0]))/absB;
717 G4double sint3 = std::sin(t3);
718 G4double cost3 = std::cos(t3);
719
720 G4double cost = -sint3*sintQ3*cospQP + cost3*tQMG;
721 G4double sint = std::sqrt((1. + cost)*(1. - cost));
722 G4double cosp = (sintQ3*cospQP*cost3 + sint3*tQMG)/sint;
723 G4double sinp = sintQ3*sinpQP/sint;
724
725 G4ThreeVector dirGamma;
726 dirGamma.set(sint*cosp, sint*sinp, cost);
727 G4LorentzVector gFourMomentum(dirGamma*gMomentum, gEnergy);
728
729 // Ingoing parent particle change
730 G4double Phi = CLHEP::twopi*randNumbs[3];
731 partDirection.set(sint3, 0., cost3);
732 fAngularGenerator->PhiRotation(partDirection, Phi);
733 kinEnergy = particleFinalEnergy;
734
735 G4double cost5 = -1. + 2.*randNumbs[6];
736 G4double phi5 = CLHEP::twopi*randNumbs[8];
737
738 G4LorentzVector eFourMomentumMQ = fAngularGenerator->eDP2(Q2, eMass2, eMass2, cost5, phi5);
739 G4LorentzVector pFourMomentumMQ = fAngularGenerator->pDP2(eMass2, eFourMomentumMQ);
740
741 G4LorentzVector eFourMomentum = eFourMomentumMQ.boost(gFourMomentum.boostVector());
742 G4LorentzVector pFourMomentum = pFourMomentumMQ.boost(gFourMomentum.boostVector());
743
744 eEnergy = eFourMomentum.t();
745 pEnergy = pFourMomentum.t();
746
747 } else if (randNumbs[7] >= W[1] && randNumbs[7] < W[2]) {
748 G4double phi5 = CLHEP::twopi*randNumbs[1];
749 G4double phi6 = CLHEP::twopi*randNumbs[2];
750 G4double muEnergyInterval = kinEnergy - 2.*eMass - particleMass;
751 particleFinalEnergy = particleMass + muEnergyInterval*randNumbs[3];
752 particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
754
755 G4double mineEnergy = eMass;
756 G4double maxeEnergy = kinEnergy - particleFinalEnergy - eMass;
757 G4double eEnergyinterval = maxeEnergy - mineEnergy;
758 eEnergy = mineEnergy + eEnergyinterval*randNumbs[4];
759
760 G4double cosp3 = 1.;
761 G4double sinp3 = 0.;
762 G4double cosp5 = std::cos(phi5);
763 G4double sinp5 = std::sin(phi5);
764 G4double cosp6 = std::cos(phi6);
765 G4double sinp6 = std::sin(phi6);
766
767 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
768 pEnergy = kinEnergy - particleFinalEnergy - eEnergy;
769 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
770
771 G4double A3 = -2.*particleMass*particleMass + 2.*kinEnergy*particleFinalEnergy;
772 G4double B3 = -2.*particleMomentum*particleFinalMomentum;
773 G4double cost3interval = G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
774 G4double expanCost3r6 = G4Exp(B3*cost3interval*randNumbs[5]);
775 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
776 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
777
778 partDirection.set(sint3, 0., cost3);
779
780 G4ThreeVector muFinalMomentumVector;
781 muFinalMomentumVector.set(particleFinalMomentum*sint3, 0., particleFinalMomentum*cost3);
782
783 G4LorentzVector muFourMomentum(particleMomentumVector, totEnergy);
784 G4LorentzVector muFinalFourMomentum(muFinalMomentumVector, particleFinalEnergy);
785 G4LorentzVector auxVec1 = muFourMomentum - muFinalFourMomentum;
786 G4double A5 = auxVec1.mag2() - 2.*eEnergy*(kinEnergy - particleFinalEnergy) +
787 2.*particleMomentumVector[2]*eMomentum - 2.*particleFinalMomentum*eMomentum*cost3;
788 G4double B5 = -2.*particleFinalMomentum*eMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5);
789 G4double absA5 = std::abs(A5);
790 G4double absB5 = std::abs(B5);
791 G4double mint5 = 0.;
792 G4double maxt5 = CLHEP::pi;
793 G4double t5interval = G4Log((absA5 + absB5*maxt5)/(absA5 + absB5*mint5))/absB5;
794 G4double argexp = absB5*t5interval*randNumbs[6] + G4Log(absA5 + absB5*mint5);
795 G4double t5 = -absA5/absB5 + G4Exp(argexp)/absB5;
796 G4double sint5 = std::sin(t5);
797 G4double cost5 = std::cos(t5);
798
799 eDirection.set(sint5*cosp5, sint5*sinp5, cost5);
800 G4ThreeVector eMomentumVector = eMomentum*eDirection;
801
802 G4ThreeVector auxVec2 = particleMomentumVector - muFinalMomentumVector - eMomentumVector;
803 G4double p1mp3mp52 = auxVec2.dot(auxVec2);
804 G4double Bp = particleFinalMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6) +
805 eMomentum*(sint5*cosp5*cosp6 + sint5*sinp5*sinp6);
806 G4double Cp = -particleMomentum + particleFinalMomentum*cost3 + eMomentum*cost5;
807 G4double A6 = p1mp3mp52 + pMomentum*pMomentum;
808 G4double B6 = 2.*pMomentum*Bp;
809 G4double C6 = 2.*pMomentum*Cp;
810 G4double mint6 = 0.;
811 G4double maxt6 = CLHEP::pi;
812 G4double absA6C6 = std::abs(A6 + C6);
813 G4double absB6 = std::abs(B6);
814 G4double t6interval = (1./(absA6C6 + absB6*mint6) - 1./(absA6C6 + absB6*maxt6))/absB6;
815 G4double t6 = (-absA6C6 + 1./(1./(absA6C6 + absB6*mint6) - absB6*t6interval*randNumbs[8]))/absB6;
816 G4double sint6 = std::sin(t6);
817 G4double cost6 = std::cos(t6);
818
819 pDirection.set(sint6*cosp6, sint6*sinp6, cost6);
820
821 } else {
822 G4double phi6 = CLHEP::twopi*randNumbs[1];
823 G4double phi5 = CLHEP::twopi*randNumbs[2];
824 G4double muFinalEnergyinterval = kinEnergy - 2.*eMass - particleMass;
825 particleFinalEnergy = particleMass + muFinalEnergyinterval*randNumbs[3];
826 particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
828
829 G4double maxpEnergy = kinEnergy - particleFinalEnergy - eMass;
830 G4double pEnergyinterval = maxpEnergy - eMass;
831 pEnergy = eMass + pEnergyinterval*randNumbs[4];
832
833 G4double cosp3 = 1.;
834 G4double sinp3 = 0.;
835 G4double cosp5 = std::cos(phi5);
836 G4double sinp5 = std::sin(phi5);
837 G4double cosp6 = std::cos(phi6);
838 G4double sinp6 = std::sin(phi6);
839
840 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
841 eEnergy = kinEnergy - particleFinalEnergy - pEnergy;
842 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
843
844 G4double A3 = -2.*particleMass*particleMass + 2.*kinEnergy*particleFinalEnergy;
845 G4double B3 = -2.*particleMomentum*particleFinalMomentum;
846 G4double cost3interval = G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
847 G4double expanCost3r6 = G4Exp(B3*cost3interval*randNumbs[5]);
848 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
849 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
850
851 partDirection.set(sint3*cosp3, sint3*sinp3, cost3);
852
853 G4ThreeVector muFinalMomentumVector;
854 muFinalMomentumVector.set(particleFinalMomentum*sint3*cosp3,
855 particleFinalMomentum*sint3*sinp3,
856 particleFinalMomentum*cost3);
857
858 G4LorentzVector muFourMomentum(particleMomentum, particleMomentumVector);
859 G4LorentzVector muFinalFourMomentum(particleFinalEnergy, muFinalMomentumVector);
860 G4LorentzVector auxVec1 = muFourMomentum - muFinalFourMomentum;
861 G4double A6 = auxVec1.mag2() -
862 2.*pEnergy*(kinEnergy - particleFinalEnergy) + 2.*particleMomentumVector[2]*pMomentum -
863 2.*particleFinalMomentum*pMomentum*cost3;
864 G4double B6 = -2.*particleFinalMomentum*pMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6);
865 G4double absA6 = std::abs(A6);
866 G4double absB6 = std::abs(B6);
867 G4double mint6 = 0.;
868 G4double maxt6 = CLHEP::pi;
869 G4double t6interval = G4Log((absA6 + absB6*maxt6)/(absA6 + absB6*mint6))/absB6;
870 G4double argexp = absB6*t6interval*randNumbs[6] + G4Log(absA6 + absB6*mint6);
871 G4double t6 = -absA6/absB6 + G4Exp(argexp)/absB6;
872 G4double sint6 = std::sin(t6);
873 G4double cost6 = std::cos(t6);
874
875 pDirection.set(sint6*cosp6, sint6*sinp6, cost6);
876 G4ThreeVector pMomentumVector = pMomentum*pDirection;
877
878 G4ThreeVector auxVec2 = particleMomentumVector - muFinalMomentumVector - pMomentumVector;
879 G4double p1mp3mp62 = auxVec2.dot(auxVec2);
880 G4double Bp = particleFinalMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5) +
881 pMomentum*(sint6*cosp6*cosp5 + sint6*sinp6*sinp5);
882 G4double Cp = -particleMomentum + particleFinalMomentum*cost3 + pMomentum*cost6;
883 G4double A5 = p1mp3mp62 + eMomentum*eMomentum;
884 G4double B5 = 2.*eMomentum*Bp;
885 G4double C5 = 2.*eMomentum*Cp;
886 G4double mint5 = 0.;
887 G4double maxt5 = CLHEP::pi;
888 G4double absA5C5 = std::abs(A5 + C5);
889 G4double absB5 = std::abs(B5);
890 G4double t5interval = (1./(absA5C5 + absB5*mint5) - 1./(absA5C5 + absB5*maxt5))/absB5;
891 G4double t5 = (-absA5C5 + 1./(1./(absA5C5 + absB5*mint5) - absB5*t5interval*randNumbs[8]))/absB5;
892 G4double sint5 = std::sin(t5);
893 G4double cost5 = std::cos(t5);
894
895 eDirection.set(sint5*cosp5, sint5*sinp5, cost5);
896 }
897 // these lines are closed temporary, will be enabled or removed after testing
898 /*
899 fAngularGenerator->Sample5DPairDirections(aDynamicParticle, eDirection, pDirection,
900 gEnergy, Q2, gMomentum,
901 particleFinalMomentum,
902 particleFinalEnergy,
903 randNumbs, W);
904 */
905 // create G4DynamicParticle object for e+e-
906 auto aParticle1 = new G4DynamicParticle(theElectron, eDirection, eEnergy);
907 auto aParticle2 = new G4DynamicParticle(thePositron, pDirection, pEnergy);
908
909 // Fill output vector
910 vdp->push_back(aParticle1);
911 vdp->push_back(aParticle2);
912
913 // if energy transfer is higher than threshold (very high by default)
914 // then stop tracking the primary particle and create a new secondary
915 if (pairEnergy > SecondaryThreshold()) {
916 fParticleChange->ProposeTrackStatus(fStopAndKill);
917 fParticleChange->SetProposedKineticEnergy(0.0);
918 auto newdp = new G4DynamicParticle(particle, muF);
919 vdp->push_back(newdp);
920 } else { // continue tracking the primary e-/e+ otherwise
921 fParticleChange->SetProposedMomentumDirection(muF.vect().unit());
922 G4double ekin = std::max(muF.e() - particleMass, 0.0);
923 fParticleChange->SetProposedKineticEnergy(ekin);
924 }
925 //G4cout << "-- G4RiGeMuPairProductionModel::SampleSecondaries done" << G4endl;
926}
G4double C(G4double temp)
G4double B(G4double temperature)
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
@ fStopAndKill
const G4double A[17]
#define C5
Hep3Vector unit() const
double dot(const Hep3Vector &) const
void set(double x, double y, double z)
HepLorentzVector & boost(double, double, double)
Hep3Vector vect() const
virtual double flat()=0
virtual void flatArray(const int size, double *vect)=0
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
G4ThreeVector GetMomentum() const
G4double GetTotalMomentum() const
G4double GetZ() const
Definition G4Element.hh:119
G4double FindScaledEnergy(G4int Z, G4double rand, G4double logTkin, G4double yymin, G4double yymax)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double SecondaryThreshold() const
#define W
Definition crc32.c:85

◆ SetLowestKineticEnergy()

void G4RiGeMuPairProductionModel::SetLowestKineticEnergy ( G4double e)
inline

Definition at line 172 of file G4RiGeMuPairProductionModel.hh.

173{
174 lowestKinEnergy = e;
175}

◆ SetParticle()

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

Definition at line 180 of file G4RiGeMuPairProductionModel.hh.

181{
182 if(nullptr == particle) {
183 particle = p;
184 particleMass = particle->GetPDGMass();
185 }
186}

Referenced by G4RiGeMuPairProductionModel(), and Initialise().

◆ StoreTables()

void G4RiGeMuPairProductionModel::StoreTables ( ) const
protected

Definition at line 963 of file G4RiGeMuPairProductionModel.cc.

964{
965 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
966 G4int Z = ZDATPAIR[iz];
967 G4Physics2DVector* pv = fElementData->GetElement2DData(Z);
968 if(nullptr == pv) {
969 DataCorrupted(Z, 1.0);
970 return;
971 }
972 std::ostringstream ss;
973 ss << "mupair/" << particle->GetParticleName() << Z << ".dat";
974 std::ofstream outfile(ss.str());
975 pv->Store(outfile);
976 }
977}
void Store(std::ofstream &fOut) const

Referenced by Initialise().

Member Data Documentation

◆ currentZ

G4int G4RiGeMuPairProductionModel::currentZ = 0
protected

◆ dy

G4double G4RiGeMuPairProductionModel::dy = 0.005
protected

Definition at line 143 of file G4RiGeMuPairProductionModel.hh.

Referenced by Initialise(), and MakeSamplingTables().

◆ emax

G4double G4RiGeMuPairProductionModel::emax
protected

◆ emin

G4double G4RiGeMuPairProductionModel::emin
protected

◆ factorForCross

G4double G4RiGeMuPairProductionModel::factorForCross
protected

◆ fParticleChange

G4ParticleChangeForLoss* G4RiGeMuPairProductionModel::fParticleChange = nullptr
protected

Definition at line 126 of file G4RiGeMuPairProductionModel.hh.

Referenced by Initialise(), and SampleSecondaries().

◆ fTableToFile

G4bool G4RiGeMuPairProductionModel::fTableToFile = false
protected

Definition at line 153 of file G4RiGeMuPairProductionModel.hh.

Referenced by Initialise().

◆ lnZ

G4double G4RiGeMuPairProductionModel::lnZ = 0.0
protected

◆ lowestKinEnergy

◆ minPairEnergy

◆ nbine

std::size_t G4RiGeMuPairProductionModel::nbine = 0
protected

Definition at line 151 of file G4RiGeMuPairProductionModel.hh.

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

◆ nbiny

std::size_t G4RiGeMuPairProductionModel::nbiny = 1000
protected

Definition at line 150 of file G4RiGeMuPairProductionModel.hh.

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

◆ NINTPAIR

const G4int G4RiGeMuPairProductionModel::NINTPAIR = 8
staticprotected

◆ nist

G4NistManager* G4RiGeMuPairProductionModel::nist = nullptr
protected

◆ nYBinPerDecade

G4int G4RiGeMuPairProductionModel::nYBinPerDecade = 4
protected

Definition at line 149 of file G4RiGeMuPairProductionModel.hh.

Referenced by Initialise().

◆ NZDATPAIR

const G4int G4RiGeMuPairProductionModel::NZDATPAIR = 5
staticprotected

◆ particle

const G4ParticleDefinition* G4RiGeMuPairProductionModel::particle = nullptr
protected

◆ particleMass

G4double G4RiGeMuPairProductionModel::particleMass = 0.0
protected

◆ randNumbs

G4double G4RiGeMuPairProductionModel::randNumbs[9]
protected

◆ sqrte

G4double G4RiGeMuPairProductionModel::sqrte
protected

◆ wgi

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

Definition at line 76 of file G4RiGeMuPairProductionModel.hh.

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

◆ xgi

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

Definition at line 71 of file G4RiGeMuPairProductionModel.hh.

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

◆ ymin

G4double G4RiGeMuPairProductionModel::ymin = -5.0
protected

◆ z13

G4double G4RiGeMuPairProductionModel::z13 = 0.0
protected

◆ z23

G4double G4RiGeMuPairProductionModel::z23 = 0.0
protected

◆ ZDATPAIR

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

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