| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // PDG Elastic cross section
30 // PDG fits, Phys.Rev. D50 (1994), 1335
31 //
32 //
33 // -------------------------------------------
34
35 #include "globals.hh"
36 #include "G4ios.hh"
37 #include "G4Log.hh"
38 #include "G4Pow.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4XPDGElastic.hh"
41 #include "G4KineticTrack.hh"
42 #include "G4ParticleDefinition.hh"
43 #include "G4DataVector.hh"
44
45 #include "G4AntiProton.hh"
46 #include "G4AntiNeutron.hh"
47 #include "G4Proton.hh"
48 #include "G4Neutron.hh"
49 #include "G4PionPlus.hh"
50 #include "G4PionMinus.hh"
51 #include "G4KaonMinus.hh"
52 #include "G4KaonPlus.hh"
53
54 const G4double G4XPDGElastic::_lowLimit = 5. *
55 const G4double G4XPDGElastic::_highLimit = DBL
56
57 // Parameters of the PDG Elastic cross-section
58 // Columns are: lower and higher fit range, X,
59 const G4int G4XPDGElastic::nPar = 7;
60 // p pi+
61 const G4double G4XPDGElastic::pPiPlusPDGFit[7]
62 // p pi-
63 const G4double G4XPDGElastic::pPiMinusPDGFit[7
64 // p K+
65 const G4double G4XPDGElastic::pKPlusPDGFit[7]
66 // p K-
67 const G4double G4XPDGElastic::pKMinusPDGFit[7]
68 // p p
69 const G4double G4XPDGElastic::ppPDGFit[7] =
70 // p pbar
71 const G4double G4XPDGElastic::ppbarPDGFit[7] =
72 // n pbar
73 const G4double G4XPDGElastic::npbarPDGFit[7] =
74
75
76 G4XPDGElastic::G4XPDGElastic()
77 {
78 const G4ParticleDefinition * proton = G4Prot
79 const G4ParticleDefinition * neutron = G4Neu
80 const G4ParticleDefinition * piPlus = G4Pion
81 const G4ParticleDefinition * piMinus = G4Pio
82 const G4ParticleDefinition * KPlus = G4KaonP
83 const G4ParticleDefinition * KMinus = G4Kaon
84 const G4ParticleDefinition * antiproton = G4
85
86 std::pair<const G4ParticleDefinition *,const
87 std::pair<const G4ParticleDefinition *,const
88 std::pair<const G4ParticleDefinition *,const
89 std::pair<const G4ParticleDefinition *,const
90 std::pair<const G4ParticleDefinition *,const
91 std::pair<const G4ParticleDefinition *,const
92 std::pair<const G4ParticleDefinition *,const
93 std::pair<const G4ParticleDefinition *,const
94 std::pair<const G4ParticleDefinition *,const
95
96 std::vector<G4double> ppData;
97 std::vector<G4double> pPiPlusData;
98 std::vector<G4double> pPiMinusData;
99 std::vector<G4double> pKPlusData;
100 std::vector<G4double> pKMinusData;
101 std::vector<G4double> ppbarData;
102 std::vector<G4double> npbarData;
103
104 G4int i;
105 for (i=0; i<2; i++)
106 {
107 ppData.push_back(ppPDGFit[i] * GeV);
108 pPiPlusData.push_back(pPiPlusPDGFit[i] *
109 pPiMinusData.push_back(pPiMinusPDGFit[i]
110 pKPlusData.push_back(pKPlusPDGFit[i] * G
111 pKMinusData.push_back(pKMinusPDGFit[i] *
112 ppbarData.push_back(ppbarPDGFit[i] * GeV
113 npbarData.push_back(npbarPDGFit[i] * GeV
114 }
115
116 for (i=2; i<nPar; i++)
117 {
118 ppData.push_back(ppPDGFit[i]);
119 pPiPlusData.push_back(pPiPlusPDGFit[i]);
120 pPiMinusData.push_back(pPiMinusPDGFit[i]
121 pKPlusData.push_back(pKPlusPDGFit[i]);
122 pKMinusData.push_back(pKMinusPDGFit[i]);
123 ppbarData.push_back(ppbarPDGFit[i]);
124 npbarData.push_back(npbarPDGFit[i]);
125 }
126
127 xMap[nn] = ppData;
128 xMap[pp] = ppData;
129 xMap[pn] = std::move(ppData);
130 xMap[piPlusp] = std::move(pPiPlusData);
131 xMap[piMinusp] = std::move(pPiMinusData);
132 xMap[KPlusp] = std::move(pKPlusData);
133 xMap[KMinusp] = std::move(pKMinusData);
134 xMap[ppbar] = std::move(ppbarData);
135 xMap[npbar] = std::move(npbarData);
136 }
137
138
139 G4XPDGElastic::~G4XPDGElastic()
140 { }
141
142
143 G4bool G4XPDGElastic::operator==(const G4XPDGE
144 {
145 return (this == (G4XPDGElastic *) &right);
146 }
147
148
149 G4bool G4XPDGElastic::operator!=(const G4XPDGE
150 {
151 return (this != (G4XPDGElastic *) &right);
152 }
153
154
155 G4double G4XPDGElastic::CrossSection(const G4K
156 {
157 // Elastic Cross-section fit, 1994 Review of
158
159 G4double sigma = 0.;
160
161 const G4ParticleDefinition* def1 = trk1.GetD
162 const G4ParticleDefinition* def2 = trk2.GetD
163
164 G4double sqrtS = (trk1.Get4Momentum() + trk2
165 G4double m_1 = def1->GetPDGMass();
166 G4double m_2 = def2->GetPDGMass();
167 G4double m_max = std::max(m_1,m_2);
168 // if (m1 > m) m = m1;
169
170 G4double pLab = 0.;
171
172 if (m_max > 0. && sqrtS > (m_1 + m_2))
173 {
174 pLab = std::sqrt( (sqrtS*sqrtS - (m_1+m_
175
176 // The PDG fit formula requires p in GeV
177
178 // Order the pair: first is the lower ma
179 std::pair<const G4ParticleDefinition *,c
180 if (def1->GetPDGMass() > def2->GetPDGMas
181 trkPair = std::pair<const G4ParticleDefiniti
182
183 std::vector<G4double> data;
184 G4double pMinFit = 0.;
185 G4double pMaxFit = 0.;
186 G4double aFit = 0.;
187 G4double bFit = 0.;
188 G4double cFit = 0.;
189 G4double dFit = 0.;
190 G4double nFit = 0.;
191
192 // Debug
193 // G4cout << "Map has " << xMap.size() <<
194 // Debug end
195
196 if (xMap.find(trkPair) != xMap.end())
197 {
198 PairDoubleMap::const_iterator iter;
199 for (iter = xMap.begin(); iter != xMap.end
200 {
201 std::pair<const G4ParticleDefinition *
202 if (thePair == trkPair)
203 {
204 data = (*iter).second;
205 pMinFit = data[0];
206 pMaxFit = data[1];
207 aFit = data[2];
208 bFit = data[3];
209 cFit = data[5];
210 dFit = data[6];
211 nFit = data[4];
212
213 if (pLab < pMinFit) return 0.0;
214 if (pLab > pMaxFit )
215 G4cout << "WARNING! G4XPDGElastic::PDG
216 << trk1.GetDefinition()->GetParticleN
217 << trk2.GetDefinition()->GetParticleN
218 << " elastic cross section: momentum
219 << pLab / GeV << " GeV outside valid
220 << pMinFit /GeV << "-" << pMaxFit / G
221 << G4endl;
222
223 pLab /= GeV;
224 if (pLab > 0.)
225 {
226 G4double logP = G4Log(pLab);
227 sigma = aFit + bFit * G4Pow::GetInst
228 sigma = sigma * millibarn;
229 }
230 }
231 }
232 }
233 else
234 {
235 G4cout << "G4XPDGElastic::CrossSection - D
236 }
237 }
238
239 if (sigma < 0.)
240 {
241 G4cout << "WARNING! G4XPDGElastic::PDGEl
242 << def1->GetParticleName() << "-" << de
243 << " elastic cross section: momentum "
244 << pLab << " GeV, negative cross sectio
245 << sigma / millibarn << " mb set to 0"
246 sigma = 0.;
247 }
248
249 return sigma;
250 }
251
252
253 G4String G4XPDGElastic::Name() const
254 {
255 G4String name = "PDGElastic ";
256 return name;
257 }
258
259
260 G4bool G4XPDGElastic::IsValid(G4double e) cons
261 {
262 G4bool answer = InLimits(e,_lowLimit,_highLi
263
264 return answer;
265 }
266
267
268
269
270
271
272
273
274
275
276
277