Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (15 years ago)
Author:
garnier
Message:

update ti head

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

Legend:

Unmodified
Added
Removed

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BGGNucleonElasticXS.cc,v 1.7 2009/11/19 19:40:45 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BGGNucleonElasticXS.cc,v 1.8 2010/10/12 06:02:41 dennis Exp $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5151#include "G4NistManager.hh"
    5252
    53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    5453
    5554G4BGGNucleonElasticXS::G4BGGNucleonElasticXS(const G4ParticleDefinition*)
    … …  
    7574}
    7675
    77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    78 
    79 G4double G4BGGNucleonElasticXS::GetIsoZACrossSection(const G4DynamicParticle* dp,
    80                                                        G4double Z,
    81                                                        G4double A,
    82                                                        G4double)
     76
     77G4double
     78G4BGGNucleonElasticXS::GetZandACrossSection(const G4DynamicParticle* dp,
     79                                            G4int Z, G4int A, G4double)
    8380{
    8481  G4double cross = 0.0;
    8582  G4double ekin = dp->GetKineticEnergy();
    86   G4int iz = G4int(Z);
    87   if(iz > 92) iz = 92;
     83  if(Z > 92) Z = 92;
    8884
    8985  if(ekin <= fLowEnergy) {
    90     cross = theCoulombFac[iz];
     86    cross = theCoulombFac[Z];
    9187    if(isProton) { cross *= CoulombFactor(ekin, A); }
    9288
    93   } else if(iz == 1) {
    94     if( A < 1.5) {
     89  } else if(Z == 1) {
     90    if( A < 2) {
    9591      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
    9692      cross = fHadron->GetElasticHadronNucleonXsc();
    … …  
    10298    }
    10399  } else if(ekin > fGlauberEnergy) {
    104     cross = theGlauberFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
     100    cross = theGlauberFac[Z]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
    105101  } else {
    106102    cross = fNucleon->GetElasticCrossSection(dp, Z, A);
    … …  
    159155  G4NistManager* nist = G4NistManager::Instance();
    160156
    161   G4double A, csup, csdn;
     157  G4double csup, csdn;
     158  G4int A;
    162159
    163160  if(verboseLevel > 0) G4cout << "### G4BGGNucleonElasticXS::Initialise for "
    … …  
    167164
    168165    G4double Z = G4double(iz);
    169     A = nist->GetAtomicMassAmu(iz);
    170 
    171     csup = fGlauber->GetElasticGlauberGribov(&dp, Z, A);
    172     csdn = fNucleon->GetElasticCrossSection(&dp, Z, A);
     166    A = G4lrint(nist->GetAtomicMassAmu(iz));
     167
     168    csup = fGlauber->GetElasticGlauberGribov(&dp, iz, A);
     169    csdn = fNucleon->GetElasticCrossSection(&dp, iz, A);
    173170
    174171    theGlauberFac[iz] = csdn/csup;
    … …  
    179176  fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
    180177  theCoulombFac[1] = fHadron->GetElasticHadronNucleonXsc();
    181   if(isProton) { theCoulombFac[1] /= CoulombFactor(fLowEnergy, 1.0); }
     178  if(isProton) { theCoulombFac[1] /= CoulombFactor(fLowEnergy, 1); }
    182179
    183180  for(G4int iz=2; iz<93; iz++) {
    184181
    185182    G4double Z = G4double(iz);
    186     A = nist->GetAtomicMassAmu(iz);
    187 
    188     theCoulombFac[iz] = fNucleon->GetElasticCrossSection(&dp, Z, A);
     183    A = G4lrint(nist->GetAtomicMassAmu(iz));
     184
     185    theCoulombFac[iz] = fNucleon->GetElasticCrossSection(&dp, iz, A);
    189186    if(isProton) { theCoulombFac[iz] /= CoulombFactor(fLowEnergy, A); }
    190187    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
    … …  
    195192//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    196193
    197 G4double G4BGGNucleonElasticXS::CoulombFactor(G4double kinEnergy, G4double A)
    198 {
    199   G4double res= 0.0;
     194G4double G4BGGNucleonElasticXS::CoulombFactor(G4double kinEnergy, G4int A)
     195{
     196  G4double res = 0.0;
    200197  if(kinEnergy <= DBL_MIN) return res;
    201   else if(A < 1.5) return kinEnergy*kinEnergy;
     198  else if(A < 2) return kinEnergy*kinEnergy;
    202199
    203200  G4double elog = std::log10(kinEnergy/GeV);
     201  G4double aa = A;
    204202
    205203  // from G4ProtonInelasticCrossSection
    206   G4double f1 = 8.0  - 8.0/A - 0.008*A;
    207   G4double f2 = 2.34 - 5.4/A - 0.0028*A;
     204  G4double f1 = 8.0  - 8.0/aa - 0.008*aa;
     205  G4double f2 = 2.34 - 5.4/aa - 0.0028*aa;
    208206
    209207  res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
    210208 
    211   f1 = 5.6 - 0.016*A;
    212   f2 = 1.37 + 1.37/A;
    213   res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
     209  f1 = 5.6 - 0.016*aa;
     210  f2 = 1.37 + 1.37/aa;
     211  res *= ( 1.0 + (0.8 + 18./aa - 0.002*aa)/(1.0 + std::exp(f1*(elog + f2))));
    214212  return res; 
    215213}
    216214
    217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    218 
    219 
     215

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BGGNucleonInelasticXS.cc,v 1.7 2009/11/19 19:40:45 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BGGNucleonInelasticXS.cc,v 1.9 2010/10/12 06:16:19 dennis Exp $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7575}
    7676
    77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    78 
    79 G4double G4BGGNucleonInelasticXS::GetIsoZACrossSection(const G4DynamicParticle* dp,
    80                                                        G4double Z,
    81                                                        G4double A,
    82                                                        G4double)
     77
     78G4double
     79G4BGGNucleonInelasticXS::GetZandACrossSection(const G4DynamicParticle* dp,
     80                                              G4int Z, G4int A, G4double)
    8381{
    8482  G4double cross = 0.0;
    8583  G4double ekin = dp->GetKineticEnergy();
    86   G4int iz = G4int(Z);
    87   if(iz > 92) iz = 92;
     84//  G4int iz = G4int(Z);
     85  if(Z > 92) Z = 92;
    8886
    8987  if(ekin <= fLowEnergy) {
    90     cross = theCoulombFac[iz]*CoulombFactor(ekin, A);
    91   } else if(iz == 1) {
    92     if( A < 1.5) {
     88    cross = theCoulombFac[Z]*CoulombFactor(ekin, A);
     89  } else if(Z == 1) {
     90    if( A < 2) {
    9391      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
    9492      cross = fHadron->GetInelasticHadronNucleonXsc();
    … …  
    10098    }
    10199  } else if(ekin > fGlauberEnergy) {
    102     cross = theGlauberFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
     100    cross = theGlauberFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
    103101  } else {
    104     cross = fNucleon->GetIsoZACrossSection(dp, Z, A);
     102    cross = fNucleon->GetZandACrossSection(dp, Z, A);
    105103  }
    106104
    … …  
    156154
    157155  G4NistManager* nist = G4NistManager::Instance();
    158   G4double A = nist->GetAtomicMassAmu(2);
     156  G4int A = G4lrint(nist->GetAtomicMassAmu(2));
    159157
    160158  G4double csup, csdn;
    … …  
    166164
    167165    G4double Z = G4double(iz);
    168     A = nist->GetAtomicMassAmu(iz);
    169 
    170     csup = fGlauber->GetInelasticGlauberGribov(&dp, Z, A);
    171     csdn = fNucleon->GetIsoZACrossSection(&dp, Z, A);
     166    A = G4lrint(nist->GetAtomicMassAmu(iz));
     167
     168    csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A);
     169//    csdn = fNucleon->GetIsoZACrossSection(&dp, Z, A);
     170    csdn = fNucleon->GetZandACrossSection(&dp, iz, A);
    172171
    173172    theGlauberFac[iz] = csdn/csup;
    … …  
    178177  fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
    179178  theCoulombFac[1] =
    180     fHadron->GetInelasticHadronNucleonXsc()/CoulombFactor(fLowEnergy,1.0);
     179    fHadron->GetInelasticHadronNucleonXsc()/CoulombFactor(fLowEnergy,1);
    181180     
    182181  for(G4int iz=2; iz<93; iz++) {
    183182
    184183    G4double Z = G4double(iz);
    185     A = nist->GetAtomicMassAmu(iz);
     184    A = G4lrint(nist->GetAtomicMassAmu(iz));
    186185
    187186    theCoulombFac[iz] =
    188       fNucleon->GetIsoZACrossSection(&dp, Z, A)/CoulombFactor(fLowEnergy,A);
     187      fNucleon->GetZandACrossSection(&dp, iz, A)/CoulombFactor(fLowEnergy,A);
    189188
    190189    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
    … …  
    193192}
    194193
    195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    196 
    197 G4double G4BGGNucleonInelasticXS::CoulombFactor(G4double kinEnergy, G4double A)
     194
     195G4double G4BGGNucleonInelasticXS::CoulombFactor(G4double kinEnergy, G4int A)
    198196{
    199197  G4double res= 0.0;
    200198  if(kinEnergy <= DBL_MIN) return res;
    201   else if (A < 1.5) return kinEnergy*kinEnergy;
     199  else if (A < 2) return kinEnergy*kinEnergy;
    202200 
    203201  G4double elog = std::log10(kinEnergy/GeV);
     202  G4double aa = A;
    204203
    205204  // from G4ProtonInelasticCrossSection
    206205  if(isProton) {
    207     G4double f1 = 8.0  - 8.0/A - 0.008*A;
    208     G4double f2 = 2.34 - 5.4/A - 0.0028*A;
     206    G4double f1 = 8.0  - 8.0/aa - 0.008*aa;
     207    G4double f2 = 2.34 - 5.4/aa - 0.0028*aa;
    209208
    210209    res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
    211210 
    212     f1 = 5.6 - 0.016*A;
    213     f2 = 1.37 + 1.37/A;
    214     res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
     211    f1 = 5.6 - 0.016*aa;
     212    f2 = 1.37 + 1.37/aa;
     213    res *= ( 1.0 + (0.8 + 18./aa - 0.002*aa)/(1.0 + std::exp(f1*(elog + f2))));
    215214  } else {
    216215
    217     G4double p3 = 0.6 + 13./A - 0.0005*A;
    218     G4double p4 = 7.2449 - 0.018242*A;
    219     G4double p5 = 1.36 + 1.8/A + 0.0005*A;
    220     G4double p6 = 1. + 200./A + 0.02*A;
    221     G4double p7 = 3.0 - (A-70.)*(A-200.)/11000.;
     216    G4double p3 = 0.6 + 13./aa - 0.0005*aa;
     217    G4double p4 = 7.2449 - 0.018242*aa;
     218    G4double p5 = 1.36 + 1.8/aa + 0.0005*aa;
     219    G4double p6 = 1. + 200./aa + 0.02*aa;
     220    G4double p7 = 3.0 - (aa-70.)*(aa-200.)/11000.;
    222221
    223222    res = (1.+p3/(1. + std::exp(p4*(elog+p5))))/(1.+std::exp(-p6*(elog+p7)));
    … …  
    227226}
    228227
    229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    230 

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BGGPionElasticXS.cc,v 1.7 2009/11/19 19:40:45 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BGGPionElasticXS.cc,v 1.8 2010/10/12 06:03:37 dennis Exp $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7777//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    7878
    79 G4double G4BGGPionElasticXS::GetIsoZACrossSection(const G4DynamicParticle* dp,
    80                                                     G4double Z,
    81                                                     G4double A,
    82                                                     G4double)
     79G4double
     80G4BGGPionElasticXS::GetZandACrossSection(const G4DynamicParticle* dp,
     81                                         G4int Z, G4int A, G4double)
    8382{
    8483  G4double cross = 0.0;
    8584  G4double ekin = dp->GetKineticEnergy();
    86   G4int iz = G4int(Z);
    87   if(iz > 92) iz = 92;
     85//  G4int iz = G4int(Z);
     86  if(Z > 92) Z = 92;
    8887
    8988  if(ekin <= fLowEnergy) {
    90     cross = theCoulombFac[iz];
     89    cross = theCoulombFac[Z];
    9190    if(isPiplus) { cross *= CoulombFactor(ekin, A); }
    9291
    93   } else if(iz == 1) {
    94     if( A < 1.5) {
     92  } else if(Z == 1) {
     93    if( A < 2) {
    9594      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
    9695      cross = fHadron->GetElasticHadronNucleonXsc();
    … …  
    102101    }
    103102  } else if(ekin > fGlauberEnergy) {
    104     cross = theGlauberFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
     103    cross = theGlauberFac[Z]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
    105104  } else {
    106105    cross = fPion->GetElasticCrossSection(dp, Z, A);
    … …  
    159158  G4NistManager* nist = G4NistManager::Instance();
    160159
    161   G4double A, csup, csdn;
     160  G4double csup, csdn;
     161  G4int A;
    162162
    163163  if(verboseLevel > 0) G4cout << "### G4BGGPionElasticXS::Initialise for "
    … …  
    167167
    168168    G4double Z = G4double(iz);
    169     A = nist->GetAtomicMassAmu(iz);
    170 
    171     csup = fGlauber->GetElasticGlauberGribov(&dp, Z, A);
    172     csdn = fPion->GetElasticCrossSection(&dp, Z, A);
     169    A = G4lrint(nist->GetAtomicMassAmu(iz));
     170
     171    csup = fGlauber->GetElasticGlauberGribov(&dp, iz, A);
     172    csdn = fPion->GetElasticCrossSection(&dp, iz, A);
    173173
    174174    theGlauberFac[iz] = csdn/csup;
    … …  
    179179  fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
    180180  theCoulombFac[1] = fHadron->GetElasticHadronNucleonXsc();
    181   if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1.0); }
     181  if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1); }
    182182
    183183  for(G4int iz=2; iz<93; iz++) {
    184184
    185185    G4double Z = G4double(iz);
    186     A = nist->GetAtomicMassAmu(iz);
    187 
    188     theCoulombFac[iz] = fPion->GetElasticCrossSection(&dp, Z, A);
     186    A = G4lrint(nist->GetAtomicMassAmu(iz));
     187
     188    theCoulombFac[iz] = fPion->GetElasticCrossSection(&dp, iz, A);
    189189    if(isPiplus) { theCoulombFac[iz] /= CoulombFactor(fLowEnergy,A); }
    190190    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A
    … …  
    193193}
    194194
    195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    196 
    197 G4double G4BGGPionElasticXS::CoulombFactor(G4double kinEnergy, G4double A)
     195
     196G4double G4BGGPionElasticXS::CoulombFactor(G4double kinEnergy, G4int A)
    198197{
    199198  G4double res= 0.0;
    200199  if(kinEnergy <= DBL_MIN) return res;
    201   else if(A < 1.5) return kinEnergy*kinEnergy;
     200  else if(A < 2) return kinEnergy*kinEnergy;
    202201 
    203202  G4double elog = std::log10(kinEnergy/GeV);
     203  G4double aa = A;
    204204
    205205  // from G4ProtonInelasticCrossSection
    206   G4double f1 = 8.0  - 8.0/A - 0.008*A;
    207   G4double f2 = 2.34 - 5.4/A - 0.0028*A;
     206  G4double f1 = 8.0  - 8.0/aa - 0.008*aa;
     207  G4double f2 = 2.34 - 5.4/aa - 0.0028*aa;
    208208
    209209  res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
    210210 
    211   f1 = 5.6 - 0.016*A;
    212   f2 = 1.37 + 1.37/A;
    213   res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
     211  f1 = 5.6 - 0.016*aa;
     212  f2 = 1.37 + 1.37/aa;
     213  res *= ( 1.0 + (0.8 + 18./aa - 0.002*aa)/(1.0 + std::exp(f1*(elog + f2))));
    214214  return res; 
    215215}
    216216
    217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    218 
    219 

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BGGPionInelasticXS.cc,v 1.8 2009/11/19 19:40:45 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BGGPionInelasticXS.cc,v 1.9 2010/10/12 06:03:49 dennis Exp $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7777//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    7878
    79 G4double G4BGGPionInelasticXS::GetIsoZACrossSection(const G4DynamicParticle* dp,
    80                                                     G4double Z,
    81                                                     G4double A,
    82                                                     G4double)
     79G4double
     80G4BGGPionInelasticXS::GetZandACrossSection(const G4DynamicParticle* dp,
     81                                           G4int Z, G4int A, G4double)
    8382{
    8483  G4double cross = 0.0;
    8584  G4double ekin = dp->GetKineticEnergy();
    86   G4int iz = G4int(Z);
    87   if(iz > 92) iz = 92;
     85//  G4int iz = G4int(Z);
     86  if(Z > 92) Z = 92;
    8887
    8988  if(ekin <= fLowEnergy) {
    90     cross = theCoulombFac[iz];
     89    cross = theCoulombFac[Z];
    9190    if(isPiplus) { cross *= CoulombFactor(ekin, A); }
    9291
    93   } else if(iz == 1) {
    94 
    95     if( A < 1.5) {
     92  } else if(Z == 1) {
     93
     94    if( A < 2) {
    9695      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
    9796      cross = fHadron->GetInelasticHadronNucleonXsc();
    … …  
    104103
    105104  } else if(ekin > fGlauberEnergy) {
    106     cross = theGlauberFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
     105    cross = theGlauberFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
    107106
    108107  } else {
    … …  
    162161  G4NistManager* nist = G4NistManager::Instance();
    163162
    164   G4double A, csup, csdn;
     163  G4double csup, csdn;
     164  G4int A;
    165165
    166166  if(verboseLevel > 0) G4cout << "### G4BGGPionInelasticXS::Initialise for "
    … …  
    172172
    173173    G4double Z = G4double(iz);
    174     A = nist->GetAtomicMassAmu(iz);
    175 
    176     csup = fGlauber->GetInelasticGlauberGribov(&dp, Z, A);
    177     csdn = fPion->GetInelasticCrossSection(&dp, Z, A);
     174    A = G4lrint(nist->GetAtomicMassAmu(iz));
     175
     176    csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A);
     177    csdn = fPion->GetInelasticCrossSection(&dp, iz, A);
    178178
    179179    theGlauberFac[iz] = csdn/csup;
    … …  
    184184  fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
    185185  theCoulombFac[1] = fHadron->GetInelasticHadronNucleonXsc();
    186   if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1.0); }
     186  if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1); }
    187187     
    188188  for(G4int iz=2; iz<93; iz++) {
    189189
    190190    G4double Z = G4double(iz);
    191     A = nist->GetAtomicMassAmu(iz);
    192 
    193     theCoulombFac[iz] = fPion->GetIsoZACrossSection(&dp, Z, A);
     191    A = G4lrint(nist->GetAtomicMassAmu(iz));
     192
     193    theCoulombFac[iz] = fPion->GetIsoZACrossSection(&dp, iz, A);
    194194    if(isPiplus) { theCoulombFac[iz] /= CoulombFactor(fLowEnergy,A); }
    195195
    … …  
    199199}
    200200
    201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    202 
    203 G4double G4BGGPionInelasticXS::CoulombFactor(G4double kinEnergy, G4double A)
     201
     202G4double G4BGGPionInelasticXS::CoulombFactor(G4double kinEnergy, G4int A)
    204203{
    205204  G4double res= 0.0;
    206205  if(kinEnergy <= DBL_MIN) return res;
    207   else if(A < 1.5) return kinEnergy*kinEnergy;
     206  else if(A < 2) return kinEnergy*kinEnergy;
    208207 
    209208  G4double elog = std::log10(kinEnergy/GeV);
     209  G4double aa = A;
    210210
    211211  // from G4ProtonInelasticCrossSection
    212   G4double f1 = 8.0  - 8.0/A - 0.008*A;
    213   G4double f2 = 2.34 - 5.4/A - 0.0028*A;
     212  G4double f1 = 8.0  - 8.0/aa - 0.008*aa;
     213  G4double f2 = 2.34 - 5.4/aa - 0.0028*aa;
    214214
    215215  res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
    216216 
    217   f1 = 5.6 - 0.016*A;
    218   f2 = 1.37 + 1.37/A;
    219   res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
     217  f1 = 5.6 - 0.016*aa;
     218  f2 = 1.37 + 1.37/aa;
     219  res *= ( 1.0 + (0.8 + 18./aa - 0.002*aa)/(1.0 + std::exp(f1*(elog + f2))));
    220220  return res; 
    221221}
    222222
    223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    224 
    225 
     223

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

    r819 r1340  
    3434// ********************************************************************
    3535//
    36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     36// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3737//
    3838// MODULE:              G4EMDissociationCrossSection.cc
    … …  
    4545// Contract:            17191/03/NL/LvH
    4646//
    47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     47// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4848//
    4949// CHANGE HISTORY
    … …  
    5656// Beta release
    5757//
    58 // 30. May 2005, J.P. Wellisch removed a compilation warning on gcc 3.4 for gant4 7.1.
    59 //
    60 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    61 ////////////////////////////////////////////////////////////////////////////////
     58// 30. May 2005, J.P. Wellisch removed a compilation warning on gcc 3.4 for
     59//               geant4 7.1.
     60//
     61// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     62//////////////////////////////////////////////////////////////////////////////
    6263//
    6364#include "G4EMDissociationCrossSection.hh"
    … …  
    6566#include "G4ParticleTable.hh"
    6667#include "G4IonTable.hh"
     68#include "G4HadTmpUtil.hh"
    6769#include "globals.hh"
    68 ////////////////////////////////////////////////////////////////////////////////
    69 //
     70
     71
    7072G4EMDissociationCrossSection::G4EMDissociationCrossSection ()
    7173{
    72 //
    7374//
    7475// This function makes use of the class which can sample the virtual photon
    … …  
    8687  xd      = 0.25;
    8788}
    88 ////////////////////////////////////////////////////////////////////////////////
     89//////////////////////////////////////////////////////////////////////////////
    8990//
    9091G4EMDissociationCrossSection::~G4EMDissociationCrossSection()
    … …  
    9293  delete thePhotonSpectrum;
    9394}
    94 ///////////////////////////////////////////////////////////////////////////////
    95 //
    96 G4bool G4EMDissociationCrossSection::IsZAApplicable
    97 (const G4DynamicParticle* theDynamicParticle, G4double /*ZZ*/, G4double AA)
     95/////////////////////////////////////////////////////////////////////////////
     96//
     97G4bool
     98G4EMDissociationCrossSection::IsIsoApplicable(const G4DynamicParticle* theDynamicParticle,
     99                                              G4int /*ZZ*/, G4int AA)
    98100{
    99101//
    … …  
    105107//
    106108  if (G4ParticleTable::GetParticleTable()->GetIonTable()->
    107     IsIon(theDynamicParticle->GetDefinition()) && AA > 1.0)
     109    IsIon(theDynamicParticle->GetDefinition()) && AA > 1)
    108110    return true;
    109111  else
    … …  
    114116  (const G4DynamicParticle* theDynamicParticle, const G4Element* theElement)
    115117{
    116   return IsZAApplicable(theDynamicParticle, 0., theElement->GetN());
     118  return IsIsoApplicable(theDynamicParticle, 0, G4lrint(theElement->GetN()));
    117119}
    118120
    … …  
    130132    G4IsotopeVector* isoVector = theElement->GetIsotopeVector();
    131133    G4double* abundVector = theElement->GetRelativeAbundanceVector();
    132     G4double ZZ;
    133     G4double AA;
     134    G4int ZZ;
     135    G4int AA;
    134136     
    135137    for (G4int i = 0; i < nIso; i++) {
    136       ZZ = G4double( (*isoVector)[i]->GetZ() );
    137       AA = G4double( (*isoVector)[i]->GetN() );
    138       sig = GetIsoZACrossSection(theDynamicParticle, ZZ, AA, temperature);
     138      ZZ = (*isoVector)[i]->GetZ();
     139      AA = (*isoVector)[i]->GetN();
     140      sig = GetZandACrossSection(theDynamicParticle, ZZ, AA, temperature);
    139141      crossSection += sig*abundVector[i];
    140142    }
    141143   
    142144  } else {
    143     crossSection =
    144       GetIsoZACrossSection(theDynamicParticle, theElement->GetZ(),
    145                            theElement->GetN(), temperature);
     145    G4int ZZ = G4lrint(theElement->GetZ());
     146    G4int AA = G4lrint(theElement->GetN());
     147    crossSection = GetZandACrossSection(theDynamicParticle, ZZ, AA, temperature);
    146148  }
    147149   
    … …  
    150152
    151153
    152 G4double G4EMDissociationCrossSection::GetIsoZACrossSection
    153   (const G4DynamicParticle *theDynamicParticle, G4double ZZ, G4double AA,
    154    G4double /*temperature*/)
    155 {
    156 //
     154G4double
     155G4EMDissociationCrossSection::GetZandACrossSection(const G4DynamicParticle *theDynamicParticle,
     156                                                   G4int ZZ, G4int AA, G4double /*temperature*/)
     157{
    157158//
    158159// Get relevant information about the projectile and target (A, Z) and
    … …  
    240241  return theCrossSectionVector;
    241242}
     243
    242244////////////////////////////////////////////////////////////////////////////////
    243245//
    … …  
    247249{
    248250//
    249 //
    250251// This is a cheaky little member function to calculate the probability of
    251252// EMD for the target in the field of the projectile ... just by reversing the
    … …  
    254255  return GetCrossSectionForProjectile (AT, ZT, AP, ZP, b, bmin);
    255256}
    256 ////////////////////////////////////////////////////////////////////////////////
    257 //
     257
     258
    258259G4double
    259   G4EMDissociationCrossSection::GetWilsonProbabilityForProtonDissociation
    260   (G4double A, G4double Z)
    261 {
    262 //
     260G4EMDissociationCrossSection::GetWilsonProbabilityForProtonDissociation(G4double A,
     261                                                                        G4double Z)
     262{
    263263//
    264264// This is a simple algorithm to choose whether a proton or neutron is ejected
    … …  
    282282  return p;
    283283}
    284 ////////////////////////////////////////////////////////////////////////////////
    285 //

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ElectroNuclearCrossSection.cc,v 1.29 2008/10/24 19:15:17 dennis Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ElectroNuclearCrossSection.cc,v 1.30 2010/10/14 05:25:22 dennis Exp $
     28// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2929//
    3030//
    … …  
    3333// The last update: M.V. Kossov, CERN/ITEP (Moscow) 17-Oct-03
    3434//
    35 //===============================================================================================
     35//============================================================================
    3636
    3737///#define debug
    … …  
    4343
    4444#include "G4ElectroNuclearCrossSection.hh"
     45#include "G4HadTmpUtil.hh"
    4546
    4647// Initialization of statics
    47 G4double  G4ElectroNuclearCrossSection::lastE=0.;  // Last used in the cross section TheEnergy
    48 G4int     G4ElectroNuclearCrossSection::lastF=0;   // Last used in the cross section TheFirstBin
    49 G4double  G4ElectroNuclearCrossSection::lastG=0.;  // Last used in the cross section TheGamma
     48// Last used in the cross section TheEnergy
     49G4double G4ElectroNuclearCrossSection::lastE=0.;
     50
     51// Last used in the cross section TheFirstBin
     52G4int G4ElectroNuclearCrossSection::lastF=0;
     53
     54// Last used in the cross section TheGamma
     55G4double G4ElectroNuclearCrossSection::lastG=0.;
     56
    5057G4double  G4ElectroNuclearCrossSection::lastH=0.;  // Last value of the High Energy A-dependence
    5158G4double* G4ElectroNuclearCrossSection::lastJ1=0;  // Pointer to the last array of the J1 function
    … …  
    5865G4double  G4ElectroNuclearCrossSection::lastSig=0.;// Last value of the Cross Section
    5966
    60 std::vector<G4double*> G4ElectroNuclearCrossSection::J1;     // Vector of pointers to the J1 tabulated functions
    61 std::vector<G4double*> G4ElectroNuclearCrossSection::J2;     // Vector of pointers to the J2 tabulated functions
    62 std::vector<G4double*> G4ElectroNuclearCrossSection::J3;     // Vector of pointers to the J3 tabulated functions
     67// Vector of pointers to the J1 tabulated functions
     68std::vector<G4double*> G4ElectroNuclearCrossSection::J1;
     69
     70// Vector of pointers to the J2 tabulated functions
     71std::vector<G4double*> G4ElectroNuclearCrossSection::J2;
     72
     73// Vector of pointers to the J3 tabulated functions
     74std::vector<G4double*> G4ElectroNuclearCrossSection::J3;
     75
    6376
    6477G4ElectroNuclearCrossSection::G4ElectroNuclearCrossSection()
    65 {
    66 }
     78{}
     79
    6780
    6881G4ElectroNuclearCrossSection::~G4ElectroNuclearCrossSection()
    … …  
    95108    G4IsotopeVector* isoVector = anEle->GetIsotopeVector();
    96109    G4double* abundVector = anEle->GetRelativeAbundanceVector();
    97     G4double ZZ;
    98     G4double AA;
     110    G4int ZZ;
     111    G4int AA;
    99112     
    100113    for (G4int i = 0; i < nIso; i++) {
    101       ZZ = G4double( (*isoVector)[i]->GetZ() );
    102       AA = G4double( (*isoVector)[i]->GetN() );
    103       sig = GetIsoZACrossSection(aPart, ZZ, AA, temperature);
     114      ZZ = (*isoVector)[i]->GetZ();
     115      AA = (*isoVector)[i]->GetN();
     116      sig = GetZandACrossSection(aPart, ZZ, AA, temperature);
    104117      xsection += sig*abundVector[i];
    105118    }
    … …  
    107120  } else {
    108121    xsection =
    109       GetIsoZACrossSection(aPart, anEle->GetZ(), anEle->GetN(),
     122      GetZandACrossSection(aPart,
     123                           G4lrint(anEle->GetZ()),
     124                           G4lrint(anEle->GetN()),
    110125                           temperature);
    111126  }
    … …  
    116131
    117132G4double
    118 G4ElectroNuclearCrossSection::GetIsoZACrossSection(const G4DynamicParticle* aPart,
    119                                               G4double ZZ, G4double AA,
    120                                               G4double /*temperature*/)
     133G4ElectroNuclearCrossSection::GetZandACrossSection(const G4DynamicParticle* aPart,
     134                                                   G4int ZZ, G4int AA,
     135                                                   G4double /*temperature*/)
    121136{
    122137  static const G4int nE=336; // !!  If you change this, change it in GetFunctions() (*.hh) !!
    123138  static const G4int mL=nE-1;
    124   static const G4double EMi=2.0612;          // Minimum tabulated Energy of the Electron
    125   static const G4double EMa=50000.;          // Maximum tabulated Energy of the Electron
    126   static const G4double lEMi=std::log(EMi);       // Minimum tabulated logarithmic Energy of the Electron
    127   static const G4double lEMa=std::log(EMa);       // Maximum tabulated logarithmic Energy of the Electron
     139  static const G4double EMi=2.0612; // Minimum tabulated Energy of the Electron
     140  static const G4double EMa=50000.; // Maximum tabulated Energy of the Electron
     141  static const G4double lEMi=std::log(EMi);  // Minimum tabulated logarithmic Energy of the Electron
     142  static const G4double lEMa=std::log(EMa);  // Maximum tabulated logarithmic Energy of the Electron
    128143  static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic step in the table for the electron Energy
    129144  static const G4double alop=1./137.036/3.14159265; //coef. for the calculated functions (Ee>50000.)
    … …  
    137152  static std::vector <G4double> colH;    // Vector of HighEnergyCoefficients (functional calculations)
    138153  // *** End of Static Definitions (Associative Memory) ***
     154
    139155  const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the electron
    140   const G4int targetAtomicNumber = static_cast<int>(AA+.499); //@@ Nat mixture (?!)
    141   const G4int targZ = static_cast<int>(ZZ+.001);
     156  const G4int targetAtomicNumber = AA;
     157  const G4int targZ = ZZ;
    142158  const G4int targN = targetAtomicNumber-targZ; // @@ Get isotops (can change initial A)
    143159  if (Energy<=EMi) return 0.;              // Energy is below the minimum energy in the table
     160
    144161  G4int PDG=aPart->GetDefinition()->GetPDGEncoding();
    145   if(  PDG == 11 || PDG == -11)            // @@ Now only for electrons, but can be fo muons
     162  if (PDG == 11 || PDG == -11)            // @@ Now only for electrons, but can be fo muons
    146163  {
    147     G4double A=targN+targZ;                // New A (can differ from G4double targetAtomicNumber)
    148     if(targN!=lastN || targZ!=lastZ)       // This nucleus was not the last used isotop
     164    G4double A = targN + targZ;           // New A (can differ from G4double targetAtomicNumber)
     165    if(targN!=lastN || targZ!=lastZ)      // This nucleus was not the last used isotop
    149166        {
    150167      lastE    = 0.;                       // New history in the electron Energy
    … …  
    185202          } // End of creation of the new set of parameters
    186203    } // End of parameters udate
     204
    187205    else if(std::abs((lastE-Energy)/Energy)<.001) return lastSig*millibarn; // Don't calc. same CS twice
    188206    // ============================== NOW Calculate the Cross Section ==========================
    … …  
    268286}
    269287
    270 // Correction function for Be,C @@ Move to header // @@@ !!! NOT used !!!
    271 //G4double G4ElectroNuclearCrossSection::LinearFit(G4double X, G4int N, const G4double* XN,
    272 //                                                 const G4double* YN)
    273 //{
    274 //  G4double Xj=XN[0];
    275 //  G4double Xh=XN[N-1];
    276 //  if(X<=Xj) return YN[0]; //-------+ !!! If you correct this, correct upper copy too!!
    277 //  else if(X>=Xh) return YN[N-1];//-|
    278 //  G4double Xp=0.; //               |
    279 //  G4int j=0;   //                  |
    280 //  while (X>Xj && j<N)//<-----------+
    281 //  {
    282 //    j++;
    283 //    Xp=Xj;
    284 //    Xj=XN[j];
    285 //  }
    286 //  return YN[j]-(Xj-X)*(YN[j]-YN[j-1])/(Xj-Xp);
    287 //}
    288288
    289289// Calculate the functions for the std::log(A)

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

    r1228 r1340  
    3333#include "G4IonTable.hh"
    3434#include "G4ParticleDefinition.hh"
    35 
    36 
    37 
    38 ////////////////////////////////////////////////////////////////////////////////
     35#include "G4HadTmpUtil.hh"
     36
     37
     38///////////////////////////////////////////////////////////////////////////////
    3939//
    4040//
    … …  
    4949}
    5050
    51 ///////////////////////////////////////////////////////////////////////////////////////
     51///////////////////////////////////////////////////////////////////////////////
    5252//
    5353//
    5454
    5555G4GGNuclNuclCrossSection::~G4GGNuclNuclCrossSection()
    56 {
    57 }
    58 
    59 
    60 ////////////////////////////////////////////////////////////////////////////////////////
     56{}
     57
     58///////////////////////////////////////////////////////////////////////////////
    6159//
    6260//
    … …  
    6563G4bool
    6664G4GGNuclNuclCrossSection::IsApplicable(const G4DynamicParticle* aDP,
    67                                           const G4Element*  anElement)
    68 {
    69   return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN());
     65                                       const G4Element*  anElement)
     66{
     67  G4int Z = G4lrint(anElement->GetZ());
     68  G4int N = G4lrint(anElement->GetN());
     69  return IsIsoApplicable(aDP, Z, N);
    7070}
    7171
    72 ////////////////////////////////////////////////////////////////////////////////////////
     72///////////////////////////////////////////////////////////////////////////////
    7373//
    7474//
    7575
    7676G4bool
    77 G4GGNuclNuclCrossSection::IsZAApplicable(const G4DynamicParticle* aDP,
    78                                             G4double Z, G4double)
    79 {
    80   G4bool applicable      = false;
    81   // G4int baryonNumber     = aDP->GetDefinition()->GetBaryonNumber();
     77G4GGNuclNuclCrossSection::IsIsoApplicable(const G4DynamicParticle* aDP,
     78                                          G4int Z, G4int)
     79{
     80  G4bool applicable = false;
    8281  G4double kineticEnergy = aDP->GetKineticEnergy();
    8382
    84   //  const G4ParticleDefinition* theParticle = aDP->GetDefinition();
    85  
    86   if ( kineticEnergy  >= fLowerLimit && Z > 1.5 ) applicable = true;
    87 
     83  if (kineticEnergy >= fLowerLimit && Z > 1) applicable = true;
    8884  return applicable;
    8985}
    9086
    91 ////////////////////////////////////////////////////////////////////////////////////////
    92 //
    93 // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
    94 // Glauber model with Gribov correction calculated in the dipole approximation on
    95 // light cone. Gaussian density helps to calculate rest integrals of the model.
    96 // [1] B.Z. Kopeliovich, nucl-th/0306044
     87///////////////////////////////////////////////////////////////////////////////
     88//
     89// Calculates total and inelastic Xsc, derives elastic as total - inelastic
     90// accordong to Glauber model with Gribov correction calculated in the dipole
     91// approximation on light cone. Gaussian density helps to calculate rest
     92// integrals of the model. [1] B.Z. Kopeliovich, nucl-th/0306044
    9793
    9894
    9995G4double G4GGNuclNuclCrossSection::
    100 GetCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement, G4double T)
    101 {
    102   return GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(), T);
    103 }
    104 
    105 ////////////////////////////////////////////////////////////////////////////////////////
    106 //
    107 // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
    108 // Glauber model with Gribov correction calculated in the dipole approximation on
    109 // light cone. Gaussian density of point-like nucleons helps to calculate rest integrals of the model.
    110 // [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above
    111 
     96GetCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement,
     97                G4double T)
     98{
     99  G4int Z = G4lrint(anElement->GetZ());
     100  G4int N = G4lrint(anElement->GetN());
     101  return GetZandACrossSection(aParticle, Z, N, T);
     102}
     103
     104///////////////////////////////////////////////////////////////////////////////
     105//
     106// Calculates total and inelastic Xsc, derives elastic as total - inelastic
     107// accordong to Glauber model with Gribov correction calculated in the dipole
     108// approximation on light cone. Gaussian density of point-like nucleons helps
     109// to calculate rest integrals of the model. [1] B.Z. Kopeliovich,
     110// nucl-th/0306044 + simplification above
    112111
    113112
    114113G4double G4GGNuclNuclCrossSection::
    115 GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double tZ, G4double tA, G4double)
    116 {
    117   G4double xsection, sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, cB, ratio;
     114GetZandACrossSection(const G4DynamicParticle* aParticle,
     115                     G4int tZ, G4int tA, G4double)
     116{
     117  G4double xsection;
     118  G4double sigma;
     119  G4double cofInelastic = 2.4;
     120  G4double cofTotal = 2.0;
     121  G4double nucleusSquare;
     122  G4double cB;
     123  G4double ratio;
    118124
    119125  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
    … …  
    132138  G4double pR = GetNucleusRadius(pA);
    133139
    134   cB = GetCoulombBarier(aParticle, tZ, tA, pR, tR);
    135 
    136   if(cB > 0.)
    137   {
     140  cB = GetCoulombBarier(aParticle, G4double(tZ), G4double(tA), pR, tR);
     141  if (cB > 0.) {
    138142
    139143    sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
    … …  
    143147
    144148    ratio = sigma/nucleusSquare;
    145 
    146149    xsection =  nucleusSquare*std::log( 1. + ratio );
    147 
    148150    fTotalXsc = xsection;
    149 
    150151    fTotalXsc *= cB;
    151152
    … …  
    153154
    154155    fInelasticXsc *= cB;
    155 
    156156    fElasticXsc   = fTotalXsc - fInelasticXsc;
    157157
    … …  
    171171 
    172172    ratio = sigma/nucleusSquare;
    173 
    174173    fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
    175174
    … …  
    186185}
    187186
    188 ////////////////////////////////////////////////////////////////////////////////////////
     187///////////////////////////////////////////////////////////////////////////////
    189188//
    190189//
    191190
    192191G4double G4GGNuclNuclCrossSection::
    193 GetCoulombBarier(const G4DynamicParticle* aParticle, G4double tZ, G4double tA, G4double pR, G4double tR)
     192GetCoulombBarier(const G4DynamicParticle* aParticle, G4double tZ, G4double tA,
     193                 G4double pR, G4double tR)
    194194{
    195195  G4double ratio;
    196196  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
    197 
    198197
    199198  G4double pTkin = aParticle->GetKineticEnergy();
    … …  
    254253
    255254  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );   // basically 2piRR
    256 
    257255  ratio = sigma/nucleusSquare;
    258 
    259 
    260   fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
    261    
     256  fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
    262257  G4double difratio = ratio/(1.+ratio);
    263258
    … …  
    272267//////////////////////////////////////////////////////////////////////////
    273268//
    274 // Return suasi-elastic/inelastic cross-section ratio
     269// Return quasi-elastic/inelastic cross-section ratio
    275270
    276271G4double G4GGNuclNuclCrossSection::
    … …  
    298293
    299294  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );   // basically 2piRR
    300 
    301295  ratio = sigma/nucleusSquare;
    302 
    303   fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
     296  fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
    304297
    305298  //  sigma = GetHNinelasticXsc(aParticle, tA, tZ);
    306299  ratio = sigma/nucleusSquare;
    307 
    308   fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
     300  fProductionXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
    309301
    310302  if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
    … …  
    315307}
    316308
    317 /////////////////////////////////////////////////////////////////////////////////////
     309///////////////////////////////////////////////////////////////////////////////
    318310//
    319311// Returns hadron-nucleon Xsc according to differnt parametrisations:
    … …  
    324316G4double
    325317G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
    326                                                   const G4Element* anElement          )
    327 {
    328   G4double At = anElement->GetN();  // number of nucleons
    329   G4double Zt = anElement->GetZ();  // number of protons
    330 
    331 
    332   return GetHadronNucleonXsc( aParticle, At, Zt );
    333 }
    334 
    335 
    336 
    337 
    338 /////////////////////////////////////////////////////////////////////////////////////
     318                                              const G4Element* anElement)
     319{
     320  G4int At = G4lrint(anElement->GetN());  // number of nucleons
     321  G4int Zt = G4lrint(anElement->GetZ());  // number of protons
     322  return GetHadronNucleonXsc(aParticle, At, Zt);
     323}
     324
     325
     326
     327
     328///////////////////////////////////////////////////////////////////////////////
    339329//
    340330// Returns hadron-nucleon Xsc according to differnt parametrisations:
    … …  
    345335G4double
    346336G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
    347                                                    G4double At,  G4double Zt       )
     337                                                   G4int At, G4int Zt)
    348338{
    349339  G4double xsection = 0.;
    350340
    351 
    352341  G4double targ_mass = G4ParticleTable::GetParticleTable()->
    353   GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
    354 
     342  GetIonTable()->GetIonMass(Zt, At);
    355343  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
    356344
    357   G4double proj_mass     = aParticle->GetMass();
     345  G4double proj_mass = aParticle->GetMass();
    358346  G4double proj_momentum = aParticle->GetMomentum().mag();
    359347  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
    … …  
    361349  sMand /= GeV*GeV;  // in GeV for parametrisation
    362350  proj_momentum /= GeV;
    363 
    364351  const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
    365  
    366352
    367353  if(pParticle == theNeutron) // as proton ???
    368354  {
    369     xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
     355    xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    370356  }
    371357  else if(pParticle == theProton)
    372358  {
    373     xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    374     // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
    375     // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
     359    xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    376360  }
    377361 
    … …  
    381365
    382366
    383 /////////////////////////////////////////////////////////////////////////////////////
     367///////////////////////////////////////////////////////////////////////////////
    384368//
    385369// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
    … …  
    388372G4double
    389373G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
    390                                                   const G4Element* anElement          )
    391 {
    392   G4double At = anElement->GetN();  // number of nucleons
    393   G4double Zt = anElement->GetZ();  // number of protons
    394 
    395 
     374                                                  const G4Element* anElement)
     375{
     376  G4int At = G4lrint(anElement->GetN());  // number of nucleons
     377  G4int Zt = G4lrint(anElement->GetZ());  // number of protons
    396378  return GetHadronNucleonXscPDG( aParticle, At, Zt );
    397379}
    398380
    399381
    400 
    401 
    402 /////////////////////////////////////////////////////////////////////////////////////
     382///////////////////////////////////////////////////////////////////////////////
    403383//
    404384// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
    … …  
    408388G4double
    409389G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
    410                                                      G4double At,  G4double Zt )
     390                                                 G4int At, G4int Zt)
    411391{
    412392  G4double xsection = 0.;
    … …  
    415395  if (Nt < 0.) Nt = 0.; 
    416396
    417 
    418397  G4double targ_mass = G4ParticleTable::GetParticleTable()->
    419   GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
    420 
     398  GetIonTable()->GetIonMass(Zt, At);
    421399  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
    422400
    423401  G4double proj_mass     = aParticle->GetMass();
    424402  G4double proj_momentum = aParticle->GetMomentum().mag();
    425 
    426403  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
    427 
    428404  sMand         /= GeV*GeV;  // in GeV for parametrisation
    429405
    … …  
    435411  G4double B    = 0.308;
    436412
    437 
    438413  const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
    439414 
    440 
    441415  if(pParticle == theNeutron) // proton-neutron fit
    442416  {
    443     xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    444                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    445     xsection  += Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    446                       + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
     417    xsection = G4double(Zt)*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
     418                  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
     419    xsection += Nt*(35.45 + B*std::pow(std::log(sMand/s0),2.)
     420             + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
    447421  }
    448422  else if(pParticle == theProton)
    449423  {
    450      
    451       xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    452                           + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
    453 
    454       xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    455                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
     424    xsection = G4double(Zt)*(35.45 + B*std::pow(std::log(sMand/s0),2.)
     425                + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
     426
     427    xsection += Nt*(35.80 + B*std::pow(std::log(sMand/s0),2.)
     428                  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    456429  }
    457430  xsection *= millibarn; // parametrised in mb
    … …  
    460433
    461434
    462 
    463 
    464 /////////////////////////////////////////////////////////////////////////////////////
     435///////////////////////////////////////////////////////////////////////////////
    465436//
    466437// Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
    … …  
    469440
    470441G4double
    471 G4GGNuclNuclCrossSection::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
     442G4GGNuclNuclCrossSection::GetHadronNucleonXscNS(G4ParticleDefinition* pParticle,
    472443                                                 G4double pTkin,
    473444                                                 G4ParticleDefinition* tParticle)
    … …  
    477448  G4double hnXsc(0);
    478449
    479 
    480   G4double targ_mass     = tParticle->GetPDGMass();
    481   G4double proj_mass     = pParticle->GetPDGMass();
     450  G4double targ_mass = tParticle->GetPDGMass();
     451  G4double proj_mass = pParticle->GetPDGMass();
    482452
    483453  G4double proj_energy   = proj_mass + pTkin;
    … …  
    497467  //  G4double eta2 = 0.458;
    498468  //  G4double B    = 0.308;
    499 
    500  
    501 
    502 
    503469 
    504470  if( proj_momentum >= 10. ) // high energy: pp = nn = np
    … …  
    506472  {
    507473    Delta = 1.;
    508 
    509     if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
    510 
    511     if( proj_momentum >= 10.)
    512     {
    513         B0 = 7.5;
    514         A0 = 100. - B0*std::log(3.0e7);
    515 
    516         xsection = A0 + B0*std::log(proj_energy) - 11
    517                   + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
    518                      0.93827*0.93827,-0.165);        //  mb
     474    if (proj_energy < 40.) Delta = 0.916+0.0021*proj_energy;
     475
     476    if (proj_momentum >= 10.) {
     477      B0 = 7.5;
     478      A0 = 100. - B0*std::log(3.0e7);
     479
     480      xsection = A0 + B0*std::log(proj_energy) - 11
     481                + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
     482                  0.93827*0.93827,-0.165);        //  mb
    519483    }
    520484  }
    … …  
    562526}
    563527
    564 /*
    565 /////////////////////////////////////////////////////////////////////////////////////
    566 //
    567 // Returns hadron-nucleon inelastic cross-section based on proper parametrisation
    568 
    569 G4double
    570 G4GGNuclNuclCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
    571                                                   const G4Element* anElement          )
    572 {
    573   G4double At = anElement->GetN();  // number of nucleons
    574   G4double Zt = anElement->GetZ();  // number of protons
    575 
    576 
    577   return GetHNinelasticXsc( aParticle, At, Zt );
    578 }
    579 
    580 /////////////////////////////////////////////////////////////////////////////////////
     528/////////////////////////////////////////////////////////////////////////////////
    581529//
    582530// Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
    583531
    584532G4double
    585 G4GGNuclNuclCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
    586                                                      G4double At,  G4double Zt )
    587 {
    588   //  G4ParticleDefinition* hadron = aParticle->GetDefinition();
    589   G4double sumInelastic, Nt = At - Zt;
    590 
    591   if(Nt < 0.) Nt = 0.;
    592  
    593   sumInelastic  = Zt*GetHadronNucleonXscNS(aParticle, theProton);
    594   sumInelastic += Nt*GetHadronNucleonXscNS(aParticle, theNeutron);   
    595  
    596   return sumInelastic;
    597 }
    598 */
    599 
    600 /////////////////////////////////////////////////////////////////////////////////////
    601 //
    602 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
    603 
    604 G4double
    605533G4GGNuclNuclCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
    606                                                      G4double At,  G4double Zt )
    607 {
    608   G4int PDGcode    = aParticle->GetDefinition()->GetPDGEncoding();
     534                                              G4int At, G4int Zt)
     535{
     536  G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
    609537  G4int absPDGcode = std::abs(PDGcode);
    610 
    611538  G4double Elab = aParticle->GetTotalEnergy();             
    612539                          // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
    … …  
    650577                           NumberOfTargetNeutrons * XelPN   );
    651578  }
     579
    652580  Xinelastic = Xtotal - Xelastic;
    653 
    654581  if(Xinelastic < 0.) Xinelastic = 0.;
    655582
    … …  
    657584}
    658585
    659 ////////////////////////////////////////////////////////////////////////////////////
    660 //
    661 //
    662 
    663 G4double
    664 G4GGNuclNuclCrossSection::GetNucleusRadius( const G4DynamicParticle* ,
    665                                                const G4Element* anElement)
    666 {
    667   G4double At       = anElement->GetN();
     586///////////////////////////////////////////////////////////////////////////////
     587//
     588//
     589
     590G4double
     591G4GGNuclNuclCrossSection::GetNucleusRadius(const G4DynamicParticle* ,
     592                                           const G4Element* anElement)
     593{
     594  G4double At = anElement->GetN();
    668595  G4double oneThird = 1.0/3.0;
    669596  G4double cubicrAt = std::pow (At, oneThird);
    670597
    671 
    672598  G4double R;  // = fRadiusConst*cubicrAt;
    673   /* 
    674   G4double tmp = std::pow( cubicrAt-1., 3.);
    675   tmp         += At;
    676   tmp         *= 0.5;
    677 
    678   if (At > 20.)   // 20.
    679   {
    680     R = fRadiusConst*std::pow (tmp, oneThird);
    681   }
    682   else
    683   {
    684     R = fRadiusConst*cubicrAt;
    685   }
    686   */
    687  
    688599  R = fRadiusConst*cubicrAt;
    689600
    690   // return R;  // !!!!
    691 
    692 
    693  
    694601  G4double meanA  = 21.;
    695 
    696602  G4double tauA1  = 40.;
    697603  G4double tauA2  = 10.;
    … …  
    715621  {
    716622    R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
    717   } 
     623  }
     624
    718625  return R;
    719  
    720 }
    721 
    722 ////////////////////////////////////////////////////////////////////////////////////
     626}
     627
     628///////////////////////////////////////////////////////////////////////////////
    723629//
    724630//
    … …  
    728634{
    729635  G4double R;
    730 
    731   // R = GetNucleusRadiusGG(At);
    732 
    733636  R = GetNucleusRadiusDE(At);
    734637
    … …  
    741644G4GGNuclNuclCrossSection::GetNucleusRadiusGG(G4double At)
    742645{
    743  
    744646  G4double oneThird = 1.0/3.0;
    745647  G4double cubicrAt = std::pow (At, oneThird);
    746648
    747 
    748   G4double R;  // = fRadiusConst*cubicrAt;
    749  
    750   /*
    751   G4double tmp = std::pow( cubicrAt-1., 3.);
    752   tmp         += At;
    753   tmp         *= 0.5;
    754 
    755   if (At > 20.)
    756   {
    757     R = fRadiusConst*std::pow (tmp, oneThird);
     649  G4double R;  // = fRadiusConst*cubicrAt; 
     650  R = fRadiusConst*cubicrAt;
     651
     652  G4double meanA = 20.;
     653  G4double tauA  = 20.;
     654
     655  if ( At > 20.)   // 20.
     656  {
     657    R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
    758658  }
    759659  else
    760660  {
    761     R = fRadiusConst*cubicrAt;
    762   }
    763   */
    764  
    765   R = fRadiusConst*cubicrAt;
    766 
    767   G4double meanA = 20.;
    768   G4double tauA  = 20.;
    769 
    770   if ( At > 20.)   // 20.
    771   {
    772     R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
    773   }
    774   else
    775   {
    776661    R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) );
    777662  }
    778663
    779664  return R;
    780 
    781665}
    782666
    … …  
    785669G4GGNuclNuclCrossSection::GetNucleusRadiusDE(G4double A)
    786670{
    787 
    788671  // algorithm from diffuse-elastic
    789672
    … …  
    797680
    798681
    799   if( A < 50. )
     682  if (A < 50.)
    800683  {
    801684    if( 10  < A && A <= 15. )     r0  = a11*( 1 - std::pow(A, -2./3.) )*fermi;   // 1.08*fermi;
    … …  
    816699
    817700
    818 ////////////////////////////////////////////////////////////////////////////////////
    819 //
    820 //
    821 
    822 G4double G4GGNuclNuclCrossSection::CalculateEcmValue( const G4double mp ,
    823                                                          const G4double mt ,
    824                                                          const G4double Plab )
     701///////////////////////////////////////////////////////////////////////////////
     702//
     703//
     704
     705G4double G4GGNuclNuclCrossSection::CalculateEcmValue(const G4double mp,
     706                                                     const G4double mt,
     707                                                     const G4double Plab)
    825708{
    826709  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
    … …  
    833716
    834717
    835 ////////////////////////////////////////////////////////////////////////////////////
    836 //
    837 //
    838 
    839 G4double G4GGNuclNuclCrossSection::CalcMandelstamS( const G4double mp ,
    840                                                        const G4double mt ,
    841                                                        const G4double Plab )
     718///////////////////////////////////////////////////////////////////////////////
     719//
     720//
     721
     722G4double G4GGNuclNuclCrossSection::CalcMandelstamS(const G4double mp,
     723                                                   const G4double mt,
     724                                                   const G4double Plab)
    842725{
    843726  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
    … …  
    847730}
    848731
    849 
    850 //
    851 //
    852 ///////////////////////////////////////////////////////////////////////////////////////
     732//
     733//
     734///////////////////////////////////////////////////////////////////////////////

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

    r819 r1340  
    6565#include "G4ParticleTable.hh"
    6666#include "G4IonTable.hh"
     67#include "G4HadTmpUtil.hh"
    6768
    6869#include <iomanip>
    69 ///////////////////////////////////////////////////////////////////////////////
    70 //
     70
     71
    7172G4GeneralSpaceNNCrossSection::G4GeneralSpaceNNCrossSection ()
    7273{
    … …  
    7980  return;
    8081}
    81 ///////////////////////////////////////////////////////////////////////////////
    82 //
     82
     83
    8384G4GeneralSpaceNNCrossSection::~G4GeneralSpaceNNCrossSection ()
    8485{
    … …  
    9293//
    9394
    94 G4bool G4GeneralSpaceNNCrossSection::IsZAApplicable
    95  (const G4DynamicParticle* theProjectile, G4double ZZ, G4double AA)
    96 {
    97   G4bool result = protonInelastic->IsZAApplicable(theProjectile, ZZ, AA);
     95G4bool G4GeneralSpaceNNCrossSection::IsIsoApplicable
     96 (const G4DynamicParticle* theProjectile, G4int ZZ, G4int AA)
     97{
     98  G4bool result = protonInelastic->IsIsoApplicable(theProjectile, ZZ, AA);
    9899  if (!result)
    99100  {
    100     result = ionProton->IsZAApplicable(theProjectile, ZZ, AA);
     101    result = ionProton->IsIsoApplicable(theProjectile, ZZ, AA);
    101102    if (!result)
    102103    {
    103       result = TripathiGeneral->IsZAApplicable(theProjectile, ZZ, AA);
     104      result = TripathiGeneral->IsIsoApplicable(theProjectile, ZZ, AA);
    104105      if (!result)
    105         result = Shen->IsZAApplicable(theProjectile, ZZ, AA);
     106        result = Shen->IsIsoApplicable(theProjectile, ZZ, AA);
    106107    }
    107108  }
    … …  
    139140    G4IsotopeVector* isoVector = theTarget->GetIsotopeVector();
    140141    G4double* abundVector = theTarget->GetRelativeAbundanceVector();
    141     G4double ZZ;
    142     G4double AA;
     142    G4int ZZ;
     143    G4int AA;
    143144     
    144145    for (G4int i = 0; i < nIso; i++) {
    145       ZZ = G4double( (*isoVector)[i]->GetZ() );
    146       AA = G4double( (*isoVector)[i]->GetN() );
    147       sig = GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     146      ZZ = (*isoVector)[i]->GetZ();
     147      AA = (*isoVector)[i]->GetN();
     148      sig = GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    148149      xsection += sig*abundVector[i];
    149150    }
    150151   
    151152  } else {
    152     xsection =
    153       GetIsoZACrossSection(theProjectile, theTarget->GetZ(), theTarget->GetN(),
    154                           theTemperature);
     153    G4int ZZ = G4lrint(theTarget->GetZ());
     154    G4int AA = G4lrint(theTarget->GetN());
     155    xsection = GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    155156  }
    156157   
    … …  
    159160
    160161
    161 G4double G4GeneralSpaceNNCrossSection::GetIsoZACrossSection
    162   (const G4DynamicParticle* theProjectile, G4double ZZ, G4double AA,
    163    G4double theTemperature)
     162G4double
     163G4GeneralSpaceNNCrossSection::GetZandACrossSection(const G4DynamicParticle* theProjectile,
     164                                                   G4int ZZ, G4int AA, G4double theTemperature)
    164165{
    165166  G4double result = 0.0;
    … …  
    186187    {
    187188      result = protonInelastic->
    188         GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     189        GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    189190      if (verboseLevel >= 2)
    190191        G4cout <<"Selecting G4ProtonInelasticCrossSection" <<G4endl;
    … …  
    193194    {
    194195      result = TripathiLight->
    195         GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     196        GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    196197      if (verboseLevel >= 2)
    197198        G4cout <<"Selecting G4TripathiLightCrossSection" <<G4endl;
    … …  
    203204    {
    204205      result = ionProton->
    205         GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     206        GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    206207      if (verboseLevel >= 2)
    207208        G4cout <<"Selecting G4IonProtonCrossSection" <<G4endl;
    … …  
    210211    {
    211212      result = TripathiLight->
    212         GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     213        GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    213214      if (verboseLevel >= 2)
    214215        G4cout <<"Selecting G4TripathiLightCrossSection" <<G4endl;
    … …  
    217218  else
    218219  {
    219     if (TripathiLight->IsZAApplicable(theProjectile, ZZ, AA))
     220    if (TripathiLight->IsIsoApplicable(theProjectile, ZZ, AA))
    220221    {
    221222      result = TripathiLight->
    222         GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     223        GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    223224      if (verboseLevel >= 2)
    224225        G4cout <<"Selecting G4TripathiLightCrossSection" <<G4endl;
    225226    }
    226     else if (TripathiGeneral->IsZAApplicable(theProjectile, ZZ, AA))
     227    else if (TripathiGeneral->IsIsoApplicable(theProjectile, ZZ, AA))
    227228    {
    228229      result = TripathiGeneral->
    229         GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     230        GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    230231      if (verboseLevel >= 2)
    231232        G4cout <<"Selecting G4TripathiCrossSection" <<G4endl;
    232233    }
    233     else if (Shen->IsZAApplicable(theProjectile, ZZ, AA))
     234    else if (Shen->IsIsoApplicable(theProjectile, ZZ, AA))
    234235    {
    235236      result = Shen->
    236         GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     237        GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    237238      if (verboseLevel >= 2)
    238239        G4cout <<"Selecting G4IonsShenCrossSection" <<G4endl;
    … …  
    247248  return result;
    248249}
    249 ///////////////////////////////////////////////////////////////////////////////
    250 //
     250

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

    r1337 r1340  
    279279                                          const G4Element*  anElement)
    280280{
    281   return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN());
     281  return IsIsoApplicable(aDP, G4lrint(anElement->GetZ()),
     282                              G4lrint(anElement->GetN()));
    282283}
    283284
    … …  
    287288
    288289G4bool
    289 G4GlauberGribovCrossSection::IsZAApplicable(const G4DynamicParticle* aDP,
    290                                             G4double Z, G4double)
     290G4GlauberGribovCrossSection::IsIsoApplicable(const G4DynamicParticle* aDP,
     291                                             G4int Z, G4int)
    291292{
    292293  G4bool applicable      = false;
    … …  
    297298 
    298299  if ( ( kineticEnergy  >= fLowerLimit &&
    299          Z > 1.5 &&      // >=  He
     300         Z > 1 &&      // >=  He
    300301       ( theParticle == theAProton   ||
    301302         theParticle == theGamma     ||
    … …  
    305306
    306307       ( kineticEnergy  >= fLowerLimit &&
    307          Z > 1.5 &&      // >=  He
     308         Z > 1 &&      // >=  He
    308309       ( theParticle == theProton    ||
    309310         theParticle == theNeutron   ||   
    … …  
    325326GetCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement, G4double T)
    326327{
    327   return GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(), T);
     328  return GetZandACrossSection(aParticle, G4lrint(anElement->GetZ()),
     329                              G4lrint(anElement->GetN()), T);
    328330}
    329331
    … …  
    338340
    339341G4double G4GlauberGribovCrossSection::
    340 GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double Z, G4double A, G4double)
     342GetZandACrossSection(const G4DynamicParticle* aParticle, G4int Z, G4int A, G4double)
    341343{
    342344  G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
    … …  
    362364  // cofInelastic = 2.0;
    363365
    364   if( A > 1.5 )
     366  if( A > 1 )
    365367  {
    366368    nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
    … …  
    420422
    421423G4double G4GlauberGribovCrossSection::
    422 GetRatioSD(const G4DynamicParticle* aParticle, G4double A, G4double Z)
     424GetRatioSD(const G4DynamicParticle* aParticle, G4int A, G4int Z)
    423425{
    424426  G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
    … …  
    462464
    463465G4double G4GlauberGribovCrossSection::
    464 GetRatioQE(const G4DynamicParticle* aParticle, G4double A, G4double Z)
     466GetRatioQE(const G4DynamicParticle* aParticle, G4int A, G4int Z)
    465467{
    466468  G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
    … …  
    510512G4double
    511513G4GlauberGribovCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
    512                                                   const G4Element* anElement          )
    513 {
    514   G4double At = anElement->GetN();  // number of nucleons
    515   G4double Zt = anElement->GetZ();  // number of protons
    516 
    517 
    518   return GetHadronNucleonXsc( aParticle, At, Zt );
     514                                                 const G4Element* anElement)
     515{
     516  G4int At = G4lrint(anElement->GetN());  // number of nucleons
     517  G4int Zt = G4lrint(anElement->GetZ());  // number of protons
     518
     519  return GetHadronNucleonXsc(aParticle, At, Zt);
    519520}
    520521
    … …  
    531532G4double
    532533G4GlauberGribovCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
    533                                                    G4double At,  G4double Zt       )
     534                                                 G4int At, G4int Zt)
    534535{
    535536  G4double xsection;
    536537
    537 
    538538  G4double targ_mass = G4ParticleTable::GetParticleTable()->
    539   GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
     539  GetIonTable()->GetIonMass(Zt, At);
     540//  GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
    540541
    541542  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
    … …  
    550551  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
    551552 
     553  G4double aa = At;
    552554
    553555  if(theParticle == theGamma)
    554556  {
    555     xsection = At*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
     557    xsection = aa*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
    556558  }
    557559  else if(theParticle == theNeutron) // as proton ???
    558560  {
    559     xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
     561    xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    560562  }
    561563  else if(theParticle == theProton)
    562564  {
    563     xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
     565    xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    564566    // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
    565567    // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
    … …  
    567569  else if(theParticle == theAProton)
    568570  {
    569     xsection = At*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
     571    xsection = aa*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
    570572  }
    571573  else if(theParticle == thePiPlus)
    572574  {
    573     xsection = At*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
     575    xsection = aa*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
    574576  }
    575577  else if(theParticle == thePiMinus)
    576578  {
    577579    // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
    578     xsection = At*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
     580    xsection = aa*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
    579581  }
    580582  else if(theParticle == theKPlus)
    581583  {
    582     xsection = At*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
     584    xsection = aa*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
    583585  }
    584586  else if(theParticle == theKMinus)
    585587  {
    586     xsection = At*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
     588    xsection = aa*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
    587589  }
    588590  else  // as proton ???
    589591  {
    590     xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
     592    xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    591593  }
    592594  xsection *= millibarn;
    … …  
    602604G4double
    603605G4GlauberGribovCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
    604                                                   const G4Element* anElement          )
    605 {
    606   G4double At = anElement->GetN();  // number of nucleons
    607   G4double Zt = anElement->GetZ();  // number of protons
    608 
    609 
    610   return GetHadronNucleonXscPDG( aParticle, At, Zt );
     606                                                    const G4Element* anElement)
     607{
     608  G4int At = G4lrint(anElement->GetN());  // number of nucleons
     609  G4int Zt = G4lrint(anElement->GetZ());  // number of protons
     610
     611  return GetHadronNucleonXscPDG(aParticle, At, Zt);
    611612}
    612613
    … …  
    622623G4double
    623624G4GlauberGribovCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
    624                                                      G4double At,  G4double Zt )
     625                                                    G4int At, G4int Zt)
    625626{
    626627  G4double xsection;
    627628
    628   G4double Nt = At-Zt;              // number of neutrons
    629 
    630   if (Nt < 0.) Nt = 0.; 
    631 
     629  G4int Nt = At-Zt;              // number of neutrons
     630  if (Nt < 0) Nt = 0;
     631 
     632  G4double zz = Zt;
     633  G4double aa = At;
     634  G4double nn = Nt;
    632635
    633636  G4double targ_mass = G4ParticleTable::GetParticleTable()->
    634   GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
     637    GetIonTable()->GetIonMass(Zt, At);
    635638
    636639  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
    … …  
    656659  if(theParticle == theNeutron) // proton-neutron fit
    657660  {
    658     xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
     661    xsection = zz*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    659662                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    660     xsection  += Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
     663    xsection  += nn*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    661664                      + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
    662665  }
    … …  
    664667  {
    665668     
    666       xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
     669      xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    667670                          + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
    668671
    669       xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
     672      xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    670673                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    671674  }
    672675  else if(theParticle == theAProton)
    673676  {
    674     xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
     677    xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    675678                          + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
    676679
    677     xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
     680    xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    678681                          + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
    679682  }
    680683  else if(theParticle == thePiPlus)
    681684  {
    682     xsection  = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
     685    xsection  = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
    683686                          + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
    684687  }
    685688  else if(theParticle == thePiMinus)
    686689  {
    687     xsection  = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
     690    xsection  = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
    688691                          + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
    689692  }
    690693  else if(theParticle == theKPlus)
    691694  {
    692     xsection  = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
     695    xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
    693696                          + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
    694697
    695     xsection += Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
     698    xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
    696699                          + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
    697700  }
    698701  else if(theParticle == theKMinus)
    699702  {
    700     xsection  = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
     703    xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
    701704                          + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
    702705
    703     xsection += Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
     706    xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
    704707                          + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
    705708  }
    706709  else if(theParticle == theSMinus)
    707710  {
    708     xsection  = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
     711    xsection  = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
    709712                          - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
    710713  }
    711714  else if(theParticle == theGamma) // modify later on
    712715  {
    713     xsection  = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
     716    xsection  = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
    714717                          + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
    715718   
    … …  
    717720  else  // as proton ???
    718721  {
    719     xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
     722    xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    720723                          + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
    721724
    722     xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
     725    xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    723726                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    724727  }
    … …  
    735738G4double
    736739G4GlauberGribovCrossSection::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
    737                                                   const G4Element* anElement          )
    738 {
    739   G4double At = anElement->GetN();  // number of nucleons
    740   G4double Zt = anElement->GetZ();  // number of protons
    741 
    742 
    743   return GetHadronNucleonXscNS( aParticle, At, Zt );
     740                                                   const G4Element* anElement)
     741{
     742  G4int At = G4lrint(anElement->GetN());  // number of nucleons
     743  G4int Zt = G4lrint(anElement->GetZ());  // number of protons
     744
     745  return GetHadronNucleonXscNS(aParticle, At, Zt);
    744746}
    745747
    … …  
    754756G4double
    755757G4GlauberGribovCrossSection::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
    756                                                      G4double At,  G4double Zt )
     758                                                   G4int At, G4int Zt)
    757759{
    758760  G4double xsection(0), Delta, A0, B0;
    … …  
    760762  G4double hnXsc(0);
    761763
    762   G4double Nt = At-Zt;              // number of neutrons
    763 
    764   if (Nt < 0.) Nt = 0.; 
    765 
     764  G4int Nt = At-Zt;              // number of neutrons
     765  if (Nt < 0) Nt = 0; 
     766
     767  G4double aa = At;
     768  G4double zz = Zt;
     769  G4double nn = Nt;
    766770
    767771  G4double targ_mass = G4ParticleTable::GetParticleTable()->
    768   GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
     772  GetIonTable()->GetIonMass(Zt, At);
    769773
    770774  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
    … …  
    810814                     0.93827*0.93827,-0.165);        //  mb
    811815      }
    812       xsection *= Zt + Nt;
     816      xsection *= zz + nn;
    813817    }
    814818    else
    … …  
    846850                 (std::pow(proj_momentum,2.50)+0.95);
    847851      }
    848       xsection = hpXsc*Zt + hnXsc*Nt;
     852      xsection = hpXsc*zz + hnXsc*nn;
    849853    }
    850854  }
    … …  
    867871                     0.93827*0.93827,-0.165);        //  mb
    868872      }
    869       xsection *= Zt + Nt;
     873      xsection *= zz + nn;
    870874    }
    871875    else
    … …  
    903907                 (std::pow(proj_momentum,2.50)+0.95);
    904908      }
    905       xsection = hpXsc*Zt + hnXsc*Nt;
     909      xsection = hpXsc*zz + hnXsc*nn;
    906910      // xsection = hpXsc*(Zt + Nt);
    907911      // xsection = hnXsc*(Zt + Nt);
    … …  
    921925    if( proj_momentum <= 1.0 )
    922926    {
    923       xsection  = Zt*(65.55 + 53.84/(proj_momentum+1.e-6)  );
     927      xsection  = zz*(65.55 + 53.84/(proj_momentum+1.e-6)  );
    924928    }
    925929    else
    926930    {
    927       xsection  = Zt*( 41.1 + 77.2*std::pow( proj_momentum, -0.68)
     931      xsection  = zz*( 41.1 + 77.2*std::pow( proj_momentum, -0.68)
    928932                       + 0.293*logP*logP - 1.82*logP );
    929933    }
    930     if ( Nt > 0.) 
    931     {
    932       xsection += Nt*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP);
     934    if ( nn > 0.) 
     935    {
     936      xsection += nn*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP);
    933937    }
    934938    else // H
    … …  
    987991      hnXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
    988992    }
    989     xsection = hpXsc*Zt + hnXsc*Nt;
     993    xsection = hpXsc*zz + hnXsc*nn;
    990994  }
    991995  else if(theParticle == thePiMinus)
    … …  
    10391043      hpXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
    10401044    }
    1041     xsection = hpXsc*Zt + hnXsc*Nt;
     1045    xsection = hpXsc*zz + hnXsc*nn;
    10421046  }
    10431047  else if(theParticle == theKPlus)
    10441048  {
    1045     xsection  = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
     1049    xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
    10461050                          + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
    10471051
    1048     xsection += Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
     1052    xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
    10491053                          + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
    10501054  }
    10511055  else if(theParticle == theKMinus)
    10521056  {
    1053     xsection  = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
     1057    xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
    10541058                          + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
    10551059
    1056     xsection += Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
     1060    xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
    10571061                          + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
    10581062  }
    10591063  else if(theParticle == theSMinus)
    10601064  {
    1061     xsection  = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
     1065    xsection  = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
    10621066                          - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
    10631067  }
    10641068  else if(theParticle == theGamma) // modify later on
    10651069  {
    1066     xsection  = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
     1070    xsection  = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
    10671071                          + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
    10681072   
    … …  
    10701074  else  // as proton ???
    10711075  {
    1072     xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
     1076    xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    10731077                          + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
    10741078
    1075     xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
     1079    xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    10761080                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    10771081  }
    … …  
    10871091G4double
    10881092G4GlauberGribovCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
    1089                                                   const G4Element* anElement          )
    1090 {
    1091   G4double At = anElement->GetN();  // number of nucleons
    1092   G4double Zt = anElement->GetZ();  // number of protons
    1093 
    1094 
    1095   return GetHNinelasticXsc( aParticle, At, Zt );
     1093                                               const G4Element* anElement)
     1094{
     1095  G4int At = G4lrint(anElement->GetN());  // number of nucleons
     1096  G4int Zt = G4lrint(anElement->GetZ());  // number of protons
     1097
     1098  return GetHNinelasticXsc(aParticle, At, Zt);
    10961099}
    10971100
    … …  
    11021105G4double
    11031106G4GlauberGribovCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
    1104                                                      G4double At,  G4double Zt )
     1107                                                     G4int At,  G4int Zt)
    11051108{
    11061109  G4ParticleDefinition* hadron = aParticle->GetDefinition();
    1107   G4double sumInelastic, Nt = At - Zt;
    1108   if(Nt < 0.) Nt = 0.;
     1110  G4double sumInelastic;
     1111  G4int Nt = At - Zt;
     1112  if(Nt < 0) Nt = 0;
    11091113 
    11101114  if( hadron == theKPlus )
    … …  
    11161120    //sumInelastic  = Zt*GetHadronNucleonXscMK(aParticle, theProton);
    11171121    // sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);   
    1118     sumInelastic  = Zt*GetHadronNucleonXscNS(aParticle, 1.0, 1.0);
    1119     sumInelastic += Nt*GetHadronNucleonXscNS(aParticle, 1.0, 0.0);   
     1122    sumInelastic  = G4double(Zt)*GetHadronNucleonXscNS(aParticle, 1, 1);
     1123    sumInelastic += G4double(Nt)*GetHadronNucleonXscNS(aParticle, 1, 0);   
    11201124  }
    11211125  return sumInelastic;
    … …  
    11291133G4double
    11301134G4GlauberGribovCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
    1131                                                      G4double At,  G4double Zt )
     1135                                                 G4int At, G4int Zt)
    11321136{
    11331137  G4int PDGcode    = aParticle->GetDefinition()->GetPDGEncoding();
    … …  
    11471151  //G4cout<<"Plab = "<<Plab<<G4endl;
    11481152
    1149   G4double NumberOfTargetProtons  = Zt;
    1150   G4double NumberOfTargetNucleons = At;
     1153  G4double NumberOfTargetProtons = G4double(Zt);
     1154  G4double NumberOfTargetNucleons = G4double(At);
    11511155  G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
    11521156
    1153   if( NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.;
     1157  if(NumberOfTargetNeutrons < 0.0) NumberOfTargetNeutrons = 0.0;
    11541158
    11551159  G4double Xtotal, Xelastic, Xinelastic;
    … …  
    11691173                         0.169*sqrLogPlab - 1.85*LogPlab;
    11701174
    1171        Xtotal          = ( NumberOfTargetProtons * XtotPP +
    1172                            NumberOfTargetNeutrons * XtotPN  );
    1173 
    1174        Xelastic        = ( NumberOfTargetProtons  * XelPP +
    1175                            NumberOfTargetNeutrons * XelPN   );
     1175       Xtotal          = (NumberOfTargetProtons * XtotPP +
     1176                          NumberOfTargetNeutrons * XtotPN);
     1177
     1178       Xelastic        = (NumberOfTargetProtons * XelPP +
     1179                          NumberOfTargetNeutrons * XelPN);
    11761180  }
    11771181  else if( PDGcode ==  211 ) //------Projectile is PionPlus -------
    … …  
    13401344
    13411345G4double
    1342 G4GlauberGribovCrossSection::GetNucleusRadius( const G4DynamicParticle* ,
    1343                                                const G4Element* anElement)
    1344 {
    1345   G4double At       = anElement->GetN();
     1346G4GlauberGribovCrossSection::GetNucleusRadius(const G4DynamicParticle* ,
     1347                                              const G4Element* anElement)
     1348{
     1349  G4int At = G4lrint(anElement->GetN());
    13461350  G4double oneThird = 1.0/3.0;
    1347   G4double cubicrAt = std::pow (At, oneThird);
     1351  G4double cubicrAt = std::pow(G4double(At), oneThird);
    13481352
    13491353  G4double R;  // = fRadiusConst*cubicrAt;
    … …  
    13771381  G4double b3 = 4.;
    13781382
    1379   if (At > 20.)   // 20.
     1383  if (At > 20)   // 20.
    13801384  {
    13811385    R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
    13821386  }
    1383   else if (At > 3.5)
     1387  else if (At > 3)
    13841388  {
    13851389    R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
    … …  
    13971401
    13981402G4double
    1399 G4GlauberGribovCrossSection::GetNucleusRadius(G4double At)
     1403G4GlauberGribovCrossSection::GetNucleusRadius(G4int At)
    14001404{
    14011405  G4double oneThird = 1.0/3.0;
    1402   G4double cubicrAt = std::pow (At, oneThird);
     1406  G4double cubicrAt = std::pow(G4double(At), oneThird);
    14031407
    14041408  G4double R;  // = fRadiusConst*cubicrAt;
    … …  
    14241428  G4double tauA  = 20.;
    14251429
    1426   if (At > 20.)   // 20.
    1427   {
    1428     R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
     1430  if (At > 20)   // 20.
     1431  {
     1432    R *= ( 0.8 + 0.2*std::exp( -(G4double(At) - meanA)/tauA) );
    14291433  }
    14301434  else
    14311435  {
    1432     R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) );
     1436    R *= ( 1.0 + 0.1*( 1. - std::exp( (G4double(At) - meanA)/tauA) ) );
    14331437  }
    14341438

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

    r1337 r1340  
    2525//
    2626//
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929//
    … …  
    12301230      G4IsotopeVector* isoVector = element->GetIsotopeVector();
    12311231      G4double* abundVector = element->GetRelativeAbundanceVector();
    1232       G4double ZZ;
    1233       G4double AA;
     1232      G4int ZZ;
     1233      G4int AA;
    12341234
    12351235      for (G4int i = 0; i < nIso; i++) {
    1236         ZZ = G4double( (*isoVector)[i]->GetZ() );
    1237         AA = G4double( (*isoVector)[i]->GetN() );
     1236        ZZ = (*isoVector)[i]->GetZ();
     1237        AA = (*isoVector)[i]->GetN();
    12381238        CalcScatteringCrossSections(particle, ZZ, AA);
    12391239        cross_section += siginelastic*abundVector[i];
    … …  
    12411241      siginelastic = cross_section;
    12421242
    1243     } else {
    1244       CalcScatteringCrossSections(particle, element->GetZ(), element->GetN());
     1243    } else {
     1244      G4int ZZ = G4lrint(element->GetZ());
     1245      G4int AA = G4lrint(element->GetN());
     1246      CalcScatteringCrossSections(particle, ZZ, AA);
    12451247    }
    12461248  }
    … …  
    12491251
    12501252
    1251 G4double G4HadronCrossSections::GetInelasticCrossSection(
    1252                           const G4DynamicParticle* particle,
    1253                           G4double ZZ, G4double AA)
     1253G4double
     1254G4HadronCrossSections::GetInelasticCrossSection(const G4DynamicParticle* particle,
     1255                                                G4int ZZ, G4int AA)
    12541256{
    12551257  prevElement = 0; // force new cross section calculation for next call of
    … …  
    12741276      G4IsotopeVector* isoVector = element->GetIsotopeVector();
    12751277      G4double* abundVector = element->GetRelativeAbundanceVector();
    1276       G4double ZZ;
    1277       G4double AA;
     1278      G4int ZZ;
     1279      G4int AA;
    12781280
    12791281      for (G4int i = 0; i < nIso; i++) {
    1280         ZZ = G4double( (*isoVector)[i]->GetZ() );
    1281         AA = G4double( (*isoVector)[i]->GetN() );
     1282        ZZ = (*isoVector)[i]->GetZ();
     1283        AA = (*isoVector)[i]->GetN();
    12821284        CalcScatteringCrossSections(particle, ZZ, AA);
    12831285        cross_section += sigelastic*abundVector[i];
    … …  
    12851287      sigelastic = cross_section;
    12861288
    1287     } else {
    1288       CalcScatteringCrossSections(particle, element->GetZ(), element->GetN());
     1289    } else {
     1290      G4int ZZ = G4lrint(element->GetZ());
     1291      G4int AA = G4lrint(element->GetN());
     1292      CalcScatteringCrossSections(particle, ZZ, AA);
    12891293    }
    12901294  }
    … …  
    12931297
    12941298
    1295 G4double G4HadronCrossSections::GetElasticCrossSection(
    1296                           const G4DynamicParticle* particle,
    1297                           G4double ZZ, G4double AA)
     1299G4double
     1300G4HadronCrossSections::GetElasticCrossSection(const G4DynamicParticle* particle,
     1301                                              G4int ZZ, G4int AA)
    12981302{
    12991303  prevElement = 0; // force new cross section calculation for next call of
    … …  
    13121316G4HadronCrossSections::CalcScatteringCrossSections(
    13131317                          const G4DynamicParticle* aParticle,
    1314                           G4double ZZ, G4double AA)
     1318                          G4int ZZ, G4int AA)
    13151319{
    13161320   G4double sigel, sigin, sigtot;
    … …  
    15831587G4HadronCrossSections::GetCaptureCrossSection(
    15841588                          const G4DynamicParticle* aParticle,
    1585                           G4double ZZ, G4double /*AA*/)
     1589                          G4int ZZ, G4int /*AA*/)
    15861590{
    15871591   if (GetParticleCode(aParticle) != 16)  { return 0.; }
    … …  
    15961600   }
    15971601
    1598    G4int izno = static_cast<G4int> (ZZ + 0.01);
     1602   G4int izno = ZZ;
    15991603   if (izno > 100) izno = 100;      // Not in GHESIG
    16001604   izno = izno - 1;      // For array indexing
    … …  
    16091613G4HadronCrossSections::GetFissionCrossSection(
    16101614                          const G4DynamicParticle* aParticle,
    1611                           G4double ZZ, G4double AA)
     1615                          G4int ZZ, G4int AA)
    16121616{
    1613    if (AA < 230.) return 0;
     1617   if (AA < 230) return 0;
    16141618
    16151619   G4double ek = aParticle->GetKineticEnergy()/GeV;     
    1616 
    1617    //   G4int i = NFISS;
    1618    //   for (G4int ii = 1; i <= NFISS; i++) {
    1619    //      if (ek < ekfiss[ii - 1]) {
    1620    //         i = ii;
    1621    //         break;
    1622    //      }
    1623    //   }
    1624    //   i = i - 1;      // For array indexing
    16251620
    16261621   G4int ie1 = 0;
    … …  
    16381633   G4int j = 4;
    16391634   if (ek <= 0.01) {
    1640       if (ZZ == 92. && std::abs(AA - 233.) < 0.5) j = 1;
    1641       else if (ZZ == 92. && std::abs(AA - 235.) < 0.5) j = 2;
    1642       else if (ZZ == 94. && std::abs(AA - 239.) < 0.5) j = 3;
     1635      if (ZZ == 92 && AA == 233) j = 1;
     1636      else if (ZZ == 92 && AA == 235) j = 2;
     1637      else if (ZZ == 94 && AA == 239) j = 3;
    16431638   }
    16441639   G4double z43ba;
    16451640   if (j == 4) {
    1646       z43ba = std::pow(ZZ, 4./3.)/AA;
     1641      z43ba = std::pow(G4double(ZZ), 4./3.)/G4double(AA);
    16471642      z43ba = std::max(-67. + 38.7*z43ba, 0.);
    16481643   }

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

    r1228 r1340  
    3333#include "G4IonTable.hh"
    3434#include "G4ParticleDefinition.hh"
    35 
    36 //////////////////////////////////////////////////////////////////////////////////////
    37 //
    38 //
     35#include "G4HadTmpUtil.hh"
    3936
    4037
    … …  
    8077
    8178G4HadronNucleonXsc::~G4HadronNucleonXsc()
    82 {
    83 }
     79{}
    8480
    8581
    … …  
    9187G4bool
    9288G4HadronNucleonXsc::IsApplicable(const G4DynamicParticle* aDP,
    93                                           const G4Element* anElement)
     89                                 const G4Element* anElement)
    9490{
    95   return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN());
     91  G4int Z = G4lrint(anElement->GetZ());
     92  G4int A = G4lrint(anElement->GetN());
     93  return IsIsoApplicable(aDP, Z, A);
    9694}
    9795
    9896////////////////////////////////////////////////////////////////////////////////////////
    9997//
    100 //
    10198
    10299G4bool
    103 G4HadronNucleonXsc::IsZAApplicable(const G4DynamicParticle* aDP,
    104                                             G4double Z, G4double)
     100G4HadronNucleonXsc::IsIsoApplicable(const G4DynamicParticle* aDP,
     101                                    G4int Z, G4int)
    105102{
    106   G4bool applicable      = false;
     103  G4bool applicable = false;
    107104  // G4int baryonNumber     = aDP->GetDefinition()->GetBaryonNumber();
    108105  G4double kineticEnergy = aDP->GetKineticEnergy();
    … …  
    111108 
    112109  if ( ( kineticEnergy  >= fLowerLimit &&
    113          Z > 1.5 &&      // >=  He
     110         Z > 1 &&      // >=  He
    114111       ( theParticle == theAProton   ||
    115112         theParticle == theGamma     ||
    … …  
    119116
    120117       ( kineticEnergy  >= 0.1*fLowerLimit &&
    121          Z > 1.5 &&      // >=  He
     118         Z > 1 &&      // >=  He
    122119       ( theParticle == theProton    ||
    123120         theParticle == theNeutron   ||   
    … …  
    127124  return applicable;
    128125}
    129 
    130 
    131126
    132127
    … …  
    210205
    211206
    212 
    213 
    214 
    215207/////////////////////////////////////////////////////////////////////////////////////
    216208//
    … …  
    242234  G4double B    = 0.308;
    243235
    244 
    245236  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
    246237
    … …  
    249240  G4bool neutron = (nucleon == theNeutron);
    250241 
    251 
    252242  if(theParticle == theNeutron) // proton-neutron fit
    253243  {
    … …  
    358348  return xsection;
    359349}
    360 
    361350
    362351
    … …  
    758747  G4bool neutron = (nucleon == theNeutron);
    759748
    760 
    761 
    762749   
    763750  if( absPDGcode > 1000 && pORn )  //------Projectile is baryon -

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

    r1055 r1340  
    3232#include "G4ParticleTable.hh"
    3333#include "G4IonTable.hh"
     34#include "G4HadTmpUtil.hh"
    3435
    3536G4double G4IonsKoxCrossSection::
    36 GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double ZZ,
    37                      G4double AA, G4double /*temperature*/)
     37GetZandACrossSection(const G4DynamicParticle* aParticle, G4int ZZ,
     38                     G4int AA, G4double /*temperature*/)
    3839{
    3940   G4double xsection = 0.0;
    4041
    4142   G4int Ap = aParticle->GetDefinition()->GetBaryonNumber();
    42    G4int Zp = int ( aParticle->GetDefinition()->GetPDGCharge() / eplus + 0.5);
     43   G4int Zp = G4int(aParticle->GetDefinition()->GetPDGCharge() / eplus + 0.5);
    4344   G4double ke_per_N = aParticle->GetKineticEnergy() / Ap;
    4445
    … …  
    4647   // if (  ke_per_N < lowerLimit ) return xsection;
    4748
    48    G4int At = int (AA + 0.5);
    49    G4int Zt = int (ZZ + 0.5 ); 
     49   G4int At = AA;
     50   G4int Zt = ZZ; 
    5051
    5152   G4double one_third = 1.0 / 3.0;
    … …  
    5455   G4double cubicrAp = std::pow ( G4double(Ap) , G4double(one_third) ); 
    5556
     57   // rc divide fermi
     58   G4double Bc = Zt * Zp / ( (rc/fermi) * (cubicrAp+cubicrAt) );
    5659
    57    G4double Bc = Zt * Zp / ( ( rc / fermi ) * (  cubicrAp + cubicrAt ) );   // rc divide fermi
    58    G4double targ_mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( Zt , At );
    59    G4double proj_mass = aParticle->GetMass(); 
     60   G4double targ_mass =
     61     G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Zt, At);
     62   G4double proj_mass = aParticle->GetMass();
    6063   G4double proj_momentum = aParticle->GetMomentum().mag();
    6164
    … …  
    6972
    7073   G4double a = 1.85;
    71    G4double Rsurf = r0 * ( a * cubicrAp * cubicrAt / ( cubicrAp + cubicrAt ) - c);
     74   G4double Rsurf = r0 * (a*cubicrAp * cubicrAt/(cubicrAp + cubicrAt) - c);
    7275   G4double D = 5.0 * ( At - 2 * Zt ) * Zp / ( Ap * At );
    7376   Rsurf = Rsurf + D * fermi;  // multiply D by fermi
    7477
    7578   G4double Rint = Rvol + Rsurf;
    76 
    7779   xsection = pi * Rint * Rint * ( 1 - Bc / ( Ecm / MeV ) );
    7880 
    7981   return xsection;
    8082}
     83
    8184
    8285G4double G4IonsKoxCrossSection::
    … …  
    9194    G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
    9295    G4double* abundVector = anElement->GetRelativeAbundanceVector();
    93     G4double ZZ;
    94     G4double AA;
     96    G4int ZZ;
     97    G4int AA;
    9598     
    9699    for (G4int i = 0; i < nIso; i++) {
    97       ZZ = G4double( (*isoVector)[i]->GetZ() );
    98       AA = G4double( (*isoVector)[i]->GetN() );
    99       sig = GetIsoZACrossSection(aParticle, ZZ, AA, temperature);
     100      ZZ = (*isoVector)[i]->GetZ();
     101      AA = (*isoVector)[i]->GetN();
     102      sig = GetZandACrossSection(aParticle, ZZ, AA, temperature);
    100103      xsection += sig*abundVector[i];
    101104    }
    102105   
    103106  } else {
    104     xsection =
    105       GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(),
    106                           temperature);
     107    G4int ZZ = G4lrint(anElement->GetZ());
     108    G4int AA = G4lrint(anElement->GetN());
     109    xsection = GetIsoZACrossSection(aParticle, ZZ, AA, temperature);
    107110  }
    108111   
    … …  
    111114
    112115
    113 G4double G4IonsKoxCrossSection::calEcm ( G4double mp , G4double mt , G4double Plab )
     116G4double
     117G4IonsKoxCrossSection::calEcm(G4double mp, G4double mt, G4double Plab)
    114118{
    115119   G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
    … …  
    121125
    122126
    123 G4double G4IonsKoxCrossSection::calCeValue( const G4double ke )
     127G4double G4IonsKoxCrossSection::calCeValue(const G4double ke)
    124128{
    125129   // Calculate c value
    … …  
    131135   G4double Ce;
    132136   G4double log10_ke = std::log10 ( ke );   
    133    if ( log10_ke > 1.5 )
     137   if (log10_ke > 1.5)
    134138   {
    135139      Ce = - 10.0 / std::pow ( G4double(log10_ke) , G4double(5) ) + 2.0;
    … …  
    137141   else
    138142   {
    139       Ce = ( - 10.0 / std::pow ( G4double(1.5) , G4double(5) ) + 2.0 ) / std::pow ( G4double(1.5) , G4double(3) ) * std::pow ( G4double(log10_ke) , G4double(3) );
     143      Ce = (-10.0/std::pow(G4double(1.5), G4double(5) ) + 2.0) /
     144           std::pow(G4double(1.5), G4double(3)) * std::pow(G4double(log10_ke), G4double(3) );
    140145
    141146   }

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

    r1055 r1340  
    3333#include "G4ParticleTable.hh"
    3434#include "G4IonTable.hh"
     35#include "G4HadTmpUtil.hh"
    3536
    3637
    3738G4double G4IonsShenCrossSection::
    38 GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double ZZ,
    39                 G4double AA, G4double /*temperature*/)
     39GetZandACrossSection(const G4DynamicParticle* aParticle, G4int ZZ,
     40                     G4int AA, G4double /*temperature*/)
    4041{
    4142   G4double xsection = 0.0;
    4243
    4344   G4int Ap = aParticle->GetDefinition()->GetBaryonNumber();
    44    G4int Zp = int ( aParticle->GetDefinition()->GetPDGCharge() / eplus + 0.5 );
     45   G4int Zp = G4int(aParticle->GetDefinition()->GetPDGCharge()/eplus + 0.5 );
    4546   G4double ke_per_N = aParticle->GetKineticEnergy() / Ap;
    4647   if ( ke_per_N > 10*GeV ) ke_per_N = 10*GeV;
    … …  
    4950   // if (  ke_per_N < lowerLimit ) return xsection;
    5051
    51    G4int At = G4int(AA);
    52    G4int Zt = G4int(ZZ);
     52   G4int At = AA;
     53   G4int Zt = ZZ;
    5354 
    5455   G4double one_third = 1.0 / 3.0;
    … …  
    6263   G4double r = Rt + Rp + 3.2;   // in fm
    6364   G4double b = 1.0;   // in MeV/fm
    64    G4double targ_mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( Zt , At );
     65   G4double targ_mass =
     66     G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Zt, At);
    6567   G4double proj_mass = aParticle->GetMass();
    6668   G4double proj_momentum = aParticle->GetMomentum().mag();
    6769
    68    G4double Ecm = calEcmValue ( proj_mass , targ_mass , proj_momentum );
     70   G4double Ecm = calEcmValue (proj_mass, targ_mass, proj_momentum);
    6971
    7072   G4double B = 1.44 * Zt * Zp / r - b * Rt * Rp / ( Rt + Rp );
    … …  
    7476   G4double c = calCeValue ( ke_per_N / MeV  ); 
    7577
    76    G4double R1 = r0 * ( cubicrAt + cubicrAp + 1.85 * cubicrAt * cubicrAp / ( cubicrAt + cubicrAp ) - c);
     78   G4double R1 = r0 * (cubicrAt + cubicrAp + 1.85*cubicrAt*cubicrAp/(cubicrAt + cubicrAp) - c);
    7779
    7880   G4double R2 = 1.0 * ( At - 2 * Zt ) * Zp / ( Ap * At );
    7981
    8082
    81    G4double R3 = 0.176 / std::pow ( G4double(Ecm) , G4double(one_third) ) * cubicrAt * cubicrAp / ( cubicrAt + cubicrAp );
     83   G4double R3 = 0.176 / std::pow(G4double(Ecm), G4double(one_third)) * cubicrAt * cubicrAp /(cubicrAt + cubicrAp);
    8284
    8385   G4double R = R1 + R2 + R3;
    8486
    8587   xsection = 10 * pi * R * R * ( 1 - B / Ecm );   
    86    xsection = xsection * millibarn;   // mulitply xsection by millibarn   
    87  
     88   xsection = xsection * millibarn;   // mulitply xsection by millibarn
     89
    8890   return xsection;
    8991}
    … …  
    101103    G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
    102104    G4double* abundVector = anElement->GetRelativeAbundanceVector();
    103     G4double ZZ;
    104     G4double AA;
     105    G4int ZZ;
     106    G4int AA;
    105107   
    106108    for (G4int i = 0; i < nIso; i++) {
    107       ZZ = G4double( (*isoVector)[i]->GetZ() );
    108       AA = G4double( (*isoVector)[i]->GetN() );
     109      ZZ = (*isoVector)[i]->GetZ();
     110      AA = (*isoVector)[i]->GetN();
    109111      sig = GetIsoZACrossSection(aParticle, ZZ, AA, temperature);
    110112      xsection += sig*abundVector[i];
    … …  
    112114 
    113115  } else {
    114     xsection =
    115       GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(),
    116                           temperature);
     116    G4int ZZ = G4lrint(anElement->GetZ());
     117    G4int AA = G4lrint(anElement->GetN());
     118    xsection = GetIsoZACrossSection(aParticle, ZZ, AA, temperature);
    117119  }
    118120 
    … …  
    121123
    122124
    123 G4double G4IonsShenCrossSection::calEcmValue( const G4double mp , const G4double mt , const G4double Plab )
     125G4double
     126G4IonsShenCrossSection::calEcmValue(const G4double mp, const G4double mt,
     127                                    const G4double Plab)
    124128{
    125129   G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
    … …  
    131135
    132136
    133 G4double G4IonsShenCrossSection::calCeValue( const G4double ke )
     137G4double G4IonsShenCrossSection::calCeValue(const G4double ke)
    134138{
    135    // Calculate c value
    136    // This value is indepenent from projectile and target particle
    137    // ke is projectile kinetic energy per nucleon in the Lab system with MeV unit
    138    // fitting function is made by T. Koi
    139    // There are no data below 30 MeV/n in Kox et al.,
     139  // Calculate c value
     140  // This value is indepenent from projectile and target particle
     141  // ke is projectile kinetic energy per nucleon in the Lab system
     142  // with MeV unit
     143  // fitting function is made by T. Koi
     144  // There are no data below 30 MeV/n in Kox et al.,
    140145
    141146   G4double Ce;
    142147   G4double log10_ke = std::log10 ( ke );   
    143    if ( log10_ke > 1.5 )
     148   if (log10_ke > 1.5)
    144149   {
    145       Ce = - 10.0 / std::pow ( G4double(log10_ke) , G4double(5) ) + 2.0;
     150     Ce = -10.0/std::pow(G4double(log10_ke), G4double(5)) + 2.0;
    146151   }
    147152   else
    148153   {
    149       Ce = ( - 10.0 / std::pow ( G4double(1.5) , G4double(5) ) + 2.0 ) / std::pow ( G4double(1.5) , G4double(3) ) * std::pow ( G4double(log10_ke) , G4double(3) );
     154     Ce = (-10.0/std::pow(G4double(1.5), G4double(5) ) + 2.0) /
     155         std::pow(G4double(1.5) , G4double(3)) * std::pow(G4double(log10_ke), G4double(3));
    150156   }
    151157   return Ce;

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

    r819 r1340  
    3131#include "G4ParticleTable.hh"
    3232#include "G4IonTable.hh"
     33#include "G4HadTmpUtil.hh"
    3334
    3435G4double G4IonsSihverCrossSection::
    35 GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double /*ZZ*/,
    36                 G4double AA, G4double /*aTemperature*/)
     36GetZandACrossSection(const G4DynamicParticle* aParticle,
     37                     G4int /*ZZ*/, G4int AA, G4double /*aTemperature*/)
    3738{
    3839   G4double xsection = 0.0;
    39 
    40    G4int At = G4int(AA);
    41    //G4int Zt = G4int(ZZ);  // not used
     40   G4int At = AA;
    4241
    4342   G4int Ap = aParticle->GetDefinition()->GetBaryonNumber();
    44    //Zp = aParticle->GetDefinition()->GetPDGCharge();   // not used 
    4543 
    4644   G4double one_third = 1.0 / 3.0;
    … …  
    4947   G4double cubicrAp = std::pow ( G4double(Ap) , G4double(one_third) ); 
    5048
    51    G4double b0 = 1.581 - 0.876 * ( 1.0 / cubicrAp + 1.0 / cubicrAt );
     49   G4double b0 = 1.581 - 0.876 * (1.0/cubicrAp + 1.0/cubicrAt);
    5250
    5351   xsection = pi * square_r0
    54             * std::pow ( G4double(cubicrAp + cubicrAt - b0 * (  1.0 / cubicrAp + 1.0 / cubicrAt ) ), G4double(2) );
     52            * std::pow(G4double(cubicrAp + cubicrAt - b0 * (1.0/cubicrAp + 1.0/cubicrAt)), G4double(2));
    5553 
    5654   return xsection;
    … …  
    6967    G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
    7068    G4double* abundVector = anElement->GetRelativeAbundanceVector();
    71     G4double ZZ;
    72     G4double AA;
     69    G4int ZZ;
     70    G4int AA;
    7371   
    7472    for (G4int i = 0; i < nIso; i++) {
    75       ZZ = G4double( (*isoVector)[i]->GetZ() );
    76       AA = G4double( (*isoVector)[i]->GetN() );
     73      ZZ = (*isoVector)[i]->GetZ();
     74      AA = (*isoVector)[i]->GetN();
    7775      sig = GetIsoZACrossSection(aParticle, ZZ, AA, temperature);
    7876      xsection += sig*abundVector[i];
    … …  
    8078 
    8179  } else {
    82     xsection =
    83       GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(),
    84                           temperature);
     80    G4int ZZ = G4lrint(anElement->GetZ());
     81    G4int AA = G4lrint(anElement->GetN());
     82    xsection = GetIsoZACrossSection(aParticle, ZZ, AA, temperature);
    8583  }
    8684   

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4NeutronCaptureXS.cc,v 1.3 2010/06/03 11:50:21 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4NeutronCaptureXS.cc,v 1.5 2010/10/15 22:36:13 dennis Exp $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4848#include <fstream>
    4949#include <sstream>
     50
    5051using namespace std;
    5152
    5253G4NeutronCaptureXS::G4NeutronCaptureXS()
     54  : emax(20*MeV),maxZ(92)
    5355{
    54   verboseLevel = 0;
    55   emax = 20.*MeV;
    56   G4cout  << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise " << G4endl;
    57   for(G4int i=0; i<93; ++i) {
    58     data[i] = 0;
     56  //  verboseLevel = 0;
     57  if(verboseLevel > 0){
     58    G4cout  << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise for Z < "
     59            << maxZ + 1 << G4endl;
    5960  }
     61  data.resize(maxZ+1, 0);
    6062  isInitialized = false;
    6163}
    … …  
    6365G4NeutronCaptureXS::~G4NeutronCaptureXS()
    6466{
    65   for(G4int i=0; i<92; ++i) {
     67  for(G4int i=0; i<=maxZ; ++i) {
    6668    delete data[i];
    6769  }
    … …  
    7678
    7779G4bool
    78 G4NeutronCaptureXS::IsZAApplicable(const G4DynamicParticle*,
    79                                       G4double /*ZZ*/, G4double /*AA*/)
     80G4NeutronCaptureXS::IsIsoApplicable(const G4DynamicParticle*,
     81                                    G4int /*ZZ*/, G4int /*AA*/)
    8082{
    8183  return false;
    … …  
    9092  G4double xs = 0.0;
    9193  G4double ekin = aParticle->GetKineticEnergy();
    92   if(ekin > emax) return xs;
     94  if(ekin > emax) { return xs; }
    9395
    9496  G4int Z = G4int(elm->GetZ());
    … …  
    99101    Initialise(Z);
    100102    pv = data[Z];
    101     if(!pv) return xs;
     103    if(!pv) { return xs; }
    102104  }
    103105
    … …  
    118120G4NeutronCaptureXS::BuildPhysicsTable(const G4ParticleDefinition& p)
    119121{
    120   G4cout << "G4NeutronCaptureXS::BuildPhysicsTable: " << G4endl;
    121   G4cout << p.GetParticleName() << G4endl;
    122   if(p.GetParticleName() != "neutron") {
    123     return;
     122  if(verboseLevel > 0){
     123    G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for "
     124           << p.GetParticleName() << G4endl;
    124125  }
    125   if(isInitialized) return;
     126  if(isInitialized || p.GetParticleName() != "neutron") { return; }
    126127  isInitialized = true;
    127128
    … …  
    140141    for(size_t i=0; i<numOfElm; ++i) {
    141142      G4int Z = G4int(((*theElmTable)[i])->GetZ());
    142       if(Z < 1) Z = 1;
    143       else if(Z > 92) Z = 92;
     143      if(Z < 1)         { Z = 1; }
     144      else if(Z > maxZ) { Z = maxZ; }
    144145      //G4cout << "Z= " << Z << G4endl;
    145146      // Initialisation
    … …  
    151152void
    152153G4NeutronCaptureXS::DumpPhysicsTable(const G4ParticleDefinition&)
    153 {
    154 }
     154{}
     155
    155156
    156157void
    157158G4NeutronCaptureXS::Initialise(G4int Z, const char* p)
    158159{
    159   if(data[Z]) return;
     160  if(data[Z]) { return; }
    160161  const char* path = p;
    161162  if(!p) {

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4NeutronElasticXS.cc,v 1.4 2010/06/03 11:50:21 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4NeutronElasticXS.cc,v 1.6 2010/10/15 22:35:58 dennis Exp $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4747#include "G4PhysicsVector.hh"
    4848#include "G4GlauberGribovCrossSection.hh"
     49#include "G4HadronNucleonXsc.hh"
    4950#include "G4NistManager.hh"
     51#include "G4Proton.hh"
    5052
    5153#include <iostream>
    5254#include <fstream>
    5355#include <sstream>
     56
    5457using namespace std;
    5558
    5659G4NeutronElasticXS::G4NeutronElasticXS()
    57 {
    58   verboseLevel = 0;
    59   G4cout  << "G4NeutronElasticXS::G4NeutronElasticXS: Initialise " << G4endl;
    60   for(G4int i=0; i<93; ++i) {
    61     data[i] = 0;
    62     coeff[i]= 1.0;
    63   }
     60  :  proton(G4Proton::Proton()), maxZ(92)
     61{
     62  //  verboseLevel = 0;
     63  if(verboseLevel > 0){
     64    G4cout  << "G4NeutronElasticXS::G4NeutronElasticXS Initialise for Z < "
     65            << maxZ + 1 << G4endl;
     66  }
     67  data.resize(maxZ+1, 0);
     68  coeff.resize(maxZ+1, 1.0);
    6469  ggXsection = new G4GlauberGribovCrossSection();
     70  fNucleon = new G4HadronNucleonXsc();
    6571  isInitialized = false;
    6672}
    … …  
    6874G4NeutronElasticXS::~G4NeutronElasticXS()
    6975{
    70   for(G4int i=0; i<92; ++i) {
     76  for(G4int i=0; i<=maxZ; ++i) {
    7177    delete data[i];
    7278  }
    … …  
    7581G4bool
    7682G4NeutronElasticXS::IsApplicable(const G4DynamicParticle*,
    77                                    const G4Element*)
     83                                 const G4Element*)
    7884{
    7985  return true;
    … …  
    8187
    8288G4bool
    83 G4NeutronElasticXS::IsZAApplicable(const G4DynamicParticle*,
    84                                      G4double /*ZZ*/, G4double /*AA*/)
     89G4NeutronElasticXS::IsIsoApplicable(const G4DynamicParticle*,
     90                                    G4int /*ZZ*/, G4int /*AA*/)
    8591{
    8692  return false;
    8793}
    88 
    8994
    9095G4double
    9196G4NeutronElasticXS::GetCrossSection(const G4DynamicParticle* aParticle,
    92                                       const G4Element* elm,
    93                                       G4double)
     97                                    const G4Element* elm, G4double)
    9498{
    9599  G4double xs = 0.0;
    … …  
    97101
    98102  G4int Z = G4int(elm->GetZ());
     103  if(Z < 1 || Z > maxZ) { return xs; }
    99104  G4PhysicsVector* pv = data[Z];
    100105  //  G4cout  << "G4NeutronElasticXS::GetCrossSection e= " << ekin << " Z= " << Z << G4endl;
    … …  
    104109    Initialise(Z);
    105110    pv = data[Z];
    106     if(!pv) return xs;
     111    if(!pv) { return xs; }
    107112  }
    108113
    109114  G4double e1 = pv->Energy(0);
    110   if(ekin <= e1) return (*pv)[0];
     115  if(ekin <= e1) { return (*pv)[0]; }
    111116
    112117  G4int n = pv->GetVectorLength() - 1;
    113118  G4double e2 = pv->Energy(n);
     119
    114120  if(ekin <= e2) {
    115121    xs = pv->Value(ekin);
     122  } else if(1 == Z) {
     123    fNucleon->GetHadronNucleonXscPDG(aParticle, proton);
     124    xs = coeff[1]*fNucleon->GetElasticHadronNucleonXsc();
    116125  } else {         
    117126    ggXsection->GetCrossSection(aParticle, elm);
    … …  
    129138G4NeutronElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
    130139{
    131   G4cout << "G4NeutronElasticXS::BuildPhysicsTable: " << G4endl;
    132   G4cout << p.GetParticleName() << G4endl;
    133   if(p.GetParticleName() != "neutron") {
    134     return;
    135   }
    136   if(isInitialized) return;
     140  if(verboseLevel > 0){
     141    G4cout << "G4NeutronElasticXS::BuildPhysicsTable for "
     142           << p.GetParticleName() << G4endl;
     143  }
     144  if(isInitialized || p.GetParticleName() != "neutron") { return; }
    137145  isInitialized = true;
    138 
    139146
    140147  // check environment variable
    … …  
    154161    for(size_t i=0; i<numOfElm; ++i) {
    155162      G4int Z = G4int(((*theElmTable)[i])->GetZ());
    156       if(Z < 1) Z = 1;
    157       else if(Z > 92) Z = 92;
     163      if(Z < 1)         { Z = 1; }
     164      else if(Z > maxZ) { Z = maxZ; }
    158165      //G4cout << "Z= " << Z << G4endl;
    159166      // Initialisation
    … …  
    173180                               const char* p)
    174181{
    175   if(data[Z]) return;
     182  if(data[Z]) { return; }
    176183  const char* path = p;
    177184  if(!p) {
    178   // check environment variable
    179   // Build the complete string identifying the file with the data set
     185    // check environment variable
     186    // Build the complete string identifying the file with the data set
    180187    path = getenv("G4NEUTRONXSDATA");
    181188    if (!path) {
    … …  
    215222    G4double sig1 = (*data[Z])[n];
    216223    dynParticle->SetKineticEnergy(emax);
    217     ggXsection->GetCrossSection(dynParticle, Elem);
    218     G4double sig2 = ggXsection->GetElasticGlauberGribovXsc();
     224    G4double sig2 = 0.0;
     225    if(1 == Z) {
     226      fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
     227      sig2 = fNucleon->GetElasticHadronNucleonXsc();
     228    } else {
     229      ggXsection->GetCrossSection(dynParticle, Elem);
     230      sig2 = ggXsection->GetElasticGlauberGribovXsc();
     231    }
    219232    if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
    220233  }

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

    r1228 r1340  
    3030
    3131#include "G4NeutronInelasticCrossSection.hh"
     32#include "G4HadTmpUtil.hh"
    3233#include "globals.hh"
    3334
    … …  
    4445    G4IsotopeVector* isoVector = anEle->GetIsotopeVector();
    4546    G4double* abundVector = anEle->GetRelativeAbundanceVector();
    46     G4double ZZ;
    47     G4double AA;
     47    G4int ZZ;
     48    G4int AA;
    4849
    4950    for (G4int i = 0; i < nIso; i++) {
    50       ZZ = G4double( (*isoVector)[i]->GetZ() );
    51       AA = G4double( (*isoVector)[i]->GetN() );
     51      ZZ = (*isoVector)[i]->GetZ();
     52      AA = (*isoVector)[i]->GetN();
    5253      psig = GetCrossSection(KE, AA, ZZ);
    5354      cross_section += psig*abundVector[i];
    … …  
    5556 
    5657  } else {
    57     cross_section = GetCrossSection(KE, anEle->GetN(), anEle->GetZ());
     58    G4int ZZ = G4lrint(anEle->GetZ());
     59    G4int AA = G4lrint(anEle->GetN());
     60    cross_section = GetCrossSection(KE, AA, ZZ);
    5861  }
    5962
    … …  
    6366   
    6467G4double G4NeutronInelasticCrossSection::
    65 GetCrossSection(G4double anEnergy, G4double atomicNumber, G4double nOfProtons)
     68GetCrossSection(G4double anEnergy, G4int AA, G4int ZZ)
    6669{
     70  G4double atomicNumber = G4double(AA);
     71  G4double nOfProtons = G4double(ZZ);
     72
    6773  if (anEnergy > 19.9*GeV )
    6874  { // constant cross section above ~20GeV.
    69     return  GetCrossSection(19.8*GeV,atomicNumber,nOfProtons);
     75    return  GetCrossSection(19.8*GeV, AA, ZZ);
    7076  }
    7177  G4double kineticEnergy = std::log10(DBL_MIN/MeV);

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4NeutronInelasticXS.cc,v 1.4 2010/06/03 11:50:21 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4NeutronInelasticXS.cc,v 1.6 2010/10/15 22:35:32 dennis Exp $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4747#include "G4PhysicsVector.hh"
    4848#include "G4GlauberGribovCrossSection.hh"
     49#include "G4HadronNucleonXsc.hh"
    4950#include "G4NistManager.hh"
     51#include "G4Proton.hh"
    5052
    5153#include <iostream>
    5254#include <fstream>
    5355#include <sstream>
     56
    5457using namespace std;
    5558
    5659G4NeutronInelasticXS::G4NeutronInelasticXS()
    57 {
    58   verboseLevel = 0;
    59   G4cout  << "G4NeutronInelasticXS::G4NeutronInelasticXS: Initialise " << G4endl;
    60   for(G4int i=0; i<93; ++i) {
    61     data[i] = 0;
    62     coeff[i]= 1.0;
    63   }
     60  : proton(G4Proton::Proton()), maxZ(92)
     61{
     62  //  verboseLevel = 0;
     63  if(verboseLevel > 0){
     64    G4cout  << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
     65            << maxZ + 1 << G4endl;
     66  }
     67  data.resize(maxZ+1, 0);
     68  coeff.resize(maxZ+1, 1.0);
    6469  ggXsection = new G4GlauberGribovCrossSection();
     70  fNucleon = new G4HadronNucleonXsc();
    6571  isInitialized = false;
    6672}
    … …  
    6874G4NeutronInelasticXS::~G4NeutronInelasticXS()
    6975{
    70   for(G4int i=0; i<92; ++i) {
     76  for(G4int i=0; i<=maxZ; ++i) {
    7177    delete data[i];
    7278  }
    … …  
    8187
    8288G4bool
    83 G4NeutronInelasticXS::IsZAApplicable(const G4DynamicParticle*,
    84                                      G4double /*ZZ*/, G4double /*AA*/)
     89G4NeutronInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
     90                                      G4int /*ZZ*/, G4int /*AA*/)
    8591{
    8692  return false;
    … …  
    97103
    98104  G4int Z = G4int(elm->GetZ());
     105  if(Z < 1 || Z > maxZ) { return xs; }
    99106  G4PhysicsVector* pv = data[Z];
    100107  //  G4cout  << "G4NeutronInelasticXS::GetCrossSection e= " << ekin << " Z= " << Z << G4endl;
    … …  
    104111    Initialise(Z);
    105112    pv = data[Z];
    106     if(!pv) return xs;
     113    if(!pv) { return xs; }
    107114  }
    108115
    109116  G4double e1 = pv->Energy(0);
    110   if(ekin <= e1) return xs;
     117  if(ekin <= e1) { return xs; }
    111118
    112119  G4int n = pv->GetVectorLength() - 1;
    … …  
    114121  if(ekin <= e2) {
    115122    xs = pv->Value(ekin);
     123  } else if(1 == Z) {
     124    fNucleon->GetHadronNucleonXscPDG(aParticle, proton);
     125    xs = coeff[1]*fNucleon->GetInelasticHadronNucleonXsc();
    116126  } else {         
    117127    ggXsection->GetCrossSection(aParticle, elm);
    … …  
    119129  }
    120130
    121   if(verboseLevel > 0){
     131  if(verboseLevel > 0) {
    122132    G4cout  << "ekin= " << ekin << ",  XSinel= " << xs << G4endl;
    123133  }
    … …  
    129139G4NeutronInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
    130140{
    131   G4cout << "G4NeutronInelasticXS::BuildPhysicsTable: " << G4endl;
    132   G4cout << p.GetParticleName() << G4endl;
    133   if(p.GetParticleName() != "neutron") {
    134     return;
    135   }
    136   if(isInitialized) return;
     141  if(verboseLevel > 0){
     142    G4cout << "G4NeutronInelasticXS::BuildPhysicsTable for "
     143           << p.GetParticleName() << G4endl;
     144  }
     145  if(isInitialized || p.GetParticleName() != "neutron") { return; }
    137146  isInitialized = true;
    138 
    139147
    140148  // check environment variable
    … …  
    154162    for(size_t i=0; i<numOfElm; ++i) {
    155163      G4int Z = G4int(((*theElmTable)[i])->GetZ());
    156       if(Z < 1) Z = 1;
    157       else if(Z > 92) Z = 92;
     164      if(Z < 1)         { Z = 1; }
     165      else if(Z > maxZ) { Z = maxZ; }
    158166      //G4cout << "Z= " << Z << G4endl;
    159167      // Initialisation
    … …  
    173181                                 const char* p)
    174182{
    175   if(data[Z]) return;
     183  if(data[Z]) { return; }
    176184  const char* path = p;
    177185  if(!p) {
    … …  
    215223    G4double sig1 = (*data[Z])[n];
    216224    dynParticle->SetKineticEnergy(emax);
    217     ggXsection->GetCrossSection(dynParticle, Elem);
    218     G4double sig2 = ggXsection->GetInelasticGlauberGribovXsc();
     225    G4double sig2 = 0.0;
     226    if(1 == Z) {
     227      fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
     228      sig2 = fNucleon->GetInelasticHadronNucleonXsc();
     229    } else {
     230      ggXsection->GetCrossSection(dynParticle, Elem);
     231      sig2 = ggXsection->GetInelasticGlauberGribovXsc();
     232    }
    219233    if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
    220234  }

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

    r1196 r1340  
    3838#include "G4Neutron.hh"
    3939#include "G4Proton.hh"
     40#include "G4HadTmpUtil.hh"
     41
    4042
    4143// Group 1: He, Be, C for 44 energies 
    … …  
    4345const G4double G4NucleonNuclearCrossSection::e1[44] =     
    4446{
    45   0.014, 0.015, 0.017, .02, 0.022, 0.025, 0.027, 0.03, 0.035, .04, 0.045, 0.05, .06, 0.07,
    46   .08, 0.09,  .1, .12, .14, .15, .16, .18, .20, .25, .30, .35, .4 , 0.5, 0.6, 0.7,  0.8, 
    47   0.9,   1, 1.5,   2,   3,   5,  7, 10,   
    48   20,   50,  100,  500, 1000
     47  0.014, 0.015, 0.017, 0.02, 0.022, 0.025, 0.027, 0.03, 0.035, 0.04,
     48  0.045, 0.05,  0.06,  0.07, 0.08,  0.09,  0.1,   0.12, 0.14,  0.15,
     49  0.16,  0.18,  0.20,  0.25, 0.30,  0.35,  0.4,   0.5,  0.6,   0.7,
     50  0.8,   0.9,   1.0,   1.5,  2.0,   3.0,   5.0,   7.0, 10.0,  20.0,
     51 50.0, 100.0, 500.0, 1000.0
    4952};
    5053
    … …  
    391394const G4double G4NucleonNuclearCrossSection::e6[46] =     
    392395{
    393   0.014, 0.015, 0.017, 0.019, .02, 0.022, 0.025, 0.027, 0.03, 0.035, .04, 0.045, 0.05,
    394   0.055, .06, 0.07, .08, 0.09, .1, .12, .14, .15, .16, .18, .20, .25, .30, .35, .4 ,
    395   0.5, 0.6, 0.7,  0.8, 0.9, 1, 1.5, 2, 3, 5, 7, 10, 20, 50, 100, 500, 1000
     396  0.014, 0.015, 0.017, 0.019, 0.02, 0.022, 0.025, 0.027, 0.03, 0.035,
     397  0.04,  0.045, 0.05,  0.055, 0.06, 0.07,  0.08,  0.09,  0.1,  0.12,
     398  0.14,  0.15,  0.16,  0.18,  0.20, 0.25,  0.30,  0.35,  0.4 , 0.5,
     399  0.6,   0.7,   0.8,   0.9,   1.0,  1.5,   2.0,   3.0,   5.0,  7.0,
     400 10.0,  20.0,  50.0, 100.0, 500.0, 1000.0
    396401};
    397402
    … …  
    443448using namespace std;
    444449
    445 //////////////////////////////////////////////////////////////////////////////////
     450///////////////////////////////////////////////////////////////////////////////
    446451//
    447452//
    … …  
    535540}
    536541
    537 ///////////////////////////////////////////////////////////////////////////////////
    538 //
     542///////////////////////////////////////////////////////////////////////////////
    539543//
    540544
    … …  
    547551////////////////////////////////////////////////////////////////////////////
    548552//
    549 //
    550553
    551554G4double G4NucleonNuclearCrossSection::
    552 GetCrossSection( const G4DynamicParticle* aParticle,
    553                  const G4Element* anElement,
    554                        G4double                )
    555 
    556 {
    557   return GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(), 0.);
     555GetCrossSection(const G4DynamicParticle* aParticle,
     556                const G4Element* anElement, G4double)
     557
     558{
     559  G4int Z = G4lrint(anElement->GetZ());
     560  G4int A = G4lrint(anElement->GetN());
     561  return GetZandACrossSection(aParticle, Z, A, 0.);
    558562}
    559563
    560564////////////////////////////////////////////////////////////////////////////
    561565//
    562 //
    563 
    564 G4bool G4NucleonNuclearCrossSection::IsApplicable(const G4DynamicParticle* aParticle,
     566
     567G4bool
     568G4NucleonNuclearCrossSection::IsApplicable(const G4DynamicParticle* aParticle,
    565569                                                  const G4Element* anElement)
    566570{
    567   return IsZAApplicable(aParticle, anElement->GetZ(), anElement->GetN());
     571  G4int Z = G4lrint(anElement->GetZ());
     572  G4int A = G4lrint(anElement->GetN());
     573  return IsIsoApplicable(aParticle, Z, A);
    568574}
    569575
    570576////////////////////////////////////////////////////////////////////////////
    571577//
    572 //
    573 
    574 G4bool G4NucleonNuclearCrossSection::IsZAApplicable(const G4DynamicParticle* aParticle,
    575                                                     G4double Z, G4double)
     578
     579G4bool
     580G4NucleonNuclearCrossSection::IsIsoApplicable(const G4DynamicParticle* aParticle,
     581                                                    G4int Z, G4int)
    576582{
    577583  G4bool result = false;
    578584  if(aParticle->GetDefinition() == theNeutron ) result = true;
    579585  if(aParticle->GetDefinition() == theProton)   result = true;
    580   if(Z < 1.5)                                   result = false;
     586  if(Z < 2)                                     result = false;
    581587  if(aParticle->GetKineticEnergy() > 999.9*GeV) result = false;
    582588  return result;
    583589}
    584590
     591
    585592////////////////////////////////////////////////////////////////////////////
    586593//
    587 //
    588594
    589595G4double G4NucleonNuclearCrossSection::
    590 GetIsoZACrossSection( const G4DynamicParticle* aParticle,
    591                       G4double  zElement, G4double, G4double  )
     596GetZandACrossSection(const G4DynamicParticle* aParticle,
     597                     G4int zElement, G4int, G4double  )
    592598{
    593599   G4double kineticEnergy = aParticle->GetKineticEnergy();
    594600 
    595601   G4double result = 0;
    596    G4int Z = G4int(zElement + 0.5);
     602   G4int Z = zElement;
    597603
    598604   // G4cout<<"Z = "<<Z<<G4endl;
    … …  
    675681
    676682/////////////////////////////////////////////////////////////////////////////
    677 //
    678683//
    679684

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

    r1337 r1340  
    2626//
    2727// The lust update: M.V. Kossov, CERN/ITEP(Moscow) 17-June-02
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     28// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2929//
    3030//
    … …  
    4040
    4141#include "G4PhotoNuclearCrossSection.hh"
     42#include "G4HadTmpUtil.hh"
    4243
    4344// Initialization of the statics
    … …  
    5960                            // Last value of the ShadowingPomeron (A-dependent)
    6061
    61 std::vector<G4double*> G4PhotoNuclearCrossSection::GDR;   // Vector of pointers to the GDRPhotonuclearCrossSection
    62 std::vector<G4double*> G4PhotoNuclearCrossSection::HEN;   // Vector of pointers to the HighEnPhotonuclearCrossSect
     62
     63// Vector of pointers to the GDRPhotonuclearCrossSection
     64std::vector<G4double*> G4PhotoNuclearCrossSection::GDR;
     65
     66// Vector of pointers to the HighEnPhotonuclearCrossSect
     67std::vector<G4double*> G4PhotoNuclearCrossSection::HEN;
     68
    6369
    6470G4PhotoNuclearCrossSection::G4PhotoNuclearCrossSection()
    65 {
    66 }
     71{}
     72
    6773
    6874G4PhotoNuclearCrossSection::~G4PhotoNuclearCrossSection()
    … …  
    9298    G4IsotopeVector* isoVector = anEle->GetIsotopeVector();
    9399    G4double* abundVector = anEle->GetRelativeAbundanceVector();
    94     G4double ZZ;
    95     G4double AA;
     100    G4int ZZ;
     101    G4int AA;
    96102 
    97103    for (G4int i = 0; i < nIso; i++) {
    98       ZZ = G4double( (*isoVector)[i]->GetZ() );
    99       AA = G4double( (*isoVector)[i]->GetN() );
     104      ZZ = (*isoVector)[i]->GetZ();
     105      AA = (*isoVector)[i]->GetN();
    100106      psig = GetIsoZACrossSection(aPart, ZZ, AA, temperature);
    101107      cross_section += psig*abundVector[i];
    … …  
    103109
    104110  } else {
    105     cross_section =
    106       GetIsoZACrossSection(aPart, anEle->GetZ(), anEle->GetN(), temperature);
     111    G4int ZZ = G4lrint(anEle->GetZ());
     112    G4int AA = G4lrint(anEle->GetN());
     113    cross_section = GetZandACrossSection(aPart, ZZ, AA, temperature);
    107114  }
    108115 
    … …  
    112119
    113120G4double
    114 G4PhotoNuclearCrossSection::GetIsoZACrossSection(const G4DynamicParticle* aPart,
    115                                                  G4double ZZ, G4double AA,
     121G4PhotoNuclearCrossSection::GetZandACrossSection(const G4DynamicParticle* aPart,
     122                                                 G4int ZZ, G4int AA,
    116123                                                 G4double /*temperature*/)
    117124{
    … …  
    146153  //
    147154  const G4double Energy = aPart->GetKineticEnergy()/MeV;
    148   const G4int targetAtomicNumber = static_cast<int>(AA+.499); //@@ Nat mixture (?!)
    149   const G4int targZ = static_cast<int>(ZZ);
     155  const G4int targetAtomicNumber = AA; //@@ Nat mixture (?!)
     156  const G4int targZ = ZZ;
    150157  const G4int targN = targetAtomicNumber-targZ;
    151158#ifdef debug
    … …  
    16091616
    16101617  static const G4double SH10[nH]={
    1611     3.918292e+00,3.904931e+00,3.893792e+00,3.886847e+00,3.886858e+00,3.897612e+00,3.924175e+00,
    1612     3.973155e+00,4.052892e+00,4.173448e+00,4.346251e+00,4.583168e+00,4.894929e+00,5.289011e+00,
    1613     5.767472e+00,6.325587e+00,6.952077e+00,7.631192e+00,8.346046e+00,9.081993e+00,9.828955e+00,
    1614     1.058224e+01,1.134205e+01,1.211228e+01,1.289914e+01,1.370979e+01,1.455140e+01,1.543028e+01,
    1615     1.635126e+01,1.731704e+01,1.832759e+01,1.937949e+01,2.046524e+01,2.157253e+01,2.268371e+01,
     1618    3.918292, 3.904931, 3.893792, 3.886847, 3.886858, 3.897612, 3.924175,
     1619    3.973155, 4.052892, 4.173448, 4.346251, 4.583168, 4.894929, 5.289011,
     1620    5.767472, 6.325587, 6.952077, 7.631192, 8.346046, 9.081993, 9.828955,
     1621   10.58224, 11.34205, 12.11228, 12.89914, 13.70979, 14.55140, 15.43028,
     1622   16.35126, 17.31704, 18.32759, 19.37949, 20.46524, 21.57253, 22.68371,
    16161623    2.377554e+01,2.481942e+01,2.578236e+01,2.662884e+01,2.732355e+01,2.783477e+01,2.813791e+01,
    16171624    2.821866e+01,2.807489e+01,2.771701e+01,2.716665e+01,2.645395e+01,2.561414e+01,2.468397e+01,

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

    r1228 r1340  
    7575 const G4double G4PiNuclearCrossSection::e2[39] = {
    7676  .02, .04, .06, .08, .10, .11, .12, .13, .14,  .15, .16, .17, .18, .20, .22,
    77   .24, .26, .28, .30, .35, .40, .45, .50, .55, .575, .60, .70, .80, .90,   1,   
     77  .24, .26, .28, .30, .35, .40, .45, .50, .55, .575, .60, .70, .80, .90,   1,
    7878    2,   3,   5,  10,  20,  50, 100, 500, 100000};
    7979
    … …  
    138138
    139139 const G4double G4PiNuclearCrossSection::na_m_t[31] = {
    140   450, 545, 705, 910, 1020, 1075, 1087, 1080, 1042, 987, 943, 885, 790, 700, 650,
    141   610, 585, 575, 585,  595,  600,  610,  556,  524, 494, 458, 445, 429, 427, 427,
    142   427};
     140  450, 545, 705, 910, 1020, 1075, 1087, 1080, 1042, 987, 943, 885, 790, 700,
     141  650, 610, 585, 575, 585,  595,  600,  610,  556,  524, 494, 458, 445, 429,
     142  427, 427, 427};
    143143
    144144 const G4double G4PiNuclearCrossSection::na_m_in[31] = {
    … …  
    156156
    157157 const G4double G4PiNuclearCrossSection::e3_1[31] = {
    158   0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.25, 0.30, 0.35, 0.40,
    159   0.45, 0.50, 0.60, 0.70, 0.80, 0.90,    1,    2,    3,    5,   10,   20,   50,  100,  500,
    160   100000};
     158  0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20,
     159  0.22, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.60, 0.70, 0.80,
     160  0.90, 1.0,  2.0,  3.0,  5.0, 10.0, 20.0, 50.0, 100.0, 500.0, 100000.0};
    161161
    162162 const G4double G4PiNuclearCrossSection::al_m_t[31] = {
    163   532, 637, 832, 1057, 1207, 1230, 1210, 1174, 1133, 1095, 1038, 970, 890, 807, 750,
    164   710, 675, 665,  670,  673,  678,  682,  618,  574,  546,  520, 507, 495, 488, 488, 
    165   488};
     163  532, 637, 832, 1057, 1207, 1230, 1210, 1174, 1133, 1095,
     164 1038, 970, 890,  807,  750,  710,  675,  665,  670,  673,
     165  678, 682, 618,  574,  546,  520,  507,  495,  488,  488,  488};
    166166
    167167 const G4double G4PiNuclearCrossSection::al_m_in[31] = {
    … …  
    171171
    172172 const G4double G4PiNuclearCrossSection::al_p_t[21] = {
    173   225, 350, 616, 945, 1122, 1175, 1157, 1128, 1088, 1045, 988, 935, 870, 787, 730,
    174   690, 660, 652, 660,  668,  678};
     173  225, 350, 616, 945, 1122, 1175, 1157, 1128, 1088, 1045,
     174  988, 935, 870, 787, 730,   690,  660,  652, 660,  668,  678};
    175175
    176176 const G4double G4PiNuclearCrossSection::al_p_in[21] = {
    … …  
    179179
    180180 const G4double G4PiNuclearCrossSection::ca_m_t[31] = {
    181   800, 980, 1240, 1460, 1570, 1600, 1580, 1535, 1475, 1425, 1375, 1295, 1200, 1083, 1000,
    182   948, 915,  895,  900,  908,  915,  922,  856,  795,  740,  705,  682,  660,  660,  660,
    183   660};
     181  800,  980, 1240, 1460, 1570, 1600, 1580, 1535, 1475, 1425,
     182 1375, 1295, 1200, 1083, 1000,  948,  915,  895,  900,  908,
     183  915,  922,  856,  795,  740,  705,  682,  660,  660,  660, 660};
    184184
    185185 const G4double G4PiNuclearCrossSection::ca_m_in[31] = {
    … …  
    222222 const G4double G4PiNuclearCrossSection::cd_m_t[34] = {3060, 3125, 3170, 3220, 3255, 3280, 3290, 3260, 3270, 3200, 3120, 3080, 3090, 2920, 2810, 2640, 2362, 2230, 2115, 2050, 2020, 2025, 2040, 2070, 2100, 1900, 1795, 1740, 1675, 1645, 1625, 1620, 1620, 1620};
    223223 const G4double G4PiNuclearCrossSection::cd_m_in[34]= {1025, 1275, 1440, 1625, 1740, 1800, 1880, 1920, 1980, 1920, 1850, 1810, 1720, 1650, 1560, 1450, 1330, 1290, 1245, 1210, 1200, 1200, 1205, 1205, 1230, 1130, 1085, 1060, 1000,  985,  975,  970,  970,  970};
    224  const G4double G4PiNuclearCrossSection::cd_p_t[28] =  {455,  780, 1170, 1700, 2120, 2400, 2600, 2720, 2820, 2840, 2800, 2760, 2720, 2640, 2560, 2450, 2252, 2130, 2035, 1985, 1970, 1975, 2005, 2035, 2070, 1880, 1795, 1740};
    225  const G4double G4PiNuclearCrossSection::cd_p_in[28] = {310,  580,  880, 1060, 1270, 1400, 1530, 1610, 1660, 1680, 1640, 1600, 1560, 1500, 1430, 1330, 1280, 1230, 1200, 1180, 1170, 1175, 1180, 1180, 1210, 1120, 1085, 1060};
     224
     225 const G4double G4PiNuclearCrossSection::cd_p_t[28] =  {
     226   455,  780, 1170, 1700, 2120, 2400, 2600, 2720, 2820, 2840,
     227  2800, 2760, 2720, 2640, 2560, 2450, 2252, 2130, 2035, 1985,
     228  1970, 1975, 2005, 2035, 2070, 1880, 1795, 1740};
     229
     230 const G4double G4PiNuclearCrossSection::cd_p_in[28] = {
     231   310,  580,  880, 1060, 1270, 1400, 1530, 1610, 1660, 1680,
     232  1640, 1600, 1560, 1500, 1430, 1330, 1280, 1230, 1200, 1180,
     233  1170, 1175, 1180, 1180, 1210, 1120, 1085, 1060};
    226234
    227235 const G4double G4PiNuclearCrossSection::e6[35] = {
    228   0.02, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.25, 
    229   0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.70, 0.80, 0.90,    1,    2,    3,    5,   10,   
    230     20,   50,  100,  500, 100000};
    231 
    232  const G4double G4PiNuclearCrossSection::sn_m_t[35] =  {3000, 3180, 3250, 3300, 3300, 3410, 3470, 3450, 3410, 3350, 3280, 3200, 3120, 3050, 2900, 2630, 2500, 2325, 2190, 2100, 2060, 2055, 2055, 2055, 2067, 2085, 2000, 1900, 1835, 1770, 1720, 1700, 1695, 1695, 1695};
    233  const G4double G4PiNuclearCrossSection::sn_m_in[35] = {1050, 1350, 1520, 1650, 1800, 1980, 2070, 2120, 2090, 2050, 1980, 1920, 1830, 1770, 1670, 1500, 1435, 1350, 1300, 1230, 1220, 1235, 1235, 1235, 1237, 1240, 1160, 1120, 1090, 1065, 1040, 1020, 1015, 1015, 1015};
     236  0.02, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.12, 0.14,
     237  0.16, 0.18, 0.20, 0.22, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50,
     238  0.55, 0.60, 0.70, 0.80, 0.90, 1.0,  2.0,  3.0,  5.0, 10.0,   
     239 20.0, 50.0, 100.0, 500.0, 100000.0};
     240
     241 const G4double G4PiNuclearCrossSection::sn_m_t[35] = {
     242  3000, 3180, 3250, 3300, 3300, 3410, 3470, 3450, 3410, 3350,
     243  3280, 3200, 3120, 3050, 2900, 2630, 2500, 2325, 2190, 2100,
     244  2060, 2055, 2055, 2055, 2067, 2085, 2000, 1900, 1835, 1770,
     245  1720, 1700, 1695, 1695, 1695};
     246
     247 const G4double G4PiNuclearCrossSection::sn_m_in[35] = {
     248  1050, 1350, 1520, 1650, 1800, 1980, 2070, 2120, 2090, 2050,
     249  1980, 1920, 1830, 1770, 1670, 1500, 1435, 1350, 1300, 1230,
     250  1220, 1235, 1235, 1235, 1237, 1240, 1160, 1120, 1090, 1065,
     251  1040, 1020, 1015, 1015, 1015};
     252
    234253 const G4double G4PiNuclearCrossSection::sn_p_t[29] =  { 465,  800, 1200, 1760, 2170, 2480, 2730, 2885, 2970, 2980, 2970, 2890, 2840, 2790, 2620, 2450, 2335, 2205, 2080, 2020, 2010, 1990, 1990, 2015, 2030, 2045, 1980, 1890, 1835};
    235254 const G4double G4PiNuclearCrossSection::sn_p_in[29] = { 315,  590,  880, 1220, 1460, 1580, 1700, 1770, 1810, 1810, 1800, 1730, 1680, 1630, 1530, 1400, 1335, 1270, 1210, 1180, 1190, 1190, 1190, 1205, 1210, 1210, 1150, 1115, 1090};
    … …  
    346365
    347366G4double G4PiNuclearCrossSection::
    348 GetIsoZACrossSection(const G4DynamicParticle* particle, G4double ZZ,
    349                      G4double /*AA*/, G4double /*temperature*/)
     367GetZandACrossSection(const G4DynamicParticle* particle, G4int ZZ,
     368                     G4int /*AA*/, G4double /*temperature*/)
    350369{
    351370  // precondition
    352   //  G4ping debug("debug_PiNuclearCrossSection");
    353371  using namespace std;
    354372  G4bool ok = false;
    … …  
    367385
    368386  G4double result = 0;
    369   G4int Z=G4lrint(ZZ);
     387  G4int Z = ZZ;
    370388  //  debug.push_back(Z);
    371389  size_t it = 0;
    372390
    373   while(it<theZ.size() && Z>theZ[it]) it++;
     391  while(it < theZ.size() && Z > theZ[it]) it++;
    374392
    375393  //  debug.push_back(theZ[it]);
    … …  
    475493 }
    476494
    477  G4double G4PiNuclearCrossSection::
    478  Interpolate(G4int Z1, G4int Z2, G4int Z, G4double x1, G4double x2)
    479  {
     495
     496G4double G4PiNuclearCrossSection::
     497Interpolate(G4int Z1, G4int Z2, G4int Z, G4double x1, G4double x2)
     498{
    480499//   Nucleon numbers obtained from G4NistManager G4 8.0
    481  static const G4double A[92]=
    482                         {1.0001, 4.0000, 6.9241, 9.0000, 10.801, 12.011, 14.004, 16.004, 19.000, 20.188,
    483                          23.000, 24.320, 27.000, 28.109, 31.000, 32.094, 35.484, 39.985, 39.135, 40.116,
    484                          45.000, 47.918, 50.998, 52.055, 55.000, 55.910, 59.000, 58.760, 63.617, 65.468,
    485                          69.798, 72.691, 75.000, 79.042, 79.986, 83.887, 85.557, 87.710, 89.000, 91.318,
    486                          93.000, 96.025, 98.000, 101.16, 103.00, 106.51, 107.96, 112.51, 114.91, 118.81,
    487                          121.86, 127.70, 127.00, 131.39, 133.00, 137.42, 139.00, 140.21, 141.00, 144.32,
    488                          145.00, 150.45, 152.04, 157.33, 159.00, 162.57, 165.00, 167.32, 169.00, 173.10,
    489                          175.03, 178.54, 181.00, 183.89, 186.25, 190.27, 192.25, 195.11, 197.00, 200.63,
    490                          204.41, 207.24, 209.00, 209.00, 210.00, 222.00, 223.00, 226.00, 227.00, 232.00,
    491                          231.00, 237.98};
     500 static const G4double A[92] = {
     501   1.0001, 4.0000, 6.9241, 9.000, 10.801, 12.011, 14.004, 16.004, 19.000,
     502  20.188, 23.000, 24.320, 27.000, 28.109, 31.000, 32.094, 35.484, 39.985,
     503  39.135, 40.116, 45.000, 47.918, 50.998, 52.055, 55.000, 55.910, 59.000,
     504  58.760, 63.617, 65.468, 69.798, 72.691, 75.000, 79.042, 79.986, 83.887,
     505  85.557, 87.710, 89.000, 91.318, 93.000, 96.025, 98.000, 101.16, 103.00,
     506 106.51, 107.96, 112.51, 114.91, 118.81, 121.86, 127.70,  127.00, 131.39,
     507 133.00, 137.42, 139.00, 140.21, 141.00, 144.32, 145.00,  150.45, 152.04,
     508 157.33, 159.00, 162.57, 165.00, 167.32, 169.00, 173.10,  175.03, 178.54,
     509 181.00, 183.89, 186.25, 190.27, 192.25, 195.11, 197.00, 200.63,  204.41,
     510 207.24, 209.00, 209.00, 210.00, 222.00, 223.00, 226.00, 227.00,  232.00,
     511 231.00, 237.98};
    492512                         
    493513 static G4bool NeedInit=true;               
    … …  
    502522 }
    503523
    504 //        for tabulated data, cross section scales with A^.75
     524// for tabulated data, cross section scales with A^.75
    505525   G4double r1 = x1 / A75[Z1-1] * A75[Z-1];
    506526   G4double r2 = x2 / A75[Z2-1] * A75[Z-1];
    507527   G4double result=0.5*(r1+r2);
    508 //       G4cout << "x1/2, z1/2 z" <<x1<<" "<<x2<<" "<<Z1<<" "<<Z2<<" "<<Z<<G4endl;
    509 //       G4cout << "res1/2 " << r1 <<" " << r2 <<" " << result<< G4endl;
     528//  G4cout << "x1/2, z1/2 z" <<x1<<" "<<x2<<" "<<Z1<<" "<<Z2<<" "<<Z<<G4endl;
     529//  G4cout << "res1/2 " << r1 <<" " << r2 <<" " << result<< G4endl;
    510530   return result;
    511  }
     531}

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

    r1228 r1340  
    3131
    3232#include "G4ProtonInelasticCrossSection.hh"
     33#include "G4HadTmpUtil.hh"
    3334#include "globals.hh"
    3435
    … …  
    4647    G4IsotopeVector* isoVector = anEle->GetIsotopeVector();
    4748    G4double* abundVector = anEle->GetRelativeAbundanceVector();
    48     G4double ZZ;
    49     G4double AA;
     49    G4int ZZ;
     50    G4int AA;
    5051 
    5152    for (G4int i = 0; i < nIso; i++) {
    52       ZZ = G4double( (*isoVector)[i]->GetZ() );
    53       AA = G4double( (*isoVector)[i]->GetN() );
     53      ZZ = (*isoVector)[i]->GetZ();
     54      AA = (*isoVector)[i]->GetN();
    5455      psig = GetCrossSection(KE, AA, ZZ);
    5556      cross_section += psig*abundVector[i];
    … …  
    5758 
    5859  } else {
    59     cross_section = GetCrossSection(KE, anEle->GetN(), anEle->GetZ());
     60    G4int ZZ = G4lrint(anEle->GetZ());
     61    G4int AA = G4lrint(anEle->GetN());
     62    cross_section = GetCrossSection(KE, AA, ZZ);
    6063  }
    6164
    … …  
    6366}
    6467
    65    
     68
    6669G4double G4ProtonInelasticCrossSection::
    67 GetCrossSection(G4double kineticEnergy, G4double atomicNumber, G4double nOfProtons)
     70GetCrossSection(G4double kineticEnergy, G4int atomicNumber, G4int nOfProtons)
    6871{   
    6972  if (kineticEnergy > 19.9*GeV )
    … …  
    7174    return  GetCrossSection(19.8*GeV,atomicNumber,nOfProtons);
    7275  }
    73   G4double nOfNeutrons = atomicNumber-nOfProtons;
     76  G4int nOfNeutrons = atomicNumber-nOfProtons;
    7477  kineticEnergy /=GeV;
    7578  G4double a = atomicNumber;
    … …  
    8083  G4double fac1=b0*(1-std::pow(a,-0.3333));
    8184  G4double fac2=1.;
    82   if(nOfNeutrons>1.5) fac2=std::log((nOfNeutrons));
     85  if(nOfNeutrons > 1) fac2=std::log((G4double(nOfNeutrons)));
    8386  G4double crossSection = 1E31*fac*fac2*(1+std::pow(a,0.3333)-fac1);
    8487

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

    r1196 r1340  
    4343
    4444G4double G4TripathiCrossSection::
    45 GetIsoZACrossSection(const G4DynamicParticle* aPart, G4double ZZ, G4double AA,
    46                 G4double /*temperature*/)
     45GetZandACrossSection(const G4DynamicParticle* aPart, G4int ZZ, G4int AA,
     46                     G4double /*temperature*/)
    4747{
    4848  G4double result = 0.;
    49  
    5049  const G4double targetAtomicNumber = AA;
    5150  const G4double nTargetProtons = ZZ;
    … …  
    157156    G4IsotopeVector* isoVector = anEle->GetIsotopeVector();
    158157    G4double* abundVector = anEle->GetRelativeAbundanceVector();
    159     G4double ZZ;
    160     G4double AA;
     158    G4int ZZ;
     159    G4int AA;
    161160     
    162161    for (G4int i = 0; i < nIso; i++) {
    163       ZZ = G4double( (*isoVector)[i]->GetZ() );
    164       AA = G4double( (*isoVector)[i]->GetN() );
    165       sig = GetIsoZACrossSection(aPart, ZZ, AA, temperature);
     162      ZZ = (*isoVector)[i]->GetZ();
     163      AA = (*isoVector)[i]->GetN();
     164      sig = GetZandACrossSection(aPart, ZZ, AA, temperature);
    166165      xsection += sig*abundVector[i];
    167166    }
    168167   
    169168  } else {
    170     xsection =
    171       GetIsoZACrossSection(aPart, anEle->GetZ(), anEle->GetN(),
    172                           temperature);
     169    G4int ZZ = G4lrint(anEle->GetZ());
     170    G4int AA = G4lrint(anEle->GetN());
     171    xsection = GetIsoZACrossSection(aPart, ZZ, AA, temperature);
    173172  }
    174173

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

    r1228 r1340  
    6363#include "G4ParticleTable.hh"
    6464#include "G4IonTable.hh"
     65#include "G4HadTmpUtil.hh"
     66
    6567
    6668///////////////////////////////////////////////////////////////////////////////
    … …  
    101103  (const G4DynamicParticle* theProjectile, const G4Element* theTarget)
    102104{
    103   return IsZAApplicable(theProjectile, theTarget->GetZ(), theTarget->GetN());
    104 }
    105 
    106 
    107 G4bool G4TripathiLightCrossSection::IsZAApplicable
    108   (const G4DynamicParticle* theProjectile, G4double ZZ, G4double AA)
     105  G4int Z = G4lrint(theTarget->GetZ());
     106  G4int A = G4lrint(theTarget->GetN());
     107  return IsIsoApplicable(theProjectile, Z, A);
     108}
     109
     110
     111G4bool
     112G4TripathiLightCrossSection::IsIsoApplicable(const G4DynamicParticle* theProjectile,
     113                                             G4int ZZ, G4int AA)
    109114{
    110115  G4bool result = false;
    … …  
    122127  return result;
    123128}
    124 ///////////////////////////////////////////////////////////////////////////////
    125 //
    126 G4double G4TripathiLightCrossSection::GetIsoZACrossSection
    127   (const G4DynamicParticle* theProjectile, G4double ZZ, G4double AA,
    128    G4double /*theTemperature*/)
     129
     130///////////////////////////////////////////////////////////////////////////////
     131//
     132G4double
     133G4TripathiLightCrossSection::GetZandACrossSection(const G4DynamicParticle* theProjectile,
     134                                                  G4int ZZ, G4int AA, G4double /*theTemperature*/)
    129135{
    130136//
    … …  
    301307          SetLowEnergyCheck(true);
    302308      G4double resultp =
    303         GetIsoZACrossSection(&slowerProjectile, ZZ, AA, 0.0);
     309        GetZandACrossSection(&slowerProjectile, ZZ, AA, 0.0);
    304310          SetLowEnergyCheck(savelowenergy);
    305311      if (resultp >result) result = 0.0;
    … …  
    322328    G4IsotopeVector* isoVector = theTarget->GetIsotopeVector();
    323329    G4double* abundVector = theTarget->GetRelativeAbundanceVector();
    324     G4double ZZ;
    325     G4double AA;
     330    G4int ZZ;
     331    G4int AA;
    326332     
    327333    for (G4int i = 0; i < nIso; i++) {
    328       ZZ = G4double( (*isoVector)[i]->GetZ() );
    329       AA = G4double( (*isoVector)[i]->GetN() );
    330       sig = GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
     334      ZZ = (*isoVector)[i]->GetZ();
     335      AA = (*isoVector)[i]->GetN();
     336      sig = GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
    331337      xsection += sig*abundVector[i];
    332338    }
    333339   
    334340  } else {
    335     xsection =
    336       GetIsoZACrossSection(theProjectile, theTarget->GetZ(), theTarget->GetN(),
    337                           theTemperature);
     341    G4int ZZ = G4lrint(theTarget->GetZ());
     342    G4int AA = G4lrint(theTarget->GetN());
     343    xsection = GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
    338344  }
    339345   
    … …  
    348354  lowEnergyCheck = aLowEnergyCheck;
    349355}
    350 ///////////////////////////////////////////////////////////////////////////////
    351 //
     356

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

    r962 r1340  
    8080//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    8181
    82 G4bool G4UElasticCrossSection::IsZAApplicable(const G4DynamicParticle* dp,
    83                                                 G4double Z, G4double A)
     82// G4bool G4UElasticCrossSection::IsZAApplicable(const G4DynamicParticle* dp,
     83//                                              G4double Z, G4double A)
     84G4bool G4UElasticCrossSection::IsIsoApplicable(const G4DynamicParticle* dp,
     85                                               G4int Z, G4int A)
    8486{
    8587  G4bool res = false;
    … …  
    9496}
    9597
    96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    97 
    98 G4double G4UElasticCrossSection::GetCrossSection(const G4DynamicParticle* dp,
    99                                                    const G4Element* elm,
    100                                                    G4double temp)
    101 {
    102   return GetIsoZACrossSection(dp, elm->GetZ(), elm->GetN(), temp);
    103 }
    104 
    105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    106 
    107 G4double G4UElasticCrossSection::GetIsoZACrossSection(const G4DynamicParticle* dp,
    108                                                         G4double Z,
    109                                                         G4double A,
    110                                                         G4double)
     98
     99G4double
     100G4UElasticCrossSection::GetCrossSection(const G4DynamicParticle* dp,
     101                                        const G4Element* elm, G4double temp)
     102{
     103  G4int Z = G4lrint(elm->GetZ());
     104  G4int N = G4lrint(elm->GetN());
     105
     106  return GetZandACrossSection(dp, Z, N, temp);
     107}
     108
     109
     110G4double
     111G4UElasticCrossSection::GetZandACrossSection(const G4DynamicParticle* dp,
     112                                             G4int Z, G4int A, G4double)
    111113{
    112114  G4double cross = 0.0;
    113115  G4double ekin = dp->GetKineticEnergy();
    114   G4int iz = G4int(Z);
    115   if(iz > 92) iz = 92;
     116//  G4int iz = G4int(Z);
     117  if(Z > 92) Z = 92;
    116118
    117119    // proton and neutron
    118120  if(fNucleon) {
    119     if(iz == 1) cross = fGheisha->GetElasticCrossSection(dp, Z, A);
     121    if(Z == 1) cross = fGheisha->GetElasticCrossSection(dp, Z, A);
    120122    else if(ekin > thEnergy) {
    121       cross = theFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
     123      cross = theFac[Z]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
    122124    } else {
    123125      cross = fNucleon->GetElasticCrossSection(dp, Z, A);
    … …  
    126128    // pions
    127129  } else if(fUPi) {
    128     if(iz == 1) cross = fGheisha->GetElasticCrossSection(dp, Z, A);
     130    if(Z == 1) cross = fGheisha->GetElasticCrossSection(dp, Z, A);
    129131    else if(ekin > thEnergy) {
    130       cross = theFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
     132      cross = theFac[Z]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
    131133    } else {
    132134      cross = fUPi->GetElasticCrossSection(dp, Z, A);
    … …  
    136138  } else {
    137139    if(hasGlauber && ekin > thEnergy) {
    138       cross = theFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
     140      cross = theFac[Z]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
    139141    } else if(fGheisha->IsApplicable(dp, Z, A)) {
    140142      cross = fGheisha->GetElasticCrossSection(dp, Z, A);
    … …  
    182184
    183185  G4NistManager* nist = G4NistManager::Instance();
    184   G4double A = nist->GetAtomicMassAmu(2);
     186  G4int A = G4lrint(nist->GetAtomicMassAmu(2));
    185187
    186188  if(fGlauber->IsZAApplicable(&dp, 2.0, A)) {
    … …  
    193195
    194196      G4double Z = G4double(iz);
    195       A = nist->GetAtomicMassAmu(iz);
    196       csup = fGlauber->GetElasticGlauberGribov(&dp, Z, A);
     197      A = G4lrint(nist->GetAtomicMassAmu(iz));
     198      csup = fGlauber->GetElasticGlauberGribov(&dp, iz, A);
    197199      // proton and neutron
    198200      if(fNucleon) {
    199         csdn = fNucleon->GetElasticCrossSection(&dp, Z, A);
     201        csdn = fNucleon->GetElasticCrossSection(&dp, iz, A);
    200202
    201203        // pions
    202204      } else if(fUPi) {
    203         csdn = fUPi->GetElasticCrossSection(&dp, Z, A);
     205        csdn = fUPi->GetElasticCrossSection(&dp, iz, A);
    204206
    205207        // other
    206       } else if(fGheisha->IsApplicable(&dp, Z, A)) {
    207         csdn = fGheisha->GetElasticCrossSection(&dp, Z, A);
     208      } else if(fGheisha->IsApplicable(&dp, iz, A)) {
     209        csdn = fGheisha->GetElasticCrossSection(&dp, iz, A);
    208210
    209211      } else {
    … …  
    217219}
    218220
    219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    220 
    221 
     221

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

    r962 r1340  
    7878}
    7979
    80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    81 
    82 G4bool G4UInelasticCrossSection::IsZAApplicable(const G4DynamicParticle* dp,
    83                                                 G4double Z, G4double A)
     80
     81G4bool
     82G4UInelasticCrossSection::IsIsoApplicable(const G4DynamicParticle* dp,
     83                                          G4int Z, G4int A)
    8484{
    8585  G4bool res = false;
    … …  
    9494}
    9595
    96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    97 
    98 G4double G4UInelasticCrossSection::GetCrossSection(const G4DynamicParticle* dp,
    99                                                    const G4Element* elm,
    100                                                    G4double temp)
    101 {
    102   return GetIsoZACrossSection(dp, elm->GetZ(), elm->GetN(), temp);
    103 }
    104 
    105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    106 
    107 G4double G4UInelasticCrossSection::GetIsoZACrossSection(const G4DynamicParticle* dp,
    108                                                         G4double Z,
    109                                                         G4double A,
    110                                                         G4double)
     96
     97G4double
     98G4UInelasticCrossSection::GetCrossSection(const G4DynamicParticle* dp,
     99                                          const G4Element* elm,
     100                                          G4double temp)
     101{
     102  G4int Z = G4lrint(elm->GetZ());
     103  G4int N = G4lrint(elm->GetN());
     104  return GetIsoZACrossSection(dp, Z, N, temp);
     105}
     106
     107
     108G4double
     109G4UInelasticCrossSection::GetZandACrossSection(const G4DynamicParticle* dp,
     110                                               G4int Z, G4int A, G4double)
    111111{
    112112  G4double cross = 0.0;
    113113  G4double ekin = dp->GetKineticEnergy();
    114   G4int iz = G4int(Z);
    115   if(iz > 92) iz = 92;
     114  //  G4int iz = G4int(Z);
     115  if(Z > 92) Z = 92;
    116116
    117117    // proton and neutron
    118118  if(fNucleon) {
    119     if(iz == 1) cross = fGheisha->GetInelasticCrossSection(dp, Z, A);
     119    if(Z == 1) cross = fGheisha->GetInelasticCrossSection(dp, Z, A);
    120120    else if(ekin > thEnergy) {
    121       cross = theFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
     121      cross = theFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
    122122    } else {
    123123      cross = fNucleon->GetIsoZACrossSection(dp, Z, A);
    … …  
    126126    // pions
    127127  } else if(fUPi) {
    128     if(iz == 1) cross = fGheisha->GetInelasticCrossSection(dp, Z, A);
     128    if(Z == 1) cross = fGheisha->GetInelasticCrossSection(dp, Z, A);
    129129    else if(ekin > thEnergy) {
    130       cross = theFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
     130      cross = theFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
    131131    } else {
    132132      cross = fUPi->GetInelasticCrossSection(dp, Z, A);
    … …  
    136136  } else {
    137137    if(hasGlauber && ekin > thEnergy) {
    138       cross = theFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
     138      cross = theFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
    139139    } else if(fGheisha->IsApplicable(dp, Z, A)){
    140140      cross = fGheisha->GetInelasticCrossSection(dp, Z, A);
    … …  
    153153}
    154154
    155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    156155
    157156void G4UInelasticCrossSection::BuildPhysicsTable(const G4ParticleDefinition& p)
    … …  
    173172
    174173
    175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    176 
    177174void G4UInelasticCrossSection::Initialise(const G4ParticleDefinition* p)
    178175{
    … …  
    182179
    183180  G4NistManager* nist = G4NistManager::Instance();
    184   G4double A = nist->GetAtomicMassAmu(2);
    185 
    186   if(fGlauber->IsZAApplicable(&dp, 2.0, A)) {
     181  G4int A = G4lrint(nist->GetAtomicMassAmu(2));
     182
     183  if(fGlauber->IsIsoApplicable(&dp, 2, A)) {
    187184    hasGlauber = true;
    188185
    … …  
    193190
    194191      G4double Z = G4double(iz);
    195       A = nist->GetAtomicMassAmu(iz);
    196       csup = fGlauber->GetInelasticGlauberGribov(&dp, Z, A);
     192      A = G4lrint(nist->GetAtomicMassAmu(iz));
     193      csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A);
    197194
    198195      // proton and neutron
    199196      if(fNucleon) {
    200         csdn = fNucleon->GetIsoZACrossSection(&dp, Z, A);
     197        csdn = fNucleon->GetIsoZACrossSection(&dp, iz, A);
    201198
    202199        // pions
    203200      } else if(fUPi) {
    204         csdn = fUPi->GetInelasticCrossSection(&dp, Z, A);
     201        csdn = fUPi->GetInelasticCrossSection(&dp, iz, A);
    205202
    206203        // other
    207       } else if(fGheisha->IsApplicable(&dp, Z, A)) {
    208         csdn = fGheisha->GetInelasticCrossSection(&dp, Z, A);
     204      } else if(fGheisha->IsApplicable(&dp, iz, A)) {
     205        csdn = fGheisha->GetInelasticCrossSection(&dp, iz, A);
    209206
    210207      } else {
    … …  
    218215}
    219216
    220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    221 
    222 

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

    r1337 r1340  
    6060}
    6161
    62 G4double G4UPiNuclearCrossSection::GetElasticCrossSection(
    63          const G4DynamicParticle* dp, G4double Z, G4double A)
     62G4double
     63G4UPiNuclearCrossSection::GetElasticCrossSection(const G4DynamicParticle* dp,
     64                                                 G4int Z, G4int A)
    6465{
    6566  G4double cross = 0.0;
    … …  
    7374}
    7475
    75 G4double G4UPiNuclearCrossSection::GetInelasticCrossSection(
    76          const G4DynamicParticle* dp, G4double Z, G4double A)
     76G4double
     77G4UPiNuclearCrossSection::GetInelasticCrossSection(const G4DynamicParticle* dp,
     78                                                   G4int Z, G4int A)
    7779{
    7880  G4double cross = 0.0;
    … …  
    8284  const G4ParticleDefinition* part = dp->GetDefinition();
    8385
    84   // Coulomb barier
     86  // Coulomb barrier
    8587  if(part == piPlus) {
    8688    if(ekin > elowest) {
    … …  
    101103
    102104G4double G4UPiNuclearCrossSection::Interpolate(
    103          G4double Z, G4double A, G4double ekin, G4PhysicsTable* table)
     105         G4int Z, G4int A, G4double ekin, G4PhysicsTable* table)
    104106{
    105107  G4double res = 0.0;
    106108  G4int idx;
    107   G4int iz = G4int(Z + 0.5);
     109  G4int iz = Z;
    108110  if(iz > 92) iz = 92;
    109111  for(idx=0; idx<NZ; idx++) {if(theZ[idx] >= iz) break;}
    … …  
    123125    G4int iz1 = theZ[idx-1];
    124126    G4double x1 = (((*table)[idx-1])->Value(ekin))*APower[iz]/APower[iz1];
    125     G4double w1 = A - theA[idx-1];
    126     G4double w2 = theA[idx] - A;
     127    G4double w1 = G4double(A) - theA[idx-1];
     128    G4double w2 = theA[idx] - G4double(A);
    127129    res = (w1*x2 + w2*x1)/(w1 + w2);
    128130  }

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VCrossSectionDataSet.cc,v 1.8 2009/01/24 11:54:47 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VCrossSectionDataSet.cc,v 1.9 2010/07/05 13:39:11 vnivanch Exp $
     27// GEANT4 tag $Name: hadr-cross-V09-03-12 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4242#include "G4VCrossSectionDataSet.hh"
    4343#include "G4CrossSectionDataSetRegistry.hh"
     44#include "G4Pow.hh"
    4445
    45 G4VCrossSectionDataSet::G4VCrossSectionDataSet() :
    46   verboseLevel(0)
     46G4VCrossSectionDataSet::G4VCrossSectionDataSet(const G4String& nam) :
     47  verboseLevel(0),minKinEnergy(0.0),maxKinEnergy(DBL_MAX),name(nam)
    4748{
    4849  G4CrossSectionDataSetRegistry::Instance()->Register(this);
    … …  
    5455}
    5556
    56 // Override this method to test for particle and isotope applicability
    57 
     57// Override thess methods to test for particle and isotope applicability
    5858G4bool
    5959G4VCrossSectionDataSet::IsZAApplicable(const G4DynamicParticle*,
    … …  
    6363}
    6464
     65G4bool
     66G4VCrossSectionDataSet::IsIsoApplicable(const G4DynamicParticle*,
     67                                        G4int /*Z*/, G4int /*A*/)
     68{
     69  return true;
     70}
    6571
    6672G4double
    … …  
    8187                                             G4double /*aTemperature*/)
    8288{
    83   return 62*std::pow(AA, 2./3.)*millibarn;
     89  return 62*G4Pow::GetInstance()->A23(AA)*millibarn;
    8490}
     91
     92G4double
     93G4VCrossSectionDataSet::GetZandACrossSection(const G4DynamicParticle*,
     94                                             G4int /*Z*/, G4int N,
     95                                             G4double /*aTemperature*/)
     96{
     97  return 62*G4Pow::GetInstance()->Z23(N)*millibarn;
     98}
Note: See TracChangeset for help on using the changeset viewer.