ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPInelasticData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPInelasticData.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 // particle_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 "G4Neutron.hh"
42 #include "G4ElementTable.hh"
43 #include "G4ParticleHPData.hh"
44 #include "G4Pow.hh"
45 
48 {
49  const char* dataDirVariable;
50  G4String particleName;
51  if( projectile == G4Neutron::Neutron() ) {
52  dataDirVariable = "G4NEUTRONHPDATA";
53  }else if( projectile == G4Proton::Proton() ) {
54  dataDirVariable = "G4PROTONHPDATA";
55  particleName = "Proton";
56  }else if( projectile == G4Deuteron::Deuteron() ) {
57  dataDirVariable = "G4DEUTERONHPDATA";
58  particleName = "Deuteron";
59  }else if( projectile == G4Triton::Triton() ) {
60  dataDirVariable = "G4TRITONHPDATA";
61  particleName = "Triton";
62  }else if( projectile == G4He3::He3() ) {
63  dataDirVariable = "G4HE3HPDATA";
64  particleName = "He3";
65  }else if( projectile == G4Alpha::Alpha() ) {
66  dataDirVariable = "G4ALPHAHPDATA";
67  particleName = "Alpha";
68  } else {
69  G4String message("G4ParticleHPInelasticData may only be called for neutron, proton, deuteron, triton, He3 or alpha, while it is called for " + projectile->GetParticleName());
70  throw G4HadronicException(__FILE__, __LINE__,message.c_str());
71  }
72  // G4cout << this << " G4ParticleHPInelasticData::G4ParticleHPInelasticData " << projectile->GetParticleName() << " DATADIR " << dataDirVariable << G4endl;//GDEB
73  G4String dataName = projectile->GetParticleName()+"HPInelasticXS";
74  dataName.at(0) = toupper(dataName.at(0)) ;
75  SetName( dataName );
76 
77  if ( !std::getenv(dataDirVariable) && !std::getenv( "G4PARTICLEHPDATA" ) ){
78  G4String message("Please setenv G4PARTICLEHPDATA (recommended) or, at least setenv " +
79  G4String(dataDirVariable) + " to point to the " + projectile->GetParticleName() + " cross-section files.");
80  throw G4HadronicException(__FILE__, __LINE__,message.c_str());
81  }
82 
83  G4String dirName;
84  if ( std::getenv(dataDirVariable) ) {
85  dirName = std::getenv(dataDirVariable);
86  } else {
87  G4String baseName = std::getenv( "G4PARTICLEHPDATA" );
88  dirName = baseName + "/" + particleName;
89  }
90  G4cout << "@@@ G4ParticleHPInelasticData instantiated for particle " << projectile->GetParticleName() << " data directory variable is " << dataDirVariable << " pointing to " << dirName << G4endl;
91 
94 
95  onFlightDB = true;
96  theCrossSections = 0;
97  theProjectile=projectile;
98 
99  theHPData = NULL;
100  instanceOfWorker = false;
101  if ( G4Threading::IsMasterThread() ) {
103  } else {
104  instanceOfWorker = true;
105  }
106  element_cache = NULL;
107  material_cache = NULL;
108  ke_cache = 0.0;
109  xs_cache = 0.0;
110 }
111 
113 {
114  if ( theCrossSections != NULL && instanceOfWorker != true ) {
116  delete theCrossSections;
117  theCrossSections = NULL;
118  }
119  if ( theHPData != NULL && instanceOfWorker != true ) {
120  delete theHPData;
121  theHPData = NULL;
122  }
123 }
124 
126  G4int /*Z*/ , G4int /*A*/ ,
127  const G4Element* /*elm*/ ,
128  const G4Material* /*mat*/ )
129 {
130  G4double eKin = dp->GetKineticEnergy();
131  if ( eKin > GetMaxKinEnergy()
132  || eKin < GetMinKinEnergy()
133  || dp->GetDefinition() != theProjectile ) return false;
134 
135  return true;
136 }
137 
139  G4int /*Z*/ , G4int /*A*/ ,
140  const G4Isotope* /*iso*/ ,
141  const G4Element* element ,
142  const G4Material* material )
143 {
144  if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
145 
146  ke_cache = dp->GetKineticEnergy();
147  element_cache = element;
149  G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
150  xs_cache = xs;
151  return xs;
152 }
153 
154 /*
155 G4bool G4ParticleHPInelasticData::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
156 {
157  G4bool result = true;
158  G4double eKin = aP->GetKineticEnergy();
159  if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
160  return result;
161 }
162 */
163 
164 //void G4ParticleHPInelasticData::BuildPhysicsTableHP(G4ParticleDefinition* projectile,const char* /* dataDirVariable */)
166 {
167  // if(&projectile!=G4Neutron::Neutron())
168  // throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
169 
170 //080428
171  if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() )
172  {
173  G4cout << "Find a flag of \"G4PHP_NEGLECT_DOPPLER\"." << G4endl;
174  G4cout << "On the fly Doppler broadening will be neglect in the cross section calculation of inelastic scattering of neutrons (<20MeV)." << G4endl;
175  onFlightDB = false;
176  }
177 
178  if ( G4Threading::IsWorkerThread() ) {
180  return;
181  } else {
182  if ( theHPData == NULL ) theHPData = G4ParticleHPData::Instance( const_cast<G4ParticleDefinition*> ( &projectile ) );
183  }
184 
185 
186 
187  size_t numberOfElements = G4Element::GetNumberOfElements();
188 // theCrossSections = new G4PhysicsTable( numberOfElements );
189 // TKDB
190  //if ( theCrossSections == 0 )
191  //{ theCrossSections = new G4PhysicsTable( numberOfElements ); }
192  if ( theCrossSections == NULL )
193  theCrossSections = new G4PhysicsTable( numberOfElements );
194  else
196 
197  // make a PhysicsVector for each element
198 
199  //G4ParticleHPData* hpData = new G4ParticleHPData(projectile); //NEW
200  static G4ThreadLocal G4ElementTable *theElementTable = 0 ;
201  if (!theElementTable) theElementTable= G4Element::GetElementTable();
202  for( size_t i=0; i<numberOfElements; ++i )
203  {
204  //NEW G4PhysicsVector* physVec = G4ParticleHPData::
205  //NEW Instance(projectile, dataDirVariable)->MakePhysicsVector((*theElementTable)[i], this);
206  //G4PhysicsVector* physVec = hpData->MakePhysicsVector((*theElementTable)[i], this);
207  G4PhysicsVector* physVec = theHPData->MakePhysicsVector((*theElementTable)[i], this);
208  theCrossSections->push_back(physVec);
209  }
210 
212 }
213 
215 {
216  if(&projectile!=theProjectile)
217  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use ParticleHP data for a wrong projectile!!!");
218 
219 //
220 // Dump element based cross section
221 // range 10e-5 eV to 20 MeV
222 // 10 point per decade
223 // in barn
224 //
225 
226  G4cout << G4endl;
227  G4cout << G4endl;
228  G4cout << "Inelastic Cross Section of Neutron HP"<< G4endl;
229  G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
230  G4cout << G4endl;
231  G4cout << "Name of Element" << G4endl;
232  G4cout << "Energy[eV] XS[barn]" << G4endl;
233  G4cout << G4endl;
234 
235  size_t numberOfElements = G4Element::GetNumberOfElements();
236  static G4ThreadLocal G4ElementTable *theElementTable = 0 ;
237  if (!theElementTable) theElementTable= G4Element::GetElementTable();
238 
239  for ( size_t i = 0 ; i < numberOfElements ; ++i )
240  {
241 
242  G4cout << (*theElementTable)[i]->GetName() << G4endl;
243 
244  G4int ie = 0;
245 
246  for ( ie = 0 ; ie < 130 ; ie++ )
247  {
248  G4double eKinetic = 1.0e-5 * G4Pow::GetInstance()->powA ( 10.0 , ie/10.0 ) *CLHEP::eV;
249  G4bool outOfRange = false;
250 
251  if ( eKinetic < 20*CLHEP::MeV )
252  {
253  G4cout << eKinetic/CLHEP::eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/CLHEP::barn << G4endl;
254  }
255 
256  }
257 
258  G4cout << G4endl;
259  }
260 
261  //G4cout << "G4ParticleHPInelasticData::DumpPhysicsTable still to be implemented"<<G4endl;
262 }
263 
264 #include "G4NucleiProperties.hh"
265 
267 GetCrossSection(const G4DynamicParticle* projectile, const G4Element*anE, G4double aT)
268 {
269  G4double result = 0;
270  G4bool outOfRange;
271  G4int index = anE->GetIndex();
272 
273  // prepare neutron
274  G4double eKinetic = projectile->GetKineticEnergy();
275 
276  if ( !onFlightDB )
277  {
278  //NEGLECT_DOPPLER
279  G4double factor = 1.0;
280  if ( eKinetic < aT * CLHEP::k_Boltzmann )
281  {
282  // below 0.1 eV neutrons
283  // Have to do some, but now just igonre.
284  // Will take care after performance check.
285  // factor = factor * targetV;
286  }
287  return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
288 
289  }
290 
291  G4ReactionProduct theNeutron( projectile->GetDefinition() );
292  theNeutron.SetMomentum( projectile->GetMomentum() );
293  theNeutron.SetKineticEnergy( eKinetic );
294 
295  // prepare thermal nucleus
296  G4Nucleus aNuc;
297  G4double eps = 0.0001;
298  G4double theA = anE->GetN();
299  G4double theZ = anE->GetZ();
300  G4double eleMass;
301  eleMass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theA+eps), static_cast<G4int>(theZ+eps) );
302 
303  G4ReactionProduct boosted;
304  G4double aXsection;
305 
306  // MC integration loop
307  G4int counter = 0;
308  G4int failCount = 0;
309  G4double buffer = 0; G4int size = G4int(std::max(10., aT/60*CLHEP::kelvin));
310  G4ThreeVector neutronVelocity = 1./theProjectile->GetPDGMass()*theNeutron.GetMomentum();
311  G4double neutronVMag = neutronVelocity.mag();
312 
313  // G4cout << " G4ParticleHPInelasticData 2 " << size << G4endl;//GDEB
314 #ifndef G4PHP_DOPPLER_LOOP_ONCE
315  while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.01*buffer) // Loop checking, 11.05.2015, T. Koi
316  {
317  if(counter) buffer = result/counter;
318  while (counter<size) // Loop checking, 11.05.2015, T. Koi
319  {
320  counter ++;
321 #endif
322  //G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus( eleMass/theProjectile->GetPDGMass(), aT );
323  //G4Nucleus::GetThermalNucleus requests normalization of mass in neutron mass
324  G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus( eleMass/G4Neutron::Neutron()->GetPDGMass(), aT );
325  boosted.Lorentz(theNeutron, aThermalNuc);
326  G4double theEkin = boosted.GetKineticEnergy();
327  aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
328  // G4cout << " G4ParticleHPInelasticData aXsection " << aXsection << " index " << index << " theEkin " << theEkin << " outOfRange " << outOfRange <<G4endl;//GDEB
329  if(aXsection <0)
330  {
331  if(failCount<1000)
332  {
333  failCount++;
334 #ifndef G4PHP_DOPPLER_LOOP_ONCE
335  counter--;
336  continue;
337 #endif
338  }
339  else
340  {
341  aXsection = 0;
342  }
343  }
344  // velocity correction.
345  G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
346  aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
347  result += aXsection;
348  }
349 #ifndef G4PHP_DOPPLER_LOOP_ONCE
350  size += size;
351  }
352  result /= counter;
353 #endif
354 
355 /*
356  // Checking impact of G4PHP_NEGLECT_DOPPLER
357  G4cout << " result " << result << " "
358  << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) << " "
359  << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) /result << G4endl;
360 */
361 // G4cout << this << " G4ParticleHPInelasticData result " << result << G4endl; //GDEB
362 
363  return result;
364 }
365 
367 {
369 }
371 {
373 }
374 void G4ParticleHPInelasticData::CrossSectionDescription(std::ostream& outFile) const
375 {
376  outFile << "Extension of High Precision cross section for inelastic reaction of proton, deuteron, triton, He3 and alpha below 20MeV\n";
377 }