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

Last change on this file was 1340, checked in by garnier, 14 years ago

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