ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LENDCrossSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LENDCrossSection.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 for LEND (Low Energy Nuclear Data)
28 // LEND is Geant4 interface for GIDI (General Interaction Data Interface)
29 // which gives a discription of nuclear and atomic reactions, such as
30 // Binary collision cross sections
31 // Particle number multiplicity distributions of reaction products
32 // Energy and angular distributions of reaction products
33 // Derived calculational constants
34 // GIDI is developped at Lawrence Livermore National Laboratory
35 // Class Description - End
36 
37 // 071025 First implementation done by T. Koi (SLAC/SCCS)
38 // 101118 Name modifications for release T. Koi (SLAC/PPA)
39 
40 #include "G4LENDCrossSection.hh"
41 #include "G4Pow.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4ElementTable.hh"
44 #include "G4HadronicException.hh"
45 
46 G4bool G4LENDCrossSection::IsIsoApplicable( const G4DynamicParticle* dp, G4int iZ , G4int iA ,
47  const G4Element* element , const G4Material* /*material*/ )
48 {
49  G4double eKin = dp->GetKineticEnergy();
50  if ( dp->GetDefinition() != proj ) return false;
51  if ( eKin > GetMaxKinEnergy() || eKin < GetMinKinEnergy() ) return false;
52 
53  //G4cout << "G4LENDCrossSection::GetIsoIsIsoApplicable this->GetName() = " << this->GetName() << ", iZ = " << iZ << ", iA = " << iA << ", allow_nat = " << allow_nat << G4endl;
54  //Check existence of target data
55  if ( element != NULL ) {
56  if ( element->GetNumberOfIsotopes() != 0 ) {
57  std::vector< const G4Isotope*> vIsotope;
58  for ( size_t i = 0 ; i != element->GetNumberOfIsotopes() ; i++ ) {
59  if ( element->GetIsotope( i )->GetN() == iA ) vIsotope.push_back( element->GetIsotope( i ) );
60  }
61  for ( size_t i = 0 ; i != vIsotope.size() ; i++ ) {
62  G4int iM = vIsotope[i]->Getm();
63  if ( get_target_from_map( lend_manager->GetNucleusEncoding( iZ , iA , iM ) ) != NULL ) return true;
64  }
65  //No isomer has data
66  //Check natural aboundance data for the element
67  if ( get_target_from_map( lend_manager->GetNucleusEncoding( iZ , 0 , 0 ) ) != NULL ) return true;
68  } else {
69  //Check for iZ and iA under assuming iM = 0
70  if ( get_target_from_map( lend_manager->GetNucleusEncoding( iZ , iA , 0 ) ) != NULL ) return true;
71  //Check natural aboundance data for the element
72  if ( get_target_from_map( lend_manager->GetNucleusEncoding( iZ , 0 , 0 ) ) != NULL ) return true;
73  }
74  } else {
75  //Check for iZ and iA under assuming iM = 0
76  if ( get_target_from_map( lend_manager->GetNucleusEncoding( iZ , iA , 0 ) ) != NULL ) return true;
77  //Check natural aboundance data for iZ
78  if ( get_target_from_map( lend_manager->GetNucleusEncoding( iZ , 0 , 0 ) ) != NULL ) return true;
79  }
80  return false;
81 }
82 
83 G4double G4LENDCrossSection::GetIsoCrossSection( const G4DynamicParticle* dp , G4int iZ , G4int iA ,
84  const G4Isotope* isotope , const G4Element* /*elment*/ , const G4Material* material )
85 {
86 
87  G4double xs = 0.0;
88  G4double ke = dp->GetKineticEnergy();
89  G4double temp = material->GetTemperature();
90  G4int iM = 0;
91  if ( isotope != NULL ) iM = isotope->Getm();
92 
93  G4GIDI_target* aTarget = get_target_from_map( lend_manager->GetNucleusEncoding( iZ , iA , iM ) );
94  if ( aTarget == NULL ) {
95  G4String message;
96  message = this->GetName();
97  message += " is unexpectedly called.";
98  //G4Exception( "G4LEND::GetIsoCrossSection(,)" , "LENDCrossSection-01" , JustWarning ,
99  G4Exception( "G4LEND::GetIsoCrossSection(,)" , "LENDCrossSection-01" , FatalException ,
100  message );
101  }
102  xs = getLENDCrossSection ( aTarget , ke , temp );
103 
104  return xs;
105 }
106 
107 
108 /*
109 G4bool G4LENDCrossSection::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
110 {
111  G4bool result = true;
112  G4double eKin = aP->GetKineticEnergy();
113  if( eKin > GetMaxKinEnergy() || aP->GetDefinition() != proj ) result = false;
114  return result;
115 }
116 */
117 
118 G4LENDCrossSection::G4LENDCrossSection( const G4String nam )
119 :G4VCrossSectionDataSet( nam )
120 {
121 
122  proj = NULL; //will be set in an inherited class
123  //default_evaluation = "endl99";
124  //default_evaluation = "ENDF.B-VII.0";
125  default_evaluation = "ENDF/BVII.1";
126 
127  allow_nat = false;
128  allow_any = false;
129 
130  SetMinKinEnergy( 0*MeV );
131  SetMaxKinEnergy( 20*MeV );
132 
133  lend_manager = G4LENDManager::GetInstance();
134 
135 }
136 
137 G4LENDCrossSection::~G4LENDCrossSection()
138 {
139 
140  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
141  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
142  {
143  delete it->second;
144  }
145 
146 }
147 
148 void G4LENDCrossSection::BuildPhysicsTable( const G4ParticleDefinition& )
149 {
150  create_used_target_map();
151 }
152 
153 void G4LENDCrossSection::DumpPhysicsTable(const G4ParticleDefinition& aP)
154 {
155 
156  if ( &aP != proj )
157  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use LEND data for particles other than neutrons!!!");
158 
159  G4cout << G4endl;
160  G4cout << "Dump Cross Sections of " << GetName() << G4endl;
161  G4cout << "(Pointwise cross-section at 300 Kelvin.)" << G4endl;
162  G4cout << G4endl;
163 
164  G4cout << "Target informaiton " << G4endl;
165 
166  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
167  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
168  {
169  G4cout
170  << "Wanted " << it->second->GetWantedEvaluation()
171  << ", Z= " << it->second->GetWantedZ()
172  << ", A= " << it->second->GetWantedA()
173  << "; Actual " << it->second->GetActualEvaluation()
174  << ", Z= " << it->second->GetActualZ()
175  << ", A= " << it->second->GetActualA()
176  << ", " << it->second->GetTarget()
177  << G4endl;
178 
179  G4int ie = 0;
180 
181  G4GIDI_target* aTarget = it->second->GetTarget();
182  G4double aT = 300;
183  for ( ie = 0 ; ie < 130 ; ie++ )
184  {
185  G4double ke = 1.0e-5 * G4Pow::GetInstance()->powA ( 10.0 , ie/10.0 ) *eV;
186 
187  if ( ke < 20*MeV )
188  {
189  G4cout << " " << GetName() << ", cross section at " << ke/eV << " [eV] = " << getLENDCrossSection ( aTarget , ke , aT )/barn << " [barn] " << G4endl;
190  }
191  }
192  G4cout << G4endl;
193 
194  }
195 
196 }
197 
198 
199 /*
200 //110810
201 //G4double G4LENDCrossSection::GetCrossSection(const G4DynamicParticle* aP , const G4Element* anElement , G4double aT)
202 G4double G4LENDCrossSection::GetCrossSection(const G4DynamicParticle* aP , int iZ , const G4Material* aMat)
203 {
204 
205 //110810
206  G4double aT = aMat->GetTemperature();
207  G4Element* anElement = lend_manager->GetNistElementBuilder()->FindOrBuildElement( iZ );
208 
209  G4double ke = aP->GetKineticEnergy();
210  G4double XS = 0.0;
211 
212  G4int numberOfIsotope = anElement->GetNumberOfIsotopes();
213 
214  if ( numberOfIsotope > 0 )
215  {
216  // User Defined Abundances
217  for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
218  {
219 
220  G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
221  G4int iA = anElement->GetIsotope( i_iso )->GetN();
222  G4double ratio = anElement->GetRelativeAbundanceVector()[i_iso];
223 
224  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
225  XS += ratio*getLENDCrossSection ( aTarget , ke , aT );
226 
227  }
228  }
229  else
230  {
231  // Natural Abundances
232  G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder();
233  G4int iZ = int ( anElement->GetZ() );
234  G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) );
235 
236  G4int Nfirst = nistElementBuild->GetNistFirstIsotopeN( iZ );
237  for ( G4int i = 0 ; i < numberOfNistIso ; i++ )
238  {
239  G4int iA = Nfirst + i;
240  G4double ratio = nistElementBuild->GetIsotopeAbundance( iZ , iA );
241  if ( ratio > 0.0 )
242  {
243  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
244  XS += ratio*getLENDCrossSection ( aTarget , ke , aT );
245  //G4cout << ke/eV << " " << iZ << " " << iMass << " " << aTarget << " " << getLENDCrossSection ( aTarget , ke , aT ) << G4endl;
246  }
247  }
248  }
249 
250  //G4cout << "XS= " << XS << G4endl;
251  return XS;
252 }
253 
254 
255 
256 //110810
257 //G4double G4LENDCrossSection::GetIsoCrossSection(const G4DynamicParticle* dp, const G4Isotope* isotope, G4double aT )
258 G4double G4LENDCrossSection::GetIsoCrossSection(const G4DynamicParticle* dp, const G4Isotope* isotope, const G4Material* aMat)
259 {
260 
261 //110810
262  G4double aT = aMat->GetTemperature();
263 
264  G4double ke = dp->GetKineticEnergy();
265 
266  G4int iZ = isotope->GetZ();
267  G4int iA = isotope->GetN();
268 
269  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
270 
271  return getLENDCrossSection ( aTarget , ke , aT );
272 
273 }
274 
275 
276 
277 //110810
278 //G4double G4LENDCrossSection::GetZandACrossSection(const G4DynamicParticle* dp, G4int iZ, G4int iA, G4double aT)
279 G4double G4LENDCrossSection::GetZandACrossSection(const G4DynamicParticle* dp, G4int iZ, G4int iA, const G4Material* aMat)
280 {
281 
282 //110810
283  G4double aT = aMat->GetTemperature();
284 
285  G4double ke = dp->GetKineticEnergy();
286 
287  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
288 
289  return getLENDCrossSection ( aTarget , ke , aT );
290 
291 }
292 */
293 
294 
295 
296 void G4LENDCrossSection::recreate_used_target_map()
297 {
298  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
299  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
300  {
301  delete it->second;
302  }
303  usedTarget_map.clear();
304 
305  create_used_target_map();
306 }
307 
308 
309 
310 void G4LENDCrossSection::create_used_target_map()
311 {
312 
313  lend_manager->RequestChangeOfVerboseLevel( verboseLevel );
314 
315  size_t numberOfElements = G4Element::GetNumberOfElements();
316  static const G4ElementTable* theElementTable = G4Element::GetElementTable();
317 
318  for ( size_t i = 0 ; i < numberOfElements ; ++i )
319  {
320 
321  const G4Element* anElement = (*theElementTable)[i];
322  G4int numberOfIsotope = anElement->GetNumberOfIsotopes();
323 
324  if ( numberOfIsotope > 0 )
325  {
326  // User Defined Abundances
327  for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
328  {
329  G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
330  G4int iA = anElement->GetIsotope( i_iso )->GetN();
331  G4int iIsomer = anElement->GetIsotope( i_iso )->Getm();
332 
333  //G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( G4Neutron::Neutron() , default_evaluation , iZ , iA );
334  G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iA , iIsomer );
335  if ( allow_nat == true ) aTarget->AllowNat();
336  if ( allow_any == true ) aTarget->AllowAny();
337  usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iA , iIsomer ) , aTarget ) );
338  }
339  }
340  else
341  {
342  // Natural Abundances
343  G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder();
344  G4int iZ = int ( anElement->GetZ() );
345  //G4cout << nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ) << G4endl;
346  G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) );
347 
348  for ( G4int ii = 0 ; ii < numberOfNistIso ; ii++ )
349  {
350  //G4cout << nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + i ) << G4endl;
351  if ( nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii ) > 0 )
352  {
353  G4int iMass = nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii;
354  //G4cout << iZ << " " << nistElementBuild->GetNistFirstIsotopeN( iZ ) + i << " " << nistElementBuild->GetIsotopeAbundance ( iZ , iMass ) << G4endl;
355  G4int iIsomer = 0;
356 
357  G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iMass );
358  if ( allow_nat == true ) aTarget->AllowNat();
359  if ( allow_any == true ) aTarget->AllowAny();
360  usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iMass , iIsomer ) , aTarget ) );
361 
362  }
363 
364  }
365  }
366  }
367  DumpLENDTargetInfo();
368 }
369 
370  // elow ehigh xs_elow xs_ehigh ke (ke < elow)
371 G4double G4LENDCrossSection::GetUltraLowEnergyExtrapolatedXS( G4double x1, G4double x2, G4double y1, G4double y2 , G4double ke )
372 {
373  //XS propotinal to 1/v at low energy -> 1/root(E)
374  //XS = a * 1/root(E) + b
375  G4double a = ( y2 - y1 ) / ( 1/std::sqrt(x2) - 1/std::sqrt(x1) );
376  G4double b = y1 - a * 1/std::sqrt(x1);
377  G4double result = a * 1/std::sqrt(ke) + b;
378  return result;
379 }
380 
381 G4GIDI_target* G4LENDCrossSection::get_target_from_map( G4int nuclear_code ) {
382  G4GIDI_target* target = NULL;
383  if ( usedTarget_map.find( nuclear_code ) != usedTarget_map.end() ) {
384  target = usedTarget_map.find( nuclear_code )->second->GetTarget();
385  }
386  return target;
387 }
388 
389 void G4LENDCrossSection::DumpLENDTargetInfo( G4bool force ) {
390 
391  if ( lend_manager->GetVerboseLevel() >= 1 || force ) {
392  if ( usedTarget_map.size() == 0 ) create_used_target_map();
393  G4cout << "Dumping UsedTarget of " << GetName() << " for " << proj->GetParticleName() << G4endl;
394  G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) " << G4endl;
395  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
396  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ ) {
397  G4cout
398  << " " << it->second->GetWantedEvaluation()
399  << ", " << it->second->GetWantedZ()
400  << ", " << it->second->GetWantedA()
401  << " -> " << it->second->GetActualEvaluation()
402  << ", " << it->second->GetActualZ()
403  << ", " << it->second->GetActualA()
404  << G4endl;
405  }
406  }
407 }
408