Ignore:
Timestamp:
May 28, 2009, 4:26:57 PM (15 years ago)
Author:
garnier
Message:

maj sur la beta de geant 4.9.3

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundProton.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 //J.M.Quesada (August 08). New source file
    27 //
    28 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse
    29 // cross section option
    30 //
    31 // JMQ (06 September 2008) Also external choice has been added for:
    32 //                      - superimposed Coulomb barrier (if useSICB=true)
     26//
     27// $Id: G4PreCompoundProton.cc,v 1.4 2009/02/11 18:06:00 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     29//
     30// -------------------------------------------------------------------
     31//
     32// GEANT4 Class file
     33//
     34//
     35// File name:     G4PreCompoundProton
     36//
     37// Author:         V.Lara
     38//
     39// Modified: 
     40// 21.08.2008 J. M. Quesada added external choice of inverse cross section option
     41// 21.08.2008 J. M. Quesada added external choice for superimposed Coulomb barrier
     42//                          (if useSICB=true)
     43//
    3344
    3445#include "G4PreCompoundProton.hh"
    3546
    36   G4ReactionProduct * G4PreCompoundProton::GetReactionProduct() const
    37   {
    38     G4ReactionProduct * theReactionProduct =
    39       new G4ReactionProduct(G4Proton::ProtonDefinition());
    40     theReactionProduct->SetMomentum(GetMomentum().vect());
    41     theReactionProduct->SetTotalEnergy(GetMomentum().e());
     47G4ReactionProduct * G4PreCompoundProton::GetReactionProduct() const
     48{
     49  G4ReactionProduct * theReactionProduct =
     50    new G4ReactionProduct(G4Proton::ProtonDefinition());
     51  theReactionProduct->SetMomentum(GetMomentum().vect());
     52  theReactionProduct->SetTotalEnergy(GetMomentum().e());
    4253#ifdef PRECOMPOUND_TEST
    43     theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     54  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    4455#endif
    45     return theReactionProduct;
    46   }
    47 
    48 
    49    G4double G4PreCompoundProton::GetRj(const G4int NumberParticles, const G4int NumberCharged)
    50   {
    51     G4double rj = 0.0;
    52     if(NumberParticles > 0) rj = static_cast<G4double>(NumberCharged)/static_cast<G4double>(NumberParticles);
    53     return rj;
    54   }
    55 
    56 
     56  return theReactionProduct;
     57}
     58
     59G4double G4PreCompoundProton::GetRj(const G4int NumberParticles, const G4int NumberCharged)
     60{
     61  G4double rj = 0.0;
     62  if(NumberParticles > 0) rj = static_cast<G4double>(NumberCharged)/static_cast<G4double>(NumberParticles);
     63  return rj;
     64}
    5765
    5866////////////////////////////////////////////////////////////////////////////////////
    … …  
    6371//OPT=3 Kalbach's parameterization
    6472//
    65  G4double G4PreCompoundProton::CrossSection(const  G4double K)
     73G4double G4PreCompoundProton::CrossSection(const  G4double K)
    6674{
    6775  //G4cout<<" In G4PreCompoundProton OPTxs="<<OPTxs<<G4endl;
    … …  
    7583  FragmentA=GetA()+GetRestA();
    7684  FragmentAthrd=std::pow(FragmentA,0.33333);
    77 
    78 
    7985
    8086  if (OPTxs==0) return GetOpt0(K);
    … …  
    94100G4double G4PreCompoundProton::GetOpt0(const  G4double K)
    95101{
    96 const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
     102  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
    97103  // cross section is now given in mb (r0 is in mm) for the sake of consistency
    98104  //with the rest of the options
    99  return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
     105  return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
    100106}
    101107//
    102108//------------
    103109//
    104   G4double G4PreCompoundProton::GetAlpha()
    105   {
    106  G4double aZ = static_cast<G4double>(GetRestZ());
    107     G4double C = 0.0;
    108     if (aZ >= 70)
    109       {
    110         C = 0.10;
    111       }
    112     else
    113       {
    114         C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
    115       }
    116     return 1.0 + C;
    117   }
     110G4double G4PreCompoundProton::GetAlpha()
     111{
     112  G4double aZ = static_cast<G4double>(GetRestZ());
     113  G4double C = 0.0;
     114  if (aZ >= 70)
     115    {
     116      C = 0.10;
     117    }
     118  else
     119    {
     120      C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
     121    }
     122  return 1.0 + C;
     123}
    118124//
    119125//-------------------
    120126// 
    121   G4double G4PreCompoundProton::GetBeta()
    122   {
     127G4double G4PreCompoundProton::GetBeta()
     128{
    123129  return -GetCoulombBarrier();
    124   }
    125 //
    126  
    127 
    128 
     130}
     131//
     132 
    129133//********************* OPT=1 : Chatterjee's cross section ************************
    130134//(fitting to cross section from Bechetti & Greenles OM potential)
    … …  
    132136G4double G4PreCompoundProton::GetOpt1(const  G4double K)
    133137{
    134  G4double Kc=K;
    135 
    136 // JMQ  xsec is set constat above limit of validity
    137  if (K>50)  Kc=50;
    138 
    139  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
    140  G4double p, p0, p1, p2,Ec,delta,q,r,ji;
    141  
    142  p0 = 15.72;
    143  p1 = 9.65;
    144  p2 = -449.0;
    145  landa0 = 0.00437;
    146  landa1 = -16.58;
    147  mu0 = 244.7;
    148  mu1 = 0.503;
    149  nu0 = 273.1;
    150  nu1 = -182.4;
    151  nu2 = -1.872; 
    152  delta=0.; 
    153 
    154  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
    155  p = p0 + p1/Ec + p2/(Ec*Ec);
    156  landa = landa0*ResidualA + landa1;
    157  mu = mu0*std::pow(ResidualA,mu1);
    158  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
    159  q = landa - nu/(Ec*Ec) - 2*p*Ec;
    160  r = mu + 2*nu/Ec + p*(Ec*Ec);
    161 
    162  ji=std::max(Kc,Ec);
    163  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
    164  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
    165  if (xs <0.0) {xs=0.0;}
    166 
    167  return xs;
    168 
    169 }
    170 
    171 
     138  G4double Kc=K;
     139
     140  // JMQ  xsec is set constat above limit of validity
     141  if (K>50)  Kc=50;
     142
     143  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
     144  G4double p, p0, p1, p2,Ec,delta,q,r,ji;
     145 
     146  p0 = 15.72;
     147  p1 = 9.65;
     148  p2 = -449.0;
     149  landa0 = 0.00437;
     150  landa1 = -16.58;
     151  mu0 = 244.7;
     152  mu1 = 0.503;
     153  nu0 = 273.1;
     154  nu1 = -182.4;
     155  nu2 = -1.872; 
     156  delta=0.; 
     157
     158  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     159  p = p0 + p1/Ec + p2/(Ec*Ec);
     160  landa = landa0*ResidualA + landa1;
     161  mu = mu0*std::pow(ResidualA,mu1);
     162  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     163  q = landa - nu/(Ec*Ec) - 2*p*Ec;
     164  r = mu + 2*nu/Ec + p*(Ec*Ec);
     165
     166  ji=std::max(Kc,Ec);
     167  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     168  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     169  if (xs <0.0) {xs=0.0;}
     170
     171  return xs;
     172
     173}
    172174
    173175//************* OPT=2 : Welisch's proton reaction cross section ************************
    … …  
    178180  G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
    179181 
    180 //This is redundant when the Coulomb  barrier is overimposed to all cross sections
    181 //It should be kept when Coulomb barrier only imposed at OPTxs=2
     182  //This is redundant when the Coulomb  barrier is overimposed to all cross sections
     183  //It should be kept when Coulomb barrier only imposed at OPTxs=2
    182184
    183185  if(!useSICB && K<=theCoulombBarrier) return xine_th=0.0;
    … …  
    316318  }   //end for E>Ec
    317319
    318   return sig;}
    319 
    320 
     320  return sig;
     321}
    321322
    322323//   ************************** end of cross sections *******************************
Note: See TracChangeset for help on using the changeset viewer.