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