Ignore:
Timestamp:
Dec 22, 2010, 3:52:27 PM (13 years ago)
Author:
garnier
Message:

geant4 tag 9.4

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1340 r1347  
    2424// ********************************************************************
    2525//
    26 // $Id: G4GEMProbability.cc,v 1.13 2010/05/19 10:21:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     26// $Id: G4GEMProbability.cc,v 1.15 2010/11/05 14:43:27 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2828//
    2929//---------------------------------------------------------------------
    … …  
    6565}
    6666
    67 G4GEMProbability:: G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin) :
     67G4GEMProbability:: G4GEMProbability(G4int anA, G4int aZ, G4double aSpin) :
    6868  theA(anA), theZ(aZ), Spin(aSpin), theCoulombBarrierPtr(0),
    6969  ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
    … …  
    9595
    9696G4double G4GEMProbability::EmissionProbability(const G4Fragment & fragment,
    97                                                const G4double MaximalKineticEnergy)
     97                                               G4double MaximalKineticEnergy)
    9898{
    9999  G4double EmissionProbability = 0.0;
    … …  
    130130
    131131G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
    132                                            const G4double MaximalKineticEnergy,
    133                                            const G4double V)
     132                                           G4double MaximalKineticEnergy,
     133                                           G4double V)
    134134
    135135// Calculate integrated probability (width) for evaporation channel
    … …  
    146146  G4double Alpha = CalcAlphaParam(fragment);
    147147  G4double Beta = CalcBetaParam(fragment);
    148  
    149  
     148   
    150149  //                       ***RESIDUAL***
    151150  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
    … …  
    159158  G4double T  = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
    160159  //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));
     160  G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0
     161        - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
    162162  //                      ***end RESIDUAL ***
    163163 
    … …  
    165165  //JMQ (September 2009) the following quantities refer to the PARENT:
    166166     
    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***
     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  //                     ***end PARENT***
    175173
    176174  G4double Width;
    … …  
    231229  //end of JMQ fix on Rb by 190709
    232230 
    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 
     231  //JMQ 160909 fix: initial level density must be calculated according to the
     232  // conditions at the initial compound nucleus
     233  // (it has been removed from previous "if" for the residual)
     234
     235  if ( U < ExCN )
     236    {
     237      //JMQ fixed bug in units
     238      //VI moved the computation here
     239      G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0
     240                                  - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
     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 x1 = std::sqrt(aCN*x);
     249      InitialLevelDensity = (pi/12.0)*std::exp(2*x1)/(x*std::sqrt(x1));
     250    }
    254251
    255252  //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
    … …  
    257254  Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
    258255   
    259 
    260256  return Width;
    261257}
    262258
    263 /*
    264 
    265 G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
    266                                            const G4double MaximalKineticEnergy,
    267                                            const G4double V)
    268 // Calculate integrated probability (width) for evaporation channel
    269 {
    270   G4double ResidualA = static_cast<G4double>(fragment.GetA() - theA);
    271   G4double ResidualZ = static_cast<G4double>(fragment.GetZ() - theZ);
    272   G4double U = fragment.GetExcitationEnergy();
    273  
    274   G4double NuclearMass = G4ParticleTable::GetParticleTable()->
    275     GetIonTable()->GetNucleusMass(theZ,theA);
    276  
    277  
    278   G4double Alpha = CalcAlphaParam(fragment);
    279   G4double Beta = CalcBetaParam(fragment);
    280  
    281  
    282   //                       ***RESIDUAL***
    283   //JMQ (September 2009) the following quantities refer to the RESIDUAL:
    284  
    285   G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection( static_cast<G4int>( ResidualA ) , static_cast<G4int>( ResidualZ ) ); 
    286  
    287   G4double a = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),
    288                                                     static_cast<G4int>(ResidualZ),MaximalKineticEnergy+V-delta0);
    289   G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
    290   G4double Ex = Ux + delta0;
    291   G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
    292   //JMQ fixed bug in units
    293   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));
    294   //                      ***end RESIDUAL ***
    295  
    296   //                       ***PARENT***
    297   //JMQ (September 2009) the following quantities refer to the PARENT:
    298      
    299   G4double deltaCN=G4PairingCorrection::GetInstance()->
    300     GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));                                     
    301   G4double aCN = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
    302                                                       static_cast<G4int>(fragment.GetZ()),U-deltaCN);
    303   G4double UxCN = (2.5 + 150.0/fragment.GetA())*MeV;
    304   G4double ExCN = UxCN + deltaCN;
    305   G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
    306   //JMQ fixed bug in units
    307   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));
    308 //                       ***end PARENT***
    309 
    310   G4double Width;
    311   G4double InitialLevelDensity;
    312   G4double t = MaximalKineticEnergy/T;
    313   if ( MaximalKineticEnergy < Ex ) {
    314     //JMQ 190709 bug in I1 fixed (T was  missing)
    315     Width = (I1(t,t)*T + (Beta+V)*I0(t))/std::exp(E0/T);
    316 //JMQ 160909 fix:  InitialLevelDensity has been taken away (different conditions for initial CN..)
    317   } else {
    318     G4double tx = Ex/T;
    319     G4double s = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
    320     G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
    321     Width = I1(t,tx)*T/std::exp(E0/T) + I3(s,sx)*std::exp(s)/(std::sqrt(2.0)*a);
    322     // For charged particles (Beta+V) = 0 beacuse Beta = -V
    323     if (theZ == 0) {
    324       Width += (Beta+V)*(I0(tx)/std::exp(E0/T) + 2.0*std::sqrt(2.0)*I2(s,sx)*std::exp(s));
    325     }
    326   }
    327  
    328   //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of hbar_planck must be used
    329   //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
    330   G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
    331  
    332   //JMQ 190709 fix on Rb and  geometrical cross sections according to Furihata's paper
    333   //                      (JAERI-Data/Code 2001-105, p6)
    334   //    G4double RN = 0.0;
    335   G4double Rb = 0.0;
    336   if (theA > 4)
    337     {
    338       //        G4double R1 = std::pow(ResidualA,1.0/3.0);
    339       //        G4double R2 = std::pow(G4double(theA),1.0/3.0);
    340       G4double Ad = std::pow(ResidualA,1.0/3.0);
    341       G4double Aj = std::pow(G4double(theA),1.0/3.0);
    342       //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
    343       Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
    344       Rb *= fermi;
    345     }
    346   else if (theA>1)
    347     {
    348       G4double Ad = std::pow(ResidualA,1.0/3.0);
    349       G4double Aj = std::pow(G4double(theA),1.0/3.0);
    350       Rb=1.5*(Aj+Ad)*fermi;
    351     }
    352   else
    353     {
    354       G4double Ad = std::pow(ResidualA,1.0/3.0);
    355       Rb = 1.5*Ad*fermi;
    356     }
    357   //   G4double GeometricalXS = pi*RN*RN*std::pow(ResidualA,2./3.);
    358   G4double GeometricalXS = pi*Rb*Rb;
    359   //end of JMQ fix on Rb by 190709
    360  
    361 
    362 
    363 //JMQ 160909 fix: initial level density must be calculated according to the
    364 // conditions at the initial compound nucleus
    365 // (it has been removed from previous "if" for the residual)
    366 
    367    if ( U < ExCN )
    368       {
    369         InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
    370       }
    371     else
    372       {
    373         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);
    374       }
    375  //
    376 
    377 
    378   //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
    379   //    Width *= std::sqrt(pi)*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
    380   Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
    381    
    382 
    383   return Width;
    384 }
    385 */
    386 
    387 G4double G4GEMProbability::I3(const G4double s, const G4double sx)
     259G4double G4GEMProbability::I3(G4double s, G4double sx)
    388260{
    389261  G4double s2 = s*s;
Note: See TracChangeset for help on using the changeset viewer.