Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/im_r_matrix/src/G4XPDGTotal.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 // 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] = std::move(ppData);
129   xMap[pn] = std::move(pnData);
130   xMap[piPlusp] = pipData;
131   xMap[piMinusp] = std::move(pipData);
132   xMap[KPlusp] = KpData;
133   xMap[KPlusn] = KnData;
134   xMap[KMinusp] = std::move(KpData);
135   xMap[KMinusn] = std::move(KnData);
136   xMap[gp] = std::move(gpData);
137   xMap[gg] = std::move(ggData);
138   xMap[nn] = std::move(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 
267