Geant4_10
G4NeutronHPFissionData.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 // 070618 fix memory leaking by T. Koi
31 // 071002 enable cross section dump by T. Koi
32 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
33 // 081124 Protect invalid read which caused run time errors by T. Koi
34 // 100729 Add safty for 0 lenght cross sections by T. Koi
35 
37 #include "G4NeutronHPManager.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4Neutron.hh"
40 #include "G4ElementTable.hh"
41 #include "G4NeutronHPData.hh"
42 #include "G4NeutronHPManager.hh"
43 
45 :G4VCrossSectionDataSet("NeutronHPFissionXS")
46 {
47  SetMinKinEnergy( 0*MeV );
48  SetMaxKinEnergy( 20*MeV );
49 
50  ke_cache = 0.0;
51  xs_cache = 0.0;
52  element_cache = NULL;
53  material_cache = NULL;
54 
55  theCrossSections = 0;
56  //BuildPhysicsTable(*G4Neutron::Neutron());
57 }
58 
60 {
61  if ( theCrossSections != NULL ) theCrossSections->clearAndDestroy();
62  delete theCrossSections;
63 }
64 
66  G4int /*Z*/ , G4int /*A*/ ,
67  const G4Element* /*elm*/ ,
68  const G4Material* /*mat*/ )
69 {
70  G4double eKin = dp->GetKineticEnergy();
71  if ( eKin > GetMaxKinEnergy()
72  || eKin < GetMinKinEnergy()
73  || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
74 
75  return true;
76 }
77 
79  G4int /*Z*/ , G4int /*A*/ ,
80  const G4Isotope* /*iso*/ ,
81  const G4Element* element ,
82  const G4Material* material )
83 {
84  if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
85 
86  ke_cache = dp->GetKineticEnergy();
87  element_cache = element;
88  material_cache = material;
89  G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
90  xs_cache = xs;
91  return xs;
92 }
93 
94 /*
95 G4bool G4NeutronHPFissionData::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
96 {
97  G4bool result = true;
98  G4double eKin = aP->GetKineticEnergy();
99  if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
100  return result;
101 }
102 */
103 
105 {
106  if(&aP!=G4Neutron::Neutron())
107  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
108  size_t numberOfElements = G4Element::GetNumberOfElements();
109  //theCrossSections = new G4PhysicsTable( numberOfElements );
110  // TKDB
111  //if ( theCrossSections == NULL ) theCrossSections = new G4PhysicsTable( numberOfElements );
112  if ( theCrossSections == NULL )
113  theCrossSections = new G4PhysicsTable( numberOfElements );
114  else
115  theCrossSections->clearAndDestroy();
116 
117  // make a PhysicsVector for each element
118 
119  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
120  for( size_t i=0; i<numberOfElements; ++i )
121  {
123  Instance()->MakePhysicsVector((*theElementTable)[i], this);
124  theCrossSections->push_back(physVec);
125  }
126 }
127 
129 {
130  if(&aP!=G4Neutron::Neutron())
131  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
132 
133 //
134 // Dump element based cross section
135 // range 10e-5 eV to 20 MeV
136 // 10 point per decade
137 // in barn
138 //
139 
140  G4cout << G4endl;
141  G4cout << G4endl;
142  G4cout << "Fission Cross Section of Neutron HP"<< G4endl;
143  G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
144  G4cout << G4endl;
145  G4cout << "Name of Element" << G4endl;
146  G4cout << "Energy[eV] XS[barn]" << G4endl;
147  G4cout << G4endl;
148 
149  size_t numberOfElements = G4Element::GetNumberOfElements();
150  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
151 
152  for ( size_t i = 0 ; i < numberOfElements ; ++i )
153  {
154 
155  G4cout << (*theElementTable)[i]->GetName() << G4endl;
156 
157  if ( (*((*theCrossSections)(i))).GetVectorLength() == 0 )
158  {
159  G4cout << "The cross-section data of the fission of this element is not available." << G4endl;
160  G4cout << G4endl;
161  continue;
162  }
163 
164  G4int ie = 0;
165 
166  for ( ie = 0 ; ie < 130 ; ie++ )
167  {
168  G4double eKinetic = 1.0e-5 * std::pow ( 10.0 , ie/10.0 ) *eV;
169  G4bool outOfRange = false;
170 
171  if ( eKinetic < 20*MeV )
172  {
173  G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
174  }
175 
176  }
177 
178  G4cout << G4endl;
179  }
180 
181  //G4cout << "G4NeutronHPFissionData::DumpPhysicsTable still to be implemented"<<G4endl;
182 }
183 
184 #include "G4NucleiProperties.hh"
185 
188 {
189  G4double result = 0;
190  if(anE->GetZ()<90) return result;
191  G4bool outOfRange;
192  G4int index = anE->GetIndex();
193 
194 // 100729 TK add safety
195 if ( ( ( *theCrossSections )( index ) )->GetVectorLength() == 0 ) return result;
196 
197  // prepare neutron
198  G4double eKinetic = aP->GetKineticEnergy();
199  G4ReactionProduct theNeutron( aP->GetDefinition() );
200  theNeutron.SetMomentum( aP->GetMomentum() );
201  theNeutron.SetKineticEnergy( eKinetic );
202 
203  // prepare thermal nucleus
204  G4Nucleus aNuc;
205  G4double eps = 0.0001;
206  G4double theA = anE->GetN();
207  G4double theZ = anE->GetZ();
208  G4double eleMass;
209  eleMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theA+eps) , static_cast<G4int>(theZ+eps) )
211 
212  G4ReactionProduct boosted;
213  G4double aXsection;
214 
215  // MC integration loop
216  G4int counter = 0;
217  G4double buffer = 0;
218  G4int size = G4int(std::max(10., aT/60*kelvin));
219  G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
220  G4double neutronVMag = neutronVelocity.mag();
221 
222  while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.01*buffer)
223  {
224  if(counter) buffer = result/counter;
225  while (counter<size)
226  {
227  counter ++;
228  G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
229  boosted.Lorentz(theNeutron, aThermalNuc);
230  G4double theEkin = boosted.GetKineticEnergy();
231  aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
232  // velocity correction.
233  G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
234  aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
235  result += aXsection;
236  }
237  size += size;
238  }
239  result /= counter;
240  return result;
241 }
242 
244 {
246 }
248 {
250 }
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
G4bool IsIsoApplicable(const G4DynamicParticle *, G4int, G4int, const G4Element *, const G4Material *)
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
string material
Definition: eplot.py:19
G4GLOB_DLL std::ostream G4cout
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, G4double aT)
static size_t GetNumberOfElements()
Definition: G4Element.cc:402
void SetMinKinEnergy(G4double value)
bool G4bool
Definition: G4Types.hh:79
size_t GetIndex() const
Definition: G4Element.hh:181
G4double GetIsoCrossSection(const G4DynamicParticle *, G4int, G4int, const G4Isotope *, const G4Element *, const G4Material *)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4NeutronHPData * Instance()
G4double GetKineticEnergy() const
void BuildPhysicsTable(const G4ParticleDefinition &)
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
void DumpPhysicsTable(const G4ParticleDefinition &)
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