ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPJENDLHEData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPJENDLHEData.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 // Class Description
27 // Cross-section data set for a high precision (based on JENDL_HE evaluated data
28 // libraries) description of elastic scattering 20 MeV ~ 3 GeV;
29 // Class Description - End
30 
31 // 15-Nov-06 First Implementation is done by T. Koi (SLAC/SCCS)
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
35 #include "G4SystemOfUnits.hh"
36 #include "G4LPhysicsFreeVector.hh"
37 #include "G4ElementTable.hh"
38 #include "G4ParticleHPData.hh"
39 #include "G4Pow.hh"
40 
42 {
43 
44  G4bool result = true;
45  G4double eKin = aP->GetKineticEnergy();
46  //if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
47  if ( eKin < 20*MeV || 3*GeV < eKin || aP->GetDefinition()!=G4Neutron::Neutron() )
48  {
49  result = false;
50  }
51 // Element Check
52  else if ( !(vElement[ anE->GetIndex() ]) ) result = false;
53 
54  return result;
55 
56 }
57 
58 
59 
61 {
62  for ( std::map< G4int , std::map< G4int , G4PhysicsVector* >* >::iterator itZ = mIsotope.begin();
63  itZ != mIsotope.end(); ++itZ ) {
64  std::map< G4int , G4PhysicsVector* >* pointer_map = itZ->second;
65  if ( pointer_map ) {
66  for ( std::map< G4int , G4PhysicsVector* >::iterator itA = pointer_map->begin();
67  itA != pointer_map->end() ; ++itA ) {
68  G4PhysicsVector* pointerPhysicsVector = itA->second;
69  if ( pointerPhysicsVector ) {
70  delete pointerPhysicsVector;
71  itA->second = NULL;
72  }
73  }
74  delete pointer_map;
75  itZ->second = NULL;
76  }
77  }
78  mIsotope.clear();
79 }
80 
81 
82 
84 :G4VCrossSectionDataSet( "JENDLHE"+reaction+"CrossSection" )
85 {
86  reactionName = reaction;
87  BuildPhysicsTable( *pd );
88 }
89 
90 
91 
93 {
94  ;
95  //delete theCrossSections;
96 }
97 
98 
99 
101 {
102 
103 // if ( &aP != G4Neutron::Neutron() )
104 // throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
106 
107  G4String baseName = std::getenv( "G4NEUTRONHPDATA" );
108  G4String dirName = baseName+"/JENDL_HE/"+particleName+"/"+reactionName ;
109  G4String aFSType = "/CrossSection/";
110  G4ParticleHPNames theNames;
111 
113 
114 // Create JENDL_HE data
115 // Create map element or isotope
116 
117  size_t numberOfElements = G4Element::GetNumberOfElements();
118  //theCrossSections = new G4PhysicsTable( numberOfElements );
119 
120  // make a PhysicsVector for each element
121 
122  static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
123  vElement.clear();
124  vElement.resize( numberOfElements );
125  for ( size_t i = 0; i < numberOfElements; ++i )
126  {
127 
128  G4Element* theElement = (*theElementTable)[i];
129  vElement[i] = false;
130 
131  // isotope
132  G4int nIso = (*theElementTable)[i]->GetNumberOfIsotopes();
133  G4int Z = static_cast<G4int> ((*theElementTable)[i]->GetZ());
134  if ( nIso!=0 )
135  {
136  G4bool found_at_least_one = false;
137  for ( G4int i1 = 0; i1 < nIso; i1++ )
138  {
139  G4int A = theElement->GetIsotope(i1)->GetN();
140 
141  if ( isThisNewIsotope( Z , A ) )
142  {
143 
144  std::stringstream ss;
145  ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
146  filename = ss.str();
147  std::fstream file;
148  file.open ( filename , std::fstream::in );
149  G4int dummy;
150  file >> dummy;
151  if ( file.good() )
152  {
153 
154  //G4cout << "Found file for Z=" << Z << ", A=" << A << ", as " << filename << G4endl;
155  found_at_least_one = true;
156 
157  // read the file
158  G4PhysicsVector* aPhysVec = readAFile ( &file );
159 
160  //Regist
161 
162  registAPhysicsVector( Z , A , aPhysVec );
163 
164  }
165  else
166  {
167  //G4cout << "No file for "<< reactionType << " Z=" << Z << ", A=" << A << G4endl;
168  }
169 
170  file.close();
171 
172  }
173  else
174  {
175  found_at_least_one = TRUE;
176  }
177  }
178 
179  if ( found_at_least_one ) vElement[i] = true;
180 
181  }
182  else
183  {
184  G4StableIsotopes theStableOnes;
185  G4int first = theStableOnes.GetFirstIsotope( Z );
186  G4bool found_at_least_one = FALSE;
187  for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(theElement->GetZ() ) ); i1++)
188  {
189  G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
190 
191  if ( isThisNewIsotope( Z , A ) )
192  {
193 
194  std::stringstream ss;
195  ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
196  filename = ss.str();
197 
198  std::fstream file;
199  file.open ( filename , std::fstream::in );
200  G4int dummy;
201  file >> dummy;
202  if ( file.good() )
203  {
204  //G4cout << "Found file for Z=" << Z << ", A=" << A << ", as " << filename << G4endl;
205  found_at_least_one = TRUE;
206  //Read the file
207 
208  G4PhysicsVector* aPhysVec = readAFile ( &file );
209 
210  //Regist the PhysicsVector
211  registAPhysicsVector( Z , A , aPhysVec );
212 
213  }
214  else
215  {
216  //G4cout << "No file for "<< reactionType << " Z=" << Z << ", A=" << A << G4endl;
217  }
218 
219  file.close();
220  }
221  else
222  {
223  found_at_least_one = TRUE;
224  }
225  }
226 
227  if ( found_at_least_one ) vElement[i] = true;
228 
229  }
230 
231  }
232 
233 }
234 
235 
236 
238 {
239  if(&aP!=G4Neutron::Neutron())
240  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
241 // G4cout << "G4ParticleHPJENDLHEData::DumpPhysicsTable still to be implemented"<<G4endl;
242 }
243 
244 
245 
248 // aTemp
249 {
250 
251  // Primary energy >20MeV
252  // Thus
253  // Not take account of Doppler broadening
254  // also
255  // Not take account of Target thermal motions
256 
257  G4double result = 0;
258 
259  G4double ek = aP->GetKineticEnergy();
260 
261  G4int nIso = anE->GetNumberOfIsotopes();
262  G4int Z = static_cast<G4int> ( anE->GetZ() );
263  if ( nIso!=0 )
264  {
265  for ( G4int i1 = 0; i1 < nIso; i1++ )
266  {
267 
268  G4int A = anE->GetIsotope(i1)->GetN();
269  G4double frac = anE->GetRelativeAbundanceVector()[ i1 ]; // This case do NOT request "*perCent".
270 
271  result += frac * getXSfromThisIsotope( Z , A , ek );
272  //G4cout << reactionType << " XS in barn " << Z << " " << A << " " << frac << " " << getXSfromThisIsotope( Z , A , ek )/barn << G4endl;
273 
274  }
275  }
276  else
277  {
278 
279  G4StableIsotopes theStableOnes;
280  G4int first = theStableOnes.GetFirstIsotope( Z );
281  for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(anE->GetZ() ) ); i1++)
282  {
283 
284  G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
285  G4double frac = theStableOnes.GetAbundance( first+i1 )*perCent; // This case request "*perCent".
286 
287  result += frac * getXSfromThisIsotope( Z , A , ek );
288  //G4cout << reactionType << " XS in barn " << Z << " " << A << " " << frac << " " << getXSfromThisIsotope( Z , A , ek )/barn << G4endl;
289 
290  }
291  }
292  return result;
293 
294 }
295 
296 
297 
299 {
300 
301  G4int dummy;
302  G4int len;
303  *file >> dummy;
304  *file >> len;
305 
306  std::vector< G4double > v_e;
307  std::vector< G4double > v_xs;
308 
309  for ( G4int i = 0 ; i < len ; i++ )
310  {
311  G4double e;
312  G4double xs;
313 
314  *file >> e;
315  *file >> xs;
316  // data are written in eV and barn.
317  v_e.push_back( e*eV );
318  v_xs.push_back( xs*barn );
319  }
320 
321  G4LPhysicsFreeVector* aPhysVec = new G4LPhysicsFreeVector( static_cast< size_t >( len ) , v_e.front() , v_e.back() );
322 
323  for ( G4int i = 0 ; i < len ; i++ )
324  {
325  aPhysVec->PutValues( static_cast< size_t >( i ) , v_e[ i ] , v_xs[ i ] );
326  }
327 
328  return aPhysVec;
329 }
330 
331 
332 
334 {
335  if ( mIsotope.find ( z ) == mIsotope.end() ) return false;
336  if ( mIsotope.find ( z ) -> second->find ( a ) == mIsotope.find ( z ) -> second->end() ) return false;
337  return true;
338 }
339 
340 
341 
343 {
344 
345  std::pair< G4int , G4PhysicsVector* > aPair = std::pair < G4int , G4PhysicsVector* > ( A , aPhysVec );
346 
347  std::map < G4int , std::map< G4int , G4PhysicsVector* >* >::iterator itm;
348  itm = mIsotope.find ( Z );
349  if ( itm != mIsotope.end() )
350  {
351  itm->second->insert ( aPair );
352  }
353  else
354  {
355  std::map< G4int , G4PhysicsVector* >* aMap = new std::map< G4int , G4PhysicsVector* >;
356  aMap->insert ( aPair );
357  mIsotope.insert( std::pair< G4int , std::map< G4int , G4PhysicsVector* >* > ( Z , aMap ) );
358  }
359 
360 }
361 
362 
363 
365 {
366 
367  G4double aXSection = 0.0;
368  G4bool outOfRange;
369 
370  G4PhysicsVector* aPhysVec;
371  if ( mIsotope.find ( Z )->second->find ( A ) != mIsotope.find ( Z )->second->end() )
372  {
373 
374  aPhysVec = mIsotope.find ( Z )->second->find ( A )->second;
375  aXSection = aPhysVec->GetValue( ek , outOfRange );
376 
377  }
378  else
379  {
380 
381  //Select closest one in the same Z
382  std::map < G4int , G4PhysicsVector* >::iterator it;
383  G4int delta0 = 99; // no mean for 99
384  for ( it = mIsotope.find ( Z )->second->begin() ; it != mIsotope.find ( Z )->second->end() ; it++ )
385  {
386  G4int delta = std::abs( A - it->first );
387  if ( delta < delta0 ) delta0 = delta;
388  }
389 
390  // Randomize of selection larger or smaller than A
391  if ( G4UniformRand() < 0.5 ) delta0 *= -1;
392  G4int A1 = A + delta0;
393  if ( mIsotope.find ( Z )->second->find ( A1 ) != mIsotope.find ( Z )->second->end() )
394  {
395  aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
396  }
397  else
398  {
399  A1 = A - delta0;
400  aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
401  }
402 
403  aXSection = aPhysVec->GetValue( ek , outOfRange );
404  // X^(2/3) factor
405  //aXSection *= std::pow ( 1.0*A/ A1 , 2.0 / 3.0 );
406  aXSection *= G4Pow::GetInstance()->A23( 1.0*A/ A1 );
407 
408  }
409 
410  return aXSection;
411 }