source: trunk/source/processes/hadronic/models/im_r_matrix/src/G4XNNElasticLowE.cc@ 1007

Last change on this file since 1007 was 962, checked in by garnier, 17 years ago

update processes
File size: 8.1 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#include "globals.hh"
27#include "G4ios.hh"
28#include "G4XNNElasticLowE.hh"
29#include "G4KineticTrack.hh"
30#include "G4ParticleDefinition.hh"
31#include "G4PhysicsVector.hh"
32#include "G4PhysicsLnVector.hh"
33#include "G4ParticleDefinition.hh"
34#include "G4Proton.hh"
35#include "G4Neutron.hh"
36
37
38const G4double G4XNNElasticLowE::_lowLimit = 0.;
39const 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
43const G4double G4XNNElasticLowE::_eMinTable = 1.8964808;
44const G4double G4XNNElasticLowE::_eStepLog = 0.01;
45
46// Cross-sections in mb
47// Units are assigned while filling the PhysicsVector
48
49const G4int G4XNNElasticLowE::tableSize = 101;
50
51const 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
69const 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
88G4XNNElasticLowE::G4XNNElasticLowE()
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()->GetParticleName()] = pp;
114 xMap[G4Neutron::NeutronDefinition()->GetParticleName()] = np;
115}
116
117
118G4XNNElasticLowE::~G4XNNElasticLowE()
119{
120 G4String name_proton = "proton";
121 G4String name_neutron = "neutron";
122 delete xMap[name_proton];
123 delete xMap[name_neutron];
124}
125
126
127G4bool G4XNNElasticLowE::operator==(const G4XNNElasticLowE &right) const
128{
129 return (this == (G4XNNElasticLowE *) &right);
130}
131
132
133G4bool G4XNNElasticLowE::operator!=(const G4XNNElasticLowE &right) const
134{
135
136 return (this != (G4XNNElasticLowE *) &right);
137}
138
139
140G4double G4XNNElasticLowE::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
141{
142 G4double sigma = 0.;
143 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
144 G4bool dummy = false;
145
146 G4String key = FindKeyParticle(trk1,trk2);
147
148 typedef std::map <G4String, G4PhysicsVector*, std::less<G4String> > StringPhysMap;
149
150 if (xMap.find(key)!= xMap.end())
151 {
152
153 StringPhysMap::const_iterator iter;
154 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
155 {
156 G4String str = (*iter).first;
157 if (str == key)
158 {
159 G4PhysicsVector* physVector = (*iter).second;
160 // G4PhysicsVector* physVector = xMap[key];
161 if (sqrtS >= _eMin && sqrtS <= _eMax)
162 {
163 sigma = physVector->GetValue(sqrtS,dummy);
164 } else if ( sqrtS < _eMin )
165 {
166 sigma = physVector->GetValue(_eMin,dummy);
167 }
168 }
169 }
170 }
171 return sigma;
172}
173
174
175void G4XNNElasticLowE::Print() const
176{
177 // Dump the pp cross-section table
178
179 G4cout << Name() << ", pp cross-section: " << G4endl;
180
181 G4bool dummy = false;
182 G4int i;
183 G4String key = G4Proton::ProtonDefinition()->GetParticleName();
184 G4PhysicsVector* pp = 0;
185
186 typedef std::map <G4String, G4PhysicsVector*, std::less<G4String> > StringPhysMap;
187 StringPhysMap::const_iterator iter;
188
189 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
190 {
191 G4String 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
212 key = G4Neutron::NeutronDefinition()->GetParticleName();
213 G4PhysicsVector* np = 0;
214 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
215 {
216 G4String 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 }
234 G4VCrossSectionSource::Print();
235}
236
237
238G4String G4XNNElasticLowE::Name() const
239{
240 G4String name("NNElasticLowE");
241 return name;
242}
243
244
245
246G4bool G4XNNElasticLowE::IsValid(G4double e) const
247{
248 G4bool answer = InLimits(e,_lowLimit,_highLimit);
249
250 return answer;
251}
252
253
Note: See TracBrowser for help on using the repository browser.