Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:
trunk/source/processes/hadronic/models/de_excitation/gem_evaporation
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4EvaporationGEMFactory.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4EvaporationGEMFactory.hh,v 1.7 2009/09/15 12:54:16 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4EvaporationGEMFactory.hh,v 1.8 2010/04/27 11:43:16 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    4141{
    4242public:
    43   G4EvaporationGEMFactory() {}
    44   virtual ~G4EvaporationGEMFactory() {}
     43  G4EvaporationGEMFactory();
     44  virtual ~G4EvaporationGEMFactory();
    4545
    4646private:
    47   G4EvaporationGEMFactory(const G4EvaporationGEMFactory & ) : G4VEvaporationFactory() {}
     47  G4EvaporationGEMFactory(const G4EvaporationGEMFactory & );
    4848  const G4EvaporationGEMFactory & operator=(const G4EvaporationGEMFactory & val);
    4949  G4bool operator==(const G4EvaporationGEMFactory & val) const;

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4GEMProbability.hh

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4GEMProbability.hh,v 1.4 2010/05/19 10:21:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
     28//
     29//---------------------------------------------------------------------
     30//
     31// Geant4 header G4GEMProbability
    2632//
    2733//
    … …  
    2935// by V. Lara (Sept 2001)
    3036//
    31 
    32 
     37// 18.05.2010 V.Ivanchenko trying to speedup the most slow method
     38//            by usage of G4Pow, integer Z and A; moved constructor,
     39//            destructor and virtual functions to source
     40//
    3341
    3442#ifndef G4GEMProbability_h
    … …  
    4250#include "G4PairingCorrection.hh"
    4351
     52class G4Pow;
     53
    4454class G4GEMProbability : public G4VEmissionProbability
    4555{
    4656public:
    47     // Only available constructor
    48     G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin) :
    49         theA(anA), theZ(aZ), Spin(aSpin), theCoulombBarrierPtr(0),
    50         ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
    51         {
    52             theEvapLDPptr = new G4EvaporationLevelDensityParameter;
    53         }
     57
     58  // Default constructor - should not be used
     59  G4GEMProbability();
     60
     61  // Only available constructor
     62  G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin);
    5463   
    55     ~G4GEMProbability()
    56         {
    57             if (theEvapLDPptr != 0) delete theEvapLDPptr;
    58         }
     64  virtual ~G4GEMProbability();
    5965
     66  inline G4int GetZ_asInt(void) const { return theZ; }
     67       
     68  inline G4int GetA_asInt(void) const { return theA;}
     69       
     70  inline G4double GetZ(void) const { return theZ; }
     71       
     72  inline G4double GetA(void) const { return theA;}
    6073
    61        
    62     G4double GetZ(void) const { return theZ; }
    63        
    64     G4double GetA(void) const { return theA;}
     74  inline G4double GetSpin(void) const { return Spin; }
    6575
    66     G4double GetSpin(void) const { return Spin; }
     76  inline G4double GetNormalization(void) const { return Normalization; }
     77   
     78  inline void SetCoulomBarrier(const G4VCoulombBarrier * aCoulombBarrierStrategy)
     79  {
     80    theCoulombBarrierPtr = aCoulombBarrierStrategy;
     81  }
    6782
    68     G4double GetNormalization(void) const { return Normalization; }
     83  inline G4double GetCoulombBarrier(const G4Fragment& fragment) const
     84  {
     85    G4double res = 0.0;
     86    if (theCoulombBarrierPtr)
     87      {
     88        G4int Acomp = fragment.GetA_asInt();
     89        G4int Zcomp = fragment.GetZ_asInt();
     90        res = theCoulombBarrierPtr->GetCoulombBarrier(Acomp-theA, Zcomp-theZ,
     91          fragment.GetExcitationEnergy()-fPairCorr->GetPairingCorrection(Acomp,Zcomp));
     92      }
     93    return res;
     94  }
    6995   
    70     void SetCoulomBarrier(const G4VCoulombBarrier * aCoulombBarrierStrategy)
    71         {
    72             theCoulombBarrierPtr = aCoulombBarrierStrategy;
    73         }
    74 
    75     G4double GetCoulombBarrier(const G4Fragment& fragment) const
    76         {
    77             if (theCoulombBarrierPtr)
    78               {
    79                 G4int Acompound = static_cast<G4int>(fragment.GetA());
    80                 G4int Zcompound = static_cast<G4int>(fragment.GetZ());
    81                 return theCoulombBarrierPtr->GetCoulombBarrier(Acompound-theA, Zcompound-theZ,
    82                                                                fragment.GetExcitationEnergy()-
    83                                                                G4PairingCorrection::GetInstance()->
    84                                                                GetPairingCorrection(Acompound,Zcompound));
    85               }
    86             else
    87               {
    88                 return 0.0;
    89               }
    90         }
    91    
    92 
    93     virtual G4double CalcAlphaParam(const G4Fragment & ) const {return 1.0;}
    94     virtual G4double CalcBetaParam(const G4Fragment & ) const {return 1.0;}
     96  virtual G4double CalcAlphaParam(const G4Fragment & ) const;
     97  virtual G4double CalcBetaParam(const G4Fragment & ) const;
    9598   
    9699protected:
    97100 
    98     void SetExcitationEnergiesPtr(std::vector<G4double> * anExcitationEnergiesPtr)
    99         {
    100             ExcitationEnergies = anExcitationEnergiesPtr;
    101         }
     101  inline void SetExcitationEnergiesPtr(std::vector<G4double> * anExcitationEnergiesPtr)
     102  {
     103    ExcitationEnergies = anExcitationEnergiesPtr;
     104  }
    102105 
    103     void SetExcitationSpinsPtr(std::vector<G4double> * anExcitationSpinsPtr)
    104         {
    105             ExcitationSpins = anExcitationSpinsPtr;
    106         }
     106  inline void SetExcitationSpinsPtr(std::vector<G4double> * anExcitationSpinsPtr)
     107  {
     108    ExcitationSpins = anExcitationSpinsPtr;
     109  }
    107110
    108     void SetExcitationLifetimesPtr(std::vector<G4double> * anExcitationLifetimesPtr)
    109         {
    110             ExcitationLifetimes = anExcitationLifetimesPtr;
    111         }
     111  inline void SetExcitationLifetimesPtr(std::vector<G4double> * anExcitationLifetimesPtr)
     112  {
     113    ExcitationLifetimes = anExcitationLifetimesPtr;
     114  }
    112115
    113     void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier)
    114         {
    115             theCoulombBarrierPtr = aCoulombBarrier;
    116         }
     116  inline void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier)
     117  {
     118    theCoulombBarrierPtr = aCoulombBarrier;
     119  }
    117120
    118  
    119     // Default constructor
    120     G4GEMProbability() {}
    121121private:
    122     // Copy constructor
    123     G4GEMProbability(const G4GEMProbability &right);
     122
     123  // Copy constructor
     124  G4GEMProbability(const G4GEMProbability &right);
    124125   
    125     const G4GEMProbability & operator=(const G4GEMProbability &right);
    126     G4bool operator==(const G4GEMProbability &right) const;
    127     G4bool operator!=(const G4GEMProbability &right) const;
     126  const G4GEMProbability & operator=(const G4GEMProbability &right);
     127  G4bool operator==(const G4GEMProbability &right) const;
     128  G4bool operator!=(const G4GEMProbability &right) const;
    128129   
    129130public:
    130     G4double EmissionProbability(const G4Fragment & fragment, const G4double anEnergy);
     131
     132  G4double EmissionProbability(const G4Fragment & fragment, const G4double anEnergy);
    131133 
    132134private:
    133135
    134     G4double CalcProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy,
    135                              const G4double V);
    136     virtual G4double CCoeficient(const G4double ) const {return 0.0;};
     136  G4double CalcProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy,
     137                           const G4double V);
    137138
     139  virtual G4double CCoeficient(const G4double ) const;
     140
     141  inline G4double I0(const G4double t);
     142  inline G4double I1(const G4double t, const G4double tx);
     143  inline G4double I2(const G4double s, const G4double sx);
     144  G4double I3(const G4double s, const G4double sx);
    138145   
    139     G4double I0(const G4double t);
    140     G4double I1(const G4double t, const G4double tx);
    141     G4double I2(const G4double s, const G4double sx);
    142     G4double I3(const G4double s, const G4double sx);
     146  // Data Members
     147
     148  G4Pow*   fG4pow;
     149  G4PairingCorrection* fPairCorr;
    143150   
    144     // Data Members
     151  G4VLevelDensityParameter * theEvapLDPptr;
     152       
     153  G4int theA;
     154  G4int theZ;
    145155   
    146     G4VLevelDensityParameter * theEvapLDPptr;
    147        
    148     G4int theA;
    149     G4int theZ;
     156  // Spin is fragment spin
     157  G4double Spin;
     158
     159  // Coulomb Barrier
     160  const G4VCoulombBarrier * theCoulombBarrierPtr;
     161 
     162  // Resonances Energy
     163  std::vector<G4double> * ExcitationEnergies;
    150164   
    151     // Spin is fragment spin
    152     G4double Spin;
     165  // Resonances Spin
     166  std::vector<G4double> * ExcitationSpins;
    153167
    154     // Coulomb Barrier
    155     const G4VCoulombBarrier * theCoulombBarrierPtr;
     168  // Resonances half lifetime
     169  std::vector<G4double> * ExcitationLifetimes;
    156170
    157    
    158     // Resonances Energy
    159     std::vector<G4double> * ExcitationEnergies;
    160    
    161     // Resonances Spin
    162     std::vector<G4double> * ExcitationSpins;
    163 
    164     // Resonances half lifetime
    165     std::vector<G4double> * ExcitationLifetimes;
    166 
    167 
    168     // Normalization
    169     G4double Normalization;
     171  // Normalization
     172  G4double Normalization;
    170173   
    171174};
    172175
     176inline G4double G4GEMProbability::I0(const G4double t)
     177{
     178  return std::exp(t) - 1.0;
     179}
     180
     181inline G4double G4GEMProbability::I1(const G4double t, const G4double tx)
     182{
     183  return (t - tx + 1.0)*std::exp(tx) - t - 1.0;
     184}
     185
     186
     187inline G4double G4GEMProbability::I2(const G4double s, const G4double sx)
     188{
     189  G4double S = 1.0/std::sqrt(s);
     190  G4double Sx = 1.0/std::sqrt(sx);
     191 
     192  G4double p1 = S*S*S*( 1.0 + S*S*( 1.5 + 3.75*S*S) );
     193  G4double p2 = Sx*Sx*Sx*( 1.0 + Sx*Sx*( 1.5 + 3.75*Sx*Sx) )*std::exp(sx-s);
     194 
     195  return p1-p2;
     196}
     197
    173198
    174199#endif

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4EvaporationGEMFactory.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4EvaporationGEMFactory.cc,v 1.9 2009/09/15 12:54:16 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4EvaporationGEMFactory.cc,v 1.10 2010/04/27 11:43:16 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    103103#include "G4PhotonEvaporation.hh"
    104104
    105 const G4EvaporationGEMFactory &
    106 G4EvaporationGEMFactory::operator=(const G4EvaporationGEMFactory & )
    107 {
    108   throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationGEMFactory::operator= meant to not be accessable.");
    109   return *this;
    110 }
    111 
    112 G4bool
    113 G4EvaporationGEMFactory::operator==(const G4EvaporationGEMFactory & ) const
    114 {
    115   throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationGEMFactory::operator== meant to not be accessable.");
    116   return false;
    117 }
    118 
    119 G4bool
    120 G4EvaporationGEMFactory::operator!=(const G4EvaporationGEMFactory & ) const
    121 {
    122   throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationGEMFactory::operator!= meant to not be accessable.");
    123   return true;
    124 }
     105G4EvaporationGEMFactory::G4EvaporationGEMFactory()
     106{}
     107 
     108G4EvaporationGEMFactory::~G4EvaporationGEMFactory()
     109{}
    125110                 
    126111std::vector<G4VEvaporationChannel*> * G4EvaporationGEMFactory::CreateChannel()
    … …  
    129114    new std::vector<G4VEvaporationChannel*>;
    130115  theChannel->reserve(68);
     116
     117  theChannel->push_back( new G4PhotonEvaporation() );  // Photon Channel
     118  theChannel->push_back( new G4CompetitiveFission() ); // Fission Channel
    131119
    132120  theChannel->push_back( new G4NeutronGEMChannel() );  // n
    … …  
    197185  theChannel->push_back( new G4Mg28GEMChannel() );     // Mg28
    198186
    199   theChannel->push_back( new G4CompetitiveFission() ); // Fission Channel
    200   theChannel->push_back( new G4PhotonEvaporation() );  // Photon Channel
    201 
    202187  return theChannel;
    203188}

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMProbability.cc

    r1196 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4GEMProbability.cc,v 1.13 2010/05/19 10:21:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
     28//
     29//---------------------------------------------------------------------
     30//
     31// Geant4 class G4GEMProbability
     32//
     33//
     34// Hadronic Process: Nuclear De-excitations
     35// by V. Lara (Sept 2001)
    2636//
    2737//
    … …  
    3545//       -InitialLeveldensity calculated according to the right conditions of the
    3646//        initial excited nucleus.
     47// J. M. Quesada 19/04/2010 fix in  emission probability calculation.
     48// V.Ivanchenko  20/04/2010 added usage of G4Pow and use more safe computation
     49// V.Ivanchenko 18/05/2010 trying to speedup the most slow method
     50//            by usage of G4Pow, integer Z and A; moved constructor,
     51//            destructor and virtual functions to source
    3752
    3853#include "G4GEMProbability.hh"
    3954#include "G4PairingCorrection.hh"
    40 
    41 
    42 
    43 G4GEMProbability::G4GEMProbability(const G4GEMProbability &) : G4VEmissionProbability()
    44 {
    45     throw G4HadronicException(__FILE__, __LINE__, "G4GEMProbability::copy_constructor meant to not be accessable");
    46 }
    47 
    48 
    49 
    50 
    51 const G4GEMProbability & G4GEMProbability::
    52 operator=(const G4GEMProbability &)
    53 {
    54   throw G4HadronicException(__FILE__, __LINE__, "G4GEMProbability::operator= meant to not be accessable");
    55   return *this;
    56 }
    57 
    58 
    59 G4bool G4GEMProbability::operator==(const G4GEMProbability &) const
    60 {
    61   return false;
    62 }
    63 
    64 G4bool G4GEMProbability::operator!=(const G4GEMProbability &) const
    65 {
    66   return true;
     55#include "G4Pow.hh"
     56#include "G4IonTable.hh"
     57
     58G4GEMProbability:: G4GEMProbability() :
     59  theA(0), theZ(0), Spin(0.0), theCoulombBarrierPtr(0),
     60  ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
     61{
     62  theEvapLDPptr = new G4EvaporationLevelDensityParameter;
     63  fG4pow = G4Pow::GetInstance();
     64  fPairCorr = G4PairingCorrection::GetInstance();
     65}
     66
     67G4GEMProbability:: G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin) :
     68  theA(anA), theZ(aZ), Spin(aSpin), theCoulombBarrierPtr(0),
     69  ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
     70{
     71  theEvapLDPptr = new G4EvaporationLevelDensityParameter;
     72  fG4pow = G4Pow::GetInstance();
     73  fPairCorr = G4PairingCorrection::GetInstance();
     74}
     75   
     76G4GEMProbability::~G4GEMProbability()
     77{
     78  delete theEvapLDPptr;
     79}
     80
     81G4double G4GEMProbability::CalcAlphaParam(const G4Fragment & ) const
     82{
     83  return 1.0;
     84}
     85
     86G4double G4GEMProbability::CalcBetaParam(const G4Fragment & ) const
     87{
     88  return 1.0;
     89}
     90
     91G4double G4GEMProbability::CCoeficient(const G4double ) const
     92{
     93  return 0.0;
    6794}
    6895
    … …  
    80107    if (ExcitationEnergies  &&  ExcitationSpins && ExcitationLifetimes) {
    81108      G4double SavedSpin = Spin;
    82       for (unsigned int i = 0; i < ExcitationEnergies->size(); i++) {
     109      for (size_t i = 0; i < ExcitationEnergies->size(); ++i) {
    83110        Spin = ExcitationSpins->operator[](i);
    84111        // substract excitation energies
    … …  
    86113        if (Tmax > 0.0) {
    87114          G4double width = CalcProbability(fragment,Tmax,CoulombBarrier);
     115          //JMQ April 2010 added condition to prevent reported crash
    88116          // update probability
    89           if (hbar_Planck*std::log(2.0)/width < ExcitationLifetimes->operator[](i)) {
     117          if (width > 0. &&
     118              hbar_Planck*fG4pow->logZ(2) < width*ExcitationLifetimes->operator[](i)) {
    90119            EmissionProbability += width;
    91120          }
    … …  
    98127  return EmissionProbability;
    99128}
     129
     130
     131G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
     132                                           const G4double MaximalKineticEnergy,
     133                                           const G4double V)
     134
     135// Calculate integrated probability (width) for evaporation channel
     136{
     137  G4int A = fragment.GetA_asInt();
     138  G4int Z = fragment.GetZ_asInt();
     139
     140  G4int ResidualA = A - theA;
     141  G4int ResidualZ = Z - theZ;
     142  G4double U = fragment.GetExcitationEnergy();
     143 
     144  G4double NuclearMass = fragment.ComputeGroundStateMass(theZ, theA);
     145 
     146  G4double Alpha = CalcAlphaParam(fragment);
     147  G4double Beta = CalcBetaParam(fragment);
     148 
     149 
     150  //                       ***RESIDUAL***
     151  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
     152 
     153  G4double delta0 = fPairCorr->GetPairingCorrection(ResidualA, ResidualZ); 
     154 
     155  G4double a = theEvapLDPptr->LevelDensityParameter(ResidualA,
     156                                                    ResidualZ,MaximalKineticEnergy+V-delta0);
     157  G4double Ux = (2.5 + 150.0/G4double(ResidualA))*MeV;
     158  G4double Ex = Ux + delta0;
     159  G4double T  = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
     160  //JMQ fixed bug in units
     161  G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
     162  //                      ***end RESIDUAL ***
     163 
     164  //                       ***PARENT***
     165  //JMQ (September 2009) the following quantities refer to the PARENT:
     166     
     167  G4double deltaCN=fPairCorr->GetPairingCorrection(A, Z);                                     
     168  G4double aCN = theEvapLDPptr->LevelDensityParameter(A, Z, U-deltaCN);
     169  G4double UxCN = (2.5 + 150.0/G4double(A))*MeV;
     170  G4double ExCN = UxCN + deltaCN;
     171  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
     172  //JMQ fixed bug in units
     173  G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
     174//                       ***end PARENT***
     175
     176  G4double Width;
     177  G4double InitialLevelDensity;
     178  G4double t = MaximalKineticEnergy/T;
     179  if ( MaximalKineticEnergy < Ex ) {
     180    //JMQ 190709 bug in I1 fixed (T was  missing)
     181    Width = (I1(t,t)*T + (Beta+V)*I0(t))/std::exp(E0/T);
     182    //JMQ 160909 fix:  InitialLevelDensity has been taken away
     183    //(different conditions for initial CN..)
     184  } else {
     185
     186    //VI minor speedup
     187    G4double expE0T = std::exp(E0/T);
     188    const G4double sqrt2 = std::sqrt(2.0);
     189
     190    G4double tx = Ex/T;
     191    G4double s = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
     192    G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
     193    Width = I1(t,tx)*T/expE0T + I3(s,sx)*std::exp(s)/(sqrt2*a);
     194    // For charged particles (Beta+V) = 0 beacuse Beta = -V
     195    if (theZ == 0) {
     196      Width += (Beta+V)*(I0(tx)/expE0T + 2.0*sqrt2*I2(s,sx)*std::exp(s));
     197    }
     198  }
     199 
     200  //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of hbar_planck must be used
     201  //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
     202  G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
     203 
     204  //JMQ 190709 fix on Rb and  geometrical cross sections according to Furihata's paper
     205  //                      (JAERI-Data/Code 2001-105, p6)
     206  //    G4double RN = 0.0;
     207  G4double Rb = 0.0;
     208  if (theA > 4)
     209    {
     210      //        G4double R1 = std::pow(ResidualA,1.0/3.0);
     211      //        G4double R2 = std::pow(G4double(theA),1.0/3.0);
     212      G4double Ad = fG4pow->Z13(ResidualA);
     213      G4double Aj = fG4pow->Z13(theA);
     214      //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
     215      Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
     216      Rb *= fermi;
     217    }
     218  else if (theA>1)
     219    {
     220      G4double Ad = fG4pow->Z13(ResidualA);
     221      G4double Aj = fG4pow->Z13(theA);
     222      Rb=1.5*(Aj+Ad)*fermi;
     223    }
     224  else
     225    {
     226      G4double Ad = fG4pow->Z13(ResidualA);
     227      Rb = 1.5*Ad*fermi;
     228    }
     229  //   G4double GeometricalXS = pi*RN*RN*std::pow(ResidualA,2./3.);
     230  G4double GeometricalXS = pi*Rb*Rb;
     231  //end of JMQ fix on Rb by 190709
     232 
     233
     234
     235//JMQ 160909 fix: initial level density must be calculated according to the
     236// conditions at the initial compound nucleus
     237// (it has been removed from previous "if" for the residual)
     238
     239   if ( U < ExCN )
     240      {
     241        InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
     242      }
     243    else
     244      {
     245        //        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
     246        //VI speedup
     247        G4double x  = U-deltaCN;
     248        G4double x2 = x*x;
     249        G4double x3 = x2*x;
     250        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*x))/std::sqrt(std::sqrt(aCN*x2*x3));
     251      }
     252 //
     253
     254
     255  //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
     256  //    Width *= std::sqrt(pi)*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
     257  Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
     258   
     259
     260  return Width;
     261}
     262
     263/*
    100264
    101265G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
    … …  
    219383  return Width;
    220384}
    221 
    222 
    223 
    224 
    225 G4double G4GEMProbability::I0(const G4double t)
    226 {
    227   G4double result = (std::exp(t) - 1.0);
    228   return result;
    229 }
    230 
    231 G4double G4GEMProbability::I1(const G4double t, const G4double tx)
    232 {
    233   G4double result = t - tx + 1.0;
    234   result *= std::exp(tx);
    235   result -= (t + 1.0);
    236   return result;
    237 }
    238 
    239 
    240 G4double G4GEMProbability::I2(const G4double s, const G4double sx)
    241 {
    242   G4double S = 1.0/std::sqrt(s);
    243   G4double Sx = 1.0/std::sqrt(sx);
    244  
    245   G4double p1 = S*S*S*( 1.0 + S*S*( 1.5 + 3.75*S*S) );
    246   G4double p2 = Sx*Sx*Sx*( 1.0 + Sx*Sx*( 1.5 + 3.75*Sx*Sx) )*std::exp(sx-s);
    247  
    248   return p1-p2;
    249 }
     385*/
    250386
    251387G4double G4GEMProbability::I3(const G4double s, const G4double sx)
Note: See TracChangeset for help on using the changeset viewer.