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