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

geant4 tag 9.4

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc

    r1315 r1347  
    4444#include "G4IonTable.hh"
    4545
    46 
    47 G4EvaporationProbability::G4EvaporationProbability(const G4EvaporationProbability &) : G4VEmissionProbability()
     46G4EvaporationProbability::G4EvaporationProbability(G4int anA, G4int aZ,
     47                                                   G4double aGamma,
     48                                                   G4VCoulombBarrier * aCoulombBarrier)
     49  : theA(anA),
     50    theZ(aZ),
     51    Gamma(aGamma),
     52    theCoulombBarrierptr(aCoulombBarrier)
     53{}
     54
     55G4EvaporationProbability::G4EvaporationProbability()
     56  : theA(0),
     57    theZ(0),
     58    Gamma(0.0),
     59    theCoulombBarrierptr(0)
     60{}
     61
     62G4EvaporationProbability::~G4EvaporationProbability()
     63{}
     64 
     65G4double
     66G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, G4double anEnergy)
    4867{
    49     throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationProbability::copy_constructor meant to not be accessable");
    50 }
    51 
    52 
    53 
    54 
    55 const G4EvaporationProbability & G4EvaporationProbability::
    56 operator=(const G4EvaporationProbability &)
     68  G4double EmissionProbability = 0.0;
     69  G4double MaximalKineticEnergy = anEnergy;
     70
     71  if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
     72    EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy);
     73
     74  }
     75  return EmissionProbability;
     76}
     77
     78////////////////////////////////////
     79
     80// Computes the integrated probability of evaporation channel
     81G4double
     82G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment,
     83                                               G4double MaximalKineticEnergy)
    5784{
    58     throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationProbability::operator= meant to not be accessable");
    59     return *this;
    60 }
    61 
    62 
    63 G4bool G4EvaporationProbability::operator==(const G4EvaporationProbability &) const
    64 {
    65     return false;
    66 }
    67 
    68 G4bool G4EvaporationProbability::operator!=(const G4EvaporationProbability &) const
    69 {
    70     return true;
    71 }
    72  
    73 G4double G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, const G4double anEnergy)
    74 {
    75     G4double EmissionProbability = 0.0;
    76     G4double MaximalKineticEnergy = anEnergy;
    77 
    78     if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
    79         EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy);
    80 
    81     }
    82     return EmissionProbability;
    83 }
    84 
    85 ////////////////////////////////////
    86 
    87 // Computes the integrated probability of evaporation channel
    88 G4double G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy)
    89 {
    90     G4double ResidualA = fragment.GetA() - theA;
    91     G4double ResidualZ = fragment.GetZ() - theZ;
    92     G4double U = fragment.GetExcitationEnergy();
     85  G4int ResidualA = fragment.GetA_asInt() - theA;
     86  G4int ResidualZ = fragment.GetZ_asInt() - theZ;
     87  G4double U = fragment.GetExcitationEnergy();
    9388   
    94  if (OPTxs==0) {
    95 
    96        
    97     G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA);
    98 
    99 
    100     G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),
    101                                                                                static_cast<G4int>(fragment.GetZ()));
    102 
    103     G4double SystemEntropy = 2.0*std::sqrt(theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
    104                                                                            static_cast<G4int>(fragment.GetZ()),U)*
    105                                       (U-delta0));
     89  if (OPTxs==0) {
     90
     91    G4double NuclearMass = fragment.ComputeGroundStateMass(theZ,theA);
     92
     93    G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(),
     94                                                      fragment.GetZ_asInt());
     95
     96    G4double SystemEntropy = 2.0*std::sqrt(
     97      theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(),fragment.GetZ_asInt(),U)*
     98      (U-delta0));
    10699                                                                 
    107 
    108     G4double RN = 1.5*fermi;
     100    const G4double RN = 1.5*fermi;
    109101
    110102    G4double Alpha = CalcAlphaParam(fragment);
    … …  
    112104       
    113105    G4double Rmax = MaximalKineticEnergy;
    114     G4double a = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),
    115                                                       static_cast<G4int>(ResidualZ),
    116                                                       Rmax);
    117     G4double GlobalFactor = static_cast<G4double>(Gamma) * (Alpha/(a*a)) *
    118         (NuclearMass*RN*RN*std::pow(ResidualA,2./3.))/
    119         (2.*pi* hbar_Planck*hbar_Planck);
     106    G4double a = theEvapLDPptr->LevelDensityParameter(ResidualA,ResidualZ,Rmax);
     107    G4double GlobalFactor = Gamma * Alpha/(a*a) *
     108        (NuclearMass*RN*RN*fG4pow->Z23(ResidualA))/
     109        (twopi* hbar_Planck*hbar_Planck);
    120110    G4double Term1 = (2.0*Beta*a-3.0)/2.0 + Rmax*a;
    121111    G4double Term2 = (2.0*Beta*a-3.0)*std::sqrt(Rmax*a) + 2.0*a*Rmax;
    122112       
    123113    G4double ExpTerm1 = 0.0;
    124     if (SystemEntropy <= 600.0) ExpTerm1 = std::exp(-SystemEntropy);
     114    if (SystemEntropy <= 600.0) { ExpTerm1 = std::exp(-SystemEntropy); }
    125115       
    126116    G4double ExpTerm2 = 2.*std::sqrt(a*Rmax) - SystemEntropy;
    127     if (ExpTerm2 > 700.0) ExpTerm2 = 700.0;
     117    if (ExpTerm2 > 700.0) { ExpTerm2 = 700.0; }
    128118    ExpTerm2 = std::exp(ExpTerm2);
    129119       
    … …  
    134124 } else if (OPTxs==1 || OPTxs==2 ||OPTxs==3 || OPTxs==4) {
    135125
    136    G4double limit;
    137    if (useSICB) limit=theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U);
    138    else limit=0.;
    139 
    140    if (MaximalKineticEnergy <= limit)  return 0.0;
    141 
    142 
    143    // if Coulomb barrier cutoff is superimposed for all cross sections the limit is the Coulomb Barrier
     126   G4double EvaporatedMass = fragment.ComputeGroundStateMass(theZ,theA);
     127   G4double ResidulalMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA);
     128   G4double limit = std::max(0.0,fragment.GetGroundStateMass()-EvaporatedMass-ResidulalMass);
     129   if (useSICB) {
     130     limit = std::max(limit,theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U));
     131   }
     132
     133   if (MaximalKineticEnergy <= limit) { return 0.0; }
     134
     135   // if Coulomb barrier cutoff is superimposed for all cross sections
     136   // then the limit is the Coulomb Barrier
    144137   G4double LowerLimit= limit;
    145138
    146    //  MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel)     
     139   //MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel)     
    147140
    148141   G4double UpperLimit = MaximalKineticEnergy;
    149142
    150 
    151143   G4double Width = IntegrateEmissionProbability(fragment,LowerLimit,UpperLimit);
    152144
    153145   return Width;
    154  } else{
     146 } else {
    155147   std::ostringstream errOs;
    156148   errOs << "Bad option for cross sections at evaporation"  <<G4endl;
    … …  
    163155
    164156G4double G4EvaporationProbability::
    165 IntegrateEmissionProbability(const G4Fragment & fragment, const G4double & Low, const G4double & Up )
     157IntegrateEmissionProbability(const G4Fragment & fragment,
     158                             const G4double & Low, const G4double & Up )
    166159{
    167 
    168160  static const G4int N = 10;
    169161  // 10-Points Gauss-Legendre abcisas and weights
    … …  
    212204//New method (OPT=1,2,3,4)
    213205
    214 G4double G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment, const G4double K)
     206G4double
     207G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment,
     208                                                           G4double K)
    215209{
    216 
     210  G4int ResidualA = fragment.GetA_asInt() - theA;
     211  G4int ResidualZ = fragment.GetZ_asInt() - theZ; 
     212  G4double U = fragment.GetExcitationEnergy();
     213  //G4cout << "### G4EvaporationProbability::ProbabilityDistributionFunction" << G4endl;
     214  //G4cout << "FragZ= " << fragment.GetZ_asInt() << " FragA= " << fragment.GetA_asInt()
     215  //     << " Z= " << theZ << "  A= " << theA << G4endl;
     216  //G4cout << "PC " << fPairCorr << "   DP " << theEvapLDPptr << G4endl;
     217
     218  // if(K <= theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U)) return 0.0;
     219
     220  G4double delta1 = fPairCorr->GetPairingCorrection(ResidualA,ResidualZ);
    217221 
    218 
    219 
    220   G4double ResidualA = fragment.GetA() - theA;
    221   G4double ResidualZ = fragment.GetZ() - theZ; 
    222   G4double U = fragment.GetExcitationEnergy();
    223 
    224   //        if(K <= theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return 0.0;   
    225 
    226   G4double delta1 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ));
    227 
    228  
    229   G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));
    230 
    231  
    232   G4double ParticleMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA);
    233 
    234   G4double theSeparationEnergy= G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),static_cast<G4int>(theZ)) +
    235     G4NucleiProperties::GetMassExcess(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ)) -
    236     G4NucleiProperties::GetMassExcess(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));
    237 
    238   G4double a0 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()),U - delta0);
    239 
    240   G4double a1 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ),U - theSeparationEnergy - delta1);
    241  
    242  
    243   G4double E0=U-delta0;
    244 
    245   G4double E1=U-theSeparationEnergy-delta1-K;
    246 
    247   if (E1<0.) return 0.;
     222  G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(),
     223                                                    fragment.GetZ_asInt());
     224
     225 
     226  G4double ParticleMass = fragment.ComputeGroundStateMass(theZ,theA);
     227  G4double ResidualMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA);
     228
     229  G4double theSeparationEnergy = ParticleMass + ResidualMass
     230    - fragment.GetGroundStateMass();
     231
     232  G4double a0 = theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(),
     233                                                     fragment.GetZ_asInt(),
     234                                                     U - delta0);
     235
     236  G4double a1 = theEvapLDPptr->LevelDensityParameter(ResidualA, ResidualZ,
     237                                                     U - theSeparationEnergy - delta1);
     238 
     239 
     240  G4double E0 = U - delta0;
     241
     242  G4double E1 = U - theSeparationEnergy - delta1 - K;
     243
     244  if (E1<0.) { return 0.; }
    248245
    249246  //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead of hbar_Planck
    250247  //Without 1/hbar_Panck remains as a width
    251   //  G4double  Prob=Gamma*ParticleMass/((pi*hbar_Planck)*(pi*hbar_Planck)*
    252   //std::exp(2*std::sqrt(a0*E0)))*K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn;
    253248
    254249  G4double Prob=Gamma*ParticleMass/((pi*hbarc)*(pi*hbarc)*std::exp(2*std::sqrt(a0*E0)))
Note: See TracChangeset for help on using the changeset viewer.