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

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

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