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

maj sur la beta de geant 4.9.3

Location:
trunk/source/processes/hadronic/models/pre_equilibrium
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/pre_equilibrium/History

    r1007 r1055  
    1414     * Please list in reverse chronological order (last date on top)
    1515     ---------------------------------------------------------------
     16
     1713-February 2009 V.Ivanchenko   hadr-pre-V09-02-03
     18---------------------------------------------------
     19G4PreCompoundXXX - changed the shape of probabilities (return back to 9.2) 
     20                   for d, t, He3, alpha (JMQ)
     21
     2212-February 2009 V.Ivanchenko   hadr-pre-V09-02-02
     23---------------------------------------------------
     24G4PreCompoundXXX - changed the shape of probabilities
     25                   for d, t, He3, alpha near the Coulomb barrier (JMQ)
     26
     2711-February 2009 V.Ivanchenko   hadr-pre-V09-02-01
     28---------------------------------------------------
     29G4PreCompoundXXX - set default Opt3 back, add decrease Coulomb barrier
     30                   for d, t, a, he3 (JMQ)
     31
     3210-February 2009 V.Ivanchenko   hadr-pre-V09-02-00
     33---------------------------------------------------
     34Some clean up of comments
     35G4PreCompoundIon - fixed probability of light ion emmision (JMQ)
     36G4PreCompoundXXX - by default Opt1 is used for d, t, a, he3,
     37                   Opt3 for n, p (JMQ)
    1638
    173909-December 2008 A.Howard   hadr-pre-V09-01-15

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.icc

    r1007 r1055  
    2424// ********************************************************************
    2525//
     26// $Id: G4PreCompoundEmission.icc,v 1.5 2009/02/10 16:01:37 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29//
     30// Author:         V.Lara
     31//
     32// Modified: 
    2633// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    2734// cross section option
    … …  
    3441  return;
    3542}
    36 
    3743
    3844inline G4double G4PreCompoundEmission::GetTotalProbability(const G4Fragment & aFragment)

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.hh

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4PreCompoundFragmentVector.hh,v 1.4 2008/09/22 10:18:36 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4PreCompoundFragmentVector.hh,v 1.6 2009/02/10 16:01:37 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// Hadronic Process: Nuclear Preequilibrium
    … …  
    3939#define G4PreCompoundFragmentVector_h 1
    4040
    41 
    4241#include "G4VPreCompoundFragment.hh"
    43 
    44 
    4542
    4643class G4PreCompoundFragmentVector

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.icc

    r962 r1055  
    2626//
    2727// $Id: G4PreCompoundFragmentVector.icc,v 1.4 2008/09/22 10:18:36 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// Hadronic Process: Nuclear Preequilibrium

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.hh

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 //J. M. Quesada (August 2008). 
    27 //Based  on previous work by V. Lara
     26//
     27// $Id: G4VPreCompoundFragment.hh,v 1.10 2009/02/10 16:01:37 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     29//
     30// J. M. Quesada (August 2008). 
     31// Based  on previous work by V. Lara
    2832//
    2933// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    … …  
    3135// JMQ (06 September 2008) Also external choice has been added for:
    3236//                      - superimposed Coulomb barrier (if useSICB=true)
    33 
    3437
    3538#ifndef G4VPreCompoundFragment_h

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.icc

    r962 r1055  
    2626//
    2727// $Id: G4VPreCompoundFragment.icc,v 1.7 2008/09/22 10:18:36 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// by V. Lara

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundAlpha.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//
     27// $Id: G4PreCompoundAlpha.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     29//
     30// -------------------------------------------------------------------
     31//
     32// GEANT4 Class file
     33//
     34//
     35// File name:     G4PreCompoundAlpha
     36//
     37// Author:         V.Lara
     38//
     39// Modified: 
     40// 21.08.2008 J. M. Quesada add choice of options 
     41// 10.02.2009 J. M. Quesada set default opt1 
     42//
    3043
    3144#include "G4PreCompoundAlpha.hh"
    3245
    33   G4ReactionProduct * G4PreCompoundAlpha::GetReactionProduct() const
    34   {
    35     G4ReactionProduct * theReactionProduct =
    36       new G4ReactionProduct(G4Alpha::AlphaDefinition());
    37     theReactionProduct->SetMomentum(GetMomentum().vect());
    38     theReactionProduct->SetTotalEnergy(GetMomentum().e());
     46G4ReactionProduct * G4PreCompoundAlpha::GetReactionProduct() const
     47{
     48  G4ReactionProduct * theReactionProduct =
     49    new G4ReactionProduct(G4Alpha::AlphaDefinition());
     50  theReactionProduct->SetMomentum(GetMomentum().vect());
     51  theReactionProduct->SetTotalEnergy(GetMomentum().e());
    3952#ifdef PRECOMPOUND_TEST
    40     theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     53  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    4154#endif
    42     return theReactionProduct;
    43   }   
    44 
    45 
    46    G4double G4PreCompoundAlpha::FactorialFactor(const G4double N, const G4double P)
    47   {
    48      return
     55  return theReactionProduct;
     56}   
     57
     58G4double G4PreCompoundAlpha::FactorialFactor(const G4double N, const G4double P)
     59{
     60  return
    4961      (N-4.0)*(P-3.0)*(
    5062                       (((N-3.0)*(P-2.0))/2.0) *(
    … …  
    5466                                                 )
    5567                       );
     68}
     69 
     70G4double G4PreCompoundAlpha::CoalescenceFactor(const G4double A)
     71{
     72  return 4096.0/(A*A*A); 
     73}   
     74
     75G4double G4PreCompoundAlpha::GetRj(const G4int NumberParticles, const G4int NumberCharged)
     76{
     77  G4double rj = 0.0;
     78  G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2)*(NumberParticles-3);
     79  if(NumberCharged >=2 && (NumberParticles-NumberCharged) >=2 ) {
     80    rj = 6.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged)*
     81                                   (NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); 
    5682  }
    57  
    58    G4double G4PreCompoundAlpha::CoalescenceFactor(const G4double A)
    59   {
    60      return 4096.0/(A*A*A); 
    61   }   
    62 
    63 
    64 
    65    G4double G4PreCompoundAlpha::GetRj(const G4int NumberParticles, const G4int NumberCharged)
    66   {
    67     G4double rj = 0.0;
    68     G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2)*(NumberParticles-3);
    69     if(NumberCharged >=2 && (NumberParticles-NumberCharged) >=2 ) rj = 6.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); 
    70  return rj;
    71   }
    72 
    73 
    74 
     83  return rj;
     84}
    7585
    7686////////////////////////////////////////////////////////////////////////////////////
    … …  
    8090//OPT=3,4 Kalbach's parameterization
    8191//
    82  G4double G4PreCompoundAlpha::CrossSection(const  G4double K)
     92G4double G4PreCompoundAlpha::CrossSection(const  G4double K)
    8393{
    8494
    … …  
    115125//----------------
    116126//
    117  G4double G4PreCompoundAlpha::GetAlpha()
    118   {
    119  G4double C = 0.0;
    120     G4double aZ = GetZ() + GetRestZ();
    121     if (aZ <= 30)
    122       {
    123         C = 0.10;
    124       }
    125     else if (aZ <= 50)
    126       {
    127         C = 0.1 + -((aZ-50.)/20.)*0.02;
    128       }
    129     else if (aZ < 70)
    130       {
    131         C = 0.08 + -((aZ-70.)/20.)*0.02;
    132       }
    133     else
    134       {
    135         C = 0.06;
    136       }
    137     return 1.0+C;
    138   }
     127G4double G4PreCompoundAlpha::GetAlpha()
     128{
     129  G4double C = 0.0;
     130  G4double aZ = GetZ() + GetRestZ();
     131  if (aZ <= 30)
     132    {
     133      C = 0.10;
     134    }
     135  else if (aZ <= 50)
     136    {
     137      C = 0.1 + -((aZ-50.)/20.)*0.02;
     138    }
     139  else if (aZ < 70)
     140    {
     141      C = 0.08 + -((aZ-70.)/20.)*0.02;
     142    }
     143  else
     144    {
     145      C = 0.06;
     146    }
     147  return 1.0+C;
     148}
    139149//
    140150//--------------------
    141151//
    142   G4double G4PreCompoundAlpha::GetBeta()
    143   {
    144       return -GetCoulombBarrier();
    145   }
     152G4double G4PreCompoundAlpha::GetBeta()
     153{
     154  return -GetCoulombBarrier();
     155}
    146156//
    147157//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    153163  G4double Kc=K;
    154164
    155   // JMQ xsec is set constat above limit of validity
     165  // JMQ xsec is set constant above limit of validity
    156166  if (K>50) Kc=50;
    157167
    … …  
    213223  G4double      ra=1.20;
    214224       
    215   ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     225  //JMQ 13/02/09 increase of reduced radius to lower the barrier
     226  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     227  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
    216228  ecsq = ec * ec;
    217229  p = p0 + p1/ec + p2/ecsq;
    … …  
    236248 
    237249  if (elab <= ec) { //start for E<Ec
    238     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
     250    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
    239251  }           //end for E<Ec
    240252  else {           //start for E>Ec
    … …  
    252264
    253265//   ************************** end of cross sections *******************************
    254 
    255 
    256 
    257 

  • 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

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

    r1007 r1055  
    2525//
    2626//
    27 // Hadronic Process: Nuclear Preequilibrium
    28 // by V. Lara
    29 
     27// $Id: G4PreCompoundEmission.cc,v 1.19 2009/02/10 16:01:37 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     29//
     30// -------------------------------------------------------------------
     31//
     32// GEANT4 Class file
     33//
     34//
     35// File name:     G4PreCompoundEmission
     36//
     37// Author:         V.Lara
     38//
     39// Modified: 
     40//
    3041
    3142#include "G4PreCompoundEmission.hh"
    … …  
    4253
    4354
    44   G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
     55G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
    4556{
    4657  return false;
    4758}
    4859
    49     G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
     60G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
    5061{
    5162  return true;
    … …  
    95106  return;
    96107}
    97 
    98 
    99108
    100109G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 //J. M. Quesada (August 2008). 
    27 //Based  on previous work by V. Lara
     26// $Id: G4PreCompoundFragment.cc,v 1.8 2009/02/10 16:01:37 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29// J. M. Quesada (August 2008). 
     30// Based  on previous work by V. Lara
    2831// JMQ (06 September 2008) Also external choice has been added for:
    2932//                      - superimposed Coulomb barrier (if useSICB=true)

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundFragmentVector.cc,v 1.7 2008/05/08 10:36:03 quesada Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PreCompoundFragmentVector.cc,v 1.11 2009/02/10 16:01:37 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2928//
    3029// Hadronic Process: Nuclear Preequilibrium

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundHe3.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//
     27// $Id: G4PreCompoundHe3.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     29//
     30// -------------------------------------------------------------------
     31//
     32// GEANT4 Class file
     33//
     34//
     35// File name:     G4PreCompoundHe3
     36//
     37// Author:         V.Lara
     38//
     39// Modified: 
     40// 21.08.2008 J. M. Quesada add choice of options 
     41// 10.02.2009 J. M. Quesada set default opt1 
     42//
    3043 
    3144#include "G4PreCompoundHe3.hh"
    3245
    33 
    34   G4ReactionProduct * G4PreCompoundHe3::GetReactionProduct() const
    35   {
    36     G4ReactionProduct * theReactionProduct =
    37       new G4ReactionProduct(G4He3::He3Definition());
    38     theReactionProduct->SetMomentum(GetMomentum().vect());
    39     theReactionProduct->SetTotalEnergy(GetMomentum().e());
     46G4ReactionProduct * G4PreCompoundHe3::GetReactionProduct() const
     47{
     48  G4ReactionProduct * theReactionProduct =
     49    new G4ReactionProduct(G4He3::He3Definition());
     50  theReactionProduct->SetMomentum(GetMomentum().vect());
     51  theReactionProduct->SetTotalEnergy(GetMomentum().e());
    4052#ifdef PRECOMPOUND_TEST
    41     theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     53  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    4254#endif
    43     return theReactionProduct;
    44   }   
    45 
    46 
    47    G4double G4PreCompoundHe3::FactorialFactor(const G4double N, const G4double P)
    48   {
    49      return
     55  return theReactionProduct;
     56}   
     57
     58G4double G4PreCompoundHe3::FactorialFactor(const G4double N, const G4double P)
     59{
     60  return
    5061      (N-3.0)*(P-2.0)*(
    5162                       (((N-2.0)*(P-1.0))/2.0) *(
    … …  
    5364                                                 )
    5465                       );
     66}
     67 
     68G4double G4PreCompoundHe3::CoalescenceFactor(const G4double A)
     69{
     70  return 243.0/(A*A);
     71}   
     72
     73G4double G4PreCompoundHe3::GetRj(const G4int NumberParticles, const G4int NumberCharged)
     74{
     75  G4double rj = 0.0;
     76  G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2);
     77  if(NumberCharged >=2 && (NumberParticles-NumberCharged) >= 1) {
     78    rj = 3.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged))
     79      / static_cast<G4double>(denominator); 
    5580  }
    56  
    57   G4double G4PreCompoundHe3::CoalescenceFactor(const G4double A)
    58   {
    59          return 243.0/(A*A);
    60   }   
    61 
    62 
    63 
    64   G4double G4PreCompoundHe3::GetRj(const G4int NumberParticles, const G4int NumberCharged)
    65   {
    66     G4double rj = 0.0;
    67     G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2);
    68     if(NumberCharged >=2 && (NumberParticles-NumberCharged) >= 1) rj = 3.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged))/static_cast<G4double>(denominator); 
    69    
    70     return rj;
    71   }
    72 
    73 
    74 
     81  return rj;
     82}
    7583
    7684////////////////////////////////////////////////////////////////////////////////////
    … …  
    8088//OPT=3,4 Kalbach's parameterization
    8189//
    82  G4double G4PreCompoundHe3::CrossSection(const  G4double K)
     90G4double G4PreCompoundHe3::CrossSection(const  G4double K)
    8391{
    8492  ResidualA=GetRestA();
    … …  
    114122//----------------
    115123//
    116   G4double G4PreCompoundHe3::GetAlpha()
    117   {
    118     G4double C = 0.0;
    119     G4double aZ = GetZ() + GetRestZ();
    120     if (aZ <= 30)
    121       {
    122         C = 0.10;
    123       }
    124     else if (aZ <= 50)
    125       {
    126         C = 0.1 + -((aZ-50.)/20.)*0.02;
    127       }
    128     else if (aZ < 70)
    129       {
    130         C = 0.08 + -((aZ-70.)/20.)*0.02;
    131       }
    132     else
    133       {
    134         C = 0.06;
    135       }
    136     return 1.0 + C*(4.0/3.0);
    137   }
     124G4double G4PreCompoundHe3::GetAlpha()
     125{
     126  G4double C = 0.0;
     127  G4double aZ = GetZ() + GetRestZ();
     128  if (aZ <= 30)
     129    {
     130      C = 0.10;
     131    }
     132  else if (aZ <= 50)
     133    {
     134      C = 0.1 + -((aZ-50.)/20.)*0.02;
     135    }
     136  else if (aZ < 70)
     137    {
     138      C = 0.08 + -((aZ-70.)/20.)*0.02;
     139    }
     140  else
     141    {
     142      C = 0.06;
     143    }
     144  return 1.0 + C*(4.0/3.0);
     145}
    138146//
    139147//--------------------
    140148//
    141    G4double G4PreCompoundHe3::GetBeta()
    142   {
    143       return -GetCoulombBarrier();
    144   }
     149G4double G4PreCompoundHe3::GetBeta()
     150{
     151  return -GetCoulombBarrier();
     152}
    145153//
    146154//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    214222  G4double      ra=0.80;
    215223       
    216   ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     224  //JMQ 13/02/09 increase of reduced radius to lower the barrier
     225  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     226  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
    217227  ecsq = ec * ec;
    218228  p = p0 + p1/ec + p2/ecsq;
    … …  
    237247 
    238248  if (elab <= ec) { //start for E<Ec
    239     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
     249    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
    240250  }           //end for E<Ec
    241251  else {           //start for E>Ec

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 //J. M. Quesada (August 2008). 
    27 //Based  on previous work by V. Lara
     26// $Id: G4PreCompoundIon.cc,v 1.16 2009/02/10 16:01:37 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29// -------------------------------------------------------------------
     30//
     31// GEANT4 Class file
     32//
     33//
     34// File name:     G4PreCompoundIon
     35//
     36// Author:         V.Lara
     37//
     38// Modified: 
     39// 10.02.2009 J. M. Quesada fixed bug in density level of light fragments 
    2840//
    2941
    … …  
    3143#include "G4PreCompoundParameters.hh"
    3244
    33  G4bool G4PreCompoundIon::IsItPossible(const G4Fragment& aFragment)
    34   {
    35     G4int pplus = aFragment.GetNumberOfCharged();   
    36     G4int pneut = aFragment.GetNumberOfParticles()-pplus;
    37     return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
    38   }
     45G4bool G4PreCompoundIon::IsItPossible(const G4Fragment& aFragment)
     46{
     47  G4int pplus = aFragment.GetNumberOfCharged();   
     48  G4int pneut = aFragment.GetNumberOfParticles()-pplus;
     49  return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
     50}
    3951
    4052G4double G4PreCompoundIon::
    … …  
    5769    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    5870
    59   G4double gj = (6.0/pi2)*GetA() *
    60     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     71  //JMQ 06/02/209  This is  THE BUG that was killing cluster emission
     72  // G4double gj = (6.0/pi2)*GetA() *
     73  //   G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     74
     75  G4double gj = g1;
    6176
    6277  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
    … …  
    7489  G4double Ej = std::max(0.0,eKin + GetBindingEnergy() -Aj);
    7590
    76 
    77  G4double pA = 1.e-25*(3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*
    78 (eKin+GetBindingEnergy()))))/(pi * r0 * r0 * std::pow(GetRestA(),2.0/3.0) )* eKin*CrossSection(eKin) /(r0*std::pow(GetRestA(),1.0/3.0)) * CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P)* GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())  ;
     91  // JMQ 10/02/09 reshaping of the formula (unnecessary std::pow elimitated)
     92  G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*
     93                (eKin+GetBindingEnergy()))))/(pi * r0 * r0 *r0* GetRestA())*
     94                eKin*CrossSection(eKin)*millibarn*
     95                CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P)*
     96                GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged());
    7997
    8098  G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0);
    … …  
    82100  G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0;
    83101
    84 
    85102  G4double Probability = pA * pB * pC;
    86103
    87104  return Probability;
    88105}
    89 
    90 
    91 
    92 
    93 
    94 

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 //J. M. Quesada (August 2008). New source file
    27 //
    28 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse
    29 // cross section option
     26//
     27// $Id: G4PreCompoundNeutron.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:     G4PreCompoundNeutron
     36//
     37// Author:         V.Lara
     38//
     39// Modified: 
     40// 21.08.2008 J. M. Quesada add choice of options 
     41// 10.02.2009 J. M. Quesada set default opt3
    3042//
     43
    3144#include "G4PreCompoundNeutron.hh"
    3245
    33 
    34   G4ReactionProduct * G4PreCompoundNeutron::GetReactionProduct() const
    35   {
    36     G4ReactionProduct * theReactionProduct =
    37       new G4ReactionProduct(G4Neutron::NeutronDefinition());
    38     theReactionProduct->SetMomentum(GetMomentum().vect());
    39     theReactionProduct->SetTotalEnergy(GetMomentum().e());
     46G4ReactionProduct * G4PreCompoundNeutron::GetReactionProduct() const
     47{
     48  G4ReactionProduct * theReactionProduct =
     49    new G4ReactionProduct(G4Neutron::NeutronDefinition());
     50  theReactionProduct->SetMomentum(GetMomentum().vect());
     51  theReactionProduct->SetTotalEnergy(GetMomentum().e());
    4052#ifdef PRECOMPOUND_TEST
    41     theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     53  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    4254#endif
    43     return theReactionProduct;
    44   }
    45 
    46 
    47    G4double G4PreCompoundNeutron::GetRj(const G4int NumberParticles, const G4int NumberCharged)
    48   {
    49     G4double rj = 0.0;
    50     if(NumberParticles > 0) rj = static_cast<G4double>(NumberParticles - NumberCharged)/static_cast<G4double>(NumberParticles);
    51     return rj;
    52   }
     55  return theReactionProduct;
     56}
     57
     58G4double G4PreCompoundNeutron::GetRj(const G4int NumberParticles, const G4int NumberCharged)
     59{
     60  G4double rj = 0.0;
     61  if(NumberParticles > 0) rj = static_cast<G4double>(NumberParticles - NumberCharged)/
     62                            static_cast<G4double>(NumberParticles);
     63  return rj;
     64}
    5365
    5466////////////////////////////////////////////////////////////////////////////////////
    … …  
    5870//OPT=3,4 Kalbach's parameterization
    5971//
    60  G4double G4PreCompoundNeutron::CrossSection(const  G4double K)
     72G4double G4PreCompoundNeutron::CrossSection(const  G4double K)
    6173{
    6274  ResidualA=GetRestA();
    … …  
    92104//-------
    93105//
    94   G4double G4PreCompoundNeutron::GetAlpha()
    95   {
    96  //   return 0.76+2.2/std::pow(GetRestA(),1.0/3.0);
     106G4double G4PreCompoundNeutron::GetAlpha()
     107{
     108  //   return 0.76+2.2/std::pow(GetRestA(),1.0/3.0);
    97109  return 0.76+2.2/ResidualAthrd;
    98   }
     110}
    99111//
    100112//------------
    101113//
    102   G4double G4PreCompoundNeutron::GetBeta()
    103   {
    104  //   return (2.12/std::pow(GetRestA(),2.0/3.0)-0.05)*MeV/GetAlpha();
    105  return (2.12/(ResidualAthrd*ResidualAthrd)-0.05)*MeV/GetAlpha();
    106   }
    107 //
    108 
    109  
    110 
     114G4double G4PreCompoundNeutron::GetBeta()
     115{
     116  //   return (2.12/std::pow(GetRestA(),2.0/3.0)-0.05)*MeV/GetAlpha();
     117  return (2.12/(ResidualAthrd*ResidualAthrd)-0.05)*MeV/GetAlpha();
     118}
     119//
    111120
    112121//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    156165  G4double b,ecut,cut,ecut2,geom,elab;
    157166
    158 //safety initialization
     167  //safety initialization
    159168  landa0=0;
    160169  landa1=0;
    … …  
    172181  spill= 1.0e+18;
    173182
    174 
    175 
    176 // PRECO xs for neutrons is choosen
     183  // PRECO xs for neutrons is choosen
    177184
    178185  p0 = -312.;
    … …  
    201208  xnulam = 1.;
    202209  etest = 32.;
    203 //          ** etest is the energy above which the rxn cross section is
    204 //          ** compared with the geometrical limit and the max taken.
    205 //          ** xnulam here is a dummy value to be used later.
    206 
    207 
     210  //          ** etest is the energy above which the rxn cross section is
     211  //          ** compared with the geometrical limit and the max taken.
     212  //          ** xnulam here is a dummy value to be used later.
    208213
    209214  a = -2.*p*ec + landa - nu/ecsq;

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 //J. M. Quesada (August 2008). 
    27 //Based  on previous work by V. Lara
     26//
     27// $Id: G4PreCompoundNucleon.cc,v 1.13 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:     G4PreCompoundNucleon
     36//
     37// Author:         V.Lara
     38//
     39// Modified: 
     40// 10.02.2009 J. M. Quesada cleanup
    2841//
    2942
    … …  
    3245
    3346G4bool G4PreCompoundNucleon::IsItPossible(const G4Fragment& aFragment)
    34   {
    35     G4int pplus = aFragment.GetNumberOfCharged();   
    36     G4int pneut = aFragment.GetNumberOfParticles()-pplus;
    37     return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
    38   }
    39 
     47{
     48  G4int pplus = aFragment.GetNumberOfCharged();   
     49  G4int pneut = aFragment.GetNumberOfParticles()-pplus;
     50  return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
     51}
    4052
    4153G4double G4PreCompoundNucleon::
    … …  
    4557  if ( !IsItPossible(aFragment) ) return 0.0;
    4658
    47 
    4859  G4double U = aFragment.GetExcitationEnergy();
    4960  G4double P = aFragment.GetNumberOfParticles();
    … …  
    5162  G4double N = P + H;
    5263 
    53 
    54 
    5564  G4double g0 = (6.0/pi2)*aFragment.GetA() *
    5665    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    … …  
    5867  G4double g1 = (6.0/pi2)*GetRestA() *
    5968    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    60 
    61 
    6269
    6370  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
    … …  
    7178
    7279
    73  G4double Probability =1.e-25* 1.0/pi2*2.0/(hbarc*hbarc*hbarc) * GetReducedMass()
    74 * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())  * eKin*CrossSection(eKin) * P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 * g1/(g0*g0);
     80  G4double Probability = 1.0/pi2*2.0/(hbarc*hbarc*hbarc) * GetReducedMass()
     81    * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) 
     82    * eKin*CrossSection(eKin)*millibarn * P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 * g1/(g0*g0);
    7583
    76 
    77 if (GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())<0.0 || CrossSection(eKin) <0)
    78 {  std::ostringstream errOs;
    79    G4cout<<"WARNING:  NEGATIVE VALUES "<<G4endl;     
    80         errOs << "Rj=" << GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())<<G4endl;
    81         errOs <<"  xsec("<<eKin<<" MeV) ="<<CrossSection(eKin)<<G4endl;
    82         errOs <<"  A="<<GetA()<<"  Z="<<GetZ()<<G4endl;
    83         throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    84               }
     84  if (GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())<0.0
     85      || CrossSection(eKin) <0) { 
     86    std::ostringstream errOs;
     87    G4cout<<"WARNING:  NEGATIVE VALUES "<<G4endl;     
     88    errOs << "Rj=" << GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())
     89          <<G4endl;
     90    errOs <<"  xsec("<<eKin<<" MeV) ="<<CrossSection(eKin)<<G4endl;
     91    errOs <<"  A="<<GetA()<<"  Z="<<GetZ()<<G4endl;
     92    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     93  }
    8594
    8695  return Probability;

  • 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 *******************************

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 //
    27 //J. M. Quesada (Feb. 08). Base on previous work by V. Lara.
    28 // New transition probabilities. Several bugs fixed.
    29 // JMQ (06 September 2008) Also external choices have been added for:
     26// $Id: G4PreCompoundTransitions.cc,v 1.20 2009/02/10 16:01:37 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29// -------------------------------------------------------------------
     30//
     31// GEANT4 Class file
     32//
     33//
     34// File name:     G4PreCompoundIon
     35//
     36// Author:         V.Lara
     37//
     38// Modified: 
     39// 16.02.2008 J. M. Quesada fixed bugs
     40// 06.09.2008 J. M. Quesada added external choices for:
    3041//                      - "never go back"  hipothesis (useNGB=true)
    3142//                      -  CEM transition probabilities (useCEMtr=true)
     43
    3244#include "G4PreCompoundTransitions.hh"
    3345#include "G4HadronicException.hh"
    … …  
    184196    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
    185197    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl;
    186     return TransitionProb1 + TransitionProb2 + TransitionProb3;}
     198    return TransitionProb1 + TransitionProb2 + TransitionProb3;
     199  }
    187200 
    188201  else  {
    … …  
    215228    return TransitionProb1 + TransitionProb2 + TransitionProb3;
    216229  }
    217 
    218 
    219230}
    220231
    … …  
    237248    }
    238249
    239   // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and in proportion to the number charges w.r.t. number of particles,  PROVIDED that there are charged particles
    240   if(deltaN < 0 && G4UniformRand() <= static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) && (result.GetNumberOfCharged() >= 1))
     250  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and in proportion
     251  // to the number charges w.r.t. number of particles,  PROVIDED that there are charged particles
     252  if(deltaN < 0 && G4UniformRand() <=
     253     static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles())
     254     && (result.GetNumberOfCharged() >= 1)) {
    241255    result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative!
    242 
    243 
    244  //JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on number of charged vs. number of particles fails
     256  }
     257
     258  // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on
     259  // number of charged vs. number of particles fails
    245260  result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
    246261  result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2);
    247262
    248  // With weight Z/A, number of charged particles is increased with +1
     263  // With weight Z/A, number of charged particles is increased with +1
    249264  if ( ( deltaN > 0 ) &&
    250265      (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/
    251                   std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
     266       std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
    252267    {
    253268      result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2);

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

    r1007 r1055  
    2525//
    2626//
    27 //J.M.Quesada (August 08). New source file
    28 //
    29 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse
    30 // cross section option
     27// $Id: G4PreCompoundTriton.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     29//
     30// -------------------------------------------------------------------
     31//
     32// GEANT4 Class file
     33//
     34//
     35// File name:     G4PreCompoundTriton
     36//
     37// Author:         V.Lara
     38//
     39// Modified: 
     40// 21.08.2008 J. M. Quesada add choice of options 
     41// 10.02.2009 J. M. Quesada set default opt1 
     42//
    3143 
    3244#include "G4PreCompoundTriton.hh"
    3345
    3446
    35   G4ReactionProduct * G4PreCompoundTriton::GetReactionProduct() const
    36   {
    37     G4ReactionProduct * theReactionProduct =
    38       new G4ReactionProduct(G4Triton::TritonDefinition());
    39     theReactionProduct->SetMomentum(GetMomentum().vect());
    40     theReactionProduct->SetTotalEnergy(GetMomentum().e());
     47G4ReactionProduct * G4PreCompoundTriton::GetReactionProduct() const
     48{
     49  G4ReactionProduct * theReactionProduct =
     50    new G4ReactionProduct(G4Triton::TritonDefinition());
     51  theReactionProduct->SetMomentum(GetMomentum().vect());
     52  theReactionProduct->SetTotalEnergy(GetMomentum().e());
    4153#ifdef PRECOMPOUND_TEST
    42     theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     54  theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    4355#endif
    44     return theReactionProduct;
    45   }   
    46 
    47    G4double G4PreCompoundTriton::FactorialFactor(const G4double N, const G4double P)
    48   {
    49       return
     56  return theReactionProduct;
     57}   
     58
     59G4double G4PreCompoundTriton::FactorialFactor(const G4double N, const G4double P)
     60{
     61  return
    5062      (N-3.0)*(P-2.0)*(
    5163                       (((N-2.0)*(P-1.0))/2.0) *(
    … …  
    5365                                                 )
    5466                       );
     67}
     68 
     69G4double G4PreCompoundTriton::CoalescenceFactor(const G4double A)
     70{
     71  return 243.0/(A*A);
     72}   
     73
     74G4double G4PreCompoundTriton::GetRj(const G4int NumberParticles, const G4int NumberCharged)
     75{
     76  G4double rj = 0.0;
     77  G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2);
     78  if(NumberCharged >= 1 && (NumberParticles-NumberCharged) >= 2) {
     79    rj = 3.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))
     80      /static_cast<G4double>(denominator);
    5581  }
    56  
    57    G4double G4PreCompoundTriton::CoalescenceFactor(const G4double A)
    58   {
    59      return 243.0/(A*A);
    60   }   
    61 
    62 
    63    G4double G4PreCompoundTriton::GetRj(const G4int NumberParticles, const G4int NumberCharged)
    64   {
    65     G4double rj = 0.0;
    66     G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2);
    67     if(NumberCharged >= 1 && (NumberParticles-NumberCharged) >= 2) rj = 3.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator);
    68     return rj;
    69   }
    70 
    71 
    72 
     82  return rj;
     83}
    7384
    7485////////////////////////////////////////////////////////////////////////////////////
    … …  
    7889//OPT=3,4 Kalbach's parameterization
    7990//
    80  G4double G4PreCompoundTriton::CrossSection(const  G4double K)
     91G4double G4PreCompoundTriton::CrossSection(const  G4double K)
    8192{
    8293  ResidualA=GetRestA();
    … …  
    111122//---------
    112123//
    113   G4double G4PreCompoundTriton::GetAlpha()
    114   {
    115     G4double C = 0.0;
    116     G4double aZ = GetZ() + GetRestZ();
    117     if (aZ >= 70)
    118       {
    119         C = 0.10;
    120       }
    121     else
    122       {
    123         C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
    124       }
    125  
    126     return 1.0 + C/3.0;
    127   }
     124G4double G4PreCompoundTriton::GetAlpha()
     125{
     126  G4double C = 0.0;
     127  G4double aZ = GetZ() + GetRestZ();
     128  if (aZ >= 70)
     129    {
     130      C = 0.10;
     131    }
     132  else
     133    {
     134      C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
     135    }
     136 
     137  return 1.0 + C/3.0;
     138}
    128139//
    129140//-------------
    130141//
    131    G4double G4PreCompoundTriton::GetBeta()
    132   {
    133       return -GetCoulombBarrier();
    134   }
     142G4double G4PreCompoundTriton::GetBeta()
     143{
     144  return -GetCoulombBarrier();
     145}
    135146//
    136147//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    183194
    184195  G4double landa, mu, nu, p , signor(1.),sig;
    185 G4double ec,ecsq,xnulam,etest(0.),a;
    186 G4double b,ecut,cut,ecut2,geom,elab;
    187 
    188 
    189  G4double     flow = 1.e-18;
    190  G4double     spill= 1.e+18;
    191 
    192 
    193  G4double     p0 = -21.45;
    194  G4double     p1 = 484.7;
    195  G4double     p2 = -1608.;
    196  G4double     landa0 = 0.0186;
    197  G4double     landa1 = -8.90;
    198  G4double     mu0 = 686.3;
    199  G4double     mu1 = 0.325;
    200  G4double     nu0 = 368.9;
    201  G4double     nu1 = -522.2;
    202  G4double     nu2 = -4.998; 
     196  G4double ec,ecsq,xnulam,etest(0.),a;
     197  G4double b,ecut,cut,ecut2,geom,elab;
     198
     199
     200  G4double     flow = 1.e-18;
     201  G4double     spill= 1.e+18;
     202
     203
     204  G4double     p0 = -21.45;
     205  G4double     p1 = 484.7;
     206  G4double     p2 = -1608.;
     207  G4double     landa0 = 0.0186;
     208  G4double     landa1 = -8.90;
     209  G4double     mu0 = 686.3;
     210  G4double     mu1 = 0.325;
     211  G4double     nu0 = 368.9;
     212  G4double     nu1 = -522.2;
     213  G4double     nu2 = -4.998; 
    203214 
    204  G4double      ra=0.80;
     215  G4double      ra=0.80;
    205216       
    206  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    207  ecsq = ec * ec;
    208  p = p0 + p1/ec + p2/ecsq;
    209  landa = landa0*ResidualA + landa1;
    210  a = std::pow(ResidualA,mu1);
    211  mu = mu0 * a;
    212  nu = a* (nu0+nu1*ec+nu2*ecsq); 
    213  xnulam = nu / landa;
    214  if (xnulam > spill) xnulam=0.;
    215  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
    216  
    217  a = -2.*p*ec + landa - nu/ecsq;
    218  b = p*ecsq + mu + 2.*nu/ec;
    219  ecut = 0.;
    220  cut = a*a - 4.*p*b;
    221  if (cut > 0.) ecut = std::sqrt(cut);
    222  ecut = (ecut-a) / (p+p);
    223  ecut2 = ecut;
    224  if (cut < 0.) ecut2 = ecut - 2.;
    225  elab = K * FragmentA / ResidualA;
    226  sig = 0.;
    227  
    228  if (elab <= ec) { //start for E<Ec
    229    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
    230  }           //end for E<Ec
    231  else {           //start for E>Ec
    232    sig = (landa*elab+mu+nu/elab) * signor;
    233    geom = 0.;
    234    if (xnulam < flow || elab < etest) return sig;
    235    geom = std::sqrt(theA*K);
    236    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
    237    geom = 31.416 * geom * geom;
    238    sig = std::max(geom,sig);
    239  }           //end for E>Ec
    240  return sig;
     217  //JMQ 13/02/09 increase of reduced radius to lower the barrier
     218  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     219  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
     220  ecsq = ec * ec;
     221  p = p0 + p1/ec + p2/ecsq;
     222  landa = landa0*ResidualA + landa1;
     223  a = std::pow(ResidualA,mu1);
     224  mu = mu0 * a;
     225  nu = a* (nu0+nu1*ec+nu2*ecsq); 
     226  xnulam = nu / landa;
     227  if (xnulam > spill) xnulam=0.;
     228  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
     229 
     230  a = -2.*p*ec + landa - nu/ecsq;
     231  b = p*ecsq + mu + 2.*nu/ec;
     232  ecut = 0.;
     233  cut = a*a - 4.*p*b;
     234  if (cut > 0.) ecut = std::sqrt(cut);
     235  ecut = (ecut-a) / (p+p);
     236  ecut2 = ecut;
     237  if (cut < 0.) ecut2 = ecut - 2.;
     238  elab = K * FragmentA / ResidualA;
     239  sig = 0.;
     240 
     241  if (elab <= ec) { //start for E<Ec
     242    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     243  }           //end for E<Ec
     244  else {           //start for E>Ec
     245    sig = (landa*elab+mu+nu/elab) * signor;
     246    geom = 0.;
     247    if (xnulam < flow || elab < etest) return sig;
     248    geom = std::sqrt(theA*K);
     249    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
     250    geom = 31.416 * geom * geom;
     251    sig = std::max(geom,sig);
     252  }           //end for E>Ec
     253  return sig;
    241254
    242255}
    243256
    244257//   ************************** end of cross sections *******************************
    245 
    246 
    247 
    248 

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
     26// $Id: G4VPreCompoundFragment.cc,v 1.12 2009/02/10 16:01:37 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2628//
    27 //J. M. Quesada (August 2008). 
    28 //Based  on previous work by V. Lara
     29// J. M. Quesada (August 2008). 
     30// Based  on previous work by V. Lara
    2931//
    3032 
Note: See TracChangeset for help on using the changeset viewer.