Geant4_10
G4XNNElasticLowE.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 #include "globals.hh"
27 #include "G4ios.hh"
28 #include "G4SystemOfUnits.hh"
29 #include "G4XNNElasticLowE.hh"
30 #include "G4KineticTrack.hh"
31 #include "G4ParticleDefinition.hh"
32 #include "G4PhysicsVector.hh"
33 #include "G4PhysicsLnVector.hh"
34 #include "G4ParticleDefinition.hh"
35 #include "G4Proton.hh"
36 #include "G4Neutron.hh"
37 
38 const G4double G4XNNElasticLowE::_lowLimit = 0.;
39 const G4double G4XNNElasticLowE::_highLimit = 3.*GeV;
40 
41 // Low energy limit of the cross-section table (in GeV)
42 // Units are assigned while filling the PhysicsVector
43 const G4double G4XNNElasticLowE::_eMinTable = 1.8964808;
44 const G4double G4XNNElasticLowE::_eStepLog = 0.01;
45 
46 // Cross-sections in mb
47 // Units are assigned while filling the PhysicsVector
48 
49 const G4int G4XNNElasticLowE::tableSize = 101;
50 
51 const G4double G4XNNElasticLowE::ppTable[101] =
52 {
53  60.00, //was 0.
54  33.48, 26.76, 25.26, 24.55, 23.94, 23.77, 23.72, 23.98,
55  25.48, 27.52, 27.72, 27.21, 25.80, 26.00, 24.32, 23.81,
56  24.37, 24.36, 23.13, 22.43, 21.71, 21.01, 20.83, 20.74,
57  20.25, 20.10, 20.59, 20.04, 20.83, 20.84, 21.07, 20.83,
58  20.79, 21.88, 21.15, 20.92, 19.00, 18.60, 17.30, 17.00,
59  16.70, 16.50, 16.20, 15.80, 15.57, 15.20, 15.00, 14.60,
60  14.20, 14.00, 13.80, 13.60, 13.40, 13.20, 13.00, 12.85,
61  12.70, 12.60, 12.50, 12.40, 12.30, 12.20, 12.10, 12.00,
62  11.90, 11.80, 11.75, 11.70, 11.64, 11.53, 11.41, 11.31,
63  11.22, 11.13, 11.05, 10.97, 10.89, 10.82, 10.75, 10.68,
64  10.61, 10.54, 10.48, 10.41, 10.35, 10.28, 10.22, 10.16,
65  10.13, 10.10, 10.08, 10.05, 10.02, 9.99, 9.96, 9.93,
66  9.90, 9.87, 9.84, 9.80
67 };
68 
69 const G4double G4XNNElasticLowE::npTable[101] =
70 {
71  1500.00, // was 0.
72  248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 37.20, 35.98,
73  35.02, 34.47, 32.48, 30.76, 29.46, 28.53, 27.84, 27.20,
74  26.53, 25.95, 25.59, 25.46, 25.00, 24.49, 24.08, 23.86,
75  23.17, 22.70, 21.88, 21.48, 20.22, 19.75, 18.97, 18.39,
76  17.98, 17.63, 17.21, 16.72, 16.68, 16.58, 16.42, 16.22,
77  15.98, 15.71, 15.42, 15.14, 14.87, 14.65, 14.44, 14.26,
78  14.10, 13.95, 13.80, 13.64, 13.47, 13.29, 13.09, 12.89,
79  12.68, 12.47, 12.27, 12.06, 11.84, 11.76, 11.69, 11.60,
80  11.50, 11.41, 11.29, 11.17, 11.06, 10.93, 10.81, 10.68,
81  10.56, 10.44, 10.33, 10.21, 10.12, 10.03, 9.96, 9.89,
82  9.83, 9.80, 9.77, 9.75, 9.74, 9.74, 9.74, 9.76,
83  9.73, 9.70, 9.68, 9.65, 9.63, 9.60, 9.57, 9.55,
84  9.52, 9.49, 9.46, 9.43
85 };
86 
87 
89 {
90  // Cross-sections are available in the range (_eMin,_eMax)
91 
92  _eMin = _eMinTable * GeV;
93  _eMax = std::exp(std::log(_eMinTable) + tableSize * _eStepLog) * GeV;
94  if (_eMin < _lowLimit)
95  throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
96  if (_highLimit > _eMax)
97  throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - High energy limit not valid");
98  G4PhysicsVector* pp = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
99 
100  _eMin = std::exp(std::log(_eMinTable)-_eStepLog)*GeV;
101  if (_eMin < _lowLimit)
102  throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
103  G4PhysicsVector* np = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
104 
105  G4int i;
106  for (i=0; i<tableSize; i++)
107  {
108  G4double value = ppTable[i] * millibarn;
109  pp->PutValue(i,value);
110  value = npTable[i] * millibarn;
111  np->PutValue(i,value);
112  }
113  xMap[G4Proton::ProtonDefinition()] = pp;
114  xMap[G4Neutron::NeutronDefinition()] = np;
115 }
116 
117 
119 {
120  delete xMap[G4Proton::ProtonDefinition()];
121  delete xMap[G4Neutron::NeutronDefinition()];
122 }
123 
124 
126 {
127  return (this == (G4XNNElasticLowE *) &right);
128 }
129 
130 
132 {
133 
134  return (this != (G4XNNElasticLowE *) &right);
135 }
136 
137 
139 {
140  G4double sigma = 0.;
141  G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
142  G4bool dummy = false;
143 
144  G4ParticleDefinition * key = FindKeyParticle(trk1,trk2);
145 
146  typedef std::map <G4ParticleDefinition *, G4PhysicsVector*, std::less<G4ParticleDefinition *> > StringPhysMap;
147 
148  if (xMap.find(key)!= xMap.end())
149  {
150 
151  StringPhysMap::const_iterator iter;
152  for (iter = xMap.begin(); iter != xMap.end(); ++iter)
153  {
154  G4ParticleDefinition * str = (*iter).first;
155  if (str == key)
156  {
157  G4PhysicsVector* physVector = (*iter).second;
158  // G4PhysicsVector* physVector = xMap[key];
159  if (sqrtS >= _eMin && sqrtS <= _eMax)
160  {
161  sigma = physVector->GetValue(sqrtS,dummy);
162  } else if ( sqrtS < _eMin )
163  {
164  sigma = physVector->GetValue(_eMin,dummy);
165  }
166  //G4cout << " sqrtS / sigma " << sqrtS/GeV << " / " <<
167  // sigma/millibarn << G4endl;
168  }
169  }
170  }
171  return sigma;
172 }
173 
174 
176 {
177  // Dump the pp cross-section table
178 
179  G4cout << Name() << ", pp cross-section: " << G4endl;
180 
181  G4bool dummy = false;
182  G4int i;
184  G4PhysicsVector* pp = 0;
185 
186  typedef std::map <G4ParticleDefinition *, G4PhysicsVector*, std::less<G4ParticleDefinition *> > StringPhysMap;
187  StringPhysMap::const_iterator iter;
188 
189  for (iter = xMap.begin(); iter != xMap.end(); ++iter)
190  {
191  G4ParticleDefinition * str = (*iter).first;
192  if (str == key)
193  {
194  pp = (*iter).second;
195  }
196  }
197 
198  if (pp != 0)
199  {
200  for (i=0; i<tableSize; i++)
201  {
202  G4double e = pp->GetLowEdgeEnergy(i);
203  G4double sigma = pp->GetValue(e,dummy) / millibarn;
204  G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
205  }
206  }
207 
208  // Dump the np cross-section table
209 
210  G4cout << Name() << ", np cross-section: " << G4endl;
211 
213  G4PhysicsVector* np = 0;
214  for (iter = xMap.begin(); iter != xMap.end(); ++iter)
215  {
216  G4ParticleDefinition * str = (*iter).first;
217  if (str == key)
218  {
219  np = (*iter).second;
220  }
221  }
222 
223  // G4PhysicsVector* np = xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
224 
225  if (np != 0)
226  {
227  for (i=0; i<tableSize; i++)
228  {
229  G4double e = np->GetLowEdgeEnergy(i);
230  G4double sigma = np->GetValue(e,dummy) / millibarn;
231  G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
232  }
233  }
235 }
236 
237 
239 {
240  G4String name("NNElasticLowE");
241  return name;
242 }
243 
244 
245 
247 {
248  G4bool answer = InLimits(e,_lowLimit,_highLimit);
249 
250  return answer;
251 }
252 
253 
G4double GetValue(G4double theEnergy, G4bool &isOutRange) const
G4bool operator==(const G4XNNElasticLowE &right) const
virtual G4double CrossSection(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
const XML_Char * name
Definition: expat.h:151
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
int millibarn
Definition: hepunit.py:40
virtual ~G4XNNElasticLowE()
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4ParticleDefinition * FindKeyParticle(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
void PutValue(size_t index, G4double theValue)
virtual void Print() const
virtual G4bool IsValid(G4double e) const
virtual void Print() const
const XML_Char int const XML_Char * value
Definition: expat.h:331
#define G4endl
Definition: G4ios.hh:61
virtual G4String Name() const
double G4double
Definition: G4Types.hh:76
G4bool operator!=(const G4XNNElasticLowE &right) const
const G4LorentzVector & Get4Momentum() const
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
G4bool InLimits(G4double e, G4double eLow, G4double eHigh) const