source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPThermalScatteringData.cc@ 1201

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

import all except CVS
File size: 11.0 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// Thermal Neutron Scattering
27// Koi, Tatsumi (SCCS/SLAC)
28//
29// Class Description
30// Cross Sections for a high precision (based on evaluated data
31// libraries) description of themal neutron scattering below 4 eV;
32// Based on Thermal neutron scattering files
33// from the evaluated nuclear data files ENDF/B-VI, Release2
34// To be used in your physics list in case you need this physics.
35// In this case you want to register an object of this class with
36// the corresponding process.
37// Class Description - End
38
39// 15-Nov-06 First implementation is done by T. Koi (SLAC/SCCS)
40// 070625 implement clearCurrentXSData to fix memory leaking by T. Koi
41
42#include "G4NeutronHPThermalScatteringData.hh"
43#include "G4Neutron.hh"
44#include "G4ElementTable.hh"
45//#include "G4NeutronHPData.hh"
46
47
48
49G4NeutronHPThermalScatteringData::G4NeutronHPThermalScatteringData()
50{
51// Upper limit of neutron energy
52 emax = 4*eV;
53
54 indexOfThermalElement.clear();
55
56 names = new G4NeutronHPThermalScatteringNames();
57
58 BuildPhysicsTable( *G4Neutron::Neutron() );
59}
60
61
62
63G4NeutronHPThermalScatteringData::~G4NeutronHPThermalScatteringData()
64{
65
66 clearCurrentXSData();
67
68 delete names;
69}
70
71
72void G4NeutronHPThermalScatteringData::clearCurrentXSData()
73{
74 std::map< G4int , std::map< G4double , G4NeutronHPVector* >* >::iterator it;
75 std::map< G4double , G4NeutronHPVector* >::iterator itt;
76
77 for ( it = coherent.begin() ; it != coherent.end() ; it++ )
78 {
79 if ( it->second != NULL )
80 {
81 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
82 {
83 delete itt->second;
84 }
85 }
86 delete it->second;
87 }
88
89 for ( it = incoherent.begin() ; it != incoherent.end() ; it++ )
90 {
91 if ( it->second != NULL )
92 {
93 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
94 {
95 delete itt->second;
96 }
97 }
98 delete it->second;
99 }
100
101 for ( it = inelastic.begin() ; it != inelastic.end() ; it++ )
102 {
103 if ( it->second != NULL )
104 {
105 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
106 {
107 delete itt->second;
108 }
109 }
110 delete it->second;
111 }
112
113 coherent.clear();
114 incoherent.clear();
115 inelastic.clear();
116
117}
118
119
120
121G4bool G4NeutronHPThermalScatteringData::IsApplicable( const G4DynamicParticle* aP , const G4Element* anEle )
122{
123 G4bool result = false;
124
125 G4double eKin = aP->GetKineticEnergy();
126 // Check energy
127 if ( eKin < emax )
128 {
129 // Check Particle Species
130 if ( aP->GetDefinition() == G4Neutron::Neutron() )
131 {
132 // anEle is one of Thermal elements
133 G4int ie = (G4int) anEle->GetIndex();
134 std::vector < G4int >::iterator it;
135 for ( it = indexOfThermalElement.begin() ; it != indexOfThermalElement.end() ; it++ )
136 {
137 if ( ie == *it ) return true;
138 }
139 }
140 }
141
142/*
143 if ( names->IsThisThermalElement ( anEle->GetName() ) )
144 {
145 // Check energy and projectile species
146 G4double eKin = aP->GetKineticEnergy();
147 if ( eKin < emax && aP->GetDefinition() == G4Neutron::Neutron() ) result = true;
148 }
149*/
150 return result;
151}
152
153
154void G4NeutronHPThermalScatteringData::BuildPhysicsTable(const G4ParticleDefinition& aP)
155{
156
157 if ( &aP != G4Neutron::Neutron() )
158 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
159
160 indexOfThermalElement.clear();
161
162 clearCurrentXSData();
163
164 static const G4ElementTable* theElementTable = G4Element::GetElementTable();
165 size_t numberOfElements = G4Element::GetNumberOfElements();
166 size_t numberOfThermalElements = 0;
167 for ( size_t i = 0 ; i < numberOfElements ; i++ )
168 {
169 if ( names->IsThisThermalElement ( (*theElementTable)[i]->GetName() ) )
170 {
171 indexOfThermalElement.push_back( i );
172 numberOfThermalElements++;
173 }
174 }
175
176 // Read Cross Section Data files
177
178 G4String dirName;
179 if ( !getenv( "G4NEUTRONHPDATA" ) )
180 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
181 G4String baseName = getenv( "G4NEUTRONHPDATA" );
182
183 dirName = baseName + "/ThermalScattering";
184
185 G4String ndl_filename;
186 G4String name;
187
188 for ( size_t i = 0 ; i < numberOfThermalElements ; i++ )
189 {
190 ndl_filename = names->GetTS_NDL_Name( (*theElementTable)[ indexOfThermalElement[ i ] ]->GetName() );
191
192 // Coherent
193 name = dirName + "/Coherent/CrossSection/" + ndl_filename;
194 std::map< G4double , G4NeutronHPVector* >* coh_amapTemp_EnergyCross = readData( name );
195 coherent.insert ( std::pair < G4int , std::map< G4double , G4NeutronHPVector* >* > ( indexOfThermalElement[ i ] , coh_amapTemp_EnergyCross ) );
196
197 // Incoherent
198 name = dirName + "/Incoherent/CrossSection/" + ndl_filename;
199 std::map< G4double , G4NeutronHPVector* >* incoh_amapTemp_EnergyCross = readData( name );
200 incoherent.insert ( std::pair < G4int , std::map< G4double , G4NeutronHPVector* >* > ( indexOfThermalElement[ i ] , incoh_amapTemp_EnergyCross ) );
201
202 // Inelastic
203 name = dirName + "/Inelastic/CrossSection/" + ndl_filename;
204 std::map< G4double , G4NeutronHPVector* >* inela_amapTemp_EnergyCross = readData( name );
205 inelastic.insert ( std::pair < G4int , std::map< G4double , G4NeutronHPVector* >* > ( indexOfThermalElement[ i ] , inela_amapTemp_EnergyCross ) );
206 }
207
208}
209
210
211
212std::map< G4double , G4NeutronHPVector* >* G4NeutronHPThermalScatteringData::readData ( G4String name )
213{
214
215 std::map< G4double , G4NeutronHPVector* >* aData = new std::map< G4double , G4NeutronHPVector* >;
216
217 std::ifstream theChannel( name.c_str() );
218
219 //G4cout << "G4NeutronHPThermalScatteringData " << name << G4endl;
220
221 G4int dummy;
222 while ( theChannel >> dummy ) // MF
223 {
224 theChannel >> dummy; // MT
225 G4double temp;
226 theChannel >> temp;
227 G4NeutronHPVector* anEnergyCross = new G4NeutronHPVector;
228 G4int nData;
229 theChannel >> nData;
230 anEnergyCross->Init ( theChannel , nData , eV , barn );
231 aData->insert ( std::pair < G4double , G4NeutronHPVector* > ( temp , anEnergyCross ) );
232 }
233 theChannel.close();
234
235 return aData;
236
237}
238
239
240
241void G4NeutronHPThermalScatteringData::DumpPhysicsTable( const G4ParticleDefinition& aP )
242{
243 if( &aP != G4Neutron::Neutron() )
244 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
245// G4cout << "G4NeutronHPThermalScatteringData::DumpPhysicsTable still to be implemented"<<G4endl;
246}
247
248//#include "G4Nucleus.hh"
249//#include "G4NucleiPropertiesTable.hh"
250//#include "G4Neutron.hh"
251//#include "G4Electron.hh"
252
253
254
255G4double G4NeutronHPThermalScatteringData::GetCrossSection( const G4DynamicParticle* aP , const G4Element*anE , G4double aT )
256{
257 G4double result = 0;
258
259 G4int iele = anE->GetIndex();
260
261 G4double Xcoh = GetX ( aP , aT , coherent.find(iele)->second );
262 G4double Xincoh = GetX ( aP , aT , incoherent.find(iele)->second );
263 G4double Xinela = GetX ( aP , aT , inelastic.find(iele)->second );
264
265 result = Xcoh + Xincoh + Xinela;
266
267 //G4cout << "G4NeutronHPThermalScatteringData::GetCrossSection Tot= " << result/barn << " Coherent= " << Xcoh/barn << " Incoherent= " << Xincoh/barn << " Inelastic= " << Xinela/barn << G4endl;
268
269 return result;
270}
271
272
273
274G4double G4NeutronHPThermalScatteringData::GetInelasticCrossSection( const G4DynamicParticle* aP , const G4Element*anE , G4double aT )
275{
276 G4double result = 0;
277 G4int iele = anE->GetIndex();
278 result = GetX ( aP , aT , inelastic.find(iele)->second );
279 return result;
280}
281
282
283
284G4double G4NeutronHPThermalScatteringData::GetCoherentCrossSection( const G4DynamicParticle* aP , const G4Element*anE , G4double aT )
285{
286 G4double result = 0;
287 G4int iele = anE->GetIndex();
288 result = GetX ( aP , aT , coherent.find(iele)->second );
289 return result;
290}
291
292
293
294G4double G4NeutronHPThermalScatteringData::GetIncoherentCrossSection( const G4DynamicParticle* aP , const G4Element*anE , G4double aT )
295{
296 G4double result = 0;
297 G4int iele = anE->GetIndex();
298 result = GetX ( aP , aT , incoherent.find(iele)->second );
299 return result;
300}
301
302
303
304
305G4double G4NeutronHPThermalScatteringData::GetX ( const G4DynamicParticle* aP, G4double aT , std::map < G4double , G4NeutronHPVector* >* amapTemp_EnergyCross )
306{
307 G4double result = 0;
308 if ( amapTemp_EnergyCross->size() == 0 ) return result;
309
310 std::map< G4double , G4NeutronHPVector* >::iterator it;
311 for ( it = amapTemp_EnergyCross->begin() ; it != amapTemp_EnergyCross->end() ; it++ )
312 {
313 if ( aT < it->first ) break;
314 }
315 if ( it == amapTemp_EnergyCross->begin() ) it++; // lower than first
316 else if ( it == amapTemp_EnergyCross->end() ) it--; // upper than last
317
318 G4double eKinetic = aP->GetKineticEnergy();
319
320 G4double TH = it->first;
321 G4double XH = it->second->GetXsec ( eKinetic );
322
323 //G4cout << "G4NeutronHPThermalScatteringData::GetX TH " << TH << " E " << eKinetic << " XH " << XH << G4endl;
324
325 it--;
326 G4double TL = it->first;
327 G4double XL = it->second->GetXsec ( eKinetic );
328
329 //G4cout << "G4NeutronHPThermalScatteringData::GetX TL " << TL << " E " << eKinetic << " XL " << XL << G4endl;
330
331 if ( TH == TL )
332 throw G4HadronicException(__FILE__, __LINE__, "Thermal Scattering Data Error!");
333
334 G4double T = aT;
335 G4double X = ( XH - XL ) / ( TH - TL ) * ( T - TL ) + XL;
336 result = X;
337
338 return result;
339}
Note: See TracBrowser for help on using the repository browser.