Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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$
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
66using namespace std;
67
68G4eeCrossSections::G4eeCrossSections()
69{
70 Initialise();
71}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74
75G4eeCrossSections::~G4eeCrossSections()
76{}
77
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79
80void 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
193G4double 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
216G4double 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
235G4double 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
256G4double 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
277G4double 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
292G4double 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
307G4double 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
319G4double 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
330G4double 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
347G4double 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
359G4double 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
367G4double 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
380G4double 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
394complex<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
402complex<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
410complex<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....
G4double
double G4double
Definition: G4Types.hh:64
G4Eta::Eta
static G4Eta * Eta()
Definition: G4Eta.cc:109
G4KaonPlus::KaonPlus
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4KaonZeroLong::KaonZeroLong
static G4KaonZeroLong * KaonZeroLong()
Definition: G4KaonZeroLong.cc:115
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4eeCrossSections::CrossSection2pi
G4double CrossSection2pi(G4double)
Definition: G4eeCrossSections.cc:193
G4eeCrossSections::~G4eeCrossSections
virtual ~G4eeCrossSections()
Definition: G4eeCrossSections.cc:75
G4eeCrossSections::CrossSection2Kcharged
G4double CrossSection2Kcharged(G4double)
Definition: G4eeCrossSections.cc:277
G4eeCrossSections::CrossSection2Kneutral
G4double CrossSection2Kneutral(G4double)
Definition: G4eeCrossSections.cc:292
G4eeCrossSections::CrossSectionEtaG
G4double CrossSectionEtaG(G4double)
Definition: G4eeCrossSections.cc:256
G4eeCrossSections::CrossSection3pi
G4double CrossSection3pi(G4double)
Definition: G4eeCrossSections.cc:216
G4eeCrossSections::CrossSectionPi0G
G4double CrossSectionPi0G(G4double)
Definition: G4eeCrossSections.cc:235
G4eeCrossSections::DpRho
std::complex< G4double > DpRho(G4double e)
Definition: G4eeCrossSections.cc:394