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

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4EvaporationChannel.cc,v 1.5 2006/06/29 20:10:27 gunter Exp $
    28 // GEANT4 tag $Name:  $
     27//J.M. Quesada (August2008). Based on:
    2928//
    3029// Hadronic Process: Nuclear De-excitations
    3130// by V. Lara (Oct 1998)
    3231//
     32// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     33// cross section option
     34// JMQ (06 September 2008) Also external choices have been added for
     35// superimposed Coulomb barrier (if useSICB is set true, by default is false)
     36
    3337
    3438#include "G4EvaporationChannel.hh"
    … …  
    3640
    3741
    38 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ,
     42
     43G4EvaporationChannel::G4EvaporationChannel(const G4int anA, const G4int aZ, const G4String & aName,
    3944                                           G4VEmissionProbability * aEmissionStrategy,
    40                                            G4VCoulombBarrier * aCoulombBarrier):
    41     A(theA),
    42     Z(theZ),
     45                                           G4VCoulombBarrier * aCoulombBarrier):
     46    G4VEvaporationChannel(aName),
     47    theA(anA),
     48    theZ(aZ),
    4349    theEvaporationProbabilityPtr(aEmissionStrategy),
    4450    theCoulombBarrierPtr(aCoulombBarrier),
    … …  
    4753{
    4854    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    49     MyOwnLevelDensity = true;
    50 }
    51 
    52 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String & aName,
    53                                            G4VEmissionProbability * aEmissionStrategy,
    54                                            G4VCoulombBarrier * aCoulombBarrier):
    55     G4VEvaporationChannel(aName),
    56     A(theA),
    57     Z(theZ),
    58     theEvaporationProbabilityPtr(aEmissionStrategy),
    59     theCoulombBarrierPtr(aCoulombBarrier),
    60     EmissionProbability(0.0),
    61     MaximalKineticEnergy(-1000.0)
    62 {
    63     theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    64     MyOwnLevelDensity = true;
    65 }
    66 
    67 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String * aName,
    68                                            G4VEmissionProbability * aEmissionStrategy,
    69                                            G4VCoulombBarrier * aCoulombBarrier):
    70     G4VEvaporationChannel(aName),
    71     A(theA),
    72     Z(theZ),
    73     theEvaporationProbabilityPtr(aEmissionStrategy),
    74     theCoulombBarrierPtr(aCoulombBarrier),
    75     EmissionProbability(0.0),
    76     MaximalKineticEnergy(-1000.0)
    77 {
    78     theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    79     MyOwnLevelDensity = true;
    80 }
    81 
     55//    MyOwnLevelDensity = true; 
     56}
    8257
    8358G4EvaporationChannel::~G4EvaporationChannel()
    8459{
    85     if (MyOwnLevelDensity) delete theLevelDensityPtr;
    86 }
    87 
    88 
    89 
     60  delete theLevelDensityPtr;
     61}
    9062
    9163G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel()
    … …  
    11284}
    11385
    114 
     86//void G4EvaporationChannel::SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)
     87//  {
     88//    if (MyOwnLevelDensity) delete theLevelDensityPtr;
     89//    theLevelDensityPtr = aLevelDensity;
     90//    MyOwnLevelDensity = false;
     91//  }
    11592
    11693void G4EvaporationChannel::Initialize(const G4Fragment & fragment)
    11794{
    118 
    119     G4int anA = static_cast<G4int>(fragment.GetA());
    120     G4int aZ = static_cast<G4int>(fragment.GetZ());
    121     AResidual = anA - A;
    122     ZResidual = aZ - Z;
    123 
    124     // Effective excitation energy
    125     G4double ExEnergy = fragment.GetExcitationEnergy() -
    126       G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
    127 
    128     // We only take into account channels which are physically allowed
    129     if (AResidual <= 0 || ZResidual <= 0 || AResidual < ZResidual ||
    130         (AResidual == ZResidual && AResidual > 1) || ExEnergy <= 0.0) {
    131         //              LevelDensityParameter = 0.0;
    132         CoulombBarrier = 0.0;
    133         //              BindingEnergy = 0.0;
    134         MaximalKineticEnergy = -1000.0*MeV;
    135         EmissionProbability = 0.0;
    136     } else {
    137         // // Get Level Density
    138         // LevelDensityParameter = theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
    139 
    140         // Coulomb Barrier calculation
    141         CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(AResidual,ZResidual,ExEnergy);
    142        
    143         // // Binding Enegy (for separate fragment from nucleus)
    144         // BindingEnergy = CalcBindingEnergy(anA,aZ);
    145 
    146         // Maximal Kinetic Energy
    147         MaximalKineticEnergy = CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()->
    148                                                         GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy);
    149                
    150         // Emission probability
    151         if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0;
    152         else {
    153             // Total emission probability for this channel
    154             EmissionProbability = theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
    155 
    156         }
     95  //for inverse cross section choice
     96  theEvaporationProbabilityPtr->SetOPTxs(OPTxs);
     97  // for superimposed Coulomb Barrier for inverse cross sections
     98  theEvaporationProbabilityPtr->UseSICB(useSICB);
     99 
     100 
     101  G4int FragmentA = static_cast<G4int>(fragment.GetA());
     102  G4int FragmentZ = static_cast<G4int>(fragment.GetZ());
     103  ResidualA = FragmentA - theA;
     104  ResidualZ = FragmentZ - theZ;
     105 
     106  //Effective excitation energy
     107  G4double ExEnergy = fragment.GetExcitationEnergy() -
     108    G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ);
     109 
     110 
     111  // Only channels which are physically allowed are taken into account
     112  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
     113      (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) {
     114    CoulombBarrier=0.0;
     115    MaximalKineticEnergy = -1000.0*MeV;
     116    EmissionProbability = 0.0;
     117  } else {
     118    CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
     119    // Maximal Kinetic Energy
     120    MaximalKineticEnergy = CalcMaximalKineticEnergy
     121      (G4ParticleTable::GetParticleTable()->
     122       GetIonTable()->GetNucleusMass(FragmentZ,FragmentA)+ExEnergy);
     123   
     124    // Emission probability
     125    // Protection for the case Tmax<V. If not set in this way we could end up in an
     126    // infinite loop in  the method GetKineticEnergy if OPTxs!=0 && useSICB=true.
     127    // Of course for OPTxs=0 we have the Coulomb barrier
     128   
     129    G4double limit;
     130    if (OPTxs==0 || (OPTxs!=0 && useSICB))
     131      limit= CoulombBarrier;
     132    else limit=0.;
     133 
     134    // The threshold for charged particle emission must be  set to 0 if Coulomb
     135    //cutoff  is included in the cross sections
     136    //if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0; 
     137    if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0;
     138    else {
     139      // Total emission probability for this channel
     140      EmissionProbability = theEvaporationProbabilityPtr->
     141        EmissionProbability(fragment, MaximalKineticEnergy);
    157142    }
    158        
    159     return;
    160 }
    161 
     143  }
     144 
     145  return;
     146}
    162147
    163148G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus)
    164149{
    165 
    166     G4double EvaporatedKineticEnergy = CalcKineticEnergy();
    167     G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
    168     G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
    169 
    170     G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
    171                                                 (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
    172  
    173     momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
    174 
    175     G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
    176     EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
    177 
    178     G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
     150  G4double EvaporatedKineticEnergy=GetKineticEnergy(theNucleus);
     151 
     152  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
     153    GetIonMass(theZ,theA);
     154  G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
     155 
     156  G4ThreeVector momentum(IsotropicVector
     157                         (std::sqrt(EvaporatedKineticEnergy*
     158                                    (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
     159 
     160  momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
     161 
     162  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
     163  EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
     164 
     165  G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum);
    179166#ifdef PRECOMPOUND_TEST
    180     EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
     167  EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
    181168#endif
    182     // ** And now the residual nucleus **
    183     G4double theExEnergy = theNucleus.GetExcitationEnergy();
    184     G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
    185         GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
    186                        static_cast<G4int>(theNucleus.GetA()));
    187     G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
    188        
    189     G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
    190     ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
    191        
    192     G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum );
     169  // ** And now the residual nucleus **
     170  G4double theExEnergy = theNucleus.GetExcitationEnergy();
     171  G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
     172    GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
     173                   static_cast<G4int>(theNucleus.GetA()));
     174  G4double ResidualEnergy = theMass +
     175    (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
     176 
     177  G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
     178  ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
     179 
     180  G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum );
     181 
    193182#ifdef PRECOMPOUND_TEST
    194     ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
     183  ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
    195184#endif
    196     G4FragmentVector * theResult = new G4FragmentVector;
    197 
     185  G4FragmentVector * theResult = new G4FragmentVector;
     186 
    198187#ifdef debug
    199     G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
    200     G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
    201     if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
    202         G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
    203         G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :"
    204                <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
    205                << " MeV" << G4endl;
    206     }
    207     if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
    208         std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
    209         std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
    210         G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
    211         G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :"
    212                <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
    213                << " MeV" << G4endl;   
    214 
    215     }
     188  G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
     189  G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
     190  if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
     191    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
     192    G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :"
     193           <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
     194           << " MeV" << G4endl;
     195  }
     196  if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
     197      std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
     198      std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
     199    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
     200    G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :"
     201           <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
     202           << " MeV" << G4endl;   
     203   
     204  }
    216205#endif
    217     theResult->push_back(EvaporatedFragment);
    218     theResult->push_back(ResidualFragment);
    219     return theResult;
     206  theResult->push_back(EvaporatedFragment);
     207  theResult->push_back(ResidualFragment);
     208  return theResult;
    220209}
    221210
    222 
    223 
    224 // G4double G4EvaporationChannel::CalcBindingEnergy(const G4int anA, const G4int aZ)
    225 // // Calculate Binding Energy for separate fragment from nucleus
    226 // {
    227 //      // Mass Excess for residual nucleus
    228 //      G4double ResNucMassExcess = G4NucleiProperties::GetNuclearMass(AResidual,ZResidual);
    229 //      // Mass Excess for fragment
    230 //      G4double FragmentMassExcess = G4NucleiProperties::GetNuclearMass(A,Z);
    231 //      // Mass Excess for Compound Nucleus
    232 //      G4double NucleusMassExcess = G4NucleiProperties::GetNuclearMass(anA,aZ);
    233 //
    234 //      return ResNucMassExcess + FragmentMassExcess - NucleusMassExcess;
    235 // }
    236 
    237 
     211/////////////////////////////////////////
     212// Calculates the maximal kinetic energy that can be carried by fragment.
    238213G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
    239     // Calculate maximal kinetic energy that can be carried by fragment.
    240 {
    241     G4double ResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( ZResidual, AResidual );
    242     G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( Z, A );
    243        
    244     G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
    245         (2.0*NucleusTotalE) -
    246         EvaporatedMass - CoulombBarrier;
    247        
    248     return T;
    249 }
    250 
    251 
    252 
    253 
    254 G4double G4EvaporationChannel::CalcKineticEnergy(void)
    255     // Samples fragment kinetic energy.
     214{
     215  G4double ResidualMass = G4ParticleTable::GetParticleTable()->
     216    GetIonTable()->GetNucleusMass( ResidualZ, ResidualA );
     217  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
     218    GetIonTable()->GetNucleusMass( theZ, theA );
     219
     220  // This is the "true" assimptotic kinetic energy (from energy conservation)   
     221  G4double Tmax = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass -               
     222                   ResidualMass*ResidualMass)/(2.0*NucleusTotalE) - EvaporatedMass;
     223 
     224  //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated
     225  //at the Coulomb barrier
     226  //IMPORTANT: meaning of Tmax differs in OPTxs=0 and OPTxs!=0
     227  //When OPTxs!=0 Tmax is the TRUE (assimptotic) maximal kinetic energy
     228 
     229  if(OPTxs==0)
     230    Tmax=Tmax- CoulombBarrier;
     231 
     232  return Tmax;
     233}
     234
     235///////////////////////////////////////////
     236//JMQ: New method for MC sampling of kinetic energy. Substitutes old CalcKineticEnergy
     237G4double G4EvaporationChannel::GetKineticEnergy(const G4Fragment & aFragment)
     238{
     239 
     240  if (OPTxs==0) {
    256241    // It uses Dostrovsky's approximation for the inverse reaction cross
    257     // in the probability for fragment emisson
    258 {
    259     if (MaximalKineticEnergy < 0.0)
    260         throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic energy is less than 0");
    261 
    262     G4double Rb = 4.0*theLevelDensityPtr->LevelDensityParameter(AResidual+A,ZResidual+Z,MaximalKineticEnergy)*
    263         MaximalKineticEnergy;
     242    // in the probability for fragment emission
     243    // MaximalKineticEnergy energy in the original version (V.Lara) was calculated at
     244    //the Coulomb barrier.
     245   
     246   
     247    if (MaximalKineticEnergy < 0.0)   
     248      throw G4HadronicException(__FILE__, __LINE__,
     249                                "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0");
     250   
     251    G4double Rb = 4.0*theLevelDensityPtr->
     252      LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)*
     253      MaximalKineticEnergy;
    264254    G4double RbSqrt = std::sqrt(Rb);
    265255    G4double PEX1 = 0.0;
    … …  
    268258    G4double FRk = 0.0;
    269259    do {
    270         G4double RandNumber = G4UniformRand();
    271         Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
    272         G4double Q1 = 1.0;
    273         G4double Q2 = 1.0;
    274         if (Z == 0) { // for emitted neutron
    275             G4double Beta = (2.12/std::pow(G4double(AResidual),2./3.) - 0.05)*MeV/
    276                 (0.76 + 2.2/std::pow(G4double(AResidual),1./3.));
    277             Q1 = 1.0 + Beta/(MaximalKineticEnergy);
    278             Q2 = Q1*std::sqrt(Q1);
    279         }
    280    
    281         FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
    282    
     260      G4double RandNumber = G4UniformRand();
     261      Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
     262      G4double Q1 = 1.0;
     263      G4double Q2 = 1.0;
     264      if (theZ == 0) { // for emitted neutron
     265        G4double Beta = (2.12/std::pow(G4double(ResidualA),2./3.) - 0.05)*MeV/
     266          (0.76 + 2.2/std::pow(G4double(ResidualA),1./3.));
     267        Q1 = 1.0 + Beta/(MaximalKineticEnergy);
     268        Q2 = Q1*std::sqrt(Q1);
     269      }
     270     
     271      FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
     272     
    283273    } while (FRk < G4UniformRand());
    284 
     274   
    285275    G4double result =  MaximalKineticEnergy * (1.0-Rk*Rk) + CoulombBarrier;
    286 
    287276    return result;
    288 }
    289  
     277  } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) {
     278   
     279    G4double V;
     280    if(useSICB) V= CoulombBarrier;
     281    else V=0.;
     282    //If Coulomb barrier is just included  in the cross sections
     283    //  G4double V=0.;
     284
     285    G4double Tmax=MaximalKineticEnergy;
     286    G4double T(0.0);
     287    G4double NormalizedProbability(1.0);
     288   
     289    // A pointer is created in order to access the distribution function.
     290    G4EvaporationProbability * G4EPtemp;
     291   
     292    if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability();
     293    else if (theA==1 && theZ==1) G4EPtemp=new G4ProtonEvaporationProbability();
     294    else if (theA==2 && theZ==1 ) G4EPtemp=new G4DeuteronEvaporationProbability();
     295    else if (theA==3 && theZ==1 ) G4EPtemp=new G4TritonEvaporationProbability();
     296    else if (theA==3 && theZ==2 ) G4EPtemp=new G4He3EvaporationProbability();
     297    else if (theA==4 && theZ==2) G4EPtemp=new G4AlphaEvaporationProbability();
     298    else {
     299      std::ostringstream errOs;
     300      errOs << "ejected particle out of range in G4EvaporationChannel"  << G4endl;
     301      throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     302    }
     303
     304      //for cross section selection and superimposed Coulom Barrier for xs
     305      G4EPtemp->SetOPTxs(OPTxs);
     306      G4EPtemp->UseSICB(useSICB);
     307
     308    do
     309      { 
     310        T=V+G4UniformRand()*(Tmax-V);
     311        NormalizedProbability=G4EPtemp->ProbabilityDistributionFunction(aFragment,T)/
     312          (this->GetEmissionProbability());
     313       
     314      }
     315    while (G4UniformRand() > NormalizedProbability);
     316   delete G4EPtemp;
     317   return T;
     318  } else{
     319    std::ostringstream errOs;
     320    errOs << "Bad option for energy sampling in evaporation"  <<G4endl;
     321    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     322  }
     323}
    290324
    291325G4ThreeVector G4EvaporationChannel::IsotropicVector(const G4double Magnitude)
    … …  
    300334                         Magnitude*CosTheta);
    301335    return Vector;
    302 }
    303 
    304 
    305 
     336            }
     337
     338
     339
     340
     341
     342   
     343
     344
     345 
     346
Note: See TracChangeset for help on using the changeset viewer.