Ignore:
Timestamp:
Apr 20, 2009, 5:54:05 PM (15 years ago)
Author:
garnier
Message:

update to geant4.9.2

Location:
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src
Files:
17 edited

Legend:

Unmodified
Added
Removed

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

    r968 r1007  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundAlpha.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    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 //
     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
    4330
    4431#include "G4PreCompoundAlpha.hh"
    4532
    46 G4ReactionProduct * G4PreCompoundAlpha::GetReactionProduct() const
    47 {
    48   G4ReactionProduct * theReactionProduct =
    49     new G4ReactionProduct(G4Alpha::AlphaDefinition());
    50   theReactionProduct->SetMomentum(GetMomentum().vect());
    51   theReactionProduct->SetTotalEnergy(GetMomentum().e());
     33  G4ReactionProduct * G4PreCompoundAlpha::GetReactionProduct() const
     34  {
     35    G4ReactionProduct * theReactionProduct =
     36      new G4ReactionProduct(G4Alpha::AlphaDefinition());
     37    theReactionProduct->SetMomentum(GetMomentum().vect());
     38    theReactionProduct->SetTotalEnergy(GetMomentum().e());
    5239#ifdef PRECOMPOUND_TEST
    53   theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     40    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    5441#endif
    55   return theReactionProduct;
    56 }   
    57 
    58 G4double G4PreCompoundAlpha::FactorialFactor(const G4double N, const G4double P)
    59 {
    60   return
     42    return theReactionProduct;
     43  }   
     44
     45
     46   G4double G4PreCompoundAlpha::FactorialFactor(const G4double N, const G4double P)
     47  {
     48     return
    6149      (N-4.0)*(P-3.0)*(
    6250                       (((N-3.0)*(P-2.0))/2.0) *(
    … …  
    6654                                                 )
    6755                       );
    68 }
    69  
    70 G4double G4PreCompoundAlpha::CoalescenceFactor(const G4double A)
    71 {
    72   return 4096.0/(A*A*A); 
    73 }   
    74 
    75 G4double 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); 
    82   }
    83   return rj;
    84 }
     56  }
     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
    8575
    8676////////////////////////////////////////////////////////////////////////////////////
    … …  
    9080//OPT=3,4 Kalbach's parameterization
    9181//
    92 G4double G4PreCompoundAlpha::CrossSection(const  G4double K)
     82 G4double G4PreCompoundAlpha::CrossSection(const  G4double K)
    9383{
    9484
    … …  
    125115//----------------
    126116//
    127 G4double 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 }
     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  }
    149139//
    150140//--------------------
    151141//
    152 G4double G4PreCompoundAlpha::GetBeta()
    153 {
    154   return -GetCoulombBarrier();
    155 }
     142  G4double G4PreCompoundAlpha::GetBeta()
     143  {
     144      return -GetCoulombBarrier();
     145  }
    156146//
    157147//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    163153  G4double Kc=K;
    164154
    165   // JMQ xsec is set constant above limit of validity
     155  // JMQ xsec is set constat above limit of validity
    166156  if (K>50) Kc=50;
    167157
    … …  
    223213  G4double      ra=1.20;
    224214       
    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);
     215  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    228216  ecsq = ec * ec;
    229217  p = p0 + p1/ec + p2/ecsq;
    … …  
    248236 
    249237  if (elab <= ec) { //start for E<Ec
    250     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     238    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
    251239  }           //end for E<Ec
    252240  else {           //start for E>Ec
    … …  
    264252
    265253//   ************************** end of cross sections *******************************
     254
     255
     256
     257

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

    r968 r1007  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PreCompoundDeuteron.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    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 //
     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
    4230
    4331#include "G4PreCompoundDeuteron.hh"
    4432
    45 G4ReactionProduct * G4PreCompoundDeuteron::GetReactionProduct() const
    46 {
    47   G4ReactionProduct * theReactionProduct =
    48     new G4ReactionProduct(G4Deuteron::DeuteronDefinition());
    49   theReactionProduct->SetMomentum(GetMomentum().vect());
    50   theReactionProduct->SetTotalEnergy(GetMomentum().e());
     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());
    5141#ifdef PRECOMPOUND_TEST
    52   theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     42    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    5343#endif
    54   return theReactionProduct;
    55 }   
    56 
    57  
    58 G4double G4PreCompoundDeuteron::FactorialFactor(const G4double N, const G4double P)
    59 {
    60   return (N-1.0)*(N-2.0)*(P-1.0)*P/2.0;
    61 }
     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;
     51  }
    6252 
    63 G4double G4PreCompoundDeuteron::CoalescenceFactor(const G4double A)
    64 {
    65   return 16.0/A;
    66 }   
    67 
    68 G4double 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);
    75   }
    76   return rj;
    77 }
     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  }
    7866
    7967////////////////////////////////////////////////////////////////////////////////////
    … …  
    8371//OPT=3,4 Kalbach's parameterization
    8472//
    85 G4double G4PreCompoundDeuteron::CrossSection(const  G4double K)
     73 G4double G4PreCompoundDeuteron::CrossSection(const  G4double K)
    8674{
    8775
    … …  
    118106//---------
    119107//
    120 G4double 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 }
     108 G4double G4PreCompoundDeuteron::GetAlpha()
     109  {
     110G4double 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  }
    134122//
    135123//---------
    136124//
    137 G4double G4PreCompoundDeuteron::GetBeta()
    138 {
    139   return -GetCoulombBarrier();
    140 }
     125  G4double G4PreCompoundDeuteron::GetBeta()
     126  {
     127      return -GetCoulombBarrier();
     128  }
    141129//
    142130//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    152140
    153141  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    154   //G4double Eo(0),epsilon1(0),epsilon2(0),discri(0);
     142//G4double Eo(0),epsilon1(0),epsilon2(0),discri(0);
    155143
    156144 
    … …  
    186174}
    187175
     176
     177
     178
    188179// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
    189180G4double G4PreCompoundDeuteron::GetOpt34(const  G4double K)
    … …  
    214205  G4double     ra=0.80;
    215206       
    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);
     207  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    219208  ecsq = ec * ec;
    220209  p = p0 + p1/ec + p2/ecsq;
    … …  
    237226  elab = K * FragmentA / ResidualA;
    238227  sig = 0.;
    239 
     228 
    240229  if (elab <= ec) { //start for E<Ec
    241     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     230    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
    242231  }           //end for E<Ec
    243232  else {           //start for E>Ec

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

    r962 r1007  
    2525//
    2626//
    27 // $Id: G4PreCompoundEmission.cc,v 1.19 2009/02/10 16:01:37 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    29 //
    30 // -------------------------------------------------------------------
    31 //
    32 // GEANT4 Class file
    33 //
    34 //
    35 // File name:     G4PreCompoundEmission
    36 //
    37 // Author:         V.Lara
    38 //
    39 // Modified: 
    40 //
     27// Hadronic Process: Nuclear Preequilibrium
     28// by V. Lara
     29
    4130
    4231#include "G4PreCompoundEmission.hh"
    … …  
    5342
    5443
    55 G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
     44  G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
    5645{
    5746  return false;
    5847}
    5948
    60 G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
     49    G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
    6150{
    6251  return true;
    … …  
    10695  return;
    10796}
     97
     98
    10899
    109100G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)

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

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

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

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

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

    r968 r1007  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundHe3.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    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 //
     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
    4330 
    4431#include "G4PreCompoundHe3.hh"
    4532
    46 G4ReactionProduct * G4PreCompoundHe3::GetReactionProduct() const
    47 {
    48   G4ReactionProduct * theReactionProduct =
    49     new G4ReactionProduct(G4He3::He3Definition());
    50   theReactionProduct->SetMomentum(GetMomentum().vect());
    51   theReactionProduct->SetTotalEnergy(GetMomentum().e());
     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());
    5240#ifdef PRECOMPOUND_TEST
    53   theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     41    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    5442#endif
    55   return theReactionProduct;
    56 }   
    57 
    58 G4double G4PreCompoundHe3::FactorialFactor(const G4double N, const G4double P)
    59 {
    60   return
     43    return theReactionProduct;
     44  }   
     45
     46
     47   G4double G4PreCompoundHe3::FactorialFactor(const G4double N, const G4double P)
     48  {
     49     return
    6150      (N-3.0)*(P-2.0)*(
    6251                       (((N-2.0)*(P-1.0))/2.0) *(
    … …  
    6453                                                 )
    6554                       );
    66 }
    67  
    68 G4double G4PreCompoundHe3::CoalescenceFactor(const G4double A)
    69 {
    70   return 243.0/(A*A);
    71 }   
    72 
    73 G4double 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); 
    80   }
    81   return rj;
    82 }
     55  }
     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
    8375
    8476////////////////////////////////////////////////////////////////////////////////////
    … …  
    8880//OPT=3,4 Kalbach's parameterization
    8981//
    90 G4double G4PreCompoundHe3::CrossSection(const  G4double K)
     82 G4double G4PreCompoundHe3::CrossSection(const  G4double K)
    9183{
    9284  ResidualA=GetRestA();
    … …  
    122114//----------------
    123115//
    124 G4double 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 }
     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  }
    146138//
    147139//--------------------
    148140//
    149 G4double G4PreCompoundHe3::GetBeta()
    150 {
    151   return -GetCoulombBarrier();
    152 }
     141   G4double G4PreCompoundHe3::GetBeta()
     142  {
     143      return -GetCoulombBarrier();
     144  }
    153145//
    154146//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    222214  G4double      ra=0.80;
    223215       
    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);
     216  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    227217  ecsq = ec * ec;
    228218  p = p0 + p1/ec + p2/ecsq;
    … …  
    247237 
    248238  if (elab <= ec) { //start for E<Ec
    249     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     239    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
    250240  }           //end for E<Ec
    251241  else {           //start for E>Ec

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

    r962 r1007  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PreCompoundIon.cc,v 1.16 2009/02/10 16:01:37 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    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 
     26//J. M. Quesada (August 2008). 
     27//Based  on previous work by V. Lara
    4028//
    4129
    … …  
    4331#include "G4PreCompoundParameters.hh"
    4432
    45 G4bool 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 }
     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  }
    5139
    5240G4double G4PreCompoundIon::
    … …  
    6957    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    7058
    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;
     59  G4double gj = (6.0/pi2)*GetA() *
     60    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    7661
    7762  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
    … …  
    8974  G4double Ej = std::max(0.0,eKin + GetBindingEnergy() -Aj);
    9075
    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());
     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())  ;
    9779
    9880  G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0);
    … …  
    10082  G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0;
    10183
     84
    10285  G4double Probability = pA * pB * pC;
    10386
    10487  return Probability;
    10588}
     89
     90
     91
     92
     93
     94

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

    r962 r1007  
    2626//
    2727// $Id: G4PreCompoundModel.cc,v 1.17 2008/12/09 14:09:59 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-02 $
    2929//
    3030// by V. Lara

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

    r968 r1007  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundNeutron.cc,v 1.4 2009/02/11 18:06:00 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    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
     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
    4230//
    43 
    4431#include "G4PreCompoundNeutron.hh"
    4532
    46 G4ReactionProduct * G4PreCompoundNeutron::GetReactionProduct() const
    47 {
    48   G4ReactionProduct * theReactionProduct =
    49     new G4ReactionProduct(G4Neutron::NeutronDefinition());
    50   theReactionProduct->SetMomentum(GetMomentum().vect());
    51   theReactionProduct->SetTotalEnergy(GetMomentum().e());
     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());
    5240#ifdef PRECOMPOUND_TEST
    53   theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     41    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    5442#endif
    55   return theReactionProduct;
    56 }
    57 
    58 G4double 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 }
     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  }
    6553
    6654////////////////////////////////////////////////////////////////////////////////////
    … …  
    7058//OPT=3,4 Kalbach's parameterization
    7159//
    72 G4double G4PreCompoundNeutron::CrossSection(const  G4double K)
     60 G4double G4PreCompoundNeutron::CrossSection(const  G4double K)
    7361{
    7462  ResidualA=GetRestA();
    … …  
    10492//-------
    10593//
    106 G4double G4PreCompoundNeutron::GetAlpha()
    107 {
    108   //   return 0.76+2.2/std::pow(GetRestA(),1.0/3.0);
     94  G4double G4PreCompoundNeutron::GetAlpha()
     95  {
     96 //   return 0.76+2.2/std::pow(GetRestA(),1.0/3.0);
    10997  return 0.76+2.2/ResidualAthrd;
    110 }
     98  }
    11199//
    112100//------------
    113101//
    114 G4double 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 //
     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
    120111
    121112//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    165156  G4double b,ecut,cut,ecut2,geom,elab;
    166157
    167   //safety initialization
     158//safety initialization
    168159  landa0=0;
    169160  landa1=0;
    … …  
    181172  spill= 1.0e+18;
    182173
    183   // PRECO xs for neutrons is choosen
     174
     175
     176// PRECO xs for neutrons is choosen
    184177
    185178  p0 = -312.;
    … …  
    208201  xnulam = 1.;
    209202  etest = 32.;
    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.
     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
    213208
    214209  a = -2.*p*ec + landa - nu/ecsq;

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

    r962 r1007  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundNucleon.cc,v 1.13 2009/02/11 18:06:00 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    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
     26//J. M. Quesada (August 2008). 
     27//Based  on previous work by V. Lara
    4128//
    4229
    … …  
    4532
    4633G4bool G4PreCompoundNucleon::IsItPossible(const G4Fragment& aFragment)
    47 {
    48   G4int pplus = aFragment.GetNumberOfCharged();   
    49   G4int pneut = aFragment.GetNumberOfParticles()-pplus;
    50   return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
    51 }
     34  {
     35    G4int pplus = aFragment.GetNumberOfCharged();   
     36    G4int pneut = aFragment.GetNumberOfParticles()-pplus;
     37    return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
     38  }
     39
    5240
    5341G4double G4PreCompoundNucleon::
    … …  
    5745  if ( !IsItPossible(aFragment) ) return 0.0;
    5846
     47
    5948  G4double U = aFragment.GetExcitationEnergy();
    6049  G4double P = aFragment.GetNumberOfParticles();
    … …  
    6251  G4double N = P + H;
    6352 
     53
     54
    6455  G4double g0 = (6.0/pi2)*aFragment.GetA() *
    6556    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    … …  
    6758  G4double g1 = (6.0/pi2)*GetRestA() *
    6859    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     60
     61
    6962
    7063  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
    … …  
    7871
    7972
    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);
     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);
    8375
    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   }
     76
     77if (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              }
    9485
    9586  return Probability;

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

    r962 r1007  
    2626//
    2727// $Id: G4PreCompoundParameters.cc,v 1.3 2006/06/29 20:59:29 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-02 $
    2929//
    3030// by V. Lara

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

    r968 r1007  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundProton.cc,v 1.4 2009/02/11 18:06:00 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    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 //
     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)
    4433
    4534#include "G4PreCompoundProton.hh"
    4635
    47 G4ReactionProduct * G4PreCompoundProton::GetReactionProduct() const
    48 {
    49   G4ReactionProduct * theReactionProduct =
    50     new G4ReactionProduct(G4Proton::ProtonDefinition());
    51   theReactionProduct->SetMomentum(GetMomentum().vect());
    52   theReactionProduct->SetTotalEnergy(GetMomentum().e());
     36  G4ReactionProduct * G4PreCompoundProton::GetReactionProduct() const
     37  {
     38    G4ReactionProduct * theReactionProduct =
     39      new G4ReactionProduct(G4Proton::ProtonDefinition());
     40    theReactionProduct->SetMomentum(GetMomentum().vect());
     41    theReactionProduct->SetTotalEnergy(GetMomentum().e());
    5342#ifdef PRECOMPOUND_TEST
    54   theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     43    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    5544#endif
    56   return theReactionProduct;
    57 }
    58 
    59 G4double 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 }
     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
    6557
    6658////////////////////////////////////////////////////////////////////////////////////
    … …  
    7163//OPT=3 Kalbach's parameterization
    7264//
    73 G4double G4PreCompoundProton::CrossSection(const  G4double K)
     65 G4double G4PreCompoundProton::CrossSection(const  G4double K)
    7466{
    7567  //G4cout<<" In G4PreCompoundProton OPTxs="<<OPTxs<<G4endl;
    … …  
    8375  FragmentA=GetA()+GetRestA();
    8476  FragmentAthrd=std::pow(FragmentA,0.33333);
     77
     78
    8579
    8680  if (OPTxs==0) return GetOpt0(K);
    … …  
    10094G4double G4PreCompoundProton::GetOpt0(const  G4double K)
    10195{
    102   const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
     96const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
    10397  // cross section is now given in mb (r0 is in mm) for the sake of consistency
    10498  //with the rest of the options
    105   return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
     99 return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K);
    106100}
    107101//
    108102//------------
    109103//
    110 G4double 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 }
     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  }
    124118//
    125119//-------------------
    126120// 
    127 G4double G4PreCompoundProton::GetBeta()
    128 {
     121  G4double G4PreCompoundProton::GetBeta()
     122  {
    129123  return -GetCoulombBarrier();
    130 }
    131 //
    132  
     124  }
     125//
     126 
     127
     128
    133129//********************* OPT=1 : Chatterjee's cross section ************************
    134130//(fitting to cross section from Bechetti & Greenles OM potential)
    … …  
    136132G4double G4PreCompoundProton::GetOpt1(const  G4double K)
    137133{
    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;
     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;
    172168
    173169}
     170
     171
    174172
    175173//************* OPT=2 : Welisch's proton reaction cross section ************************
    … …  
    180178  G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
    181179 
    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
     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
    184182
    185183  if(!useSICB && K<=theCoulombBarrier) return xine_th=0.0;
    … …  
    318316  }   //end for E>Ec
    319317
    320   return sig;
    321 }
     318  return sig;}
     319
     320
    322321
    323322//   ************************** end of cross sections *******************************

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

    r962 r1007  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PreCompoundTransitions.cc,v 1.20 2009/02/10 16:01:37 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    2826//
    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:
     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:
    4130//                      - "never go back"  hipothesis (useNGB=true)
    4231//                      -  CEM transition probabilities (useCEMtr=true)
    43 
    4432#include "G4PreCompoundTransitions.hh"
    4533#include "G4HadronicException.hh"
    … …  
    196184    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
    197185    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl;
    198     return TransitionProb1 + TransitionProb2 + TransitionProb3;
    199   }
     186    return TransitionProb1 + TransitionProb2 + TransitionProb3;}
    200187 
    201188  else  {
    … …  
    228215    return TransitionProb1 + TransitionProb2 + TransitionProb3;
    229216  }
     217
     218
    230219}
    231220
    … …  
    248237    }
    249238
    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)) {
     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))
    255241    result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative!
    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
     242
     243
     244 //JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on number of charged vs. number of particles fails
    260245  result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
    261246  result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2);
    262247
    263   // With weight Z/A, number of charged particles is increased with +1
     248 // With weight Z/A, number of charged particles is increased with +1
    264249  if ( ( deltaN > 0 ) &&
    265250      (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/
    266        std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
     251                  std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
    267252    {
    268253      result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2);

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

    r968 r1007  
    2525//
    2626//
    27 // $Id: G4PreCompoundTriton.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    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 //
     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
    4331 
    4432#include "G4PreCompoundTriton.hh"
    4533
    4634
    47 G4ReactionProduct * G4PreCompoundTriton::GetReactionProduct() const
    48 {
    49   G4ReactionProduct * theReactionProduct =
    50     new G4ReactionProduct(G4Triton::TritonDefinition());
    51   theReactionProduct->SetMomentum(GetMomentum().vect());
    52   theReactionProduct->SetTotalEnergy(GetMomentum().e());
     35  G4ReactionProduct * G4PreCompoundTriton::GetReactionProduct() const
     36  {
     37    G4ReactionProduct * theReactionProduct =
     38      new G4ReactionProduct(G4Triton::TritonDefinition());
     39    theReactionProduct->SetMomentum(GetMomentum().vect());
     40    theReactionProduct->SetTotalEnergy(GetMomentum().e());
    5341#ifdef PRECOMPOUND_TEST
    54   theReactionProduct->SetCreatorModel("G4PrecompoundModel");
     42    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
    5543#endif
    56   return theReactionProduct;
    57 }   
    58 
    59 G4double G4PreCompoundTriton::FactorialFactor(const G4double N, const G4double P)
    60 {
    61   return
     44    return theReactionProduct;
     45  }   
     46
     47   G4double G4PreCompoundTriton::FactorialFactor(const G4double N, const G4double P)
     48  {
     49      return
    6250      (N-3.0)*(P-2.0)*(
    6351                       (((N-2.0)*(P-1.0))/2.0) *(
    … …  
    6553                                                 )
    6654                       );
    67 }
     55  }
    6856 
    69 G4double G4PreCompoundTriton::CoalescenceFactor(const G4double A)
    70 {
    71   return 243.0/(A*A);
    72 }   
    73 
    74 G4double 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);
    81   }
    82   return rj;
    83 }
     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
    8473
    8574////////////////////////////////////////////////////////////////////////////////////
    … …  
    8978//OPT=3,4 Kalbach's parameterization
    9079//
    91 G4double G4PreCompoundTriton::CrossSection(const  G4double K)
     80 G4double G4PreCompoundTriton::CrossSection(const  G4double K)
    9281{
    9382  ResidualA=GetRestA();
    … …  
    122111//---------
    123112//
    124 G4double 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 }
     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  }
    139128//
    140129//-------------
    141130//
    142 G4double G4PreCompoundTriton::GetBeta()
    143 {
    144   return -GetCoulombBarrier();
    145 }
     131   G4double G4PreCompoundTriton::GetBeta()
     132  {
     133      return -GetCoulombBarrier();
     134  }
    146135//
    147136//********************* OPT=1,2 : Chatterjee's cross section ************************
    … …  
    194183
    195184  G4double landa, mu, nu, p , signor(1.),sig;
    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; 
     185G4double ec,ecsq,xnulam,etest(0.),a;
     186G4double 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; 
    214203 
    215   G4double      ra=0.80;
     204 G4double      ra=0.80;
    216205       
    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;
     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;
    254241
    255242}
    256243
    257244//   ************************** end of cross sections *******************************
     245
     246
     247
     248

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

    r962 r1007  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VPreCompoundFragment.cc,v 1.12 2009/02/10 16:01:37 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02-ref-02 $
    2826//
    29 // J. M. Quesada (August 2008). 
    30 // Based  on previous work by V. Lara
     27//J. M. Quesada (August 2008). 
     28//Based  on previous work by V. Lara
    3129//
    3230 

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

    r962 r1007  
    2626//
    2727// $Id: G4VPreCompoundIon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-02 $
    2929//
    3030// by V. Lara

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

    r962 r1007  
    2727
    2828// $Id: G4VPreCompoundNucleon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $
    29 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     29// GEANT4 tag $Name: geant4-09-02 $
    3030
    3131//
Note: See TracChangeset for help on using the changeset viewer.