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