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/G4PreCompoundDeuteron.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
     26// $Id: G4PreCompoundDeuteron.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29// -------------------------------------------------------------------
     30//
     31// GEANT4 Class file
     32//
     33//
     34// File name:     G4PreCompoundDeuteron
     35//
     36// Author:         V.Lara
     37//
     38// Modified: 
     39// 21.08.2008 J. M. Quesada add choice of options 
     40// 10.02.2009 J. M. Quesada set default opt1 
     41//
    3042
    3143#include "G4PreCompoundDeuteron.hh"
    3244
    33 
    34 
    35   G4ReactionProduct * G4PreCompoundDeuteron::GetReactionProduct() const
    36   {
    37     G4ReactionProduct * theReactionProduct =
    38       new G4ReactionProduct(G4Deuteron::DeuteronDefinition());
    39     theReactionProduct->SetMomentum(GetMomentum().vect());
    40     theReactionProduct->SetTotalEnergy(GetMomentum().e());
     45G4ReactionProduct * G4PreCompoundDeuteron::GetReactionProduct() const
     46{
     47  G4ReactionProduct * theReactionProduct =
     48    new G4ReactionProduct(G4Deuteron::DeuteronDefinition());
     49  theReactionProduct->SetMomentum(GetMomentum().vect());
     50  theReactionProduct->SetTotalEnergy(GetMomentum().e());
    4151#ifdef PRECOMPOUND_TEST
    42     theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     52  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    4353#endif
    44     return theReactionProduct;
    45   }   
    46 
    47  
    48    G4double G4PreCompoundDeuteron::FactorialFactor(const G4double N, const G4double P)
    49   {
    50       return (N-1.0)*(N-2.0)*(P-1.0)*P/2.0;
     54  return theReactionProduct;
     55}   
     56
     57 
     58G4double G4PreCompoundDeuteron::FactorialFactor(const G4double N, const G4double P)
     59{
     60  return (N-1.0)*(N-2.0)*(P-1.0)*P/2.0;
     61}
     62 
     63G4double G4PreCompoundDeuteron::CoalescenceFactor(const G4double A)
     64{
     65  return 16.0/A;
     66}   
     67
     68G4double G4PreCompoundDeuteron::GetRj(const G4int NumberParticles, const G4int NumberCharged)
     69{
     70  G4double rj = 0.0;
     71  G4double denominator = NumberParticles*(NumberParticles-1);
     72  if(NumberCharged >=1 && (NumberParticles-NumberCharged) >=1) {
     73    rj = 2.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged))
     74      / static_cast<G4double>(denominator);
    5175  }
    52  
    53    G4double G4PreCompoundDeuteron::CoalescenceFactor(const G4double A)
    54   {
    55     return 16.0/A;
    56   }   
    57 
    58   G4double G4PreCompoundDeuteron::GetRj(const G4int NumberParticles, const G4int NumberCharged)
    59   {
    60     G4double rj = 0.0;
    61     G4double denominator = NumberParticles*(NumberParticles-1);
    62    if(NumberCharged >=1 && (NumberParticles-NumberCharged) >=1) rj = 2.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged))/static_cast<G4double>(denominator);
    63 
    64     return rj;
    65   }
     76  return rj;
     77}
    6678
    6779////////////////////////////////////////////////////////////////////////////////////
    … …  
    7183//OPT=3,4 Kalbach's parameterization
    7284//
    73  G4double G4PreCompoundDeuteron::CrossSection(const  G4double K)
     85G4double G4PreCompoundDeuteron::CrossSection(const  G4double K)
    7486{
    7587
    … …  
    106118//---------
    107119//
    108  G4double G4PreCompoundDeuteron::GetAlpha()
    109   {
    110 G4double C = 0.0;
    111     G4double aZ = GetZ() + GetRestZ();
    112     if (aZ >= 70)
    113       {
    114         C = 0.10;
    115       }
    116     else
    117       {
    118         C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
    119       }
    120     return 1.0 + C/2.0;
    121   }
     120G4double G4PreCompoundDeuteron::GetAlpha()
     121{
     122  G4double C = 0.0;
     123  G4double aZ = GetZ() + GetRestZ();
     124  if (aZ >= 70)
     125    {
     126      C = 0.10;
     127    }
     128  else
     129    {
     130      C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
     131    }
     132  return 1.0 + C/2.0;
     133}
    122134//
    123135//---------
    124136//
    125   G4double G4PreCompoundDeuteron::GetBeta()
    126   {
    127       return -GetCoulombBarrier();
    128   }
     137G4double G4PreCompoundDeuteron::GetBeta()
     138{
     139  return -GetCoulombBarrier();
     140}
    129141//
    130142//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    140152
    141153  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    142 //G4double Eo(0),epsilon1(0),epsilon2(0),discri(0);
     154  //G4double Eo(0),epsilon1(0),epsilon2(0),discri(0);
    143155
    144156 
    … …  
    174186}
    175187
    176 
    177 
    178 
    179188// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
    180189G4double G4PreCompoundDeuteron::GetOpt34(const  G4double K)
    … …  
    205214  G4double     ra=0.80;
    206215       
    207   ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     216  //JMQ 13/02/09 increase of reduced radius to lower the barrier
     217  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     218  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
    208219  ecsq = ec * ec;
    209220  p = p0 + p1/ec + p2/ecsq;
    … …  
    226237  elab = K * FragmentA / ResidualA;
    227238  sig = 0.;
    228  
     239
    229240  if (elab <= ec) { //start for E<Ec
    230     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
     241    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
    231242  }           //end for E<Ec
    232243  else {           //start for E>Ec
Note: See TracChangeset for help on using the changeset viewer.