ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizedPairProductionCrossSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PolarizedPairProductionCrossSection.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: G4PolarizedPairProductionCrossSection
33 //
34 // Author: Andreas Schaelicke on the base of Karim Laihems code
35 //
36 // Creation date: 16.08.2006
37 //
38 
39 #include "G4PolarizedPairProductionCrossSection.hh"
40 #include "G4PhysicalConstants.hh"
41 
42 G4bool G4PolarizedPairProductionCrossSection::scrnInitialized=false;
43 G4double G4PolarizedPairProductionCrossSection::SCRN [3][20];
44 // screening function lookup table;
45 
46 void G4PolarizedPairProductionCrossSection::InitializeMe()
47 {
48  if (!scrnInitialized) {
49  SCRN [1][1]= 0.5 ; SCRN [2][1] = 0.0145;
50  SCRN [1][2]= 1.0 ; SCRN [2][2] = 0.0490;
51  SCRN [1][3]= 2.0 ; SCRN [2][3] = 0.1400;
52  SCRN [1][4]= 4.0 ; SCRN [2][4] = 0.3312;
53  SCRN [1][5]= 8.0 ; SCRN [2][5] = 0.6758;
54  SCRN [1][6]= 15.0 ; SCRN [2][6] = 1.126;
55  SCRN [1][7]= 20.0 ; SCRN [2][7] = 1.367;
56  SCRN [1][8]= 25.0 ; SCRN [2][8] = 1.564;
57  SCRN [1][9]= 30.0 ; SCRN [2][9] = 1.731;
58  SCRN [1][10]= 35.0 ; SCRN [2][10]= 1.875;
59  SCRN [1][11]= 40.0 ; SCRN [2][11]= 2.001;
60  SCRN [1][12]= 45.0 ; SCRN [2][12]= 2.114;
61  SCRN [1][13]= 50.0 ; SCRN [2][13]= 2.216;
62  SCRN [1][14]= 60.0 ; SCRN [2][14]= 2.393;
63  SCRN [1][15]= 70.0 ; SCRN [2][15]= 2.545;
64  SCRN [1][16]= 80.0 ; SCRN [2][16]= 2.676;
65  SCRN [1][17]= 90.0 ; SCRN [2][17]= 2.793;
66  SCRN [1][18]= 100.0 ; SCRN [2][18]= 2.897;
67  SCRN [1][19]= 120.0 ; SCRN [2][19]= 3.078;
68 
69  scrnInitialized=true;
70  }
71 }
72 
73 G4PolarizedPairProductionCrossSection::G4PolarizedPairProductionCrossSection()
74 {
75  InitializeMe();
76 }
77 
78 
79 void G4PolarizedPairProductionCrossSection::Initialize(G4double aGammaE,
80  G4double aLept0E,
81  G4double sintheta,
82  const G4StokesVector & beamPol,
83  const G4StokesVector & /*p1*/,
84  G4int /*flag*/)
85 {
86  G4double aLept1E = aGammaE - aLept0E;
87 
88  G4double Stokes_P3 = beamPol.z() ;
89  // **************************************************************************
90 
91  G4double m0_c2 = electron_mass_c2;
92  G4double Lept0E = aLept0E/m0_c2+1., Lept0E2 = Lept0E * Lept0E ;
93  G4double GammaE = aGammaE/m0_c2;
94  G4double Lept1E = aLept1E/m0_c2-1., Lept1E2 = Lept1E * Lept1E ;
95 
96 
97  // const G4Element* theSelectedElement = theModel->SelectedAtom();
98 
99  // ******* Gamma Transvers Momentum
100 
101  G4double TMom = std::sqrt(Lept0E2 -1.)* sintheta, u = TMom , u2 =u * u ;
102  G4double Xsi = 1./(1.+u2) , Xsi2 = Xsi * Xsi ;
103 
104  // G4double theZ = theSelectedElement->GetZ();
105 
106  // G4double fCoul = theSelectedElement->GetfCoulomb();
107  G4double delta = 12. * std::pow(theZ, 1./3.) * Lept0E * Lept1E * Xsi / (121. * GammaE);
108  G4double GG=0.;
109 
110  if(delta < 0.5) {
111  GG = std::log(2.* Lept0E * Lept1E / GammaE) - 2. - fCoul;
112  }
113  else if ( delta < 120.) {
114  for (G4int j=2; j<=19; j++) {
115  if(SCRN[1][j] >= delta) {
116  GG =std::log(2. * Lept0E * Lept1E / GammaE) - 2. - fCoul
117  -(SCRN[2][j-1]+(delta-SCRN[1][j-1])*(SCRN[2][j]-SCRN[2][j-1])/(SCRN[1][j]-SCRN[1][j-1]));
118  break;
119  }
120  }
121  }
122  else {
123  G4double alpha_sc = (111. * std::pow(theZ, -1./3.)) / Xsi;
124  GG = std::log(alpha_sc)- 2. - fCoul;
125  }
126 
127  if(GG<-1.) GG=-1.; // *KL* do we need this ?!
128 
129 
130  G4double I_Lepton = (Lept0E2 + Lept1E2)*(3+2*GG) + 2. * Lept0E * Lept1E * (1. + 4. * u2 * Xsi2 * GG);
131 
132  // G4double D_Lepton1 = -8 * Lept0E * Lept1E * u2 * Xsi2 * GG / I_Lepton;
133 
134  G4double L_Lepton1 = GammaE * ((Lept0E - Lept1E) * (3. + 2. * GG)+2 * Lept1E * (1. + 4. * u2 * Xsi2 * GG))/I_Lepton;
135 
136  G4double T_Lepton1 = 4. * GammaE * Lept1E * Xsi * u * (1. - 2. * Xsi) * GG / I_Lepton ;
137 
138 
139  G4double Stokes_S1 = (Stokes_P3 * T_Lepton1) ;
140  G4double Stokes_S2 = 0.;
141  G4double Stokes_S3 = (Stokes_P3 * L_Lepton1) ;
142 
143 
144  theFinalElectronPolarization.setX(Stokes_S1);
145  theFinalElectronPolarization.setY(Stokes_S2);
146  theFinalElectronPolarization.setZ(Stokes_S3);
147 
148  if(theFinalElectronPolarization.mag2()>1.) {
149  G4cout<<" WARNING in pol-conv theFinalElectronPolarization \n";
150  G4cout
151  <<"\t"<<theFinalElectronPolarization
152  <<"\t GG\t"<<GG
153  <<"\t delta\t"<<delta
154  <<G4endl;
155  theFinalElectronPolarization.setX(0.);
156  theFinalElectronPolarization.setY(0.);
157  theFinalElectronPolarization.setZ(Stokes_S3);
158  if(Stokes_S3>1.) theFinalElectronPolarization.setZ(1.);
159  }
160 
161 
162  G4double L_Lepton2 = GammaE * ((Lept1E - Lept0E) * (3. + 2. * GG)+2 * Lept0E * (1. + 4. * u2 * Xsi2 * GG))/I_Lepton;
163 
164  G4double T_Lepton2 = 4. * GammaE * Lept0E * Xsi * u * (1. - 2. * Xsi) * GG / I_Lepton ;
165 
166  G4double Stokes_SS1 = (Stokes_P3 * T_Lepton2) ;
167  G4double Stokes_SS2 = 0.;
168  G4double Stokes_SS3 = (Stokes_P3 * L_Lepton2) ;
169 
170  theFinalPositronPolarization.SetPhoton();
171 
172  theFinalPositronPolarization.setX(Stokes_SS1);
173  theFinalPositronPolarization.setY(Stokes_SS2);
174  theFinalPositronPolarization.setZ(Stokes_SS3);
175 
176  if(theFinalPositronPolarization.mag2()>1.) {
177  G4cout<<" WARNING in pol-conv theFinalPositronPolarization \n";
178  G4cout
179  <<"\t"<<theFinalPositronPolarization
180  <<"\t GG\t"<<GG
181  <<"\t delta\t"<<delta
182  <<G4endl;
183  }
184 }
185 
186 G4double G4PolarizedPairProductionCrossSection::XSection(const G4StokesVector & /*pol2*/,
187  const G4StokesVector & /*pol3*/)
188 {
189  G4cout<<"ERROR dummy routine G4PolarizedPairProductionCrossSection::XSection called \n";
190  return 0.;
191 }
192 
193  // return expected mean polarisation
194 G4StokesVector G4PolarizedPairProductionCrossSection::GetPol2()
195 {
196  // electron/positron
197  return theFinalElectronPolarization;
198 }
199 G4StokesVector G4PolarizedPairProductionCrossSection::GetPol3()
200 {
201  // photon
202  return theFinalPositronPolarization;;
203 }
204 
205