Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (15 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/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

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