Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

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

Legend:

Unmodified
Added
Removed

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

    r819 r962  
    2525//
    2626//
    27 // Hadronic Process: Nuclear Preequilibrium
    28 // by V. Lara
    29 
     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//
    3041
    3142#include "G4PreCompoundEmission.hh"
    … …  
    4253
    4354
    44   G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
     55G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
    4556{
    4657  return false;
    4758}
    4859
    49     G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
     60G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
    5061{
    5162  return true;
    … …  
    95106  return;
    96107}
    97 
    98 
    99108
    100109G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)
    … …  
    163172  // Update nucleus parameters:
    164173  // --------------------------
     174
    165175  // Number of excitons
    166176  aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
    … …  
    175185  // Charge
    176186  aFragment.SetZ(theFragment->GetRestZ());
    177    
     187
    178188   
    179189  // Perform Lorentz boosts
    … …  
    275285}
    276286
    277 
    278287G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g,
    279288                                    const G4double E, const G4double Ef) const
     289{
     290       
     291  G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g);
     292  G4double alpha = (p*p+h*h)/(2.0*g);
     293 
     294  if ( (E-alpha) < 0 ) return 0;
     295
     296  G4double factp=factorial(p);
     297
     298  G4double facth=factorial(h);
     299
     300  G4double factph=factorial(p+h-1);
     301 
     302  G4double logConst =  (p+h)*std::log(g) - std::log (factph) - std::log(factp) - std::log(facth);
     303
     304// initialise values using j=0
     305
     306  G4double t1=1;
     307  G4double t2=1;
     308  G4double logt3=(p+h-1) * std::log(E-Aph);
     309  G4double tot = std::exp( logt3 + logConst );
     310
     311// and now sum rest of terms
     312  G4int j(1); 
     313  while ( (j <= h) && ((E - alpha - j*Ef) > 0.0) )
     314    {
     315          t1 *= -1.;
     316          t2 *= (h+1-j)/j;
     317          logt3 = (p+h-1) * std::log( E - j*Ef - Aph) + logConst;
     318          G4double t3 = std::exp(logt3);
     319          tot += t1*t2*t3;
     320          j++;
     321    }
     322       
     323  return tot;
     324}
     325
     326
     327
     328G4double G4PreCompoundEmission::factorial(G4double a) const
    280329{
    281330  // Values of factorial function from 0 to 60
    … …  
    349398  //      fact[n] = fact[n-1]*static_cast<G4double>(n);
    350399  //    }
    351        
    352   G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g);
    353   G4double alpha = (p*p+h*h)/(2.0*g);
    354 
    355   G4double tot = 0.0;
    356   for (G4int j = 0; j <= h; j++)
    357     {
    358       if (E-alpha-static_cast<G4double>(j)*Ef > 0.0)
    359         {
    360           G4double t1 = std::pow(-1.0, static_cast<G4double>(j));
    361           G4double t2 = fact[static_cast<G4int>(h)]/ (fact[static_cast<G4int>(h)-j]*fact[j]);
    362           G4double t3 = E - static_cast<G4double>(j)*Ef - Aph;
    363           t3 = std::pow(t3,p+h-1);
    364           tot += t1*t2*t3;
    365         }
    366       else
     400  G4double result(0.0);
     401  G4int ia = static_cast<G4int>(a);
     402  if (ia < factablesize)
     403    {
     404      result = fact[ia];
     405    }
     406  else
     407    {
     408      result = fact[factablesize-1];
     409      for (G4int n = factablesize; n < ia+1; ++n)
    367410        {
    368           break;
     411          result *= static_cast<G4double>(n);
    369412        }
    370413    }
    371414   
    372 
    373   G4double factp(0.0);
    374   G4int ph = static_cast<G4int>(p);
    375   if (ph < factablesize)
    376     {
    377       factp = fact[ph];
    378     }
    379   else
    380     {
    381       factp = fact[factablesize-1];
    382       for (G4int n = factablesize; n < ph+1; ++n)
    383         {
    384           factp *= static_cast<G4double>(n);
    385         }
    386     }
    387 
    388   G4double facth(0.0);
    389   ph = static_cast<G4int>(h);
    390   if (ph < factablesize)
    391     {
    392       facth = fact[ph];
    393     }
    394   else
    395     {
    396       facth = fact[factablesize-1];
    397       for (G4int n = factablesize; n < ph+1; ++n)
    398         {
    399           facth *= static_cast<G4double>(n);
    400         }
    401     }
    402   G4double factph(0.0);
    403   ph = static_cast<G4int>(p+h)-1;
    404   if (ph < factablesize)
    405     {
    406       factph = fact[ph];
    407     }
    408   else
    409     {
    410       factph = fact[factablesize-1];
    411       for (G4int n = factablesize; n < ph+1; ++n)
    412         {
    413           factph *= static_cast<G4double>(n);
    414         }
    415     }
    416  
    417   tot *= std::pow(g,p+h)/factph;
    418   tot /= factp;
    419   tot /= facth;
    420    
    421   return tot;
    422 }
    423 
     415    return result;
     416}
    424417
    425418

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

    r819 r962  
    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 $
    2628//
    27 // by V. Lara
    28  
     29// J. M. Quesada (August 2008). 
     30// Based  on previous work by V. Lara
     31// JMQ (06 September 2008) Also external choice has been added for:
     32//                      - superimposed Coulomb barrier (if useSICB=true)
     33//
    2934#include "G4PreCompoundFragment.hh"
    3035
    … …  
    4146                      const G4String & aName):
    4247  G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName)
    43 {}
     48{
     49}
    4450
    4551
    … …  
    7177CalcEmissionProbability(const G4Fragment & aFragment)
    7278{
    73   if (GetEnergyThreshold() <= 0.0)
     79// If  theCoulombBarrier effect is included in the emission probabilities
     80//if (GetMaximalKineticEnergy() <= 0.0)
     81  G4double limit;
     82  if(OPTxs==0 ||  useSICB) limit= theCoulombBarrier;
     83  else limit=0.;
     84  if (GetMaximalKineticEnergy() <= limit)
    7485    {
    7586      theEmissionProbability = 0.0;
    7687      return 0.0;
    7788  }   
    78   // Coulomb barrier is the lower limit
    79   // of integration over kinetic energy
    80   G4double LowerLimit = theCoulombBarrier;
    81  
    82   // Excitation energy of nucleus after fragment emission is the upper limit
    83   // of integration over kinetic energy
    84   G4double UpperLimit = this->GetMaximalKineticEnergy();
     89// If  theCoulombBarrier effect is included in the emission probabilities
     90//  G4double LowerLimit = 0.;
     91// Coulomb barrier is the lower limit
     92// of integration over kinetic energy
     93  G4double LowerLimit = limit;
     94
     95// Excitation energy of nucleus after fragment emission is the upper limit of integration over kinetic energy
     96  G4double UpperLimit = GetMaximalKineticEnergy();
    8597 
    8698  theEmissionProbability =
    8799    IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment);
    88  
    89100  return theEmissionProbability;
    90101}
    … …  
    130141    }
    131142  Total  *= (Up-Low)/2.0;
    132 
    133143  return Total;
    134144}
    … …  
    140150GetKineticEnergy(const G4Fragment & aFragment)
    141151{
    142   G4double V = this->GetCoulombBarrier();
    143   G4double Tmax =  this->GetMaximalKineticEnergy() - V;
    144  
     152
     153//      G4double V = this->GetCoulombBarrier();// alternative way for accessing the Coulomb barrier
     154//                                             //should be equivalent (in fact it is)
     155  G4double V;
     156  if(OPTxs==0 || useSICB) V= theCoulombBarrier;//let's keep this way for consistency with CalcEmissionProbability method
     157  else V=0.;
     158
     159  G4double Tmax =  GetMaximalKineticEnergy() ; 
    145160  G4double T(0.0);
    146161  G4double NormalizedProbability(1.0);
    147162  do
    148163    {
    149       T = V + G4UniformRand()*Tmax;
    150       NormalizedProbability = this->ProbabilityDistributionFunction(T,aFragment)/
    151         this->GetEmissionProbability();
    152      
    153     }
    154   while (G4UniformRand() > NormalizedProbability);
    155  
     164      T =V+ G4UniformRand()*(Tmax-V);
     165      NormalizedProbability = ProbabilityDistributionFunction(T,aFragment)/GetEmissionProbability();     
     166    }   while (G4UniformRand() > NormalizedProbability); 
    156167  return T;
    157168}

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundFragmentVector.cc,v 1.5 2006/06/29 20:59:21 gunter Exp $
    28 // GEANT4 tag $Name:  $
     26// $Id: G4PreCompoundFragmentVector.cc,v 1.11 2009/02/10 16:01:37 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2928//
    3029// Hadronic Process: Nuclear Preequilibrium
    … …  
    8685  for (i = theChannels->begin(); i != theChannels->end(); ++i) {
    8786    accumulation += (*i)->GetEmissionProbability();
     87
    8888    running.push_back(accumulation);
    8989  }
    … …  
    126126          (*i)->IncrementStage();
    127127        }
    128     } 
     128    }
    129129
    130130  return theChannels->operator[](ChosenChannel);

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // by V. Lara
     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 
     40//
    2741
    2842#include "G4PreCompoundIon.hh"
    2943#include "G4PreCompoundParameters.hh"
    30 //#include "G4EvaporationLevelDensityParameter.hh"
    3144
     45G4bool 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}
    3251
    3352G4double G4PreCompoundIon::
    … …  
    4362  G4double H = aFragment.GetNumberOfHoles();
    4463  G4double N = P + H;
    45  
    46   // g = (6.0/pi2)*a*A
    47   //  G4EvaporationLevelDensityParameter theLDP;
     64
    4865  G4double g0 = (6.0/pi2)*aFragment.GetA() *
    4966    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    50   //    theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);
     67 
    5168  G4double g1 = (6.0/pi2)*GetRestA() *
    5269    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    53   //    theLDP.LevelDensityParameter(GetRestA(),GetRestZ(),U);
    54   //no longer used:  G4double gj = (6.0/pi2)*GetA() *
    55   //  -----"-----    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    56   //    theLDP.LevelDensityParameter(GetA(),GetZ(),U);
    5770
    58   G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
    59   G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1);
    60   G4double Aj = GetA()*(GetA()+1.0)/4.0;
     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;
     76
     77  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
     78
     79  G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1);
     80
     81  G4double Aj = GetA()*(GetA()+1.0)/4.0/gj;
     82
    6183
    6284  G4double E0 = std::max(0.0,U - A0);
    6385  if (E0 == 0.0) return 0.0;
    64   //  G4double E1 = std::max(0.0,GetMaximalKineticEnergy() + GetCoulombBarrier() - eKin - A1);
    65   G4double E1 = (std::max(0.0,GetMaximalKineticEnergy() - eKin - A1))/g1; //JMQ fix
    66   //  G4double Ej = std::max(0.0,U + GetBindingEnergy() - Aj);
    67   G4double Ej = std::max(0.0,eKin + GetBindingEnergy() - Aj); //JMQ fix
    6886
    69   G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*(eKin+GetBindingEnergy()))))*
    70       GetAlpha()*(eKin+GetBeta())/(r0*std::pow(GetRestA(),1.0/3.0)) *
    71       CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P);
     87  G4double E1 = (std::max(0.0,GetMaximalKineticEnergy() - eKin - A1));
     88
     89  G4double Ej = std::max(0.0,eKin + GetBindingEnergy() -Aj);
     90
     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());
    7297
    7398  G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0);
     99 
     100  G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0;
    74101
    75   //  G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0;
    76   G4double pC = std::pow((g1*Ej)/(g0*E0),GetA()-1.0)*(g1/g0)/E0; //JMQ fix
    77 
    78   G4double Probability = pA * pB * pC * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged());
     102  G4double Probability = pA * pB * pC;
    79103
    80104  return Probability;

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4PreCompoundModel.cc,v 1.11 2007/10/11 14:19:36 ahoward Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4PreCompoundModel.cc,v 1.17 2008/12/09 14:09:59 ahoward Exp $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// by V. Lara
     31//
     32//J. M. Quesada (Apr.08). Several changes. Soft cut-off switched off.
     33//(May. 08). Protection against non-physical preeq. transitional regime has
     34// been set
     35//
     36// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     37// cross section option
     38// JMQ (06 September 2008) Also external choices have been added for:
     39//                      - superimposed Coulomb barrier (useSICB=true)
     40//                      - "never go back"  hipothesis (useNGB=true)
     41//                      - soft cutoff from preeq. to equlibrium (useSCO=true)
     42//                      - CEM transition probabilities (useCEMtr=true) 
     43
    3144
    3245#include "G4PreCompoundModel.hh"
    … …  
    3447#include "G4PreCompoundTransitions.hh"
    3548#include "G4GNASHTransitions.hh"
     49#include "G4ParticleDefinition.hh"
    3650
    3751
    … …  
    160174  // Copy of the initial state
    161175  G4Fragment aFragment(theInitialState);
     176
     177  if (aFragment.GetExcitationEnergy() < 10*eV)
     178    {
     179      // Perform Equilibrium Emission
     180      PerformEquilibriumEmission(aFragment,Result);
     181      return Result;
     182    }
    162183 
    163184  if (aFragment.GetA() < 5) {
    … …  
    175196  if (useHETCEmission) aEmission.SetHETCModel();
    176197  aEmission.SetUp(theInitialState);
    177 
     198 
     199  //for cross section options
     200 
     201  if (OPTxs!= 0 && OPTxs!=1 && OPTxs !=2 && OPTxs !=3 && OPTxs !=4  ) {
     202    std::ostringstream errOs;
     203    errOs << "BAD CROSS SECTION OPTION in G4PreCompoundModel.cc !!"  <<G4endl;
     204    throw G4HadronicException(__FILE__, __LINE__, errOs.str());}
     205  else aEmission.SetOPTxs(OPTxs);
     206 
     207  //for the choice of superimposed Coulomb Barrier for inverse cross sections
     208 
     209   aEmission.UseSICB(useSICB);
     210 
     211 
     212  //----------
     213 
    178214  G4VPreCompoundTransitions * aTransition = 0;
    179215  if (useGNASHTransition)
    … …  
    184220    {
    185221      aTransition = new G4PreCompoundTransitions;
    186     }
    187 
    188 
     222      // for the choice of "never go back" hypothesis and CEM transition probabilities 
     223      if (useNGB) aTransition->UseNGB(useNGB);
     224      if (useCEMtr) aTransition->UseCEMtr(useCEMtr);
     225    }
     226 
    189227  // Main loop. It is performed until equilibrium deexcitation.
     228  //G4int fragment=0;
     229 
    190230  for (;;) {
     231   
     232    //fragment++;
     233    //G4cout<<"-------------------------------------------------------------------"<<G4endl;
     234    //G4cout<<"Fragment number .. "<<fragment<<G4endl;
     235   
    191236    // Initialize fragment according with the nucleus parameters
    192237    aEmission.Initialize(aFragment);
    193 
    194     // Equilibrium exciton number
    195     //    G4EvaporationLevelDensityParameter theLDP;
    196     //    G4double g = (6.0/pi2)*aFragment.GetA()*
    197     //      theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),
    198     //                             aFragment.GetExcitationEnergy());
     238   
     239   
    199240   
    200241    G4double g = (6.0/pi2)*aFragment.GetA()*
    201242      G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    202     G4int EquilibriumExcitonNumber = static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy())
    203                                                         + 0.5);
    204    
    205     // Loop for transitions, it is performed while there are preequilibrium transitions.
     243   
     244   
     245   
     246   
     247    G4int EquilibriumExcitonNumber = static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy())+ 0.5);
     248//   
     249//    G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
     250//
     251// J. M. Quesada (Jan. 08)  equilibrium hole number could be used as preeq.- evap. delimiter (IAEA report)
     252//    G4int EquilibriumHoleNumber = static_cast<G4int>(0.2*std::sqrt(g*aFragment.GetExcitationEnergy())+ 0.5);
     253   
     254// Loop for transitions, it is performed while there are preequilibrium transitions.
    206255    G4bool ThereIsTransition = false;
     256   
     257    //        G4cout<<"----------------------------------------"<<G4endl;
     258    //        G4double NP=aFragment.GetNumberOfParticles();
     259    //        G4double NH=aFragment.GetNumberOfHoles();
     260    //        G4double NE=aFragment.GetNumberOfExcitons();
     261    //        G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
     262    //        G4cout<<"N. excitons="<<NE<<"  N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
     263       
     264   
     265    //G4int transition=0;
    207266    do
    208267      {
     268        //transition++;
     269        //G4cout<<"transition number .."<<transition<<G4endl;
     270        //G4cout<<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
    209271        G4bool go_ahead = false;
     272        // soft cutoff criterium as an "ad-hoc" solution to force increase in  evaporation 
     273        //       G4double test = static_cast<G4double>(aFragment.GetNumberOfHoles());
    210274        G4double test = static_cast<G4double>(aFragment.GetNumberOfExcitons());
    211         if (test < EquilibriumExcitonNumber)
    212           {
    213             test /= static_cast<G4double>(EquilibriumExcitonNumber);
    214             test -= 1.0;
    215             test = test*test;
    216             test /= 0.32;
    217             test = 1.0 - std::exp(-test);
    218             go_ahead = (G4UniformRand() < test);
    219           }
    220      
     275       
     276       
     277        if (test < EquilibriumExcitonNumber) go_ahead=true;
     278        //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
     279        if (useSCO) {
     280          if (test < EquilibriumExcitonNumber)
     281            //  if (test < EquilibriumHoleNumber)
     282            {
     283              test /= static_cast<G4double>(EquilibriumExcitonNumber);
     284              //     test /= static_cast<G4double>(EquilibriumHoleNumber);
     285              test -= 1.0;
     286              test = test*test;
     287              test /= 0.32;
     288              test = 1.0 - std::exp(-test);
     289              go_ahead = (G4UniformRand() < test);
     290             
     291            }
     292        }
     293       
     294        //JMQ: WARNING:  CalculateProbability MUST be called prior to Get methods !! (O values would be returned otherwise)
     295        G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);
     296        G4double P1=aTransition->GetTransitionProb1();
     297        G4double P2=aTransition->GetTransitionProb2();
     298        G4double P3=aTransition->GetTransitionProb3();
     299        //       G4cout<<"P1="<<P1<<" P2="<<P2<<"  P3="<<P3<<G4endl;
     300       
     301       
     302        //J.M. Quesada (May. 08). Physical criterium (lamdas)  PREVAILS over approximation (critical exciton number)
     303        if(P1<=(P2+P3)) go_ahead=false;
     304       
    221305        if (go_ahead &&  aFragment.GetA() > 4)
    222           {
    223             //          if (aFragment.GetNumberOfParticles() < 1) {
    224             //            aFragment.SetNumberOfHoles(aFragment.GetNumberOfHoles()+1);
    225             //            aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()+1);       
    226             //          }
    227                                
     306          {                             
    228307            G4double TotalEmissionProbability = aEmission.GetTotalProbability(aFragment);
    229        
     308            //
     309            //  G4cout<<"TotalEmissionProbability="<<TotalEmissionProbability<<G4endl;
     310            //
    230311            // Check if number of excitons is greater than 0
    231312            // else perform equilibrium emission
    … …  
    242323           
    243324            //      G4PreCompoundTransitions aTransition(aFragment);
    244        
     325           
     326            //J.M.Quesada (May 08) this has already been done in order to decide what to do (preeq-eq)
    245327            // Sum of transition probabilities
    246             G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);
    247        
     328            //      G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);
     329           
    248330            // Sum of all probabilities
    249331            G4double TotalProbability = TotalEmissionProbability + TotalTransitionProbability;
    250        
     332           
    251333            // Select subprocess
    252334            if (G4UniformRand() > TotalEmissionProbability/TotalProbability)
    253335              {
    254336                // It will be transition to state with a new number of excitons
    255                 ThereIsTransition = true;
    256                
     337                ThereIsTransition = true;               
    257338                // Perform the transition
    258339                aFragment = aTransition->PerformTransition(aFragment);
    … …  
    262343                // It will be fragment emission
    263344                ThereIsTransition = false;
    264 
    265                 // Perform the emission and Add emitted fragment to Result
    266345                Result->push_back(aEmission.PerformEmission(aFragment));
    267346              }
    … …  
    288367{
    289368  G4ReactionProductVector * theEquilibriumResult;
     369
    290370  theEquilibriumResult = GetExcitationHandler()->BreakItUp(aFragment);
    291371 
    … …  
    353433
    354434#endif
     435
     436
     437

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // by V. Lara
     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
     41//
    2742
    2843#include "G4PreCompoundNucleon.hh"
    2944#include "G4PreCompoundParameters.hh"
    3045
     46G4bool 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}
    3152
    3253G4double G4PreCompoundNucleon::
    … …  
    3657  if ( !IsItPossible(aFragment) ) return 0.0;
    3758
    38   const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
    39 
    4059  G4double U = aFragment.GetExcitationEnergy();
    4160  G4double P = aFragment.GetNumberOfParticles();
    … …  
    4362  G4double N = P + H;
    4463 
    45   // g = (6.0/pi2)*a*A
    46   //  G4EvaporationLevelDensityParameter theLDP;
    4764  G4double g0 = (6.0/pi2)*aFragment.GetA() *
    4865    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    49   //    theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);
     66 
    5067  G4double g1 = (6.0/pi2)*GetRestA() *
    5168    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    52     //    theLDP.LevelDensityParameter(static_cast<G4int>(GetRestA()),static_cast<G4int>(GetRestZ()),U);
    5369
    5470  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
    55   G4double A1 = (A0*g0 - P/2.0)/g1;
     71
     72  G4double A1 = (A0 - P/2.0)/g1;
    5673 
    5774  G4double E0 = std::max(0.0,U - A0);
    … …  
    6178
    6279
    63   G4double Probability = 2.0/(hbarc*hbarc*hbarc) * GetReducedMass() *
    64       GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) * r0 * r0 * std::pow(GetRestA(),2.0/3.0) * GetAlpha() * (eKin + GetBeta()) *
    65       P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 *
    66       g1/(g0*g0);
     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);
     83
     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  }
    6794
    6895  return Probability;
    6996}
     97
     98
     99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112

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

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

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

    r819 r962  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundTransitions.cc,v 1.13 2007/07/23 12:48:54 ahoward Exp $
    28 // GEANT4 tag $Name:  $
    29 //
    30 // by V. Lara
     26// $Id: G4PreCompoundTransitions.cc,v 1.20 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// 16.02.2008 J. M. Quesada fixed bugs
     40// 06.09.2008 J. M. Quesada added external choices for:
     41//                      - "never go back"  hipothesis (useNGB=true)
     42//                      -  CEM transition probabilities (useCEMtr=true)
    3143
    3244#include "G4PreCompoundTransitions.hh"
    … …  
    5567CalculateProbability(const G4Fragment & aFragment)
    5668{
     69  //G4cout<<"In G4PreCompoundTransitions.cc  useNGB="<<useNGB<<G4endl;
     70  //G4cout<<"In G4PreCompoundTransitions.cc  useCEMtr="<<useCEMtr<<G4endl;
     71
    5772  // Fermi energy
    5873  const G4double FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
    … …  
    7489  G4double Z = static_cast<G4double>(aFragment.GetZ());
    7590  G4double U = aFragment.GetExcitationEnergy();
    76 
    77   // Relative Energy (T_{rel})
    78   G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N;
    79  
    80   // Sample kind of nucleon-projectile
    81   G4bool ChargedNucleon(false);
    82   G4double chtest = 0.5;
    83   if (P > 0) chtest = aFragment.GetNumberOfCharged()/P;
    84   if (G4UniformRand() < chtest) ChargedNucleon = true;
    85 
    86   // Relative Velocity:
    87   // <V_{rel}>^2
    88   G4double RelativeVelocitySqr(0.0);
    89   if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2;
    90   else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2;
    91 
    92   // <V_{rel}>
    93   G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
    94 
    95   // Proton-Proton Cross Section
    96   G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn;
    97   // Proton-Neutron Cross Section
    98   G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn;
    99  
    100   // Averaged Cross Section: \sigma(V_{rel})
    101   //  G4double AveragedXSection = (ppXSection+npXSection)/2.0;
    102   G4double AveragedXSection(0.0);
    103   if (ChargedNucleon)
    104     {
    105       AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0);
    106     }
    107   else
    108     {
    109       AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0);
    110     }
    111 
    112   // Fermi relative energy ratio
    113   G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
    114 
    115   // This factor is introduced to take into account the Pauli principle
    116   G4double PauliFactor = 1.0 - (7.0/5.0)*FermiRelRatio;
    117   if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
    118 
    119   // Interaction volume
    120   G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
    121 
    122   // Transition probability for \Delta n = +2
    123   //  TransitionProb1 = 0.00332*AveragedXSection*PauliFactor*std::sqrt(RelativeEnergy)/
    124   //    std::pow(1.2 + 1.0/(4.7*RelativeVelocity), 3.0);
    125   TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint;
    126   if (TransitionProb1 < 0.0) TransitionProb1 = 0.0;
    127 
    128   // g = (6.0/pi2)*aA;
    129   //  G4double a = theLDP.LevelDensityParameter(A,Z,U-G4PairingCorrection::GetInstance()->GetPairingCorrection(A,Z));
    130   G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    131   // GE = g*E where E is Excitation Energy
    132   G4double GE = (6.0/pi2)*a*A*U;
    133 
    134 
    135   // F(p,h) = 0.25*(p^2 + h^2 + p - h) - 0.5*h
    136   G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
    137   G4bool NeverGoBack(false);
    138   //AH  if (U-Fph < 0.0) NeverGoBack = true;
    139   if (GE-Fph < 0.0) NeverGoBack = true;
    140   // F(p+1,h+1)
    141   G4double Fph1 = Fph + N/2.0;
    142   // (n+1)/n ((g*E - F(p,h))/(g*E - F(p+1,h+1)))^(n+1)
    143   G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0);
    144 
    145 
    146   if (NeverGoBack)
    147     {
     91 
     92  if(U<10*eV) return 0.0;
     93 
     94  //J. M. Quesada (Feb. 08) new physics
     95  // OPT=1 Transitions are calculated according to Gudima's paper (original in G4PreCompound from VL)
     96  // OPT=2 Transitions are calculated according to Gupta's formulae
     97  //
     98 
     99 
     100 
     101  if (useCEMtr){
     102
     103   
     104    // Relative Energy (T_{rel})
     105    G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N;
     106   
     107    // Sample kind of nucleon-projectile
     108    G4bool ChargedNucleon(false);
     109    G4double chtest = 0.5;
     110    if (P > 0) chtest = aFragment.GetNumberOfCharged()/P;
     111    if (G4UniformRand() < chtest) ChargedNucleon = true;
     112   
     113    // Relative Velocity:
     114    // <V_{rel}>^2
     115    G4double RelativeVelocitySqr(0.0);
     116    if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2;
     117    else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2;
     118   
     119    // <V_{rel}>
     120    G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
     121   
     122    // Proton-Proton Cross Section
     123    G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn;
     124    // Proton-Neutron Cross Section
     125    G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn;
     126   
     127    // Averaged Cross Section: \sigma(V_{rel})
     128    //  G4double AveragedXSection = (ppXSection+npXSection)/2.0;
     129    G4double AveragedXSection(0.0);
     130    if (ChargedNucleon)
     131      {
     132        //JMQ: small bug fixed
     133        //      AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0);
     134        AveragedXSection = ((Z-1.0) * ppXSection + (A-Z) * npXSection) / (A-1.0);
     135      }
     136    else
     137      {
     138        AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0);
     139      }
     140   
     141    // Fermi relative energy ratio
     142    G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
     143   
     144    // This factor is introduced to take into account the Pauli principle
     145    G4double PauliFactor = 1.0 - (7.0/5.0)*FermiRelRatio;
     146    if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
     147   
     148    // Interaction volume
     149    //  G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
     150    G4double xx=2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity);
     151    G4double Vint = (4.0/3.0)*pi*xx*xx*xx;
     152   
     153    // Transition probability for \Delta n = +2
     154   
     155    TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint;
     156    if (TransitionProb1 < 0.0) TransitionProb1 = 0.0;
     157   
     158    G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     159    // GE = g*E where E is Excitation Energy
     160    G4double GE = (6.0/pi2)*a*A*U;
     161   
     162    G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
     163   
     164    //G4bool NeverGoBack(false);
     165    G4bool NeverGoBack;
     166    if(useNGB)  NeverGoBack=true;
     167    else NeverGoBack=false;
     168   
     169   
     170    //JMQ/AH  bug fixed: if (U-Fph < 0.0) NeverGoBack = true;
     171    if (GE-Fph < 0.0) NeverGoBack = true;
     172   
     173    // F(p+1,h+1)
     174    G4double Fph1 = Fph + N/2.0;
     175   
     176    G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0);
     177   
     178   
     179    if (NeverGoBack)
     180      {
    148181      TransitionProb2 = 0.0;
    149182      TransitionProb3 = 0.0;
    150     }
    151   else
    152     {
    153       // Transition probability for \Delta n = -2 (at F(p,h) = 0)
    154       //  TransitionProb2 = max(0, (TransitionProb1*P*H*(P+H+1.0)*(P+H-2.0))/(GE*GE));
    155       //  TransitionProb2 = (TransitionProb1*P*H*(P+H+1.0)*(P+H-2.0))/(GE*GE);
    156       TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph));
     183      }
     184    else
     185      {
     186        // Transition probability for \Delta n = -2 (at F(p,h) = 0)
     187        TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph));
     188        if (TransitionProb2 < 0.0) TransitionProb2 = 0.0;
     189       
     190        // Transition probability for \Delta n = 0 (at F(p,h) = 0)
     191        TransitionProb3 = TransitionProb1* ((N+1.0)/N) * ProbFactor  * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph);
     192        if (TransitionProb3 < 0.0) TransitionProb3 = 0.0;
     193      }
     194   
     195    //  G4cout<<"U = "<<U<<G4endl;
     196    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
     197    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl;
     198    return TransitionProb1 + TransitionProb2 + TransitionProb3;
     199  }
     200 
     201  else  {
     202    //JMQ: Transition probabilities from Gupta's work
     203   
     204    G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     205    // GE = g*E where E is Excitation Energy
     206    G4double GE = (6.0/pi2)*a*A*U;
     207   
     208    G4double Kmfp=2.;
     209     
     210   
     211    TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
     212    if (TransitionProb1 < 0.0) TransitionProb1 = 0.0;
     213   
     214    if (useNGB){
     215      TransitionProb2=0.;
     216      TransitionProb3=0.;
     217    }
     218    else{       
     219      if (N<=1) TransitionProb2=0. ;   
     220      else  TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);     
    157221      if (TransitionProb2 < 0.0) TransitionProb2 = 0.0;
    158      
    159      
    160       // Transition probability for \Delta n = 0 (at F(p,h) = 0)
    161       //  TransitionProb3 = TransitionProb1*(P+H+1.0)*(P*(P-1.0)+4.0*P*H+H*(H-1.0))/((P+H)*GE);
    162       TransitionProb3 = TransitionProb1 * ProbFactor * ((N+1.0)/N) * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph);
    163       if (TransitionProb3 < 0.0) TransitionProb3 = 0.0;
    164     }
    165  
    166   return TransitionProb1 + TransitionProb2 + TransitionProb3;
     222      TransitionProb3=0.;
     223    }
     224   
     225      //  G4cout<<"U = "<<U<<G4endl;
     226    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
     227    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl;
     228    return TransitionProb1 + TransitionProb2 + TransitionProb3;
     229  }
    167230}
    168231
    … …  
    185248    }
    186249
    187   // 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
    188   if(deltaN < 0 && G4UniformRand() <= static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) && (result.GetNumberOfCharged() >= 1))
     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)) {
    189255    result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative!
    190 
    191   // result.SetNumberOfParticles was here
    192   // result.SetNumberOfHoles was here
    193   // the following lines have to be before SetNumberOfCharged, otherwise the check on number of charged vs. number of particles fails
     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
    194260  result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
    195261  result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2);
    196262
    197   // With weight Z/A, number of charged particles is decreased on +1
    198   //  if ((deltaN > 0 || result.GetNumberOfCharged() > 0) && // AH/JMQ check is now in initialize within G4VPreCompoundFragment
     263  // With weight Z/A, number of charged particles is increased with +1
    199264  if ( ( deltaN > 0 ) &&
    200265      (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/
    201                   std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
     266       std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
    202267    {
    203268      result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2);
    … …  
    210275    }
    211276 
    212   // moved from above to make code more readable
    213   //  result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
    214   //  result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2);
    215 
    216277  return result;
    217278}

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

    r819 r962  
    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 $
    2628//
    27 // $Id: G4VPreCompoundFragment.cc,v 1.4 2007/07/23 09:56:40 ahoward Exp $
    28 // GEANT4 tag $Name:  $
     29// J. M. Quesada (August 2008). 
     30// Based  on previous work by V. Lara
    2931//
    30 // by V. Lara
    3132 
    3233#include "G4VPreCompoundFragment.hh"
    … …  
    142143  if ((theRestNucleusA < theRestNucleusZ) ||
    143144      (theRestNucleusA < theA) ||
    144       (theRestNucleusZ < theZ) ||
    145       (aFragment.GetNumberOfCharged() < theZ)) // AH last argument from JMQ
     145      (theRestNucleusZ < theZ))
    146146    {
    147147      // In order to be sure that emission probability will be 0.
    … …  
    155155    GetCoulombBarrier(static_cast<G4int>(theRestNucleusA),static_cast<G4int>(theRestNucleusZ),
    156156                      aFragment.GetExcitationEnergy());
    157  
     157
     158
    158159  // Compute Binding Energies for fragments
    159160  // (needed to separate a fragment from the nucleus)
    … …  
    164165 
    165166  // Compute Maximal Kinetic Energy which can be carried by fragments after separation
     167  // This is the true (assimptotic) maximal kinetic energy
    166168  G4double m = aFragment.GetMomentum().m();
    167169  G4double rm = GetRestNuclearMass();
    168170  G4double em = GetNuclearMass();
    169171  theMaximalKineticEnergy = ((m - rm)*(m + rm) + em*em)/(2.0*m) - em;
    170   // - theCoulombBarrier;
     172 
    171173 
    172174  return;

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4VPreCompoundIon.cc,v 1.5 2007/07/23 09:56:40 ahoward Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4VPreCompoundIon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// by V. Lara
    31 
    32 #include "G4VPreCompoundIon.hh"
    33 #include "G4PreCompoundParameters.hh"
    34 //#include "G4EvaporationLevelDensityParameter.hh"
    35 
    36 
    37 G4double G4VPreCompoundIon::
    38 ProbabilityDistributionFunction(const G4double eKin,
    39                                 const G4Fragment& aFragment)
    40 {
    41   if ( !IsItPossible(aFragment) ) return 0.0;
    42 
    43   const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
    44 
    45   G4double U = aFragment.GetExcitationEnergy();
    46   G4double P = aFragment.GetNumberOfParticles();
    47   G4double H = aFragment.GetNumberOfHoles();
    48   G4double N = P + H;
    49 
    50   // g = (6.0/pi2)*a*A
    51   //  G4EvaporationLevelDensityParameter theLDP;
    52   G4double g0 = (6.0/pi2)*aFragment.GetA() *
    53     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    54   //    theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);
    55   G4double g1 = (6.0/pi2)*GetRestA() *
    56     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    57   //    theLDP.LevelDensityParameter(GetRestA(),GetRestZ(),U);
    58   G4double gj = (6.0/pi2)*GetA() *
    59     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    60   //    theLDP.LevelDensityParameter(GetA(),GetZ(),U);
    61 
    62   G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
    63   G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1);
    64   G4double Aj = GetA()*(GetA()+1.0)/4.0;
    65 
    66   G4double E0 = std::max(0.0,U - A0);
    67   if (E0 == 0.0) return 0.0;
    68   G4double E1 = std::max(0.0,GetMaximalKineticEnergy() + GetCoulombBarrier() - eKin - A1);
    69   G4double Ej = std::max(0.0,U + GetBindingEnergy() - Aj);
    70 
    71   G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*(eKin+GetBindingEnergy()))))*
    72       GetAlpha()*(eKin+GetBeta())/(r0*std::pow(GetRestA(),1.0/3.0)) *
    73       CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P);
    74 
    75   G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0);
    76 
    77   G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0;
    78 
    79   G4double Probability = pA * pB * pC;
    80 
    81   return Probability;
    82 }
     31//
     32//J.M. Quesada (Apr. 2008) DUMMY abstract base class for ions
     33// it is not ihherited by anything
     34// Suppressed
    8335
    8436
    … …  
    8840
    8941
     42
     43
     44
     45
     46

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4VPreCompoundNucleon.cc,v 1.5 2007/07/23 09:56:40 ahoward Exp $
    28 // GEANT4 tag $Name:  $
     27
     28// $Id: G4VPreCompoundNucleon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $
     29// GEANT4 tag $Name: geant4-09-02-ref-02 $
     30
    2931//
    3032// by V. Lara
    31 
    32 #include "G4VPreCompoundNucleon.hh"
    33 #include "G4PreCompoundParameters.hh"
    34 //#include "G4EvaporationLevelDensityParameter.hh"
     33//
     34//J.M. Quesada (Apr. 2008) DUMMY abstract base class for nucleons.
     35// it is not ihherited by anything
     36// Suppressed
    3537
    3638
    37 G4double G4VPreCompoundNucleon::
    38 ProbabilityDistributionFunction(const G4double eKin,
    39                                 const G4Fragment& aFragment)
    40 {
    41   if ( !IsItPossible(aFragment) ) return 0.0;
    42 
    43   const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
    44 
    45   G4double U = aFragment.GetExcitationEnergy();
    46   G4double P = aFragment.GetNumberOfParticles();
    47   G4double H = aFragment.GetNumberOfHoles();
    48   G4double N = P + H;
    49 
    50   // g = (6.0/pi2)*a*A
    51   //  G4EvaporationLevelDensityParameter theLDP;
    52   G4double g0 = (6.0/pi2)*aFragment.GetA() *
    53     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    54   //    theLDP.LevelDensityParameter(G4int(aFragment.GetA()),G4int(aFragment.GetZ()),U);
    55   G4double g1 = (6.0/pi2)*GetRestA() *
    56     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
    57     //    theLDP.LevelDensityParameter(G4int(GetRestA()),G4int(GetRestZ()),U);
    58 
    59 
    60   G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
    61   G4double A1 = (A0*g0 - P/2.0)/g1;
    62 
    63   G4double E0 = std::max(0.0,U - A0);
    64   if (E0 == 0.0) return 0.0;
    65   G4double E1 = std::max(0.0,U - eKin - GetBindingEnergy() - A1);
    66 
    67   G4double Probability = 2.0/(hbarc*hbarc*hbarc) * GetReducedMass() *
    68       GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) * r0 * r0 * std::pow(GetRestA(),2.0/3.0) * GetAlpha() * (eKin + GetBeta()) *
    69       P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 *
    70       g1/(g0*g0);
    71 
    72   return Probability;
    73 }
    74 
    75 
Note: See TracChangeset for help on using the changeset viewer.