source: trunk/source/processes/hadronic/models/im_r_matrix/src/G4XPDGTotal.cc @ 1334

Last change on this file since 1334 was 1228, checked in by garnier, 15 years ago

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