Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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