Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4NeutronHPThermalScattering.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// Thermal Neutron Scattering
27// Koi, Tatsumi (SLAC/SCCS)
28//
29// Class Description:
30//
31// Final State Generators for a high precision (based on evaluated data
32// libraries) description of themal neutron scattering below 4 eV;
33// Based on Thermal neutron scattering files
34// from the evaluated nuclear data files ENDF/B-VI, Release2
35// To be used in your physics list in case you need this physics.
36// In this case you want to register an object of this class with
37// the corresponding process.
38
39
40// 070625 Fix memory leaking at destructor by T. Koi
41// 081201 Fix memory leaking at destructor by T. Koi
42// 100729 Add model name in constructor Problem #1116
43
44#include "G4NeutronHPThermalScattering.hh"
45#include "G4SystemOfUnits.hh"
46#include "G4Neutron.hh"
47#include "G4ElementTable.hh"
48
49G4NeutronHPThermalScattering::G4NeutronHPThermalScattering()
50 :G4HadronicInteraction("NeutronHPThermalScattering")
51{
52
53 theHPElastic = new G4NeutronHPElastic();
54
55 SetMinEnergy( 0.*eV );
56 SetMaxEnergy( 4*eV );
57 theXSection = new G4NeutronHPThermalScatteringData();
58 theXSection->BuildPhysicsTable( *(G4Neutron::Neutron()) );
59
60 buildPhysicsTable();
61
62}
63
64
65
66G4NeutronHPThermalScattering::~G4NeutronHPThermalScattering()
67{
68
69 for ( std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >::iterator it = incoherentFSs.begin() ; it != incoherentFSs.end() ; it++ )
70 {
71 std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
72 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
73 {
74 std::vector< E_isoAng* >::iterator ittt;
75 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
76 {
77 delete *ittt;
78 }
79 delete itt->second;
80 }
81 delete it->second;
82 }
83
84 for ( std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it = coherentFSs.begin() ; it != coherentFSs.end() ; it++ )
85 {
86 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt;
87 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
88 {
89 std::vector < std::pair< G4double , G4double >* >::iterator ittt;
90 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
91 {
92 delete *ittt;
93 }
94 delete itt->second;
95 }
96 delete it->second;
97 }
98
99 for ( std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >::iterator it = inelasticFSs.begin() ; it != inelasticFSs.end() ; it++ )
100 {
101 std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt;
102 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
103 {
104 std::vector < E_P_E_isoAng* >::iterator ittt;
105 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
106 {
107 std::vector < E_isoAng* >::iterator it4;
108 for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ )
109 {
110 delete *it4;
111 }
112 delete *ittt;
113 }
114 delete itt->second;
115 }
116 delete it->second;
117 }
118
119 delete theHPElastic;
120 delete theXSection;
121}
122
123
124
125std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* G4NeutronHPThermalScattering::readACoherentFSDATA( G4String name )
126{
127
128 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* aCoherentFSDATA = new std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >;
129
130 std::ifstream theChannel( name.c_str() );
131
132 std::vector< G4double > vBraggE;
133
134 G4int dummy;
135 while ( theChannel >> dummy ) // MF
136 {
137 theChannel >> dummy; // MT
138 G4double temp;
139 theChannel >> temp;
140 std::vector < std::pair< G4double , G4double >* >* anBragE_P = new std::vector < std::pair< G4double , G4double >* >;
141
142 G4int n;
143 theChannel >> n;
144 for ( G4int i = 0 ; i < n ; i++ )
145 {
146 G4double Ei;
147 G4double Pi;
148 if ( aCoherentFSDATA->size() == 0 )
149 {
150 theChannel >> Ei;
151 vBraggE.push_back( Ei );
152 }
153 else
154 {
155 Ei = vBraggE[ i ];
156 }
157 theChannel >> Pi;
158 anBragE_P->push_back ( new std::pair < G4double , G4double > ( Ei , Pi ) );
159 //G4cout << "Coherent Elastic " << Ei << " " << Pi << G4endl;
160 }
161 aCoherentFSDATA->insert ( std::pair < G4double , std::vector < std::pair< G4double , G4double >* >* > ( temp , anBragE_P ) );
162 }
163
164 return aCoherentFSDATA;
165}
166
167
168
169std::map < G4double , std::vector < E_P_E_isoAng* >* >* G4NeutronHPThermalScattering::readAnInelasticFSDATA ( G4String name )
170{
171 std::map < G4double , std::vector < E_P_E_isoAng* >* >* anT_E_P_E_isoAng = new std::map < G4double , std::vector < E_P_E_isoAng* >* >;
172
173 std::ifstream theChannel( name.c_str() );
174
175 G4int dummy;
176 while ( theChannel >> dummy ) // MF
177 {
178 theChannel >> dummy; // MT
179 G4double temp;
180 theChannel >> temp;
181 std::vector < E_P_E_isoAng* >* vE_P_E_isoAng = new std::vector < E_P_E_isoAng* >;
182 G4int n;
183 theChannel >> n;
184 for ( G4int i = 0 ; i < n ; i++ )
185 {
186 vE_P_E_isoAng->push_back ( readAnE_P_E_isoAng ( &theChannel ) );
187 }
188 anT_E_P_E_isoAng->insert ( std::pair < G4double , std::vector < E_P_E_isoAng* >* > ( temp , vE_P_E_isoAng ) );
189 }
190 theChannel.close();
191
192 return anT_E_P_E_isoAng;
193}
194
195
196
197E_P_E_isoAng* G4NeutronHPThermalScattering::readAnE_P_E_isoAng( std::ifstream* file )
198{
199 E_P_E_isoAng* aData = new E_P_E_isoAng;
200
201 G4double dummy;
202 G4double energy;
203 G4int nep , nl;
204 *file >> dummy;
205 *file >> energy;
206 aData->energy = energy*eV;
207 *file >> dummy;
208 *file >> dummy;
209 *file >> nep;
210 *file >> nl;
211 aData->n = nep/nl;
212 for ( G4int i = 0 ; i < aData->n ; i++ )
213 {
214 G4double prob;
215 E_isoAng* anE_isoAng = new E_isoAng;
216 aData->vE_isoAngle.push_back( anE_isoAng );
217 *file >> energy;
218 anE_isoAng->energy = energy*eV;
219 anE_isoAng->n = nl - 2;
220 anE_isoAng->isoAngle.resize( anE_isoAng->n );
221 *file >> prob;
222 aData->prob.push_back( prob );
223 //G4cout << "G4NeutronHPThermalScattering inelastic " << energy/eV << " " << i << " " << prob << " " << aData->prob[ i ] << G4endl;
224 for ( G4int j = 0 ; j < anE_isoAng->n ; j++ )
225 {
226 G4double x;
227 *file >> x;
228 anE_isoAng->isoAngle[j] = x ;
229 //G4cout << "G4NeutronHPThermalScattering inelastic " << x << anE_isoAng->isoAngle[j] << G4endl;
230 }
231 }
232
233 // Calcuate sum_of_provXdEs
234 G4double total = 0;
235 for ( G4int i = 0 ; i < aData->n - 1 ; i++ )
236 {
237 G4double E_L = aData->vE_isoAngle[i]->energy/eV;
238 G4double E_H = aData->vE_isoAngle[i+1]->energy/eV;
239 G4double dE = E_H - E_L;
240 total += ( ( aData->prob[i] ) * dE );
241 }
242 aData->sum_of_probXdEs = total;
243
244 return aData;
245}
246
247
248
249std::map < G4double , std::vector < E_isoAng* >* >* G4NeutronHPThermalScattering::readAnIncoherentFSDATA ( G4String name )
250{
251 std::map < G4double , std::vector < E_isoAng* >* >* T_E = new std::map < G4double , std::vector < E_isoAng* >* >;
252
253 std::ifstream theChannel( name.c_str() );
254
255 G4int dummy;
256 while ( theChannel >> dummy ) // MF
257 {
258 theChannel >> dummy; // MT
259 G4double temp;
260 theChannel >> temp;
261 std::vector < E_isoAng* >* vE_isoAng = new std::vector < E_isoAng* >;
262 G4int n;
263 theChannel >> n;
264 for ( G4int i = 0 ; i < n ; i++ )
265 vE_isoAng->push_back ( readAnE_isoAng( &theChannel ) );
266 T_E->insert ( std::pair < G4double , std::vector < E_isoAng* >* > ( temp , vE_isoAng ) );
267 }
268 theChannel.close();
269
270 return T_E;
271}
272
273
274
275E_isoAng* G4NeutronHPThermalScattering::readAnE_isoAng( std::ifstream* file )
276{
277 E_isoAng* aData = new E_isoAng;
278
279 G4double dummy;
280 G4double energy;
281 G4int n;
282 *file >> dummy;
283 *file >> energy;
284 *file >> dummy;
285 *file >> dummy;
286 *file >> n;
287 *file >> dummy;
288 aData->energy = energy*eV;
289 aData->n = n-2;
290 aData->isoAngle.resize( n );
291
292 *file >> dummy;
293 *file >> dummy;
294 for ( G4int i = 0 ; i < aData->n ; i++ )
295 *file >> aData->isoAngle[i];
296
297 return aData;
298}
299
300
301
302G4HadFinalState* G4NeutronHPThermalScattering::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus )
303{
304
305// Select Element > Reaction >
306
307 const G4Material * theMaterial = aTrack.GetMaterial();
308 G4double aTemp = theMaterial->GetTemperature();
309 G4int n = theMaterial->GetNumberOfElements();
310 //static const G4ElementTable* theElementTable = G4Element::GetElementTable();
311
312 G4bool findThermalElement = false;
313 G4int ielement;
314 const G4Element* theElement = NULL;
315 for ( G4int i = 0; i < n ; i++ )
316 {
317 theElement = theMaterial->GetElement(i);
318 //Select target element
319 if ( aNucleus.GetZ_asInt() == (G4int)(theElement->GetZ() + 0.5 ) )
320 {
321 //Check Applicability of Thermal Scattering
322 if ( getTS_ID( NULL , theElement ) != -1 )
323 {
324 ielement = getTS_ID( NULL , theElement );
325 findThermalElement = true;
326 break;
327 }
328 else if ( getTS_ID( theMaterial , theElement ) != -1 )
329 {
330 ielement = getTS_ID( theMaterial , theElement );
331 findThermalElement = true;
332 break;
333 }
334 }
335 }
336
337 if ( findThermalElement == true )
338 {
339
340// Select Reaction (Inelastic, coherent, incoherent)
341
342 G4ParticleDefinition* pd = const_cast< G4ParticleDefinition* >( aTrack.GetDefinition() );
343 G4DynamicParticle* dp = new G4DynamicParticle ( pd , aTrack.Get4Momentum() );
344 G4double total = theXSection->GetCrossSection( dp , theElement , theMaterial );
345 G4double inelastic = theXSection->GetInelasticCrossSection( dp , theElement , theMaterial );
346
347
348 G4double random = G4UniformRand();
349 if ( random <= inelastic/total )
350 {
351 // Inelastic
352
353 // T_L and T_H
354 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator it;
355 std::vector<G4double> v_temp;
356 v_temp.clear();
357 for ( it = inelasticFSs.find( ielement )->second->begin() ; it != inelasticFSs.find( ielement )->second->end() ; it++ )
358 {
359 v_temp.push_back( it->first );
360 }
361
362// T_L T_H
363 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
364//
365// For T_L aNEP_EPM_TL and T_H aNEP_EPM_TH
366//
367 std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = 0;
368 std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = 0;
369
370 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
371 {
372 vNEP_EPM_TL = inelasticFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
373 vNEP_EPM_TH = inelasticFSs.find( ielement )->second->find ( tempLH.second/kelvin )->second;
374 }
375 else if ( tempLH.first == 0.0 )
376 {
377 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
378 itm = inelasticFSs.find( ielement )->second->begin();
379 vNEP_EPM_TL = itm->second;
380 itm++;
381 vNEP_EPM_TH = itm->second;
382 }
383 else if ( tempLH.second == 0.0 )
384 {
385 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
386 itm = inelasticFSs.find( ielement )->second->end();
387 itm--;
388 vNEP_EPM_TH = itm->second;
389 itm--;
390 vNEP_EPM_TL = itm->second;
391 }
392
393//
394
395 G4double rand_for_sE = G4UniformRand();
396
397 std::pair< G4double , E_isoAng > TL = create_sE_and_EPM_from_pE_and_vE_P_E_isoAng ( rand_for_sE , aTrack.GetKineticEnergy() , vNEP_EPM_TL );
398 std::pair< G4double , E_isoAng > TH = create_sE_and_EPM_from_pE_and_vE_P_E_isoAng ( rand_for_sE , aTrack.GetKineticEnergy() , vNEP_EPM_TH );
399
400 G4double sE;
401 sE = get_linear_interpolated ( aTemp , std::pair < G4double , G4double > ( tempLH.first , TL.first ) , std::pair < G4double , G4double > ( tempLH.second , TH.first ) );
402 E_isoAng anE_isoAng;
403 if ( TL.second.n == TH.second.n )
404 {
405 anE_isoAng.energy = sE;
406 anE_isoAng.n = TL.second.n;
407 for ( G4int i=0 ; i < anE_isoAng.n ; i++ )
408 {
409 G4double angle;
410 angle = get_linear_interpolated ( aTemp , std::pair< G4double , G4double > ( tempLH.first , TL.second.isoAngle[ i ] ) , std::pair< G4double , G4double > ( tempLH.second , TH.second.isoAngle[ i ] ) );
411 anE_isoAng.isoAngle.push_back( angle );
412 }
413 }
414 else
415 {
416 std::cout << "Do not Suuport yet." << std::endl;
417 }
418
419 //set
420 theParticleChange.SetEnergyChange( sE );
421 G4double mu = getMu( &anE_isoAng );
422 theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
423
424 }
425 //else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , (*theElementTable)[ ielement ] , aTemp ) ) / total )
426 else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , theElement , theMaterial ) ) / total )
427 {
428 // Coherent Elastic
429
430 G4double E = aTrack.GetKineticEnergy();
431
432 // T_L and T_H
433 std::map < G4double , std::vector< std::pair< G4double , G4double >* >* >::iterator it;
434 std::vector<G4double> v_temp;
435 v_temp.clear();
436 for ( it = coherentFSs.find( ielement )->second->begin() ; it != coherentFSs.find( ielement )->second->end() ; it++ )
437 {
438 v_temp.push_back( it->first );
439 }
440
441// T_L T_H
442 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
443//
444//
445// For T_L anEPM_TL and T_H anEPM_TH
446//
447 std::vector< std::pair< G4double , G4double >* >* pvE_p_TL = 0;
448 std::vector< std::pair< G4double , G4double >* >* pvE_p_TH = 0;
449
450 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
451 {
452 pvE_p_TL = coherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
453 pvE_p_TH = coherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
454 }
455 else if ( tempLH.first == 0.0 )
456 {
457 pvE_p_TL = coherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second;
458 pvE_p_TH = coherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second;
459 }
460 else if ( tempLH.second == 0.0 )
461 {
462 pvE_p_TL = coherentFSs.find( ielement )->second->find ( v_temp.back() )->second;
463 std::vector< G4double >::iterator itv;
464 itv = v_temp.end();
465 itv--;
466 itv--;
467 pvE_p_TL = coherentFSs.find( ielement )->second->find ( *itv )->second;
468 }
469
470
471 std::vector< G4double > vE_T;
472 std::vector< G4double > vp_T;
473
474 G4int n1 = pvE_p_TL->size();
475 //G4int n2 = pvE_p_TH->size();
476
477 for ( G4int i=1 ; i < n1 ; i++ )
478 {
479 if ( (*pvE_p_TL)[i]->first != (*pvE_p_TH)[i]->first ) abort();
480 vE_T.push_back ( (*pvE_p_TL)[i]->first );
481 vp_T.push_back ( get_linear_interpolated ( aTemp , std::pair< G4double , G4double > ( tempLH.first , (*pvE_p_TL)[i]->second ) , std::pair< G4double , G4double > ( tempLH.second , (*pvE_p_TL)[i]->second ) ) );
482 }
483
484 G4int j = 0;
485 for ( G4int i = 1 ; i < n ; i++ )
486 {
487 if ( E/eV < vE_T[ i ] )
488 {
489 j = i-1;
490 break;
491 }
492 }
493
494 G4double rand_for_mu = G4UniformRand();
495
496 G4int k = 0;
497 for ( G4int i = 1 ; i < j ; i++ )
498 {
499 G4double Pi = vp_T[ i ] / vp_T[ j ];
500 if ( rand_for_mu < Pi )
501 {
502 k = i-1;
503 break;
504 }
505 }
506
507 //G4double Ei = vE_T[ j ];
508 G4double Ei = vE_T[ k ];
509
510 G4double mu = 1 - 2 * Ei / (E/eV) ;
511 //111102
512 if ( mu < -1.0 ) mu = -1.0;
513
514 theParticleChange.SetEnergyChange( E );
515 theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
516
517
518 }
519 else
520 {
521 // InCoherent Elastic
522
523 // T_L and T_H
524 std::map < G4double , std::vector < E_isoAng* >* >::iterator it;
525 std::vector<G4double> v_temp;
526 v_temp.clear();
527 for ( it = incoherentFSs.find( ielement )->second->begin() ; it != incoherentFSs.find( ielement )->second->end() ; it++ )
528 {
529 v_temp.push_back( it->first );
530 }
531
532// T_L T_H
533 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
534
535//
536// For T_L anEPM_TL and T_H anEPM_TH
537//
538
539 E_isoAng anEPM_TL_E;
540 E_isoAng anEPM_TH_E;
541
542 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
543 {
544 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second );
545 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( tempLH.second/kelvin )->second );
546 }
547 else if ( tempLH.first == 0.0 )
548 {
549 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second );
550 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second );
551 }
552 else if ( tempLH.second == 0.0 )
553 {
554 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp.back() )->second );
555 std::vector< G4double >::iterator itv;
556 itv = v_temp.end();
557 itv--;
558 itv--;
559 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( *itv )->second );
560 }
561
562 // E_isoAng for aTemp and aTrack.GetKineticEnergy()
563 E_isoAng anEPM_T_E;
564
565 if ( anEPM_TL_E.n == anEPM_TH_E.n )
566 {
567 anEPM_T_E.n = anEPM_TL_E.n;
568 for ( G4int i=0 ; i < anEPM_TL_E.n ; i++ )
569 {
570 G4double angle;
571 angle = get_linear_interpolated ( aTemp , std::pair< G4double , G4double > ( tempLH.first , anEPM_TL_E.isoAngle[ i ] ) , std::pair< G4double , G4double > ( tempLH.second , anEPM_TH_E.isoAngle[ i ] ) );
572 anEPM_T_E.isoAngle.push_back( angle );
573 }
574 }
575 else
576 {
577 std::cout << "Do not Suuport yet." << std::endl;
578 }
579
580 // Decide mu
581 G4double mu = getMu ( &anEPM_T_E );
582
583 // Set Final State
584 theParticleChange.SetEnergyChange( aTrack.GetKineticEnergy() ); // No energy change in Elastic
585 theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
586
587 }
588 delete dp;
589
590 return &theParticleChange;
591
592 }
593 else
594 {
595 // Not thermal element
596 // Neutron HP will handle
597 return theHPElastic -> ApplyYourself( aTrack, aNucleus );
598 }
599
600}
601
602
603
604G4double G4NeutronHPThermalScattering::getMu( E_isoAng* anEPM )
605{
606
607 G4double random = G4UniformRand();
608 G4double result = 0.0;
609
610 G4int in = int ( random * ( (*anEPM).n ) );
611
612 if ( in != 0 )
613 {
614 G4double mu_l = (*anEPM).isoAngle[ in-1 ];
615 G4double mu_h = (*anEPM).isoAngle[ in ];
616 result = ( mu_h - mu_l ) * ( random * ( (*anEPM).n ) - in ) + mu_l;
617 }
618 else
619 {
620 G4double x = random * (*anEPM).n;
621 G4double D = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) + ( 1 - (*anEPM).isoAngle[ (*anEPM).n - 1 ] );
622 G4double ratio = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) / D;
623 if ( x <= ratio )
624 {
625 G4double mu_l = -1;
626 G4double mu_h = (*anEPM).isoAngle[ 0 ];
627 result = ( mu_h - mu_l ) * x + mu_l;
628 }
629 else
630 {
631 G4double mu_l = (*anEPM).isoAngle[ (*anEPM).n - 1 ];
632 G4double mu_h = 1;
633 result = ( mu_h - mu_l ) * x + mu_l;
634 }
635 }
636 return result;
637}
638
639
640
641std::pair < G4double , G4double > G4NeutronHPThermalScattering::find_LH ( G4double x , std::vector< G4double >* aVector )
642{
643 G4double L = 0.0;
644 G4double H = 0.0;
645 std::vector< G4double >::iterator it;
646 for ( it = aVector->begin() ; it != aVector->end() ; it++ )
647 {
648 if ( x <= *it )
649 {
650 H = *it;
651 if ( it != aVector->begin() )
652 {
653 it--;
654 L = *it;
655 }
656 else
657 {
658 L = 0.0;
659 }
660 break;
661 }
662 }
663 if ( H == 0.0 )
664 L = aVector->back();
665
666 return std::pair < G4double , G4double > ( L , H );
667}
668
669
670
671G4double G4NeutronHPThermalScattering::get_linear_interpolated ( G4double x , std::pair< G4double , G4double > Low , std::pair< G4double , G4double > High )
672{
673 G4double y=0.0;
674 if ( High.first - Low.first != 0 )
675 y = ( High.second - Low.second ) / ( High.first - Low.first ) * ( x - Low.first ) + Low.second;
676 else
677 std::cout << "G4NeutronHPThermalScattering liner interpolation err!!" << std::endl;
678
679 return y;
680}
681
682
683
684E_isoAng G4NeutronHPThermalScattering::create_E_isoAng_from_energy ( G4double energy , std::vector< E_isoAng* >* vEPM )
685{
686 E_isoAng anEPM_T_E;
687
688 std::vector< E_isoAng* >::iterator iv;
689
690 std::vector< G4double > v_e;
691 v_e.clear();
692 for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
693 v_e.push_back ( (*iv)->energy );
694
695 std::pair < G4double , G4double > energyLH = find_LH ( energy , &v_e );
696 //std::cout << " " << energy/eV << " " << energyLH.first/eV << " " << energyLH.second/eV << std::endl;
697
698 E_isoAng* panEPM_T_EL=0;
699 E_isoAng* panEPM_T_EH=0;
700
701 if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
702 {
703 for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
704 {
705 if ( energyLH.first == (*iv)->energy )
706 break;
707 }
708 panEPM_T_EL = *iv;
709 iv++;
710 panEPM_T_EH = *iv;
711 }
712 else if ( energyLH.first == 0.0 )
713 {
714 panEPM_T_EL = (*vEPM)[0];
715 panEPM_T_EH = (*vEPM)[1];
716 }
717 else if ( energyLH.second == 0.0 )
718 {
719 panEPM_T_EH = (*vEPM).back();
720 iv = vEPM->end();
721 iv--;
722 iv--;
723 panEPM_T_EL = *iv;
724 }
725
726 if ( panEPM_T_EL->n == panEPM_T_EH->n )
727 {
728 anEPM_T_E.energy = energy;
729 anEPM_T_E.n = panEPM_T_EL->n;
730
731 for ( G4int i=0 ; i < panEPM_T_EL->n ; i++ )
732 {
733 G4double angle;
734 angle = get_linear_interpolated ( energy , std::pair< G4double , G4double > ( energyLH.first , panEPM_T_EL->isoAngle[ i ] ) , std::pair< G4double , G4double > ( energyLH.second , panEPM_T_EH->isoAngle[ i ] ) );
735 anEPM_T_E.isoAngle.push_back( angle );
736 }
737 }
738 else
739 {
740 G4cout << "G4NeutronHPThermalScattering Do not Suuport yet." << G4endl;
741 }
742
743
744 return anEPM_T_E;
745}
746
747
748
749G4double G4NeutronHPThermalScattering::get_secondary_energy_from_E_P_E_isoAng ( G4double random , E_P_E_isoAng* anE_P_E_isoAng )
750{
751
752 G4double secondary_energy = 0.0;
753
754 G4int n = anE_P_E_isoAng->n;
755 G4double sum_p = 0.0; // sum_p_H
756 G4double sum_p_L = 0.0;
757
758 G4double total=0.0;
759
760/*
761 delete for speed up
762 for ( G4int i = 0 ; i < n-1 ; i++ )
763 {
764 G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
765 G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
766 G4double dE = E_H - E_L;
767 total += ( ( anE_P_E_isoAng->prob[i] ) * dE );
768 }
769
770 if ( std::abs( total - anE_P_E_isoAng->sum_of_probXdEs ) > 1.0e-14 ) std::cout << total - anE_P_E_isoAng->sum_of_probXdEs << std::endl;
771*/
772 total = anE_P_E_isoAng->sum_of_probXdEs;
773
774 for ( G4int i = 0 ; i < n-1 ; i++ )
775 {
776 G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
777 G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
778 G4double dE = E_H - E_L;
779 sum_p += ( ( anE_P_E_isoAng->prob[i] ) * dE );
780
781 if ( random <= sum_p/total )
782 {
783 secondary_energy = get_linear_interpolated ( random , std::pair < G4double , G4double > ( sum_p_L/total , E_L ) , std::pair < G4double , G4double > ( sum_p/total , E_H ) );
784 secondary_energy = secondary_energy*eV; //need eV
785 break;
786 }
787 sum_p_L = sum_p;
788 }
789
790 return secondary_energy;
791}
792
793
794
795std::pair< G4double , E_isoAng > G4NeutronHPThermalScattering::create_sE_and_EPM_from_pE_and_vE_P_E_isoAng ( G4double rand_for_sE , G4double pE , std::vector < E_P_E_isoAng* >* vNEP_EPM )
796{
797
798 std::map< G4double , G4int > map_energy;
799 map_energy.clear();
800 std::vector< G4double > v_energy;
801 v_energy.clear();
802 std::vector< E_P_E_isoAng* >::iterator itv;
803 G4int i = 0;
804 for ( itv = vNEP_EPM->begin(); itv != vNEP_EPM->end(); itv++ )
805 {
806 v_energy.push_back( (*itv)->energy );
807 map_energy.insert( std::pair < G4double , G4int > ( (*itv)->energy , i ) );
808 i++;
809 }
810
811 std::pair < G4double , G4double > energyLH = find_LH ( pE , &v_energy );
812
813 E_P_E_isoAng* pE_P_E_isoAng_EL = 0;
814 E_P_E_isoAng* pE_P_E_isoAng_EH = 0;
815
816 if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
817 {
818 pE_P_E_isoAng_EL = (*vNEP_EPM)[ map_energy.find ( energyLH.first )->second ];
819 pE_P_E_isoAng_EH = (*vNEP_EPM)[ map_energy.find ( energyLH.second )->second ];
820 }
821 else if ( energyLH.first == 0.0 )
822 {
823 pE_P_E_isoAng_EL = (*vNEP_EPM)[ 0 ];
824 pE_P_E_isoAng_EH = (*vNEP_EPM)[ 1 ];
825 }
826 if ( energyLH.second == 0.0 )
827 {
828 pE_P_E_isoAng_EH = (*vNEP_EPM).back();
829 itv = vNEP_EPM->end();
830 itv--;
831 itv--;
832 pE_P_E_isoAng_EL = *itv;
833 }
834
835
836 G4double sE;
837 G4double sE_L;
838 G4double sE_H;
839
840
841 sE_L = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EL );
842 sE_H = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EH );
843
844 sE = get_linear_interpolated ( pE , std::pair < G4double , G4double > ( energyLH.first , sE_L ) , std::pair < G4double , G4double > ( energyLH.second , sE_H ) );
845
846
847 E_isoAng E_isoAng_L = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EL->vE_isoAngle) );
848 E_isoAng E_isoAng_H = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EH->vE_isoAngle) );
849
850 E_isoAng anE_isoAng;
851 if ( E_isoAng_L.n == E_isoAng_H.n )
852 {
853 anE_isoAng.n = E_isoAng_L.n;
854 for ( G4int j=0 ; j < anE_isoAng.n ; j++ )
855 {
856 G4double angle;
857 angle = get_linear_interpolated ( sE , std::pair< G4double , G4double > ( sE_L , E_isoAng_L.isoAngle[ j ] ) , std::pair< G4double , G4double > ( sE_H , E_isoAng_H.isoAngle[ j ] ) );
858 anE_isoAng.isoAngle.push_back( angle );
859 }
860 }
861 else
862 {
863 std::cout << "Do not Suuport yet." << std::endl;
864 }
865
866
867
868 return std::pair< G4double , E_isoAng >( sE , anE_isoAng);
869}
870
871void G4NeutronHPThermalScattering::buildPhysicsTable()
872{
873
874 dic.clear();
875 std::map < G4String , G4int > co_dic;
876
877 //Searching Nist Materials
878 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
879 size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
880 for ( size_t i = 0 ; i < numberOfMaterials ; i++ )
881 {
882 G4Material* material = (*theMaterialTable)[i];
883 size_t numberOfElements = material->GetNumberOfElements();
884 for ( size_t j = 0 ; j < numberOfElements ; j++ )
885 {
886 const G4Element* element = material->GetElement(j);
887 if ( names.IsThisThermalElement ( material->GetName() , element->GetName() ) )
888 {
889 G4int ts_ID_of_this_geometry;
890 G4String ts_ndl_name = names.GetTS_NDL_Name( material->GetName() , element->GetName() );
891 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
892 {
893 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
894 }
895 else
896 {
897 ts_ID_of_this_geometry = co_dic.size();
898 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
899 }
900
901 //G4cout << "Neutron HP Thermal Scattering: Registering a material-element pair of "
902 // << material->GetName() << " " << element->GetName()
903 // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
904
905 dic.insert( std::pair < std::pair < G4Material* , const G4Element* > , G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
906 }
907 }
908 }
909
910 //Searching TS Elements
911 static const G4ElementTable* theElementTable = G4Element::GetElementTable();
912 size_t numberOfElements = G4Element::GetNumberOfElements();
913 //size_t numberOfThermalElements = 0;
914 for ( size_t i = 0 ; i < numberOfElements ; i++ )
915 {
916 const G4Element* element = (*theElementTable)[i];
917 if ( names.IsThisThermalElement ( element->GetName() ) )
918 {
919 if ( names.IsThisThermalElement ( element->GetName() ) )
920 {
921 G4int ts_ID_of_this_geometry;
922 G4String ts_ndl_name = names.GetTS_NDL_Name( element->GetName() );
923 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
924 {
925 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
926 }
927 else
928 {
929 ts_ID_of_this_geometry = co_dic.size();
930 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
931 }
932
933 //G4cout << "Neutron HP Thermal Scattering: Registering an element of "
934 // << material->GetName() << " " << element->GetName()
935 // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
936
937 dic.insert( std::pair < std::pair < const G4Material* , const G4Element* > , G4int > ( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) , ts_ID_of_this_geometry ) );
938 }
939 }
940 }
941
942 G4cout << G4endl;
943 G4cout << "Neutron HP Thermal Scattering: Following material-element pairs or elements are registered." << G4endl;
944 for ( std::map < std::pair < const G4Material* , const G4Element* > , G4int >::iterator it = dic.begin() ; it != dic.end() ; it++ )
945 {
946 if ( it->first.first != NULL )
947 {
948 G4cout << "Material " << it->first.first->GetName() << " - Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
949 }
950 else
951 {
952 G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
953 }
954 }
955 G4cout << G4endl;
956
957 // Read Cross Section Data files
958
959 G4String dirName;
960 if ( !getenv( "G4NEUTRONHPDATA" ) )
961 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
962 dirName = getenv( "G4NEUTRONHPDATA" );
963
964 //dirName = baseName + "/ThermalScattering";
965
966 G4String name;
967
968 for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
969 {
970 G4String tsndlName = it->first;
971 G4int ts_ID = it->second;
972
973 // Coherent
974 G4String fsName = "/ThermalScattering/Coherent/FS/";
975 G4String fileName = dirName + fsName + tsndlName;
976 coherentFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* > ( ts_ID , readACoherentFSDATA( fileName ) ) );
977
978 // incoherent elastic
979 fsName = "/ThermalScattering/Incoherent/FS/";
980 fileName = dirName + fsName + tsndlName;
981 incoherentFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < E_isoAng* >* >* > ( ts_ID , readAnIncoherentFSDATA( fileName ) ) );
982
983 // inelastic
984 fsName = "/ThermalScattering/Inelastic/FS/";
985 fileName = dirName + fsName + tsndlName;
986 inelasticFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* > ( ts_ID , readAnInelasticFSDATA( fileName ) ) );
987 }
988}
989
990
991G4int G4NeutronHPThermalScattering::getTS_ID ( const G4Material* material , const G4Element* element )
992{
993 G4int result = -1;
994 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
995 result = dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
996 return result;
997}
998
999const std::pair<G4double, G4double> G4NeutronHPThermalScattering::GetFatalEnergyCheckLevels() const
1000{
1001 //return std::pair<G4double, G4double>(10*perCent,10*GeV);
1002 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
1003}
G4ElementTable
std::vector< G4Element * > G4ElementTable
Definition: G4ElementTable.hh:43
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:43
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Element::GetNumberOfElements
static size_t GetNumberOfElements()
Definition: G4Element.cc:406
G4Element::GetName
const G4String & GetName() const
Definition: G4Element.hh:127
G4Element::GetElementTable
static const G4ElementTable * GetElementTable()
Definition: G4Element.cc:399
G4HadFinalState::SetEnergyChange
void SetEnergyChange(G4double anEnergy)
Definition: G4HadFinalState.cc:42
G4HadFinalState::SetMomentumChange
void SetMomentumChange(const G4ThreeVector &aV)
Definition: G4HadFinalState.cc:54
G4HadProjectile::GetMaterial
const G4Material * GetMaterial() const
Definition: G4HadProjectile.hh:76
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:86
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:90
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:103
G4Material::GetMaterialTable
static const G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:562
G4Material::GetNumberOfMaterials
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:569
G4Material::GetTemperature
G4double GetTemperature() const
Definition: G4Material.hh:181
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:201
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:177
G4NeutronHPThermalScatteringData::GetInelasticCrossSection
G4double GetInelasticCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4NeutronHPThermalScatteringData.cc:441
G4NeutronHPThermalScatteringData::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4NeutronHPThermalScatteringData.cc:204
G4NeutronHPThermalScatteringData::GetCrossSection
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4NeutronHPThermalScatteringData.cc:419
G4NeutronHPThermalScatteringData::GetCoherentCrossSection
G4double GetCoherentCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4NeutronHPThermalScatteringData.cc:450
G4NeutronHPThermalScatteringNames::GetTS_NDL_Name
G4String GetTS_NDL_Name(G4String nameG4Element)
Definition: G4NeutronHPThermalScatteringNames.hh:55
G4NeutronHPThermalScattering::GetFatalEnergyCheckLevels
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
Definition: G4NeutronHPThermalScattering.cc:999
G4NeutronHPThermalScattering::ApplyYourself
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4NeutronHPThermalScattering.cc:302
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4String::first
G4int first(char) const
E_P_E_isoAng
Definition: G4NeutronHPThermalScattering.hh:56
E_P_E_isoAng::prob
std::vector< G4double > prob
Definition: G4NeutronHPThermalScattering.hh:59
E_P_E_isoAng::vE_isoAngle
std::vector< E_isoAng * > vE_isoAngle
Definition: G4NeutronHPThermalScattering.hh:60
E_P_E_isoAng::n
G4int n
Definition: G4NeutronHPThermalScattering.hh:58
E_P_E_isoAng::sum_of_probXdEs
G4double sum_of_probXdEs
Definition: G4NeutronHPThermalScattering.hh:61
E_P_E_isoAng::energy
G4double energy
Definition: G4NeutronHPThermalScattering.hh:57
E_isoAng
Definition: G4NeutronHPThermalScattering.hh:49
E_isoAng::isoAngle
std::vector< G4double > isoAngle
Definition: G4NeutronHPThermalScattering.hh:52
E_isoAng::energy
G4double energy
Definition: G4NeutronHPThermalScattering.hh:50
E_isoAng::n
G4int n
Definition: G4NeutronHPThermalScattering.hh:51
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83