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

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

import all except CVS
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 delete xMap[G4Proton::ProtonDefinition()->GetParticleName()];
121 delete xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
122}
123
124
125G4bool G4XNNElasticLowE::operator==(const G4XNNElasticLowE &right) const
126{
127 return (this == (G4XNNElasticLowE *) &right);
128}
129
130
131G4bool G4XNNElasticLowE::operator!=(const G4XNNElasticLowE &right) const
132{
133
134 return (this != (G4XNNElasticLowE *) &right);
135}
136
137
138G4double G4XNNElasticLowE::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
139{
140 G4double sigma = 0.;
141 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
142 G4bool dummy = false;
143
144 G4String key = FindKeyParticle(trk1,trk2);
145
146 typedef std::map <G4String, G4PhysicsVector*, std::less<G4String> > 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 G4String 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 }
167 }
168 }
169 return sigma;
170}
171
172
173void G4XNNElasticLowE::Print() const
174{
175 // Dump the pp cross-section table
176
177 G4cout << Name() << ", pp cross-section: " << G4endl;
178
179 G4bool dummy = false;
180 G4int i;
181 G4String key = G4Proton::ProtonDefinition()->GetParticleName();
182 G4PhysicsVector* pp = 0;
183
184 typedef std::map <G4String, G4PhysicsVector*, std::less<G4String> > StringPhysMap;
185 StringPhysMap::const_iterator iter;
186
187 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
188 {
189 G4String str = (*iter).first;
190 if (str == key)
191 {
192 pp = (*iter).second;
193 }
194 }
195
196 if (pp != 0)
197 {
198 for (i=0; i<tableSize; i++)
199 {
200 G4double e = pp->GetLowEdgeEnergy(i);
201 G4double sigma = pp->GetValue(e,dummy) / millibarn;
202 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
203 }
204 }
205
206 // Dump the np cross-section table
207
208 G4cout << Name() << ", np cross-section: " << G4endl;
209
210 key = G4Neutron::NeutronDefinition()->GetParticleName();
211 G4PhysicsVector* np = 0;
212 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
213 {
214 G4String str = (*iter).first;
215 if (str == key)
216 {
217 np = (*iter).second;
218 }
219 }
220
221 // G4PhysicsVector* np = xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
222
223 if (np != 0)
224 {
225 for (i=0; i<tableSize; i++)
226 {
227 G4double e = np->GetLowEdgeEnergy(i);
228 G4double sigma = np->GetValue(e,dummy) / millibarn;
229 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
230 }
231 }
232 G4VCrossSectionSource::Print();
233}
234
235
236G4String G4XNNElasticLowE::Name() const
237{
238 G4String name("NNElasticLowE");
239 return name;
240}
241
242
243
244G4bool G4XNNElasticLowE::IsValid(G4double e) const
245{
246 G4bool answer = InLimits(e,_lowLimit,_highLimit);
247
248 return answer;
249}
250
251
Note: See TracBrowser for help on using the repository browser.