source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPJENDLHEData.cc@ 1036

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

import all except CVS
File size: 11.7 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// Class Description
27// Cross-section data set for a high precision (based on JENDL_HE evaluated data
28// libraries) description of elastic scattering 20 MeV ~ 3 GeV;
29// Class Description - End
30
31// 15-Nov-06 First Implementation is done by T. Koi (SLAC/SCCS)
32
33#include "G4NeutronHPJENDLHEData.hh"
34#include "G4LPhysicsFreeVector.hh"
35#include "G4ElementTable.hh"
36#include "G4NeutronHPData.hh"
37
38G4bool G4NeutronHPJENDLHEData::IsApplicable(const G4DynamicParticle*aP, const G4Element* anE)
39{
40
41 G4bool result = true;
42 G4double eKin = aP->GetKineticEnergy();
43 //if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
44 if ( eKin < 20*MeV || 3*GeV < eKin || aP->GetDefinition()!=G4Neutron::Neutron() )
45 {
46 result = false;
47 }
48// Element Check
49 else if ( !(vElement[ anE->GetIndex() ]) ) result = false;
50
51 return result;
52
53}
54
55
56
57G4NeutronHPJENDLHEData::G4NeutronHPJENDLHEData()
58{
59 ;
60}
61
62
63
64G4NeutronHPJENDLHEData::G4NeutronHPJENDLHEData( G4String reaction , G4ParticleDefinition* pd )
65{
66 reactionName = reaction;
67 BuildPhysicsTable( *pd );
68}
69
70
71
72G4NeutronHPJENDLHEData::~G4NeutronHPJENDLHEData()
73{
74 ;
75 //delete theCrossSections;
76}
77
78
79
80void G4NeutronHPJENDLHEData::BuildPhysicsTable( const G4ParticleDefinition& aP )
81{
82
83// if ( &aP != G4Neutron::Neutron() )
84// throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
85 particleName = aP.GetParticleName();
86
87 G4String baseName = getenv( "G4NEUTRONHPDATA" );
88 G4String dirName = baseName+"/JENDL_HE/"+particleName+"/"+reactionName ;
89 G4String aFSType = "/CrossSection/";
90 G4NeutronHPNames theNames;
91
92 G4String filename;
93
94// Create JENDL_HE data
95// Create map element or isotope
96
97 size_t numberOfElements = G4Element::GetNumberOfElements();
98 //theCrossSections = new G4PhysicsTable( numberOfElements );
99
100 // make a PhysicsVector for each element
101
102 static const G4ElementTable *theElementTable = G4Element::GetElementTable();
103 vElement.clear();
104 vElement.resize( numberOfElements );
105 for ( size_t i = 0; i < numberOfElements; ++i )
106 {
107
108 G4Element* theElement = (*theElementTable)[i];
109 vElement[i] = false;
110
111 // isotope
112 G4int nIso = (*theElementTable)[i]->GetNumberOfIsotopes();
113 G4int Z = static_cast<G4int> ((*theElementTable)[i]->GetZ());
114 if ( nIso!=0 )
115 {
116 G4bool found_at_least_one = false;
117 for ( G4int i1 = 0; i1 < nIso; i1++ )
118 {
119 G4int A = theElement->GetIsotope(i1)->GetN();
120
121 if ( isThisNewIsotope( Z , A ) )
122 {
123
124 std::stringstream ss;
125 ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
126 filename = ss.str();
127 std::fstream file;
128 file.open ( filename , std::fstream::in );
129 G4int dummy;
130 file >> dummy;
131 if ( file.good() )
132 {
133
134 //G4cout << "Found file for Z=" << Z << ", A=" << A << ", as " << filename << G4endl;
135 found_at_least_one = true;
136
137 // read the file
138 G4PhysicsVector* aPhysVec = readAFile ( &file );
139
140 //Regist
141
142 registAPhysicsVector( Z , A , aPhysVec );
143
144 }
145 else
146 {
147 //G4cout << "No file for "<< reactionType << " Z=" << Z << ", A=" << A << G4endl;
148 }
149
150 file.close();
151
152 }
153 else
154 {
155 found_at_least_one = TRUE;
156 }
157 }
158
159 if ( found_at_least_one ) vElement[i] = true;
160
161 }
162 else
163 {
164 G4StableIsotopes theStableOnes;
165 G4int first = theStableOnes.GetFirstIsotope( Z );
166 G4bool found_at_least_one = FALSE;
167 for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(theElement->GetZ() ) ); i1++)
168 {
169 G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
170
171 if ( isThisNewIsotope( Z , A ) )
172 {
173
174 std::stringstream ss;
175 ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
176 filename = ss.str();
177
178 std::fstream file;
179 file.open ( filename , std::fstream::in );
180 G4int dummy;
181 file >> dummy;
182 if ( file.good() )
183 {
184 //G4cout << "Found file for Z=" << Z << ", A=" << A << ", as " << filename << G4endl;
185 found_at_least_one = TRUE;
186 //Read the file
187
188 G4PhysicsVector* aPhysVec = readAFile ( &file );
189
190 //Regist the PhysicsVector
191 registAPhysicsVector( Z , A , aPhysVec );
192
193 }
194 else
195 {
196 //G4cout << "No file for "<< reactionType << " Z=" << Z << ", A=" << A << G4endl;
197 }
198
199 file.close();
200 }
201 else
202 {
203 found_at_least_one = TRUE;
204 }
205 }
206
207 if ( found_at_least_one ) vElement[i] = true;
208
209 }
210
211 }
212
213}
214
215
216
217void G4NeutronHPJENDLHEData::DumpPhysicsTable(const G4ParticleDefinition& aP)
218{
219 if(&aP!=G4Neutron::Neutron())
220 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
221// G4cout << "G4NeutronHPJENDLHEData::DumpPhysicsTable still to be implemented"<<G4endl;
222}
223
224
225
226G4double G4NeutronHPJENDLHEData::
227GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double )
228// aTemp
229{
230
231 // Primary energy >20MeV
232 // Thus
233 // Not take account of Doppler broadening
234 // also
235 // Not take account of Target thermal motions
236
237 G4double result = 0;
238
239 G4double ek = aP->GetKineticEnergy();
240
241 G4int nIso = anE->GetNumberOfIsotopes();
242 G4int Z = static_cast<G4int> ( anE->GetZ() );
243 if ( nIso!=0 )
244 {
245 for ( G4int i1 = 0; i1 < nIso; i1++ )
246 {
247
248 G4int A = anE->GetIsotope(i1)->GetN();
249 G4double frac = anE->GetRelativeAbundanceVector()[ i1 ]; // This case do NOT request "*perCent".
250
251 result += frac * getXSfromThisIsotope( Z , A , ek );
252 //G4cout << reactionType << " XS in barn " << Z << " " << A << " " << frac << " " << getXSfromThisIsotope( Z , A , ek )/barn << G4endl;
253
254 }
255 }
256 else
257 {
258
259 G4StableIsotopes theStableOnes;
260 G4int first = theStableOnes.GetFirstIsotope( Z );
261 for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(anE->GetZ() ) ); i1++)
262 {
263
264 G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
265 G4double frac = theStableOnes.GetAbundance( first+i1 )*perCent; // This case request "*perCent".
266
267 result += frac * getXSfromThisIsotope( Z , A , ek );
268 //G4cout << reactionType << " XS in barn " << Z << " " << A << " " << frac << " " << getXSfromThisIsotope( Z , A , ek )/barn << G4endl;
269
270 }
271 }
272 return result;
273
274}
275
276
277
278G4PhysicsVector* G4NeutronHPJENDLHEData::readAFile ( std::fstream* file )
279{
280
281 G4int dummy;
282 G4int len;
283 *file >> dummy;
284 *file >> len;
285
286 std::vector< G4double > v_e;
287 std::vector< G4double > v_xs;
288
289 for ( G4int i = 0 ; i < len ; i++ )
290 {
291 G4double e;
292 G4double xs;
293
294 *file >> e;
295 *file >> xs;
296 // data are written in eV and barn.
297 v_e.push_back( e*eV );
298 v_xs.push_back( xs*barn );
299 }
300
301 G4LPhysicsFreeVector* aPhysVec = new G4LPhysicsFreeVector( static_cast< size_t >( len ) , v_e.front() , v_e.back() );
302
303 for ( G4int i = 0 ; i < len ; i++ )
304 {
305 aPhysVec->PutValues( static_cast< size_t >( i ) , v_e[ i ] , v_xs[ i ] );
306 }
307
308 return aPhysVec;
309}
310
311
312
313G4bool G4NeutronHPJENDLHEData::isThisInMap( G4int z , G4int a )
314{
315 if ( mIsotope.find ( z ) == mIsotope.end() ) return false;
316 if ( mIsotope.find ( z ) -> second->find ( a ) == mIsotope.find ( z ) -> second->end() ) return false;
317 return true;
318}
319
320
321
322void G4NeutronHPJENDLHEData::registAPhysicsVector( G4int Z , G4int A , G4PhysicsVector* aPhysVec )
323{
324
325 std::pair< G4int , G4PhysicsVector* > aPair = std::pair < G4int , G4PhysicsVector* > ( A , aPhysVec );
326
327 std::map < G4int , std::map< G4int , G4PhysicsVector* >* >::iterator itm;
328 itm = mIsotope.find ( Z );
329 if ( itm != mIsotope.end() )
330 {
331 itm->second->insert ( aPair );
332 }
333 else
334 {
335 std::map< G4int , G4PhysicsVector* >* aMap = new std::map< G4int , G4PhysicsVector* >;
336 aMap->insert ( aPair );
337 mIsotope.insert( std::pair< G4int , std::map< G4int , G4PhysicsVector* >* > ( Z , aMap ) );
338 }
339
340}
341
342
343
344G4double G4NeutronHPJENDLHEData::getXSfromThisIsotope( G4int Z , G4int A , G4double ek )
345{
346
347 G4double aXSection = 0.0;
348 G4bool outOfRange;
349
350 G4PhysicsVector* aPhysVec;
351 if ( mIsotope.find ( Z )->second->find ( A ) != mIsotope.find ( Z )->second->end() )
352 {
353
354 aPhysVec = mIsotope.find ( Z )->second->find ( A )->second;
355 aXSection = aPhysVec->GetValue( ek , outOfRange );
356
357 }
358 else
359 {
360
361 //Select closest one in the same Z
362 std::map < G4int , G4PhysicsVector* >::iterator it;
363 G4int delta0 = 99; // no mean for 99
364 for ( it = mIsotope.find ( Z )->second->begin() ; it != mIsotope.find ( Z )->second->end() ; it++ )
365 {
366 G4int delta = std::abs( A - it->first );
367 if ( delta < delta0 ) delta0 = delta;
368 }
369
370 // Randomize of selection larger or smaller than A
371 if ( G4UniformRand() < 0.5 ) delta0 *= -1;
372 G4int A1 = A + delta0;
373 if ( mIsotope.find ( Z )->second->find ( A1 ) != mIsotope.find ( Z )->second->end() )
374 {
375 aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
376 }
377 else
378 {
379 A1 = A - delta0;
380 aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
381 }
382
383 aXSection = aPhysVec->GetValue( ek , outOfRange );
384 // X^(2/3) factor
385 aXSection *= std::pow ( 1.0*A/ A1 , 2.0 / 3.0 );
386
387 }
388
389 return aXSection;
390}
Note: See TracBrowser for help on using the repository browser.