ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BetheHeitler5DModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BetheHeitler5DModel.cc
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: G4BetheHeitler5DModel.cc
33 //
34 // Authors:
35 // Igor Semeniouk and Denis Bernard,
36 // LLR, Ecole polytechnique & CNRS/IN2P3, 91128 Palaiseau, France
37 //
38 // Acknowledgement of the support of the French National Research Agency
39 // (ANR-13-BS05-0002).
40 //
41 // Reference: Nucl. Instrum. Meth. A 899 (2018) 85 (arXiv:1802.08253 [hep-ph])
42 // Nucl. Instrum. Meth., A 936 (2019) 290
43 //
44 // Class Description:
45 //
46 // Generates the conversion of a high-energy photon to an e+e- pair, either in the field of an
47 // atomic electron (triplet) or nucleus (nuclear).
48 // Samples the five-dimensional (5D) differential cross-section analytical expression:
49 // . Non polarized conversion:
50 // H.A. Bethe, W. Heitler, Proc. R. Soc. Lond. Ser. A 146 (1934) 83.
51 // . Polarized conversion:
52 // T. H. Berlin and L. Madansky, Phys. Rev. 78 (1950) 623,
53 // M. M. May, Phys. Rev. 84 (1951) 265,
54 // J. M. Jauch and F. Rohrlich, The theory of photons and electrons, 1976.
55 //
56 // All the above expressions are named "Bethe-Heitler" here.
57 //
58 // Bethe & Heitler, put in Feynman diagram parlance, compute only the two dominant diagrams of
59 // the first order Born development, which is an excellent approximation for nuclear conversion
60 // and for high-energy triplet conversion.
61 //
62 // Only the linear polarisation of the incoming photon takes part in these expressions.
63 // The circular polarisation of the incoming photon does not (take part) and no polarisation
64 // is transfered to the final leptons.
65 //
66 // In case conversion takes place in the field of an isolated nucleus or electron, the bare
67 // Bethe-Heitler expression is used.
68 //
69 // In case the nucleus or the electron are part of an atom, the screening of the target field
70 // by the other electrons of the atom is described by a simple form factor, function of q2:
71 // . nuclear: N.F. Mott, H.S.W. Massey, The Theory of Atomic Collisions, 1934.
72 // . triplet: J.A. Wheeler and W.E. Lamb, Phys. Rev. 55 (1939) 858.
73 //
74 // The nuclear form factor that affects the probability of very large-q2 events, is not considered.
75 //
76 // In principle the code is valid from threshold, that is from 2 * m_e c^2 for nuclear and from
77 // 4 * m_e c^2 for triplet, up to infinity, while in pratice the divergence of the differential
78 // cross section at small q2 and, at high-energy, at small polar angle, make it break down at
79 // some point that depends on machine precision.
80 //
81 // Very-high-energy (above a few tens of TeV) LPM suppression effects in the normalized differential
82 // cross-section are not considered.
83 //
84 // The 5D differential cross section is sampled without any high-energy nor small
85 // angle approximation(s).
86 // The generation is strictly energy-momentum conserving when all particles in the final state
87 // are taken into account, that is, including the recoiling target.
88 // (In contrast with the BH expressions taken at face values, for which the electron energy is
89 // taken to be EMinus = GammaEnergy - EPlus)
90 //
91 // Tests include the examination of 1D distributions: see TestEm15
92 //
93 // Total cross sections are not computed (we inherit from other classes).
94 // We just convert a photon on a target when asked to do so.
95 //
96 // Pure nuclear, pure triplet and 1/Z triplet/nuclear mixture can be generated.
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 "G4MuonPlus.hh"
109 #include "G4MuonMinus.hh"
110 #include "G4IonTable.hh"
111 #include "G4NucleiProperties.hh"
112 
113 #include "Randomize.hh"
114 #include "G4ParticleChangeForGamma.hh"
115 #include "G4Pow.hh"
116 #include "G4Log.hh"
117 #include "G4Exp.hh"
118 
119 #include "G4LorentzVector.hh"
120 #include "G4ThreeVector.hh"
121 #include "G4RotationMatrix.hh"
122 
123 #include <cassert>
124 
125 // // Q : Use enum G4EmProcessSubType hire ?
126 // enum G45DConversionMode
127 // {
128 // kEPair, kMuPair
129 // };
130 
131 const G4int kEPair = 0;
132 const G4int kMuPair = 1;
133 
134 
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
136 
137 G4BetheHeitler5DModel::G4BetheHeitler5DModel(const G4ParticleDefinition* pd,
138  const G4String& nam)
139  : G4PairProductionRelModel(pd, nam),fVerbose(1),fConversionType(0),
140  iraw(false),
141  fLepton1(G4Electron::Definition()),fLepton2(G4Positron::Definition()),
142  fConvMode(kEPair),
143  fTheMuPlus(G4MuonPlus::Definition()),fTheMuMinus(G4MuonMinus::Definition())
144 {
145  theIonTable = G4IonTable::GetIonTable();
146  //Q: Do we need this on Model
147  SetLowEnergyLimit(2*fTheElectron->GetPDGMass());
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 
152 G4BetheHeitler5DModel::~G4BetheHeitler5DModel()
153 {}
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156 
157 void G4BetheHeitler5DModel::Initialise(const G4ParticleDefinition* part,
158  const G4DataVector& vec)
159 {
160  G4PairProductionRelModel::Initialise(part, vec);
161 
162  G4EmParameters* theManager = G4EmParameters::Instance();
163  // place to initialise model parameters
164  // Verbosity levels: ( Can redefine as needed, but some consideration )
165  // 0 = nothing
166  // > 2 print results
167  // > 3 print rejection warning from transformation (fix bug from gammaray .. )
168  // > 4 print photon direction & polarisation
169  fVerbose = theManager->Verbose();
170  fConversionType = theManager->GetConversionType();
172  // iraw :
173  // true : isolated electron or nucleus.
174  // false : inside atom -> screening form factor
175  iraw = theManager->OnIsolated();
176  // G4cout << "BH5DModel::Initialise verbose " << fVerbose
177  // << " isolated " << iraw << " ctype "<< fConversionType << G4endl;
178 
179  //Q: Do we need this on Model
180  // The Leptons defined via SetLeptonPair(..) method
181  SetLowEnergyLimit(2*CLHEP::electron_mass_c2);
182 
183  if (fConvMode == kEPair) {
184  assert(fLepton1->GetPDGEncoding() == fTheElectron->GetPDGEncoding()) ;
185  if (fVerbose > 3)
186  G4cout << "BH5DModel::Initialise conversion to e+ e-" << G4endl;
187  }
188 
189  if (fConvMode == kMuPair) {
190  assert(fLepton1->GetPDGEncoding() == fTheMuMinus->GetPDGEncoding()) ;
191  if (fVerbose > 3)
192  G4cout << "BH5DModel::Initialise conversion to mu+ mu-" << G4endl;
193  }
194 }
195 
196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
197 
198 void G4BetheHeitler5DModel::SetLeptonPair(const G4ParticleDefinition* p1,
199  const G4ParticleDefinition* p2)
200 {
201  // Lepton1 - nagative charged particle
202  if ( p1->GetPDGEncoding() < 0 ){
203  if ( p1->GetPDGEncoding() ==
204  G4Positron::Definition()->GetPDGEncoding() ) {
205  SetConversionMode(kEPair);
206  fLepton1 = p2;
207  fLepton2 = p1;
208  // if (fVerbose)
209  G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to e+ e-"
210  << G4endl;
211  } else if ( p1->GetPDGEncoding() ==
212  G4MuonPlus::Definition()->GetPDGEncoding() ) {
213  SetConversionMode(kMuPair);
214  fLepton1 = p2;
215  fLepton2 = p1;
216  // if (fVerbose)
217  G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to mu+ mu-"
218  << G4endl;
219  } else {
220  // Exception
221  G4ExceptionDescription ed;
222  ed << "Model not applicable to particle(s) "
223  << p1->GetParticleName() << ", "
224  << p2->GetParticleName();
225  G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0002",
226  FatalException, ed);
227  }
228  } else {
229  if ( p1->GetPDGEncoding() ==
230  G4Electron::Definition()->GetPDGEncoding() ) {
231  SetConversionMode(kEPair);
232  fLepton1 = p1;
233  fLepton2 = p2;
234  // if (fVerbose)
235  G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to e+ e-"
236  << G4endl;
237  } else if ( p1->GetPDGEncoding() ==
238  G4MuonMinus::Definition()->GetPDGEncoding() ) {
239  SetConversionMode(kMuPair);
240  fLepton1 = p1;
241  fLepton2 = p2;
242  // if (fVerbose)
243  G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to mu+ mu-"
244  << G4endl;
245  } else {
246  // Exception
247  G4ExceptionDescription ed;
248  ed << "Model not applicable to particle(s) "
249  << p1->GetParticleName() << ", "
250  << p2->GetParticleName();
251  G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0002",
252  FatalException, ed);
253  }
254  }
255  if ( fLepton1->GetPDGEncoding() != fLepton2->GetAntiPDGEncoding() ) {
256  G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0007",
257  FatalErrorInArgument, "pair must be particle, antiparticle ");
258  G4cerr << "BH5DModel::SetLeptonPair BAD paricle/anti particle pair"
259  << fLepton1->GetParticleName() << ", "
260  << fLepton2->GetParticleName() << G4endl;
261  }
262 }
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265 
266 G4double G4BetheHeitler5DModel::MaxDiffCrossSection(const G4double* par,
267  G4double Z,
268  G4double e,
269  G4double loge) const
270 {
271  const G4double Q = e/par[9];
272  return par[0] * G4Exp((par[2]+loge*par[4])*loge)
273  / (par[1]+ G4Exp(par[3]*loge)+G4Exp(par[5]*loge))
274  * (1+par[7]*G4Exp(par[8]*G4Log(Z))*Q/(1+Q));
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278 
279 void
280 G4BetheHeitler5DModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
281  const G4MaterialCutsCouple* couple,
282  const G4DynamicParticle* aDynamicGamma,
283  G4double, G4double)
284 {
285  // MeV
286  static const G4double ElectronMass = CLHEP::electron_mass_c2;
287 
288  const G4double LeptonMass = fLepton1->GetPDGMass();
289  const G4double LeptonMass2 = LeptonMass*LeptonMass;
290 
291  static const G4double alpha0 = CLHEP::fine_structure_const;
292  // mm
293  static const G4double r0 = CLHEP::classic_electr_radius;
294  // mbarn
295  static const G4double r02 = r0*r0*1.e+25;
296  static const G4double twoPi = CLHEP::twopi;
297  static const G4double factor = alpha0 * r02 / (twoPi*twoPi);
298  // static const G4double factor1 = pow((6.0 * pi),(1.0/3.0))/(8.*alpha0*ElectronMass);
299  static const G4double factor1 = 2.66134007899/(8.*alpha0*ElectronMass);
300  //
301  G4double PairInvMassMin = 2.*LeptonMass;
302  G4double TrThreshold = 2.0 * ( (LeptonMass2)/ElectronMass + LeptonMass);
303 
304  //
305  static const G4double nu[2][10] = {
306  //electron
307  { 0.0227436, 0.0582046, 3.0322675, 2.8275065, -0.0034004,
308  1.1212766, 1.8989468, 68.3492750, 0.0211186, 14.4},
309  //muon
310  {0.67810E-06, 0.86037E+05, 2.0008395, 1.6739719, -0.0057279,
311  1.4222, 0.0, 263230.0, 0.0521, 51.1338}
312  };
313  static const G4double tr[2][10] = {
314  //electron
315  { 0.0332350, 4.3942537, 2.8515925, 2.6351695, -0.0031510,
316  1.5737305, 1.8104647, 20.6434021, -0.0272586, 28.9},
317  //muon
318  {0.10382E-03, 0.14408E+17, 4.1368679, 3.2662121, -0.0163091,
319  0.0000, 0.0, 0.0, 0.0000, 1.0000}
320  };
321  //
322  static const G4double para[2][3][2] = {
323  //electron
324  { {11., -16.},{-1.17, -2.95},{-2., -0.5} },
325  //muon
326  { {17.5, 1.},{-1.17, -2.95},{2., 6.} }
327  };
328  //
329  static const G4double correctionIndex = 1.4;
330  //
331  const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
332  // Protection, Will not be true tot cross section = 0
333  if ( GammaEnergy <= PairInvMassMin) { return; }
334 
335  const G4double GammaEnergy2 = GammaEnergy*GammaEnergy;
336 
338  const G4ParticleMomentum GammaDirection =
339  aDynamicGamma->GetMomentumDirection();
340  G4ThreeVector GammaPolarization = aDynamicGamma->GetPolarization();
341 
342  // The protection polarization perpendicular to the direction vector,
343  // as it done in G4LivermorePolarizedGammaConversionModel,
344  // assuming Direction is unitary vector
345  // (projection to plane) p_proj = p - (p o d)/(d o d) x d
346  if ( GammaPolarization.howOrthogonal(GammaDirection) != 0) {
347  GammaPolarization -= GammaPolarization.dot(GammaDirection) * GammaDirection;
348  }
349  // End of Protection
350  //
351  const G4double GammaPolarizationMag = GammaPolarization.mag();
352 
354  // target element
355  // select randomly one element constituting the material
356  const G4Element* anElement = SelectTargetAtom(couple, fTheGamma, GammaEnergy,
357  aDynamicGamma->GetLogKineticEnergy() );
358  // Atomic number
359  const G4int Z = anElement->GetZasInt();
360  const G4int A = SelectIsotopeNumber(anElement);
361  const G4double iZ13 = 1./anElement->GetIonisation()->GetZ3();
362  const G4double targetMass = G4NucleiProperties::GetNuclearMass(A, Z);
363 
364  const G4double NuThreshold = 2.0 * ( (LeptonMass2)/targetMass + LeptonMass);
365  // No conversion possible below nuclear threshold
366  if ( GammaEnergy <= NuThreshold) { return; }
367 
368  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
369 
370  // itriplet : true -- triplet, false -- nuclear.
371  G4bool itriplet = false;
372  if (fConversionType == 1) {
373  itriplet = false;
374  } else if (fConversionType == 2) {
375  itriplet = true;
376  if ( GammaEnergy <= TrThreshold ) return;
377  } else if ( GammaEnergy > TrThreshold ) {
378  // choose triplet or nuclear from a triplet/nuclear=1/Z
379  // total cross section ratio.
380  // approximate at low energies !
381  if(rndmEngine->flat()*(Z+1) < 1.) {
382  itriplet = true;
383  }
384  }
385 
386  //
387  const G4double RecoilMass = itriplet ? ElectronMass : targetMass;
388  const G4double RecoilMass2 = RecoilMass*RecoilMass;
389  const G4double sCMS = 2.*RecoilMass*GammaEnergy + RecoilMass2;
390  const G4double sCMSPlusRM2 = sCMS + RecoilMass2;
391  const G4double sqrts = std::sqrt(sCMS);
392  const G4double isqrts2 = 1./(2.*sqrts);
393  //
394  const G4double PairInvMassMax = sqrts-RecoilMass;
395  const G4double PairInvMassRange = PairInvMassMax/PairInvMassMin;
396  const G4double lnPairInvMassRange = G4Log(PairInvMassRange);
397 
398  // initial state. Defines z axis of "0" frame as along photon propagation.
399  // Since CMS(0., 0., GammaEnergy, GammaEnergy+RecoilMass) set some constants
400  const G4double betaCMS = G4LorentzVector(0.0,0.0,GammaEnergy,GammaEnergy+RecoilMass).beta();
401 
402  // maximum value of pdf
403  const G4double EffectiveZ = iraw ? 0.5 : Z;
404  const G4double Threshold = itriplet ? TrThreshold : NuThreshold;
405  const G4double AvailableEnergy = GammaEnergy - Threshold;
406  const G4double LogAvailableEnergy = G4Log(AvailableEnergy);
407  //
408  const G4double MaxDiffCross = itriplet
409  ? MaxDiffCrossSection(tr[fConvMode],
410  EffectiveZ, AvailableEnergy, LogAvailableEnergy)
411  : MaxDiffCrossSection(nu[fConvMode],
412  EffectiveZ, AvailableEnergy, LogAvailableEnergy);
413  //
414  // 50% safety marging factor
415  const G4double ymax = 1.5 * MaxDiffCross;
416  // x1 bounds
417  const G4double xu1 = (LogAvailableEnergy > para[fConvMode][2][0])
418  ? para[fConvMode][0][0] +
419  para[fConvMode][1][0]*LogAvailableEnergy
420  : para[fConvMode][0][0] +
421  para[fConvMode][2][0]*para[fConvMode][1][0];
422  const G4double xl1 = (LogAvailableEnergy > para[fConvMode][2][1])
423  ? para[fConvMode][0][1] +
424  para[fConvMode][1][1]*LogAvailableEnergy
425  : para[fConvMode][0][1] +
426  para[fConvMode][2][1]*para[fConvMode][1][1];
427  //
428  G4LorentzVector Recoil;
429  G4LorentzVector LeptonPlus;
430  G4LorentzVector LeptonMinus;
431  G4double pdf = 0.;
432 
433  G4double rndmv6[6];
434  // START Sampling
435  do {
436 
437  rndmEngine->flatArray(6, rndmv6);
438 
440  // pdf pow(x,c) with c = 1.4
441  // integral y = pow(x,(c+1))/(c+1) @ x = 1 => y = 1 /(1+c)
442  // invCdf exp( log(y /* *( c + 1.0 )/ (c + 1.0 ) */ ) /( c + 1.0) )
444  const G4double X1 =
445  G4Exp(G4Log(rndmv6[0])/(correctionIndex + 1.0));
446 
447  const G4double x0 = G4Exp(xl1 + (xu1 - xl1)*rndmv6[1]);
448  const G4double dum0 = 1./(1.+x0);
449  const G4double cosTheta = (x0-1.)*dum0;
450  const G4double sinTheta = std::sqrt(4.*x0)*dum0;
451 
452  const G4double PairInvMass = PairInvMassMin*G4Exp(X1*X1*lnPairInvMassRange);
453 
454  // G4double rndmv3[3];
455  // rndmEngine->flatArray(3, rndmv3);
456 
457  // cos and sin theta-lepton
458  const G4double cosThetaLept = std::cos(pi*rndmv6[2]);
459  // sin(ThetaLept) is always in [0,+1] if ThetaLept is in [0,pi]
460  const G4double sinThetaLept = std::sqrt((1.-cosThetaLept)*(1.+cosThetaLept));
461  // cos and sin phi-lepton
462  const G4double cosPhiLept = std::cos(twoPi*rndmv6[3]-pi);
463  // sin(PhiLept) is in [-1,0] if PhiLept in [-pi,0) and
464  // is in [0,+1] if PhiLept in [0,+pi]
465  const G4double sinPhiLept = std::copysign(std::sqrt((1.-cosPhiLept)*(1.+cosPhiLept)),rndmv6[3]-0.5);
466  // cos and sin phi
467  const G4double cosPhi = std::cos(twoPi*rndmv6[4]-pi);
468  const G4double sinPhi = std::copysign(std::sqrt((1.-cosPhi)*(1.+cosPhi)),rndmv6[4]-0.5);
469 
471  // frames:
472  // 3 : the laboratory Lorentz frame, Geant4 axes definition
473  // 0 : the laboratory Lorentz frame, axes along photon direction and polarisation
474  // 1 : the center-of-mass Lorentz frame
475  // 2 : the pair Lorentz frame
477 
478  // in the center-of-mass frame
479 
480  const G4double RecEnergyCMS = (sCMSPlusRM2-PairInvMass*PairInvMass)*isqrts2;
481  const G4double LeptonEnergy2 = PairInvMass*0.5;
482 
483  // New way of calucaltion thePRecoil to avoid underflow
484  G4double abp = std::max((2.0*GammaEnergy*RecoilMass -
485  PairInvMass*PairInvMass + 2.0*PairInvMass*RecoilMass)*
486  (2.0*GammaEnergy*RecoilMass -
487  PairInvMass*PairInvMass - 2.0*PairInvMass*RecoilMass),0.0);
488 
489  G4double thePRecoil = std::sqrt(abp) * isqrts2;
490 
491  // back to the center-of-mass frame
492  Recoil.set( thePRecoil*sinTheta*cosPhi,
493  thePRecoil*sinTheta*sinPhi,
494  thePRecoil*cosTheta,
495  RecEnergyCMS);
496 
497  // in the pair frame
498  const G4double thePLepton = std::sqrt( (LeptonEnergy2-LeptonMass)
499  *(LeptonEnergy2+LeptonMass));
500 
501  LeptonPlus.set(thePLepton*sinThetaLept*cosPhiLept,
502  thePLepton*sinThetaLept*sinPhiLept,
503  thePLepton*cosThetaLept,
504  LeptonEnergy2);
505 
506  LeptonMinus.set(-LeptonPlus.x(),
507  -LeptonPlus.y(),
508  -LeptonPlus.z(),
509  LeptonEnergy2);
510 
511 
512  // Normalisation of final state phase space:
513  // Section 47 of Particle Data Group, Chin. Phys. C, 40, 100001 (2016)
514  // const G4double Norme = Recoil1.vect().mag() * LeptonPlus2.vect().mag();
515  const G4double Norme = Recoil.vect().mag() * LeptonPlus.vect().mag();
516 
517  // e+, e- to CMS frame from pair frame
518 
519  // boost vector from Pair to CMS
520  const G4ThreeVector pair2cms =
521  G4LorentzVector( -Recoil.x(), -Recoil.y(), -Recoil.z(),
522  sqrts-RecEnergyCMS).boostVector();
523 
524  LeptonPlus.boost(pair2cms);
525  LeptonMinus.boost(pair2cms);
526 
527  // back to the laboratory frame (make use of the CMS(0,0,Eg,Eg+RM)) form
528 
529  Recoil.boostZ(betaCMS);
530  LeptonPlus.boostZ(betaCMS);
531  LeptonMinus.boostZ(betaCMS);
532 
533  // Jacobian factors
534  const G4double Jacob0 = x0*dum0*dum0;
535  const G4double Jacob1 = 2.*X1*lnPairInvMassRange*PairInvMass;
536  const G4double Jacob2 = std::abs(sinThetaLept);
537 
538  const G4double EPlus = LeptonPlus.t();
539  const G4double PPlus = LeptonPlus.vect().mag();
540  const G4double sinThetaPlus = LeptonPlus.vect().perp()/PPlus;
541  const G4double cosThetaPlus = LeptonPlus.vect().cosTheta();
542 
543  const G4double pPX = LeptonPlus.x();
544  const G4double pPY = LeptonPlus.y();
545  const G4double dum1 = 1./std::sqrt( pPX*pPX + pPY*pPY );
546  const G4double cosPhiPlus = pPX*dum1;
547  const G4double sinPhiPlus = pPY*dum1;
548 
549  // denominators:
550  // the two cancelling leading terms for forward emission at high energy, removed
551  const G4double elMassCTP = LeptonMass*cosThetaPlus;
552  const G4double ePlusSTP = EPlus*sinThetaPlus;
553  const G4double DPlus = (elMassCTP*elMassCTP + ePlusSTP*ePlusSTP)
554  /(EPlus + PPlus*cosThetaPlus);
555 
556  const G4double EMinus = LeptonMinus.t();
557  const G4double PMinus = LeptonMinus.vect().mag();
558  const G4double sinThetaMinus = LeptonMinus.vect().perp()/PMinus;
559  const G4double cosThetaMinus = LeptonMinus.vect().cosTheta();
560 
561  const G4double ePX = LeptonMinus.x();
562  const G4double ePY = LeptonMinus.y();
563  const G4double dum2 = 1./std::sqrt( ePX*ePX + ePY*ePY );
564  const G4double cosPhiMinus = ePX*dum2;
565  const G4double sinPhiMinus = ePY*dum2;
566 
567  const G4double elMassCTM = LeptonMass*cosThetaMinus;
568  const G4double eMinSTM = EMinus*sinThetaMinus;
569  const G4double DMinus = (elMassCTM*elMassCTM + eMinSTM*eMinSTM)
570  /(EMinus + PMinus*cosThetaMinus);
571 
572  // cos(phiMinus-PhiPlus)
573  const G4double cosdPhi = cosPhiPlus*cosPhiMinus + sinPhiPlus*sinPhiMinus;
574  const G4double PRec = Recoil.vect().mag();
575  const G4double q2 = PRec*PRec;
576 
577  const G4double BigPhi = -LeptonMass2 / (GammaEnergy*GammaEnergy2 * q2*q2);
578 
579  G4double FormFactor = 1.;
580  if (!iraw) {
581  if (itriplet) {
582  const G4double qun = factor1*iZ13*iZ13;
583  const G4double nun = qun * PRec;
584  if (nun < 1.) {
585  FormFactor = (nun < 0.01) ? (13.8-55.4*std::sqrt(nun))*nun
586  : std::sqrt(1-(nun-1)*(nun-1));
587  } // else FormFactor = 1 by default
588  } else {
589  const G4double dum3 = 217.*PRec*iZ13;
590  const G4double AFF = 1./(1. + dum3*dum3);
591  FormFactor = (1.-AFF)*(1-AFF);
592  }
593  } // else FormFactor = 1 by default
594 
595  G4double betheheitler;
596  if (GammaPolarizationMag==0.) {
597  const G4double pPlusSTP = PPlus*sinThetaPlus;
598  const G4double pMinusSTM = PMinus*sinThetaMinus;
599  const G4double pPlusSTPperDP = pPlusSTP/DPlus;
600  const G4double pMinusSTMperDM = pMinusSTM/DMinus;
601  const G4double dunpol = BigPhi*(
602  pPlusSTPperDP *pPlusSTPperDP *(4.*EMinus*EMinus-q2)
603  + pMinusSTMperDM*pMinusSTMperDM*(4.*EPlus*EPlus - q2)
604  + 2.*pPlusSTPperDP*pMinusSTMperDM*cosdPhi
605  *(4.*EPlus*EMinus + q2 - 2.*GammaEnergy2)
606  - 2.*GammaEnergy2*(pPlusSTP*pPlusSTP+pMinusSTM*pMinusSTM)/(DMinus*DPlus));
607  betheheitler = dunpol * factor;
608  } else {
609  const G4double pPlusSTP = PPlus*sinThetaPlus;
610  const G4double pMinusSTM = PMinus*sinThetaMinus;
611  const G4double pPlusSTPCPPperDP = pPlusSTP*cosPhiPlus/DPlus;
612  const G4double pMinusSTMCPMperDM = pMinusSTM*cosPhiMinus/DMinus;
613  const G4double caa = 2.*(EPlus*pMinusSTMCPMperDM+EMinus*pPlusSTPCPPperDP);
614  const G4double cbb = pMinusSTMCPMperDM-pPlusSTPCPPperDP;
615  const G4double ccc = (pPlusSTP*pPlusSTP + pMinusSTM*pMinusSTM
616  +2.*pPlusSTP*pMinusSTM*cosdPhi)/ (DMinus*DPlus);
617  const G4double dtot= 2.*BigPhi*( caa*caa - q2*cbb*cbb - GammaEnergy2*ccc);
618  betheheitler = dtot * factor;
619  }
620  //
621  const G4double cross = Norme * Jacob0 * Jacob1 * Jacob2 * betheheitler
622  * FormFactor * RecoilMass / sqrts;
623  pdf = cross * (xu1 - xl1) / G4Exp(correctionIndex*G4Log(X1)); // cond1;
624  } while ( pdf < ymax * rndmv6[5] );
625  // END of Sampling
626 
627  if ( fVerbose > 2 ) {
628  G4double recul = std::sqrt(Recoil.x()*Recoil.x()+Recoil.y()*Recoil.y()
629  +Recoil.z()*Recoil.z());
630  G4cout << "BetheHeitler5DModel GammaEnergy= " << GammaEnergy
631  << " PDF= " << pdf << " ymax= " << ymax
632  << " recul= " << recul << G4endl;
633  }
634 
635  // back to Geant4 system
636 
637  if ( fVerbose > 4 ) {
638  G4cout << "BetheHeitler5DModel GammaDirection " << GammaDirection << G4endl;
639  G4cout << "BetheHeitler5DModel GammaPolarization " << GammaPolarization << G4endl;
640  G4cout << "BetheHeitler5DModel GammaEnergy " << GammaEnergy << G4endl;
641  G4cout << "BetheHeitler5DModel Conv "
642  << (itriplet ? "triplet" : "nucl") << G4endl;
643  }
644 
645  if (GammaPolarizationMag == 0.0) {
646  // set polarization axis orthohonal to direction
647  GammaPolarization = GammaDirection.orthogonal().unit();
648  } else {
649  // GammaPolarization not a unit vector
650  GammaPolarization /= GammaPolarizationMag;
651  }
652 
653  // The unit norm vector that is orthogonal to the two others
654  G4ThreeVector yGrec = GammaDirection.cross(GammaPolarization);
655 
656  // rotation from gamma ref. sys. to World
657  G4RotationMatrix GtoW(GammaPolarization,yGrec,GammaDirection);
658 
659  Recoil.transform(GtoW);
660  LeptonPlus.transform(GtoW);
661  LeptonMinus.transform(GtoW);
662 
663  if ( fVerbose > 2 ) {
664  G4cout << "BetheHeitler5DModel Recoil " << Recoil.x() << " " << Recoil.y() << " " << Recoil.z()
665  << " " << Recoil.t() << " " << G4endl;
666  G4cout << "BetheHeitler5DModel LeptonPlus " << LeptonPlus.x() << " " << LeptonPlus.y() << " "
667  << LeptonPlus.z() << " " << LeptonPlus.t() << " " << G4endl;
668  G4cout << "BetheHeitler5DModel LeptonMinus " << LeptonMinus.x() << " " << LeptonMinus.y() << " "
669  << LeptonMinus.z() << " " << LeptonMinus.t() << " " << G4endl;
670  }
671 
672  // Create secondaries
673  G4DynamicParticle* aParticle1 = new G4DynamicParticle(fLepton1,LeptonMinus);
674  G4DynamicParticle* aParticle2 = new G4DynamicParticle(fLepton2,LeptonPlus);
675 
676  // create G4DynamicParticle object for the particle3 ( recoil )
677  G4ParticleDefinition* RecoilPart;
678  if (itriplet) {
679  // triplet
680  RecoilPart = fTheElectron;
681  } else{
682  RecoilPart = theIonTable->GetIon(Z, A, 0);
683  }
684  G4DynamicParticle* aParticle3 = new G4DynamicParticle(RecoilPart,Recoil);
685 
686  // Fill output vector
687  fvect->push_back(aParticle1);
688  fvect->push_back(aParticle2);
689  fvect->push_back(aParticle3);
690 
691  // kill incident photon
692  fParticleChange->SetProposedKineticEnergy(0.);
693  fParticleChange->ProposeTrackStatus(fStopAndKill);
694 }
695 
696 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....