| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4RiGeMuPairProductionModel
33 //
34 // Authors: Girardo Depaola & Ricardo Pacheco
35 //
36 // Creation date: 29.10.2024
37 //
38 //
39 // -------------------------------------------------------------------
40 //
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
68
69 const G4int G4RiGeMuPairProductionModel::ZDATPAIR[] = {1, 4, 13, 29, 92};
70
71 const G4double G4RiGeMuPairProductionModel::xgi[] = {
72 0.0198550717512320, 0.1016667612931865, 0.2372337950418355, 0.4082826787521750,
73 0.5917173212478250, 0.7627662049581645, 0.8983332387068135, 0.9801449282487680
74 };
75
76 const G4double G4RiGeMuPairProductionModel::wgi[] = {
77 0.0506142681451880, 0.1111905172266870, 0.1568533229389435, 0.1813418916891810,
78 0.1813418916891810, 0.1568533229389435, 0.1111905172266870, 0.0506142681451880
79 };
80
81 namespace
82 {
83 G4Mutex theRiGeMuPairMutex = G4MUTEX_INITIALIZER;
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........oooOO0OOooo........oooOO0OOooo....
93
94 G4RiGeMuPairProductionModel::G4RiGeMuPairProductionModel(const G4ParticleDefinition* p)
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 {
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, p->GetPDGMass()*8.0);
111 }
112 emin = lowestKinEnergy;
113 emax = emin*10000.;
114 fAngularGenerator = new G4RiGeAngularGenerator();
115 SetAngularDistribution(fAngularGenerator);
116 for (G4int i=0; i<9; ++i) { randNumbs[i] = 0.0; }
117 }
118
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120
121 G4double
122 G4RiGeMuPairProductionModel::MinPrimaryEnergy(const G4Material*,
123 const G4ParticleDefinition*,
124 G4double cut)
125 {
126 return std::max(lowestKinEnergy, cut);
127 }
128
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130
131 void G4RiGeMuPairProductionModel::Initialise(const G4ParticleDefinition* p,
132 const G4DataVector& cuts)
133 {
134 SetParticle(p);
135
136 if (nullptr == fParticleChange) {
137 fParticleChange = GetParticleChangeForLoss();
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
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 should not work
170 if(lowestKinEnergy >= HighEnergyLimit()) { return; }
171
172 if (p == particle) {
173 auto data = G4ElementDataRegistry::Instance();
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 }
182 G4bool useDataFile = G4EmParameters::Instance()->RetrieveMuDataFromFile();
183 if (useDataFile) { useDataFile = RetrieveTables(); }
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........oooOO0OOooo........oooOO0OOooo....
195
196 void G4RiGeMuPairProductionModel::InitialiseLocal(const G4ParticleDefinition* p,
197 G4VEmModel* masterModel)
198 {
199 if(p == particle && lowestKinEnergy < HighEnergyLimit()) {
200 SetElementSelectors(masterModel->GetElementSelectors());
201 }
202 }
203
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
205
206 G4double
207 G4RiGeMuPairProductionModel::ComputeDEDXPerVolume(const G4Material* material,
208 const G4ParticleDefinition*,
209 G4double kineticEnergy,
210 G4double cutEnergy)
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 }
230
231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
232
233 G4double G4RiGeMuPairProductionModel::ComputMuPairLoss(G4double Z, G4double tkin,
234 G4double cutEnergy,
235 G4double tmax)
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-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 }
262
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
264
265 G4double
266 G4RiGeMuPairProductionModel::ComputeMicroscopicCrossSection(G4double tkin,
267 G4double Z,
268 G4double cutEnergy)
269 {
270 G4double cross = 0.;
271 G4double tmax = MaxSecondaryEnergyForElement(tkin, Z);
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 }
294
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
296
297 G4double G4RiGeMuPairProductionModel::ComputeDMicroscopicCrossSection(
298 G4double tkin,
299 G4double Z,
300 G4double pairEnergy)
301 // Calculates the differential (D) microscopic cross section
302 // using the cross section formula of R.P. Kokoulin (18/01/98)
303 // Code modified by R.P. Kokoulin, V.N. Ivanchenko (27/01/04)
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 }
430
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
432
433 G4double
434 G4RiGeMuPairProductionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
435 G4double kineticEnergy,
436 G4double Z, G4double,
437 G4double cutEnergy,
438 G4double maxEnergy)
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 }
454
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456
457 void G4RiGeMuPairProductionModel::MakeSamplingTables()
458 {
459 G4double factore = G4Exp(G4Log(emax/emin)/G4double(nbine));
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 }
523
524 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
525
526 void G4RiGeMuPairProductionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
527 const G4MaterialCutsCouple* couple,
528 const G4DynamicParticle* aDynamicParticle,
529 G4double tmin,
530 G4double tmax)
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 particleMomentum = aDynamicParticle->GetTotalMomentum();
538 G4ThreeVector particleMomentumVector = aDynamicParticle->GetMomentum();
539 G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
540
541 G4double minQ2 = 4.*eMass2;
542 G4double maxQ2 = (kinEnergy - particleMass)*(kinEnergy - particleMass);
543 G4double intervalQ2 = G4Log(maxQ2/minQ2);
544
545 // Square invariant of mass of the pair
546 G4double Q2 = minQ2*G4Exp(intervalQ2*randNumbs[4]);
547
548 G4double mingEnergy = std::sqrt(Q2);
549 G4double maxgEnergy = kinEnergy - particleMass;
550 G4double intervalgEnergy = maxgEnergy - mingEnergy;
551
552 // Energy of virtual gamma
553 G4double gEnergy = mingEnergy + intervalgEnergy*randNumbs[5];
554
555 // Momentum module of the virtual gamma
556 G4double gMomentum = std::sqrt(gEnergy*gEnergy - Q2);
557
558 // Energy and momentum module of the outgoing parent particle
559 G4double particleFinalEnergy = kinEnergy - gEnergy;
560 G4double particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
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 << "------- G4RiGeMuPairProductionModel::SampleSecondaries E(MeV)= "
569 // << kinEnergy << " "
570 // << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
571
572 // select randomly one element constituing the material
573 const G4Element* anElement = SelectRandomAtom(couple,particle,kinEnergy);
574
575 // define interval of energy transfer
576 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy,
577 anElement->GetZ());
578 G4double maxEnergy = std::min(tmax, maxPairEnergy);
579 G4double minEnergy = std::max(tmin, minPairEnergy);
580
581 if (minEnergy >= maxEnergy) { return; }
582
583 //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
584 // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
585 // << " ymin= " << ymin << " dy= " << dy << G4endl;
586
587 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
588
589 G4double coeff = G4Log(minPairEnergy/kinEnergy)/ymin;
590
591 // compute limits
592 G4double yymin = G4Log(minEnergy/kinEnergy)/coeff;
593 G4double yymax = G4Log(maxEnergy/kinEnergy)/coeff;
594
595 //G4cout << "yymin= " << yymin << " yymax= " << yymax << G4endl;
596
597 // units should not be used, bacause table was built without
598 G4double logTkin = G4Log(kinEnergy/CLHEP::MeV);
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= " << iz2 << G4endl;
620 do {
621 ++count;
622 // sampling using only one random number
623 G4double rand = rndmEngine->flat();
624
625 G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
626 if(iz1 != iz2) {
627 G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
628 G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1]);
629 G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2]);
630 //G4cout << count << ". x= " << x << " x2= " << x2
631 // << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
632 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
633 }
634 //G4cout << "x= " << x << " coeff= " << coeff << G4endl;
635 pairEnergy = kinEnergy*G4Exp(x*coeff);
636
637 // Loop checking, 30-Oct-2024, Vladimir Ivanchenko
638 } while((pairEnergy < minEnergy || pairEnergy > maxEnergy) && 50 > count);
639
640 //G4cout << "## pairEnergy(GeV)= " << pairEnergy/GeV
641 // << " Etot(GeV)= " << totalEnergy/GeV << G4endl;
642 rndmEngine->flatArray(9, randNumbs);
643 G4double phi3 = CLHEP::twopi*randNumbs[0];
644 fAngularGenerator->PhiRotation(partDirection, phi3);
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*particleMomentum);
653 G4double tginterval = G4Log((A1 + B1)/(A1 - B1))/B1;
654
655 G4double costg = (-A1 + (A1 - B1)*G4Exp(B1*tginterval*randNumbs[1]))/B1;
656 G4double sintg = std::sqrt((1.0 - costg)*(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, costg);
663 G4LorentzVector gFourMomentum(gEnergy, dirGamma*gMomentum);
664
665 G4double Ap = particleMomentum*particleMomentum +
666 particleFinalMomentum*particleFinalMomentum + gMomentum*gMomentum;
667 G4double A = Ap - 2.*particleMomentum*gMomentum*costg;
668 G4double B = 2.*particleMomentum*gMomentum*sintg*cospg;
669 G4double C = 2.*particleFinalMomentum*gMomentum*costg -
670 2.*particleMomentum*particleFinalMomentum;
671 G4double absB = std::abs(B);
672 G4double t1interval = (1./(A + C + absB*mint3) - 1./(A + C + absB*maxt3))/absB;
673 G4double t1 = (-(A + C) + 1./(1./(A + C + absB*mint3) - absB*t1interval*randNumbs[0]))/absB;
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(partDirection, Phi);
681 kinEnergy = particleFinalEnergy;
682
683 G4double cost5 = -1. + 2.*randNumbs[6];
684 G4double phi5 = CLHEP::twopi*randNumbs[8];
685
686 G4LorentzVector eFourMomentumMQ = fAngularGenerator->eDP2(Q2, eMass2, eMass2, cost5, phi5);
687 G4LorentzVector pFourMomentumMQ = fAngularGenerator->pDP2(eMass2, eFourMomentumMQ);
688
689 G4LorentzVector eFourMomentum = eFourMomentumMQ.boost(gFourMomentum.boostVector());
690 G4LorentzVector pFourMomentum = pFourMomentumMQ.boost(gFourMomentum.boostVector());
691
692 eEnergy = eFourMomentum.t();
693 pEnergy = pFourMomentum.t();
694
695 } else if (randNumbs[7] >= W[0] && randNumbs[7] < W[1]) {
696 G4double A3 = Q2 + 2.*gEnergy*particleFinalEnergy;
697 G4double B3 = -2.*gMomentum*particleFinalMomentum;
698
699 G4double tQ3interval = G4Log((A3 + B3)/(A3 - B3))/B3;
700 G4double tQMG = (-A3 + (A3 - B3)*G4Exp(B3*tQ3interval*randNumbs[0]))/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*particleMomentum +
708 particleFinalMomentum*particleFinalMomentum + gMomentum*gMomentum;
709 G4double A = Ap + 2.*particleFinalMomentum*gMomentum*tQMG;
710 G4double B = -2.*particleMomentum*gMomentum*sintQ3*cospQP;
711 G4double C = -2.*particleMomentum*gMomentum*tQMG - 2.*particleMomentum*particleFinalMomentum;
712
713 G4double absB = std::abs(B);
714 G4double t3interval = (1./(A + C + absB*mint3) - 1./(A + C + absB*maxt3))/absB;
715 G4double t3 = (-(A + C) + 1./(1./(A + C + absB*mint3) - absB*t3interval*randNumbs[0]))/absB;
716 G4double sint3 = std::sin(t3);
717 G4double cost3 = std::cos(t3);
718
719 G4double cost = -sint3*sintQ3*cospQP + cost3*tQMG;
720 G4double sint = std::sqrt((1. + cost)*(1. - cost));
721 G4double cosp = (sintQ3*cospQP*cost3 + sint3*tQMG)/sint;
722 G4double sinp = sintQ3*sinpQP/sint;
723
724 G4ThreeVector dirGamma;
725 dirGamma.set(sint*cosp, sint*sinp, cost);
726 G4LorentzVector gFourMomentum(gEnergy, dirGamma*gMomentum);
727
728 // Ingoing parent particle change
729 G4double Phi = CLHEP::twopi*randNumbs[3];
730 partDirection.set(sint3, 0., cost3);
731 fAngularGenerator->PhiRotation(partDirection, Phi);
732 kinEnergy = particleFinalEnergy;
733
734 G4double cost5 = -1. + 2.*randNumbs[6];
735 G4double phi5 = CLHEP::twopi*randNumbs[8];
736
737 G4LorentzVector eFourMomentumMQ = fAngularGenerator->eDP2(Q2, eMass2, eMass2, cost5, phi5);
738 G4LorentzVector pFourMomentumMQ = fAngularGenerator->pDP2(eMass2, eFourMomentumMQ);
739
740 G4LorentzVector eFourMomentum = eFourMomentumMQ.boost(gFourMomentum.boostVector());
741 G4LorentzVector pFourMomentum = pFourMomentumMQ.boost(gFourMomentum.boostVector());
742
743 eEnergy = eFourMomentum.t();
744 pEnergy = pFourMomentum.t();
745
746 } else if (randNumbs[7] >= W[1] && randNumbs[7] < W[2]) {
747 G4double phi5 = CLHEP::twopi*randNumbs[1];
748 G4double phi6 = CLHEP::twopi*randNumbs[2];
749 G4double muEnergyInterval = kinEnergy - 2.*eMass - particleMass;
750 particleFinalEnergy = particleMass + muEnergyInterval*randNumbs[3];
751 particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
752 particleMass*particleMass);
753
754 G4double mineEnergy = eMass;
755 G4double maxeEnergy = kinEnergy - particleFinalEnergy - eMass;
756 G4double eEnergyinterval = maxeEnergy - mineEnergy;
757 eEnergy = mineEnergy + eEnergyinterval*randNumbs[4];
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*eEnergy - eMass*eMass);
767 pEnergy = kinEnergy - particleFinalEnergy - eEnergy;
768 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
769
770 G4double A3 = -2.*particleMass*particleMass + 2.*kinEnergy*particleFinalEnergy;
771 G4double B3 = -2.*particleMomentum*particleFinalMomentum;
772 G4double cost3interval = G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
773 G4double expanCost3r6 = G4Exp(B3*cost3interval*randNumbs[5]);
774 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
775 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
776
777 partDirection.set(sint3, 0., cost3);
778
779 G4ThreeVector muFinalMomentumVector;
780 muFinalMomentumVector.set(particleFinalMomentum*sint3, 0., particleFinalMomentum*cost3);
781
782 G4LorentzVector muFourMomentum(particleMomentum, particleMomentumVector);
783 G4LorentzVector muFinalFourMomentum(particleFinalEnergy, muFinalMomentumVector);
784 G4LorentzVector auxVec1 = muFourMomentum - muFinalFourMomentum;
785 G4double A5 = auxVec1.mag2() - 2.*eEnergy*(kinEnergy - particleFinalEnergy) +
786 2.*particleMomentumVector[2]*eMomentum - 2.*particleFinalMomentum*eMomentum*cost3;
787 G4double B5 = -2.*particleFinalMomentum*eMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5);
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*maxt5)/(absA5 + absB5*mint5))/absB5;
793 G4double argexp = absB5*t5interval*randNumbs[6] + G4Log(absA5 + absB5*mint5);
794 G4double t5 = -absA5/absB5 + G4Exp(argexp)/absB5;
795 G4double sint5 = std::sin(t5);
796 G4double cost5 = std::cos(t5);
797
798 eDirection.set(sint5*cosp5, sint5*sinp5, cost5);
799 G4ThreeVector eMomentumVector = eMomentum*eDirection;
800
801 G4ThreeVector auxVec2 = particleMomentumVector - muFinalMomentumVector - eMomentumVector;
802 G4double p1mp3mp52 = auxVec2.dot(auxVec2);
803 G4double Bp = particleFinalMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6) +
804 eMomentum*(sint5*cosp5*cosp6 + sint5*sinp5*sinp6);
805 G4double Cp = -particleMomentum + particleFinalMomentum*cost3 + eMomentum*cost5;
806 G4double A6 = p1mp3mp52 + pMomentum*pMomentum;
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*mint6) - 1./(absA6C6 + absB6*maxt6))/absB6;
814 G4double t6 = (-absA6C6 + 1./(1./(absA6C6 + absB6*mint6) - absB6*t6interval*randNumbs[8]))/absB6;
815 G4double sint6 = std::sin(t6);
816 G4double cost6 = std::cos(t6);
817
818 pDirection.set(sint6*cosp6, sint6*sinp6, cost6);
819
820 } else {
821 G4double phi6 = CLHEP::twopi*randNumbs[1];
822 G4double phi5 = CLHEP::twopi*randNumbs[2];
823 G4double muFinalEnergyinterval = kinEnergy - 2.*eMass - particleMass;
824 particleFinalEnergy = particleMass + muFinalEnergyinterval*randNumbs[3];
825 particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
826 particleMass*particleMass);
827
828 G4double maxpEnergy = kinEnergy - particleFinalEnergy - eMass;
829 G4double pEnergyinterval = maxpEnergy - eMass;
830 pEnergy = eMass + pEnergyinterval*randNumbs[4];
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*pEnergy - eMass*eMass);
840 eEnergy = kinEnergy - particleFinalEnergy - pEnergy;
841 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
842
843 G4double A3 = -2.*particleMass*particleMass + 2.*kinEnergy*particleFinalEnergy;
844 G4double B3 = -2.*particleMomentum*particleFinalMomentum;
845 G4double cost3interval = G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
846 G4double expanCost3r6 = G4Exp(B3*cost3interval*randNumbs[5]);
847 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
848 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
849
850 partDirection.set(sint3*cosp3, sint3*sinp3, cost3);
851
852 G4ThreeVector muFinalMomentumVector;
853 muFinalMomentumVector.set(particleFinalMomentum*sint3*cosp3,
854 particleFinalMomentum*sint3*sinp3,
855 particleFinalMomentum*cost3);
856
857 G4LorentzVector muFourMomentum(particleMomentum, particleMomentumVector);
858 G4LorentzVector muFinalFourMomentum(particleFinalEnergy, muFinalMomentumVector);
859 G4LorentzVector auxVec1 = muFourMomentum - muFinalFourMomentum;
860 G4double A6 = auxVec1.mag2() -
861 2.*pEnergy*(kinEnergy - particleFinalEnergy) + 2.*particleMomentumVector[2]*pMomentum -
862 2.*particleFinalMomentum*pMomentum*cost3;
863 G4double B6 = -2.*particleFinalMomentum*pMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6);
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*maxt6)/(absA6 + absB6*mint6))/absB6;
869 G4double argexp = absB6*t6interval*randNumbs[6] + G4Log(absA6 + absB6*mint6);
870 G4double t6 = -absA6/absB6 + G4Exp(argexp)/absB6;
871 G4double sint6 = std::sin(t6);
872 G4double cost6 = std::cos(t6);
873
874 pDirection.set(sint6*cosp6, sint6*sinp6, cost6);
875 G4ThreeVector pMomentumVector = pMomentum*pDirection;
876
877 G4ThreeVector auxVec2 = particleMomentumVector - muFinalMomentumVector - pMomentumVector;
878 G4double p1mp3mp62 = auxVec2.dot(auxVec2);
879 G4double Bp = particleFinalMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5) +
880 pMomentum*(sint6*cosp6*cosp5 + sint6*sinp6*sinp5);
881 G4double Cp = -particleMomentum + particleFinalMomentum*cost3 + pMomentum*cost6;
882 G4double A5 = p1mp3mp62 + eMomentum*eMomentum;
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*mint5) - 1./(absA5C5 + absB5*maxt5))/absB5;
890 G4double t5 = (-absA5C5 + 1./(1./(absA5C5 + absB5*mint5) - absB5*t5interval*randNumbs[8]))/absB5;
891 G4double sint5 = std::sin(t5);
892 G4double cost5 = std::cos(t5);
893
894 eDirection.set(sint5*cosp5, sint5*sinp5, cost5);
895 }
896
897 fAngularGenerator->Sample5DPairDirections(aDynamicParticle, eDirection, pDirection,
898 gEnergy, Q2, gMomentum,
899 particleFinalMomentum,
900 particleFinalEnergy,
901 randNumbs, W);
902
903 // create G4DynamicParticle object for e+e-
904 auto aParticle1 = new G4DynamicParticle(theElectron, eDirection, eEnergy);
905 auto aParticle2 = new G4DynamicParticle(thePositron, pDirection, pEnergy);
906
907 // Fill output vector
908 vdp->push_back(aParticle1);
909 vdp->push_back(aParticle2);
910
911 // if energy transfer is higher than threshold (very high by default)
912 // then stop tracking the primary particle and create a new secondary
913 if (pairEnergy > SecondaryThreshold()) {
914 fParticleChange->ProposeTrackStatus(fStopAndKill);
915 fParticleChange->SetProposedKineticEnergy(0.0);
916 auto newdp = new G4DynamicParticle(particle, muF);
917 vdp->push_back(newdp);
918 } else { // continue tracking the primary e-/e+ otherwise
919 fParticleChange->SetProposedMomentumDirection(muF.vect().unit());
920 G4double ekin = std::max(muF.e() - particleMass, 0.0);
921 fParticleChange->SetProposedKineticEnergy(ekin);
922 }
923 //G4cout << "-- G4RiGeMuPairProductionModel::SampleSecondaries done" << G4endl;
924 }
925
926 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
927
928 G4double
929 G4RiGeMuPairProductionModel::FindScaledEnergy(G4int iz, G4double rand,
930 G4double logTkin,
931 G4double yymin, G4double yymax)
932 {
933 G4double res = yymin;
934 G4Physics2DVector* pv = fElementData->GetElement2DData(iz);
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], logTkin); }
940 else { res = pv->FindLinearX((pmin + rand*(pmax - pmin))/p0, logTkin); }
941 } else {
942 DataCorrupted(ZDATPAIR[iz], logTkin);
943 }
944 return res;
945 }
946
947 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
948
949 void G4RiGeMuPairProductionModel::DataCorrupted(G4int Z, G4double logTkin) const
950 {
951 G4ExceptionDescription ed;
952 ed << "G4ElementData is not properly initialized Z= " << Z
953 << " Ekin(MeV)= " << G4Exp(logTkin)
954 << " IsMasterThread= " << IsMaster()
955 << " Model " << GetName();
956 G4Exception("G4RiGeMuPairProductionModel::()", "em0033", FatalException, ed, "");
957 }
958
959 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
960
961 void G4RiGeMuPairProductionModel::StoreTables() const
962 {
963 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
964 G4int Z = ZDATPAIR[iz];
965 G4Physics2DVector* pv = fElementData->GetElement2DData(Z);
966 if(nullptr == pv) {
967 DataCorrupted(Z, 1.0);
968 return;
969 }
970 std::ostringstream ss;
971 ss << "mupair/" << particle->GetParticleName() << Z << ".dat";
972 std::ofstream outfile(ss.str());
973 pv->Store(outfile);
974 }
975 }
976
977 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
978
979 G4bool G4RiGeMuPairProductionModel::RetrieveTables()
980 {
981 for (G4int iz=0; iz<NZDATPAIR; ++iz) {
982 G4double Z = ZDATPAIR[iz];
983 G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
984 std::ostringstream ss;
985 ss << G4EmParameters::Instance()->GetDirLEDATA() << "/mupair/"
986 << particle->GetParticleName() << Z << ".dat";
987 std::ifstream infile(ss.str(), std::ios::in);
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........oooOO0OOooo........oooOO0OOooo......
998