ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPThermalScatteringData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPThermalScatteringData.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 (SCCS/SLAC)
28 //
29 // Class Description
30 // Cross Sections for a high precision (based on evaluated data
31 // libraries) description of themal neutron scattering below 4 eV;
32 // Based on Thermal neutron scattering files
33 // from the evaluated nuclear data files ENDF/B-VI, Release2
34 // To be used in your physics list in case you need this physics.
35 // In this case you want to register an object of this class with
36 // the corresponding process.
37 // Class Description - End
38 
39 // 15-Nov-06 First implementation is done by T. Koi (SLAC/SCCS)
40 // 070625 implement clearCurrentXSData to fix memory leaking by T. Koi
41 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
42 //
43 
44 #include <list>
45 #include <algorithm>
46 
48 #include "G4ParticleHPManager.hh"
49 
50 #include "G4SystemOfUnits.hh"
51 #include "G4Neutron.hh"
52 #include "G4ElementTable.hh"
53 
54 #include "G4Threading.hh"
55 
57 :G4VCrossSectionDataSet("NeutronHPThermalScatteringData")
58 ,coherent(NULL)
59 ,incoherent(NULL)
60 ,inelastic(NULL)
61 {
62 // Upper limit of neutron energy
63  emax = 4*eV;
64  SetMinKinEnergy( 0*MeV );
66 
67  ke_cache = 0.0;
68  xs_cache = 0.0;
69  element_cache = NULL;
70  material_cache = NULL;
71 
72  indexOfThermalElement.clear();
73 
75 }
76 
78 {
79 
81 
82  delete names;
83 }
84 
86  G4int /*Z*/ , G4int /*A*/ ,
87  const G4Element* element ,
88  const G4Material* material )
89 {
90  G4double eKin = dp->GetKineticEnergy();
91  if ( eKin > 4.0*eV //GetMaxKinEnergy()
92  || eKin < 0 //GetMinKinEnergy()
93  || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
94 
95  if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) ) != dic.end()
96  || dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() ) return true;
97 
98  return false;
99 
100 // return IsApplicable( dp , element );
101 /*
102  G4double eKin = dp->GetKineticEnergy();
103  if ( eKin > 4.0*eV //GetMaxKinEnergy()
104  || eKin < 0 //GetMinKinEnergy()
105  || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
106  return true;
107 */
108 }
109 
111  G4int /*Z*/ , G4int /*A*/ ,
112  const G4Isotope* /*iso*/ ,
113  const G4Element* element ,
114  const G4Material* material )
115 {
116  //if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
117 
118  ke_cache = dp->GetKineticEnergy();
119  element_cache = element;
121  //G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
122  G4double xs = GetCrossSection( dp , element , material );
123  xs_cache = xs;
124  return xs;
125  //return GetCrossSection( dp , element , material->GetTemperature() );
126 }
127 
129 {
130  std::map< G4int , std::map< G4double , G4ParticleHPVector* >* >::iterator it;
131  std::map< G4double , G4ParticleHPVector* >::iterator itt;
132 
133  if ( coherent != NULL ) {
134  for ( it = coherent->begin() ; it != coherent->end() ; it++ )
135  {
136  if ( it->second != NULL )
137  {
138  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
139  {
140  delete itt->second;
141  }
142  }
143  delete it->second;
144  }
145  coherent->clear();
146  }
147 
148  if ( incoherent != NULL ) {
149  for ( it = incoherent->begin() ; it != incoherent->end() ; it++ )
150  {
151  if ( it->second != NULL )
152  {
153  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
154  {
155  delete itt->second;
156  }
157  }
158  delete it->second;
159  }
160  incoherent->clear();
161  }
162 
163  if ( inelastic != NULL ) {
164  for ( it = inelastic->begin() ; it != inelastic->end() ; it++ )
165  {
166  if ( it->second != NULL )
167  {
168  for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
169  {
170  delete itt->second;
171  }
172  }
173  delete it->second;
174  }
175  inelastic->clear();
176  }
177 
178 }
179 
180 
181 
183 {
184  G4bool result = false;
185 
186  G4double eKin = aP->GetKineticEnergy();
187  // Check energy
188  if ( eKin < emax )
189  {
190  // Check Particle Species
191  if ( aP->GetDefinition() == G4Neutron::Neutron() )
192  {
193  // anEle is one of Thermal elements
194  G4int ie = (G4int) anEle->GetIndex();
195  std::vector < G4int >::iterator it;
196  for ( it = indexOfThermalElement.begin() ; it != indexOfThermalElement.end() ; it++ )
197  {
198  if ( ie == *it ) return true;
199  }
200  }
201  }
202 
203 /*
204  if ( names->IsThisThermalElement ( anEle->GetName() ) )
205  {
206  // Check energy and projectile species
207  G4double eKin = aP->GetKineticEnergy();
208  if ( eKin < emax && aP->GetDefinition() == G4Neutron::Neutron() ) result = true;
209  }
210 */
211  return result;
212 }
213 
214 
216 {
217 
218  if ( &aP != G4Neutron::Neutron() )
219  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
220 
221  //std::map < std::pair < G4Material* , const G4Element* > , G4int > dic;
222  //
223  dic.clear();
225 
226  std::map < G4String , G4int > co_dic;
227 
228  //Searching Nist Materials
229  static G4ThreadLocal G4MaterialTable* theMaterialTable = 0 ; if (!theMaterialTable) theMaterialTable= G4Material::GetMaterialTable();
230  size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
231  for ( size_t i = 0 ; i < numberOfMaterials ; i++ )
232  {
233  G4Material* material = (*theMaterialTable)[i];
234  size_t numberOfElements = material->GetNumberOfElements();
235  for ( size_t j = 0 ; j < numberOfElements ; j++ )
236  {
237  const G4Element* element = material->GetElement(j);
238  if ( names->IsThisThermalElement ( material->GetName() , element->GetName() ) )
239  {
240  G4int ts_ID_of_this_geometry;
241  G4String ts_ndl_name = names->GetTS_NDL_Name( material->GetName() , element->GetName() );
242  if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
243  {
244  ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
245  }
246  else
247  {
248  ts_ID_of_this_geometry = co_dic.size();
249  co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
250  }
251 
252  //G4cout << "Neutron HP Thermal Scattering Data : Registering a material-element pair of "
253  // << material->GetName() << " " << element->GetName()
254  // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
255 
256  dic.insert( std::pair < std::pair < G4Material* , const G4Element* > , G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
257  }
258  }
259  }
260 
261  //Searching TS Elements
262  static G4ThreadLocal G4ElementTable* theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
263  size_t numberOfElements = G4Element::GetNumberOfElements();
264  //size_t numberOfThermalElements = 0;
265  for ( size_t i = 0 ; i < numberOfElements ; i++ )
266  {
267  const G4Element* element = (*theElementTable)[i];
268  if ( names->IsThisThermalElement ( element->GetName() ) )
269  {
270  if ( names->IsThisThermalElement ( element->GetName() ) )
271  {
272  G4int ts_ID_of_this_geometry;
273  G4String ts_ndl_name = names->GetTS_NDL_Name( element->GetName() );
274  if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
275  {
276  ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
277  }
278  else
279  {
280  ts_ID_of_this_geometry = co_dic.size();
281  co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
282  }
283 
284  //G4cout << "Neutron HP Thermal Scattering: Registering an element of "
285  // << material->GetName() << " " << element->GetName()
286  // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
287 
288  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 ) );
289  }
290  }
291  }
292 
293  G4cout << G4endl;
294  G4cout << "Neutron HP Thermal Scattering Data: Following material-element pairs and/or elements are registered." << G4endl;
295  for ( std::map < std::pair < const G4Material* , const G4Element* > , G4int >::iterator it = dic.begin() ; it != dic.end() ; it++ )
296  {
297  if ( it->first.first != NULL )
298  {
299  G4cout << "Material " << it->first.first->GetName() << " - Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
300  }
301  else
302  {
303  G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
304  }
305  }
306  G4cout << G4endl;
307 
308 
309  //G4cout << "Neutron HP Thermal Scattering Data: Following NDL thermal scattering files are assigned to the internal thermal scattering id." << G4endl;
310  //for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
311  //{
312  // G4cout << "NDL file name " << it->first << ", internal thermal scattering id " << it->second << G4endl;
313  //}
314 
316 
320 
321  if ( G4Threading::IsMasterThread() ) {
322 
323  if ( coherent == NULL ) coherent = new std::map< G4int , std::map< G4double , G4ParticleHPVector* >* >;
324  if ( incoherent == NULL ) incoherent = new std::map< G4int , std::map< G4double , G4ParticleHPVector* >* >;
325  if ( inelastic == NULL ) inelastic = new std::map< G4int , std::map< G4double , G4ParticleHPVector* >* >;
326 
327 
328  // Read Cross Section Data files
329 
330  G4String dirName;
331  if ( !std::getenv( "G4NEUTRONHPDATA" ) )
332  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
333  G4String baseName = std::getenv( "G4NEUTRONHPDATA" );
334 
335  dirName = baseName + "/ThermalScattering";
336 
337  G4String ndl_filename;
338  G4String full_name;
339 
340  for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
341  {
342 
343  ndl_filename = it->first;
344  G4int ts_ID = it->second;
345 
346  // Coherent
347  full_name = dirName + "/Coherent/CrossSection/" + ndl_filename;
348  std::map< G4double , G4ParticleHPVector* >* coh_amapTemp_EnergyCross = readData( full_name );
349  coherent->insert ( std::pair < G4int , std::map< G4double , G4ParticleHPVector* >* > ( ts_ID , coh_amapTemp_EnergyCross ) );
350 
351  // Incoherent
352  full_name = dirName + "/Incoherent/CrossSection/" + ndl_filename;
353  std::map< G4double , G4ParticleHPVector* >* incoh_amapTemp_EnergyCross = readData( full_name );
354  incoherent->insert ( std::pair < G4int , std::map< G4double , G4ParticleHPVector* >* > ( ts_ID , incoh_amapTemp_EnergyCross ) );
355 
356  // Inelastic
357  full_name = dirName + "/Inelastic/CrossSection/" + ndl_filename;
358  std::map< G4double , G4ParticleHPVector* >* inela_amapTemp_EnergyCross = readData( full_name );
359  inelastic->insert ( std::pair < G4int , std::map< G4double , G4ParticleHPVector* >* > ( ts_ID , inela_amapTemp_EnergyCross ) );
360 
361  }
365  }
366 }
367 
368 
369 
370 std::map< G4double , G4ParticleHPVector* >* G4ParticleHPThermalScatteringData::readData ( G4String full_name )
371 {
372 
373  std::map< G4double , G4ParticleHPVector* >* aData = new std::map< G4double , G4ParticleHPVector* >;
374 
375  //std::ifstream theChannel( full_name.c_str() );
376  std::istringstream theChannel;
377  G4ParticleHPManager::GetInstance()->GetDataStream(full_name,theChannel);
378 
379  //G4cout << "G4ParticleHPThermalScatteringData " << name << G4endl;
380 
381  G4int dummy;
382  while ( theChannel >> dummy ) // MF // Loop checking, 11.05.2015, T. Koi
383  {
384  theChannel >> dummy; // MT
385  G4double temp;
386  theChannel >> temp;
387  G4ParticleHPVector* anEnergyCross = new G4ParticleHPVector;
388  G4int nData;
389  theChannel >> nData;
390  anEnergyCross->Init ( theChannel , nData , eV , barn );
391  aData->insert ( std::pair < G4double , G4ParticleHPVector* > ( temp , anEnergyCross ) );
392  }
393  //theChannel.close();
394 
395  return aData;
396 
397 }
398 
399 
400 
402 {
403  if( &aP != G4Neutron::Neutron() )
404  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
405 // G4cout << "G4ParticleHPThermalScatteringData::DumpPhysicsTable still to be implemented"<<G4endl;
406 }
407 
408 //#include "G4Nucleus.hh"
409 //#include "G4NucleiPropertiesTable.hh"
410 //#include "G4Neutron.hh"
411 //#include "G4Electron.hh"
412 
413 
414 
415 /*
416 G4double G4ParticleHPThermalScatteringData::GetCrossSection( const G4DynamicParticle* aP , const G4Element*anE , G4double aT )
417 {
418 
419  G4double result = 0;
420  const G4Material* aM = NULL;
421 
422  G4int iele = anE->GetIndex();
423 
424  if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , anE ) ) != dic.end() )
425  {
426  iele = dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , anE ) )->second;
427  }
428  else if ( dic.find( std::pair < const G4Material* , const G4Element* > ( aM , anE ) ) != dic.end() )
429  {
430  iele = dic.find( std::pair < const G4Material* , const G4Element* > ( aM , anE ) )->second;
431  }
432  else
433  {
434  return result;
435  }
436 
437  G4double Xcoh = GetX ( aP , aT , coherent.find(iele)->second );
438  G4double Xincoh = GetX ( aP , aT , incoherent.find(iele)->second );
439  G4double Xinela = GetX ( aP , aT , inelastic.find(iele)->second );
440 
441  result = Xcoh + Xincoh + Xinela;
442 
443  //G4cout << "G4ParticleHPThermalScatteringData::GetCrossSection Tot= " << result/barn << " Coherent= " << Xcoh/barn << " Incoherent= " << Xincoh/barn << " Inelastic= " << Xinela/barn << G4endl;
444 
445  return result;
446 
447 }
448 */
449 
451 {
452  G4double result = 0;
453 
454  G4int ts_id =getTS_ID( aM , anE );
455 
456  if ( ts_id == -1 ) return result;
457 
458  G4double aT = aM->GetTemperature();
459 
460  G4double Xcoh = GetX ( aP , aT , coherent->find(ts_id)->second );
461  G4double Xincoh = GetX ( aP , aT , incoherent->find(ts_id)->second );
462  G4double Xinela = GetX ( aP , aT , inelastic->find(ts_id)->second );
463 
464  result = Xcoh + Xincoh + Xinela;
465 
466  //G4cout << "G4ParticleHPThermalScatteringData::GetCrossSection Tot= " << result/barn << " Coherent= " << Xcoh/barn << " Incoherent= " << Xincoh/barn << " Inelastic= " << Xinela/barn << G4endl;
467 
468  return result;
469 }
470 
471 
473 {
474  G4double result = 0;
475  G4int ts_id = getTS_ID( aM , anE );
476  G4double aT = aM->GetTemperature();
477  result = GetX ( aP , aT , inelastic->find( ts_id )->second );
478  return result;
479 }
480 
482 {
483  G4double result = 0;
484  G4int ts_id = getTS_ID( aM , anE );
485  G4double aT = aM->GetTemperature();
486  result = GetX ( aP , aT , coherent->find( ts_id )->second );
487  return result;
488 }
489 
491 {
492  G4double result = 0;
493  G4int ts_id = getTS_ID( aM , anE );
494  G4double aT = aM->GetTemperature();
495  result = GetX ( aP , aT , incoherent->find( ts_id )->second );
496  return result;
497 }
498 
499 
500 
502 {
503  G4int result = -1;
504  if ( dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) ) != dic.end() )
505  return dic.find( std::pair < const G4Material* , const G4Element* > ( (G4Material*)NULL , element ) )->second;
506  if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
507  return dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
508  return result;
509 }
510 
511 
512 
513 
514 G4double G4ParticleHPThermalScatteringData::GetX ( const G4DynamicParticle* aP, G4double aT , std::map < G4double , G4ParticleHPVector* >* amapTemp_EnergyCross )
515 {
516 
517  G4double result = 0;
518  if ( amapTemp_EnergyCross->size() == 0 ) return result;
519 
520 
521  G4double eKinetic = aP->GetKineticEnergy();
522 
523  if ( amapTemp_EnergyCross->size() == 1 ) {
524  if ( std::fabs ( aT - amapTemp_EnergyCross->begin()->first ) / amapTemp_EnergyCross->begin()->first > 0.1 ) {
525  G4cout << "G4ParticleHPThermalScatteringData:: The temperature of material ("
526  << aT/kelvin << "K) is different more than 10% from temperature of thermal scattering file expected ("
527  << amapTemp_EnergyCross->begin()->first << "K). Result may not be reliable."
528  << G4endl;
529  }
530  result = amapTemp_EnergyCross->begin()->second->GetXsec ( eKinetic );
531  return result;
532  }
533 
534  std::map< G4double , G4ParticleHPVector* >::iterator it;
535  for ( it = amapTemp_EnergyCross->begin() ; it != amapTemp_EnergyCross->end() ; it++ ) {
536  if ( aT < it->first ) break;
537  }
538  //if ( it == amapTemp_EnergyCross->begin() && it != amapTemp_EnergyCross->end() ) it++; // lower than the first
539  //if ( it != amapTemp_EnergyCross->begin() && it == amapTemp_EnergyCross->end() ) it--; // upper than the last
540  if ( it == amapTemp_EnergyCross->begin() ) {
541  it++; // lower than the first
542  } else if ( it == amapTemp_EnergyCross->end() ) {
543  it--; // upper than the last
544  }
545 
546  G4double TH = it->first;
547  G4double XH = it->second->GetXsec ( eKinetic );
548 
549  //G4cout << "G4ParticleHPThermalScatteringData::GetX TH " << TH << " E " << eKinetic << " XH " << XH << G4endl;
550 
551  if ( it != amapTemp_EnergyCross->begin() ) it--;
552  G4double TL = it->first;
553  G4double XL = it->second->GetXsec ( eKinetic );
554 
555  //G4cout << "G4ParticleHPThermalScatteringData::GetX TL " << TL << " E " << eKinetic << " XL " << XL << G4endl;
556 
557  if ( TH == TL )
558  throw G4HadronicException(__FILE__, __LINE__, "Thermal Scattering Data Error!");
559 
560  G4double T = aT;
561  G4double X = ( XH - XL ) / ( TH - TL ) * ( T - TL ) + XL;
562  result = X;
563 
564  return result;
565 }
566 
567 
569 {
570  names->AddThermalElement( nameG4Element , filename );
571 }
573 {
574  outFile << "High Precision cross data based on thermal scattering data in evaluated nuclear data libraries for neutrons below 5eV on specific materials\n" ;
575 }