Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4XPDGElastic.cc
Go to the documentation of this file.
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 // $Id: G4XPDGElastic.cc,v 1.4 2010-03-12 15:45:18 gunter Exp $ //
28 // -------------------------------------------------------------------
29 //
30 // PDG Elastic cross section
31 // PDG fits, Phys.Rev. D50 (1994), 1335
32 //
33 //
34 // -------------------------------------------------------------------
35 
36 #include "globals.hh"
37 #include "G4ios.hh"
38 #include "G4Log.hh"
39 #include "G4Pow.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4XPDGElastic.hh"
42 #include "G4KineticTrack.hh"
43 #include "G4ParticleDefinition.hh"
44 #include "G4DataVector.hh"
45 
46 #include "G4AntiProton.hh"
47 #include "G4AntiNeutron.hh"
48 #include "G4Proton.hh"
49 #include "G4Neutron.hh"
50 #include "G4PionPlus.hh"
51 #include "G4PionMinus.hh"
52 #include "G4KaonMinus.hh"
53 #include "G4KaonPlus.hh"
54 
55 const G4double G4XPDGElastic::_lowLimit = 5. * GeV;
56 const G4double G4XPDGElastic::_highLimit = DBL_MAX;
57 
58 // Parameters of the PDG Elastic cross-section fit (Rev. Particle Properties, 1998)
59 // Columns are: lower and higher fit range, X, Y1, Y2
60 const G4int G4XPDGElastic::nPar = 7;
61 // p pi+
62 const G4double G4XPDGElastic::pPiPlusPDGFit[7] = { 2., 200., 0., 11.4, -0.4, 0.079, 0. };
63 // p pi-
64 const G4double G4XPDGElastic::pPiMinusPDGFit[7] = { 2., 360., 1.76, 11.2, -0.64, 0.043, 0. };
65 // p K+
66 const G4double G4XPDGElastic::pKPlusPDGFit[7] = { 2., 175., 5.0, 8.1, -1.8, 0.16, -1.3 };
67 // p K-
68 const G4double G4XPDGElastic::pKMinusPDGFit[7] = { 2., 175., 7.3, 0., 0., 0.29, -2.40 };
69 // p p
70 const G4double G4XPDGElastic::ppPDGFit[7] = { 2., 2100., 11.9, 26.9, -1.21, 0.169, -1.85 };
71 // p pbar
72 const G4double G4XPDGElastic::ppbarPDGFit[7] = { 5., 1730000., 10.2, 52.7, -1.16, 0.125, -1.28 };
73 // n pbar
74 const G4double G4XPDGElastic::npbarPDGFit[7] = { 1.1, 5.55, 36.5, 0., 0., 0., -11.9 };
75 
76 
78 {
86 
87  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> pp(proton,proton);
88  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> pn(proton,neutron);
89  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> piPlusp(piPlus,proton);
90  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> piMinusp(piMinus,proton);
91  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> KPlusp(KPlus,proton);
92  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> KMinusp(KMinus,proton);
93  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> nn(neutron,neutron);
94  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> ppbar(proton,antiproton);
95  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> npbar(antiproton,neutron);
96 
97  std::vector<G4double> ppData;
98  std::vector<G4double> pPiPlusData;
99  std::vector<G4double> pPiMinusData;
100  std::vector<G4double> pKPlusData;
101  std::vector<G4double> pKMinusData;
102  std::vector<G4double> ppbarData;
103  std::vector<G4double> npbarData;
104 
105  G4int i;
106  for (i=0; i<2; i++)
107  {
108  ppData.push_back(ppPDGFit[i] * GeV);
109  pPiPlusData.push_back(pPiPlusPDGFit[i] * GeV);
110  pPiMinusData.push_back(pPiMinusPDGFit[i] * GeV);
111  pKPlusData.push_back(pKPlusPDGFit[i] * GeV);
112  pKMinusData.push_back(pKMinusPDGFit[i] * GeV);
113  ppbarData.push_back(ppbarPDGFit[i] * GeV);
114  npbarData.push_back(npbarPDGFit[i] * GeV);
115  }
116 
117  for (i=2; i<nPar; i++)
118  {
119  ppData.push_back(ppPDGFit[i]);
120  pPiPlusData.push_back(pPiPlusPDGFit[i]);
121  pPiMinusData.push_back(pPiMinusPDGFit[i]);
122  pKPlusData.push_back(pKPlusPDGFit[i]);
123  pKMinusData.push_back(pKMinusPDGFit[i]);
124  ppbarData.push_back(ppbarPDGFit[i]);
125  npbarData.push_back(npbarPDGFit[i]);
126  }
127 
128  xMap[nn] = ppData;
129  xMap[pp] = ppData;
130  xMap[pn] = ppData;
131  xMap[piPlusp] = pPiPlusData;
132  xMap[piMinusp] = pPiMinusData;
133  xMap[KPlusp] = pKPlusData;
134  xMap[KMinusp] = pKMinusData;
135  xMap[ppbar] = ppbarData;
136  xMap[npbar] = npbarData;
137 }
138 
139 
141 { }
142 
143 
145 {
146  return (this == (G4XPDGElastic *) &right);
147 }
148 
149 
151 {
152  return (this != (G4XPDGElastic *) &right);
153 }
154 
155 
157 {
158  // Elastic Cross-section fit, 1994 Review of Particle Properties, (1994), 1
159 
160  G4double sigma = 0.;
161 
162  const G4ParticleDefinition* def1 = trk1.GetDefinition();
163  const G4ParticleDefinition* def2 = trk2.GetDefinition();
164 
165  G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
166  G4double m_1 = def1->GetPDGMass();
167  G4double m_2 = def2->GetPDGMass();
168  G4double m_max = std::max(m_1,m_2);
169  // if (m1 > m) m = m1;
170 
171  G4double pLab = 0.;
172 
173  if (m_max > 0. && sqrtS > (m_1 + m_2))
174  {
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);
176 
177  // The PDG fit formula requires p in GeV/c
178 
179  // Order the pair: first is the lower mass particle, second is the higher mass one
180  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> trkPair(def1,def2);
181  if (def1->GetPDGMass() > def2->GetPDGMass())
182  trkPair = std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *>(def2,def1);
183 
184  std::vector<G4double> data;
185  G4double pMinFit = 0.;
186  G4double pMaxFit = 0.;
187  G4double aFit = 0.;
188  G4double bFit = 0.;
189  G4double cFit = 0.;
190  G4double dFit = 0.;
191  G4double nFit = 0.;
192 
193  // Debug
194 // G4cout << "Map has " << xMap.size() << " elements" << G4endl;
195  // Debug end
196 
197  if (xMap.find(trkPair) != xMap.end())
198  {
199  PairDoubleMap::const_iterator iter;
200  for (iter = xMap.begin(); iter != xMap.end(); ++iter)
201  {
202  std::pair<const G4ParticleDefinition *,const G4ParticleDefinition *> thePair = (*iter).first;
203  if (thePair == trkPair)
204  {
205  data = (*iter).second;
206  pMinFit = data[0];
207  pMaxFit = data[1];
208  aFit = data[2];
209  bFit = data[3];
210  cFit = data[5];
211  dFit = data[6];
212  nFit = data[4];
213 
214  if (pLab < pMinFit) return 0.0;
215  if (pLab > pMaxFit )
216  G4cout << "WARNING! G4XPDGElastic::PDGElastic "
217  << trk1.GetDefinition()->GetParticleName() << "-"
218  << trk2.GetDefinition()->GetParticleName()
219  << " elastic cross section: momentum "
220  << pLab / GeV << " GeV outside valid fit range "
221  << pMinFit /GeV << "-" << pMaxFit / GeV
222  << G4endl;
223 
224  pLab /= GeV;
225  if (pLab > 0.)
226  {
227  G4double logP = G4Log(pLab);
228  sigma = aFit + bFit * G4Pow::GetInstance()->powA(pLab, nFit) + cFit * logP*logP + dFit * logP;
229  sigma = sigma * millibarn;
230  }
231  }
232  }
233  }
234  else
235  {
236  G4cout << "G4XPDGElastic::CrossSection - Data not found in Map" << G4endl;
237  }
238  }
239 
240  if (sigma < 0.)
241  {
242  G4cout << "WARNING! G4XPDGElastic::PDGElastic "
243  << def1->GetParticleName() << "-" << def2->GetParticleName()
244  << " elastic cross section: momentum "
245  << pLab << " GeV, negative cross section "
246  << sigma / millibarn << " mb set to 0" << G4endl;
247  sigma = 0.;
248  }
249 
250  return sigma;
251 }
252 
253 
255 {
256  G4String name = "PDGElastic ";
257  return name;
258 }
259 
260 
262 {
263  G4bool answer = InLimits(e,_lowLimit,_highLimit);
264 
265  return answer;
266 }
267 
268 
269 
270 
271 
272 
273 
274 
275 
276 
277 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
const XML_Char * name
Definition: expat.h:151
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
static G4KaonPlus * KaonPlusDefinition()
Definition: G4KaonPlus.cc:108
virtual G4String Name() const
static G4KaonMinus * KaonMinusDefinition()
Definition: G4KaonMinus.cc:108
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
static G4AntiProton * AntiProtonDefinition()
Definition: G4AntiProton.cc:88
int G4int
Definition: G4Types.hh:78
G4bool operator!=(const G4XPDGElastic &right) const
const G4String & GetParticleName() const
const XML_Char const XML_Char * data
Definition: expat.h:268
G4GLOB_DLL std::ostream G4cout
static G4PionMinus * PionMinusDefinition()
Definition: G4PionMinus.cc:93
bool G4bool
Definition: G4Types.hh:79
static G4PionPlus * PionPlusDefinition()
Definition: G4PionPlus.cc:93
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
G4bool operator==(const G4XPDGElastic &right) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
virtual G4double CrossSection(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
const G4LorentzVector & Get4Momentum() const
#define DBL_MAX
Definition: templates.hh:83
virtual G4bool IsValid(G4double e) const
const G4ParticleDefinition * GetDefinition() const
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
virtual ~G4XPDGElastic()
static constexpr double millibarn
Definition: G4SIunits.hh:106
G4bool InLimits(G4double e, G4double eLow, G4double eHigh) const