ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPCaptureData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPCaptureData.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 // 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 ) {
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;
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
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(std::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 }