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