source: trunk/source/processes/hadronic/models/im_r_matrix/src/G4XPDGElastic.cc@ 962

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

import all except CVS
File size: 9.1 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// $Id: G4XPDGElastic.cc,v 1.3 2006/06/29 20:42:42 gunter Exp $ //
28// -------------------------------------------------------------------
29//
30// PDG Elastic cross section
31// PDG fits, Phys.Rev. D50 (1994), 1335
32//
33//
34// -------------------------------------------------------------------
35
36#include "globals.hh"
37#include "G4ios.hh"
38#include "G4XPDGElastic.hh"
39#include "G4KineticTrack.hh"
40#include "G4ParticleDefinition.hh"
41#include "G4DataVector.hh"
42
43#include "G4AntiProton.hh"
44#include "G4AntiNeutron.hh"
45#include "G4Proton.hh"
46#include "G4Neutron.hh"
47#include "G4PionPlus.hh"
48#include "G4PionMinus.hh"
49#include "G4KaonMinus.hh"
50#include "G4KaonPlus.hh"
51
52const G4double G4XPDGElastic::_lowLimit = 5. * GeV;
53const G4double G4XPDGElastic::_highLimit = DBL_MAX;
54
55// Parameters of the PDG Elastic cross-section fit (Rev. Particle Properties, 1998)
56// Columns are: lower and higher fit range, X, Y1, Y2
57const G4int G4XPDGElastic::nPar = 7;
58// p pi+
59const G4double G4XPDGElastic::pPiPlusPDGFit[7] = { 2., 200., 0., 11.4, -0.4, 0.079, 0. };
60// p pi-
61const G4double G4XPDGElastic::pPiMinusPDGFit[7] = { 2., 360., 1.76, 11.2, -0.64, 0.043, 0. };
62// p K+
63const G4double G4XPDGElastic::pKPlusPDGFit[7] = { 2., 175., 5.0, 8.1, -1.8, 0.16, -1.3 };
64// p K-
65const G4double G4XPDGElastic::pKMinusPDGFit[7] = { 2., 175., 7.3, 0., 0., 0.29, -2.40 };
66// p p
67const G4double G4XPDGElastic::ppPDGFit[7] = { 2., 2100., 11.9, 26.9, -1.21, 0.169, -1.85 };
68// p pbar
69const G4double G4XPDGElastic::ppbarPDGFit[7] = { 5., 1730000., 10.2, 52.7, -1.16, 0.125, -1.28 };
70// n pbar
71const G4double G4XPDGElastic::npbarPDGFit[7] = { 1.1, 5.55, 36.5, 0., 0., 0., -11.9 };
72
73
74G4XPDGElastic::G4XPDGElastic()
75{
76 G4String proton = G4Proton::ProtonDefinition()->GetParticleName();
77 G4String neutron = G4Neutron::NeutronDefinition()->GetParticleName();
78 G4String piPlus = G4PionPlus::PionPlusDefinition()->GetParticleName();
79 G4String piMinus = G4PionMinus::PionMinusDefinition()->GetParticleName();
80 G4String KPlus = G4KaonPlus::KaonPlusDefinition()->GetParticleName();
81 G4String KMinus = G4KaonMinus::KaonMinusDefinition()->GetParticleName();
82 G4String antiproton = G4AntiProton::AntiProtonDefinition()->GetParticleName();
83
84 std::pair<G4String,G4String> pp(proton,proton);
85 std::pair<G4String,G4String> pn(proton,neutron);
86 std::pair<G4String,G4String> piPlusp(piPlus,proton);
87 std::pair<G4String,G4String> piMinusp(piMinus,proton);
88 std::pair<G4String,G4String> KPlusp(KPlus,proton);
89 std::pair<G4String,G4String> KMinusp(KMinus,proton);
90 std::pair<G4String,G4String> nn(neutron,neutron);
91 std::pair<G4String,G4String> ppbar(proton,antiproton);
92 std::pair<G4String,G4String> npbar(antiproton,neutron);
93
94 std::vector<G4double> ppData;
95 std::vector<G4double> pPiPlusData;
96 std::vector<G4double> pPiMinusData;
97 std::vector<G4double> pKPlusData;
98 std::vector<G4double> pKMinusData;
99 std::vector<G4double> ppbarData;
100 std::vector<G4double> npbarData;
101
102 G4int i;
103 for (i=0; i<2; i++)
104 {
105 ppData.push_back(ppPDGFit[i] * GeV);
106 pPiPlusData.push_back(pPiPlusPDGFit[i] * GeV);
107 pPiMinusData.push_back(pPiMinusPDGFit[i] * GeV);
108 pKPlusData.push_back(pKPlusPDGFit[i] * GeV);
109 pKMinusData.push_back(pKMinusPDGFit[i] * GeV);
110 ppbarData.push_back(ppbarPDGFit[i] * GeV);
111 npbarData.push_back(npbarPDGFit[i] * GeV);
112 }
113
114 for (i=2; i<nPar; i++)
115 {
116 ppData.push_back(ppPDGFit[i]);
117 pPiPlusData.push_back(pPiPlusPDGFit[i]);
118 pPiMinusData.push_back(pPiMinusPDGFit[i]);
119 pKPlusData.push_back(pKPlusPDGFit[i]);
120 pKMinusData.push_back(pKMinusPDGFit[i]);
121 ppbarData.push_back(ppbarPDGFit[i]);
122 npbarData.push_back(npbarPDGFit[i]);
123 }
124
125 xMap[nn] = ppData;
126 xMap[pp] = ppData;
127 xMap[pn] = ppData;
128 xMap[piPlusp] = pPiPlusData;
129 xMap[piMinusp] = pPiMinusData;
130 xMap[KPlusp] = pKPlusData;
131 xMap[KMinusp] = pKMinusData;
132 xMap[ppbar] = ppbarData;
133 xMap[npbar] = npbarData;
134}
135
136
137G4XPDGElastic::~G4XPDGElastic()
138{ }
139
140
141G4bool G4XPDGElastic::operator==(const G4XPDGElastic &right) const
142{
143 return (this == (G4XPDGElastic *) &right);
144}
145
146
147G4bool G4XPDGElastic::operator!=(const G4XPDGElastic &right) const
148{
149 return (this != (G4XPDGElastic *) &right);
150}
151
152
153G4double G4XPDGElastic::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
154{
155 // Elastic Cross-section fit, 1994 Review of Particle Properties, (1994), 1
156
157 G4double sigma = 0.;
158
159 G4ParticleDefinition* def1 = trk1.GetDefinition();
160 G4ParticleDefinition* def2 = trk2.GetDefinition();
161 G4String name1 = def1->GetParticleName();
162 G4String name2 = def2->GetParticleName();
163
164 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
165 G4double m1 = def1->GetPDGMass();
166 G4double m2 = def2->GetPDGMass();
167 G4double m = std::max(m1,m2);
168 // if (m1 > m) m = m1;
169
170 G4double pLab = 0.;
171
172 if (m > 0. && sqrtS > (m1 + m2))
173 {
174 pLab = std::sqrt( (sqrtS*sqrtS - (m1+m2)*(m1+m2) ) * (sqrtS*sqrtS - (m1-m2)*(m1-m2)) ) / (2*m);
175
176 // The PDG fit formula requires p in GeV/c
177
178 G4double enc1 = def1->GetPDGEncoding();
179 G4double enc2 = def2->GetPDGEncoding();
180 G4double coeff = -1.;
181 if ( (enc1 < 0 && enc2 >0) || (enc2 < 0 && enc1 >0) ) coeff = 1.;
182
183 // Order the pair: first is the lower mass particle, second is the higher mass one
184 std::pair<G4String,G4String> trkPair(def1->GetParticleName(),def2->GetParticleName());
185 if (def1->GetPDGMass() > def2->GetPDGMass())
186 trkPair = std::pair<G4String,G4String>(def2->GetParticleName(),def1->GetParticleName());
187
188 std::vector<G4double> data;
189 G4double pMinFit = 0.;
190 G4double pMaxFit = 0.;
191 G4double aFit = 0.;
192 G4double bFit = 0.;
193 G4double cFit = 0.;
194 G4double dFit = 0.;
195 G4double nFit = 0.;
196
197 // Debug
198// G4cout << "Map has " << xMap.size() << " elements" << G4endl;
199 // Debug end
200
201 if (xMap.find(trkPair) != xMap.end())
202 {
203 PairDoubleMap::const_iterator iter;
204 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
205 {
206 std::pair<G4String,G4String> thePair = (*iter).first;
207 if (thePair == trkPair)
208 {
209 data = (*iter).second;
210 pMinFit = data[0];
211 pMaxFit = data[1];
212 aFit = data[2];
213 bFit = data[3];
214 cFit = data[5];
215 dFit = data[6];
216 nFit = data[4];
217
218 if (pLab < pMinFit) return 0.0;
219 if (pLab > pMaxFit )
220 G4cout << "WARNING! G4XPDGElastic::PDGElastic "
221 << trk1.GetDefinition()->GetParticleName() << "-"
222 << trk2.GetDefinition()->GetParticleName()
223 << " elastic cross section: momentum "
224 << pLab / GeV << " GeV outside valid fit range "
225 << pMinFit /GeV << "-" << pMaxFit / GeV
226 << G4endl;
227
228 pLab /= GeV;
229 if (pLab > 0.)
230 {
231 G4double logP = std::log(pLab);
232 sigma = aFit + bFit * std::pow(pLab, nFit) + cFit * logP*logP + dFit * logP;
233 sigma = sigma * millibarn;
234 }
235 }
236 }
237 }
238 else
239 {
240 G4cout << "G4XPDGElastic::CrossSection - Data not found in Map" << G4endl;
241 }
242 }
243
244 if (sigma < 0.)
245 {
246 G4cout << "WARNING! G4XPDGElastic::PDGElastic "
247 << name1 << "-" << name2
248 << " elastic cross section: momentum "
249 << pLab << " GeV, negative cross section "
250 << sigma / millibarn << " mb set to 0" << G4endl;
251 sigma = 0.;
252 }
253
254 return sigma;
255}
256
257
258G4String G4XPDGElastic::Name() const
259{
260 G4String name = "PDGElastic ";
261 return name;
262}
263
264
265G4bool G4XPDGElastic::IsValid(G4double e) const
266{
267 G4bool answer = InLimits(e,_lowLimit,_highLimit);
268
269 return answer;
270}
271
272
273
274
275
276
277
278
279
280
281
Note: See TracBrowser for help on using the repository browser.