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

Last change on this file since 836 was 819, checked in by garnier, 16 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.