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

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

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