Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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