ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eeCrossSections.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eeCrossSections.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 header file
30 //
31 //
32 // File name: G4eeCrossSections
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 25.10.2003
37 //
38 // Modifications:
39 // 10.07.2008 Updated for PDG Jour. Physics, G33, 1 (2006)
40 //
41 // -------------------------------------------------------------------
42 //
43 
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
48 #include "G4eeCrossSections.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4PionPlus.hh"
52 #include "G4PionMinus.hh"
53 #include "G4PionZero.hh"
54 #include "G4Eta.hh"
55 #include "G4KaonPlus.hh"
56 #include "G4KaonMinus.hh"
57 #include "G4KaonZeroLong.hh"
58 #include "G4PhysicsLinearVector.hh"
59 
60 #include <iostream>
61 #include <fstream>
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
65 using namespace std;
66 
67 G4eeCrossSections::G4eeCrossSections()
68 {
69  Initialise();
70 }
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
74 G4eeCrossSections::~G4eeCrossSections()
75 {}
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 
79 void G4eeCrossSections::Initialise()
80 {
81  MsPi = G4PionPlus::PionPlus()->GetPDGMass();
82  MsPi0= G4PionZero::PionZero()->GetPDGMass();
83  MsEta= G4Eta::Eta()->GetPDGMass();
84  MsEtap=957.78*MeV;
85  MsKs = G4KaonZeroLong::KaonZeroLong()->GetPDGMass();
86  MsKc = G4KaonPlus::KaonPlus()->GetPDGMass();
87  MsRho= 775.5*MeV;
88  MsOm = 782.62*MeV;
89  MsF0 = 980.0*MeV;
90  MsA0 = 984.7*MeV;
91  MsPhi= 1019.46*MeV;
92  MsK892 = 891.66*MeV;
93  MsK0892 = 896.0*MeV;
94  GRho = 149.4*MeV;
95  GOm = 8.49*MeV;
96  GPhi = 4.26*MeV;
97  GK892 = 50.8*MeV;
98  GK0892 = 50.3*MeV;
99  PhRho = 0.0;
100  PhOm = 0.0;
101  PhPhi = 155.0*degree;
102  PhRhoPi = 186.0*degree;
103 
104  BrRhoPiG = 4.5e-4;
105  BrRhoPi0G= 6.8e-4;
106  BrRhoEtaG= 2.95e-4;
107  BrRhoEe = 4.7e-5;
108  BrOm3Pi = 0.891;
109  BrOmPi0G= 0.089;
110  BrOmEtaG= 4.9e-4;
111  BrOm2Pi = 0.017;
112  PhOm2Pi = 90.0;
113  BrOmEe = 7.18e-5;
114  BrPhi2Kc = 0.492;
115  BrPhiKsKl= 0.34;
116  BrPhi3Pi = 0.153;
117  BrPhiPi0G= 1.25e-3;
118  BrPhiEtaG= 1.301e-2;
119  BrPhi2Pi = 7.3e-5;
120  PhPhi2Pi = -20.0*degree;
121  BrPhiEe = 2.97e-4;
122 
123  MsRho3 = MsRho*MsRho*MsRho;
124  MsOm3 = MsOm*MsOm*MsOm;
125  MsPhi3 = MsPhi*MsPhi*MsPhi;
126 
127  MeVnb = 3.8938e+11*nanobarn;
128  Alpha = fine_structure_const;
129 
130  AOmRho = 3.0;
131  ARhoPRho = 0.72;
132  cterm=0.;
133  mssig = 600.*MeV;
134  gsig = 500.*MeV;
135  brsigpipi = 1.;
136 
137  msrho1450 = 1459.*MeV;
138  msrho1700 = 1688.8*MeV;
139  grho1450 = 171.*MeV;
140  grho1700 = 161.*MeV;
141  arhoompi0 = 1.;
142  arho1450ompi0 = 1.;
143  arho1700ompi0 = 1.;
144  phrhoompi0 = 0.;
145  phrho1450ompi0 = pi;
146  phrho1700ompi0 = 0.;
147  aomrhopi0 = 1.;
148  phomrhopi0 = 0.;
149  arhopi0pi0g = 0.;
150  aompi0pi0g = 0.;
151  phrhopi0pi0g = 0.;
152  phompi0pi0g = 0.;
153  brrho1450ompi0 = 0.02;
154  brrho1450pipi = 0.50;
155  brrho1700ompi0 = 1.0;
156  brrho1700pipi = 0.02;
157  aphirhopi0 = 1.;
158  phphirhopi0 = pi;
159  arhosigg = 0.;
160  phrhosigg = 0.;
161  aomsigg = 0.;
162  phomsigg = 0.;
163 
164  G4String w0, w1, w2;
165  ph3p = 0;
166 
167  /*
168  G4double emin, emax;
169  G4int nbins;
170  const G4String fname = "wrhopi.wid";
171  ifstream fi(fname.c_str());
172  fi >> w0 >> nbins >> w1 >> emin >> w2 >> emax;
173  emin *= MeV;
174  emax *= MeV;
175  ph3p = new G4PhysicsLinearVector(emin,emax,nbins);
176  G4int nlines = nbins/5;
177  G4double s0, s1, s2, s3, s4;
178  for(G4int i=0; i<nlines; i++) {
179  fi >> s0 >> s1 >> s2 >> s3 >> s4;
180  ph3p->PutValue(5*i, s0);
181  ph3p->PutValue(5*i + 1, s1);
182  ph3p->PutValue(5*i + 2, s2);
183  ph3p->PutValue(5*i + 3, s3);
184  ph3p->PutValue(5*i + 4, s4);
185  }
186  fi.close();
187  */
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
191 
192 G4double G4eeCrossSections::CrossSection2pi(G4double e)
193 {
194  complex<G4double> xr(cos(PhRho),sin(PhRho));
195  complex<G4double> xo(cos(PhOm2Pi),sin(PhOm2Pi));
196  complex<G4double> xf(cos(PhPhi2Pi),sin(PhPhi2Pi));
197 
198  G4double s_inv = e*e;
199  complex<G4double> drho = DpRho(e);
200  complex<G4double> dom = DpOm(e);
201  complex<G4double> dphi = DpPhi(e);
202 
203  complex<G4double> amp =
204  sqrt(Width2p(s_inv,MsRho,GRho,1.0,MsPi)*MsRho3*BrRhoEe*GRho)*xr/drho
205  + sqrt(Width2p(s_inv,MsOm,GOm,BrOm2Pi,MsPi)*MsOm3*BrOmEe*GOm)*xo/dom
206  + sqrt(Width2p(s_inv,MsPhi,GPhi,BrPhi2Pi,MsPi)*MsPhi3*BrPhiEe*GPhi)*xf/dphi;
207 
208  G4double cross = 12.0*pi*MeVnb*norm(amp)/(e*s_inv);
209 
210  return cross;
211 }
212 
213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
214 
215 G4double G4eeCrossSections::CrossSection3pi(G4double e)
216 {
217  complex<G4double> xf(cos(PhPhi2Pi),sin(PhPhi));
218 
219  G4double s_inv = e*e;
220  complex<G4double> dom = DpOm(e);
221  complex<G4double> dphi = DpPhi(e);
222 
223  complex<G4double> amp =
224  sqrt(Width3p(s_inv,MsOm,GOm,BrOm3Pi)*MsOm3*BrOmEe*GOm)/dom
225  + sqrt(Width3p(s_inv,MsPhi,GPhi,BrPhi3Pi)*MsPhi3*BrPhiEe*GPhi)*xf/dphi;
226 
227  G4double cross = 12.0*pi*MeVnb*norm(amp)/(e*s_inv);
228 
229  return cross;
230 }
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 
234 G4double G4eeCrossSections::CrossSectionPi0G(G4double e)
235 {
236  complex<G4double> xf(cos(PhPhi),sin(PhPhi));
237 
238  G4double s_inv = e*e;
239  complex<G4double> drho = DpRho(e);
240  complex<G4double> dom = DpOm(e);
241  complex<G4double> dphi = DpPhi(e);
242 
243  complex<G4double> amp =
244  sqrt(WidthPg(s_inv,MsRho,GRho,BrRhoPi0G,MsPi0)*MsRho3*BrRhoEe*GRho)/drho
245  + sqrt(WidthPg(s_inv,MsOm,GOm,BrOmPi0G,MsPi0)*MsOm3*BrOmEe*GOm)/dom
246  + sqrt(WidthPg(s_inv,MsPhi,GPhi,BrPhiPi0G,MsPi0)*MsPhi3*BrPhiEe*GPhi)*xf/dphi;
247 
248  G4double cross = 12.0*pi*MeVnb*norm(amp)/(e*s_inv);
249 
250  return cross;
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254 
255 G4double G4eeCrossSections::CrossSectionEtaG(G4double e)
256 {
257  complex<G4double> xf(cos(PhPhi),sin(PhPhi));
258 
259  G4double s_inv = e*e;
260  complex<G4double> drho = DpRho(e);
261  complex<G4double> dom = DpOm(e);
262  complex<G4double> dphi = DpPhi(e);
263 
264  complex<G4double> amp =
265  sqrt(WidthPg(s_inv,MsRho,GRho,BrRhoEtaG,MsEta)*MsRho3*BrRhoEe*GRho)/drho
266  + sqrt(WidthPg(s_inv,MsOm,GOm,BrOmEtaG,MsEta)*MsOm3*BrOmEe*GOm)/dom
267  + sqrt(WidthPg(s_inv,MsPhi,GPhi,BrPhiEtaG,MsEta)*MsPhi3*BrPhiEe*GPhi)*xf/dphi;
268 
269  G4double cross = 12.0*pi*MeVnb*norm(amp)/(e*s_inv);
270 
271  return cross;
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
275 
276 G4double G4eeCrossSections::CrossSection2Kcharged(G4double e)
277 {
278  G4double s_inv = e*e;
279  complex<G4double> dphi = DpPhi(e);
280 
281  complex<G4double> amp =
282  sqrt(Width2p(s_inv,MsPhi,GPhi,BrPhi2Kc,MsKc)*MsPhi3*BrPhiEe*GPhi)/dphi;
283 
284  G4double cross = 12.0*pi*MeVnb*norm(amp)/(e*s_inv);
285 
286  return cross;
287 }
288 
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
290 
291 G4double G4eeCrossSections::CrossSection2Kneutral(G4double e)
292 {
293  G4double s_inv = e*e;
294  complex<G4double> dphi = DpPhi(e);
295 
296  complex<G4double> amp =
297  sqrt(Width2p(s_inv,MsPhi,GPhi,BrPhiKsKl,MsKs)*MsPhi3*BrPhiEe*GPhi)/dphi;
298 
299  G4double cross = 12.0*pi*MeVnb*norm(amp)/(e*s_inv);
300 
301  return cross;
302 }
303 
304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
305 
306 G4double G4eeCrossSections::Width2p(G4double s_inv, G4double mres,
307  G4double gconst, G4double br, G4double mp)
308 {
309  G4double mp2 = 4.0*mp*mp;
310  G4double s0 = mres*mres;
311  G4double f = (s_inv - mp2)/(s0 - mp2);
312  if(f < 0.0) f = 0.0;
313  return gconst*br*sqrt(f)*f*s0/s_inv;
314 }
315 
316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
317 
318 G4double G4eeCrossSections::Width3p(G4double s_inv, G4double mres,
319  G4double gconst, G4double br)
320 {
321  G4double w = PhaseSpace3p(sqrt(s_inv));
322  G4double w0= PhaseSpace3p(mres);
323  G4double x = gconst*br*w/w0;
324  return x;
325 }
326 
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
328 
329 G4double G4eeCrossSections::PhaseSpace3p(G4double e)
330 {
331  // E.A.Kuraev, Z.K.Silagadze.
332  // Once more about the omega->3 pi contact term.
333  // Yadernaya Phisica, 1995, V58, N9, p.1678-1694.
334 
335  // G4bool b;
336  // G4double x = ph3p->GetValue(e, b);
337  G4double x = 1.0;
338  G4double emev = e/MeV;
339  G4double y = 414.12/emev;
340  x *= pow(e/MsOm, 5.0) * pow(emev*0.1, 3.0)*(1.0 - y*y);
341  return x;
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 
346 G4double G4eeCrossSections::WidthPg(G4double s_inv, G4double mres,
347  G4double gconst, G4double br, G4double mp)
348 {
349  G4double mp2 = mp*mp;
350  G4double s0 = mres*mres;
351  G4double f = (s_inv - mp2)*mres/((s0 - mp2)*sqrt(s_inv));
352  if(f < 0.0) f = 0.0;
353  return gconst*br*f*f*f;
354 }
355 
356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
357 
358 G4double G4eeCrossSections::WidthRho(G4double e)
359 {
360  G4double w = Width2p(e*e, MsRho, GRho, 1.0, MsPi);
361  return w;
362 }
363 
364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
365 
366 G4double G4eeCrossSections::WidthOm(G4double e)
367 {
368  G4double s_inv = e*e;
369  G4double w = (Width3p(s_inv, MsOm, GOm, BrOm3Pi) +
370  WidthPg(s_inv, MsOm, GOm, BrOmPi0G, MsPi0) +
371  WidthPg(s_inv, MsOm, GOm, BrOmEtaG, MsEta) +
372  Width2p(s_inv, MsOm, GOm, BrOm2Pi, MsPi)) /
373  (BrOm3Pi+BrOmPi0G+BrOmEtaG+BrOm2Pi);
374  return w;
375 }
376 
377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
378 
379 G4double G4eeCrossSections::WidthPhi(G4double e)
380 {
381  G4double s_inv = e*e;
382  G4double w = (Width3p(s_inv, MsPhi, GPhi, BrPhi3Pi) +
383  WidthPg(s_inv, MsPhi, GPhi, BrPhiPi0G, MsPi0) +
384  WidthPg(s_inv, MsPhi, GPhi, BrPhiEtaG, MsEta) +
385  Width2p(s_inv, MsPhi, GPhi, BrPhi2Kc, MsKc) +
386  Width2p(s_inv, MsPhi, GPhi, BrPhiKsKl, MsKs)) /
387  (BrPhi3Pi+BrPhiPi0G+BrPhiEtaG+BrPhi2Kc+BrPhiKsKl);
388  return w;
389 }
390 
391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
392 
393 complex<G4double> G4eeCrossSections::DpRho(G4double e)
394 {
395  complex<G4double> d(MsRho*MsRho - e*e, -e*WidthRho(e));
396  return d;
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 
401 complex<G4double> G4eeCrossSections::DpOm(G4double e)
402 {
403  complex<G4double> d(MsOm*MsOm - e*e, -e*WidthOm(e));
404  return d;
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408 
409 complex<G4double> G4eeCrossSections::DpPhi(G4double e)
410 {
411  complex<G4double> d(MsPhi*MsPhi - e*e, -e*WidthPhi(e));
412  return d;
413 }
414 
415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....