| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4RiGeMuPairProductionModel
33 //
34 // Authors: Girardo Depaola & Ricardo Pa
35 //
36 // Creation date: 29.10.2024
37 //
38 //
39 // -------------------------------------------
40 //
41 //....oooOO0OOooo........oooOO0OOooo........oo
42 //....oooOO0OOooo........oooOO0OOooo........oo
43
44 #include "G4RiGeMuPairProductionModel.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4EmParameters.hh"
48 #include "G4Electron.hh"
49 #include "G4Positron.hh"
50 #include "G4MuonMinus.hh"
51 #include "G4MuonPlus.hh"
52 #include "Randomize.hh"
53 #include "G4Material.hh"
54 #include "G4Element.hh"
55 #include "G4ElementVector.hh"
56 #include "G4ElementDataRegistry.hh"
57 #include "G4ProductionCutsTable.hh"
58 #include "G4ParticleChangeForLoss.hh"
59 #include "G4RiGeAngularGenerator.hh"
60 #include "G4Log.hh"
61 #include "G4Exp.hh"
62 #include "G4AutoLock.hh"
63
64 #include <iostream>
65 #include <fstream>
66
67 //....oooOO0OOooo........oooOO0OOooo........oo
68
69 const G4int G4RiGeMuPairProductionModel::ZDATP
70
71 const G4double G4RiGeMuPairProductionModel::xg
72 0.0198550717512320, 0.1016667612931865, 0.
73 0.5917173212478250, 0.7627662049581645, 0.
74 };
75
76 const G4double G4RiGeMuPairProductionModel::wg
77 0.0506142681451880, 0.1111905172266870, 0.
78 0.1813418916891810, 0.1568533229389435, 0.
79 };
80
81 namespace
82 {
83 G4Mutex theRiGeMuPairMutex = G4MUTEX_INITIAL
84
85 const G4double ak1 = 6.9;
86 const G4double ak2 = 1.0;
87
88 // Channel weights
89 const G4double W[3] = {0.25, 0.5, 0.75};
90 }
91
92 //....oooOO0OOooo........oooOO0OOooo........oo
93
94 G4RiGeMuPairProductionModel::G4RiGeMuPairProdu
95 : G4VEmModel("muPairProdRiGe"),
96 factorForCross(CLHEP::fine_structure_const
97 CLHEP::classic_electr_radius*CLHEP::cla
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 {
103 nist = G4NistManager::Instance();
104
105 theElectron = G4Electron::Electron();
106 thePositron = G4Positron::Positron();
107
108 if (nullptr != p) {
109 SetParticle(p);
110 lowestKinEnergy = std::max(lowestKinEnergy
111 }
112 emin = lowestKinEnergy;
113 emax = emin*10000.;
114 fAngularGenerator = new G4RiGeAngularGenerat
115 SetAngularDistribution(fAngularGenerator);
116 for (G4int i=0; i<9; ++i) { randNumbs[i] = 0
117 }
118
119 //....oooOO0OOooo........oooOO0OOooo........oo
120
121 G4double
122 G4RiGeMuPairProductionModel::MinPrimaryEnergy(
123
124
125 {
126 return std::max(lowestKinEnergy, cut);
127 }
128
129 //....oooOO0OOooo........oooOO0OOooo........oo
130
131 void G4RiGeMuPairProductionModel::Initialise(c
132 c
133 {
134 SetParticle(p);
135
136 if (nullptr == fParticleChange) {
137 fParticleChange = GetParticleChangeForLoss
138
139 // define scale of internal table for each
140 if (0 == nbine) {
141 emin = std::max(lowestKinEnergy, LowEner
142 emax = std::max(HighEnergyLimit(), emin*
143 nbine = std::size_t(nYBinPerDecade*std::
144 if(nbine < 3) { nbine = 3; }
145
146 ymin = G4Log(minPairEnergy/emin);
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 s
170 if(lowestKinEnergy >= HighEnergyLimit()) { r
171
172 if (p == particle) {
173 auto data = G4ElementDataRegistry::Instanc
174 fElementData = data->GetElementDataByName(
175 if (nullptr == fElementData) {
176 G4AutoLock l(&theRiGeMuPairMutex);
177 fElementData = data->GetElementDataByNam
178 if (nullptr == fElementData) {
179 fElementData = new G4ElementData(NZDAT
180 fElementData->SetName(dataName);
181 }
182 G4bool useDataFile = G4EmParameters::Ins
183 if (useDataFile) { useDataFile = Retrie
184 if (!useDataFile) { MakeSamplingTables()
185 if (fTableToFile) { StoreTables(); }
186 l.unlock();
187 }
188 if (IsMaster()) {
189 InitialiseElementSelectors(p, cuts);
190 }
191 }
192 }
193
194 //....oooOO0OOooo........oooOO0OOooo........oo
195
196 void G4RiGeMuPairProductionModel::InitialiseLo
197
198 {
199 if(p == particle && lowestKinEnergy < HighEn
200 SetElementSelectors(masterModel->GetElemen
201 }
202 }
203
204 //....oooOO0OOooo........oooOO0OOooo........oo
205
206 G4double
207 G4RiGeMuPairProductionModel::ComputeDEDXPerVol
208
209
210
211 {
212 G4double dedx = 0.0;
213 if (cutEnergy <= minPairEnergy || kineticEne
214 { return dedx; }
215
216 const G4ElementVector* theElementVector = ma
217 const G4double* theAtomicNumDensityVector =
218 material->G
219
220 // loop for elements in the material
221 for (std::size_t i=0; i<material->GetNumberO
222 G4double Z = (*theElementVector)[i]->GetZ
223 G4double tmax = MaxSecondaryEnergyForElem
224 G4double loss = ComputMuPairLoss(Z, kinet
225 dedx += loss*theAtomicNumDensityVector[i]
226 }
227 dedx = std::max(dedx, 0.0);
228 return dedx;
229 }
230
231 //....oooOO0OOooo........oooOO0OOooo........oo
232
233 G4double G4RiGeMuPairProductionModel::ComputMu
234
235
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)
244 G4double aaa = G4Log(minPairEnergy);
245 G4double bbb = G4Log(cut);
246
247 G4int kkk = std::min(std::max(G4lrint((bbb-a
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*ComputeDMicroscopi
255 }
256 x += hhh;
257 }
258 loss *= hhh;
259 loss = std::max(loss, 0.0);
260 return loss;
261 }
262
263 //....oooOO0OOooo........oooOO0OOooo........oo
264
265 G4double
266 G4RiGeMuPairProductionModel::ComputeMicroscopi
267 G4double Z,
268 G4double cutEnergy)
269 {
270 G4double cross = 0.;
271 G4double tmax = MaxSecondaryEnergyForElement
272 G4double cut = std::max(cutEnergy, minPairE
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-a
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]*ComputeDMicroscopicCr
286 }
287 x += hhh;
288 }
289
290 cross *= hhh;
291 cross = std::max(cross, 0.0);
292 return cross;
293 }
294
295 //....oooOO0OOooo........oooOO0OOooo........oo
296
297 G4double G4RiGeMuPairProductionModel::ComputeD
298 G4d
299 G4d
300 G4d
301 // Calculates the differential (D) microscopi
302 // using the cross section formula of R.P. Kok
303 // Code modified by R.P. Kokoulin, V.N. Ivanch
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 - pairEn
317
318 if (residEnergy <= 0.75*sqrte*z13*particleMa
319 return 0.0;
320
321 G4double a0 = 1.0 / (totalEnergy * residEner
322 G4double alf = 4.0 * electron_mass_c2 / pair
323 G4double rt = std::sqrt(1.0 - alf);
324 G4double delta = 6.0 * particleMass * partic
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::ele
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 =
338 else { bbb = bbbtf; g1 = g1tf; g2 =
339
340 G4double zeta = 0.0;
341 G4double z1exp = totalEnergy / (particleMass
342
343 // 35.221047195922 is the root of zeta1(x) =
344 // condition below is the same as zeta1 > 0.
345 if (z1exp > 35.221047195922)
346 {
347 G4double z2exp = totalEnergy / (particleMa
348 zeta = (0.073 * G4Log(z1exp) - 0.26) / (0.
349 }
350
351 G4double z2 = Z*(Z+zeta);
352 G4double screen0 = 2.*electron_mass_c2*sqrte
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
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 =
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) *
381 G4double yed = b62*G4Log(3.0 + xii[i]) + (
382
383 G4double ymu = b62 * (1.0 + rho2[i]) + 6.0
384 G4double ymd = (b40 + 3.0)*(1.0 + rho2[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 + xi
398 (1.0 - rho2[i] - beta)/xi1[i] - (3.0 + rho2[
399 } else {
400 be[i] = 0.5*(3.0 - rho2[i] + 2.0*beta*(1
401 }
402
403 if(xi[i] >= 0.001) {
404 G4double a10 = (1.0 + 2.0 * beta) * (1.0
405 bm[i] = ((1.0 + rho2[i])*(1.0 + 1.5 * be
406 xi[i] * (1.0 - rho2[i] - beta)
407 } else {
408 bm[i] = 0.5*(5.0 - rho2[i] + beta * (3.0
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 - rh
416 G4double ale = G4Log(bbb/z13*std::sqrt(xi1
417 G4double cre = 0.5*G4Log(1. + 2.25*z23*xi1
418
419 G4double fe = (ale-cre)*be[i];
420 fe = std::max(fe, 0.0);
421
422 G4double alm_crm = G4Log(bbb*massratio/(1.
423 G4double fm = std::max(alm_crm*bm[i], 0.0)
424
425 sum += wgi[i]*(1.0 + rho[i])*(fe + fm);
426 }
427
428 return -tmn*sum*factorForCross*z2*residEnerg
429 }
430
431 //....oooOO0OOooo........oooOO0OOooo........oo
432
433 G4double
434 G4RiGeMuPairProductionModel::ComputeCrossSecti
435
436
437
438
439 {
440 G4double cross = 0.0;
441 if (kineticEnergy <= lowestKinEnergy) { retu
442
443 G4double maxPairEnergy = MaxSecondaryEnergyF
444 G4double tmax = std::min(maxEnergy, maxPairE
445 G4double cut = std::max(cutEnergy, minPairE
446 if (cut >= tmax) { return cross; }
447
448 cross = ComputeMicroscopicCrossSection(kinet
449 if(tmax < kineticEnergy) {
450 cross -= ComputeMicroscopicCrossSection(ki
451 }
452 return cross;
453 }
454
455 //....oooOO0OOooo........oooOO0OOooo........oo
456
457 void G4RiGeMuPairProductionModel::MakeSampling
458 {
459 G4double factore = G4Exp(G4Log(emax/emin)/G4
460
461 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
462
463 G4double Z = ZDATPAIR[iz];
464 G4Physics2DVector* pv = new G4Physics2DVec
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 = MaxSecondaryEne
471 /*
472 G4cout << "it= " << it << " E= " << kinE
473 << " " << particle->GetParticleN
474 << " maxE= " << maxPairEnergy <<
475 << " ymin= " << ymin << G4endl;
476 */
477 G4double coef = G4Log(minPairEnergy/kinE
478 G4double ymax = G4Log(maxPairEnergy/kinE
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= "
487 << " mE= " << maxPairEnergy << "
488 << " dy= " << dy << " c= " << co
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*(
500
501 // not multiplied by interval, becau
502 // will be used only for sampling
503 //G4cout << "i= " << i << " x= " <<
504 // << " Egamma= " << ep << G
505 xSec += ep*ComputeDMicroscopicCrossS
506
507 // last bin before the kinematic lim
508 } else if(i == imax) {
509 G4double ep = kinEnergy*G4Exp(coef*(
510 xSec += ep*fac*ComputeDMicroscopicCr
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 }
523
524 //....oooOO0OOooo........oooOO0OOooo........oo
525
526 void G4RiGeMuPairProductionModel::SampleSecond
527
528
529
530
531 {
532 G4double eMass = CLHEP::electron_mass_c2;
533 G4double eMass2 = eMass*eMass;
534
535 // Energy and momentum of the pramary partic
536 G4double kinEnergy = aDynamicParticle->GetKi
537 G4double particleMomentum = aDynamicParticle
538 G4ThreeVector particleMomentumVector = aDyna
539 G4ThreeVector partDirection = aDynamicPartic
540
541 G4double minQ2 = 4.*eMass2;
542 G4double maxQ2 = (kinEnergy - particleMass)*
543 G4double intervalQ2 = G4Log(maxQ2/minQ2);
544
545 // Square invariant of mass of the pair
546 G4double Q2 = minQ2*G4Exp(intervalQ2*randNum
547
548 G4double mingEnergy = std::sqrt(Q2);
549 G4double maxgEnergy = kinEnergy - particleMa
550 G4double intervalgEnergy = maxgEnergy - ming
551
552 // Energy of virtual gamma
553 G4double gEnergy = mingEnergy + intervalgEne
554
555 // Momentum module of the virtual gamma
556 G4double gMomentum = std::sqrt(gEnergy*gEner
557
558 // Energy and momentum module of the outgoin
559 G4double particleFinalEnergy = kinEnergy - g
560 G4double particleFinalMomentum = std::sqrt(p
561 particleMass*particleMass);
562
563 G4double mint3 = 0.;
564 G4double maxt3 = CLHEP::pi;
565 G4double Cmin = std::cos(maxt3);
566 G4double Cmax = std::cos(mint3);
567
568 //G4cout << "------- G4RiGeMuPairProductionM
569 // << kinEnergy << " "
570 // << aDynamicParticle->GetDefinitio
571
572 // select randomly one element constituing t
573 const G4Element* anElement = SelectRandomAto
574
575 // define interval of energy transfer
576 G4double maxPairEnergy = MaxSecondaryEnergyF
577
578 G4double maxEnergy = std::min(tmax, maxPairE
579 G4double minEnergy = std::max(tmin, minPairE
580
581 if (minEnergy >= maxEnergy) { return; }
582
583 //G4cout << "emin= " << minEnergy << " emax=
584 // << " minPair= " << minPairEnergy << " max
585 // << " ymin= " << ymin << " dy= " << dy
586
587 CLHEP::HepRandomEngine* rndmEngine = G4Rando
588
589 G4double coeff = G4Log(minPairEnergy/kinEner
590
591 // compute limits
592 G4double yymin = G4Log(minEnergy/kinEnergy)/
593 G4double yymax = G4Log(maxEnergy/kinEnergy)/
594
595 //G4cout << "yymin= " << yymin << " yymax=
596
597 // units should not be used, bacause table w
598 G4double logTkin = G4Log(kinEnergy/CLHEP::Me
599
600 // sample e-e+ energy, pair energy first
601
602 // select sample table via Z
603 G4int iz1(0), iz2(0);
604 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
605 if(currentZ == ZDATPAIR[iz]) {
606 iz1 = iz2 = iz;
607 break;
608 } else if(currentZ < ZDATPAIR[iz]) {
609 iz2 = iz;
610 if(iz > 0) { iz1 = iz-1; }
611 else { iz1 = iz2; }
612 break;
613 }
614 }
615 if (0 == iz1) { iz1 = iz2 = NZDATPAIR-1; }
616
617 G4double pairEnergy = 0.0;
618 G4int count = 0;
619 //G4cout << "start loop Z1= " << iz1 << " Z2
620 do {
621 ++count;
622 // sampling using only one random number
623 G4double rand = rndmEngine->flat();
624
625 G4double x = FindScaledEnergy(iz1, rand, l
626 if(iz1 != iz2) {
627 G4double x2 = FindScaledEnergy(iz2, rand
628 G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1
629 G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2
630 //G4cout << count << ". x= " << x << "
631 // << " Z1= " << iz1 << " Z2
632 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
633 }
634 //G4cout << "x= " << x << " coeff= " << c
635 pairEnergy = kinEnergy*G4Exp(x*coeff);
636
637 // Loop checking, 30-Oct-2024, Vladimir Iv
638 } while((pairEnergy < minEnergy || pairEnerg
639
640 //G4cout << "## pairEnergy(GeV)= " << pairEn
641 // << " Etot(GeV)= " << totalEnergy/
642 rndmEngine->flatArray(9, randNumbs);
643 G4double phi3 = CLHEP::twopi*randNumbs[0];
644 fAngularGenerator->PhiRotation(partDirection
645
646 G4LorentzVector muF;
647 G4ThreeVector eDirection, pDirection;
648 G4double eEnergy, pEnergy;
649
650 if (randNumbs[7] < W[0]) {
651 G4double A1 = -(Q2 - 2.*kinEnergy*gEnergy)
652 G4double B1 = -(2.*gMomentum*particleMomen
653 G4double tginterval = G4Log((A1 + B1)/(A1
654
655 G4double costg = (-A1 + (A1 - B1)*G4Exp(B1
656 G4double sintg = std::sqrt((1.0 - costg)*(
657 G4double phig = CLHEP::twopi*randNumbs[2]
658 G4double sinpg = std::sin(phig);
659 G4double cospg = std::cos(phig);
660
661 G4ThreeVector dirGamma;
662 dirGamma.set(sintg*cospg, sintg*sinpg, cos
663 G4LorentzVector gFourMomentum(gEnergy, dir
664
665 G4double Ap = particleMomentum*particleMom
666 particleFinalMomentum*particleFinalMomen
667 G4double A = Ap - 2.*particleMomentum*gMom
668 G4double B = 2.*particleMomentum*gMomentum
669 G4double C = 2.*particleFinalMomentum*gMom
670 2.*particleMomentum*particleFinalMomentu
671 G4double absB = std::abs(B);
672 G4double t1interval = (1./(A + C + absB*mi
673 G4double t1 = (-(A + C) + 1./(1./(A + C +
674 G4double sint1 = std::sin(t1);
675 G4double cost1 = std::cos(t1);
676
677 // Ingoing parent particle change
678 G4double Phi = CLHEP::twopi*randNumbs[3];
679 partDirection.set(sint1, 0., cost1);
680 fAngularGenerator->PhiRotation(partDirecti
681 kinEnergy = particleFinalEnergy;
682
683 G4double cost5 = -1. + 2.*randNumbs[6];
684 G4double phi5 = CLHEP::twopi*randNumbs[8];
685
686 G4LorentzVector eFourMomentumMQ = fAngular
687 G4LorentzVector pFourMomentumMQ = fAngular
688
689 G4LorentzVector eFourMomentum = eFourMomen
690 G4LorentzVector pFourMomentum = pFourMomen
691
692 eEnergy = eFourMomentum.t();
693 pEnergy = pFourMomentum.t();
694
695 } else if (randNumbs[7] >= W[0] && randNumbs
696 G4double A3 = Q2 + 2.*gEnergy*particleFina
697 G4double B3 = -2.*gMomentum*particleFinalM
698
699 G4double tQ3interval = G4Log((A3 + B3)/(A3
700 G4double tQMG = (-A3 + (A3 - B3)*G4Exp(B3*
701 G4double phiQP = CLHEP::twopi*randNumbs[2]
702
703 G4double sintQ3 = std::sqrt(1. - tQMG*tQMG
704 G4double cospQP = std::cos(phiQP);
705 G4double sinpQP = std::sin(phiQP);
706
707 G4double Ap = particleMomentum*particleMom
708 particleFinalMomentum*particleFinalMomen
709 G4double A = Ap + 2.*particleFinalMomentum
710 G4double B = -2.*particleMomentum*gMomentu
711 G4double C = -2.*particleMomentum*gMomentu
712
713 G4double absB = std::abs(B);
714 G4double t3interval = (1./(A + C + absB*mi
715 G4double t3 = (-(A + C) + 1./(1./(A + C +
716 G4double sint3 = std::sin(t3);
717 G4double cost3 = std::cos(t3);
718
719 G4double cost = -sint3*sintQ3*cospQP + cos
720 G4double sint = std::sqrt((1. + cost)*(1.
721 G4double cosp = (sintQ3*cospQP*cost3 + sin
722 G4double sinp = sintQ3*sinpQP/sint;
723
724 G4ThreeVector dirGamma;
725 dirGamma.set(sint*cosp, sint*sinp, cost);
726 G4LorentzVector gFourMomentum(gEnergy, dir
727
728 // Ingoing parent particle change
729 G4double Phi = CLHEP::twopi*randNumbs[3];
730 partDirection.set(sint3, 0., cost3);
731 fAngularGenerator->PhiRotation(partDirecti
732 kinEnergy = particleFinalEnergy;
733
734 G4double cost5 = -1. + 2.*randNumbs[6];
735 G4double phi5 = CLHEP::twopi*randNumbs[8];
736
737 G4LorentzVector eFourMomentumMQ = fAngular
738 G4LorentzVector pFourMomentumMQ = fAngular
739
740 G4LorentzVector eFourMomentum = eFourMomen
741 G4LorentzVector pFourMomentum = pFourMomen
742
743 eEnergy = eFourMomentum.t();
744 pEnergy = pFourMomentum.t();
745
746 } else if (randNumbs[7] >= W[1] && randNumbs
747 G4double phi5 = CLHEP::twopi*randNumbs[1];
748 G4double phi6 = CLHEP::twopi*randNumbs[2];
749 G4double muEnergyInterval = kinEnergy - 2.
750 particleFinalEnergy = particleMass + muEne
751 particleFinalMomentum = std::sqrt(particle
752 particleMass*particleMass);
753
754 G4double mineEnergy = eMass;
755 G4double maxeEnergy = kinEnergy - particle
756 G4double eEnergyinterval = maxeEnergy - mi
757 eEnergy = mineEnergy + eEnergyinterval*ran
758
759 G4double cosp3 = 1.;
760 G4double sinp3 = 0.;
761 G4double cosp5 = std::cos(phi5);
762 G4double sinp5 = std::sin(phi5);
763 G4double cosp6 = std::cos(phi6);
764 G4double sinp6 = std::sin(phi6);
765
766 G4double eMomentum = std::sqrt(eEnergy*eEn
767 pEnergy = kinEnergy - particleFinalEnergy
768 G4double pMomentum = std::sqrt(pEnergy*pEn
769
770 G4double A3 = -2.*particleMass*particleMas
771 G4double B3 = -2.*particleMomentum*particl
772 G4double cost3interval = G4Log((A3 + B3*Cm
773 G4double expanCost3r6 = G4Exp(B3*cost3inte
774 G4double cost3 = A3*(expanCost3r6 - 1.)/B3
775 G4double sint3 = std::sqrt((1. - cost3)*(1
776
777 partDirection.set(sint3, 0., cost3);
778
779 G4ThreeVector muFinalMomentumVector;
780 muFinalMomentumVector.set(particleFinalMom
781
782 G4LorentzVector muFourMomentum(particleMom
783 G4LorentzVector muFinalFourMomentum(partic
784 G4LorentzVector auxVec1 = muFourMomentum -
785 G4double A5 = auxVec1.mag2() - 2.*eEnergy*
786 2.*particleMomentumVector[2]*eMomentum -
787 G4double B5 = -2.*particleFinalMomentum*eM
788 G4double absA5 = std::abs(A5);
789 G4double absB5 = std::abs(B5);
790 G4double mint5 = 0.;
791 G4double maxt5 = CLHEP::pi;
792 G4double t5interval = G4Log((absA5 + absB5
793 G4double argexp = absB5*t5interval*randNum
794 G4double t5 = -absA5/absB5 + G4Exp(argexp)
795 G4double sint5 = std::sin(t5);
796 G4double cost5 = std::cos(t5);
797
798 eDirection.set(sint5*cosp5, sint5*sinp5, c
799 G4ThreeVector eMomentumVector = eMomentum*
800
801 G4ThreeVector auxVec2 = particleMomentumVe
802 G4double p1mp3mp52 = auxVec2.dot(auxVec2);
803 G4double Bp = particleFinalMomentum*(sint3
804 eMomentum*(sint5*cosp5*cosp6 + sint5*sin
805 G4double Cp = -particleMomentum + particle
806 G4double A6 = p1mp3mp52 + pMomentum*pMomen
807 G4double B6 = 2.*pMomentum*Bp;
808 G4double C6 = 2.*pMomentum*Cp;
809 G4double mint6 = 0.;
810 G4double maxt6 = CLHEP::pi;
811 G4double absA6C6 = std::abs(A6 + C6);
812 G4double absB6 = std::abs(B6);
813 G4double t6interval = (1./(absA6C6 + absB6
814 G4double t6 = (-absA6C6 + 1./(1./(absA6C6
815 G4double sint6 = std::sin(t6);
816 G4double cost6 = std::cos(t6);
817
818 pDirection.set(sint6*cosp6, sint6*sinp6, c
819
820 } else {
821 G4double phi6 = CLHEP::twopi*randNumbs[1];
822 G4double phi5 = CLHEP::twopi*randNumbs[2];
823 G4double muFinalEnergyinterval = kinEnergy
824 particleFinalEnergy = particleMass + muFin
825 particleFinalMomentum = std::sqrt(particle
826 particleMass*particleMass);
827
828 G4double maxpEnergy = kinEnergy - particle
829 G4double pEnergyinterval = maxpEnergy - eM
830 pEnergy = eMass + pEnergyinterval*randNumb
831
832 G4double cosp3 = 1.;
833 G4double sinp3 = 0.;
834 G4double cosp5 = std::cos(phi5);
835 G4double sinp5 = std::sin(phi5);
836 G4double cosp6 = std::cos(phi6);
837 G4double sinp6 = std::sin(phi6);
838
839 G4double pMomentum = std::sqrt(pEnergy*pEn
840 eEnergy = kinEnergy - particleFinalEnergy
841 G4double eMomentum = std::sqrt(eEnergy*eEn
842
843 G4double A3 = -2.*particleMass*particleMas
844 G4double B3 = -2.*particleMomentum*particl
845 G4double cost3interval = G4Log((A3 + B3*Cm
846 G4double expanCost3r6 = G4Exp(B3*cost3inte
847 G4double cost3 = A3*(expanCost3r6 - 1.)/B3
848 G4double sint3 = std::sqrt((1. - cost3)*(1
849
850 partDirection.set(sint3*cosp3, sint3*sinp3
851
852 G4ThreeVector muFinalMomentumVector;
853 muFinalMomentumVector.set(particleFinalMom
854 particleFinalMomentum*sint3*sinp3,
855 particleFinalMomentum*cost3);
856
857 G4LorentzVector muFourMomentum(particleMom
858 G4LorentzVector muFinalFourMomentum(partic
859 G4LorentzVector auxVec1 = muFourMomentum -
860 G4double A6 = auxVec1.mag2() -
861 2.*pEnergy*(kinEnergy - particleFinalEne
862 2.*particleFinalMomentum*pMomentum*cost3
863 G4double B6 = -2.*particleFinalMomentum*pM
864 G4double absA6 = std::abs(A6);
865 G4double absB6 = std::abs(B6);
866 G4double mint6 = 0.;
867 G4double maxt6 = CLHEP::pi;
868 G4double t6interval = G4Log((absA6 + absB6
869 G4double argexp = absB6*t6interval*randNum
870 G4double t6 = -absA6/absB6 + G4Exp(argexp)
871 G4double sint6 = std::sin(t6);
872 G4double cost6 = std::cos(t6);
873
874 pDirection.set(sint6*cosp6, sint6*sinp6, c
875 G4ThreeVector pMomentumVector = pMomentum*
876
877 G4ThreeVector auxVec2 = particleMomentumVe
878 G4double p1mp3mp62 = auxVec2.dot(auxVec2);
879 G4double Bp = particleFinalMomentum*(sint3
880 pMomentum*(sint6*cosp6*cosp5 + sint6*sin
881 G4double Cp = -particleMomentum + particle
882 G4double A5 = p1mp3mp62 + eMomentum*eMomen
883 G4double B5 = 2.*eMomentum*Bp;
884 G4double C5 = 2.*eMomentum*Cp;
885 G4double mint5 = 0.;
886 G4double maxt5 = CLHEP::pi;
887 G4double absA5C5 = std::abs(A5 + C5);
888 G4double absB5 = std::abs(B5);
889 G4double t5interval = (1./(absA5C5 + absB5
890 G4double t5 = (-absA5C5 + 1./(1./(absA5C5
891 G4double sint5 = std::sin(t5);
892 G4double cost5 = std::cos(t5);
893
894 eDirection.set(sint5*cosp5, sint5*sinp5, c
895 }
896
897 fAngularGenerator->Sample5DPairDirections(aD
898 gEnergy, Q2, gMomentum,
899 particleFinalMomentum,
900 particleFinalEnergy,
901 randNumbs, W);
902
903 // create G4DynamicParticle object for e+e-
904 auto aParticle1 = new G4DynamicParticle(theE
905 auto aParticle2 = new G4DynamicParticle(theP
906
907 // Fill output vector
908 vdp->push_back(aParticle1);
909 vdp->push_back(aParticle2);
910
911 // if energy transfer is higher than thresho
912 // then stop tracking the primary particle a
913 if (pairEnergy > SecondaryThreshold()) {
914 fParticleChange->ProposeTrackStatus(fStopA
915 fParticleChange->SetProposedKineticEnergy(
916 auto newdp = new G4DynamicParticle(particl
917 vdp->push_back(newdp);
918 } else { // continue tracking the primary e-
919 fParticleChange->SetProposedMomentumDirect
920 G4double ekin = std::max(muF.e() - particl
921 fParticleChange->SetProposedKineticEnergy(
922 }
923 //G4cout << "-- G4RiGeMuPairProductionModel:
924 }
925
926 //....oooOO0OOooo........oooOO0OOooo........oo
927
928 G4double
929 G4RiGeMuPairProductionModel::FindScaledEnergy(
930 G4double logTkin,
931 G4double yymin, G4double yymax
932 {
933 G4double res = yymin;
934 G4Physics2DVector* pv = fElementData->GetEle
935 if (nullptr != pv) {
936 G4double pmin = pv->Value(yymin, logTkin);
937 G4double pmax = pv->Value(yymax, logTkin);
938 G4double p0 = pv->Value(0.0, logTkin);
939 if(p0 <= 0.0) { DataCorrupted(ZDATPAIR[iz]
940 else { res = pv->FindLinearX((pmin + rand*
941 } else {
942 DataCorrupted(ZDATPAIR[iz], logTkin);
943 }
944 return res;
945 }
946
947 //....oooOO0OOooo........oooOO0OOooo........oo
948
949 void G4RiGeMuPairProductionModel::DataCorrupte
950 {
951 G4ExceptionDescription ed;
952 ed << "G4ElementData is not properly initial
953 << " Ekin(MeV)= " << G4Exp(logTkin)
954 << " IsMasterThread= " << IsMaster()
955 << " Model " << GetName();
956 G4Exception("G4RiGeMuPairProductionModel::()
957 }
958
959 //....oooOO0OOooo........oooOO0OOooo........oo
960
961 void G4RiGeMuPairProductionModel::StoreTables(
962 {
963 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
964 G4int Z = ZDATPAIR[iz];
965 G4Physics2DVector* pv = fElementData->GetE
966 if(nullptr == pv) {
967 DataCorrupted(Z, 1.0);
968 return;
969 }
970 std::ostringstream ss;
971 ss << "mupair/" << particle->GetParticleNa
972 std::ofstream outfile(ss.str());
973 pv->Store(outfile);
974 }
975 }
976
977 //....oooOO0OOooo........oooOO0OOooo........oo
978
979 G4bool G4RiGeMuPairProductionModel::RetrieveTa
980 {
981 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
982 G4double Z = ZDATPAIR[iz];
983 G4Physics2DVector* pv = new G4Physics2DVec
984 std::ostringstream ss;
985 ss << G4EmParameters::Instance()->GetDirLE
986 << particle->GetParticleName() << Z <<
987 std::ifstream infile(ss.str(), std::ios::i
988 if(!pv->Retrieve(infile)) {
989 delete pv;
990 return false;
991 }
992 fElementData->InitialiseForElement(iz, pv)
993 }
994 return true;
995 }
996
997 //....oooOO0OOooo........oooOO0OOooo........oo
998