Geant4_10
G4NeutronHPInelasticData.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 //
38 #include "G4NeutronHPManager.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4Neutron.hh"
42 #include "G4ElementTable.hh"
43 #include "G4NeutronHPData.hh"
44 #include "G4NeutronHPManager.hh"
45 
47 :G4VCrossSectionDataSet("NeutronHPInelasticXS")
48 {
49 
50  SetMinKinEnergy( 0*MeV );
51  SetMaxKinEnergy( 20*MeV );
52 
53  ke_cache = 0.0;
54  xs_cache = 0.0;
55  element_cache = NULL;
56  material_cache = NULL;
57 
58  onFlightDB = true;
59  theCrossSections = 0;
60  //BuildPhysicsTable(*G4Neutron::Neutron());
61 }
62 
64 {
65  if ( theCrossSections != 0 ) theCrossSections->clearAndDestroy();
66  delete theCrossSections;
67 }
68 
70  G4int /*Z*/ , G4int /*A*/ ,
71  const G4Element* /*elm*/ ,
72  const G4Material* /*mat*/ )
73 {
74  G4double eKin = dp->GetKineticEnergy();
75  if ( eKin > GetMaxKinEnergy()
76  || eKin < GetMinKinEnergy()
77  || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
78 
79  return true;
80 }
81 
83  G4int /*Z*/ , G4int /*A*/ ,
84  const G4Isotope* /*iso*/ ,
85  const G4Element* element ,
86  const G4Material* material )
87 {
88  if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
89 
90  ke_cache = dp->GetKineticEnergy();
91  element_cache = element;
92  material_cache = material;
93  G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
94  xs_cache = xs;
95  return xs;
96  //return GetCrossSection( dp , element , material->GetTemperature() );
97 }
98 
99 /*
100 G4bool G4NeutronHPInelasticData::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
101 {
102  G4bool result = true;
103  G4double eKin = aP->GetKineticEnergy();
104  if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
105  return result;
106 }
107 */
108 
110 {
111  if(&aP!=G4Neutron::Neutron())
112  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
113 
114 //080428
115  if ( getenv( "G4NEUTRONHP_NEGLECT_DOPPLER" ) )
116  {
117  G4cout << "Find environment variable of \"G4NEUTRONHP_NEGLECT_DOPPLER\"." << G4endl;
118  G4cout << "On the fly Doppler broadening will be neglect in the cross section calculation of inelastic scattering of neutrons (<20MeV)." << G4endl;
119  onFlightDB = false;
120  }
121 
122  size_t numberOfElements = G4Element::GetNumberOfElements();
123 // theCrossSections = new G4PhysicsTable( numberOfElements );
124 // TKDB
125  //if ( theCrossSections == 0 )
126  //{ theCrossSections = new G4PhysicsTable( numberOfElements ); }
127  if ( theCrossSections == NULL )
128  theCrossSections = new G4PhysicsTable( numberOfElements );
129  else
130  theCrossSections->clearAndDestroy();
131 
132  // make a PhysicsVector for each element
133 
134  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
135  for( size_t i=0; i<numberOfElements; ++i )
136  {
138  Instance()->MakePhysicsVector((*theElementTable)[i], this);
139  theCrossSections->push_back(physVec);
140  }
141 }
142 
144 {
145  if(&aP!=G4Neutron::Neutron())
146  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
147 
148 //
149 // Dump element based cross section
150 // range 10e-5 eV to 20 MeV
151 // 10 point per decade
152 // in barn
153 //
154 
155  G4cout << G4endl;
156  G4cout << G4endl;
157  G4cout << "Inelastic Cross Section of Neutron HP"<< G4endl;
158  G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
159  G4cout << G4endl;
160  G4cout << "Name of Element" << G4endl;
161  G4cout << "Energy[eV] XS[barn]" << G4endl;
162  G4cout << G4endl;
163 
164  size_t numberOfElements = G4Element::GetNumberOfElements();
165  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
166 
167  for ( size_t i = 0 ; i < numberOfElements ; ++i )
168  {
169 
170  G4cout << (*theElementTable)[i]->GetName() << G4endl;
171 
172  G4int ie = 0;
173 
174  for ( ie = 0 ; ie < 130 ; ie++ )
175  {
176  G4double eKinetic = 1.0e-5 * std::pow ( 10.0 , ie/10.0 ) *eV;
177  G4bool outOfRange = false;
178 
179  if ( eKinetic < 20*MeV )
180  {
181  G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
182  }
183 
184  }
185 
186  G4cout << G4endl;
187  }
188 
189  //G4cout << "G4NeutronHPInelasticData::DumpPhysicsTable still to be implemented"<<G4endl;
190 }
191 
192 #include "G4NucleiProperties.hh"
193 
196 {
197  G4double result = 0;
198  G4bool outOfRange;
199  G4int index = anE->GetIndex();
200 
201  // prepare neutron
202  G4double eKinetic = aP->GetKineticEnergy();
203 
204  // T. K.
205 //if ( getenv( "G4NEUTRONHP_NEGLECT_DOPPLER" ) )
206 //080428
207  if ( !onFlightDB )
208  {
209  G4double factor = 1.0;
210  if ( eKinetic < aT * k_Boltzmann )
211  {
212  // below 0.1 eV neutrons
213  // Have to do some, but now just igonre.
214  // Will take care after performance check.
215  // factor = factor * targetV;
216  }
217  return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
218  }
219 
220  G4ReactionProduct theNeutron( aP->GetDefinition() );
221  theNeutron.SetMomentum( aP->GetMomentum() );
222  theNeutron.SetKineticEnergy( eKinetic );
223 
224  // prepare thermal nucleus
225  G4Nucleus aNuc;
226  G4double eps = 0.0001;
227  G4double theA = anE->GetN();
228  G4double theZ = anE->GetZ();
229  G4double eleMass;
230  eleMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theA+eps), static_cast<G4int>(theZ+eps))
232 
233  G4ReactionProduct boosted;
234  G4double aXsection;
235 
236  // MC integration loop
237  G4int counter = 0;
238  G4int failCount = 0;
239  G4double buffer = 0;
240  G4int size = G4int(std::max(10., aT/60*kelvin));
241  G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
242  G4double neutronVMag = neutronVelocity.mag();
243 
244  while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.01*buffer)
245  {
246  if(counter) buffer = result/counter;
247  while (counter<size)
248  {
249  counter ++;
250  G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
251  boosted.Lorentz(theNeutron, aThermalNuc);
252  G4double theEkin = boosted.GetKineticEnergy();
253  aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
254  if(aXsection <0)
255  {
256  if(failCount<1000)
257  {
258  failCount++;
259  counter--;
260  continue;
261  }
262  else
263  {
264  aXsection = 0;
265  }
266  }
267  // velocity correction.
268  G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
269  aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
270  result += aXsection;
271  }
272  size += size;
273  }
274  result /= counter;
275 /*
276  // Checking impact of G4NEUTRONHP_NEGLECT_DOPPLER
277  G4cout << " result " << result << " "
278  << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) << " "
279  << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) /result << G4endl;
280 */
281  return result;
282 }
283 
285 {
287 }
289 {
291 }
G4PhysicsVector * MakePhysicsVector(G4Element *thE, G4NeutronHPFissionData *theP)
static G4double GetNuclearMass(const G4double A, const G4double Z)
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4double GetKineticEnergy() const
Int_t index
Definition: macro.C:9
G4double GetN() const
Definition: G4Element.hh:134
static G4NeutronHPManager * GetInstance()
void SetMomentum(const G4double x, const G4double y, const G4double z)
void push_back(G4PhysicsVector *)
G4double GetZ() const
Definition: G4Element.hh:131
G4double G4NeutronHPJENDLHEData::G4double result
#define buffer
Definition: xmlparse.cc:611
G4ParticleDefinition * GetDefinition() const
#define G4ThreadLocal
Definition: tls.hh:52
G4ReactionProduct GetThermalNucleus(G4double aMass, G4double temp=-1) const
Definition: G4Nucleus.cc:130
int G4int
Definition: G4Types.hh:78
void BuildPhysicsTable(const G4ParticleDefinition &)
G4double GetIsoCrossSection(const G4DynamicParticle *, G4int, G4int, const G4Isotope *, const G4Element *, const G4Material *)
string material
Definition: eplot.py:19
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, G4double aT)
G4GLOB_DLL std::ostream G4cout
float k_Boltzmann
Definition: hepunit.py:299
static size_t GetNumberOfElements()
Definition: G4Element.cc:402
void SetMinKinEnergy(G4double value)
bool G4bool
Definition: G4Types.hh:79
void DumpPhysicsTable(const G4ParticleDefinition &)
size_t GetIndex() const
Definition: G4Element.hh:181
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4NeutronHPData * Instance()
G4double GetKineticEnergy() const
void SetVerboseLevel(G4int i)
G4double GetPDGMass() const
void SetMaxKinEnergy(G4double value)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4bool IsIsoApplicable(const G4DynamicParticle *, G4int, G4int, const G4Element *, const G4Material *)
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:395
double mag() const
G4double GetMass() const
void clearAndDestroy()
G4ThreeVector GetMomentum() const