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