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

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 9.3 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.4 2010/03/12 15:45:18 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 G4ParticleDefinition * proton = G4Proton::ProtonDefinition();
77 G4ParticleDefinition * neutron = G4Neutron::NeutronDefinition();
78 G4ParticleDefinition * piPlus = G4PionPlus::PionPlusDefinition();
79 G4ParticleDefinition * piMinus = G4PionMinus::PionMinusDefinition();
80 G4ParticleDefinition * KPlus = G4KaonPlus::KaonPlusDefinition();
81 G4ParticleDefinition * KMinus = G4KaonMinus::KaonMinusDefinition();
82 G4ParticleDefinition * antiproton = G4AntiProton::AntiProtonDefinition();
83
84 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> pp(proton,proton);
85 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> pn(proton,neutron);
86 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> piPlusp(piPlus,proton);
87 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> piMinusp(piMinus,proton);
88 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> KPlusp(KPlus,proton);
89 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> KMinusp(KMinus,proton);
90 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> nn(neutron,neutron);
91 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> ppbar(proton,antiproton);
92 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> 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
162 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
163 G4double m1 = def1->GetPDGMass();
164 G4double m2 = def2->GetPDGMass();
165 G4double m = std::max(m1,m2);
166 // if (m1 > m) m = m1;
167
168 G4double pLab = 0.;
169
170 if (m > 0. && sqrtS > (m1 + m2))
171 {
172 pLab = std::sqrt( (sqrtS*sqrtS - (m1+m2)*(m1+m2) ) * (sqrtS*sqrtS - (m1-m2)*(m1-m2)) ) / (2*m);
173
174 // The PDG fit formula requires p in GeV/c
175
176 G4double enc1 = def1->GetPDGEncoding();
177 G4double enc2 = def2->GetPDGEncoding();
178 G4double coeff = -1.;
179 if ( (enc1 < 0 && enc2 >0) || (enc2 < 0 && enc1 >0) ) coeff = 1.;
180
181 // Order the pair: first is the lower mass particle, second is the higher mass one
182 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> trkPair(def1,def2);
183 if (def1->GetPDGMass() > def2->GetPDGMass())
184 trkPair = std::pair<G4ParticleDefinition *,G4ParticleDefinition *>(def2,def1);
185
186 std::vector<G4double> data;
187 G4double pMinFit = 0.;
188 G4double pMaxFit = 0.;
189 G4double aFit = 0.;
190 G4double bFit = 0.;
191 G4double cFit = 0.;
192 G4double dFit = 0.;
193 G4double nFit = 0.;
194
195 // Debug
196// G4cout << "Map has " << xMap.size() << " elements" << G4endl;
197 // Debug end
198
199 if (xMap.find(trkPair) != xMap.end())
200 {
201 PairDoubleMap::const_iterator iter;
202 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
203 {
204 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> thePair = (*iter).first;
205 if (thePair == trkPair)
206 {
207 data = (*iter).second;
208 pMinFit = data[0];
209 pMaxFit = data[1];
210 aFit = data[2];
211 bFit = data[3];
212 cFit = data[5];
213 dFit = data[6];
214 nFit = data[4];
215
216 if (pLab < pMinFit) return 0.0;
217 if (pLab > pMaxFit )
218 G4cout << "WARNING! G4XPDGElastic::PDGElastic "
219 << trk1.GetDefinition()->GetParticleName() << "-"
220 << trk2.GetDefinition()->GetParticleName()
221 << " elastic cross section: momentum "
222 << pLab / GeV << " GeV outside valid fit range "
223 << pMinFit /GeV << "-" << pMaxFit / GeV
224 << G4endl;
225
226 pLab /= GeV;
227 if (pLab > 0.)
228 {
229 G4double logP = std::log(pLab);
230 sigma = aFit + bFit * std::pow(pLab, nFit) + cFit * logP*logP + dFit * logP;
231 sigma = sigma * millibarn;
232 }
233 }
234 }
235 }
236 else
237 {
238 G4cout << "G4XPDGElastic::CrossSection - Data not found in Map" << G4endl;
239 }
240 }
241
242 if (sigma < 0.)
243 {
244 G4cout << "WARNING! G4XPDGElastic::PDGElastic "
245 << def1->GetParticleName() << "-" << def2->GetParticleName()
246 << " elastic cross section: momentum "
247 << pLab << " GeV, negative cross section "
248 << sigma / millibarn << " mb set to 0" << G4endl;
249 sigma = 0.;
250 }
251
252 return sigma;
253}
254
255
256G4String G4XPDGElastic::Name() const
257{
258 G4String name = "PDGElastic ";
259 return name;
260}
261
262
263G4bool G4XPDGElastic::IsValid(G4double e) const
264{
265 G4bool answer = InLimits(e,_lowLimit,_highLimit);
266
267 return answer;
268}
269
270
271
272
273
274
275
276
277
278
279
Note: See TracBrowser for help on using the repository browser.