| 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: G4MuPairProductionModel
33 //
34 // Author: Vladimir Ivanchenko on base
35 //
36 // Creation date: 24.06.2002
37 //
38 // Modifications:
39 //
40 // 04-12-02 Change G4DynamicParticle construct
41 // 23-12-02 Change interface in order to move
42 // 24-01-03 Fix for compounds (V.Ivanchenko)
43 // 27-01-03 Make models region aware (V.Ivanch
44 // 13-02-03 Add model (V.Ivanchenko)
45 // 06-06-03 Fix in cross section calculation f
46 // 20-10-03 2*xi in ComputeDDMicroscopicCrossS
47 // 8 integration points in ComputeDMi
48 // 12-01-04 Take min cut of e- and e+ not its
49 // 10-02-04 Update parameterisation using R.Ko
50 // 28-04-04 For complex materials repeat calcu
51 // material (V.Ivanchenko)
52 // 01-11-04 Fix bug inside ComputeDMicroscopic
53 // 08-04-05 Major optimisation of internal int
54 // 03-08-05 Add SetParticle method (V.Ivantche
55 // 23-10-05 Add protection in sampling of e+e-
56 // low cuts (V.Ivantchenko)
57 // 13-02-06 Add ComputeCrossSectionPerAtom (mm
58 // 24-04-07 Add protection in SelectRandomAtom
59 // 12-05-06 Updated sampling (use cut) in Sele
60 // 11-10-07 Add ignoreCut flag (V.Ivanchenko)
61
62 //
63 // Class Description:
64 //
65 //
66 // -------------------------------------------
67 //
68 //....oooOO0OOooo........oooOO0OOooo........oo
69 //....oooOO0OOooo........oooOO0OOooo........oo
70
71 #include "G4MuPairProductionModel.hh"
72 #include "G4PhysicalConstants.hh"
73 #include "G4SystemOfUnits.hh"
74 #include "G4EmParameters.hh"
75 #include "G4Electron.hh"
76 #include "G4Positron.hh"
77 #include "G4MuonMinus.hh"
78 #include "G4MuonPlus.hh"
79 #include "Randomize.hh"
80 #include "G4Material.hh"
81 #include "G4Element.hh"
82 #include "G4ElementVector.hh"
83 #include "G4ElementDataRegistry.hh"
84 #include "G4ProductionCutsTable.hh"
85 #include "G4ParticleChangeForLoss.hh"
86 #include "G4ModifiedMephi.hh"
87 #include "G4Log.hh"
88 #include "G4Exp.hh"
89 #include "G4AutoLock.hh"
90
91 #include <iostream>
92 #include <fstream>
93
94 //....oooOO0OOooo........oooOO0OOooo........oo
95
96 const G4int G4MuPairProductionModel::ZDATPAIR[
97
98 const G4double G4MuPairProductionModel::xgi[]
99 0.0198550717512320, 0.1016667612931865, 0.
100 0.5917173212478250, 0.7627662049581645, 0.
101 };
102
103 const G4double G4MuPairProductionModel::wgi[]
104 0.0506142681451880, 0.1111905172266870, 0.
105 0.1813418916891810, 0.1568533229389435, 0.
106 };
107
108 namespace
109 {
110 G4Mutex theMuPairMutex = G4MUTEX_INITIALIZER
111
112 const G4double ak1 = 6.9;
113 const G4double ak2 = 1.0;
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oo
117
118 G4MuPairProductionModel::G4MuPairProductionMod
119
120 : G4VEmModel(nam),
121 factorForCross(CLHEP::fine_structure_const*C
122 CLHEP::classic_electr_radius*CLHEP::class
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 {
128 nist = G4NistManager::Instance();
129
130 theElectron = G4Electron::Electron();
131 thePositron = G4Positron::Positron();
132
133 if(nullptr != p) {
134 SetParticle(p);
135 lowestKinEnergy = std::max(lowestKinEnergy
136 }
137 emin = lowestKinEnergy;
138 emax = emin*10000.;
139 SetAngularDistribution(new G4ModifiedMephi()
140 }
141
142 //....oooOO0OOooo........oooOO0OOooo........oo
143
144 G4double G4MuPairProductionModel::MinPrimaryEn
145
146
147 {
148 return std::max(lowestKinEnergy, cut);
149 }
150
151 //....oooOO0OOooo........oooOO0OOooo........oo
152
153 void G4MuPairProductionModel::Initialise(const
154 const
155 {
156 SetParticle(p);
157
158 if (nullptr == fParticleChange) {
159 fParticleChange = GetParticleChangeForLoss
160
161 // define scale of internal table for each
162 if (0 == nbine) {
163 emin = std::max(lowestKinEnergy, LowEner
164 emax = std::max(HighEnergyLimit(), emin*
165 nbine = std::size_t(nYBinPerDecade*std::
166 if(nbine < 3) { nbine = 3; }
167
168 ymin = G4Log(minPairEnergy/emin);
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 s
192 if(lowestKinEnergy >= HighEnergyLimit()) { r
193
194 if (p == particle) {
195 fElementData =
196 G4ElementDataRegistry::Instance()->GetEl
197 if (nullptr == fElementData) {
198 G4AutoLock l(&theMuPairMutex);
199 fElementData =
200 G4ElementDataRegistry::Instance()->GetElemen
201 if (nullptr == fElementData) {
202 fElementData = new G4ElementData(NZDATPAIR);
203 fElementData->SetName(dataName);
204 }
205 G4bool useDataFile = G4EmParameters::Ins
206 if (useDataFile) { useDataFile = Retrie
207 if (!useDataFile) { MakeSamplingTables()
208 if (fTableToFile) { StoreTables(); }
209 l.unlock();
210 }
211 if (IsMaster()) {
212 InitialiseElementSelectors(p, cuts);
213 }
214 }
215 }
216
217 //....oooOO0OOooo........oooOO0OOooo........oo
218
219 void G4MuPairProductionModel::InitialiseLocal(
220
221 {
222 if(p == particle && lowestKinEnergy < HighEn
223 SetElementSelectors(masterModel->GetElemen
224 }
225 }
226
227 //....oooOO0OOooo........oooOO0OOooo........oo
228
229 G4double G4MuPairProductionModel::ComputeDEDXP
230
231
232
233
234 {
235 G4double dedx = 0.0;
236 if (cutEnergy <= minPairEnergy || kineticEne
237 { return dedx; }
238
239 const G4ElementVector* theElementVector = ma
240 const G4double* theAtomicNumDensityVector =
241 material->G
242
243 // loop for elements in the material
244 for (std::size_t i=0; i<material->GetNumberO
245 G4double Z = (*theElementVector)[i]->GetZ
246 G4double tmax = MaxSecondaryEnergyForElem
247 G4double loss = ComputMuPairLoss(Z, kinet
248 dedx += loss*theAtomicNumDensityVector[i]
249 }
250 dedx = std::max(dedx, 0.0);
251 return dedx;
252 }
253
254 //....oooOO0OOooo........oooOO0OOooo........oo
255
256 G4double G4MuPairProductionModel::ComputMuPair
257
258
259
260 {
261 G4double loss = 0.0;
262
263 G4double cut = std::min(cutEnergy, tmax);
264 if(cut <= minPairEnergy) { return loss; }
265
266 // calculate the rectricted loss
267 // numerical integration in log(PairEnergy)
268 G4double aaa = G4Log(minPairEnergy);
269 G4double bbb = G4Log(cut);
270
271 G4int kkk = std::min(std::max(G4lrint((bbb-a
272 G4double hhh = (bbb-aaa)/kkk;
273 G4double x = aaa;
274
275 for (G4int l=0 ; l<kkk; ++l) {
276 for (G4int ll=0; ll<NINTPAIR; ++ll) {
277 G4double ep = G4Exp(x+xgi[ll]*hhh);
278 loss += wgi[ll]*ep*ep*ComputeDMicroscopi
279 }
280 x += hhh;
281 }
282 loss *= hhh;
283 loss = std::max(loss, 0.0);
284 return loss;
285 }
286
287 //....oooOO0OOooo........oooOO0OOooo........oo
288
289 G4double G4MuPairProductionModel::ComputeMicro
290 G4d
291 G4d
292 G4d
293 {
294 G4double cross = 0.;
295 G4double tmax = MaxSecondaryEnergyForElement
296 G4double cut = std::max(cutEnergy, minPairE
297 if (tmax <= cut) { return cross; }
298
299 G4double aaa = G4Log(cut);
300 G4double bbb = G4Log(tmax);
301 G4int kkk = std::min(std::max(G4lrint((bbb-a
302
303 G4double hhh = (bbb-aaa)/(kkk);
304 G4double x = aaa;
305
306 for (G4int l=0; l<kkk; ++l) {
307 for (G4int i=0; i<NINTPAIR; ++i) {
308 G4double ep = G4Exp(x + xgi[i]*hhh);
309 cross += ep*wgi[i]*ComputeDMicroscopicCr
310 }
311 x += hhh;
312 }
313
314 cross *= hhh;
315 cross = std::max(cross, 0.0);
316 return cross;
317 }
318
319 //....oooOO0OOooo........oooOO0OOooo........oo
320
321 G4double G4MuPairProductionModel::ComputeDMicr
322 G4d
323 G4d
324 G4d
325 // Calculates the differential (D) microscopi
326 // using the cross section formula of R.P. Kok
327 // Code modified by R.P. Kokoulin, V.N. Ivanch
328 {
329 static const G4double bbbtf= 183. ;
330 static const G4double bbbh = 202.4 ;
331 static const G4double g1tf = 1.95e-5 ;
332 static const G4double g2tf = 5.3e-5 ;
333 static const G4double g1h = 4.4e-5 ;
334 static const G4double g2h = 4.8e-5 ;
335
336 if (pairEnergy <= minPairEnergy)
337 return 0.0;
338
339 G4double totalEnergy = tkin + particleMass;
340 G4double residEnergy = totalEnergy - pairEn
341
342 if (residEnergy <= 0.75*sqrte*z13*particleMa
343 return 0.0;
344
345 G4double a0 = 1.0 / (totalEnergy * residEner
346 G4double alf = 4.0 * electron_mass_c2 / pair
347 G4double rt = std::sqrt(1.0 - alf);
348 G4double delta = 6.0 * particleMass * partic
349 G4double tmnexp = alf/(1.0 + rt) + delta*rt;
350
351 if(tmnexp >= 1.0) { return 0.0; }
352
353 G4double tmn = G4Log(tmnexp);
354
355 G4double massratio = particleMass/CLHEP::ele
356 G4double massratio2 = massratio*massratio;
357 G4double inv_massratio2 = 1.0 / massratio2;
358
359 // zeta calculation
360 G4double bbb,g1,g2;
361 if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 =
362 else { bbb = bbbtf; g1 = g1tf; g2 =
363
364 G4double zeta = 0.0;
365 G4double z1exp = totalEnergy / (particleMass
366
367 // 35.221047195922 is the root of zeta1(x) =
368 // condition below is the same as zeta1 > 0.
369 if (z1exp > 35.221047195922)
370 {
371 G4double z2exp = totalEnergy / (particleMa
372 zeta = (0.073 * G4Log(z1exp) - 0.26) / (0.
373 }
374
375 G4double z2 = Z*(Z+zeta);
376 G4double screen0 = 2.*electron_mass_c2*sqrte
377 G4double beta = 0.5*pairEnergy*pairEnergy*a0
378 G4double xi0 = 0.5*massratio2*beta;
379
380 // Gaussian integration in ln(1-ro) ( with 8
381 G4double rho[NINTPAIR];
382 G4double rho2[NINTPAIR];
383 G4double xi[NINTPAIR];
384 G4double xi1[NINTPAIR];
385 G4double xii[NINTPAIR];
386
387 for (G4int i = 0; i < NINTPAIR; ++i)
388 {
389 rho[i] = G4Exp(tmn*xgi[i]) - 1.0; // rho =
390 rho2[i] = rho[i] * rho[i];
391 xi[i] = xi0*(1.0-rho2[i]);
392 xi1[i] = 1.0 + xi[i];
393 xii[i] = 1.0 / xi[i];
394 }
395
396 G4double ye1[NINTPAIR];
397 G4double ym1[NINTPAIR];
398
399 G4double b40 = 4.0 * beta;
400 G4double b62 = 6.0 * beta + 2.0;
401
402 for (G4int i = 0; i < NINTPAIR; ++i)
403 {
404 G4double yeu = (b40 + 5.0) + (b40 - 1.0) *
405 G4double yed = b62*G4Log(3.0 + xii[i]) + (
406
407 G4double ymu = b62 * (1.0 + rho2[i]) + 6.0
408 G4double ymd = (b40 + 3.0)*(1.0 + rho2[i])
409 + 2.0 - 3.0 * rho2[i];
410
411 ye1[i] = 1.0 + yeu / yed;
412 ym1[i] = 1.0 + ymu / ymd;
413 }
414
415 G4double be[NINTPAIR];
416 G4double bm[NINTPAIR];
417
418 for(G4int i = 0; i < NINTPAIR; ++i) {
419 if(xi[i] <= 1000.0) {
420 be[i] = ((2.0 + rho2[i])*(1.0 + beta) +
421 xi[i]*(3.0 + rho2[i]))*G4Log(1.0 + xi
422 (1.0 - rho2[i] - beta)/xi1[i] - (3.0 + rho2[
423 } else {
424 be[i] = 0.5*(3.0 - rho2[i] + 2.0*beta*(1
425 }
426
427 if(xi[i] >= 0.001) {
428 G4double a10 = (1.0 + 2.0 * beta) * (1.0
429 bm[i] = ((1.0 + rho2[i])*(1.0 + 1.5 * be
430 xi[i] * (1.0 - rho2[i] - beta)
431 } else {
432 bm[i] = 0.5*(5.0 - rho2[i] + beta * (3.0
433 }
434 }
435
436 G4double sum = 0.0;
437
438 for (G4int i = 0; i < NINTPAIR; ++i) {
439 G4double screen = screen0*xi1[i]/(1.0 - rh
440 G4double ale = G4Log(bbb/z13*std::sqrt(xi1
441 G4double cre = 0.5*G4Log(1. + 2.25*z23*xi1
442
443 G4double fe = (ale-cre)*be[i];
444 fe = std::max(fe, 0.0);
445
446 G4double alm_crm = G4Log(bbb*massratio/(1.
447 G4double fm = std::max(alm_crm*bm[i], 0.0)
448
449 sum += wgi[i]*(1.0 + rho[i])*(fe + fm);
450 }
451
452 return -tmn*sum*factorForCross*z2*residEnerg
453 }
454
455 //....oooOO0OOooo........oooOO0OOooo........oo
456
457 G4double G4MuPairProductionModel::ComputeCross
458 con
459
460
461
462
463 {
464 G4double cross = 0.0;
465 if (kineticEnergy <= lowestKinEnergy) { retu
466
467 G4double maxPairEnergy = MaxSecondaryEnergyF
468 G4double tmax = std::min(maxEnergy, maxPairE
469 G4double cut = std::max(cutEnergy, minPairE
470 if (cut >= tmax) { return cross; }
471
472 cross = ComputeMicroscopicCrossSection(kinet
473 if(tmax < kineticEnergy) {
474 cross -= ComputeMicroscopicCrossSection(ki
475 }
476 return cross;
477 }
478
479 //....oooOO0OOooo........oooOO0OOooo........oo
480
481 void G4MuPairProductionModel::MakeSamplingTabl
482 {
483 G4double factore = G4Exp(G4Log(emax/emin)/G4
484
485 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
486
487 G4double Z = ZDATPAIR[iz];
488 G4Physics2DVector* pv = new G4Physics2DVec
489 G4double kinEnergy = emin;
490
491 for (std::size_t it=0; it<=nbine; ++it) {
492
493 pv->PutY(it, G4Log(kinEnergy/CLHEP::MeV)
494 G4double maxPairEnergy = MaxSecondaryEne
495 /*
496 G4cout << "it= " << it << " E= " << kinE
497 << " " << particle->GetParticleN
498 << " maxE= " << maxPairEnergy <<
499 << " ymin= " << ymin << G4endl;
500 */
501 G4double coef = G4Log(minPairEnergy/kinE
502 G4double ymax = G4Log(maxPairEnergy/kinE
503 G4double fac = (ymax - ymin)/dy;
504 std::size_t imax = (std::size_t)fac;
505 fac -= (G4double)imax;
506
507 G4double xSec = 0.0;
508 G4double x = ymin;
509 /*
510 G4cout << "Z= " << currentZ << " z13= "
511 << " mE= " << maxPairEnergy << "
512 << " dy= " << dy << " c= " << co
513 */
514 // start from zero
515 pv->PutValue(0, it, 0.0);
516 if(0 == it) { pv->PutX(nbiny, 0.0); }
517
518 for (std::size_t i=0; i<nbiny; ++i) {
519
520 if(0 == it) { pv->PutX(i, x); }
521
522 if(i < imax) {
523 G4double ep = kinEnergy*G4Exp(coef*(
524
525 // not multiplied by interval, becau
526 // will be used only for sampling
527 //G4cout << "i= " << i << " x= " <<
528 // << " Egamma= " << ep << G
529 xSec += ep*ComputeDMicroscopicCrossS
530
531 // last bin before the kinematic lim
532 } else if(i == imax) {
533 G4double ep = kinEnergy*G4Exp(coef*(
534 xSec += ep*fac*ComputeDMicroscopicCr
535 }
536 pv->PutValue(i + 1, it, xSec);
537 x += dy;
538 }
539 kinEnergy *= factore;
540
541 // to avoid precision lost
542 if(it+1 == nbine) { kinEnergy = emax; }
543 }
544 fElementData->InitialiseForElement(iz, pv)
545 }
546 }
547
548 //....oooOO0OOooo........oooOO0OOooo........oo
549
550 void G4MuPairProductionModel::SampleSecondarie
551 std::vector<G4Dy
552 const G4Material
553 const G4DynamicP
554 G4double tmin,
555 G4double tmax)
556 {
557 G4double kinEnergy = aDynamicParticle->GetKi
558 //G4cout << "------- G4MuPairProductionModel
559 // << kinEnergy << " "
560 // << aDynamicParticle->GetDefinitio
561 G4double totalEnergy = kinEnergy + particl
562 G4double totalMomentum =
563 std::sqrt(kinEnergy*(kinEnergy + 2.0*parti
564
565 G4ThreeVector partDirection = aDynamicPartic
566
567 // select randomly one element constituing t
568 const G4Element* anElement = SelectRandomAto
569
570 // define interval of energy transfer
571 G4double maxPairEnergy = MaxSecondaryEnergyF
572
573 G4double maxEnergy = std::min(tmax, maxPairE
574 G4double minEnergy = std::max(tmin, minPairE
575
576 if (minEnergy >= maxEnergy) { return; }
577 //G4cout << "emin= " << minEnergy << " emax=
578 // << " minPair= " << minPairEnergy << " max
579 // << " ymin= " << ymin << " dy= " << dy
580
581 G4double coeff = G4Log(minPairEnergy/kinEner
582
583 // compute limits
584 G4double yymin = G4Log(minEnergy/kinEnergy)/
585 G4double yymax = G4Log(maxEnergy/kinEnergy)/
586
587 //G4cout << "yymin= " << yymin << " yymax=
588
589 // units should not be used, bacause table w
590 G4double logTkin = G4Log(kinEnergy/CLHEP::Me
591
592 // sample e-e+ energy, pair energy first
593
594 // select sample table via Z
595 G4int iz1(0), iz2(0);
596 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
597 if(currentZ == ZDATPAIR[iz]) {
598 iz1 = iz2 = iz;
599 break;
600 } else if(currentZ < ZDATPAIR[iz]) {
601 iz2 = iz;
602 if(iz > 0) { iz1 = iz-1; }
603 else { iz1 = iz2; }
604 break;
605 }
606 }
607 if (0 == iz1) { iz1 = iz2 = NZDATPAIR-1; }
608
609 G4double pairEnergy = 0.0;
610 G4int count = 0;
611 //G4cout << "start loop Z1= " << iz1 << " Z2
612 do {
613 ++count;
614 // sampling using only one random number
615 G4double rand = G4UniformRand();
616
617 G4double x = FindScaledEnergy(iz1, rand, l
618 if(iz1 != iz2) {
619 G4double x2 = FindScaledEnergy(iz2, rand
620 G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1
621 G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2
622 //G4cout << count << ". x= " << x << "
623 // << " Z1= " << iz1 << " Z2
624 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
625 }
626 //G4cout << "x= " << x << " coeff= " << c
627 pairEnergy = kinEnergy*G4Exp(x*coeff);
628
629 // Loop checking, 03-Aug-2015, Vladimir Iv
630 } while((pairEnergy < minEnergy || pairEnerg
631
632 //G4cout << "## pairEnergy(GeV)= " << pairEn
633 // << " Etot(GeV)= " << totalEnergy/
634
635 // sample r=(E+-E-)/pairEnergy ( uniformly
636 G4double rmax =
637 (1.-6.*particleMass*particleMass/(totalEne
638 *std::sqrt(1.-minPairEnergy/pairEnergy);
639 G4double r = rmax * (-1.+2.*G4UniformRand())
640
641 // compute energies from pairEnergy,r
642 G4double eEnergy = (1.-r)*pairEnergy*0.5;
643 G4double pEnergy = pairEnergy - eEnergy;
644
645 // Sample angles
646 G4ThreeVector eDirection, pDirection;
647 //
648 GetAngularDistribution()->SamplePairDirectio
649
650
651 // create G4DynamicParticle object for e+e-
652 eEnergy = std::max(eEnergy - CLHEP::electron
653 pEnergy = std::max(pEnergy - CLHEP::electron
654 auto aParticle1 = new G4DynamicParticle(theE
655 auto aParticle2 = new G4DynamicParticle(theP
656 // Fill output vector
657 vdp->push_back(aParticle1);
658 vdp->push_back(aParticle2);
659
660 // primary change
661 kinEnergy -= pairEnergy;
662 partDirection *= totalMomentum;
663 partDirection -= (aParticle1->GetMomentum()
664 partDirection = partDirection.unit();
665
666 // if energy transfer is higher than thresho
667 // then stop tracking the primary particle a
668 if (pairEnergy > SecondaryThreshold()) {
669 fParticleChange->ProposeTrackStatus(fStopA
670 fParticleChange->SetProposedKineticEnergy(
671 auto newdp = new G4DynamicParticle(particl
672 vdp->push_back(newdp);
673 } else { // continue tracking the primary e-
674 fParticleChange->SetProposedMomentumDirect
675 fParticleChange->SetProposedKineticEnergy(
676 }
677 //G4cout << "-- G4MuPairProductionModel::Sam
678 }
679
680 //....oooOO0OOooo........oooOO0OOooo........oo
681
682 G4double
683 G4MuPairProductionModel::FindScaledEnergy(G4in
684 G4double logTkin,
685 G4double yymin, G4double yymax)
686 {
687 G4double res = yymin;
688 G4Physics2DVector* pv = fElementData->GetEle
689 if (nullptr != pv) {
690 G4double pmin = pv->Value(yymin, logTkin);
691 G4double pmax = pv->Value(yymax, logTkin);
692 G4double p0 = pv->Value(0.0, logTkin);
693 if(p0 <= 0.0) { DataCorrupted(ZDATPAIR[iz]
694 else { res = pv->FindLinearX((pmin + rand*
695 } else {
696 DataCorrupted(ZDATPAIR[iz], logTkin);
697 }
698 return res;
699 }
700
701 //....oooOO0OOooo........oooOO0OOooo........oo
702
703 void G4MuPairProductionModel::DataCorrupted(G4
704 {
705 G4ExceptionDescription ed;
706 ed << "G4ElementData is not properly initial
707 << " Ekin(MeV)= " << G4Exp(logTkin)
708 << " IsMasterThread= " << IsMaster()
709 << " Model " << GetName();
710 G4Exception("G4MuPairProductionModel::()", "
711 }
712
713 //....oooOO0OOooo........oooOO0OOooo........oo
714
715 void G4MuPairProductionModel::StoreTables() co
716 {
717 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
718 G4int Z = ZDATPAIR[iz];
719 G4Physics2DVector* pv = fElementData->GetE
720 if(nullptr == pv) {
721 DataCorrupted(Z, 1.0);
722 return;
723 }
724 std::ostringstream ss;
725 ss << "mupair/" << particle->GetParticleNa
726 std::ofstream outfile(ss.str());
727 pv->Store(outfile);
728 }
729 }
730
731 //....oooOO0OOooo........oooOO0OOooo........oo
732
733 G4bool G4MuPairProductionModel::RetrieveTables
734 {
735 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
736 G4double Z = ZDATPAIR[iz];
737 G4Physics2DVector* pv = new G4Physics2DVec
738 std::ostringstream ss;
739 ss << G4EmParameters::Instance()->GetDirLE
740 << particle->GetParticleName() << Z <<
741 std::ifstream infile(ss.str(), std::ios::i
742 if(!pv->Retrieve(infile)) {
743 delete pv;
744 return false;
745 }
746 fElementData->InitialiseForElement(iz, pv)
747 }
748 return true;
749 }
750
751 //....oooOO0OOooo........oooOO0OOooo........oo
752