Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/G4LENDCrossSection.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 ( G4int i = 0 ; i != (G4int)element->GetNumberOfIsotopes() ; ++i ) {
 59             if ( element->GetIsotope( i )->GetN() == iA ) vIsotope.push_back( element->GetIsotope( i ) ); 
 60          }
 61          for ( std::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    std::size_t numberOfElements = G4Element::GetNumberOfElements();
316    static const G4ElementTable* theElementTable = G4Element::GetElementTable();
317 
318    for ( std::size_t i = 0 ; i < numberOfElements ; ++i )
319    {
320 
321       const G4Element* anElement = (*theElementTable)[i];
322       G4int numberOfIsotope = (G4int)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 
409