source: trunk/source/processes/hadronic/models/im_r_matrix/src/G4XPDGTotal.cc@ 1007

Last change on this file since 1007 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 8.9 KB
Line 
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// PDG fits, 1998 Review of Particle Properties, European Phys. J. 3(1998), 1
28// -------------------------------------------------------------------
29
30#include "globals.hh"
31#include "G4ios.hh"
32#include "G4XPDGTotal.hh"
33#include "G4KineticTrack.hh"
34#include "G4ParticleDefinition.hh"
35#include "G4DataVector.hh"
36#include "G4AntiProton.hh"
37#include "G4AntiNeutron.hh"
38#include "G4Proton.hh"
39#include "G4Neutron.hh"
40#include "G4PionPlus.hh"
41#include "G4PionMinus.hh"
42#include "G4Gamma.hh"
43#include "G4KaonMinus.hh"
44#include "G4KaonPlus.hh"
45
46const G4double G4XPDGTotal::_lowLimit = 3. * GeV; // 2.99999 * GeV;
47const G4double G4XPDGTotal::_highLimit = DBL_MAX;
48
49// Parameters of the PDG total cross-section fit (Rev. Particle Properties, 1998)
50// Columns are: lower and higher fit range, X, Y1, Y2
51const G4int G4XPDGTotal::nFit = 5;
52// p p
53const G4double G4XPDGTotal::ppPDGFit[5] = { 3., 40000., 18.256, 60.19, 33.43 };
54// n p
55const G4double G4XPDGTotal::npPDGFit[5] = { 3., 40., 18.256, 61.14, 29.80 };
56// pi p
57const G4double G4XPDGTotal::pipPDGFit[5] = { 3., 40., 11.568, 27.55, 5.62 };
58// K p
59const G4double G4XPDGTotal::KpPDGFit[5] = { 3., 40., 10.376, 15.57, 13.19 };
60// K n
61const G4double G4XPDGTotal::KnPDGFit[5] = { 3., 40., 10.376, 14.29, 7.38 };
62// gamma p
63const G4double G4XPDGTotal::gammapPDGFit[5] = { 3., 300., 0.0577, 0.1171, 0. };
64//gamma gamma
65const G4double G4XPDGTotal::gammagammaPDGFit[5] = { 3., 300., 0.000156, 0.00032, 0. };
66
67
68G4XPDGTotal::G4XPDGTotal()
69{
70 std::pair<G4String,G4String> pp(G4Proton::ProtonDefinition()->GetParticleName(),
71 G4Proton::ProtonDefinition()->GetParticleName());
72 std::pair<G4String,G4String> pn(G4Proton::ProtonDefinition()->GetParticleName(),
73 G4Neutron::NeutronDefinition()->GetParticleName());
74 std::pair<G4String,G4String> piPlusp(G4PionPlus::PionPlusDefinition()->GetParticleName(),
75 G4Proton::ProtonDefinition()->GetParticleName());
76 std::pair<G4String,G4String> piMinusp(G4PionMinus::PionMinusDefinition()->GetParticleName(),
77 G4Proton::ProtonDefinition()->GetParticleName());
78 std::pair<G4String,G4String> KPlusp(G4KaonPlus::KaonPlusDefinition()->GetParticleName(),
79 G4Proton::ProtonDefinition()->GetParticleName());
80 std::pair<G4String,G4String> KPlusn(G4KaonPlus::KaonPlusDefinition()->GetParticleName(),
81 G4Neutron::NeutronDefinition()->GetParticleName());
82 std::pair<G4String,G4String> KMinusp(G4KaonMinus::KaonMinusDefinition()->GetParticleName(),
83 G4Proton::ProtonDefinition()->GetParticleName());
84 std::pair<G4String,G4String> KMinusn(G4KaonMinus::KaonMinusDefinition()->GetParticleName(),
85 G4Neutron::NeutronDefinition()->GetParticleName());
86 std::pair<G4String,G4String> gp(G4Gamma::GammaDefinition()->GetParticleName(),
87 G4Proton::ProtonDefinition()->GetParticleName());
88 std::pair<G4String,G4String> gg(G4Gamma::GammaDefinition()->GetParticleName(),
89 G4Gamma::GammaDefinition()->GetParticleName());
90 std::pair<G4String,G4String> nn(G4Neutron::NeutronDefinition()->GetParticleName(),
91 G4Neutron::NeutronDefinition()->GetParticleName());
92
93 std::vector<G4double> nnData;
94 std::vector<G4double> ppData;
95 std::vector<G4double> pnData;
96 std::vector<G4double> pipData;
97 std::vector<G4double> KpData;
98 std::vector<G4double> KnData;
99 std::vector<G4double> gpData;
100 std::vector<G4double> ggData;
101
102 G4int i;
103 for (i=0; i<2; i++)
104 {
105 nnData.push_back(ppPDGFit[i] * GeV);
106 ppData.push_back(ppPDGFit[i] * GeV);
107 pnData.push_back(npPDGFit[i] * GeV);
108 pipData.push_back(pipPDGFit[i] * GeV);
109 KpData.push_back(KpPDGFit[i] * GeV);
110 KnData.push_back(KnPDGFit[i] * GeV);
111 gpData.push_back(gammapPDGFit[i] * GeV);
112 ggData.push_back(gammagammaPDGFit[i] * GeV);
113 }
114 for (i=2; i<nFit; i++)
115 {
116 nnData.push_back(ppPDGFit[i]);
117 ppData.push_back(ppPDGFit[i]);
118 pnData.push_back(npPDGFit[i]);
119 pipData.push_back(pipPDGFit[i]);
120 KpData.push_back(KpPDGFit[i]);
121 KnData.push_back(KnPDGFit[i]);
122 gpData.push_back(gammapPDGFit[i]);
123 ggData.push_back(gammagammaPDGFit[i]);
124 }
125
126 xMap[pp] = ppData;
127 xMap[pn] = pnData;
128 xMap[piPlusp] = pipData;
129 xMap[piMinusp] = pipData;
130 xMap[KPlusp] = KpData;
131 xMap[KPlusn] = KnData;
132 xMap[KMinusp] = KpData;
133 xMap[KMinusn] = KnData;
134 xMap[gp] = gpData;
135 xMap[gg] = ggData;
136 xMap[nn] = nnData;
137}
138
139
140G4XPDGTotal::~G4XPDGTotal()
141{ }
142
143
144G4bool G4XPDGTotal::operator==(const G4XPDGTotal &right) const
145{
146 return (this == (G4XPDGTotal *) &right);
147}
148
149
150G4bool G4XPDGTotal::operator!=(const G4XPDGTotal &right) const
151{
152 return (this != (G4XPDGTotal *) &right);
153}
154
155
156G4double G4XPDGTotal::CrossSection(const G4KineticTrack& trk1,
157 const G4KineticTrack& trk2) const
158{
159 G4double sigma = 0.;
160
161 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
162
163 G4ParticleDefinition* def1 = trk1.GetDefinition();
164 G4ParticleDefinition* def2 = trk2.GetDefinition();
165
166 G4double enc1 = def1->GetPDGEncoding();
167 G4double enc2 = def2->GetPDGEncoding();
168 G4double coeff = -1.;
169 if ( (enc1 < 0 && enc2 >0) || (enc2 < 0 && enc1 >0) ) coeff = 1.;
170
171 // Order the pair: first is the lower mass particle, second is the higher mass one
172 std::pair<G4String,G4String> trkPair(def1->GetParticleName(),def2->GetParticleName());
173
174 if (def1->GetPDGMass() > def2->GetPDGMass())
175 trkPair = std::pair<G4String,G4String>(def2->GetParticleName(),def1->GetParticleName());
176
177 std::vector<G4double> data;
178
179 if (xMap.find(trkPair) != xMap.end())
180 {
181
182 PairDoubleMap::const_iterator iter;
183 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
184 {
185 std::pair<G4String,G4String> thePair = (*iter).first;
186 if (thePair == trkPair)
187 {
188 data = (*iter).second;
189
190 G4double eMinFit = data[0];
191 G4double eMaxFit = data[1];
192 G4double xFit = data[2];
193 G4double y1Fit = data[3];
194 G4double y2Fit = data[4];
195
196 // Total Cross-section fit, 1998 Review of Particle Properties, European Phys. J. 3(1998), 1
197
198 // Parameters from the PDG fit
199 const G4double epsilon = 0.095;
200 const G4double eta1 = -0.34;
201 const G4double eta2 = -0.55;
202
203 if (sqrtS < eMinFit || sqrtS > eMaxFit)
204 {
205 G4cout << "WARNING! G4XPDGTotal::PDGTotal extrapolating cross section at "
206 << sqrtS / GeV
207 << " GeV outside the PDG fit range "
208 << eMinFit / GeV << " - " << eMaxFit / GeV << " GeV " << G4endl;
209 }
210
211 G4double s = (sqrtS * sqrtS) / (GeV*GeV);
212
213 sigma = ( (xFit * std::pow(s,epsilon)) +
214 (y1Fit * std::pow(s,eta1)) +
215 (coeff * y2Fit * std::pow(s,eta2)) ) * millibarn;
216
217 if (sigma < 0.)
218 {
219 G4String name1 = def1->GetParticleName();
220 G4String name2 = def2->GetParticleName();
221 G4cout << "WARNING! G4XPDGTotal::PDGTotal "
222 << name1 << "-" << name2
223 << " total cross section: Ecm "
224 << sqrtS / GeV << " GeV, negative cross section "
225 << sigma / millibarn << " mb set to 0" << G4endl;
226 sigma = 0.;
227 }
228 }
229 }
230 }
231 return sigma;
232}
233
234
235G4String G4XPDGTotal::Name() const
236{
237 G4String name = "PDGTotal ";
238 return name;
239}
240
241
242G4bool G4XPDGTotal::IsValid(G4double e) const
243{
244 G4bool answer = InLimits(e,_lowLimit,_highLimit);
245
246 return answer;
247}
248
249
250G4double G4XPDGTotal::PDGTotal(G4double ,G4double ) const
251{
252 return 0.;
253}
254
255
256
257
258
259
260
261
262
263
264
Note: See TracBrowser for help on using the repository browser.