ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4hPairProductionModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4hPairProductionModel.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: G4hPairProductionModel
33 //
34 // Author: Vladimir Ivanchenko on base of G4MuPairProductionModel
35 //
36 // Creation date: 28.02.2008
37 //
38 // Modifications:
39 //
40 
41 //
42 // Class Description:
43 //
44 //
45 // -------------------------------------------------------------------
46 //
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 
50 #include "G4hPairProductionModel.hh"
51 #include "G4PhysicalConstants.hh"
52 #include "G4Log.hh"
53 #include "G4Exp.hh"
54 
55 using namespace std;
56 
57 G4hPairProductionModel::G4hPairProductionModel(const G4ParticleDefinition* p,
58  const G4String& nam)
59  : G4MuPairProductionModel(p, nam)
60 {}
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 
64 G4hPairProductionModel::~G4hPairProductionModel()
65 {}
66 
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 
69 G4double G4hPairProductionModel::ComputeDMicroscopicCrossSection(
70  G4double tkin,
71  G4double Z,
72  G4double pairEnergy)
73 // differential cross section
74 {
75  static const G4double bbbtf= 183. ;
76  static const G4double bbbh = 202.4 ;
77  static const G4double g1tf = 1.95e-5 ;
78  static const G4double g2tf = 5.3e-5 ;
79  static const G4double g1h = 4.4e-5 ;
80  static const G4double g2h = 4.8e-5 ;
81 
82  G4double totalEnergy = tkin + particleMass;
83  G4double residEnergy = totalEnergy - pairEnergy;
84  G4double massratio = particleMass/electron_mass_c2 ;
85  G4double massratio2 = massratio*massratio ;
86  G4double cross = 0.;
87 
88  G4double c3 = 0.75*sqrte*particleMass;
89  if (residEnergy <= c3*z13) { return cross; }
90 
91  G4double c7 = 4.*CLHEP::electron_mass_c2;
92  G4double c8 = 6.*particleMass*particleMass;
93  G4double alf = c7/pairEnergy;
94  G4double a3 = 1. - alf;
95  if (a3 <= 0.) { return cross; }
96 
97  // zeta calculation
98  G4double bbb,g1,g2;
99  if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
100  else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
101 
102  G4double zeta = 0.;
103  G4double zeta1 =
104  0.073*G4Log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26;
105  if ( zeta1 > 0.)
106  {
107  G4double zeta2 =
108  0.058*G4Log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14;
109  zeta = zeta1/zeta2 ;
110  }
111 
112  G4double z2 = Z*(Z+zeta);
113  G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
114  G4double a0 = totalEnergy*residEnergy;
115  G4double a1 = pairEnergy*pairEnergy/a0;
116  G4double bet = 0.5*a1;
117  G4double xi0 = 0.25*massratio2*a1;
118  G4double del = c8/a0;
119 
120  G4double rta3 = sqrt(a3);
121  G4double tmnexp = alf/(1. + rta3) + del*rta3;
122  if(tmnexp >= 1.0) { return cross; }
123 
124  G4double tmn = G4Log(tmnexp);
125  G4double sum = 0.;
126 
127  // Gaussian integration in ln(1-ro) ( with 8 points)
128  for (G4int i=0; i<8; i++)
129  {
130  G4double a4 = G4Exp(tmn*xgi[i]); // a4 = (1.-asymmetry)
131  G4double a5 = a4*(2.-a4) ;
132  G4double a6 = 1.-a5 ;
133  G4double a7 = 1.+a6 ;
134  G4double a9 = 3.+a6 ;
135  G4double xi = xi0*a5 ;
136  G4double xii = 1./xi ;
137  G4double xi1 = 1.+xi ;
138  G4double screen = screen0*xi1/a5 ;
139  G4double yeu = 5.-a6+4.*bet*a7 ;
140  G4double yed = 2.*(1.+3.*bet)*G4Log(3.+xii)-a6-a1*(2.-a6) ;
141  G4double ye1 = 1.+yeu/yed ;
142  G4double ale=G4Log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
143  G4double cre = 0.5*G4Log(1.+2.25*z23*xi1*ye1/massratio2) ;
144  G4double be;
145 
146  if (xi <= 1.e3) {
147  be = ((2.+a6)*(1.+bet)+xi*a9)*G4Log(1.+xii)+(a5-bet)/xi1-a9;
148  } else {
149  be = (3.-a6+a1*a7)/(2.*xi);
150  }
151  G4double fe = (ale-cre)*be;
152  if ( fe < 0.) { fe = 0.; }
153 
154  G4double ymu = 4.+a6 +3.*bet*a7 ;
155  G4double ymd = a7*(1.5+a1)*G4Log(3.+xi)+1.-1.5*a6 ;
156  G4double ym1 = 1.+ymu/ymd ;
157  G4double alm_crm = G4Log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
158  G4double a10,bm;
159  if ( xi >= 1.e-3)
160  {
161  a10 = (1.+a1)*a5 ;
162  bm = (a7*(1.+1.5*bet)-a10*xii)*G4Log(xi1)+xi*(a5-bet)/xi1+a10;
163  } else {
164  bm = (5.-a6+bet*a9)*(xi/2.);
165  }
166 
167  G4double fm = alm_crm*bm;
168  if ( fm < 0.) { fm = 0.; }
169 
170  sum += wgi[i]*a4*(fe+fm/massratio2);
171  }
172 
173  cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
174 
175  return cross;
176 }
177 
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......