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: G4BetheHeitler5DModel.cc 33 // 34 // Authors: 35 // Igor Semeniouk and Denis Bernard, 36 // LLR, Ecole polytechnique & CNRS/IN2P3, 9112 37 // 38 // Acknowledgement of the support of the Frenc 39 // (ANR-13-BS05-0002). 40 // 41 // Reference: Nucl. Instrum. Meth. A 899 (2018 42 // Nucl. Instrum. Meth., A 936 (201 43 // 44 // Class Description: 45 // 46 // Generates the conversion of a high-energy p 47 // atomic electron (triplet) or nucleus (nucle 48 // Samples the five-dimensional (5D) different 49 // . Non polarized conversion: 50 // H.A. Bethe, W. Heitler, Proc. R. Soc. Lon 51 // . Polarized conversion: 52 // T. H. Berlin and L. Madansky, Phys. Rev. 53 // M. M. May, Phys. Rev. 84 (1951) 265, 54 // J. M. Jauch and F. Rohrlich, The theory o 55 // 56 // All the above expressions are named "Bethe- 57 // 58 // Bethe & Heitler, put in Feynman diagram par 59 // the first order Born development, which is 60 // and for high-energy triplet conversion. 61 // 62 // Only the linear polarisation of the incomin 63 // The circular polarisation of the incoming p 64 // is transfered to the final leptons. 65 // 66 // In case conversion takes place in the field 67 // Bethe-Heitler expression is used. 68 // 69 // In case the nucleus or the electron are par 70 // by the other electrons of the atom is descr 71 // . nuclear: N.F. Mott, H.S.W. Massey, The Th 72 // . triplet: J.A. Wheeler and W.E. Lamb, Phys 73 // 74 // The nuclear form factor that affects the pr 75 // 76 // In principle the code is valid from thresho 77 // 4 * m_e c^2 for triplet, up to infinity, wh 78 // cross section at small q2 and, at high-ener 79 // some point that depends on machine precisio 80 // 81 // Very-high-energy (above a few tens of TeV) 82 // cross-section are not considered. 83 // 84 // The 5D differential cross section is sample 85 // angle approximation(s). 86 // The generation is strictly energy-momentum 87 // are taken into account, that is, including 88 // (In contrast with the BH expressions taken 89 // taken to be EMinus = GammaEnergy - EPlus) 90 // 91 // Tests include the examination of 1D distrib 92 // 93 // Total cross sections are not computed (we i 94 // We just convert a photon on a target when a 95 // 96 // Pure nuclear, pure triplet and 1/Z triplet/ 97 // 98 // ------------------------------------------- 99 100 #include "G4BetheHeitler5DModel.hh" 101 #include "G4EmParameters.hh" 102 103 #include "G4PhysicalConstants.hh" 104 #include "G4SystemOfUnits.hh" 105 #include "G4Electron.hh" 106 #include "G4Positron.hh" 107 #include "G4Gamma.hh" 108 #include "G4IonTable.hh" 109 #include "G4NucleiProperties.hh" 110 111 #include "Randomize.hh" 112 #include "G4ParticleChangeForGamma.hh" 113 #include "G4Pow.hh" 114 #include "G4Log.hh" 115 #include "G4Exp.hh" 116 117 #include "G4LorentzVector.hh" 118 #include "G4ThreeVector.hh" 119 #include "G4RotationMatrix.hh" 120 121 #include <cassert> 122 123 const G4int kEPair = 0; 124 const G4int kMuPair = 1; 125 126 127 //....oooOO0OOooo........oooOO0OOooo........oo 128 129 G4BetheHeitler5DModel::G4BetheHeitler5DModel(c 130 c 131 : G4PairProductionRelModel(pd, nam), 132 fLepton1(G4Electron::Definition()),fLepton 133 fTheMuPlus(nullptr),fTheMuMinus(nullptr), 134 fVerbose(1), 135 fConversionType(0), 136 fConvMode(kEPair), 137 iraw(false) 138 { 139 theIonTable = G4IonTable::GetIonTable(); 140 //Q: Do we need this on Model 141 SetLowEnergyLimit(2*fTheElectron->GetPDGMass 142 } 143 144 //....oooOO0OOooo........oooOO0OOooo........oo 145 146 G4BetheHeitler5DModel::~G4BetheHeitler5DModel( 147 148 //....oooOO0OOooo........oooOO0OOooo........oo 149 150 void G4BetheHeitler5DModel::Initialise(const G 151 const G4DataVector& vec) 152 { 153 G4PairProductionRelModel::Initialise(part, v 154 155 G4EmParameters* theManager = G4EmParameters: 156 // place to initialise model parameters 157 // Verbosity levels: ( Can redefine as neede 158 // 0 = nothing 159 // > 2 print results 160 // > 3 print rejection warning from transfor 161 // > 4 print photon direction & polarisation 162 fVerbose = theManager->Verbose(); 163 fConversionType = theManager->GetConversionT 164 //////////////////////////////////////////// 165 // iraw : 166 // true : isolated electron or nucleus 167 // false : inside atom -> screening for 168 iraw = theManager->OnIsolated(); 169 // G4cout << "BH5DModel::Initialise verbose 170 // << " isolated " << iraw << " ctype "<< 171 172 //Q: Do we need this on Model 173 // The Leptons defined via SetLeptonPair(..) 174 SetLowEnergyLimit(2*CLHEP::electron_mass_c2) 175 } 176 177 //....oooOO0OOooo........oooOO0OOooo........oo 178 179 void G4BetheHeitler5DModel::SetLeptonPair(cons 180 const G4ParticleDefinition* p2) 181 { 182 G4int pdg1 = p1->GetPDGEncoding(); 183 G4int pdg2 = p2->GetPDGEncoding(); 184 G4int pdg = std::abs(pdg1); 185 if ( pdg1 != -pdg2 || (pdg != 11 && pdg != 1 186 G4ExceptionDescription ed; 187 ed << " Wrong pair of leptons: " << p1->Ge 188 << " and " << p1->GetParticleName(); 189 G4Exception("G4BetheHeitler5DModel::SetLep 190 FatalErrorInArgument, ed, ""); 191 } else { 192 if ( pdg == 11 ) { 193 SetConversionMode(kEPair); 194 if( pdg1 == 11 ) { 195 fLepton1 = p1; 196 fLepton2 = p2; 197 } else { 198 fLepton1 = p2; 199 fLepton2 = p1; 200 } 201 if (fVerbose > 0) 202 G4cout << "G4BetheHeitler5DModel::SetLeptonP 203 << G4endl; 204 } else { 205 SetConversionMode(kMuPair); 206 if( pdg1 == 13 ) { 207 fLepton1 = p1; 208 fLepton2 = p2; 209 } else { 210 fLepton1 = p2; 211 fLepton2 = p1; 212 } 213 fTheMuPlus = fLepton2; 214 fTheMuMinus= fLepton1; 215 if (fVerbose > 0) 216 G4cout << "G4BetheHeitler5DModel::SetLeptonP 217 << G4endl; 218 } 219 } 220 } 221 222 //....oooOO0OOooo........oooOO0OOooo........oo 223 224 G4double G4BetheHeitler5DModel::MaxDiffCrossSe 225 226 227 228 { 229 const G4double Q = e/par[9]; 230 return par[0] * G4Exp((par[2]+loge*par[4])*l 231 / (par[1]+ G4Exp(par[3]*loge)+G4Exp(p 232 * (1+par[7]*G4Exp(par[8]*G4Log(Z))*Q/ 233 } 234 235 //....oooOO0OOooo........oooOO0OOooo........oo 236 237 void 238 G4BetheHeitler5DModel::SampleSecondaries(std:: 239 const 240 const 241 G4dou 242 { 243 // MeV 244 static const G4double ElectronMass = CLHEP 245 246 const G4double LeptonMass = fLepton1->GetPDG 247 const G4double LeptonMass2 = LeptonMass*Lep 248 249 static const G4double alpha0 = CLHEP 250 // mm 251 static const G4double r0 = CLHEP 252 // mbarn 253 static const G4double r02 = r0*r0 254 static const G4double twoPi = CLHEP 255 static const G4double factor = alpha 256 // static const G4double factor1 = p 257 static const G4double factor1 = 2.661 258 // 259 G4double PairInvMassMin = 2.*LeptonMass; 260 G4double TrThreshold = 2.0 * ( (LeptonMass2 261 262 // 263 static const G4double nu[2][10] = { 264 //electron 265 { 0.0227436, 0.0582046, 3.0322675, 2.82750 266 1.1212766, 1.8989468, 68.3492750, 0.0211 267 //muon 268 {0.67810E-06, 0.86037E+05, 2.0008395, 1.67 269 1.4222, 0.0, 263230.0, 0.0521, 51.1338} 270 }; 271 static const G4double tr[2][10] = { 272 //electron 273 { 0.0332350, 4.3942537, 2.8515925, 2.6351 274 1.5737305, 1.8104647, 20.6434021, -0.027 275 //muon 276 {0.10382E-03, 0.14408E+17, 4.1368679, 3.26 277 0.0000, 0.0, 0.0, 0.0000, 1.0000} 278 }; 279 // 280 static const G4double para[2][3][2] = { 281 //electron 282 { {11., -16.},{-1.17, -2.95},{-2., -0.5} } 283 //muon 284 { {17.5, 1.},{-1.17, -2.95},{2., 6.} } 285 }; 286 // 287 static const G4double correctionIndex = 1.4; 288 // 289 const G4double GammaEnergy = aDynamicGamma- 290 // Protection, Will not be true tot cross se 291 if ( GammaEnergy <= PairInvMassMin) { return 292 293 const G4double GammaEnergy2 = GammaEnergy*Ga 294 295 //////////////////////////////////////////// 296 const G4ParticleMomentum GammaDirection = 297 aDynamicGamma->GetMomentumDirection(); 298 G4ThreeVector GammaPolarization = aDynamicGa 299 300 // The protection polarization perpendicular 301 // as it done in G4LivermorePolarizedGammaCo 302 // assuming Direction is unitary vector 303 // (projection to plane) p_proj = p - (p o 304 if ( GammaPolarization.howOrthogonal(GammaDi 305 GammaPolarization -= GammaPolarization.dot 306 } 307 // End of Protection 308 // 309 const G4double GammaPolarizationMag = GammaP 310 311 //////////////////////////////////////////// 312 // target element 313 // select randomly one element constituting 314 const G4Element* anElement = SelectTargetAt 315 aDyna 316 // Atomic number 317 const G4int Z = anElement->GetZasInt() 318 const G4int A = SelectIsotopeNumber(an 319 const G4double iZ13 = 1./anElement->GetIonis 320 const G4double targetMass = G4NucleiProperti 321 322 const G4double NuThreshold = 2.0 * ( (Lept 323 // No conversion possible below nuclear thre 324 if ( GammaEnergy <= NuThreshold) { return; } 325 326 CLHEP::HepRandomEngine* rndmEngine = G4Rando 327 328 // itriplet : true -- triplet, false -- nucl 329 G4bool itriplet = false; 330 if (fConversionType == 1) { 331 itriplet = false; 332 } else if (fConversionType == 2) { 333 itriplet = true; 334 if ( GammaEnergy <= TrThreshold ) return; 335 } else if ( GammaEnergy > TrThreshold ) { 336 // choose triplet or nuclear from a triple 337 // total cross section ratio. 338 // approximate at low energies ! 339 if(rndmEngine->flat()*(Z+1) < 1.) { 340 itriplet = true; 341 } 342 } 343 344 // 345 const G4double RecoilMass = itriplet ? Elec 346 const G4double RecoilMass2 = RecoilMass*Reco 347 const G4double sCMS = 2.*RecoilMass*G 348 const G4double sCMSPlusRM2 = sCMS + RecoilMa 349 const G4double sqrts = std::sqrt(sCMS) 350 const G4double isqrts2 = 1./(2.*sqrts); 351 // 352 const G4double PairInvMassMax = sqrts-Reco 353 const G4double PairInvMassRange = PairInvMas 354 const G4double lnPairInvMassRange = G4Log(Pa 355 356 // initial state. Defines z axis of "0" fram 357 // Since CMS(0., 0., GammaEnergy, GammaEnerg 358 const G4double betaCMS = G4LorentzVector(0.0 359 360 // maximum value of pdf 361 const G4double EffectiveZ = iraw ? 0.5 : Z; 362 const G4double Threshold = itriplet ? TrThr 363 const G4double AvailableEnergy = GammaEne 364 const G4double LogAvailableEnergy = G4Log(Av 365 // 366 const G4double MaxDiffCross = itriplet 367 ? MaxDiffCrossSection(tr[fConvMode], 368 EffectiveZ, AvailableEnergy, LogAvaila 369 : MaxDiffCrossSection(nu[fConvMode], 370 EffectiveZ, AvailableEnergy, LogAva 371 // 372 // 50% safety marging factor 373 const G4double ymax = 1.5 * MaxDiffCross; 374 // x1 bounds 375 const G4double xu1 = (LogAvailableEnergy > 376 ? para[fConvMode][0][0] + 377 para[fConvMode][1][0]*LogAvailableEner 378 : para[fConvMode][0][0] 379 para[fConvMode][2][0]*para[fConvMode][ 380 const G4double xl1 = (LogAvailableEnergy > 381 ? para[fConvMode][0][1] 382 para[fConvMode][1][1]*LogAvailableEner 383 : para[fConvMode][0][1] 384 para[fConvMode][2][1]*para[fConvMode][ 385 // 386 G4LorentzVector Recoil; 387 G4LorentzVector LeptonPlus; 388 G4LorentzVector LeptonMinus; 389 G4double pdf = 0.; 390 391 G4double rndmv6[6] = {0.0}; 392 const G4double corrFac = 1.0/(correctionInde 393 const G4double expLowLim = -20.; 394 const G4double logLowLim = G4Exp(expLowLim/c 395 G4double z0, z1, z2, x0, x1; 396 G4double betheheitler, sinTheta, cosTheta, d 397 // START Sampling 398 do { 399 400 rndmEngine->flatArray(6, rndmv6); 401 402 ////////////////////////////////////////// 403 // pdf pow(x,c) with c = 1.4 404 // integral y = pow(x,(c+1))/(c+1) @ x = 1 405 // invCdf exp( log(y /* *( c + 1.0 )/ (c + 406 ////////////////////////////////////////// 407 408 z0 = (rndmv6[0] > logLowLim) ? G4Log(rndmv 409 G4double X1 = (z0 > expLowLim) ? G4Exp(z0) 410 z1 = xl1 + (xu1 - xl1)*rndmv6[1]; 411 if (z1 > expLowLim) { 412 x0 = G4Exp(z1); 413 dum0 = 1.0/(1.0 + x0); 414 x1 = dum0*x0; 415 cosTheta = -1.0 + 2.0*x1; 416 sinTheta = 2*std::sqrt(x1*(1.0 - x1)); 417 } else { 418 x0 = 0.0; 419 dum0 = 1.0; 420 cosTheta = -1.0; 421 sinTheta = 0.0; 422 } 423 424 z2 = X1*X1*lnPairInvMassRange; 425 const G4double PairInvMass = PairInvMassMi 426 427 // cos and sin theta-lepton 428 const G4double cosThetaLept = std::cos(pi* 429 // sin(ThetaLept) is always in [0,+1] if T 430 const G4double sinThetaLept = std::sqrt((1 431 // cos and sin phi-lepton 432 const G4double cosPhiLept = std::cos(two 433 // sin(PhiLept) is in [-1,0] if PhiLept in 434 // is in [0,+1] if PhiLept in 435 const G4double sinPhiLept = std::copysig 436 // cos and sin phi 437 const G4double cosPhi = std::cos(two 438 const G4double sinPhi = std::copysig 439 440 ////////////////////////////////////////// 441 // frames: 442 // 3 : the laboratory Lorentz frame, Geant 443 // 0 : the laboratory Lorentz frame, axes 444 // 1 : the center-of-mass Lorentz frame 445 // 2 : the pair Lorentz frame 446 ////////////////////////////////////////// 447 448 // in the center-of-mass frame 449 450 const G4double RecEnergyCMS = (sCMSPlusRM 451 const G4double LeptonEnergy2 = PairInvMass 452 453 // New way of calucaltion thePRecoil to av 454 G4double abp = std::max((2.0*GammaEnergy*R 455 PairInvMass*PairInvMass + 2.0*PairI 456 (2.0*GammaEnergy*R 457 PairInvMass*PairInvMass - 2.0*PairI 458 459 G4double thePRecoil = std::sqrt(abp) * isq 460 461 // back to the center-of-mass frame 462 Recoil.set( thePRecoil*sinTheta*cosPhi, 463 thePRecoil*sinTheta*sinPhi, 464 thePRecoil*cosTheta, 465 RecEnergyCMS); 466 467 // in the pair frame 468 const G4double thePLepton = std::sqrt( 469 * 470 471 LeptonPlus.set(thePLepton*sinThetaLept*cos 472 thePLepton*sinThetaLept*sinPhiLept, 473 thePLepton*cosThetaLept, 474 LeptonEnergy2); 475 476 LeptonMinus.set(-LeptonPlus.x(), 477 -LeptonPlus.y(), 478 -LeptonPlus.z(), 479 LeptonEnergy2); 480 481 482 // Normalisation of final state phase spac 483 // Section 47 of Particle Data Group, Chin 484 // const G4double Norme = Recoil1.vect( 485 const G4double Norme = Recoil.vect().mag() 486 487 // e+, e- to CMS frame from pair frame 488 489 // boost vector from Pair to CMS 490 const G4ThreeVector pair2cms = 491 G4LorentzVector( -Recoil.x(), -Recoil.y(), 492 sqrts-RecEnergyCMS).boostVector(); 493 494 LeptonPlus.boost(pair2cms); 495 LeptonMinus.boost(pair2cms); 496 497 // back to the laboratory frame (make use 498 499 Recoil.boostZ(betaCMS); 500 LeptonPlus.boostZ(betaCMS); 501 LeptonMinus.boostZ(betaCMS); 502 503 // Jacobian factors 504 const G4double Jacob0 = x0*dum0*dum0; 505 const G4double Jacob1 = 2.*X1*lnPairInvMas 506 const G4double Jacob2 = std::abs(sinThetaL 507 508 // there is no probability to have a lepto 509 // X and Y components of momentum may be z 510 const G4double EPlus = LeptonPlus.t(); 511 const G4double PPlus = LeptonPlus.vect().m 512 const G4double pPX = LeptonPlus.x(); 513 const G4double pPY = LeptonPlus.y(); 514 const G4double pPZ = LeptonPlus.z(); 515 G4double sinPhiPlus = 1.0; 516 G4double cosPhiPlus = 0.0; 517 G4double sinThetaPlus = 0.0; 518 G4double cosThetaPlus = pPZ/PPlus; 519 if (cosThetaPlus < 1.0 && cosThetaPlus > - 520 sinThetaPlus = std::sqrt((1.0 - cosTheta 521 sinPhiPlus = pPY/(PPlus*sinThetaPlus); 522 cosPhiPlus = pPX/(PPlus*sinThetaPlus); 523 } 524 525 // denominators: 526 // the two cancelling leading terms for fo 527 const G4double elMassCTP = LeptonMass*cosT 528 const G4double ePlusSTP = EPlus*sinThetaP 529 const G4double DPlus = (elMassCTP*elMa 530 /(EPlus + PPlus* 531 532 // there is no probability to have a lepto 533 // X and Y components of momentum may be z 534 const G4double EMinus = LeptonMinus.t(); 535 const G4double PMinus = LeptonMinus.vect() 536 const G4double ePX = LeptonMinus.x(); 537 const G4double ePY = LeptonMinus.y(); 538 const G4double ePZ = LeptonMinus.z(); 539 G4double sinPhiMinus = 0.0; 540 G4double cosPhiMinus = 1.0; 541 G4double sinThetaMinus = 0.0; 542 G4double cosThetaMinus = ePZ/PMinus; 543 if (cosThetaMinus < 1.0 && cosThetaMinus > 544 sinThetaMinus = std::sqrt((1.0 - cosThet 545 sinPhiMinus = ePY/(PMinus*sinThetaMinus) 546 cosPhiMinus = ePX/(PMinus*sinThetaMinus) 547 } 548 549 const G4double elMassCTM = LeptonMass*cosT 550 const G4double eMinSTM = EMinus*sinTheta 551 const G4double DMinus = (elMassCTM*elMa 552 /(EMinus + PMinu 553 554 // cos(phiMinus-PhiPlus) 555 const G4double cosdPhi = cosPhiPlus*cosPhi 556 const G4double PRec = Recoil.vect().mag 557 const G4double q2 = PRec*PRec; 558 559 const G4double BigPhi = -LeptonMass2 / (G 560 561 G4double FormFactor = 1.; 562 if (!iraw) { 563 if (itriplet) { 564 const G4double qun = factor1*iZ13*iZ13; 565 const G4double nun = qun * PRec; 566 if (nun < 1.) { 567 FormFactor = (nun < 0.01) ? (13.8-5 568 : std::sq 569 } // else FormFactor = 1 by default 570 } else { 571 const G4double dum3 = 217.*PRec*iZ13; 572 const G4double AFF = 1./(1. + dum3*dum3); 573 FormFactor = (1.-AFF)*(1-AFF); 574 } 575 } // else FormFactor = 1 by default 576 577 if (GammaPolarizationMag==0.) { 578 const G4double pPlusSTP = PPlus*sinThe 579 const G4double pMinusSTM = PMinus*sinTh 580 const G4double pPlusSTPperDP = pPlusSTP 581 const G4double pMinusSTMperDM = pMinusST 582 const G4double dunpol = BigPhi*( 583 pPlusSTPperDP *pPlusSTPperDP 584 + pMinusSTMperDM*pMinusSTMperD 585 + 2.*pPlusSTPperDP*pMinusSTMpe 586 *(4.*EPlus*EMinus + q2 - 2 587 - 2.*GammaEnergy2*(pPlusSTP*pP 588 betheheitler = dunpol * factor; 589 } else { 590 const G4double pPlusSTP = PPlus*sinThet 591 const G4double pMinusSTM = PMinus*sinThe 592 const G4double pPlusSTPCPPperDP = pPlus 593 const G4double pMinusSTMCPMperDM = pMinu 594 const G4double caa = 2.*(EPlus*pMinusSTM 595 const G4double cbb = pMinusSTMCPMperDM-p 596 const G4double ccc = (pPlusSTP*pPlusSTP 597 +2.*pPlusSTP*pMinusS 598 const G4double dtot= 2.*BigPhi*( caa*caa 599 betheheitler = dtot * factor; 600 } 601 // 602 const G4double cross = Norme * Jacob0 * J 603 * FormFactor * Recoi 604 pdf = cross * (xu1 - xl1) / G4Exp(correcti 605 } while ( pdf < ymax * rndmv6[5] ); 606 // END of Sampling 607 608 if ( fVerbose > 2 ) { 609 G4double recul = std::sqrt(Recoil.x()*Reco 610 +Recoil.z()*Reco 611 G4cout << "BetheHeitler5DModel GammaEnergy 612 << " PDF= " << pdf << " ymax= " << ymax 613 << " recul= " << recul << G4endl; 614 } 615 616 // back to Geant4 system 617 618 if ( fVerbose > 4 ) { 619 G4cout << "BetheHeitler5DModel GammaDirect 620 G4cout << "BetheHeitler5DModel GammaPolari 621 G4cout << "BetheHeitler5DModel GammaEnergy 622 G4cout << "BetheHeitler5DModel Conv " 623 << (itriplet ? "triplet" : "nucl") << G4e 624 } 625 626 if (GammaPolarizationMag == 0.0) { 627 // set polarization axis orthohonal to dir 628 GammaPolarization = GammaDirection.orthogo 629 } else { 630 // GammaPolarization not a unit vector 631 GammaPolarization /= GammaPolarizationMag; 632 } 633 634 // The unit norm vector that is orthogonal t 635 G4ThreeVector yGrec = GammaDirection.cross(G 636 637 // rotation from gamma ref. sys. to World 638 G4RotationMatrix GtoW(GammaPolarization,yGre 639 640 Recoil.transform(GtoW); 641 LeptonPlus.transform(GtoW); 642 LeptonMinus.transform(GtoW); 643 644 if ( fVerbose > 2 ) { 645 G4cout << "BetheHeitler5DModel Recoil " << 646 << " " << Recoil.t() << " " << G4endl; 647 G4cout << "BetheHeitler5DModel LeptonPlus 648 << LeptonPlus.z() << " " << LeptonPlus.t( 649 G4cout << "BetheHeitler5DModel LeptonMinus 650 << LeptonMinus.z() << " " << LeptonMinus. 651 } 652 653 // Create secondaries 654 auto aParticle1 = new G4DynamicParticle(fLep 655 auto aParticle2 = new G4DynamicParticle(fLep 656 657 // create G4DynamicParticle object for the p 658 G4ParticleDefinition* RecoilPart; 659 if (itriplet) { 660 // triplet 661 RecoilPart = fTheElectron; 662 } else{ 663 RecoilPart = theIonTable->GetIon(Z, A, 0); 664 } 665 auto aParticle3 = new G4DynamicParticle(Reco 666 667 // Fill output vector 668 fvect->push_back(aParticle1); 669 fvect->push_back(aParticle2); 670 fvect->push_back(aParticle3); 671 672 // kill incident photon 673 fParticleChange->SetProposedKineticEnergy(0. 674 fParticleChange->ProposeTrackStatus(fStopAnd 675 } 676 677 //....oooOO0OOooo........oooOO0OOooo........oo 678