ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4XPDGElastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4XPDGElastic.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 // -------------------------------------------------------------------
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. * GeV;
55 const G4double G4XPDGElastic::_highLimit = DBL_MAX;
56 
57 // Parameters of the PDG Elastic cross-section fit (Rev. Particle Properties, 1998)
58 // Columns are: lower and higher fit range, X, Y1, Y2
59 const G4int G4XPDGElastic::nPar = 7;
60 // p pi+
61 const G4double G4XPDGElastic::pPiPlusPDGFit[7] = { 2., 200., 0., 11.4, -0.4, 0.079, 0. };
62 // p pi-
63 const G4double G4XPDGElastic::pPiMinusPDGFit[7] = { 2., 360., 1.76, 11.2, -0.64, 0.043, 0. };
64 // p K+
65 const G4double G4XPDGElastic::pKPlusPDGFit[7] = { 2., 175., 5.0, 8.1, -1.8, 0.16, -1.3 };
66 // p K-
67 const G4double G4XPDGElastic::pKMinusPDGFit[7] = { 2., 175., 7.3, 0., 0., 0.29, -2.40 };
68 // p p
69 const G4double G4XPDGElastic::ppPDGFit[7] = { 2., 2100., 11.9, 26.9, -1.21, 0.169, -1.85 };
70 // p pbar
71 const G4double G4XPDGElastic::ppbarPDGFit[7] = { 5., 1730000., 10.2, 52.7, -1.16, 0.125, -1.28 };
72 // n pbar
73 const G4double G4XPDGElastic::npbarPDGFit[7] = { 1.1, 5.55, 36.5, 0., 0., 0., -11.9 };
74 
75 
76 G4XPDGElastic::G4XPDGElastic()
77 {
78  const G4ParticleDefinition * proton = G4Proton::ProtonDefinition();
79  const G4ParticleDefinition * neutron = G4Neutron::NeutronDefinition();
80  const G4ParticleDefinition * piPlus = G4PionPlus::PionPlusDefinition();
81  const G4ParticleDefinition * piMinus = G4PionMinus::PionMinusDefinition();
82  const G4ParticleDefinition * KPlus = G4KaonPlus::KaonPlusDefinition();
83  const G4ParticleDefinition * KMinus = G4KaonMinus::KaonMinusDefinition();
84  const G4ParticleDefinition * antiproton = G4AntiProton::AntiProtonDefinition();
85 
86  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> pp(proton,proton);
87  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> pn(proton,neutron);
88  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> piPlusp(piPlus,proton);
89  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> piMinusp(piMinus,proton);
90  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> KPlusp(KPlus,proton);
91  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> KMinusp(KMinus,proton);
92  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> nn(neutron,neutron);
93  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> ppbar(proton,antiproton);
94  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> npbar(antiproton,neutron);
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] * GeV);
109  pPiMinusData.push_back(pPiMinusPDGFit[i] * GeV);
110  pKPlusData.push_back(pKPlusPDGFit[i] * GeV);
111  pKMinusData.push_back(pKMinusPDGFit[i] * GeV);
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] = ppData;
130  xMap[piPlusp] = pPiPlusData;
131  xMap[piMinusp] = pPiMinusData;
132  xMap[KPlusp] = pKPlusData;
133  xMap[KMinusp] = pKMinusData;
134  xMap[ppbar] = ppbarData;
135  xMap[npbar] = npbarData;
136 }
137 
138 
139 G4XPDGElastic::~G4XPDGElastic()
140 { }
141 
142 
143 G4bool G4XPDGElastic::operator==(const G4XPDGElastic &right) const
144 {
145  return (this == (G4XPDGElastic *) &right);
146 }
147 
148 
149 G4bool G4XPDGElastic::operator!=(const G4XPDGElastic &right) const
150 {
151  return (this != (G4XPDGElastic *) &right);
152 }
153 
154 
155 G4double G4XPDGElastic::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
156 {
157  // Elastic Cross-section fit, 1994 Review of Particle Properties, (1994), 1
158 
159  G4double sigma = 0.;
160 
161  const G4ParticleDefinition* def1 = trk1.GetDefinition();
162  const G4ParticleDefinition* def2 = trk2.GetDefinition();
163 
164  G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
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_2)*(m_1+m_2) ) * (sqrtS*sqrtS - (m_1-m_2)*(m_1-m_2)) ) / (2*m_max);
175 
176  // The PDG fit formula requires p in GeV/c
177 
178  // Order the pair: first is the lower mass particle, second is the higher mass one
179  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> trkPair(def1,def2);
180  if (def1->GetPDGMass() > def2->GetPDGMass())
181  trkPair = std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *>(def2,def1);
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() << " elements" << G4endl;
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(); ++iter)
200  {
201  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> thePair = (*iter).first;
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::PDGElastic "
216  << trk1.GetDefinition()->GetParticleName() << "-"
217  << trk2.GetDefinition()->GetParticleName()
218  << " elastic cross section: momentum "
219  << pLab / GeV << " GeV outside valid fit range "
220  << pMinFit /GeV << "-" << pMaxFit / GeV
221  << G4endl;
222 
223  pLab /= GeV;
224  if (pLab > 0.)
225  {
226  G4double logP = G4Log(pLab);
227  sigma = aFit + bFit * G4Pow::GetInstance()->powA(pLab, nFit) + cFit * logP*logP + dFit * logP;
228  sigma = sigma * millibarn;
229  }
230  }
231  }
232  }
233  else
234  {
235  G4cout << "G4XPDGElastic::CrossSection - Data not found in Map" << G4endl;
236  }
237  }
238 
239  if (sigma < 0.)
240  {
241  G4cout << "WARNING! G4XPDGElastic::PDGElastic "
242  << def1->GetParticleName() << "-" << def2->GetParticleName()
243  << " elastic cross section: momentum "
244  << pLab << " GeV, negative cross section "
245  << sigma / millibarn << " mb set to 0" << G4endl;
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) const
261 {
262  G4bool answer = InLimits(e,_lowLimit,_highLimit);
263 
264  return answer;
265 }
266 
267 
268 
269 
270 
271 
272 
273 
274 
275 
276