Geant4_10
G4LENDCrossSection.cc
Go to the documentation of this file.
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 "G4SystemOfUnits.hh"
42 #include "G4ElementTable.hh"
43 #include "G4HadronicException.hh"
44 
45 //TK110811
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  return true;
54 }
55 
57  const G4Isotope* isotope , const G4Element* /*elment*/ , const G4Material* material )
58 {
59 
60  G4double xs = 0.0;
61  G4double ke = dp->GetKineticEnergy();
62  G4double temp = material->GetTemperature();
63  G4int iM = isotope->Getm();
64 
65  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA , iM ) )->second->GetTarget();
66 
67  xs = getLENDCrossSection ( aTarget , ke , temp );
68 
69  return xs;
70 }
71 
72 
73 /*
74 G4bool G4LENDCrossSection::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
75 {
76  G4bool result = true;
77  G4double eKin = aP->GetKineticEnergy();
78  if( eKin > GetMaxKinEnergy() || aP->GetDefinition() != proj ) result = false;
79  return result;
80 }
81 */
82 
85 {
86 
87  proj = NULL; //will be set in an inherited class
88  //default_evaluation = "endl99";
89  default_evaluation = "ENDF.B-VII.0";
90 
91  allow_nat = false;
92  allow_any = false;
93 
94  SetMinKinEnergy( 0*MeV );
95  SetMaxKinEnergy( 20*MeV );
96 
97  lend_manager = G4LENDManager::GetInstance();
98 
99 }
100 
102 {
103 
104  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
105  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
106  {
107  delete it->second;
108  }
109 
110 }
111 
113 {
115 }
116 
118 {
119 
120  if ( &aP != proj )
121  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use LEND data for particles other than neutrons!!!");
122 
123  G4cout << G4endl;
124  G4cout << "Dump Cross Sections of " << GetName() << G4endl;
125  G4cout << "(Pointwise cross-section at 300 Kelvin.)" << G4endl;
126  G4cout << G4endl;
127 
128  G4cout << "Target informaiton " << G4endl;
129 
130  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
131  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
132  {
133  G4cout
134  << "Wanted " << it->second->GetWantedEvaluation()
135  << ", Z= " << it->second->GetWantedZ()
136  << ", A= " << it->second->GetWantedA()
137  << "; Actual " << it->second->GetActualEvaluation()
138  << ", Z= " << it->second->GetActualZ()
139  << ", A= " << it->second->GetActualA()
140  << ", " << it->second->GetTarget()
141  << G4endl;
142 
143  G4int ie = 0;
144 
145  G4GIDI_target* aTarget = it->second->GetTarget();
146  G4double aT = 300;
147  for ( ie = 0 ; ie < 130 ; ie++ )
148  {
149  G4double ke = 1.0e-5 * std::pow ( 10.0 , ie/10.0 ) *eV;
150 
151  if ( ke < 20*MeV )
152  {
153  G4cout << " " << GetName() << ", cross section at " << ke/eV << " [eV] = " << getLENDCrossSection ( aTarget , ke , aT )/barn << " [barn] " << G4endl;
154  }
155  }
156  G4cout << G4endl;
157 
158  }
159 
160 }
161 
162 
163 /*
164 //110810
165 //G4double G4LENDCrossSection::GetCrossSection(const G4DynamicParticle* aP , const G4Element* anElement , G4double aT)
166 G4double G4LENDCrossSection::GetCrossSection(const G4DynamicParticle* aP , int iZ , const G4Material* aMat)
167 {
168 
169 //110810
170  G4double aT = aMat->GetTemperature();
171  G4Element* anElement = lend_manager->GetNistElementBuilder()->FindOrBuildElement( iZ );
172 
173  G4double ke = aP->GetKineticEnergy();
174  G4double XS = 0.0;
175 
176  G4int numberOfIsotope = anElement->GetNumberOfIsotopes();
177 
178  if ( numberOfIsotope > 0 )
179  {
180  // User Defined Abundances
181  for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
182  {
183 
184  G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
185  G4int iA = anElement->GetIsotope( i_iso )->GetN();
186  G4double ratio = anElement->GetRelativeAbundanceVector()[i_iso];
187 
188  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
189  XS += ratio*getLENDCrossSection ( aTarget , ke , aT );
190 
191  }
192  }
193  else
194  {
195  // Natural Abundances
196  G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder();
197  G4int iZ = int ( anElement->GetZ() );
198  G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) );
199 
200  G4int Nfirst = nistElementBuild->GetNistFirstIsotopeN( iZ );
201  for ( G4int i = 0 ; i < numberOfNistIso ; i++ )
202  {
203  G4int iA = Nfirst + i;
204  G4double ratio = nistElementBuild->GetIsotopeAbundance( iZ , iA );
205  if ( ratio > 0.0 )
206  {
207  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
208  XS += ratio*getLENDCrossSection ( aTarget , ke , aT );
209  //G4cout << ke/eV << " " << iZ << " " << iMass << " " << aTarget << " " << getLENDCrossSection ( aTarget , ke , aT ) << G4endl;
210  }
211  }
212  }
213 
214  //G4cout << "XS= " << XS << G4endl;
215  return XS;
216 }
217 
218 
219 
220 //110810
221 //G4double G4LENDCrossSection::GetIsoCrossSection(const G4DynamicParticle* dp, const G4Isotope* isotope, G4double aT )
222 G4double G4LENDCrossSection::GetIsoCrossSection(const G4DynamicParticle* dp, const G4Isotope* isotope, const G4Material* aMat)
223 {
224 
225 //110810
226  G4double aT = aMat->GetTemperature();
227 
228  G4double ke = dp->GetKineticEnergy();
229 
230  G4int iZ = isotope->GetZ();
231  G4int iA = isotope->GetN();
232 
233  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
234 
235  return getLENDCrossSection ( aTarget , ke , aT );
236 
237 }
238 
239 
240 
241 //110810
242 //G4double G4LENDCrossSection::GetZandACrossSection(const G4DynamicParticle* dp, G4int iZ, G4int iA, G4double aT)
243 G4double G4LENDCrossSection::GetZandACrossSection(const G4DynamicParticle* dp, G4int iZ, G4int iA, const G4Material* aMat)
244 {
245 
246 //110810
247  G4double aT = aMat->GetTemperature();
248 
249  G4double ke = dp->GetKineticEnergy();
250 
251  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
252 
253  return getLENDCrossSection ( aTarget , ke , aT );
254 
255 }
256 */
257 
258 
259 
260 void G4LENDCrossSection::recreate_used_target_map()
261 {
262  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
263  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
264  {
265  delete it->second;
266  }
267  usedTarget_map.clear();
268 
270 }
271 
272 
273 
275 {
276 
277  lend_manager->RequestChangeOfVerboseLevel( verboseLevel );
278 
279  size_t numberOfElements = G4Element::GetNumberOfElements();
280  static const G4ElementTable* theElementTable = G4Element::GetElementTable();
281 
282  for ( size_t i = 0 ; i < numberOfElements ; ++i )
283  {
284 
285  const G4Element* anElement = (*theElementTable)[i];
286  G4int numberOfIsotope = anElement->GetNumberOfIsotopes();
287 
288  if ( numberOfIsotope > 0 )
289  {
290  // User Defined Abundances
291  for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
292  {
293  G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
294  G4int iA = anElement->GetIsotope( i_iso )->GetN();
295  G4int iIsomer = anElement->GetIsotope( i_iso )->Getm();
296 
297  //G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( G4Neutron::Neutron() , default_evaluation , iZ , iA );
298  G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iA );
299  if ( allow_nat == true ) aTarget->AllowNat();
300  if ( allow_any == true ) aTarget->AllowAny();
301  usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iA , iIsomer ) , aTarget ) );
302  }
303  }
304  else
305  {
306  // Natural Abundances
307  G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder();
308  G4int iZ = int ( anElement->GetZ() );
309  //G4cout << nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ) << G4endl;
310  G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) );
311 
312  for ( G4int ii = 0 ; ii < numberOfNistIso ; ii++ )
313  {
314  //G4cout << nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + i ) << G4endl;
315  if ( nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii ) > 0 )
316  {
317  G4int iMass = nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii;
318  //G4cout << iZ << " " << nistElementBuild->GetNistFirstIsotopeN( iZ ) + i << " " << nistElementBuild->GetIsotopeAbundance ( iZ , iMass ) << G4endl;
319  G4int iIsomer = 0;
320 
321  G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iMass );
322  if ( allow_nat == true ) aTarget->AllowNat();
323  if ( allow_any == true ) aTarget->AllowAny();
324  usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iMass , iIsomer ) , aTarget ) );
325 
326  }
327 
328  }
329  }
330  }
331 
332  G4cout << "Dump UsedTarget for " << GetName() << G4endl;
333  G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) , Pointer of Target" << G4endl;
334  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
335  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
336  {
337  G4cout
338  << " " << it->second->GetWantedEvaluation()
339  << ", " << it->second->GetWantedZ()
340  << ", " << it->second->GetWantedA()
341  << " -> " << it->second->GetActualEvaluation()
342  << ", " << it->second->GetActualZ()
343  << ", " << it->second->GetActualA()
344  << ", " << it->second->GetTarget()
345  << G4endl;
346  }
347 
348 }
349 
350  // elow ehigh xs_elow xs_ehigh ke (ke < elow)
352 {
353  //XS propotinal to 1/v at low energy -> 1/root(E)
354  //XS = a * 1/root(E) + b
355  G4double a = ( y2 - y1 ) / ( 1/std::sqrt(x2) - 1/std::sqrt(x1) );
356  G4double b = y1 - a * 1/std::sqrt(x1);
357  G4double result = a * 1/std::sqrt(ke) + b;
358  return result;
359 }
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:158
Double_t y2[nxs]
Definition: Style.C:21
virtual G4double getLENDCrossSection(G4GIDI_target *, G4double, G4double)
void DumpPhysicsTable(const G4ParticleDefinition &)
G4double GetUltraLowEnergyExtrapolatedXS(G4double, G4double, G4double, G4double, G4double)
Double_t y1[nxs]
Definition: Style.C:20
Double_t x2[nxs]
Definition: Style.C:19
tuple a
Definition: test.py:11
G4double GetKineticEnergy() const
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4NistElementBuilder * GetNistElementBuilder()
G4double GetZ() const
Definition: G4Element.hh:131
G4double G4NeutronHPJENDLHEData::G4double result
G4double GetIsotopeAbundance(G4int Z, G4int N) const
G4int GetNumberOfNistIsotopes(G4int Z) const
G4ParticleDefinition * GetDefinition() const
G4ParticleDefinition * proj
const G4String & GetName() const
int G4int
Definition: G4Types.hh:78
G4LENDCrossSection(const G4String name="")
G4bool IsIsoApplicable(const G4DynamicParticle *, G4int, G4int, const G4Element *, const G4Material *)
G4int GetNucleusEncoding(G4int iZ, G4int iA, G4int iM)
string material
Definition: eplot.py:19
tuple b
Definition: test.py:12
G4int GetN() const
Definition: G4Isotope.hh:94
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:402
G4int Getm() const
Definition: G4Isotope.hh:100
void SetMinKinEnergy(G4double value)
bool G4bool
Definition: G4Types.hh:79
G4int GetNistFirstIsotopeN(G4int Z) const
Double_t x1[nxs]
Definition: Style.C:18
G4int GetZ() const
Definition: G4Isotope.hh:91
static G4LENDManager * GetInstance()
void SetMaxKinEnergy(G4double value)
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:169
G4double GetTemperature() const
Definition: G4Material.hh:180
#define G4endl
Definition: G4ios.hh:61
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
G4bool RequestChangeOfVerboseLevel(G4int)
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:395
G4double GetIsoCrossSection(const G4DynamicParticle *, G4int, G4int, const G4Isotope *, const G4Element *, const G4Material *)
void BuildPhysicsTable(const G4ParticleDefinition &)