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/pre_equilibrium/exciton_model/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PreCompoundAlpha.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PreCompoundAlpha.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3939// Modified: 
    4040// 21.08.2008 J. M. Quesada add choice of options 
    41 // 10.02.2009 J. M. Quesada set default opt1 
    4241//
    4342
    … …  
    243242  ecut = (ecut-a) / (p+p);
    244243  ecut2 = ecut;
    245   if (cut < 0.) ecut2 = ecut - 2.;
     244//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     245//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     246//  if (cut < 0.) ecut2 = ecut - 2.;
     247  if (cut < 0.) ecut2 = ecut;
    246248  elab = K * FragmentA / ResidualA;
    247249  sig = 0.;

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

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PreCompoundDeuteron.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4PreCompoundDeuteron.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    3838// Modified: 
    3939// 21.08.2008 J. M. Quesada add choice of options 
    40 // 10.02.2009 J. M. Quesada set default opt1 
    4140//
    4241
    … …  
    234233  ecut = (ecut-a) / (p+p);
    235234  ecut2 = ecut;
    236   if (cut < 0.) ecut2 = ecut - 2.;
     235//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     236//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     237//  if (cut < 0.) ecut2 = ecut - 2.;
     238  if (cut < 0.) ecut2 = ecut;
    237239  elab = K * FragmentA / ResidualA;
    238240  sig = 0.;

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PreCompoundEmission.cc,v 1.22 2009/11/13 17:40:14 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PreCompoundEmission.cc,v 1.28 2010/02/25 10:27:36 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3838//
    3939// Modified: 
     40// 15.01.2010 J.M.Quesada  added protection against unphysical values of parameter an
     41// 19.01.2010 V.Ivanchenko simplified computation of parameter an, sample cosTheta
     42//                         instead of theta; protect all calls to sqrt
    4043//
    4144
    … …  
    4649#include "G4HETCEmissionFactory.hh"
    4750
    48 const G4PreCompoundEmission & G4PreCompoundEmission::operator=(const G4PreCompoundEmission &)
    49 {
    50   throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundEmission::operator= meant to not be accessable");
     51const G4PreCompoundEmission &
     52G4PreCompoundEmission::operator=(const G4PreCompoundEmission &)
     53{
     54  throw G4HadronicException(__FILE__, __LINE__,
     55                            "G4PreCompoundEmission::operator= meant to not be accessable");
    5156  return *this;
    5257}
    … …  
    166171  // check that Excitation energy is >= 0
    167172  G4double anU = RestMomentum.m()-theFragment->GetRestNuclearMass();
    168   if (anU < 0.0) throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundModel::DeExcite: Excitation energy less than 0!");
    169    
    170    
    171    
     173  if (anU < 0.0) {
     174    throw G4HadronicException(__FILE__, __LINE__,
     175                              "G4PreCompoundModel::DeExcite: Excitation energy less than 0!");
     176  }
     177     
    172178  // Update nucleus parameters:
    173179  // --------------------------
    … …  
    204210}
    205211
    206 
    207 G4ThreeVector G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment * theFragment,
    208                                                         const G4Fragment& aFragment,
    209                                                          const G4double KineticEnergyOfEmittedFragment) const
     212G4ThreeVector
     213G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment * theFragment,
     214                                          const G4Fragment& aFragment,
     215                                           const G4double kinEnergyOfEmittedFrag) const
    210216{
    211217  G4double p = aFragment.GetNumberOfParticles();
    212218  G4double h = aFragment.GetNumberOfHoles();
    213219  G4double U = aFragment.GetExcitationEnergy();
    214        
    215   // Kinetic Energy of emitted fragment
    216   // G4double KineticEnergyOfEmittedFragment = theFragment->GetKineticEnergy(aFragment);
     220
     221  G4double ekin = std::max(0.0, kinEnergyOfEmittedFrag);
    217222       
    218223  // Emission particle separation energy
    219224  G4double Bemission = theFragment->GetBindingEnergy();
    220        
     225
    221226  // Fermi energy
    222227  G4double Ef = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
    … …  
    225230  //  G4EvaporationLevelDensityParameter theLDP;
    226231  //  G4double g = (6.0/pi2)*aFragment.GetA()*
    227   //    theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);
    228   G4double g = (6.0/pi2)*aFragment.GetA()*
    229     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     232
     233  G4double g = (6.0/pi2)*aFragment.GetA()
     234    *G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    230235       
    231236  // Average exciton energy relative to bottom of nuclear well
    … …  
    236241  //  G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission;
    237242
    238 
    239 
    240243  G4double w_num = rho(p+1,h,g,Uf,Ef);
    241244  G4double w_den = rho(p,h,g,Uf,Ef);
    … …  
    251254 
    252255
    253   G4double zeta = std::max(1.0,9.3/std::sqrt(KineticEnergyOfEmittedFragment/MeV));
    254        
    255   G4double an = 3.0*std::sqrt((ProjEnergy+Ef)*(KineticEnergyOfEmittedFragment+Bemission+Ef));
    256   if (aFragment.GetNumberOfExcitons() == 1)
    257     {
    258       an /= (zeta*2.0*aFragment.GetNumberOfExcitons()*Eav);
    259     }
    260   else
    261     {
    262       an /= (zeta*(aFragment.GetNumberOfExcitons()-1.0)*Eav);
    263     }
     256  // VI + JMQ 19/01/2010 update computation of the parameter an
     257  //
     258  G4double an = 0.0;
     259  G4double Eeff = ekin + Bemission + Ef;
     260  if(ekin > DBL_MIN && Eeff > DBL_MIN) {
     261
     262    G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV));
     263 
     264    an = 3.0*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav);
     265
     266    G4int ne = aFragment.GetNumberOfExcitons() - 1;
     267    if ( ne > 1 ) { an /= (G4double)ne; }
    264268                       
    265                        
    266 //    G4double expan=std::exp(an);
    267 //    thetaold = std::acos(std::log(expan-random*(expan-1.0/expan))/an);
    268 //      exp(an) becomes large for large an
    269 //      1/exp(an) becomes large for large negative an
    270 //       take the large value out of the log depending on the sign of an to avoid FP exceptions
     269    // protection of exponent
     270    if ( an > 10. ) { an = 10.; }
     271  }
     272
     273  // sample cosine of theta and not theta as in old versions 
     274  G4double random = G4UniformRand();
     275  G4double cost;
    271276 
    272   G4double random=G4UniformRand();
    273   G4double exp2an=0;
    274   G4double theta;
    275   if (an > 0.) {
    276     if (an < 25.) { exp2an = std::exp(-2*an); } // we subtract from 1, exp(-50)~1e-21, so will not
    277                                            //  change numerical result
    278     theta = std::acos(1+ std::log(1-random*(1-exp2an))/an);
    279   } else if ( an < 0.) {
    280     if ( an > -25.) { exp2an = std::exp(2*an); } // similar to above, except we compare to rndm*1
    281     theta = std::acos(std::log(exp2an-random*(exp2an-1))/an - 1.);
    282   } else {   // an==0 now.     
    283     theta=std::acos(1.-2*random);
    284   }
    285  
    286   if (std::abs(an) < 50 )
    287   {
    288      
    289   }
     277  if(an < 0.1) { cost = 1. - 2*random; }
     278  else {
     279    G4double exp2an = std::exp(-2*an);
     280    cost = 1. + std::log(1-random*(1-exp2an))/an;
     281    if(cost > 1.) { cost = 1.; }
     282    else if(cost < -1.) {cost = -1.; }
     283  } 
    290284
    291285  G4double phi = twopi*G4UniformRand();
    292286 
    293287  // Calculate the momentum magnitude of emitted fragment       
    294   G4double EmittedMass = theFragment->GetNuclearMass();
    295   G4double pmag = std::sqrt(KineticEnergyOfEmittedFragment*(KineticEnergyOfEmittedFragment+2.0*EmittedMass));
    296  
    297  
    298   G4double sinTheta = std::sin(theta);
    299   //  G4double cosTheta = std::sqrt(1.0-sinTheta*sinTheta);
    300   G4double cosTheta = std::cos(theta);
    301 
    302   G4ThreeVector momentum(pmag*std::cos(phi)*sinTheta,pmag*std::sin(phi)*sinTheta,pmag*cosTheta);
     288  G4double pmag = std::sqrt(ekin*(ekin + 2.0*theFragment->GetNuclearMass()));
     289 
     290  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
     291
     292  G4ThreeVector momentum(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost);
    303293  // theta is the angle wrt the incident direction
    304294  momentum.rotateUz(theIncidentDirection);
    … …  
    309299G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g,
    310300                                    const G4double E, const G4double Ef) const
    311 {
    312        
    313   G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g);
    314   G4double alpha = (p*p+h*h)/(2.0*g);
    315  
    316   if ( (E-alpha) < 0 ) return 0;
     301{       
     302  // 25.02.2010 V.Ivanchenko added more protections
     303  G4double Aph   = (p*p + h*h + p - 3.0*h)/(4.0*g);
     304  //  G4double alpha = (p*p + h*h)/(2.0*g);
     305 
     306  if ( E - Aph < 0.0) { return 0.0; }
    317307 
    318308  G4double logConst =  (p+h)*std::log(g) - logfactorial(p+h-1) - logfactorial(p) - logfactorial(h);
    319309
    320 // initialise values using j=0
     310  // initialise values using j=0
    321311
    322312  G4double t1=1;
    323313  G4double t2=1;
    324   G4double logt3=(p+h-1) * std::log(E-Aph);
    325   G4double tot = std::exp( logt3 + logConst );
    326 
    327 // and now sum rest of terms
    328   G4int j(1); 
    329   while ( (j <= h) && ((E - alpha - j*Ef) > 0.0) )
    330     {
    331           t1 *= -1.;
    332           t2 *= (h+1-j)/j;
    333           logt3 = (p+h-1) * std::log( E - j*Ef - Aph) + logConst;
    334           G4double t3 = std::exp(logt3);
    335           tot += t1*t2*t3;
    336           j++;
     314  G4double logt3=(p+h-1) * std::log(E-Aph) + logConst;
     315  const G4double logmax = 200.;
     316  if(logt3 > logmax) { logt3 = logmax; }
     317  G4double tot = std::exp( logt3 );
     318
     319  // and now sum rest of terms
     320  // 25.02.2010 V.Ivanchenko change while to for loop and cleanup
     321  G4double Eeff = E - Aph;
     322  for(G4int j=1; j<=h; ++j)
     323    {
     324      Eeff -= Ef;
     325      if(Eeff < 0.0) { break; }
     326      t1 *= -1.;
     327      t2 *= (G4double)(h+1-j)/(G4double)j;
     328      logt3 = (p+h-1) * std::log( Eeff) + logConst;
     329      if(logt3 > logmax) { logt3 = logmax; }
     330      tot += t1*t2*std::exp(logt3);
    337331    }
    338332       
    339333  return tot;
    340334}
    341 
    342 
    343335
    344336G4double G4PreCompoundEmission::factorial(G4double a) const

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PreCompoundHe3.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PreCompoundHe3.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3939// Modified: 
    4040// 21.08.2008 J. M. Quesada add choice of options 
    41 // 10.02.2009 J. M. Quesada set default opt1 
    4241//
    4342 
    … …  
    242241  ecut = (ecut-a) / (p+p);
    243242  ecut2 = ecut;
    244   if (cut < 0.) ecut2 = ecut - 2.;
     243//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     244//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     245//  if (cut < 0.) ecut2 = ecut - 2.;
     246  if (cut < 0.) ecut2 = ecut;
    245247  elab = K * FragmentA / ResidualA;
    246248  sig = 0.;

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PreCompoundProton.cc,v 1.4 2009/02/11 18:06:00 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PreCompoundProton.cc,v 1.5 2010/04/09 14:06:17 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    6868//OPT=0 Dostrovski's parameterization
    6969//OPT=1 Chatterjee's paramaterization
    70 //OPT=2 Wellisch's parametarization
     70//OPT=2,4 Wellisch's parametarization
    7171//OPT=3 Kalbach's parameterization
    7272//
    … …  
    224224  //     ** p from  becchetti and greenlees (but modified with sub-barrier
    225225  //     ** correction function and xp2 changed from -449)
    226   // JMQ (june 2008) : refinement of proton cross section for light systems
    227   //
     226
    228227  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
    229228  G4double p, p0, p1, p2;
    … …  
    281280  ecut = (ecut-a) / (p+p);
    282281  ecut2 = ecut;
    283   if (cut < 0.) ecut2 = ecut - 2.;
     282//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     283//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     284//  if (cut < 0.) ecut2 = ecut - 2.;
     285  if (cut < 0.) ecut2 = ecut;
    284286  elab = K * FragmentA / ResidualA;
    285287  sig = 0.;
    … …  
    290292    signor2 = 1. + std::exp(signor2);
    291293    sig = sig / signor2;
    292     // signor2 is empirical global corr'n at low elab for protons in PRECO, not enough for light targets
    293     //  refinement for proton cross section
    294     if (ResidualZ<=26)
    295       sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));
    296    
    297   }              //end for E<Ec
     294  }              //end for E<=Ec
    298295  else{           //start for  E>Ec
    299296    sig = (landa*elab+mu+nu/elab) * signor;
    300    
    301     //  refinement for proton cross section
    302     if ( ResidualZ<=26 && elab <=(afit*ResidualZ+bfit)*ec)
    303       sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));   
    304     //
    305    
    306297    geom = 0.;
    307298   

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PreCompoundTriton.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PreCompoundTriton.cc,v 1.6 2010/04/09 14:06:17 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3939// Modified: 
    4040// 21.08.2008 J. M. Quesada add choice of options 
    41 // 10.02.2009 J. M. Quesada set default opt1 
    4241//
    4342 
    … …  
    235234  ecut = (ecut-a) / (p+p);
    236235  ecut2 = ecut;
    237   if (cut < 0.) ecut2 = ecut - 2.;
     236//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     237//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     238//  if (cut < 0.) ecut2 = ecut - 2.;
     239  if (cut < 0.) ecut2 = ecut;
    238240  elab = K * FragmentA / ResidualA;
    239241  sig = 0.;
Note: See TracChangeset for help on using the changeset viewer.