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