Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPElasticData.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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 070523 add neglecting doppler broadening on the fly. T. Koi
31 // 070613 fix memory leaking by T. Koi
32 // 071002 enable cross section dump by T. Koi
33 // 080428 change checking point of "neglecting doppler broadening" flag
34 // from GetCrossSection to BuildPhysicsTable by T. Koi
35 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
36 //
37 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
38 //
40 #include "G4ParticleHPManager.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4Neutron.hh"
44 #include "G4ElementTable.hh"
45 #include "G4ParticleHPData.hh"
46 #include "G4ParticleHPManager.hh"
47 #include "G4Pow.hh"
48 
50 :G4VCrossSectionDataSet("NeutronHPElasticXS")
51 {
52  SetMinKinEnergy( 0*MeV );
53  SetMaxKinEnergy( 20*MeV );
54 
55  theCrossSections = 0;
56  onFlightDB = true;
57  instanceOfWorker = false;
59  instanceOfWorker = true;
60  }
61  element_cache = NULL;
62  material_cache = NULL;
63  ke_cache = 0.0;
64  xs_cache = 0.0;
65 // BuildPhysicsTable( *G4Neutron::Neutron() );
66 }
67 
69 {
70  if ( theCrossSections != NULL && instanceOfWorker != true ) {
71  theCrossSections->clearAndDestroy();
72  delete theCrossSections;
73  theCrossSections = NULL;
74  }
75 }
76 
78  G4int /*Z*/ , G4int /*A*/ ,
79  const G4Element* /*elm*/ ,
80  const G4Material* /*mat*/ )
81 {
82 
83  G4double eKin = dp->GetKineticEnergy();
84  if ( eKin > GetMaxKinEnergy()
85  || eKin < GetMinKinEnergy()
86  || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
87 
88  return true;
89 }
90 
92  G4int /*Z*/ , G4int /*A*/ ,
93  const G4Isotope* /*iso*/ ,
94  const G4Element* element ,
95  const G4Material* material )
96 {
97  if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
98 
99  ke_cache = dp->GetKineticEnergy();
100  element_cache = element;
101  material_cache = material;
102  G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
103  xs_cache = xs;
104  return xs;
105 }
106 
107 /*
108 G4bool G4ParticleHPElasticData::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
109 {
110  G4bool result = true;
111  G4double eKin = aP->GetKineticEnergy();
112  if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
113  return result;
114 }
115 */
116 
118 {
119 
120  if(&aP!=G4Neutron::Neutron())
121  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
122 
123 //080428
124  if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() )
125  {
126  G4cout << "Find a flag of \"G4NEUTRONHP_NEGLECT_DOPPLER\"." << G4endl;
127  G4cout << "On the fly Doppler broadening will be neglect in the cross section calculation of elastic scattering of neutrons (<20MeV)." << G4endl;
128  onFlightDB = false;
129  }
130 
131  if ( G4Threading::IsWorkerThread() ) {
133  return;
134  }
135 
136  size_t numberOfElements = G4Element::GetNumberOfElements();
137 // TKDB
138  //if ( theCrossSections == 0 ) theCrossSections = new G4PhysicsTable( numberOfElements );
139  if ( theCrossSections == NULL )
140  theCrossSections = new G4PhysicsTable( numberOfElements );
141  else
142  theCrossSections->clearAndDestroy();
143 
144  // make a PhysicsVector for each element
145 
146  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
147  for( size_t i=0; i<numberOfElements; ++i )
148  {
150  Instance(G4Neutron::Neutron())->MakePhysicsVector((*theElementTable)[i], this);
151  theCrossSections->push_back(physVec);
152  }
153 
155 }
156 
158 {
159  if(&aP!=G4Neutron::Neutron())
160  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
161 
162 //
163 // Dump element based cross section
164 // range 10e-5 eV to 20 MeV
165 // 10 point per decade
166 // in barn
167 //
168 
169  G4cout << G4endl;
170  G4cout << G4endl;
171  G4cout << "Elastic Cross Section of Neutron HP"<< G4endl;
172  G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
173  G4cout << G4endl;
174  G4cout << "Name of Element" << G4endl;
175  G4cout << "Energy[eV] XS[barn]" << G4endl;
176  G4cout << G4endl;
177 
178  size_t numberOfElements = G4Element::GetNumberOfElements();
179  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
180 
181  for ( size_t i = 0 ; i < numberOfElements ; ++i )
182  {
183 
184  G4cout << (*theElementTable)[i]->GetName() << G4endl;
185 
186  G4int ie = 0;
187 
188  for ( ie = 0 ; ie < 130 ; ie++ )
189  {
190  G4double eKinetic = 1.0e-5 * G4Pow::GetInstance()->powA ( 10.0 , ie/10.0 ) *eV;
191  G4bool outOfRange = false;
192 
193  if ( eKinetic < 20*MeV )
194  {
195  G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
196  }
197 
198  }
199 
200  G4cout << G4endl;
201  }
202 
203 
204 // G4cout << "G4ParticleHPElasticData::DumpPhysicsTable still to be implemented"<<G4endl;
205 }
206 
207 #include "G4Nucleus.hh"
208 #include "G4NucleiProperties.hh"
209 #include "G4Neutron.hh"
210 #include "G4Electron.hh"
211 
214 {
215  G4double result = 0;
216  G4bool outOfRange;
217  G4int index = anE->GetIndex();
218 
219  // prepare neutron
220  G4double eKinetic = aP->GetKineticEnergy();
221 
222  if ( !onFlightDB )
223  {
224  //NEGLECT_DOPPLER_B.
225  G4double factor = 1.0;
226  if ( eKinetic < aT * k_Boltzmann )
227  {
228  // below 0.1 eV neutrons
229  // Have to do some, but now just igonre.
230  // Will take care after performance check.
231  // factor = factor * targetV;
232  }
233  return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
234  }
235 
236  G4ReactionProduct theNeutron( aP->GetDefinition() );
237  theNeutron.SetMomentum( aP->GetMomentum() );
238  theNeutron.SetKineticEnergy( eKinetic );
239 
240  // prepare thermal nucleus
241  G4Nucleus aNuc;
242  G4double eps = 0.0001;
243  G4double theA = anE->GetN();
244  G4double theZ = anE->GetZ();
245  G4double eleMass;
246 
247 
248  eleMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theA+eps) , static_cast<G4int>(theZ+eps) )
250 
251  G4ReactionProduct boosted;
252  G4double aXsection;
253 
254  // MC integration loop
255  G4int counter = 0;
256  G4double buffer = 0;
257  G4int size = G4int(std::max(10., aT/60*kelvin));
258  G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
259  G4double neutronVMag = neutronVelocity.mag();
260 
261  while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.03*buffer) // Loop checking, 11.05.2015, T. Koi
262  {
263  if(counter) buffer = result/counter;
264  while (counter<size) // Loop checking, 11.05.2015, T. Koi
265  {
266  counter ++;
267  G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
268  boosted.Lorentz(theNeutron, aThermalNuc);
269  G4double theEkin = boosted.GetKineticEnergy();
270  aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
271  // velocity correction.
272  G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
273  aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
274  result += aXsection;
275  }
276  size += size;
277  }
278  result /= counter;
279 /*
280  // Checking impact of G4NEUTRONHP_NEGLECT_DOPPLER
281  G4cout << " result " << result << " "
282  << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) << " "
283  << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) /result << G4endl;
284 */
285  return result;
286 }
287 
290 {
292 }
293 
295 SetVerboseLevel( G4int newValue )
296 {
298 }
299 void G4ParticleHPElasticData::CrossSectionDescription(std::ostream& outFile) const
300 {
301  outFile << "High Precision cross data based on Evaluated Nuclear Data Files (ENDF) for elastic reaction of neutrons below 20MeV\n" ;
302 }
G4double G4ParticleHPJENDLHEData::G4double result
static G4ParticleHPManager * GetInstance()
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
static G4double GetNuclearMass(const G4double A, const G4double Z)
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4double GetKineticEnergy() const
G4double GetIsoCrossSection(const G4DynamicParticle *, G4int, G4int, const G4Isotope *, const G4Element *, const G4Material *)
void BuildPhysicsTable(const G4ParticleDefinition &)
G4double GetN() const
Definition: G4Element.hh:135
void SetMomentum(const G4double x, const G4double y, const G4double z)
void push_back(G4PhysicsVector *)
G4double GetZ() const
Definition: G4Element.hh:131
void DumpPhysicsTable(const G4ParticleDefinition &)
#define buffer
Definition: xmlparse.cc:628
G4bool IsIsoApplicable(const G4DynamicParticle *, G4int, G4int, const G4Element *, const G4Material *)
static const G4double eps
G4ParticleDefinition * GetDefinition() const
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, G4double aT)
#define G4ThreadLocal
Definition: tls.hh:89
G4ReactionProduct GetThermalNucleus(G4double aMass, G4double temp=-1) const
Definition: G4Nucleus.cc:143
int G4int
Definition: G4Types.hh:78
virtual void CrossSectionDescription(std::ostream &) const
string material
Definition: eplot.py:19
void RegisterElasticCrossSections(G4PhysicsTable *val)
G4GLOB_DLL std::ostream G4cout
float k_Boltzmann
Definition: hepunit.py:299
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
G4PhysicsTable * GetElasticCrossSections()
void SetMinKinEnergy(G4double value)
bool G4bool
Definition: G4Types.hh:79
size_t GetIndex() const
Definition: G4Element.hh:182
static constexpr double eV
Definition: G4SIunits.hh:215
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static constexpr double kelvin
Definition: G4SIunits.hh:281
G4bool IsWorkerThread()
Definition: G4Threading.cc:145
G4double GetKineticEnergy() const
G4PhysicsVector * MakePhysicsVector(G4Element *thE, G4ParticleHPFissionData *theP)
G4double GetPDGMass() const
void SetMaxKinEnergy(G4double value)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4ParticleHPData * Instance(G4ParticleDefinition *projectile)
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
static constexpr double barn
Definition: G4SIunits.hh:105
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
double mag() const
G4double GetMass() const
void clearAndDestroy()
G4ThreeVector GetMomentum() const