Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (16 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/hadronic/cross_sections/src
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/cross_sections/src/G4BGGNucleonElasticXS.cc

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BGGNucleonElasticXS.cc,v 1.1 2007/03/13 15:19:30 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4BGGNucleonElasticXS.cc,v 1.3 2008/12/01 16:50:23 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646#include "G4GlauberGribovCrossSection.hh"
    4747#include "G4NucleonNuclearCrossSection.hh"
     48#include "G4HadronNucleonXsc.hh"
    4849#include "G4Proton.hh"
    4950#include "G4Neutron.hh"
    … …  
    5253//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    5354
    54 G4BGGNucleonElasticXS::G4BGGNucleonElasticXS(const G4ParticleDefinition* p)
     55G4BGGNucleonElasticXS::G4BGGNucleonElasticXS(const G4ParticleDefinition*)
    5556{
    5657  verboseLevel = 0;
    57   thEnergy     = 100.*GeV;
    58   if(p == G4Proton::Proton() || p == G4Neutron::Neutron()) {
    59     fNucleon = new G4NucleonNuclearCrossSection();
    60     fGlauber = new G4GlauberGribovCrossSection();
    61     particle = p;
    62     Initialise();
    63   } else {
    64     fNucleon = 0;
    65     fGlauber = 0;
    66     particle = 0;
    67     if(p) G4cout << "### G4BGGNucleonElasticXS WARNING: is not applicable to "
    68                  << p->GetParticleName()
    69                  << G4endl;
    70     else  G4cout << "### G4BGGNucleonElasticXS WARNING: particle is not defined "
    71                  << G4endl;
    72   }
     58  fGlauberEnergy = 91.*GeV;
     59  fLowEnergy = 20.*MeV;
     60  fNucleon = 0;
     61  fGlauber = 0;
     62  fHadron  = 0;
     63  particle = 0;
     64  isProton = false;
     65  isInitialized = false;
    7366}
    7467
    … …  
    7972  delete fGlauber;
    8073  delete fNucleon;
     74  delete fHadron;
    8175}
    8276
    … …  
    9084  G4double cross = 0.0;
    9185  G4double ekin = dp->GetKineticEnergy();
    92   G4int iz = G4int(Z + 0.5);
     86  G4int iz = G4int(Z);
    9387  if(iz > 92) iz = 92;
    9488
    95   if(ekin > thEnergy) {
    96     cross = theFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
     89  if(ekin <= fLowEnergy) {
     90    cross = theCoulombFac[iz];
     91    if(isProton) { cross *= CoulombFactor(ekin, A); }
     92
     93  } else if(iz == 1) {
     94    if( A < 1.5) {
     95      //fHadron->GetHadronNucleonXscPDG(dp, G4Proton::Proton());
     96      //fHadron->GetHadronNucleonXscEL(dp, G4Proton::Proton());
     97      //fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
     98      fHadron->GetHadronNucleonXscMK(dp, G4Proton::Proton());
     99      cross = fHadron->GetElasticHadronNucleonXsc();
     100    } else {
     101      fHadron->GetHadronNucleonXscMK(dp, G4Proton::Proton());
     102      cross = fHadron->GetElasticHadronNucleonXsc();
     103      fHadron->GetHadronNucleonXscMK(dp, G4Neutron::Neutron());
     104      cross += fHadron->GetElasticHadronNucleonXsc();
     105    }
     106  } else if(ekin > fGlauberEnergy) {
     107    cross = theGlauberFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
    97108  } else {
    98     cross = fNucleon->GetIsoZACrossSection(dp, Z, A);
     109    cross = fNucleon->GetElasticCrossSection(dp, Z, A);
    99110  }
    100111
    … …  
    112123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    113124
    114 void G4BGGNucleonElasticXS::BuildPhysicsTable(const G4ParticleDefinition&)
    115 {
     125void G4BGGNucleonElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
     126{
     127  if(&p == G4Proton::Proton() || &p == G4Neutron::Neutron()) {
     128    particle = &p;
     129    Initialise();
     130  } else {
     131    G4cout << "### G4BGGNucleonElasticXS WARNING: is not applicable to "
     132           << p.GetParticleName()
     133           << G4endl;
     134  }
    116135}
    117136
    … …  
    127146void G4BGGNucleonElasticXS::Initialise()
    128147{
     148  if(isInitialized) return;
     149  isInitialized = true;
     150
     151  fNucleon = new G4NucleonNuclearCrossSection();
     152  fGlauber = new G4GlauberGribovCrossSection();
     153  fHadron  = new G4HadronNucleonXsc();
     154  fNucleon->BuildPhysicsTable(*particle);
     155  fGlauber->BuildPhysicsTable(*particle);
     156  if(particle == G4Proton::Proton()) isProton = true;
     157
    129158  G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
    130159  G4ThreeVector mom(0.0,0.0,1.0);
    131   G4DynamicParticle dp(part, mom, thEnergy);
     160  G4DynamicParticle dp(part, mom, fGlauberEnergy);
    132161
    133162  G4NistManager* nist = G4NistManager::Instance();
    134   G4double A = nist->GetAtomicMassAmu(2);
    135 
    136   G4double csup, csdn;
     163
     164  G4double A, csup, csdn;
    137165
    138166  if(verboseLevel > 0) G4cout << "### G4BGGNucleonElasticXS::Initialise for "
    … …  
    145173
    146174    csup = fGlauber->GetElasticGlauberGribov(&dp, Z, A);
    147     csdn = fNucleon->GetIsoZACrossSection(&dp, Z, A);
    148 
    149     theFac[iz] = csdn/csup;
     175    csdn = fNucleon->GetElasticCrossSection(&dp, Z, A);
     176
     177    theGlauberFac[iz] = csdn/csup;
    150178    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
    151                                 << " factor= " << theFac[iz] << G4endl;
    152   }
    153 }
    154 
    155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    156 
    157 
     179                                << " factor= " << theGlauberFac[iz] << G4endl;
     180  }
     181  dp.SetKineticEnergy(fLowEnergy);
     182  fHadron->GetHadronNucleonXscMK(&dp, G4Proton::Proton());
     183  theCoulombFac[1] = fHadron->GetElasticHadronNucleonXsc();
     184  if(isProton) { theCoulombFac[1] /= CoulombFactor(fLowEnergy, 1.0); }
     185
     186  for(G4int iz=2; iz<93; iz++) {
     187
     188    G4double Z = G4double(iz);
     189    A = nist->GetAtomicMassAmu(iz);
     190
     191    theCoulombFac[iz] = fNucleon->GetElasticCrossSection(&dp, Z, A);
     192    if(isProton) { theCoulombFac[iz] /= CoulombFactor(fLowEnergy, A); }
     193    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
     194                                << " factor= " << theCoulombFac[iz] << G4endl;
     195  }
     196}
     197
     198//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     199
     200G4double G4BGGNucleonElasticXS::CoulombFactor(G4double kinEnergy, G4double A)
     201{
     202  G4double res= 0.0;
     203  if(kinEnergy <= DBL_MIN) return res;
     204  else if(A < 1.5) return kinEnergy*kinEnergy;
     205
     206  G4double elog = std::log10(kinEnergy/GeV);
     207
     208  // from G4ProtonInelasticCrossSection
     209  G4double f1 = 8.0  - 8.0/A - 0.008*A;
     210  G4double f2 = 2.34 - 5.4/A - 0.0028*A;
     211
     212  res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
     213 
     214  f1 = 5.6 - 0.016*A;
     215  f2 = 1.37 + 1.37/A;
     216  res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
     217  return res; 
     218}
     219
     220//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     221
     222

  • trunk/source/processes/hadronic/cross_sections/src/G4BGGNucleonInelasticXS.cc

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BGGNucleonInelasticXS.cc,v 1.1 2007/03/13 15:19:30 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4BGGNucleonInelasticXS.cc,v 1.3 2008/12/01 16:50:23 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646#include "G4GlauberGribovCrossSection.hh"
    4747#include "G4NucleonNuclearCrossSection.hh"
     48#include "G4HadronNucleonXsc.hh"
    4849#include "G4Proton.hh"
    4950#include "G4Neutron.hh"
    … …  
    5556{
    5657  verboseLevel = 0;
    57   thEnergy     = 100.*GeV;
    58   if(p == G4Proton::Proton() || p == G4Neutron::Neutron()) {
    59     fNucleon = new G4NucleonNuclearCrossSection();
    60     fGlauber = new G4GlauberGribovCrossSection();
    61     particle = p;
    62     Initialise();
    63   } else {
    64     fNucleon = 0;
    65     fGlauber = 0;
    66     particle = 0;
    67     if(p) G4cout << "### G4BGGNucleonInelasticXS WARNING: is not applicable to "
    68                  << p->GetParticleName()
    69                  << G4endl;
    70     else  G4cout << "### G4BGGNucleonInelasticXS WARNING: particle is not defined "
    71                  << G4endl;
    72   }
     58  fGlauberEnergy = 91.*GeV;
     59  fLowEnergy = 20.*MeV;
     60  fNucleon = 0;
     61  fGlauber = 0;
     62  fHadron  = 0;
     63  particle = p;
     64  isProton = false;
     65  isInitialized = false;
    7366}
    7467
    … …  
    7972  delete fGlauber;
    8073  delete fNucleon;
     74  delete fHadron;
    8175}
    8276
    … …  
    9084  G4double cross = 0.0;
    9185  G4double ekin = dp->GetKineticEnergy();
    92   G4int iz = G4int(Z + 0.5);
     86  G4int iz = G4int(Z);
    9387  if(iz > 92) iz = 92;
    9488
    95   if(ekin > thEnergy) {
    96     cross = theFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
     89  if(ekin <= fLowEnergy) {
     90    cross = theCoulombFac[iz]*CoulombFactor(ekin, A);
     91  } else if(iz == 1) {
     92    if( A < 1.5) {
     93      //fHadron->GetHadronNucleonXscPDG(dp, G4Proton::Proton());
     94      //fHadron->GetHadronNucleonXscEL(dp, G4Proton::Proton());
     95      //fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
     96      //fHadron->GetHadronNucleonXscVU(dp, G4Proton::Proton());
     97      fHadron->GetHadronNucleonXscMK(dp, G4Proton::Proton());
     98      cross = fHadron->GetInelasticHadronNucleonXsc();
     99    } else {
     100      fHadron->GetHadronNucleonXscMK(dp, G4Proton::Proton());
     101      cross = fHadron->GetInelasticHadronNucleonXsc();
     102      fHadron->GetHadronNucleonXscMK(dp, G4Neutron::Neutron());
     103      cross += fHadron->GetInelasticHadronNucleonXsc();
     104    }
     105  } else if(ekin > fGlauberEnergy) {
     106    cross = theGlauberFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
    97107  } else {
    98108    cross = fNucleon->GetIsoZACrossSection(dp, Z, A);
    … …  
    112122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    113123
    114 void G4BGGNucleonInelasticXS::BuildPhysicsTable(const G4ParticleDefinition&)
    115 {
     124void G4BGGNucleonInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
     125{
     126  if(&p == G4Proton::Proton() || &p == G4Neutron::Neutron()) {
     127    particle = &p;
     128    Initialise();
     129  } else {
     130    G4cout << "### G4BGGNucleonInelasticXS WARNING: is not applicable to "
     131           << p.GetParticleName()
     132           << G4endl;
     133  }
    116134}
    117135
    … …  
    127145void G4BGGNucleonInelasticXS::Initialise()
    128146{
     147  if(isInitialized) return;
     148  isInitialized = true;
     149
     150  fNucleon = new G4NucleonNuclearCrossSection();
     151  fGlauber = new G4GlauberGribovCrossSection();
     152  fHadron  = new G4HadronNucleonXsc();
     153  fNucleon->BuildPhysicsTable(*particle);
     154  fGlauber->BuildPhysicsTable(*particle);
     155  if(particle == G4Proton::Proton()) isProton = true;
     156
    129157  G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
    130158  G4ThreeVector mom(0.0,0.0,1.0);
    131   G4DynamicParticle dp(part, mom, thEnergy);
     159  G4DynamicParticle dp(part, mom, fGlauberEnergy);
    132160
    133161  G4NistManager* nist = G4NistManager::Instance();
    … …  
    147175    csdn = fNucleon->GetIsoZACrossSection(&dp, Z, A);
    148176
    149     theFac[iz] = csdn/csup;
     177    theGlauberFac[iz] = csdn/csup;
    150178    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
    151                                 << " factor= " << theFac[iz] << G4endl;
    152   }
    153 }
    154 
    155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    156 
    157 
     179                                << " factor= " << theGlauberFac[iz] << G4endl;
     180  }
     181  dp.SetKineticEnergy(fLowEnergy);
     182  fHadron->GetHadronNucleonXscMK(&dp, G4Proton::Proton());
     183  theCoulombFac[1] =
     184    fHadron->GetInelasticHadronNucleonXsc()/CoulombFactor(fLowEnergy,1.0);
     185     
     186  for(G4int iz=2; iz<93; iz++) {
     187
     188    G4double Z = G4double(iz);
     189    A = nist->GetAtomicMassAmu(iz);
     190
     191    theCoulombFac[iz] =
     192      fNucleon->GetIsoZACrossSection(&dp, Z, A)/CoulombFactor(fLowEnergy,A);
     193
     194    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
     195                                << " factor= " << theCoulombFac[iz] << G4endl;
     196  }
     197}
     198
     199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     200
     201G4double G4BGGNucleonInelasticXS::CoulombFactor(G4double kinEnergy, G4double A)
     202{
     203  G4double res= 0.0;
     204  if(kinEnergy <= DBL_MIN) return res;
     205  else if (A < 1.5) return kinEnergy*kinEnergy;
     206 
     207  G4double elog = std::log10(kinEnergy/GeV);
     208
     209  // from G4ProtonInelasticCrossSection
     210  if(isProton) {
     211    G4double f1 = 8.0  - 8.0/A - 0.008*A;
     212    G4double f2 = 2.34 - 5.4/A - 0.0028*A;
     213
     214    res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
     215 
     216    f1 = 5.6 - 0.016*A;
     217    f2 = 1.37 + 1.37/A;
     218    res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
     219  } else {
     220
     221    G4double p3 = 0.6 + 13./A - 0.0005*A;
     222    G4double p4 = 7.2449 - 0.018242*A;
     223    G4double p5 = 1.36 + 1.8/A + 0.0005*A;
     224    G4double p6 = 1. + 200./A + 0.02*A;
     225    G4double p7 = 3.0 - (A-70.)*(A-200.)/11000.;
     226
     227    res = (1.+p3/(1. + std::exp(p4*(elog+p5))))/(1.+std::exp(-p6*(elog+p7)));
     228
     229  }
     230  return res; 
     231}
     232
     233//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     234

  • trunk/source/processes/hadronic/cross_sections/src/G4BGGPionElasticXS.cc

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BGGPionElasticXS.cc,v 1.1 2007/03/13 15:19:30 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4BGGPionElasticXS.cc,v 1.3 2008/12/01 16:50:23 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646#include "G4GlauberGribovCrossSection.hh"
    4747#include "G4UPiNuclearCrossSection.hh"
     48#include "G4HadronNucleonXsc.hh"
    4849#include "G4PionPlus.hh"
    4950#include "G4PionMinus.hh"
    … …  
    5253//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    5354
    54 G4BGGPionElasticXS::G4BGGPionElasticXS(const G4ParticleDefinition* p)
     55G4BGGPionElasticXS::G4BGGPionElasticXS(const G4ParticleDefinition*)
    5556{
    5657  verboseLevel = 0;
    57   thEnergy     = 100.*GeV;
    58   if(p == G4PionPlus::PionPlus() || p == G4PionMinus::PionMinus()) {
    59     fPion = new G4UPiNuclearCrossSection();
    60     fGlauber = new G4GlauberGribovCrossSection();
    61     particle = p;
    62     Initialise();
    63   } else {
    64     fPion = 0;
    65     fGlauber = 0;
    66     particle = 0;
    67     if(p) G4cout << "### G4BGGPionElasticXS WARNING: is not applicable to "
    68                  << p->GetParticleName()
    69                  << G4endl;
    70     else  G4cout << "### G4BGGPionElasticXS WARNING: particle is not defined "
    71                  << G4endl;
    72   }
     58  fGlauberEnergy = 91.*GeV;
     59  fLowEnergy = 20.*MeV;
     60  fPion = 0;
     61  fGlauber = 0;
     62  fHadron  = 0;
     63  particle = 0;
     64  isPiplus = false;
     65  isInitialized = false;
    7366}
    7467
    … …  
    7972  delete fGlauber;
    8073  delete fPion;
     74  delete fHadron;
    8175}
    8276
    … …  
    9084  G4double cross = 0.0;
    9185  G4double ekin = dp->GetKineticEnergy();
    92   G4int iz = G4int(Z + 0.5);
     86  G4int iz = G4int(Z);
    9387  if(iz > 92) iz = 92;
    9488
    95   if(ekin > thEnergy) {
    96     cross = theFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
     89  if(ekin <= fLowEnergy) {
     90    cross = theCoulombFac[iz];
     91    if(isPiplus) { cross *= CoulombFactor(ekin, A); }
     92
     93  } else if(iz == 1) {
     94    if( A < 1.5) {
     95      //fHadron->GetHadronNucleonXscPDG(dp, G4Proton::Proton());
     96      //fHadron->GetHadronNucleonXscEL(dp, G4Proton::Proton());
     97      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
     98      //fHadron->GetHadronNucleonXscVU(dp, G4Proton::Proton());
     99      //fHadron->GetHadronNucleonXscMK(dp, G4Proton::Proton());
     100      cross = fHadron->GetElasticHadronNucleonXsc();
     101    } else {
     102      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
     103      cross = fHadron->GetElasticHadronNucleonXsc();
     104      fHadron->GetHadronNucleonXscNS(dp, G4Neutron::Neutron());
     105      cross += fHadron->GetElasticHadronNucleonXsc();
     106    }
     107  } else if(ekin > fGlauberEnergy) {
     108    cross = theGlauberFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
    97109  } else {
    98110    cross = fPion->GetElasticCrossSection(dp, Z, A);
    … …  
    112124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    113125
    114 void G4BGGPionElasticXS::BuildPhysicsTable(const G4ParticleDefinition&)
    115 {
     126void G4BGGPionElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
     127{
     128  if(&p == G4PionPlus::PionPlus() || &p == G4PionMinus::PionMinus()) {
     129    particle = &p;
     130    Initialise();
     131  } else {
     132    G4cout << "### G4BGGPionElasticXS WARNING: is not applicable to "
     133           << p.GetParticleName()
     134           << G4endl;
     135  }
    116136}
    117137
    … …  
    127147void G4BGGPionElasticXS::Initialise()
    128148{
     149  if(isInitialized) return;
     150  isInitialized = true;
     151
     152  fPion = new G4UPiNuclearCrossSection();
     153  fGlauber = new G4GlauberGribovCrossSection();
     154  fHadron  = new G4HadronNucleonXsc();
     155  fPion->BuildPhysicsTable(*particle);
     156  fGlauber->BuildPhysicsTable(*particle);
     157  if(particle == G4PionPlus::PionPlus()) isPiplus = true;
     158
    129159  G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
    130160  G4ThreeVector mom(0.0,0.0,1.0);
    131   G4DynamicParticle dp(part, mom, thEnergy);
     161  G4DynamicParticle dp(part, mom, fGlauberEnergy);
    132162
    133163  G4NistManager* nist = G4NistManager::Instance();
    134   G4double A = nist->GetAtomicMassAmu(2);
    135 
    136   G4double csup, csdn;
     164
     165  G4double A, csup, csdn;
    137166
    138167  if(verboseLevel > 0) G4cout << "### G4BGGPionElasticXS::Initialise for "
    … …  
    147176    csdn = fPion->GetElasticCrossSection(&dp, Z, A);
    148177
    149     theFac[iz] = csdn/csup;
     178    theGlauberFac[iz] = csdn/csup;
    150179    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
    151                                 << " factor= " << theFac[iz] << G4endl;
    152   }
    153 }
    154 
    155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    156 
    157 
     180                                << " factor= " << theGlauberFac[iz] << G4endl;
     181  }
     182  dp.SetKineticEnergy(fLowEnergy);
     183  fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
     184  theCoulombFac[1] = fHadron->GetElasticHadronNucleonXsc();
     185  if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1.0); }
     186
     187  for(G4int iz=2; iz<93; iz++) {
     188
     189    G4double Z = G4double(iz);
     190    A = nist->GetAtomicMassAmu(iz);
     191
     192    theCoulombFac[iz] = fPion->GetElasticCrossSection(&dp, Z, A);
     193    if(isPiplus) { theCoulombFac[iz] /= CoulombFactor(fLowEnergy,A); }
     194    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
     195                                << " factor= " << theCoulombFac[iz] << G4endl;
     196  }
     197}
     198
     199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     200
     201G4double G4BGGPionElasticXS::CoulombFactor(G4double kinEnergy, G4double A)
     202{
     203  G4double res= 0.0;
     204  if(kinEnergy <= DBL_MIN) return res;
     205  else if(A < 1.5) return kinEnergy*kinEnergy;
     206 
     207  G4double elog = std::log10(kinEnergy/GeV);
     208
     209  // from G4ProtonInelasticCrossSection
     210  G4double f1 = 8.0  - 8.0/A - 0.008*A;
     211  G4double f2 = 2.34 - 5.4/A - 0.0028*A;
     212
     213  res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
     214 
     215  f1 = 5.6 - 0.016*A;
     216  f2 = 1.37 + 1.37/A;
     217  res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
     218  return res; 
     219}
     220
     221//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     222
     223

  • trunk/source/processes/hadronic/cross_sections/src/G4BGGPionInelasticXS.cc

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BGGPionInelasticXS.cc,v 1.1 2007/03/13 15:19:30 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4BGGPionInelasticXS.cc,v 1.3 2008/12/01 16:50:23 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646#include "G4GlauberGribovCrossSection.hh"
    4747#include "G4UPiNuclearCrossSection.hh"
     48#include "G4HadronNucleonXsc.hh"
    4849#include "G4PionPlus.hh"
    4950#include "G4PionMinus.hh"
    … …  
    5556{
    5657  verboseLevel = 0;
    57   thEnergy     = 100.*GeV;
    58   if(p == G4PionPlus::PionPlus() || p == G4PionMinus::PionMinus()) {
    59     fPion = new G4UPiNuclearCrossSection();
    60     fGlauber = new G4GlauberGribovCrossSection();
    61     particle = p;
    62     Initialise();
    63   } else {
    64     fPion = 0;
    65     fGlauber = 0;
    66     particle = 0;
    67     if(p) G4cout << "### G4BGGPionInelasticXS WARNING: is not applicable to "
    68                  << p->GetParticleName()
    69                  << G4endl;
    70     else  G4cout << "### G4BGGPionInelasticXS WARNING: particle is not defined "
    71                  << G4endl;
    72   }
     58  fGlauberEnergy = 91.*GeV;
     59  fLowEnergy = 20.*MeV;
     60  fPion = 0;
     61  fGlauber = 0;
     62  fHadron  = 0;
     63  particle = p;
     64  isPiplus = false;
     65  isInitialized = false;
    7366}
    7467
    … …  
    7972  delete fGlauber;
    8073  delete fPion;
     74  delete fHadron;
    8175}
    8276
    … …  
    9084  G4double cross = 0.0;
    9185  G4double ekin = dp->GetKineticEnergy();
    92   G4int iz = G4int(Z + 0.5);
     86  G4int iz = G4int(Z);
    9387  if(iz > 92) iz = 92;
    9488
    95   if(ekin > thEnergy) {
    96     cross = theFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
     89  if(ekin <= fLowEnergy) {
     90    cross = theCoulombFac[iz];
     91    if(isPiplus) { cross *= CoulombFactor(ekin, A); }
     92
     93  } else if(iz == 1) {
     94
     95    if( A < 1.5) {
     96      //fHadron->GetHadronNucleonXscPDG(dp, G4Proton::Proton());
     97      //fHadron->GetHadronNucleonXscEL(dp, G4Proton::Proton());
     98      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
     99      //fHadron->GetHadronNucleonXscVU(dp, G4Proton::Proton());
     100      //fHadron->GetHadronNucleonXscMK(dp, G4Proton::Proton());
     101      cross = fHadron->GetInelasticHadronNucleonXsc();
     102    } else {
     103      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
     104      cross = fHadron->GetInelasticHadronNucleonXsc();
     105      fHadron->GetHadronNucleonXscNS(dp, G4Neutron::Neutron());
     106      cross += fHadron->GetInelasticHadronNucleonXsc();
     107    }
     108
     109  } else if(ekin > fGlauberEnergy) {
     110    cross = theGlauberFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
     111
    97112  } else {
    98     cross = fPion->GetInelasticCrossSection(dp, Z, A);
     113    cross = fPion->GetIsoZACrossSection(dp, Z, A);
    99114  }
    100115
    … …  
    112127//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    113128
    114 void G4BGGPionInelasticXS::BuildPhysicsTable(const G4ParticleDefinition&)
    115 {
     129void G4BGGPionInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
     130{
     131  if(&p == G4PionPlus::PionPlus() || &p == G4PionMinus::PionMinus()) {
     132    particle = &p;
     133    Initialise();
     134  } else {
     135    G4cout << "### G4BGGPionInelasticXS WARNING: is not applicable to "
     136           << p.GetParticleName()
     137           << G4endl;
     138  }
    116139}
    117140
    … …  
    127150void G4BGGPionInelasticXS::Initialise()
    128151{
     152  if(isInitialized) return;
     153  isInitialized = true;
     154
     155  fPion = new G4UPiNuclearCrossSection();
     156  fGlauber = new G4GlauberGribovCrossSection();
     157  fHadron  = new G4HadronNucleonXsc();
     158  fPion->BuildPhysicsTable(*particle);
     159  fGlauber->BuildPhysicsTable(*particle);
     160  if(particle == G4PionPlus::PionPlus()) isPiplus = true;
     161
    129162  G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
    130163  G4ThreeVector mom(0.0,0.0,1.0);
    131   G4DynamicParticle dp(part, mom, thEnergy);
     164  G4DynamicParticle dp(part, mom, fGlauberEnergy);
    132165
    133166  G4NistManager* nist = G4NistManager::Instance();
    134   G4double A = nist->GetAtomicMassAmu(2);
    135 
    136   G4double csup, csdn;
     167
     168  G4double A, csup, csdn;
    137169
    138170  if(verboseLevel > 0) G4cout << "### G4BGGPionInelasticXS::Initialise for "
    139                               << particle->GetParticleName() << G4endl;
     171                              << particle->GetParticleName()
     172                              << " isPiplus: " << isPiplus
     173                              << G4endl;
    140174
    141175  for(G4int iz=2; iz<93; iz++) {
    … …  
    147181    csdn = fPion->GetInelasticCrossSection(&dp, Z, A);
    148182
    149     theFac[iz] = csdn/csup;
     183    theGlauberFac[iz] = csdn/csup;
    150184    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
    151                                 << " factor= " << theFac[iz] << G4endl;
    152   }
    153 }
    154 
    155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    156 
    157 
     185                                << " factor= " << theGlauberFac[iz] << G4endl;
     186  }
     187  dp.SetKineticEnergy(fLowEnergy);
     188  fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
     189  theCoulombFac[1] = fHadron->GetInelasticHadronNucleonXsc();
     190  if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1.0); }
     191     
     192  for(G4int iz=2; iz<93; iz++) {
     193
     194    G4double Z = G4double(iz);
     195    A = nist->GetAtomicMassAmu(iz);
     196
     197    theCoulombFac[iz] = fPion->GetIsoZACrossSection(&dp, Z, A);
     198    if(isPiplus) { theCoulombFac[iz] /= CoulombFactor(fLowEnergy,A); }
     199
     200    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
     201                                << " factor= " << theCoulombFac[iz] << G4endl;
     202  }
     203}
     204
     205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     206
     207G4double G4BGGPionInelasticXS::CoulombFactor(G4double kinEnergy, G4double A)
     208{
     209  G4double res= 0.0;
     210  if(kinEnergy <= DBL_MIN) return res;
     211  else if(A < 1.5) return kinEnergy*kinEnergy;
     212 
     213  G4double elog = std::log10(kinEnergy/GeV);
     214
     215  // from G4ProtonInelasticCrossSection
     216  G4double f1 = 8.0  - 8.0/A - 0.008*A;
     217  G4double f2 = 2.34 - 5.4/A - 0.0028*A;
     218
     219  res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
     220 
     221  f1 = 5.6 - 0.016*A;
     222  f2 = 1.37 + 1.37/A;
     223  res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
     224  return res; 
     225}
     226
     227//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     228
     229

  • trunk/source/processes/hadronic/cross_sections/src/G4CrossSectionDataStore.cc

    r819 r962  
    2424// ********************************************************************
    2525//
     26// $Id: G4CrossSectionDataStore.cc,v 1.16 2009/01/24 11:54:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
     28//
     29// -------------------------------------------------------------------
     30//
     31// GEANT4 Class file
     32//
     33//
     34// File name:     G4CrossSectionDataStore
     35//
     36//
     37// Modifications:
     38// 23.01.2009 V.Ivanchenko add destruction of data sets
     39//
    2640//
    2741
    … …  
    3145#include "G4HadTmpUtil.hh"
    3246#include "Randomize.hh"
    33 
     47#include "G4Nucleus.hh"
     48
     49G4CrossSectionDataStore::G4CrossSectionDataStore() :
     50  NDataSetList(0), verboseLevel(0)
     51{}
     52
     53G4CrossSectionDataStore::~G4CrossSectionDataStore()
     54{}
    3455
    3556G4double
    … …  
    140161
    141162
    142 std::pair<G4double, G4double>
    143 G4CrossSectionDataStore::SelectRandomIsotope(const G4DynamicParticle* particle,
    144                                              const G4Material* aMaterial)
     163G4Element* G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* particle,
     164                                                const G4Material* aMaterial,
     165                                                G4Nucleus& target)
     166{
     167  G4double aTemp = aMaterial->GetTemperature();
     168  const G4int nElements = aMaterial->GetNumberOfElements();
     169  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
     170  G4Element* anElement = (*theElementVector)[0];
     171  G4VCrossSectionDataSet* inCharge;
     172  G4int i;
     173
     174  // compounds
     175  if(1 < nElements) {
     176    G4double* xsec = new G4double [nElements];
     177    const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume();
     178    G4double cross = 0.0;
     179    for(i=0; i<nElements; i++) {
     180      anElement= (*theElementVector)[i];
     181      inCharge = whichDataSetInCharge(particle, anElement);
     182      cross   += theAtomsPerVolumeVector[i]*
     183        inCharge->GetCrossSection(particle, anElement, aTemp);
     184      xsec[i]  = cross;
     185    }
     186    cross *= G4UniformRand();
     187
     188    for(i=0; i<nElements; i++) {
     189      if( cross <=  xsec[i] ) {
     190        anElement = (*theElementVector)[i];
     191        break;
     192      }
     193    }
     194    delete [] xsec;
     195  }
     196
     197  // element have been selected
     198  inCharge = whichDataSetInCharge(particle, anElement);
     199  G4double ZZ = anElement->GetZ();
     200  G4double AA;
     201
     202  // Collect abundance weighted cross sections and A values for each isotope
     203  // in each element
     204 
     205  const G4int nIsoPerElement = anElement->GetNumberOfIsotopes();
     206
     207  // user-defined isotope abundances
     208  if (0 < nIsoPerElement) {
     209    G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
     210    AA = G4double((*isoVector)[0]->GetN());
     211    if(1 < nIsoPerElement) {
     212
     213      G4double* xsec  = new G4double [nIsoPerElement];
     214      G4double iso_xs = 0.0;
     215      G4double cross  = 0.0;
     216
     217      G4double* abundVector = anElement->GetRelativeAbundanceVector();
     218      G4bool elementXS = false;
     219      for (i = 0; i<nIsoPerElement; i++) {
     220        if (inCharge->IsZAApplicable(particle, ZZ, G4double((*isoVector)[i]->GetN()))) {
     221          iso_xs = inCharge->GetIsoCrossSection(particle, (*isoVector)[i], aTemp);
     222        } else if (elementXS == false) {
     223          iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
     224          elementXS = true;
     225        }
     226
     227        cross  += abundVector[i]*iso_xs;
     228        xsec[i] = cross;
     229      }
     230      cross *= G4UniformRand();
     231      for (i = 0; i<nIsoPerElement; i++) {
     232        if(cross <= xsec[i]) {
     233          AA = G4double((*isoVector)[i]->GetN());
     234          break;
     235        }
     236      }
     237      delete [] xsec;
     238    }
     239    // natural abundances
     240  } else {
     241
     242    G4StableIsotopes theDefaultIsotopes; 
     243    G4int Z = G4int(ZZ + 0.5);
     244    const G4int nIso = theDefaultIsotopes.GetNumberOfIsotopes(Z);
     245    G4int index = theDefaultIsotopes.GetFirstIsotope(Z);
     246    AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index));
     247   
     248    if(1 < nIso) {
     249
     250      G4double* xsec  = new G4double [nIso];
     251      G4double iso_xs = 0.0;
     252      G4double cross  = 0.0;
     253      G4bool elementXS= false;
     254
     255      for (i = 0; i<nIso; i++) {
     256        AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+i));
     257        if (inCharge->IsZAApplicable(particle, ZZ, AA )) {
     258          iso_xs = inCharge->GetIsoZACrossSection(particle, ZZ, AA, aTemp);
     259        } else if (elementXS == false) {
     260          iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
     261          elementXS = true;
     262        }
     263        cross  += theDefaultIsotopes.GetAbundance(index+i)*iso_xs;
     264        xsec[i] = cross;
     265      }
     266      cross *= G4UniformRand();
     267      for (i = 0; i<nIso; i++) {
     268        if(cross <= xsec[i]) {
     269          AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+i));
     270          break;
     271        }
     272      }
     273      delete [] xsec;
     274    }
     275  }
     276  //G4cout << "XS: " << particle->GetDefinition()->GetParticleName()
     277  //     << " e(GeV)= " << particle->GetKineticEnergy()/GeV
     278  //     << " in " << aMaterial->GetName()
     279  //     << " ZZ= " << ZZ << " AA= " << AA << "  " << anElement->GetName() << G4endl;
     280
     281  target.SetParameters(AA, ZZ);
     282  return anElement;
     283}
     284
     285/*
     286G4Element* G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* particle,
     287                                                const G4Material* aMaterial,
     288                                                G4Nucleus& target)
    145289{
    146290  static G4StableIsotopes theDefaultIsotopes;  // natural abundances and
    … …  
    229373  // cross section is zero
    230374   
     375  anElement = (*theElementVector)[0];
     376  G4double ZZ = anElement->GetZ();
     377  G4double AA = AvaluesPerElement[0][0];
     378  if (crossSectionTotal != 0.) {
     379    G4double random = G4UniformRand();
     380    for(G4int i=0; i < nElements; i++) {
     381      if(i!=0) runningSum[i] += runningSum[i-1];
     382      if(random <= runningSum[i]/crossSectionTotal) {
     383        anElement = (*theElementVector)[i];
     384        ZZ = anElement->GetZ();
     385
     386        // Compare random number to running sum over isotope xc to choose A
     387
     388        G4int nIso = awicsPerElement[i].size();
     389        G4double* running = new G4double[nIso];
     390        for (G4int j=0; j < nIso; j++) {
     391          running[j] = awicsPerElement[i][j];
     392          if(j!=0) running[j] += running[j-1];
     393        }
     394
     395        G4double trial = G4UniformRand();
     396        for (G4int j=0; j < nIso; j++) {
     397          AA = AvaluesPerElement[i][j];
     398          if (trial <= running[j]/running[nIso-1]) break;
     399        }
     400        delete [] running;
     401        break;
     402      }
     403    }
     404  }
     405
     406  //G4cout << "XS: " << particle->GetDefinition()->GetParticleName()
     407  //     << " e(GeV)= " << particle->GetKineticEnergy()/GeV
     408  //     << " in " << aMaterial->GetName()
     409  //     << " ZZ= " << ZZ << " AA= " << AA << "  " << anElement->GetName() << G4endl;
     410
     411  target.SetParameters(AA, ZZ);
     412  return anElement;
     413}
     414
     415std::pair<G4double, G4double>
     416G4CrossSectionDataStore::SelectRandomIsotope(const G4DynamicParticle* particle,
     417                                             const G4Material* aMaterial)
     418{
     419  static G4StableIsotopes theDefaultIsotopes;  // natural abundances and
     420                                               // stable isotopes
     421  G4double aTemp = aMaterial->GetTemperature();
     422  G4int nElements = aMaterial->GetNumberOfElements();
     423  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
     424  const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume();
     425  std::vector<std::vector<G4double> > awicsPerElement;
     426  std::vector<std::vector<G4double> > AvaluesPerElement;
     427  G4Element* anElement;
     428
     429  // Collect abundance weighted cross sections and A values for each isotope
     430  // in each element
     431 
     432  for (G4int i = 0; i < nElements; i++) {
     433    anElement = (*theElementVector)[i];
     434    G4int nIsoPerElement = anElement->GetNumberOfIsotopes();
     435    std::vector<G4double> isoholder;
     436    std::vector<G4double> aholder;
     437    G4double iso_xs = DBL_MIN;
     438
     439    if (nIsoPerElement) { // user-defined isotope abundances
     440      G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
     441      G4double* abundVector = anElement->GetRelativeAbundanceVector();
     442      G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement);
     443      G4bool elementXS = false;
     444      for (G4int j = 0; j < nIsoPerElement; j++) {
     445        if (inCharge->IsZAApplicable(particle, (*isoVector)[j]->GetZ(),
     446                                               (*isoVector)[j]->GetN() ) ) {
     447          iso_xs = inCharge->GetIsoCrossSection(particle, (*isoVector)[j], aTemp);
     448        } else if (elementXS == false) {
     449          iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
     450          elementXS = true;
     451        }
     452
     453        isoholder.push_back(abundVector[j]*iso_xs);
     454        aholder.push_back(G4double((*isoVector)[j]->GetN()));
     455      }
     456
     457    } else { // natural abundances
     458      G4int ZZ = G4lrint(anElement->GetZ());
     459      nIsoPerElement = theDefaultIsotopes.GetNumberOfIsotopes(ZZ);
     460      G4int index = theDefaultIsotopes.GetFirstIsotope(ZZ);
     461      G4double AA;
     462      G4double abundance;
     463
     464      G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement);
     465      G4bool elementXS = false;
     466
     467      for (G4int j = 0; j < nIsoPerElement; j++) {
     468        AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+j));
     469        aholder.push_back(AA);
     470        if (inCharge->IsZAApplicable(particle, G4double(ZZ), AA) ) {
     471          iso_xs = inCharge->GetIsoZACrossSection(particle, G4double(ZZ), AA, aTemp);
     472        } else if (elementXS == false) {
     473          iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
     474          elementXS = true;
     475        }
     476        abundance = theDefaultIsotopes.GetAbundance(index+j)/100.0;
     477        isoholder.push_back(abundance*iso_xs);
     478      }
     479    }
     480
     481    awicsPerElement.push_back(isoholder);
     482    AvaluesPerElement.push_back(aholder);
     483  }
     484
     485  // Calculate running sums for isotope selection
     486
     487  G4double crossSectionTotal = 0;
     488  G4double xSectionPerElement;
     489  std::vector<G4double> runningSum;
     490
     491  for (G4int i=0; i < nElements; i++) {
     492    xSectionPerElement = 0;
     493    for (G4int j=0; j < G4int(awicsPerElement[i].size()); j++)
     494                     xSectionPerElement += awicsPerElement[i][j];
     495    runningSum.push_back(theAtomsPerVolumeVector[i]*xSectionPerElement);
     496    crossSectionTotal += runningSum[i];
     497  }
     498
     499  // Compare random number to running sum over element xc to choose Z
     500
     501  // Initialize Z and A to first element and first isotope in case
     502  // cross section is zero
     503   
    231504  G4double ZZ = (*theElementVector)[0]->GetZ();
    232505  G4double AA = AvaluesPerElement[0][0];
    … …  
    259532  return std::pair<G4double, G4double>(ZZ, AA);
    260533}
    261 
     534*/
    262535
    263536void
    264537G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* aDataSet)
    265538{
    266    if (NDataSetList == NDataSetMax) {
    267       G4cout << "WARNING: G4CrossSectionDataStore::AddDataSet: "<<G4endl;
    268       G4cout << "         reached maximum number of data sets";
    269       G4cout << "         data set not added !!!!!!!!!!!!!!!!";
    270       return;
    271    }
    272    DataSetList[NDataSetList] = aDataSet;
    273    NDataSetList++;
    274 }
    275 
     539  DataSetList.push_back(aDataSet);
     540  NDataSetList++;
     541}
    276542
    277543void
    … …  
    279545BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
    280546{
    281    if (NDataSetList == 0)
    282    {
    283      G4Exception("G4CrossSectionDataStore", "007", FatalException,
    284                  "BuildPhysicsTable: no data sets registered");
    285      return;
    286    }
    287    for (G4int i = NDataSetList-1; i >= 0; i--) {
     547  if(NDataSetList > 0) {
     548    for (G4int i=0; i<NDataSetList; i++) {
    288549      DataSetList[i]->BuildPhysicsTable(aParticleType);
    289    }
     550    }
     551  }
    290552}
    291553
    … …  
    295557DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
    296558{
    297    if (NDataSetList == 0) {
    298       G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: no data sets registered"<<G4endl;
    299       return;
    300    }
    301    for (G4int i = NDataSetList-1; i >= 0; i--) {
    302       DataSetList[i]->DumpPhysicsTable(aParticleType);
    303    }
    304 }
     559  if (NDataSetList == 0) {
     560    G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: "
     561           << " no data sets registered"<<G4endl;
     562    return;
     563  }
     564  for (G4int i = NDataSetList-1; i >= 0; i--) {
     565    DataSetList[i]->DumpPhysicsTable(aParticleType);
     566  }
     567}

  • trunk/source/processes/hadronic/cross_sections/src/G4ElectroNuclearCrossSection.cc

    r819 r962  
    2525//
    2626//
    27 // $Id: G4ElectroNuclearCrossSection.cc,v 1.28 2008/01/17 10:07:11 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4ElectroNuclearCrossSection.cc,v 1.29 2008/10/24 19:15:17 dennis Exp $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030//
    … …  
    250250  //G4double mT= G4QPDGCode(111).GetNuclMass(Z,N,0);
    251251  G4double mT= 0.;
    252   if(G4NucleiPropertiesTable::IsInTable(Z,A)) mT=G4NucleiProperties::GetNuclearMass(A,Z);
     252  if(G4NucleiProperties::IsInStableTable(A,Z)) mT = G4NucleiProperties::GetNuclearMass(A,Z);
    253253  // If it is not in the Table of Stable Nuclei, then the Threshold=inf
    254254  else return infEn;             
    … …  
    256256  G4double mP= infEn;
    257257  //if(Z) mP= G4QPDGCode(111).GetNuclMass(Z-1,N,0);
    258   if(Z&&G4NucleiPropertiesTable::IsInTable(Z-1,A-1)) mP=G4NucleiProperties::GetNuclearMass(A-1,Z-1);
     258  if(Z&&G4NucleiProperties::IsInStableTable(A-1,Z-1)) mP = G4NucleiProperties::GetNuclearMass(A-1,Z-1);
    259259  else return infEn;
    260260  G4double mN= infEn;
    261261  //if(N) mN= G4QPDGCode(111).GetNuclMass(Z,N-1,0);
    262   if(N&&G4NucleiPropertiesTable::IsInTable(Z,A-1)) mN=G4NucleiProperties::GetNuclearMass(A-1,Z);
     262  if(N&&G4NucleiProperties::IsInStableTable(A-1,Z)) mN = G4NucleiProperties::GetNuclearMass(A-1,Z);
    263263  else return infEn;
    264264  G4double dP= mP+mProt-mT;

  • trunk/source/processes/hadronic/cross_sections/src/G4GlauberGribovCrossSection.cc

    r819 r962  
    4040//
    4141
     42const G4double G4GlauberGribovCrossSection::fNeutronBarCorrectionTot[93] = {
     43
     441.0, 1.0,     1.118517e+00, 1.082002e+00, 1.116171e+00, 1.078747e+00, 1.061315e+00,
     451.058205e+00, 1.082663e+00, 1.068500e+00, 1.076912e+00, 1.083475e+00, 1.079117e+00,
     461.071856e+00, 1.071990e+00, 1.073774e+00, 1.079356e+00, 1.081314e+00, 1.082056e+00,
     471.090772e+00, 1.096776e+00, 1.095828e+00, 1.097678e+00, 1.099157e+00, 1.103677e+00,
     481.105132e+00, 1.109806e+00, 1.110816e+00, 1.117378e+00, 1.115165e+00, 1.115710e+00,
     491.111855e+00, 1.110482e+00, 1.110112e+00, 1.106676e+00, 1.108706e+00, 1.105549e+00,
     501.106318e+00, 1.106242e+00, 1.107672e+00, 1.107342e+00, 1.108119e+00, 1.106655e+00,
     511.102588e+00, 1.096657e+00, 1.092920e+00, 1.086629e+00, 1.083592e+00, 1.076030e+00,
     521.083777e+00, 1.089460e+00, 1.086545e+00, 1.079924e+00, 1.082218e+00, 1.077798e+00,
     531.077062e+00, 1.072825e+00, 1.072241e+00, 1.072104e+00, 1.072490e+00, 1.069829e+00,
     541.070398e+00, 1.065458e+00, 1.064968e+00, 1.060524e+00, 1.060048e+00, 1.057620e+00,
     551.056428e+00, 1.055366e+00, 1.055017e+00, 1.052304e+00, 1.051767e+00, 1.049728e+00,
     561.048745e+00, 1.047399e+00, 1.045876e+00, 1.042972e+00, 1.041824e+00, 1.039993e+00,
     571.039021e+00, 1.036627e+00, 1.034176e+00, 1.032526e+00, 1.033633e+00, 1.036107e+00,
     581.037803e+00, 1.031266e+00, 1.032991e+00, 1.033284e+00, 1.035015e+00, 1.033945e+00,
     591.037075e+00, 1.034721e+00
     60
     61};
     62
     63const G4double G4GlauberGribovCrossSection::fNeutronBarCorrectionIn[93] = {
     64
     651.0, 1.0,     1.167421e+00, 1.156250e+00, 1.205364e+00, 1.154225e+00, 1.120391e+00,
     661.124632e+00, 1.129460e+00, 1.107863e+00, 1.102152e+00, 1.104593e+00, 1.100285e+00,
     671.098450e+00, 1.092677e+00, 1.101124e+00, 1.106461e+00, 1.115049e+00, 1.123903e+00,
     681.126661e+00, 1.131259e+00, 1.133949e+00, 1.134185e+00, 1.133767e+00, 1.132813e+00,
     691.131515e+00, 1.130338e+00, 1.134171e+00, 1.139206e+00, 1.141474e+00, 1.142189e+00,
     701.140725e+00, 1.140100e+00, 1.139848e+00, 1.137674e+00, 1.138645e+00, 1.136339e+00,
     711.136439e+00, 1.135946e+00, 1.136431e+00, 1.135702e+00, 1.135703e+00, 1.134113e+00,
     721.131935e+00, 1.128381e+00, 1.126373e+00, 1.122453e+00, 1.120908e+00, 1.115953e+00,
     731.115947e+00, 1.114426e+00, 1.111749e+00, 1.106207e+00, 1.107494e+00, 1.103622e+00,
     741.102576e+00, 1.098816e+00, 1.097889e+00, 1.097306e+00, 1.097130e+00, 1.094578e+00,
     751.094552e+00, 1.090222e+00, 1.089358e+00, 1.085409e+00, 1.084560e+00, 1.082182e+00,
     761.080773e+00, 1.079464e+00, 1.078724e+00, 1.076121e+00, 1.075235e+00, 1.073159e+00,
     771.071920e+00, 1.070395e+00, 1.069503e+00, 1.067525e+00, 1.066919e+00, 1.065779e+00,
     781.065319e+00, 1.063730e+00, 1.062092e+00, 1.061085e+00, 1.059908e+00, 1.059815e+00,
     791.059109e+00, 1.051920e+00, 1.051258e+00, 1.049473e+00, 1.048823e+00, 1.045984e+00,
     801.046435e+00, 1.042614e+00
     81
     82};
     83
     84const G4double G4GlauberGribovCrossSection::fProtonBarCorrectionTot[93] = {
     85
     861.0, 1.0,     
     871.118515e+00, 1.082000e+00, 1.116169e+00, 1.078745e+00, 1.061313e+00, 1.058203e+00,
     881.082661e+00, 1.068498e+00, 1.076910e+00, 1.083474e+00, 1.079115e+00, 1.071854e+00,
     891.071988e+00, 1.073772e+00, 1.079355e+00, 1.081312e+00, 1.082054e+00, 1.090770e+00,
     901.096774e+00, 1.095827e+00, 1.097677e+00, 1.099156e+00, 1.103676e+00, 1.105130e+00,
     911.109805e+00, 1.110814e+00, 1.117377e+00, 1.115163e+00, 1.115708e+00, 1.111853e+00,
     921.110480e+00, 1.110111e+00, 1.106674e+00, 1.108705e+00, 1.105548e+00, 1.106317e+00,
     931.106241e+00, 1.107671e+00, 1.107341e+00, 1.108118e+00, 1.106654e+00, 1.102586e+00,
     941.096655e+00, 1.092918e+00, 1.086628e+00, 1.083590e+00, 1.076028e+00, 1.083776e+00,
     951.089458e+00, 1.086543e+00, 1.079923e+00, 1.082216e+00, 1.077797e+00, 1.077061e+00,
     961.072824e+00, 1.072239e+00, 1.072103e+00, 1.072488e+00, 1.069828e+00, 1.070396e+00,
     971.065456e+00, 1.064966e+00, 1.060523e+00, 1.060047e+00, 1.057618e+00, 1.056427e+00,
     981.055365e+00, 1.055016e+00, 1.052303e+00, 1.051766e+00, 1.049727e+00, 1.048743e+00,
     991.047397e+00, 1.045875e+00, 1.042971e+00, 1.041823e+00, 1.039992e+00, 1.039019e+00,
     1001.036626e+00, 1.034175e+00, 1.032525e+00, 1.033632e+00, 1.036106e+00, 1.037802e+00,
     1011.031265e+00, 1.032990e+00, 1.033283e+00, 1.035014e+00, 1.033944e+00, 1.037074e+00,
     1021.034720e+00
     103
     104};
     105
     106const G4double G4GlauberGribovCrossSection::fProtonBarCorrectionIn[93] = {
     107
     1081.0, 1.0,     
     1091.167419e+00, 1.156248e+00, 1.205362e+00, 1.154224e+00, 1.120390e+00, 1.124630e+00,
     1101.129459e+00, 1.107861e+00, 1.102151e+00, 1.104591e+00, 1.100284e+00, 1.098449e+00,
     1111.092675e+00, 1.101122e+00, 1.106460e+00, 1.115048e+00, 1.123902e+00, 1.126659e+00,
     1121.131258e+00, 1.133948e+00, 1.134183e+00, 1.133766e+00, 1.132812e+00, 1.131514e+00,
     1131.130337e+00, 1.134170e+00, 1.139205e+00, 1.141472e+00, 1.142188e+00, 1.140724e+00,
     1141.140099e+00, 1.139847e+00, 1.137672e+00, 1.138644e+00, 1.136338e+00, 1.136438e+00,
     1151.135945e+00, 1.136429e+00, 1.135701e+00, 1.135702e+00, 1.134112e+00, 1.131934e+00,
     1161.128380e+00, 1.126371e+00, 1.122452e+00, 1.120907e+00, 1.115952e+00, 1.115946e+00,
     1171.114425e+00, 1.111748e+00, 1.106205e+00, 1.107493e+00, 1.103621e+00, 1.102575e+00,
     1181.098815e+00, 1.097888e+00, 1.097305e+00, 1.097129e+00, 1.094577e+00, 1.094551e+00,
     1191.090221e+00, 1.089357e+00, 1.085408e+00, 1.084559e+00, 1.082181e+00, 1.080772e+00,
     1201.079463e+00, 1.078723e+00, 1.076120e+00, 1.075234e+00, 1.073158e+00, 1.071919e+00,
     1211.070394e+00, 1.069502e+00, 1.067524e+00, 1.066918e+00, 1.065778e+00, 1.065318e+00,
     1221.063729e+00, 1.062091e+00, 1.061084e+00, 1.059907e+00, 1.059814e+00, 1.059108e+00,
     1231.051919e+00, 1.051257e+00, 1.049472e+00, 1.048822e+00, 1.045983e+00, 1.046434e+00,
     1241.042613e+00
     125
     126};
     127
     128
     129const G4double G4GlauberGribovCrossSection::fPionPlusBarCorrectionTot[93] = {
     130
     1311.0, 1.0,     
     1321.075927e+00, 1.074407e+00, 1.126098e+00, 1.100127e+00, 1.089742e+00, 1.083536e+00,
     1331.089988e+00, 1.103566e+00, 1.096922e+00, 1.126573e+00, 1.132734e+00, 1.136512e+00,
     1341.136629e+00, 1.133086e+00, 1.132428e+00, 1.129299e+00, 1.125622e+00, 1.126992e+00,
     1351.127840e+00, 1.162670e+00, 1.160392e+00, 1.157864e+00, 1.157227e+00, 1.154627e+00,
     1361.192555e+00, 1.197243e+00, 1.197911e+00, 1.200326e+00, 1.220053e+00, 1.215019e+00,
     1371.211703e+00, 1.209080e+00, 1.204248e+00, 1.203328e+00, 1.198671e+00, 1.196840e+00,
     1381.194392e+00, 1.193037e+00, 1.190408e+00, 1.188583e+00, 1.206127e+00, 1.210028e+00,
     1391.206434e+00, 1.204456e+00, 1.200547e+00, 1.199058e+00, 1.200174e+00, 1.200276e+00,
     1401.198912e+00, 1.213048e+00, 1.207160e+00, 1.208020e+00, 1.203814e+00, 1.202380e+00,
     1411.198306e+00, 1.197002e+00, 1.196027e+00, 1.195449e+00, 1.192563e+00, 1.192135e+00,
     1421.187556e+00, 1.186308e+00, 1.182124e+00, 1.180900e+00, 1.178224e+00, 1.176471e+00,
     1431.174811e+00, 1.173702e+00, 1.170827e+00, 1.169581e+00, 1.167205e+00, 1.165626e+00,
     1441.180244e+00, 1.177626e+00, 1.175121e+00, 1.173903e+00, 1.172192e+00, 1.171128e+00,
     1451.168997e+00, 1.166826e+00, 1.164130e+00, 1.165412e+00, 1.165504e+00, 1.165020e+00,
     1461.158462e+00, 1.158014e+00, 1.156519e+00, 1.156081e+00, 1.153602e+00, 1.154190e+00,
     1471.152974e+00
     148 
     149};
     150
     151const G4double G4GlauberGribovCrossSection::fPionPlusBarCorrectionIn[93] = {
     152
     1531.0, 1.0,   
     1541.140246e+00, 1.097872e+00, 1.104301e+00, 1.068722e+00, 1.044495e+00, 1.062622e+00,
     1551.047987e+00, 1.037032e+00, 1.035686e+00, 1.042870e+00, 1.052222e+00, 1.065100e+00,
     1561.070480e+00, 1.078286e+00, 1.081488e+00, 1.089713e+00, 1.099105e+00, 1.098003e+00,
     1571.102175e+00, 1.117707e+00, 1.121734e+00, 1.125229e+00, 1.126457e+00, 1.128905e+00,
     1581.137312e+00, 1.126263e+00, 1.126459e+00, 1.115191e+00, 1.116986e+00, 1.117184e+00,
     1591.117037e+00, 1.116777e+00, 1.115858e+00, 1.115745e+00, 1.114489e+00, 1.113993e+00,
     1601.113226e+00, 1.112818e+00, 1.111890e+00, 1.111238e+00, 1.111209e+00, 1.111775e+00,
     1611.110256e+00, 1.109414e+00, 1.107647e+00, 1.106980e+00, 1.106096e+00, 1.107331e+00,
     1621.107849e+00, 1.106407e+00, 1.103426e+00, 1.103896e+00, 1.101756e+00, 1.101031e+00,
     1631.098915e+00, 1.098260e+00, 1.097768e+00, 1.097487e+00, 1.095964e+00, 1.095773e+00,
     1641.093348e+00, 1.092687e+00, 1.090465e+00, 1.089821e+00, 1.088394e+00, 1.087462e+00,
     1651.086571e+00, 1.085997e+00, 1.084451e+00, 1.083798e+00, 1.082513e+00, 1.081670e+00,
     1661.080735e+00, 1.075659e+00, 1.074341e+00, 1.073689e+00, 1.072787e+00, 1.072237e+00,
     1671.071107e+00, 1.069955e+00, 1.064856e+00, 1.065873e+00, 1.065938e+00, 1.065694e+00,
     1681.062192e+00, 1.061967e+00, 1.061180e+00, 1.060960e+00, 1.059646e+00, 1.059975e+00,
     1691.059658e+00
     170 
     171};
     172
     173
     174const G4double G4GlauberGribovCrossSection::fPionMinusBarCorrectionTot[93] = {
     175
     1761.0, 1.0,     
     1771.075927e+00, 1.077959e+00, 1.129145e+00, 1.102088e+00, 1.089765e+00, 1.083542e+00,
     1781.089995e+00, 1.104895e+00, 1.097154e+00, 1.127663e+00, 1.133063e+00, 1.137425e+00,
     1791.136724e+00, 1.133859e+00, 1.132498e+00, 1.130276e+00, 1.127896e+00, 1.127656e+00,
     1801.127905e+00, 1.164210e+00, 1.162259e+00, 1.160075e+00, 1.158978e+00, 1.156649e+00,
     1811.194157e+00, 1.199177e+00, 1.198983e+00, 1.202325e+00, 1.221967e+00, 1.217548e+00,
     1821.214389e+00, 1.211760e+00, 1.207335e+00, 1.206081e+00, 1.201766e+00, 1.199779e+00,
     1831.197283e+00, 1.195706e+00, 1.193071e+00, 1.191115e+00, 1.208838e+00, 1.212681e+00,
     1841.209235e+00, 1.207163e+00, 1.203451e+00, 1.201807e+00, 1.203283e+00, 1.203388e+00,
     1851.202244e+00, 1.216509e+00, 1.211066e+00, 1.211504e+00, 1.207539e+00, 1.205991e+00,
     1861.202143e+00, 1.200724e+00, 1.199595e+00, 1.198815e+00, 1.196025e+00, 1.195390e+00,
     1871.191137e+00, 1.189791e+00, 1.185888e+00, 1.184575e+00, 1.181996e+00, 1.180229e+00,
     1881.178545e+00, 1.177355e+00, 1.174616e+00, 1.173312e+00, 1.171016e+00, 1.169424e+00,
     1891.184120e+00, 1.181478e+00, 1.179085e+00, 1.177817e+00, 1.176124e+00, 1.175003e+00,
     1901.172947e+00, 1.170858e+00, 1.168170e+00, 1.169397e+00, 1.169304e+00, 1.168706e+00,
     1911.162774e+00, 1.162217e+00, 1.160740e+00, 1.160196e+00, 1.157857e+00, 1.158220e+00,
     1921.157267e+00
     193};
     194
     195
     196const G4double G4GlauberGribovCrossSection::fPionMinusBarCorrectionIn[93] = {
     197
     1981.0, 1.0,   
     1991.140246e+00, 1.100898e+00, 1.106773e+00, 1.070289e+00, 1.044514e+00, 1.062628e+00,
     2001.047992e+00, 1.038041e+00, 1.035862e+00, 1.043679e+00, 1.052466e+00, 1.065780e+00,
     2011.070551e+00, 1.078869e+00, 1.081541e+00, 1.090455e+00, 1.100847e+00, 1.098511e+00,
     2021.102226e+00, 1.118865e+00, 1.123143e+00, 1.126904e+00, 1.127785e+00, 1.130444e+00,
     2031.138502e+00, 1.127678e+00, 1.127244e+00, 1.116634e+00, 1.118347e+00, 1.118988e+00,
     2041.118957e+00, 1.118696e+00, 1.118074e+00, 1.117722e+00, 1.116717e+00, 1.116111e+00,
     2051.115311e+00, 1.114745e+00, 1.113814e+00, 1.113069e+00, 1.113141e+00, 1.113660e+00,
     2061.112249e+00, 1.111343e+00, 1.109718e+00, 1.108942e+00, 1.108310e+00, 1.109549e+00,
     2071.110227e+00, 1.108846e+00, 1.106183e+00, 1.106354e+00, 1.104388e+00, 1.103583e+00,
     2081.101632e+00, 1.100896e+00, 1.100296e+00, 1.099873e+00, 1.098420e+00, 1.098082e+00,
     2091.095892e+00, 1.095162e+00, 1.093144e+00, 1.092438e+00, 1.091083e+00, 1.090142e+00,
     2101.089236e+00, 1.088604e+00, 1.087159e+00, 1.086465e+00, 1.085239e+00, 1.084388e+00,
     2111.083473e+00, 1.078373e+00, 1.077136e+00, 1.076450e+00, 1.075561e+00, 1.074973e+00,
     2121.073898e+00, 1.072806e+00, 1.067706e+00, 1.068684e+00, 1.068618e+00, 1.068294e+00,
     2131.065241e+00, 1.064939e+00, 1.064166e+00, 1.063872e+00, 1.062659e+00, 1.062828e+00,
     2141.062699e+00
     215
     216};
     217
     218
     219
     220
     221////////////////////////////////////////////////////////////////////////////////
     222//
     223//
    42224
    43225G4GlauberGribovCrossSection::G4GlauberGribovCrossSection()
    44 : fUpperLimit( 10000 * GeV ),
     226: fUpperLimit( 100000 * GeV ),
    45227  fLowerLimit( 3 * GeV ),
    46228  fRadiusConst( 1.08*fermi )  // 1.1, 1.3 ?
    … …  
    121303         theParticle == theSMinus)      )    || 
    122304
    123        ( kineticEnergy  >= 0.1*fLowerLimit &&
     305       ( kineticEnergy  >= fLowerLimit &&
    124306         Z > 1.5 &&      // >=  He
    125307       ( theParticle == theProton    ||
    … …  
    185367  xsection =  nucleusSquare*std::log( 1. + ratio );
    186368
     369  xsection *= GetParticleBarCorTot(theParticle, Z);
     370
    187371  fTotalXsc = xsection;
    188372
    … …  
    190374
    191375  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
     376
     377  fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);
    192378
    193379  fElasticXsc   = fTotalXsc - fInelasticXsc;

  • trunk/source/processes/hadronic/cross_sections/src/G4HadronCaptureDataSet.cc

    r819 r962  
    2626//
    2727// $Id: G4HadronCaptureDataSet.cc,v 1.8 2006/06/29 19:57:35 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030//

  • trunk/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc

    r819 r962  
    2525//
    2626//
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929//

  • trunk/source/processes/hadronic/cross_sections/src/G4HadronElasticDataSet.cc

    r819 r962  
    2626//
    2727// $Id: G4HadronElasticDataSet.cc,v 1.8 2006/06/29 19:57:39 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030//

  • trunk/source/processes/hadronic/cross_sections/src/G4HadronFissionDataSet.cc

    r819 r962  
    2626//
    2727// $Id: G4HadronFissionDataSet.cc,v 1.8 2006/06/29 19:57:41 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030//

  • trunk/source/processes/hadronic/cross_sections/src/G4HadronInelasticDataSet.cc

    r819 r962  
    2626//
    2727// $Id: G4HadronInelasticDataSet.cc,v 1.8 2006/06/29 19:57:43 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030//

  • trunk/source/processes/hadronic/cross_sections/src/G4PhotoNuclearCrossSection.cc

    r819 r962  
    2626//
    2727// The lust update: M.V. Kossov, CERN/ITEP(Moscow) 17-June-02
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030//
    … …  
    284284  //G4double mT= G4QPDGCode(111).GetNuclMass(Z,N,0);
    285285  G4double mT= 0.;
    286   if(G4NucleiPropertiesTable::IsInTable(Z,A))
     286  if(G4NucleiProperties::IsInStableTable(A,Z))
    287287    mT=G4NucleiProperties::GetNuclearMass(A,Z);
    288288  else
    … …  
    297297  //if(Z) mP= G4QPDGCode(111).GetNuclMass(Z-1,N,0);
    298298
    299   if(Z && G4NucleiPropertiesTable::IsInTable(Z-1,A-1))
     299  if(Z && G4NucleiProperties::IsInStableTable(A-1,Z-1))
    300300  {
    301301    mP = G4NucleiProperties::GetNuclearMass(A-1,Z-1);
    … …  
    310310  G4double mN= infEn;
    311311  //if(N) mN= G4QPDGCode(111).GetNuclMass(Z,N-1,0);
    312   if(N&&G4NucleiPropertiesTable::IsInTable(Z,A-1))
     312  if(N&&G4NucleiProperties::IsInStableTable(A-1,Z))
    313313    mN=G4NucleiProperties::GetNuclearMass(A-1,Z);
    314314  /*

  • trunk/source/processes/hadronic/cross_sections/src/G4PiData.cc

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // GEANT4 tag $Name: $
     26// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2727//
    2828// --------------------------------------------------------------------

  • trunk/source/processes/hadronic/cross_sections/src/G4PiNuclearCrossSection.cc

    r819 r962  
    2727 #include "G4HadronicException.hh"
    2828 #include "G4HadTmpUtil.hh"
    29  #include "G4ping.hh"
     29// #include "G4ping.hh"
    3030
    3131 // by J.P Wellisch, Sun Sep 15 2002.
    … …  
    350350{
    351351  // precondition
    352   G4ping debug("debug_PiNuclearCrossSection");
     352  //  G4ping debug("debug_PiNuclearCrossSection");
    353353  using namespace std;
    354354  G4bool ok = false;
    … …  
    368368  G4double result = 0;
    369369  G4int Z=G4lrint(ZZ);
    370   debug.push_back(Z);
     370  //  debug.push_back(Z);
    371371  size_t it = 0;
    372372
    373373  while(it<theZ.size() && Z>theZ[it]) it++;
    374374
    375   debug.push_back(theZ[it]);
    376   debug.push_back(kineticEnergy);
     375  //  debug.push_back(theZ[it]);
     376  //  debug.push_back(kineticEnergy);
    377377
    378378  if(Z > theZ[it])
    … …  
    390390      fTotalXsc = thePimData[it]->TotalXSection(kineticEnergy);
    391391
    392       debug.push_back("D1 ");
    393       debug.push_back(result);
    394       debug.push_back(fTotalXsc);
     392      //      debug.push_back("D1 ");
     393      //      debug.push_back(result);
     394      //      debug.push_back(fTotalXsc);
    395395    }
    396396    else
    … …  
    406406      fTotalXsc = Interpolate(Z1, Z2, Z, xt1, xt2);
    407407
    408       debug.push_back("D2 ");
    409       debug.push_back(x1);
    410       debug.push_back(x2);
    411       debug.push_back(xt1);
    412       debug.push_back(xt2);
    413       debug.push_back(Z1);
    414       debug.push_back(Z2);
    415       debug.push_back(result);
    416       debug.push_back(fTotalXsc);
     408      //      debug.push_back("D2 ");
     409      //      debug.push_back(x1);
     410      //      debug.push_back(x2);
     411      //      debug.push_back(xt1);
     412      //      debug.push_back(xt2);
     413      //      debug.push_back(Z1);
     414      //      debug.push_back(Z2);
     415      //      debug.push_back(result);
     416      //      debug.push_back(fTotalXsc);
    417417    }
    418418  }
    … …  
    430430      fTotalXsc = theData->operator[](it)->TotalXSection(kineticEnergy);
    431431
    432       debug.push_back("D3 ");
    433       debug.push_back(result);
    434       debug.push_back(fTotalXsc);
     432      //      debug.push_back("D3 ");
     433      //      debug.push_back(result);
     434      //      debug.push_back(fTotalXsc);
    435435    }
    436436    else
    … …  
    456456      fTotalXsc = Interpolate(Z1, Z2, Z, xt1, xt2);
    457457
    458       debug.push_back("D4 ");
    459       debug.push_back(x1);
    460       debug.push_back(xt1);
    461       debug.push_back(x2);
    462       debug.push_back(xt2);
    463       debug.push_back(Z1);
    464       debug.push_back(Z2);
    465       debug.push_back(result);
    466       debug.push_back(fTotalXsc);
     458      //      debug.push_back("D4 ");
     459      //      debug.push_back(x1);
     460      //      debug.push_back(xt1);
     461      //      debug.push_back(x2);
     462      //      debug.push_back(xt2);
     463      //      debug.push_back(Z1);
     464      //      debug.push_back(Z2);
     465      //      debug.push_back(result);
     466      //      debug.push_back(fTotalXsc);
    467467    }
    468468  }
    469   debug.dump();
     469  //  debug.dump();
    470470
    471471  fElasticXsc = fTotalXsc - result;

  • trunk/source/processes/hadronic/cross_sections/src/G4UElasticCrossSection.cc

    r819 r962  
    5050//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    5151
    52 G4UElasticCrossSection::G4UElasticCrossSection(const G4ParticleDefinition* p)
    53 {
     52G4UElasticCrossSection::G4UElasticCrossSection(const G4ParticleDefinition*)
     53{
     54  verboseLevel = 0;
    5455  hasGlauber = false;
    55   thEnergy   = 100.*GeV;
     56  thEnergy   = 90.*GeV;
    5657  fGlauber   = new G4GlauberGribovCrossSection();
    5758  fGheisha   = G4HadronCrossSections::Instance();
    5859  fNucleon   = 0;
    5960  fUPi       = 0;
    60   if(p == G4Proton::Proton() || p == G4Neutron::Neutron())
    61     fNucleon = new G4NucleonNuclearCrossSection();
    62   else if(p == G4PionPlus::PionPlus() || p == G4PionMinus::PionMinus())
    63     fUPi = new G4UPiNuclearCrossSection();
    64   verboseLevel = 0;
    65   Initialise(p);
    6661}
    6762
    … …  
    117112  G4double cross = 0.0;
    118113  G4double ekin = dp->GetKineticEnergy();
    119   G4int iz = G4int(Z + 0.5);
     114  G4int iz = G4int(Z);
    120115  if(iz > 92) iz = 92;
    121116
    … …  
    160155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    161156
    162 void G4UElasticCrossSection::BuildPhysicsTable(const G4ParticleDefinition&)
    163 {
     157void G4UElasticCrossSection::BuildPhysicsTable(const G4ParticleDefinition& p)
     158{
     159  if(&p == G4Proton::Proton() || &p == G4Neutron::Neutron()) {
     160    fNucleon = new G4NucleonNuclearCrossSection();
     161  } else if(&p == G4PionPlus::PionPlus() || &p == G4PionMinus::PionMinus()) {
     162    fUPi = new G4UPiNuclearCrossSection();
     163  }
     164  Initialise(&p);
    164165}
    165166

  • trunk/source/processes/hadronic/cross_sections/src/G4UInelasticCrossSection.cc

    r819 r962  
    5050//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    5151
    52 G4UInelasticCrossSection::G4UInelasticCrossSection(const G4ParticleDefinition* p)
    53 {
     52G4UInelasticCrossSection::G4UInelasticCrossSection(const G4ParticleDefinition*)
     53{
     54  verboseLevel = 0;
    5455  hasGlauber = false;
    55   thEnergy   = 100.*GeV;
     56  thEnergy   = 90.*GeV;
    5657  fGlauber   = new G4GlauberGribovCrossSection();
    5758  fGheisha   = G4HadronCrossSections::Instance();
    5859  fNucleon   = 0;
    5960  fUPi       = 0;
    60   if(p == G4Proton::Proton() || p == G4Neutron::Neutron())
    61     fNucleon = new G4NucleonNuclearCrossSection();
    62   else if(p == G4PionPlus::PionPlus() || p == G4PionMinus::PionMinus())
    63     fUPi = new G4UPiNuclearCrossSection();
    64   verboseLevel = 0;
    65   Initialise(p);
    6661}
    6762
    … …  
    117112  G4double cross = 0.0;
    118113  G4double ekin = dp->GetKineticEnergy();
    119   G4int iz = G4int(Z + 0.5);
     114  G4int iz = G4int(Z);
    120115  if(iz > 92) iz = 92;
    121116
    … …  
    160155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    161156
    162 void G4UInelasticCrossSection::BuildPhysicsTable(const G4ParticleDefinition&)
    163 {
     157void G4UInelasticCrossSection::BuildPhysicsTable(const G4ParticleDefinition& p)
     158{
     159  if(&p == G4Proton::Proton() || &p == G4Neutron::Neutron()) {
     160    fNucleon = new G4NucleonNuclearCrossSection();
     161  } else if(&p == G4PionPlus::PionPlus() || &p == G4PionMinus::PionMinus()) {
     162    fUPi = new G4UPiNuclearCrossSection();
     163  }
     164  Initialise(&p);
    164165}
    165166

  • trunk/source/processes/hadronic/cross_sections/src/G4UPiNuclearCrossSection.cc

    r819 r962  
    9898         G4double Z, G4double A, G4double ekin, G4PhysicsTable* table)
    9999{
     100  G4double res = 0.0;
    100101  G4bool b;
    101102  G4int idx;
    102103  G4int iz = G4int(Z + 0.5);
    103   for(idx=0; idx<NZ-1; idx++) {if(theZ[idx] >= iz) break;}
    104 
    105   G4double x = (((*table)[idx])->GetValue(ekin, b));
    106 
    107   // one of elements in the table
    108   if(iz >= theZ[idx]) {
    109     G4double A1 = theA[idx];
    110     x *= std::pow(A/A1,aPower);
     104  if(iz > 92) iz = 92;
     105  for(idx=0; idx<NZ; idx++) {if(theZ[idx] >= iz) break;}
     106  if(idx >= NZ) idx = NZ - 1;
     107  G4int iz2 = theZ[idx];
     108
     109  G4double x2 = (((*table)[idx])->GetValue(ekin, b))*APower[iz]/APower[iz2];
     110
     111  // use only one Z
     112  if(iz >= theZ[idx] || idx == 0) {
     113    res = x2;
    111114
    112115  // Interpolation between Z
    113116  } else {
    114117
    115     G4double x2 = x;
    116     G4double x1 = x;
    117     if(idx == 0) {
    118       idx++;
    119       x2 = (((*table)[idx])->GetValue(ekin, b));
    120     } else {
    121       x1 = (((*table)[idx-1])->GetValue(ekin, b));
    122     }
    123     G4double A1 = theA[idx-1];
    124     G4double A2 = theA[idx];
    125     G4double w1 = A - A1;
    126     G4double w2 = A2 - A;
    127     G4double y1 = x1*std::pow(A/A1,aPower);
    128     if(w1 <= 0.0) x = y1;
    129     else {
    130       G4double y2 = x2*std::pow(A/A2,aPower);
    131       if(w2 <= 0.0) x = y2;
    132       else  x = (w2*y1 + w1*y2)/(w1 + w2);
    133     }
    134   }
    135   return x;
     118    G4int iz1 = theZ[idx-1];
     119    G4double x1 = (((*table)[idx-1])->GetValue(ekin, b))*APower[iz]/APower[iz1];
     120    G4double w1 = A - theA[idx-1];
     121    G4double w2 = theA[idx] - A;
     122    res = (w1*x1 + w2*x2)/(w1 + w2);
     123  }
     124  return res;
    136125}
    137126
    … …  
    143132{
    144133  G4LPhysicsFreeVector* pvin = new G4LPhysicsFreeVector(n,e[0]*GeV,e[n-1]*GeV);
     134  //pvin->SetSpline(true);
    145135  G4LPhysicsFreeVector* pvel = new G4LPhysicsFreeVector(n,e[0]*GeV,e[n-1]*GeV);
     136  //pvel->SetSpline(true);
    146137  for(G4int i=0; i<n; i++) {
    147138    pvin->PutValues(i,e[i]*GeV,in[i]*millibarn);
    … …  
    188179
    189180  G4NistManager* nist = G4NistManager::Instance();
    190   for(G4int i=0; i<n; i++) {
     181  G4int i;
     182  for(i=0; i<n; i++) {
    191183    theZ.push_back(iz[i]);
    192184    theA.push_back(nist->GetAtomicMassAmu(iz[i]));
     185  }
     186  for(i=1; i<92; i++) {
     187    APower[i] = std::pow(nist->GetAtomicMassAmu(i),aPower);
    193188  }
    194189

  • trunk/source/processes/hadronic/cross_sections/src/G4VCrossSectionDataSet.cc

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VCrossSectionDataSet.cc,v 1.7 2006/12/29 02:05:48 dennis Exp $
    27 //
     26// $Id: G4VCrossSectionDataSet.cc,v 1.8 2009/01/24 11:54:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
     28//
     29// -------------------------------------------------------------------
     30//
     31// GEANT4 Class file
     32//
     33//
     34// File name:    G4VCrossSectionDataSet
     35//
     36// Author  F.W. Jones, TRIUMF, 20-JAN-97
     37//
     38// Modifications:
     39// 23.01.2009 V.Ivanchenko move constructor and destructor to source
     40//
    2841
    2942#include "G4VCrossSectionDataSet.hh"
     43#include "G4CrossSectionDataSetRegistry.hh"
     44
     45G4VCrossSectionDataSet::G4VCrossSectionDataSet() :
     46  verboseLevel(0)
     47{
     48  G4CrossSectionDataSetRegistry::Instance()->Register(this);
     49}
     50
     51G4VCrossSectionDataSet::~G4VCrossSectionDataSet()
     52{
     53  G4CrossSectionDataSetRegistry::Instance()->DeRegister(this);
     54}
    3055
    3156// Override this method to test for particle and isotope applicability
Note: See TracChangeset for help on using the changeset viewer.