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

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

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