ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPThermalScattering.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPThermalScattering.cc
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 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
44 //
48 #include "G4ParticleHPElastic.hh"
49 #include "G4ParticleHPManager.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4Neutron.hh"
52 #include "G4ElementTable.hh"
53 #include "G4MaterialTable.hh"
54 #include "G4Threading.hh"
55 
57  :G4HadronicInteraction("NeutronHPThermalScattering")
58 ,coherentFSs(NULL)
59 ,incoherentFSs(NULL)
60 ,inelasticFSs(NULL)
61 {
63 
64  SetMinEnergy( 0.*eV );
65  SetMaxEnergy( 4*eV );
67 
68  //sizeOfMaterialTable = G4Material::GetMaterialTable()->size();
69  //buildPhysicsTable();
70  nMaterial = 0;
71  nElement = 0;
72 }
73 
74 
75 
77 {
78 
79 /*
80  for ( std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >::iterator it = incoherentFSs->begin() ; it != incoherentFSs->end() ; it++ )
81  {
82  std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
83  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
84  {
85  std::vector< E_isoAng* >::iterator ittt;
86  for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
87  {
88  delete *ittt;
89  }
90  delete itt->second;
91  }
92  delete it->second;
93  }
94 
95  for ( std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it = coherentFSs->begin() ; it != coherentFSs->end() ; it++ )
96  {
97  std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt;
98  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
99  {
100  std::vector < std::pair< G4double , G4double >* >::iterator ittt;
101  for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
102  {
103  delete *ittt;
104  }
105  delete itt->second;
106  }
107  delete it->second;
108  }
109 
110  for ( std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >::iterator it = inelasticFSs->begin() ; it != inelasticFSs->end() ; it++ )
111  {
112  std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt;
113  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
114  {
115  std::vector < E_P_E_isoAng* >::iterator ittt;
116  for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
117  {
118  std::vector < E_isoAng* >::iterator it4;
119  for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ )
120  {
121  delete *it4;
122  }
123  delete *ittt;
124  }
125  delete itt->second;
126  }
127  delete it->second;
128  }
129 */
130 
131  delete theHPElastic;
132  //TKDB 160506
133  //delete theXSection;
134 }
135 
137 
138 if ( incoherentFSs != NULL ) {
139  for ( std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >::iterator it = incoherentFSs->begin() ; it != incoherentFSs->end() ; it++ )
140  {
141  std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
142  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
143  {
144  std::vector< E_isoAng* >::iterator ittt;
145  for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
146  {
147  delete *ittt;
148  }
149  delete itt->second;
150  }
151  delete it->second;
152  }
153 }
154 
155 if ( coherentFSs != NULL ) {
156  for ( std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it = coherentFSs->begin() ; it != coherentFSs->end() ; it++ )
157  {
158  std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt;
159  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
160  {
161  std::vector < std::pair< G4double , G4double >* >::iterator ittt;
162  for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
163  {
164  delete *ittt;
165  }
166  delete itt->second;
167  }
168  delete it->second;
169  }
170 }
171 
172 if ( inelasticFSs != NULL ) {
173  for ( std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >::iterator it = inelasticFSs->begin() ; it != inelasticFSs->end() ; it++ )
174  {
175  std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt;
176  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
177  {
178  std::vector < E_P_E_isoAng* >::iterator ittt;
179  for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
180  {
181  std::vector < E_isoAng* >::iterator it4;
182  for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ )
183  {
184  delete *it4;
185  }
186  delete *ittt;
187  }
188  delete itt->second;
189  }
190  delete it->second;
191  }
192 }
193 
194  incoherentFSs = NULL;
195  coherentFSs = NULL;
196  inelasticFSs = NULL;
197 
198 }
199 
200 
201 
204  theHPElastic->BuildPhysicsTable( particle );
205 }
206 
207 
208 
209 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* G4ParticleHPThermalScattering::readACoherentFSDATA( G4String name )
210 {
211 
212  std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* aCoherentFSDATA = new std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >;
213 
214  //std::ifstream theChannel( name.c_str() );
215  std::istringstream theChannel(std::ios::in);
217 
218  std::vector< G4double > vBraggE;
219 
220  G4int dummy;
221  while ( theChannel >> dummy ) // MF // Loop checking, 11.05.2015, T. Koi
222  {
223  theChannel >> dummy; // MT
224  G4double temp;
225  theChannel >> temp;
226  std::vector < std::pair< G4double , G4double >* >* anBragE_P = new std::vector < std::pair< G4double , G4double >* >;
227 
228  G4int n;
229  theChannel >> n;
230  for ( G4int i = 0 ; i < n ; i++ )
231  {
232  G4double Ei;
233  G4double Pi;
234  if ( aCoherentFSDATA->size() == 0 )
235  {
236  theChannel >> Ei;
237  vBraggE.push_back( Ei );
238  }
239  else
240  {
241  Ei = vBraggE[ i ];
242  }
243  theChannel >> Pi;
244  anBragE_P->push_back ( new std::pair < G4double , G4double > ( Ei , Pi ) );
245  //G4cout << "Coherent Elastic " << Ei << " " << Pi << G4endl;
246  }
247  aCoherentFSDATA->insert ( std::pair < G4double , std::vector < std::pair< G4double , G4double >* >* > ( temp , anBragE_P ) );
248  }
249 
250  return aCoherentFSDATA;
251 }
252 
253 
254 
255 std::map < G4double , std::vector < E_P_E_isoAng* >* >* G4ParticleHPThermalScattering::readAnInelasticFSDATA ( G4String name )
256 {
257  std::map < G4double , std::vector < E_P_E_isoAng* >* >* anT_E_P_E_isoAng = new std::map < G4double , std::vector < E_P_E_isoAng* >* >;
258 
259  //std::ifstream theChannel( name.c_str() );
260  std::istringstream theChannel(std::ios::in);
262 
263  G4int dummy;
264  while ( theChannel >> dummy ) // MF // Loop checking, 11.05.2015, T. Koi
265  {
266  theChannel >> dummy; // MT
267  G4double temp;
268  theChannel >> temp;
269  std::vector < E_P_E_isoAng* >* vE_P_E_isoAng = new std::vector < E_P_E_isoAng* >;
270  G4int n;
271  theChannel >> n;
272  for ( G4int i = 0 ; i < n ; i++ )
273  {
274  vE_P_E_isoAng->push_back ( readAnE_P_E_isoAng ( &theChannel ) );
275  }
276  anT_E_P_E_isoAng->insert ( std::pair < G4double , std::vector < E_P_E_isoAng* >* > ( temp , vE_P_E_isoAng ) );
277  }
278  //theChannel.close();
279 
280  return anT_E_P_E_isoAng;
281 }
282 
283 
284 
286 {
287  E_P_E_isoAng* aData = new E_P_E_isoAng;
288 
289  G4double dummy;
291  G4int nep , nl;
292  *file >> dummy;
293  *file >> energy;
294  aData->energy = energy*eV;
295  *file >> dummy;
296  *file >> dummy;
297  *file >> nep;
298  *file >> nl;
299  aData->n = nep/nl;
300  for ( G4int i = 0 ; i < aData->n ; i++ )
301  {
302  G4double prob;
303  E_isoAng* anE_isoAng = new E_isoAng;
304  aData->vE_isoAngle.push_back( anE_isoAng );
305  *file >> energy;
306  anE_isoAng->energy = energy*eV;
307  anE_isoAng->n = nl - 2;
308  anE_isoAng->isoAngle.resize( anE_isoAng->n );
309  *file >> prob;
310  aData->prob.push_back( prob );
311  //G4cout << "G4ParticleHPThermalScattering inelastic " << energy/eV << " " << i << " " << prob << " " << aData->prob[ i ] << G4endl;
312  for ( G4int j = 0 ; j < anE_isoAng->n ; j++ )
313  {
314  G4double x;
315  *file >> x;
316  anE_isoAng->isoAngle[j] = x ;
317  //G4cout << "G4ParticleHPThermalScattering inelastic " << x << anE_isoAng->isoAngle[j] << G4endl;
318  }
319  }
320 
321  // Calcuate sum_of_provXdEs
322  G4double total = 0;
323  for ( G4int i = 0 ; i < aData->n - 1 ; i++ )
324  {
325  G4double E_L = aData->vE_isoAngle[i]->energy/eV;
326  G4double E_H = aData->vE_isoAngle[i+1]->energy/eV;
327  G4double dE = E_H - E_L;
328  total += ( ( aData->prob[i] ) * dE );
329  }
330  aData->sum_of_probXdEs = total;
331 
332  return aData;
333 }
334 
335 
336 
337 std::map < G4double , std::vector < E_isoAng* >* >* G4ParticleHPThermalScattering::readAnIncoherentFSDATA ( G4String name )
338 {
339  std::map < G4double , std::vector < E_isoAng* >* >* T_E = new std::map < G4double , std::vector < E_isoAng* >* >;
340 
341  //std::ifstream theChannel( name.c_str() );
342  std::istringstream theChannel(std::ios::in);
344 
345  G4int dummy;
346  while ( theChannel >> dummy ) // MF // Loop checking, 11.05.2015, T. Koi
347  {
348  theChannel >> dummy; // MT
349  G4double temp;
350  theChannel >> temp;
351  std::vector < E_isoAng* >* vE_isoAng = new std::vector < E_isoAng* >;
352  G4int n;
353  theChannel >> n;
354  for ( G4int i = 0 ; i < n ; i++ )
355  vE_isoAng->push_back ( readAnE_isoAng( &theChannel ) );
356  T_E->insert ( std::pair < G4double , std::vector < E_isoAng* >* > ( temp , vE_isoAng ) );
357  }
358  //theChannel.close();
359 
360  return T_E;
361 }
362 
363 
364 
366 {
367  E_isoAng* aData = new E_isoAng;
368 
369  G4double dummy;
371  G4int n;
372  *file >> dummy;
373  *file >> energy;
374  *file >> dummy;
375  *file >> dummy;
376  *file >> n;
377  *file >> dummy;
378  aData->energy = energy*eV;
379  aData->n = n-2;
380  aData->isoAngle.resize( n );
381 
382  *file >> dummy;
383  *file >> dummy;
384  for ( G4int i = 0 ; i < aData->n ; i++ )
385  *file >> aData->isoAngle[i];
386 
387  return aData;
388 }
389 
390 
391 
393 {
394 
395 /*
396  //Trick for dynamically generated materials
397  if ( sizeOfMaterialTable != G4Material::GetMaterialTable()->size() ) {
398  sizeOfMaterialTable = G4Material::GetMaterialTable()->size();
399  buildPhysicsTable();
400  theXSection->BuildPhysicsTable( *aTrack.GetDefinition() );
401  }
402 */
403 // Select Element > Reaction >
404 
405  const G4Material * theMaterial = aTrack.GetMaterial();
406  G4double aTemp = theMaterial->GetTemperature();
407  G4int n = theMaterial->GetNumberOfElements();
408  //static const G4ElementTable* theElementTable = G4Element::GetElementTable();
409 
410  G4bool findThermalElement = false;
411  G4int ielement;
412  const G4Element* theElement = NULL;
413  for ( G4int i = 0; i < n ; i++ )
414  {
415  theElement = theMaterial->GetElement(i);
416  //Select target element
417  if ( aNucleus.GetZ_asInt() == (G4int)(theElement->GetZ() + 0.5 ) )
418  {
419  //Check Applicability of Thermal Scattering
420  if ( getTS_ID( NULL , theElement ) != -1 )
421  {
422  ielement = getTS_ID( NULL , theElement );
423  findThermalElement = true;
424  break;
425  }
426  else if ( getTS_ID( theMaterial , theElement ) != -1 )
427  {
428  ielement = getTS_ID( theMaterial , theElement );
429  findThermalElement = true;
430  break;
431  }
432  }
433  }
434 
435  if ( findThermalElement == true )
436  {
437 
438 // Select Reaction (Inelastic, coherent, incoherent)
439 
440  const G4ParticleDefinition* pd = aTrack.GetDefinition();
441  G4DynamicParticle* dp = new G4DynamicParticle ( pd , aTrack.Get4Momentum() );
442  G4double total = theXSection->GetCrossSection( dp , theElement , theMaterial );
443  G4double inelastic = theXSection->GetInelasticCrossSection( dp , theElement , theMaterial );
444 
445 
446  G4double random = G4UniformRand();
447  if ( random <= inelastic/total )
448  {
449  // Inelastic
450 
451  // T_L and T_H
452  std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator it;
453  std::vector<G4double> v_temp;
454  v_temp.clear();
455  for ( it = inelasticFSs->find( ielement )->second->begin() ; it != inelasticFSs->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 // For T_L aNEP_EPM_TL and T_H aNEP_EPM_TH
464 //
465  std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = 0;
466  std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = 0;
467 
468  if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
469  {
470  vNEP_EPM_TL = inelasticFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second;
471  vNEP_EPM_TH = inelasticFSs->find( ielement )->second->find ( tempLH.second/kelvin )->second;
472  }
473  else if ( tempLH.first == 0.0 )
474  {
475  std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
476  itm = inelasticFSs->find( ielement )->second->begin();
477  vNEP_EPM_TL = itm->second;
478  itm++;
479  vNEP_EPM_TH = itm->second;
480  tempLH.first = tempLH.second;
481  tempLH.second = itm->first;
482  }
483  else if ( tempLH.second == 0.0 )
484  {
485  std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
486  itm = inelasticFSs->find( ielement )->second->end();
487  itm--;
488  vNEP_EPM_TH = itm->second;
489  itm--;
490  vNEP_EPM_TL = itm->second;
491  tempLH.second = tempLH.first;
492  tempLH.first = itm->first;
493  }
494 
495  G4double rand_for_sE = G4UniformRand();
496 
497  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 );
498  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 );
499 
500  G4double sE;
501  sE = get_linear_interpolated ( aTemp , std::pair < G4double , G4double > ( tempLH.first , TL.first ) , std::pair < G4double , G4double > ( tempLH.second , TH.first ) );
502 
503  G4double mu=1.0;
504  E_isoAng anE_isoAng;
505  if ( TL.second.n == TH.second.n )
506  {
507  anE_isoAng.energy = sE;
508  anE_isoAng.n = TL.second.n;
509  for ( G4int i=0 ; i < anE_isoAng.n ; i++ )
510  {
511  G4double angle;
512  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 ] ) );
513  anE_isoAng.isoAngle.push_back( angle );
514  }
515  mu = getMu( &anE_isoAng );
516 
517  } else {
518  //TL.second.n != TH.second.n
519  throw G4HadronicException(__FILE__, __LINE__, "A problem is found in Thermal Scattering Data! Do not yet supported");
520  }
521 
522  //set
524  theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
525 
526  }
527  //else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , (*theElementTable)[ ielement ] , aTemp ) ) / total )
528  else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , theElement , theMaterial ) ) / total )
529  {
530  // Coherent Elastic
531 
532  G4double E = aTrack.GetKineticEnergy();
533 
534  // T_L and T_H
535  std::map < G4double , std::vector< std::pair< G4double , G4double >* >* >::iterator it;
536  std::vector<G4double> v_temp;
537  v_temp.clear();
538  for ( it = coherentFSs->find( ielement )->second->begin() ; it != coherentFSs->find( ielement )->second->end() ; it++ )
539  {
540  v_temp.push_back( it->first );
541  }
542 
543 // T_L T_H
544  std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
545 //
546 //
547 // For T_L anEPM_TL and T_H anEPM_TH
548 //
549  std::vector< std::pair< G4double , G4double >* >* pvE_p_TL = NULL;
550  std::vector< std::pair< G4double , G4double >* >* pvE_p_TH = NULL;
551 
552  if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
553  {
554  pvE_p_TL = coherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second;
555  pvE_p_TH = coherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second;
556  }
557  else if ( tempLH.first == 0.0 )
558  {
559  pvE_p_TL = coherentFSs->find( ielement )->second->find ( v_temp[ 0 ] )->second;
560  pvE_p_TH = coherentFSs->find( ielement )->second->find ( v_temp[ 1 ] )->second;
561  tempLH.first = tempLH.second;
562  tempLH.second = v_temp[ 1 ];
563  }
564  else if ( tempLH.second == 0.0 )
565  {
566  pvE_p_TH = coherentFSs->find( ielement )->second->find ( v_temp.back() )->second;
567  std::vector< G4double >::iterator itv;
568  itv = v_temp.end();
569  itv--;
570  itv--;
571  pvE_p_TL = coherentFSs->find( ielement )->second->find ( *itv )->second;
572  tempLH.second = tempLH.first;
573  tempLH.first = *itv;
574  }
575  else
576  {
577  //tempLH.first == 0.0 && tempLH.second
578  throw G4HadronicException(__FILE__, __LINE__, "A problem is found in Thermal Scattering Data! Unexpected temperature values in data");
579  }
580 
581  std::vector< G4double > vE_T;
582  std::vector< G4double > vp_T;
583 
584  G4int n1 = pvE_p_TL->size();
585  //G4int n2 = pvE_p_TH->size();
586 
587  //171005 fix bug, contribution from H.N. TRAN@CEA
588  for ( G4int i=0 ; i < n1 ; i++ )
589  {
590  if ( (*pvE_p_TL)[i]->first != (*pvE_p_TH)[i]->first ) throw G4HadronicException(__FILE__, __LINE__, "A problem is found in Thermal Scattering Data!");
591  vE_T.push_back ( (*pvE_p_TL)[i]->first );
592  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 ) ) );
593  }
594 
595  G4int j = 0;
596  for ( G4int i = 1 ; i < n1 ; i++ )
597  {
598  if ( E/eV < vE_T[ i ] )
599  {
600  j = i-1;
601  break;
602  }
603  }
604 
605  G4double rand_for_mu = G4UniformRand();
606 
607  G4int k = 0;
608  for ( G4int i = 0 ; i <= j ; i++ )
609  {
610  G4double Pi = vp_T[ i ] / vp_T[ j ];
611  if ( rand_for_mu < Pi )
612  {
613  k = i;
614  break;
615  }
616  }
617 
618  G4double Ei = vE_T[ k ];
619 
620  G4double mu = 1 - 2 * Ei / (E/eV) ;
621  //111102
622  if ( mu < -1.0 ) mu = -1.0;
623  //G4cout << "E= " << E/eV << ", Ei= " << Ei << ", mu= " << mu << G4endl;
624 
626  theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
627 
628  }
629  else
630  {
631  // InCoherent Elastic
632 
633  // T_L and T_H
634  std::map < G4double , std::vector < E_isoAng* >* >::iterator it;
635  std::vector<G4double> v_temp;
636  v_temp.clear();
637  for ( it = incoherentFSs->find( ielement )->second->begin() ; it != incoherentFSs->find( ielement )->second->end() ; it++ )
638  {
639  v_temp.push_back( it->first );
640  }
641 
642 // T_L T_H
643  std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
644 
645 //
646 // For T_L anEPM_TL and T_H anEPM_TH
647 //
648 
649  E_isoAng anEPM_TL_E;
650  E_isoAng anEPM_TH_E;
651 
652  if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) {
653  //Interpolate TL and TH
654  anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second );
655  anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( tempLH.second/kelvin )->second );
656  } else if ( tempLH.first == 0.0 ) {
657  //Extrapolate T0 and T1
658  anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( v_temp[ 0 ] )->second );
659  anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( v_temp[ 1 ] )->second );
660  tempLH.first = tempLH.second;
661  tempLH.second = v_temp[ 1 ];
662  } else if ( tempLH.second == 0.0 ) {
663  //Extrapolate Tmax-1 and Tmax
664  anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( v_temp.back() )->second );
665  std::vector< G4double >::iterator itv;
666  itv = v_temp.end();
667  itv--;
668  itv--;
669  anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( *itv )->second );
670  tempLH.second = tempLH.first;
671  tempLH.first = *itv;
672  }
673 
674  // E_isoAng for aTemp and aTrack.GetKineticEnergy()
675  G4double mu=1.0;
676  E_isoAng anEPM_T_E;
677 
678  if ( anEPM_TL_E.n == anEPM_TH_E.n )
679  {
680  anEPM_T_E.n = anEPM_TL_E.n;
681  for ( G4int i=0 ; i < anEPM_TL_E.n ; i++ )
682  {
683  G4double angle;
684  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 ] ) );
685  anEPM_T_E.isoAngle.push_back( angle );
686  }
687  mu = getMu ( &anEPM_T_E );
688 
689  } else {
690  // anEPM_TL_E.n != anEPM_TH_E.n
691  throw G4HadronicException(__FILE__, __LINE__, "A problem is found in Thermal Scattering Data! Do not yet supported");
692  }
693 
694  // Set Final State
695  theParticleChange.SetEnergyChange( aTrack.GetKineticEnergy() ); // No energy change in Elastic
696  theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
697 
698  }
699  delete dp;
700 
701  return &theParticleChange;
702 
703  }
704  else
705  {
706  // Not thermal element
707  // Neutron HP will handle
708  return theHPElastic -> ApplyYourself( aTrack, aNucleus );
709  }
710 
711 }
712 
713 
714 
716 {
717 
718  G4double random = G4UniformRand();
719  G4double result = 0.0;
720 
721  G4int in = int ( random * ( (*anEPM).n ) );
722 
723  if ( in != 0 )
724  {
725  G4double mu_l = (*anEPM).isoAngle[ in-1 ];
726  G4double mu_h = (*anEPM).isoAngle[ in ];
727  result = ( mu_h - mu_l ) * ( random * ( (*anEPM).n ) - in ) + mu_l;
728  }
729  else
730  {
731  G4double x = random * (*anEPM).n;
732  //Bugzilla 1971
733  G4double ratio = 0.5;
735  if ( x <= ratio )
736  {
737  G4double mu_l = -1;
738  G4double mu_h = (*anEPM).isoAngle[ 0 ];
739  result = ( mu_h - mu_l ) * xx + mu_l;
740  }
741  else
742  {
743  G4double mu_l = (*anEPM).isoAngle[ (*anEPM).n - 1 ];
744  G4double mu_h = 1;
745  result = ( mu_h - mu_l ) * xx + mu_l;
746  }
747  }
748  return result;
749 }
750 
751 
752 
753 std::pair < G4double , G4double > G4ParticleHPThermalScattering::find_LH ( G4double x , std::vector< G4double >* aVector )
754 {
755  G4double LL = 0.0;
756  G4double H = 0.0;
757 
758  // v->size() == 1 --> LL=H=v(0)
759  if ( aVector->size() == 1 ) {
760  LL = aVector->front();
761  H = aVector->front();
762  } else {
763  // 1) temp < v(0) -> LL=0.0 H=v(0)
764  // 2) v(i-1) < temp <= v(i) -> LL=v(i-1) H=v(i)
765  // 3) v(imax) < temp -> LL=v(imax) H=0.0
766  for ( std::vector< G4double >::iterator
767  it = aVector->begin() ; it != aVector->end() ; it++ ) {
768  if ( x <= *it ) {
769  H = *it;
770  if ( it != aVector->begin() ) {
771  // 2)
772  it--;
773  LL = *it;
774  } else {
775  // 1)
776  LL = 0.0;
777  }
778  break;
779  }
780  }
781  // 3)
782  if ( H == 0.0 ) LL = aVector->back();
783  }
784 
785  return std::pair < G4double , G4double > ( LL , H );
786 }
787 
788 
789 
790 G4double G4ParticleHPThermalScattering::get_linear_interpolated ( G4double x , std::pair< G4double , G4double > Low , std::pair< G4double , G4double > High )
791 {
792  G4double y=0.0;
793  if ( High.first - Low.first != 0 ) {
794  y = ( High.second - Low.second ) / ( High.first - Low.first ) * ( x - Low.first ) + Low.second;
795  } else {
796  if ( High.second == Low.second ) {
797  y = High.second;
798  } else {
799  G4cout << "G4ParticleHPThermalScattering liner interpolation err!!" << G4endl;
800  }
801  }
802 
803  return y;
804 }
805 
806 
807 
808 E_isoAng
810  std::vector<E_isoAng*>* vEPM)
811 {
812  E_isoAng anEPM_T_E;
813  std::vector<E_isoAng*>::iterator iv;
814 
815  std::vector<G4double> v_e;
816  v_e.clear();
817  for (iv = vEPM->begin(); iv != vEPM->end(); iv++)
818  v_e.push_back( (*iv)->energy );
819 
820  std::pair<G4double, G4double> energyLH = find_LH(energy, &v_e);
821  //G4cout << " " << energy/eV << " " << energyLH.first/eV << " " << energyLH.second/eV << G4endl;
822 
823  E_isoAng* panEPM_T_EL = 0;
824  E_isoAng* panEPM_T_EH = 0;
825 
826  if (energyLH.first != 0.0 && energyLH.second != 0.0) {
827  for (iv = vEPM->begin(); iv != vEPM->end(); iv++) {
828  if (energyLH.first == (*iv)->energy) {
829  panEPM_T_EL = *iv;
830  iv++;
831  panEPM_T_EH = *iv;
832  break;
833  }
834  }
835 
836  } else if (energyLH.first == 0.0) {
837  panEPM_T_EL = (*vEPM)[0];
838  panEPM_T_EH = (*vEPM)[1];
839 
840  } else if (energyLH.second == 0.0) {
841  panEPM_T_EH = (*vEPM).back();
842  iv = vEPM->end();
843  iv--;
844  iv--;
845  panEPM_T_EL = *iv;
846  }
847 
848  if (panEPM_T_EL != 0 && panEPM_T_EH != 0) {
849  //checking isoAng has proper values or not
850  // Inelastic/FS, the first and last entries of *vEPM has all zero values.
851  if ( !(check_E_isoAng(panEPM_T_EL) ) ) panEPM_T_EL = panEPM_T_EH;
852  if ( !(check_E_isoAng(panEPM_T_EH) ) ) panEPM_T_EH = panEPM_T_EL;
853 
854  if (panEPM_T_EL->n == panEPM_T_EH->n) {
855  anEPM_T_E.energy = energy;
856  anEPM_T_E.n = panEPM_T_EL->n;
857 
858  for (G4int i=0; i < panEPM_T_EL->n; i++) {
859  G4double angle;
860  angle = get_linear_interpolated(energy, std::pair<G4double,G4double>(energyLH.first, panEPM_T_EL->isoAngle[i] ),
861  std::pair<G4double,G4double>(energyLH.second, panEPM_T_EH->isoAngle[i] ) );
862  anEPM_T_E.isoAngle.push_back(angle);
863  }
864 
865  } else {
866  G4Exception("G4ParticleHPThermalScattering::create_E_isoAng_from_energy",
867  "NotSupported", JustWarning,
868  "G4ParticleHPThermalScattering does not support yet EL->n != EH->n.");
869  }
870 
871  } else {
872  G4Exception("G4ParticleHPThermalScattering::create_E_isoAng_from_energy",
873  "HAD_THERM_000", FatalException,
874  "Pointer panEPM_T_EL or panEPM_T_EH is zero");
875  }
876 
877  return anEPM_T_E;
878 }
879 
880 
881 
883 {
884 
885  G4double secondary_energy = 0.0;
886 
887  G4int n = anE_P_E_isoAng->n;
888  G4double sum_p = 0.0; // sum_p_H
889  G4double sum_p_L = 0.0;
890 
891  G4double total=0.0;
892 
893 /*
894  delete for speed up
895  for ( G4int i = 0 ; i < n-1 ; i++ )
896  {
897  G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
898  G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
899  G4double dE = E_H - E_L;
900  total += ( ( anE_P_E_isoAng->prob[i] ) * dE );
901  }
902 
903  if ( std::abs( total - anE_P_E_isoAng->sum_of_probXdEs ) > 1.0e-14 ) G4cout << total - anE_P_E_isoAng->sum_of_probXdEs << G4endl;
904 */
905  total = anE_P_E_isoAng->sum_of_probXdEs;
906 
907  for ( G4int i = 0 ; i < n-1 ; i++ )
908  {
909  G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
910  G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
911  G4double dE = E_H - E_L;
912  sum_p += ( ( anE_P_E_isoAng->prob[i] ) * dE );
913 
914  if ( random <= sum_p/total )
915  {
916  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 ) );
917  secondary_energy = secondary_energy*eV; //need eV
918  break;
919  }
920  sum_p_L = sum_p;
921  }
922 
923  return secondary_energy;
924 }
925 
926 
927 
928 std::pair< G4double , E_isoAng > G4ParticleHPThermalScattering::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 )
929 {
930 
931  std::map< G4double , G4int > map_energy;
932  map_energy.clear();
933  std::vector< G4double > v_energy;
934  v_energy.clear();
935  std::vector< E_P_E_isoAng* >::iterator itv;
936  G4int i = 0;
937  for ( itv = vNEP_EPM->begin(); itv != vNEP_EPM->end(); itv++ )
938  {
939  v_energy.push_back( (*itv)->energy );
940  map_energy.insert( std::pair < G4double , G4int > ( (*itv)->energy , i ) );
941  i++;
942  }
943 
944  std::pair < G4double , G4double > energyLH = find_LH ( pE , &v_energy );
945 
946  E_P_E_isoAng* pE_P_E_isoAng_EL = 0;
947  E_P_E_isoAng* pE_P_E_isoAng_EH = 0;
948 
949  if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
950  {
951  pE_P_E_isoAng_EL = (*vNEP_EPM)[ map_energy.find ( energyLH.first )->second ];
952  pE_P_E_isoAng_EH = (*vNEP_EPM)[ map_energy.find ( energyLH.second )->second ];
953  }
954  else if ( energyLH.first == 0.0 )
955  {
956  pE_P_E_isoAng_EL = (*vNEP_EPM)[ 0 ];
957  pE_P_E_isoAng_EH = (*vNEP_EPM)[ 1 ];
958  }
959  if ( energyLH.second == 0.0 )
960  {
961  pE_P_E_isoAng_EH = (*vNEP_EPM).back();
962  itv = vNEP_EPM->end();
963  itv--;
964  itv--;
965  pE_P_E_isoAng_EL = *itv;
966  }
967 
968 
969  G4double sE;
970  G4double sE_L;
971  G4double sE_H;
972 
973 
974  sE_L = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EL );
975  sE_H = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EH );
976 
977  sE = get_linear_interpolated ( pE , std::pair < G4double , G4double > ( energyLH.first , sE_L ) , std::pair < G4double , G4double > ( energyLH.second , sE_H ) );
978 
979 
980  E_isoAng E_isoAng_L = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EL->vE_isoAngle) );
981  E_isoAng E_isoAng_H = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EH->vE_isoAngle) );
982 
983  E_isoAng anE_isoAng;
984  //For defeating warning message from compiler
985  anE_isoAng.n = 1;
986  anE_isoAng.energy = sE; //never used
987  if ( E_isoAng_L.n == E_isoAng_H.n )
988  {
989  anE_isoAng.n = E_isoAng_L.n;
990  for ( G4int j=0 ; j < anE_isoAng.n ; j++ )
991  {
992  G4double angle;
993  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 ] ) );
994  anE_isoAng.isoAngle.push_back( angle );
995  }
996  }
997  else
998  {
999  //G4cout << "Do not Suuport yet." << G4endl;
1000  throw G4HadronicException(__FILE__, __LINE__, "Unexpected values!");
1001  }
1002 
1003 
1004 
1005  return std::pair< G4double , E_isoAng >( sE , anE_isoAng);
1006 }
1007 
1009 {
1010 
1011  //Is rebuild of physics table a necessity
1012  if ( nMaterial == G4Material::GetMaterialTable()->size() && nElement == G4Element::GetElementTable()->size() ) {
1013  return;
1014  } else {
1016  nElement = G4Element::GetElementTable()->size();
1017  }
1018 
1019  dic.clear();
1020  std::map < G4String , G4int > co_dic;
1021 
1022  //Searching Nist Materials
1023  static G4ThreadLocal G4MaterialTable* theMaterialTable = 0 ; if (!theMaterialTable) theMaterialTable= G4Material::GetMaterialTable();
1024  size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
1025  for ( size_t i = 0 ; i < numberOfMaterials ; i++ )
1026  {
1027  G4Material* material = (*theMaterialTable)[i];
1028  size_t numberOfElements = material->GetNumberOfElements();
1029  for ( size_t j = 0 ; j < numberOfElements ; j++ )
1030  {
1031  const G4Element* element = material->GetElement(j);
1032  if ( names.IsThisThermalElement ( material->GetName() , element->GetName() ) )
1033  {
1034  G4int ts_ID_of_this_geometry;
1035  G4String ts_ndl_name = names.GetTS_NDL_Name( material->GetName() , element->GetName() );
1036  if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
1037  {
1038  ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
1039  }
1040  else
1041  {
1042  ts_ID_of_this_geometry = co_dic.size();
1043  co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
1044  }
1045 
1046  //G4cout << "Neutron HP Thermal Scattering: Registering a material-element pair of "
1047  // << material->GetName() << " " << element->GetName()
1048  // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
1049 
1050  dic.insert( std::pair < std::pair < G4Material* , const G4Element* > , G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
1051  }
1052  }
1053  }
1054 
1055  //Searching TS Elements
1056  static G4ThreadLocal G4ElementTable* theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
1057  size_t numberOfElements = G4Element::GetNumberOfElements();
1058  //size_t numberOfThermalElements = 0;
1059  for ( size_t i = 0 ; i < numberOfElements ; i++ )
1060  {
1061  const G4Element* element = (*theElementTable)[i];
1062  if ( names.IsThisThermalElement ( element->GetName() ) )
1063  {
1064  if ( names.IsThisThermalElement ( element->GetName() ) )
1065  {
1066  G4int ts_ID_of_this_geometry;
1067  G4String ts_ndl_name = names.GetTS_NDL_Name( element->GetName() );
1068  if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
1069  {
1070  ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
1071  }
1072  else
1073  {
1074  ts_ID_of_this_geometry = co_dic.size();
1075  co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
1076  }
1077 
1078  //G4cout << "Neutron HP Thermal Scattering: Registering an element of "
1079  // << material->GetName() << " " << element->GetName()
1080  // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
1081 
1082  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 ) );
1083  }
1084  }
1085  }
1086 
1087  G4cout << G4endl;
1088  G4cout << "Neutron HP Thermal Scattering: Following material-element pairs or elements are registered." << G4endl;
1089  for ( std::map < std::pair < const G4Material* , const G4Element* > , G4int >::iterator it = dic.begin() ; it != dic.end() ; it++ )
1090  {
1091  if ( it->first.first != NULL )
1092  {
1093  G4cout << "Material " << it->first.first->GetName() << " - Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
1094  }
1095  else
1096  {
1097  G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
1098  }
1099  }
1100  G4cout << G4endl;
1101 
1102  // Read Cross Section Data files
1103 
1108 
1109  if ( G4Threading::IsMasterThread() ) {
1110 
1112 
1113  if ( coherentFSs == NULL ) coherentFSs = new std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >;
1114  if ( incoherentFSs == NULL ) incoherentFSs = new std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >;
1115  if ( inelasticFSs == NULL ) inelasticFSs = new std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >;
1116 
1117  G4String dirName;
1118  if ( !std::getenv( "G4NEUTRONHPDATA" ) )
1119  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
1120  dirName = std::getenv( "G4NEUTRONHPDATA" );
1121 
1122  //G4String name;
1123 
1124  for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
1125  {
1126  G4String tsndlName = it->first;
1127  G4int ts_ID = it->second;
1128 
1129  // Coherent
1130  G4String fsName = "/ThermalScattering/Coherent/FS/";
1131  G4String fileName = dirName + fsName + tsndlName;
1132  coherentFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* > ( ts_ID , readACoherentFSDATA( fileName ) ) );
1133 
1134  // incoherent elastic
1135  fsName = "/ThermalScattering/Incoherent/FS/";
1136  fileName = dirName + fsName + tsndlName;
1137  incoherentFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < E_isoAng* >* >* > ( ts_ID , readAnIncoherentFSDATA( fileName ) ) );
1138 
1139  // inelastic
1140  fsName = "/ThermalScattering/Inelastic/FS/";
1141  fileName = dirName + fsName + tsndlName;
1142  inelasticFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* > ( ts_ID , readAnInelasticFSDATA( fileName ) ) );
1143  }
1144 
1148  }
1149 
1151 }
1152 
1153 
1155 {
1156  G4int result = -1;
1157  if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
1158  result = dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
1159  return result;
1160 }
1161 
1162 const std::pair<G4double, G4double> G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels() const
1163 {
1164  //return std::pair<G4double, G4double>(10*perCent,10*GeV);
1165  return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
1166 }
1167 
1169 {
1170  names.AddThermalElement( nameG4Element , filename );
1171  theXSection->AddUserThermalScatteringFile( nameG4Element , filename );
1173 }
1174 
1175 
1177 {
1178  G4bool result=false;
1179 
1180  G4int n = anE_IsoAng->n;
1181  G4double sum=0.0;
1182  for ( G4int i = 0 ; i < n ; i++ ) {
1183  sum += anE_IsoAng->isoAngle[ i ];
1184  }
1185  if ( sum != 0.0 ) result = true;
1186 
1187  return result;
1188 }
1189 
1190 void G4ParticleHPThermalScattering::ModelDescription(std::ostream& outFile) const
1191 {
1192  outFile << "High Precision model based on thermal scattering data in\n"
1193  << "evaluated nuclear data libraries for neutrons below 5eV\n"
1194  << "on specific materials\n";
1195 }