| 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