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