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