ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPFissionData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPFissionData.cc
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. Ko
35 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 //
38 #include "G4ParticleHPManager.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4Neutron.hh"
42 #include "G4ElementTable.hh"
43 #include "G4ParticleHPData.hh"
44 #include "G4ParticleHPManager.hh"
45 #include "G4Pow.hh"
46 
48 :G4VCrossSectionDataSet("NeutronHPFissionXS")
49 {
50  SetMinKinEnergy( 0*MeV );
51  SetMaxKinEnergy( 20*MeV );
52 
53  theCrossSections = 0;
54  onFlightDB = true;
55  instanceOfWorker = false;
57  instanceOfWorker = true;
58  }
59  element_cache = NULL;
60  material_cache = NULL;
61  ke_cache = 0.0;
62  xs_cache = 0.0;
63  //BuildPhysicsTable(*G4Neutron::Neutron());
64 }
65 
67 {
68  if ( theCrossSections != NULL && instanceOfWorker != true ) {
70  delete theCrossSections;
71  theCrossSections = NULL;
72  }
73 }
74 
76  G4int /*Z*/ , G4int /*A*/ ,
77  const G4Element* /*elm*/ ,
78  const G4Material* /*mat*/ )
79 {
80  G4double eKin = dp->GetKineticEnergy();
81  if ( eKin > GetMaxKinEnergy()
82  || eKin < GetMinKinEnergy()
83  || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
84 
85  return true;
86 }
87 
89  G4int /*Z*/ , G4int /*A*/ ,
90  const G4Isotope* /*iso*/ ,
91  const G4Element* element ,
92  const G4Material* material )
93 {
94  if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
95 
96  ke_cache = dp->GetKineticEnergy();
97  element_cache = element;
99  G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
100  xs_cache = xs;
101  return xs;
102 }
103 
104 /*
105 G4bool G4ParticleHPFissionData::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
106 {
107  G4bool result = true;
108  G4double eKin = aP->GetKineticEnergy();
109  if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
110  return result;
111 }
112 */
113 
115 {
116 
117  if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
118  G4cout << "Find a flag of \"G4NEUTRONHP_NEGLECT_DOPPLER\"." << G4endl;
119  G4cout << "On the fly Doppler broadening will be neglect in the cross section calculation of fission reaction of neutrons (<20MeV)." << G4endl;
120  onFlightDB = false;
121  }
122 
123  if(&aP!=G4Neutron::Neutron())
124  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
125 
126  if ( G4Threading::IsWorkerThread() ) {
128  return;
129  }
130 
131  size_t numberOfElements = G4Element::GetNumberOfElements();
132  //theCrossSections = new G4PhysicsTable( numberOfElements );
133  // TKDB
134  //if ( theCrossSections == NULL ) theCrossSections = new G4PhysicsTable( numberOfElements );
135  if ( theCrossSections == NULL )
136  theCrossSections = new G4PhysicsTable( numberOfElements );
137  else
139 
140  // make a PhysicsVector for each element
141 
142  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
143  for( size_t i=0; i<numberOfElements; ++i )
144  {
146  Instance(G4Neutron::Neutron())->MakePhysicsVector((*theElementTable)[i], this);
147  theCrossSections->push_back(physVec);
148  }
149 
151 }
152 
154 {
155  if(&aP!=G4Neutron::Neutron())
156  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
157 
158 //
159 // Dump element based cross section
160 // range 10e-5 eV to 20 MeV
161 // 10 point per decade
162 // in barn
163 //
164 
165  G4cout << G4endl;
166  G4cout << G4endl;
167  G4cout << "Fission Cross Section of Neutron HP"<< G4endl;
168  G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
169  G4cout << G4endl;
170  G4cout << "Name of Element" << G4endl;
171  G4cout << "Energy[eV] XS[barn]" << G4endl;
172  G4cout << G4endl;
173 
174  size_t numberOfElements = G4Element::GetNumberOfElements();
175  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
176 
177  for ( size_t i = 0 ; i < numberOfElements ; ++i )
178  {
179 
180  G4cout << (*theElementTable)[i]->GetName() << G4endl;
181 
182  if ( (*((*theCrossSections)(i))).GetVectorLength() == 0 )
183  {
184  G4cout << "The cross-section data of the fission of this element is not available." << G4endl;
185  G4cout << G4endl;
186  continue;
187  }
188 
189  G4int ie = 0;
190 
191  for ( ie = 0 ; ie < 130 ; ie++ )
192  {
193  G4double eKinetic = 1.0e-5 * G4Pow::GetInstance()->powA ( 10.0 , ie/10.0 ) *eV;
194  G4bool outOfRange = false;
195 
196  if ( eKinetic < 20*MeV )
197  {
198  G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
199  }
200 
201  }
202 
203  G4cout << G4endl;
204  }
205 
206  //G4cout << "G4ParticleHPFissionData::DumpPhysicsTable still to be implemented"<<G4endl;
207 }
208 
209 #include "G4NucleiProperties.hh"
210 
213 {
214  G4double result = 0;
215  if(anE->GetZ()<88) return result;
216  G4bool outOfRange;
217  G4int index = anE->GetIndex();
218 
219 // 100729 TK add safety
220 if ( ( ( *theCrossSections )( index ) )->GetVectorLength() == 0 ) return result;
221 
222  // prepare neutron
223  G4double eKinetic = aP->GetKineticEnergy();
224  G4ReactionProduct theNeutronRP( aP->GetDefinition() );
225  theNeutronRP.SetMomentum( aP->GetMomentum() );
226  theNeutronRP.SetKineticEnergy( eKinetic );
227 
228  if ( !onFlightDB ) {
229  //NEGLECT_DOPPLER
230  G4double factor = 1.0;
231  if ( eKinetic < aT * k_Boltzmann ) {
232  // below 0.1 eV neutrons
233  // Have to do some, but now just igonre.
234  // Will take care after performance check.
235  // factor = factor * targetV;
236  }
237  return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
238  }
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  eleMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theA+eps) , static_cast<G4int>(theZ+eps) )
248 
249  G4ReactionProduct boosted;
250  G4double aXsection;
251 
252  // MC integration loop
253  G4int counter = 0;
254  G4double buffer = 0;
255  G4int size = G4int(std::max(10., aT/60*kelvin));
256  G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutronRP.GetMomentum();
257  G4double neutronVMag = neutronVelocity.mag();
258 
259  while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.01*buffer) // Loop checking, 11.05.2015, T. Koi
260  {
261  if(counter) buffer = result/counter;
262  while (counter<size) // Loop checking, 11.05.2015, T. Koi
263  {
264  counter ++;
265  G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
266  boosted.Lorentz(theNeutronRP, aThermalNuc);
267  G4double theEkin = boosted.GetKineticEnergy();
268  aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
269  // velocity correction.
270  G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
271  aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
272  result += aXsection;
273  }
274  size += size;
275  }
276  result /= counter;
277  return result;
278 }
279 
281 {
283 }
285 {
287 }
288 void G4ParticleHPFissionData::CrossSectionDescription(std::ostream& outFile) const
289 {
290  outFile << "High Precision cross data based on Evaluated Nuclear Data Files (ENDF) for induced fission reaction of neutrons below 20MeV\n" ;
291 }