Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/im_r_matrix/src/G4XPDGElastic.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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] = std::move(ppData);
130   xMap[piPlusp] = std::move(pPiPlusData);
131   xMap[piMinusp] = std::move(pPiMinusData);
132   xMap[KPlusp] = std::move(pKPlusData);
133   xMap[KMinusp] = std::move(pKMinusData);
134   xMap[ppbar] = std::move(ppbarData);
135   xMap[npbar] = std::move(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 
277