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