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

geant4 tag 9.4

Location:
trunk/source/processes/hadronic/models/de_excitation/fission/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4CompetitiveFission.cc,v 1.12 2010/05/03 16:49:19 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4CompetitiveFission.cc,v 1.14 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    3939#include "G4PairingCorrection.hh"
    4040#include "G4ParticleMomentum.hh"
     41#include "G4Pow.hh"
    4142
    4243G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission")
    … …  
    5758}
    5859
    59 G4CompetitiveFission::G4CompetitiveFission(const G4CompetitiveFission &) : G4VEvaporationChannel()
    60 {
    61 }
    62 
    6360G4CompetitiveFission::~G4CompetitiveFission()
    6461{
    … …  
    7067}
    7168
    72 const G4CompetitiveFission & G4CompetitiveFission::operator=(const G4CompetitiveFission &)
    73 {
    74     throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::operator= meant to not be accessable");
    75     return *this;
    76 }
    77 
    78 G4bool G4CompetitiveFission::operator==(const G4CompetitiveFission &right) const
    79 {
    80     return (this == (G4CompetitiveFission *) &right);
    81 }
    82 
    83 G4bool G4CompetitiveFission::operator!=(const G4CompetitiveFission &right) const
    84 {
    85     return (this != (G4CompetitiveFission *) &right);
    86 }
    87 
    88 
    89 
    90 
    9169void G4CompetitiveFission::Initialize(const G4Fragment & fragment)
    9270{
    93     G4int anA = static_cast<G4int>(fragment.GetA());
    94     G4int aZ = static_cast<G4int>(fragment.GetZ());
    95     G4double ExEnergy = fragment.GetExcitationEnergy() -
    96            G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(anA,aZ);
    97  
    98 
    99     // Saddle point excitation energy ---> A = 65
    100     // Fission is excluded for A < 65
    101     if (anA >= 65 && ExEnergy > 0.0) {
    102         FissionBarrier = theFissionBarrierPtr->FissionBarrier(anA,aZ,ExEnergy);
    103         MaximalKineticEnergy = ExEnergy - FissionBarrier;
    104         LevelDensityParameter = theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
    105         FissionProbability = theFissionProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
    106     }
    107     else {
    108         MaximalKineticEnergy = -1000.0*MeV;
    109         LevelDensityParameter = 0.0;
    110         FissionProbability = 0.0;
    111     }
    112 
    113     return;
    114 }
    115 
    116 
     71  G4int anA = fragment.GetA_asInt();
     72  G4int aZ  = fragment.GetZ_asInt();
     73  G4double ExEnergy = fragment.GetExcitationEnergy() -
     74    G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(anA,aZ);
     75 
     76
     77  // Saddle point excitation energy ---> A = 65
     78  // Fission is excluded for A < 65
     79  if (anA >= 65 && ExEnergy > 0.0) {
     80    FissionBarrier = theFissionBarrierPtr->FissionBarrier(anA,aZ,ExEnergy);
     81    MaximalKineticEnergy = ExEnergy - FissionBarrier;
     82    LevelDensityParameter =
     83      theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
     84    FissionProbability =
     85      theFissionProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
     86    }
     87  else {
     88    MaximalKineticEnergy = -1000.0*MeV;
     89    LevelDensityParameter = 0.0;
     90    FissionProbability = 0.0;
     91  }
     92}
    11793
    11894G4FragmentVector * G4CompetitiveFission::BreakUp(const G4Fragment & theNucleus)
    11995{
    120     // Nucleus data
    121     // Atomic number of nucleus
    122     G4int A = static_cast<G4int>(theNucleus.GetA());
    123     // Charge of nucleus
    124     G4int Z = static_cast<G4int>(theNucleus.GetZ());
    125     //   Excitation energy (in MeV)
    126     G4double U = theNucleus.GetExcitationEnergy() -
    127         G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
    128     // Check that U > 0
    129     if (U <= 0.0) {
    130         G4FragmentVector * theResult = new  G4FragmentVector;
    131         theResult->push_back(new G4Fragment(theNucleus));
    132         return theResult;
    133     }
    134 
    135 
    136     // Atomic Mass of Nucleus (in MeV)
    137     G4double M = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
    138 
    139     // Nucleus Momentum
    140     G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum();
    141 
    142     // Calculate fission parameters
    143     G4FissionParameters theParameters(A,Z,U,FissionBarrier);
    144  
    145     // First fragment
    146     G4int A1 = 0;
    147     G4int Z1 = 0;
    148     G4double M1 = 0.0;
    149 
    150     // Second fragment
    151     G4int A2 = 0;
    152     G4int Z2 = 0;
    153     G4double M2 = 0.0;
    154 
    155     G4double FragmentsExcitationEnergy = 0.0;
    156     G4double FragmentsKineticEnergy = 0.0;
    157 
    158     //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation
    159     G4double FissionPairingEnergy=G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
    160 
    161     G4int Trials = 0;
    162     do {
    163 
    164         // First fragment
    165         A1 = FissionAtomicNumber(A,theParameters);
    166         Z1 = FissionCharge(A,Z,A1);
    167         M1 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z1,A1);
    168 
    169 
    170         // Second Fragment
    171         A2 = A - A1;
    172         Z2 = Z - Z1;
    173         if (A2 < 1 || Z2 < 0)
    174             throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakUp: Can't define second fragment! ");
    175 
    176         M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2);
    177 
    178         // Check that fragment masses are less or equal than total energy
    179         if (M1 + M2 > theNucleusMomentum.e())
    180             throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy");
    181 
    182         // Maximal Kinetic Energy (available energy for fragments)
    183         G4double Tmax = M + U - M1 - M2;
    184 
    185         FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
    186                                                        A1, Z1,
    187                                                        A2, Z2,
    188                                                        U , Tmax,
    189                                                        theParameters);
     96  // Nucleus data
     97  // Atomic number of nucleus
     98  G4int A = theNucleus.GetA_asInt();
     99  // Charge of nucleus
     100  G4int Z = theNucleus.GetZ_asInt();
     101  //   Excitation energy (in MeV)
     102  G4double U = theNucleus.GetExcitationEnergy() -
     103    G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
     104  // Check that U > 0
     105  if (U <= 0.0) {
     106    G4FragmentVector * theResult = new  G4FragmentVector;
     107    theResult->push_back(new G4Fragment(theNucleus));
     108    return theResult;
     109  }
     110
     111  // Atomic Mass of Nucleus (in MeV)
     112  G4double M = theNucleus.GetGroundStateMass();
     113
     114  // Nucleus Momentum
     115  G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum();
     116
     117  // Calculate fission parameters
     118  G4FissionParameters theParameters(A,Z,U,FissionBarrier);
     119 
     120  // First fragment
     121  G4int A1 = 0;
     122  G4int Z1 = 0;
     123  G4double M1 = 0.0;
     124
     125  // Second fragment
     126  G4int A2 = 0;
     127  G4int Z2 = 0;
     128  G4double M2 = 0.0;
     129
     130  G4double FragmentsExcitationEnergy = 0.0;
     131  G4double FragmentsKineticEnergy = 0.0;
     132
     133  //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation
     134  G4double FissionPairingEnergy=
     135    G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
     136
     137  G4int Trials = 0;
     138  do {
     139
     140    // First fragment
     141    A1 = FissionAtomicNumber(A,theParameters);
     142    Z1 = FissionCharge(A,Z,A1);
     143    M1 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z1,A1);
     144
     145    // Second Fragment
     146    A2 = A - A1;
     147    Z2 = Z - Z1;
     148    if (A2 < 1 || Z2 < 0) {
     149      throw G4HadronicException(__FILE__, __LINE__,
     150        "G4CompetitiveFission::BreakUp: Can't define second fragment! ");
     151    }
     152    M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2);
     153
     154    // Check that fragment masses are less or equal than total energy
     155    if (M1 + M2 > theNucleusMomentum.e()) {
     156      throw G4HadronicException(__FILE__, __LINE__,
     157        "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy");
     158    }
     159    // Maximal Kinetic Energy (available energy for fragments)
     160    G4double Tmax = M + U - M1 - M2;
     161
     162    FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
     163                                                   A1, Z1,
     164                                                   A2, Z2,
     165                                                   U , Tmax,
     166                                                   theParameters);
    190167   
    191         // Excitation Energy
    192         //      FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
    193         // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
    194         // fragments carry the fission pairing energy in form of
    195         //excitation energy
    196 
    197         FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy+FissionPairingEnergy;
    198 
    199     } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
     168    // Excitation Energy
     169    //  FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
     170    // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
     171    // fragments carry the fission pairing energy in form of
     172    //excitation energy
     173
     174    FragmentsExcitationEnergy =
     175      Tmax - FragmentsKineticEnergy+FissionPairingEnergy;
     176
     177  } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
    200178   
    201  
    202 
    203     if (FragmentsExcitationEnergy <= 0.0)
    204         throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
    205  
    206 
    207     // while (FragmentsExcitationEnergy < 0 && Trials < 100);
    208  
    209     // Fragment 1
    210     G4double U1 = FragmentsExcitationEnergy * (static_cast<G4double>(A1)/static_cast<G4double>(A));
     179  if (FragmentsExcitationEnergy <= 0.0) {
     180    throw G4HadronicException(__FILE__, __LINE__,
     181      "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
     182  }
     183
     184  // while (FragmentsExcitationEnergy < 0 && Trials < 100);
     185 
     186  // Fragment 1
     187  G4double U1 = FragmentsExcitationEnergy * A1/static_cast<G4double>(A);
    211188    // Fragment 2
    212     G4double U2 = FragmentsExcitationEnergy * (static_cast<G4double>(A2)/static_cast<G4double>(A));
    213 
    214 
    215     //JMQ  04/03/09 Full relativistic calculation is performed
    216     //
    217     G4double Fragment1KineticEnergy=(FragmentsKineticEnergy*(FragmentsKineticEnergy+2*(M2+U2)))/(2*(M1+U1+M2+U2+FragmentsKineticEnergy));
    218     G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt(Fragment1KineticEnergy*(Fragment1KineticEnergy+2*(M1+U1)))));
    219     G4ParticleMomentum Momentum2(-Momentum1);
    220     G4LorentzVector FourMomentum1(Momentum1,std::sqrt(Momentum1.mag2()+(M1+U1)*(M1+U1)));
    221     G4LorentzVector FourMomentum2(Momentum2,std::sqrt(Momentum2.mag2()+(M2+U2)*(M2+U2)));
    222 
    223     //JMQ 04/03/09 now we do Lorentz boosts (instead of Galileo boosts)
    224     FourMomentum1.boost(theNucleusMomentum.boostVector());
    225     FourMomentum2.boost(theNucleusMomentum.boostVector());
     189  G4double U2 = FragmentsExcitationEnergy * A2/static_cast<G4double>(A);
     190
     191  //JMQ  04/03/09 Full relativistic calculation is performed
     192  //
     193  G4double Fragment1KineticEnergy=
     194    (FragmentsKineticEnergy*(FragmentsKineticEnergy+2*(M2+U2)))
     195    /(2*(M1+U1+M2+U2+FragmentsKineticEnergy));
     196  G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt(Fragment1KineticEnergy*(Fragment1KineticEnergy+2*(M1+U1)))));
     197  G4ParticleMomentum Momentum2(-Momentum1);
     198  G4LorentzVector FourMomentum1(Momentum1,std::sqrt(Momentum1.mag2()+(M1+U1)*(M1+U1)));
     199  G4LorentzVector FourMomentum2(Momentum2,std::sqrt(Momentum2.mag2()+(M2+U2)*(M2+U2)));
     200
     201  //JMQ 04/03/09 now we do Lorentz boosts (instead of Galileo boosts)
     202  FourMomentum1.boost(theNucleusMomentum.boostVector());
     203  FourMomentum2.boost(theNucleusMomentum.boostVector());
    226204   
    227     //////////JMQ 04/03: Old version calculation is commented
    228     // There was vioation of energy momentum conservation
    229 
    230     //    G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) /
    231     //                          ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy);
    232 
    233     //G4ParticleMomentum momentum1 = IsotropicVector( Pmax );
    234     //  G4ParticleMomentum momentum2( -momentum1 );
    235 
    236     // Perform a Galileo boost for fragments
    237     //    momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));
    238     //    momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));
    239 
    240 
    241     // Create 4-momentum for first fragment
    242     // Warning!! Energy conservation is broken
    243     //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
    244     //    G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1)));
    245 
    246     // Create 4-momentum for second fragment
    247     // Warning!! Energy conservation is broken
    248     //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
    249     //    G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2)));
    250 
    251     //////////
    252 
    253     // Create Fragments
    254     G4Fragment * Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
    255     if (!Fragment1) throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakItUp: Can't create Fragment1! ");
    256     G4Fragment * Fragment2 = new G4Fragment( A2, Z2, FourMomentum2);
    257     if (!Fragment2) throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakItUp: Can't create Fragment2! ");
    258 
    259 #ifdef PRECOMPOUND_TEST
    260     Fragment1->SetCreatorModel(G4String("G4CompetitiveFission"));
    261     Fragment2->SetCreatorModel(G4String("G4CompetitiveFission"));
     205  //////////JMQ 04/03: Old version calculation is commented
     206  // There was vioation of energy momentum conservation
     207
     208  //    G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) /
     209  //                            ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy);
     210
     211  //G4ParticleMomentum momentum1 = IsotropicVector( Pmax );
     212  //  G4ParticleMomentum momentum2( -momentum1 );
     213
     214  // Perform a Galileo boost for fragments
     215  //    momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));
     216  //    momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));
     217
     218
     219  // Create 4-momentum for first fragment
     220  // Warning!! Energy conservation is broken
     221  //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
     222  //    G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1)));
     223
     224  // Create 4-momentum for second fragment
     225  // Warning!! Energy conservation is broken
     226  //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
     227  //    G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2)));
     228
     229  //////////
     230
     231  // Create Fragments
     232  G4Fragment * Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
     233  G4Fragment * Fragment2 = new G4Fragment( A2, Z2, FourMomentum2);
     234
     235  // Create Fragment Vector
     236  G4FragmentVector * theResult = new G4FragmentVector;
     237
     238  theResult->push_back(Fragment1);
     239  theResult->push_back(Fragment2);
     240
     241#ifdef debug
     242  CheckConservation(theNucleus,theResult);
    262243#endif
    263     // Create Fragment Vector
    264     G4FragmentVector * theResult = new G4FragmentVector;
    265 
    266     theResult->push_back(Fragment1);
    267     theResult->push_back(Fragment2);
    268 
    269 #ifdef debug
    270     CheckConservation(theNucleus,theResult);
    271 #endif
    272 
    273     return theResult;
    274 }
    275 
    276 
    277 
    278 G4int G4CompetitiveFission::FissionAtomicNumber(const G4int A, const G4FissionParameters & theParam)
    279     // Calculates the atomic number of a fission product
    280 {
    281 
    282     // For Simplicity reading code
    283     const G4double A1 = theParam.GetA1();
    284     const G4double A2 = theParam.GetA2();
    285     const G4double As = theParam.GetAs();
    286 //    const G4double Sigma1 = theParam.GetSigma1();
    287     const G4double Sigma2 = theParam.GetSigma2();
    288     const G4double SigmaS = theParam.GetSigmaS();
    289     const G4double w = theParam.GetW();
    290 
    291  
    292 //    G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) +
    293 //      std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
    294 
    295 //    G4double FsymA1A2 = std::exp(-((As-(A1+A2))*(As-(A1+A2)))/(2.0*SigmaS*SigmaS));
    296 
    297 
    298     G4double C2A = A2 + 3.72*Sigma2;
    299     G4double C2S = As + 3.72*SigmaS;
    300  
    301     G4double C2 = 0.0;
    302     if (w > 1000.0 ) C2 = C2S;
    303     else if (w < 0.001) C2 = C2A;
    304     else C2 =  std::max(C2A,C2S);
    305 
    306     G4double C1 = A-C2;
    307     if (C1 < 30.0) {
    308         C2 = A-30.0;
    309         C1 = 30.0;
    310     }
    311 
    312     G4double Am1 = (As + A1)/2.0;
    313     G4double Am2 = (A1 + A2)/2.0;
    314 
    315     // Get Mass distributions as sum of symmetric and asymmetric Gasussians
    316     G4double Mass1 = MassDistribution(As,A,theParam);
    317     G4double Mass2 = MassDistribution(Am1,A,theParam);
    318     G4double Mass3 = MassDistribution(A1,A,theParam);
    319     G4double Mass4 = MassDistribution(Am2,A,theParam);
    320     G4double Mass5 = MassDistribution(A2,A,theParam);
    321     // get maximal value among Mass1,...,Mass5
    322     G4double MassMax = Mass1;
    323     if (Mass2 > MassMax) MassMax = Mass2;
    324     if (Mass3 > MassMax) MassMax = Mass3;
    325     if (Mass4 > MassMax) MassMax = Mass4;
    326     if (Mass5 > MassMax) MassMax = Mass5;
    327 
    328     // Sample a fragment mass number, which lies between C1 and C2
    329     G4double m;
    330     G4double Pm;
    331     do {
    332         m = C1+G4UniformRand()*(C2-C1);
    333         Pm = MassDistribution(m,A,theParam);
    334     } while (G4UniformRand() > Pm/MassMax);
    335 
    336     return static_cast<G4int>(m+0.5);
    337 }
    338 
    339 
    340 
    341 
    342 
    343 
    344 G4double G4CompetitiveFission::MassDistribution(const G4double x, const G4double A,
    345                                                 const G4FissionParameters & theParam)
    346     // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
    347     // which consist of symmetric and asymmetric sum of gaussians components.
    348 {
    349     G4double Xsym = std::exp(-0.5*(x-theParam.GetAs())*(x-theParam.GetAs())/
    350                         (theParam.GetSigmaS()*theParam.GetSigmaS()));
    351 
    352     G4double Xasym = std::exp(-0.5*(x-theParam.GetA2())*(x-theParam.GetA2())/
    353                          (theParam.GetSigma2()*theParam.GetSigma2())) +
    354         std::exp(-0.5*(x-(A-theParam.GetA2()))*(x-(A-theParam.GetA2()))/
    355             (theParam.GetSigma2()*theParam.GetSigma2())) +
    356         0.5*std::exp(-0.5*(x-theParam.GetA1())*(x-theParam.GetA1())/
    357                 (theParam.GetSigma1()*theParam.GetSigma1())) +
    358         0.5*std::exp(-0.5*(x-(A-theParam.GetA1()))*(x-(A-theParam.GetA1()))/
    359                 (theParam.GetSigma1()*theParam.GetSigma1()));
    360 
    361     if (theParam.GetW() > 1000) return Xsym;
    362     else if (theParam.GetW() < 0.001) return Xasym;
    363     else return theParam.GetW()*Xsym+Xasym;
    364 }
    365 
    366 
    367 
    368 
    369 G4int G4CompetitiveFission::FissionCharge(const G4double A,
    370                                           const G4double Z,
    371                                           const G4double Af)
    372     // Calculates the charge of a fission product for a given atomic number Af
    373 {
    374     const G4double sigma = 0.6;
    375     G4double DeltaZ = 0.0;
    376     if (Af >= 134.0) DeltaZ = -0.45;                    //                      134 <= Af
    377     else if (Af <= (A-134.0)) DeltaZ = 0.45;             // Af <= (A-134)
    378     else DeltaZ = -0.45*(Af-(A/2.0))/(134.0-(A/2.0));   //       (A-134) < Af < 134
    379 
    380     G4double Zmean = (Af/A)*Z + DeltaZ;
     244
     245  return theResult;
     246}
     247
     248G4int
     249G4CompetitiveFission::FissionAtomicNumber(G4int A,
     250                                          const G4FissionParameters & theParam)
     251  // Calculates the atomic number of a fission product
     252{
     253
     254  // For Simplicity reading code
     255  const G4double A1 = theParam.GetA1();
     256  const G4double A2 = theParam.GetA2();
     257  const G4double As = theParam.GetAs();
     258  //    const G4double Sigma1 = theParam.GetSigma1();
     259  const G4double Sigma2 = theParam.GetSigma2();
     260  const G4double SigmaS = theParam.GetSigmaS();
     261  const G4double w = theParam.GetW();
     262 
     263  //    G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) +
     264  //    std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
     265
     266  //    G4double FsymA1A2 = std::exp(-((As-(A1+A2))*(As-(A1+A2)))/(2.0*SigmaS*SigmaS));
     267
     268  G4double C2A = A2 + 3.72*Sigma2;
     269  G4double C2S = As + 3.72*SigmaS;
     270 
     271  G4double C2 = 0.0;
     272  if (w > 1000.0 ) C2 = C2S;
     273  else if (w < 0.001) C2 = C2A;
     274  else C2 =  std::max(C2A,C2S);
     275
     276  G4double C1 = A-C2;
     277  if (C1 < 30.0) {
     278    C2 = A-30.0;
     279    C1 = 30.0;
     280  }
     281
     282  G4double Am1 = (As + A1)/2.0;
     283  G4double Am2 = (A1 + A2)/2.0;
     284
     285  // Get Mass distributions as sum of symmetric and asymmetric Gasussians
     286  G4double Mass1 = MassDistribution(As,A,theParam);
     287  G4double Mass2 = MassDistribution(Am1,A,theParam);
     288  G4double Mass3 = MassDistribution(A1,A,theParam);
     289  G4double Mass4 = MassDistribution(Am2,A,theParam);
     290  G4double Mass5 = MassDistribution(A2,A,theParam);
     291  // get maximal value among Mass1,...,Mass5
     292  G4double MassMax = Mass1;
     293  if (Mass2 > MassMax) MassMax = Mass2;
     294  if (Mass3 > MassMax) MassMax = Mass3;
     295  if (Mass4 > MassMax) MassMax = Mass4;
     296  if (Mass5 > MassMax) MassMax = Mass5;
     297
     298  // Sample a fragment mass number, which lies between C1 and C2
     299  G4double m;
     300  G4double Pm;
     301  do {
     302    m = C1+G4UniformRand()*(C2-C1);
     303    Pm = MassDistribution(m,A,theParam);
     304  } while (MassMax*G4UniformRand() > Pm);
     305
     306  return static_cast<G4int>(m+0.5);
     307}
     308
     309G4double
     310G4CompetitiveFission::MassDistribution(G4double x, G4double A,
     311                                       const G4FissionParameters & theParam)
     312  // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
     313  // which consist of symmetric and asymmetric sum of gaussians components.
     314{
     315  G4double Xsym = std::exp(-0.5*(x-theParam.GetAs())*(x-theParam.GetAs())/
     316                           (theParam.GetSigmaS()*theParam.GetSigmaS()));
     317
     318  G4double Xasym = std::exp(-0.5*(x-theParam.GetA2())*(x-theParam.GetA2())/
     319                            (theParam.GetSigma2()*theParam.GetSigma2())) +
     320    std::exp(-0.5*(x-(A-theParam.GetA2()))*(x-(A-theParam.GetA2()))/
     321             (theParam.GetSigma2()*theParam.GetSigma2())) +
     322    0.5*std::exp(-0.5*(x-theParam.GetA1())*(x-theParam.GetA1())/
     323                 (theParam.GetSigma1()*theParam.GetSigma1())) +
     324    0.5*std::exp(-0.5*(x-(A-theParam.GetA1()))*(x-(A-theParam.GetA1()))/
     325                 (theParam.GetSigma1()*theParam.GetSigma1()));
     326
     327  if (theParam.GetW() > 1000) return Xsym;
     328  else if (theParam.GetW() < 0.001) return Xasym;
     329  else return theParam.GetW()*Xsym+Xasym;
     330}
     331
     332G4int G4CompetitiveFission::FissionCharge(G4double A, G4double Z,
     333                                          G4double Af)
     334  // Calculates the charge of a fission product for a given atomic number Af
     335{
     336  const G4double sigma = 0.6;
     337  G4double DeltaZ = 0.0;
     338  if (Af >= 134.0) DeltaZ = -0.45;                    //                      134 <= Af
     339  else if (Af <= (A-134.0)) DeltaZ = 0.45;             // Af <= (A-134)
     340  else DeltaZ = -0.45*(Af-(A/2.0))/(134.0-(A/2.0));   //       (A-134) < Af < 134
     341
     342  G4double Zmean = (Af/A)*Z + DeltaZ;
    381343 
    382     G4double theZ;
    383     do {
    384         theZ = G4RandGauss::shoot(Zmean,sigma);
    385     } while (theZ  < 1.0 || theZ > (Z-1.0) || theZ > Af);
    386     //  return static_cast<G4int>(theZ+0.5);
    387     return static_cast<G4int>(theZ+0.5);
    388 }
    389 
    390 
    391 
    392 
    393 G4double G4CompetitiveFission::FissionKineticEnergy(const G4double A, const G4double Z,
    394                                                     const G4double Af1, const G4double /*Zf1*/,
    395                                                     const G4double Af2, const G4double /*Zf2*/,
    396                                                     const G4double /*U*/, const G4double Tmax,
    397                                                     const G4FissionParameters & theParam)
    398     // Gives the kinetic energy of fission products
    399 {
    400     // Find maximal value of A for fragments
    401     G4double AfMax = std::max(Af1,Af2);
    402     if (AfMax < (A/2.0)) AfMax = A - AfMax;
    403 
    404     // Weights for symmetric and asymmetric components
    405     G4double Pas;
    406     if (theParam.GetW() > 1000) Pas = 0.0;
    407     else {
    408         G4double P1 = 0.5*std::exp(-0.5*(AfMax-theParam.GetA1())*(AfMax-theParam.GetA1())/
    409                               (theParam.GetSigma1()*theParam.GetSigma1()));
    410 
    411         G4double P2 = std::exp(-0.5*(AfMax-theParam.GetA2())*(AfMax-theParam.GetA2())/
    412                           (theParam.GetSigma2()*theParam.GetSigma2()));
    413 
    414         Pas = P1+P2;
    415     }
    416 
    417 
    418     G4double Ps;
    419     if (theParam.GetW() < 0.001) Ps = 0.0;
    420     else
    421         Ps = theParam.GetW()*std::exp(-0.5*(AfMax-theParam.GetAs())*(AfMax-theParam.GetAs())/
    422                                  (theParam.GetSigmaS()*theParam.GetSigmaS()));
    423  
    424 
    425 
    426     G4double Psy = Ps/(Pas+Ps);
    427 
    428 
    429     // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
    430     G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
    431     G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
    432     G4double Xas = PPas / (PPas+PPsy);
    433     G4double Xsy = PPsy / (PPas+PPsy);
    434 
    435 
    436     // Average kinetic energy for symmetric and asymmetric components
    437     G4double Eaverage = 0.1071*MeV*(Z*Z)/std::pow(A,1.0/3.0) + 22.2*MeV;
    438 
    439 
    440     // Compute maximal average kinetic energy of fragments and Energy Dispersion (sqrt)
    441     G4double TaverageAfMax;
    442     G4double ESigma;
    443     // Select randomly fission mode (symmetric or asymmetric)
    444     if (G4UniformRand() > Psy) { // Asymmetric Mode
    445         G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1();
    446         G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1();
    447         G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2();
    448         G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2();
    449         // scale factor
    450         G4double ScaleFactor = 0.5*theParam.GetSigma1()*(AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
    451             theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22));
    452         // Compute average kinetic energy for fragment with AfMax
    453         TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * AsymmetricRatio(A,AfMax);
    454         ESigma = 10.0*MeV; // MeV
    455 
    456     } else { // Symmetric Mode
    457         G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
    458         // scale factor
    459         G4double ScaleFactor = theParam.GetW()*theParam.GetSigmaS()*SymmetricRatio(A,As0);
    460         // Compute average kinetic energy for fragment with AfMax
    461         TaverageAfMax = (Eaverage - 12.5*MeV*Xas) * (PPsy/ScaleFactor) * SymmetricRatio(A,AfMax);
    462         ESigma = 8.0*MeV;
    463     }
    464 
    465 
    466     // Select randomly, in accordance with Gaussian distribution, fragment kinetic energy
    467     G4double KineticEnergy;
    468     G4int i = 0;
    469     do {
    470         KineticEnergy = G4RandGauss::shoot(TaverageAfMax,ESigma);
    471         if (i++ > 100) return Eaverage;
    472     } while (KineticEnergy < Eaverage-3.72*ESigma ||
    473              KineticEnergy > Eaverage+3.72*ESigma ||
    474              KineticEnergy > Tmax);
    475 
    476     return KineticEnergy;
    477 }
    478 
    479 
    480 
    481 
    482 
    483 G4double G4CompetitiveFission::AsymmetricRatio(const G4double A,const G4double A11)
    484 {
    485     const G4double B1 = 23.5;
    486     const G4double A00 = 134.0;
    487     return Ratio(A,A11,B1,A00);
    488 }
    489 
    490 G4double G4CompetitiveFission::SymmetricRatio(const G4double A,const G4double A11)
    491 {
    492     const G4double B1 = 5.32;
    493     const G4double A00 = A/2.0;
    494     return Ratio(A,A11,B1,A00);
    495 }
    496 
    497 G4double G4CompetitiveFission::Ratio(const G4double A,const G4double A11,
    498                                      const G4double B1,const G4double A00)
    499 {
    500     if (A == 0) throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::Ratio: A == 0!");
    501     if (A11 >= A/2.0 && A11 <= (A00+10.0)) return 1.0-B1*((A11-A00)/A)*((A11-A00)/A);
    502     else return 1.0-B1*(10.0/A)*(10.0/A)-2.0*(10.0/A)*B1*((A11-A00-10.0)/A);
    503 }
    504 
    505 
    506 
    507 
     344  G4double theZ;
     345  do {
     346    theZ = G4RandGauss::shoot(Zmean,sigma);
     347  } while (theZ  < 1.0 || theZ > (Z-1.0) || theZ > Af);
     348  //  return static_cast<G4int>(theZ+0.5);
     349  return static_cast<G4int>(theZ+0.5);
     350}
     351
     352G4double
     353G4CompetitiveFission::FissionKineticEnergy(G4int A, G4int Z,
     354                                           G4double Af1, G4double /*Zf1*/,
     355                                           G4double Af2, G4double /*Zf2*/,
     356                                           G4double /*U*/, G4double Tmax,
     357                                           const G4FissionParameters & theParam)
     358  // Gives the kinetic energy of fission products
     359{
     360  // Find maximal value of A for fragments
     361  G4double AfMax = std::max(Af1,Af2);
     362  if (AfMax < (A/2.0)) AfMax = A - AfMax;
     363
     364  // Weights for symmetric and asymmetric components
     365  G4double Pas;
     366  if (theParam.GetW() > 1000) Pas = 0.0;
     367  else {
     368    G4double P1 = 0.5*std::exp(-0.5*(AfMax-theParam.GetA1())*(AfMax-theParam.GetA1())/
     369                               (theParam.GetSigma1()*theParam.GetSigma1()));
     370
     371    G4double P2 = std::exp(-0.5*(AfMax-theParam.GetA2())*(AfMax-theParam.GetA2())/
     372                           (theParam.GetSigma2()*theParam.GetSigma2()));
     373
     374    Pas = P1+P2;
     375  }
     376
     377  G4double Ps;
     378  if (theParam.GetW() < 0.001) Ps = 0.0;
     379  else {
     380    Ps = theParam.GetW()*std::exp(-0.5*(AfMax-theParam.GetAs())*(AfMax-theParam.GetAs())/
     381                                  (theParam.GetSigmaS()*theParam.GetSigmaS()));
     382  }
     383  G4double Psy = Ps/(Pas+Ps);
     384
     385  // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
     386  G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
     387  G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
     388  G4double Xas = PPas / (PPas+PPsy);
     389  G4double Xsy = PPsy / (PPas+PPsy);
     390
     391  // Average kinetic energy for symmetric and asymmetric components
     392  G4double Eaverage = 0.1071*MeV*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2*MeV;
     393
     394
     395  // Compute maximal average kinetic energy of fragments and Energy Dispersion (sqrt)
     396  G4double TaverageAfMax;
     397  G4double ESigma;
     398  // Select randomly fission mode (symmetric or asymmetric)
     399  if (G4UniformRand() > Psy) { // Asymmetric Mode
     400    G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1();
     401    G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1();
     402    G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2();
     403    G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2();
     404    // scale factor
     405    G4double ScaleFactor = 0.5*theParam.GetSigma1()*(AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
     406      theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22));
     407    // Compute average kinetic energy for fragment with AfMax
     408    TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * AsymmetricRatio(A,AfMax);
     409    ESigma = 10.0*MeV; // MeV
     410
     411  } else { // Symmetric Mode
     412    G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
     413    // scale factor
     414    G4double ScaleFactor = theParam.GetW()*theParam.GetSigmaS()*SymmetricRatio(A,As0);
     415    // Compute average kinetic energy for fragment with AfMax
     416    TaverageAfMax = (Eaverage - 12.5*MeV*Xas) * (PPsy/ScaleFactor) * SymmetricRatio(A,AfMax);
     417    ESigma = 8.0*MeV;
     418  }
     419
     420
     421  // Select randomly, in accordance with Gaussian distribution, fragment kinetic energy
     422  G4double KineticEnergy;
     423  G4int i = 0;
     424  do {
     425    KineticEnergy = G4RandGauss::shoot(TaverageAfMax,ESigma);
     426    if (i++ > 100) return Eaverage;
     427  } while (KineticEnergy < Eaverage-3.72*ESigma ||
     428           KineticEnergy > Eaverage+3.72*ESigma ||
     429           KineticEnergy > Tmax);
     430 
     431  return KineticEnergy;
     432}
     433
     434G4double G4CompetitiveFission::AsymmetricRatio(G4int A, G4double A11)
     435{
     436  const G4double B1 = 23.5;
     437  const G4double A00 = 134.0;
     438  return Ratio(G4double(A),A11,B1,A00);
     439}
     440
     441G4double G4CompetitiveFission::SymmetricRatio(G4int A, G4double A11)
     442{
     443  const G4double B1 = 5.32;
     444  const G4double A00 = A/2.0;
     445  return Ratio(G4double(A),A11,B1,A00);
     446}
     447
     448G4double G4CompetitiveFission::Ratio(G4double A, G4double A11,
     449                                     G4double B1, G4double A00)
     450{
     451  if (A == 0.0) {
     452    throw G4HadronicException(__FILE__, __LINE__,
     453                              "G4CompetitiveFission::Ratio: A == 0!");
     454  }
     455  if (A11 >= A/2.0 && A11 <= (A00+10.0)) {
     456    return 1.0-B1*((A11-A00)/A)*((A11-A00)/A);
     457  } else {
     458    return 1.0-B1*(10.0/A)*(10.0/A)-2.0*(10.0/A)*B1*((A11-A00-10.0)/A);
     459  }
     460}
    508461
    509462G4ThreeVector G4CompetitiveFission::IsotropicVector(const G4double Magnitude)
    510     // Samples a isotropic random vectorwith a magnitud given by Magnitude.
    511     // By default Magnitude = 1.0
    512 {
    513     G4double CosTheta = 1.0 - 2.0*G4UniformRand();
    514     G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
    515     G4double Phi = twopi*G4UniformRand();
    516     G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
    517                          Magnitude*std::sin(Phi)*SinTheta,
    518                          Magnitude*CosTheta);
    519     return Vector;
    520 }
    521 
     463  // Samples a isotropic random vectorwith a magnitud given by Magnitude.
     464  // By default Magnitude = 1.0
     465{
     466  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
     467  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
     468  G4double Phi = twopi*G4UniformRand();
     469  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
     470                       Magnitude*std::sin(Phi)*SinTheta,
     471                       Magnitude*CosTheta);
     472  return Vector;
     473}
    522474
    523475#ifdef debug

  • trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionBarrier.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4FissionBarrier.cc,v 1.7 2009/12/16 16:50:07 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4FissionBarrier.cc,v 1.8 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Oct 1998)
    32 
     32//
     33// 17.11.2010 V.Ivanchenko cleanup and add usage of G4Pow
    3334
    3435#include "G4FissionBarrier.hh"
     36#include "G4Pow.hh"
    3537
    36 G4FissionBarrier::G4FissionBarrier(const G4FissionBarrier & ) : G4VFissionBarrier()
     38G4FissionBarrier::G4FissionBarrier()
     39{}
     40
     41G4FissionBarrier::~G4FissionBarrier()
     42{}
     43
     44G4double
     45G4FissionBarrier::FissionBarrier(G4int A, G4int Z, G4double U)
     46  // Compute fission barrier according with Barashenkov's
     47  // prescription for A >= 65
    3748{
    38     throw G4HadronicException(__FILE__, __LINE__, "G4FissionBarrier::copy_constructor meant to not be accessable.");
     49  if (A >= 65) {
     50    return BarashenkovFissionBarrier(A,Z)/(1.0 + std::sqrt(U/(2.0*A)));
     51  } else { return 100.0*GeV; }
    3952}
    4053
     54G4double
     55G4FissionBarrier::BarashenkovFissionBarrier(G4int A, G4int Z)
     56  // Calculates Fission Barrier heights
     57{
     58  // Liquid drop model parameters for
     59  // surface energy of a spherical nucleus
     60  const G4double aSurf = 17.9439*MeV;
     61  // and coulomb energy
     62  const G4double aCoul = 0.7053*MeV;
     63  const G4double k = 1.7826;
     64  G4int N = A - Z;
    4165
    42 const G4FissionBarrier & G4FissionBarrier::operator=(const G4FissionBarrier & )
    43 {
    44     throw G4HadronicException(__FILE__, __LINE__, "G4FissionBarrier::operator= meant to not be accessable.");
    45     return *this;
     66  // fissibility parameter
     67  G4double x = (aCoul/(2.0*aSurf))*(Z*Z)/static_cast<G4double>(A);
     68  x /= (1.0 - k*(N-Z)*(N-Z)/static_cast<G4double>(A*A));
     69 
     70  // Liquid drop model part of Fission Barrier
     71  G4double BF0 = aSurf*G4Pow::GetInstance()->Z23(A);
     72  if (x <= 2./3.) { BF0 *= 0.38*(3./4.-x); }
     73  else { BF0 *= 0.83*(1. - x)*(1. - x)*(1. - x); }
     74
     75  //
     76  G4double D = 1.248*MeV;
     77  D *= (N - 2*(N/2) + Z - 2*(Z/2));
     78
     79  return BF0 + D - SellPlusPairingCorrection(Z,N);
    4680}
    4781
    48 G4bool G4FissionBarrier::operator==(const G4FissionBarrier & ) const
    49 {
    50     return false;
    51 }
    52 
    53 G4bool G4FissionBarrier::operator!=(const G4FissionBarrier & ) const
    54 {
    55     return true;
    56 }
    57 
    58 
    59 
    60 G4double G4FissionBarrier::FissionBarrier(const G4int A, const G4int Z, const G4double U)
    61     // Compute fission barrier according with Barashenkov's prescription for A >= 65
    62 {
    63     if (A >= 65) return BarashenkovFissionBarrier(A,Z)/(1.0 + std::sqrt(U/(2.0*static_cast<G4double>(A))));
    64     else return 100.0*GeV;
    65 }
    66 
    67 
    68 G4double G4FissionBarrier::BarashenkovFissionBarrier(const G4int A, const G4int Z)
    69     // Calculates Fission Barrier heights
    70 {
    71     // Liquid drop model parameters for
    72     // surface energy of a spherical nucleus
    73     const G4double aSurf = 17.9439*MeV;
    74     // and coulomb energy
    75     const G4double aCoul = 0.7053*MeV;
    76     const G4int N = A - Z;
    77     const G4double k = 1.7826;
    78 
    79     // fissibility parameter
    80     G4double x = (aCoul/(2.0*aSurf))*(static_cast<G4double>(Z)*static_cast<G4double>(Z))/static_cast<G4double>(A);
    81     x /= (1.0 - k*(static_cast<G4double>(N-Z)/static_cast<G4double>(A))*
    82           (static_cast<G4double>(N-Z)/static_cast<G4double>(A)));
    83  
    84     // Liquid drop model part of Fission Barrier
    85     G4double BF0 = aSurf*std::pow(static_cast<G4double>(A),2./3.);
    86     if (x <= 2./3.) BF0 *= 0.38*(3./4.-x);
    87     else BF0 *= 0.83*(1. - x)*(1. - x)*(1. - x);
    88 
    89 
    90     //
    91     G4double D = 1.248*MeV;
    92     D *= (static_cast<G4double>(N)-2.0*(N/2)) + (static_cast<G4double>(Z)-2.0*(Z/2));
    93 
    94     return BF0 + D - SellPlusPairingCorrection(Z,N);
    95 }
    96 

  • trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionLevelDensityParameter.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4FissionLevelDensityParameter.cc,v 1.7 2009/11/21 18:05:26 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4FissionLevelDensityParameter.cc,v 1.8 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    3939
    4040
    41 G4FissionLevelDensityParameter::
    42 G4FissionLevelDensityParameter(const G4FissionLevelDensityParameter &) : G4VLevelDensityParameter()
    43 {
    44     throw G4HadronicException(__FILE__, __LINE__, "G4FissionLevelDensityParameter::copy_constructor meant to not be accessable");
    45 }
     41G4FissionLevelDensityParameter::G4FissionLevelDensityParameter()
     42{}
    4643
    47 
    48 const G4FissionLevelDensityParameter & G4FissionLevelDensityParameter::
    49 operator=(const G4FissionLevelDensityParameter &)
    50 {
    51     throw G4HadronicException(__FILE__, __LINE__, "G4FissionLevelDensityParameter::operator= meant to not be accessable");
    52     return *this;
    53 }
    54 
    55 
    56 G4bool G4FissionLevelDensityParameter::
    57 operator==(const G4FissionLevelDensityParameter &) const
    58 {
    59     return false;
    60 }
    61 
    62 G4bool G4FissionLevelDensityParameter::
    63 operator!=(const G4FissionLevelDensityParameter &) const
    64 {
    65     return true;
    66 }
     44G4FissionLevelDensityParameter::~G4FissionLevelDensityParameter()
     45{}
    6746
    6847
    6948G4double G4FissionLevelDensityParameter::
    70 LevelDensityParameter(const G4int A,const G4int Z,const G4double U) const
     49LevelDensityParameter(G4int A, G4int Z, G4double U) const
    7150{
    7251  G4double EvapLDP =

  • trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionParameters.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4FissionParameters.cc,v 1.7 2009/11/19 10:30:49 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4FissionParameters.cc,v 1.8 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    4141const G4double G4FissionParameters::A2 = 141.0;
    4242
    43 
    44 G4FissionParameters::G4FissionParameters(const G4int A, const G4int Z, const G4double ExEnergy,
    45                                          const G4double FissionBarrier)
     43G4FissionParameters::G4FissionParameters(G4int A, G4int Z, G4double ExEnergy,
     44                                         G4double FissionBarrier)
    4645{
    47     G4double U = ExEnergy;
     46  G4double U = ExEnergy;
     47 
     48  As = A/2.0;
    4849   
    49     As = A/2.0;
     50  if (A <= 235) Sigma2 = 5.6;  // MeV
     51  else Sigma2 = 5.6 + 0.096*(A-235); // MeV
    5052   
    51     if (A <= 235) Sigma2 = 5.6;  // MeV
    52     else Sigma2 = 5.6 + 0.096*(A-235); // MeV
     53  Sigma1 = 0.5*Sigma2; // MeV
    5354   
    54     Sigma1 = 0.5*Sigma2; // MeV
     55  SigmaS = std::exp(0.00553*U/MeV + 2.1386); // MeV
    5556   
    56     SigmaS = std::exp(0.00553*U/MeV + 2.1386); // MeV
     57  //JMQ 310509
     58  //    if (SigmaS > 20.0) SigmaS = 20.0;
     59  //   SigmaS*=1.3;
     60  //JMQ 301009: retuning (after CEM transition prob.have been chosen as default)
     61  SigmaS*=0.8;
     62  //
    5763   
    58     //JMQ 310509
    59     //    if (SigmaS > 20.0) SigmaS = 20.0;
    60     //   SigmaS*=1.3;
    61     //JMQ 301009: retuning (after CEM transition prob.have been chosen as default)
    62     SigmaS*=0.8;
    63     //
     64  G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) +
     65    std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
    6466   
    65     G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) +
    66       std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
    67    
    68     G4double FsymA1A2 = std::exp(-((As-(A1+A2)/2.0)*(As-(A1+A2)/2.0))/(2.0*SigmaS*SigmaS));
     67  G4double FsymA1A2 = std::exp(-((As-(A1+A2)/2.0)*(As-(A1+A2)/2.0))
     68                               /(2.0*SigmaS*SigmaS));
    6969   
    7070   
    71     G4double wa = 0.0;
    72     w = 0.0;
    73     if (Z >= 90) {         // Z >= 90
    74       if (U <= 16.25) wa = std::exp(0.5385*U/MeV-9.9564);  // U <= 16.25 MeV
    75       else wa = std::exp(0.09197*U/MeV-2.7003);            // U  > 16.25 MeV
    76     } else if (Z == 89) {  // Z == 89
    77       wa = std::exp(0.09197*U-1.0808);
    78     } else if (Z >= 82) {  //  82 <= Z <= 88
    79       G4double X = FissionBarrier - 7.5*MeV;
    80       if (X < 0.0) X = 0.0;
    81       wa = std::exp(0.09197*(U-X)/MeV-1.0808);
    82     } else {               // Z < 82
    83       w = 1001.0;
    84     }
     71  G4double wa = 0.0;
     72  w = 0.0;
     73  if (Z >= 90) {         // Z >= 90
     74    if (U <= 16.25) wa = std::exp(0.5385*U/MeV-9.9564);  // U <= 16.25 MeV
     75    else wa = std::exp(0.09197*U/MeV-2.7003);            // U  > 16.25 MeV
     76  } else if (Z == 89) {  // Z == 89
     77    wa = std::exp(0.09197*U-1.0808);
     78  } else if (Z >= 82) {  //  82 <= Z <= 88
     79    G4double X = FissionBarrier - 7.5*MeV;
     80    if (X < 0.0) X = 0.0;
     81    wa = std::exp(0.09197*(U-X)/MeV-1.0808);
     82  } else {               // Z < 82
     83    w = 1001.0;
     84  }
    8585   
    86     if (w == 0.0) {
    87       G4double w1 = std::max(1.03*wa - FasymAsym, 0.0001);
    88       G4double w2 = std::max(1.0 - FsymA1A2*wa,   0.0001);
     86  if (w == 0.0) {
     87    G4double w1 = std::max(1.03*wa - FasymAsym, 0.0001);
     88    G4double w2 = std::max(1.0 - FsymA1A2*wa,   0.0001);
     89   
     90    w = w1/w2;
    8991     
    90       w = w1/w2;
    91      
    92       if (82 <= Z && Z < 89 && A < 227)  w *= std::exp(0.3*(227-A));
    93     }
     92    if (82 <= Z && Z < 89 && A < 227)  w *= std::exp(0.3*(227-A));
     93  }
    9494   
    9595}
    9696
     97G4FissionParameters::~G4FissionParameters()
     98{}
    9799
    98 G4FissionParameters::G4FissionParameters(const G4FissionParameters &)
    99 {
    100     throw G4HadronicException(__FILE__, __LINE__, "G4FissionParameters::copy_constructor meant to not be accessable");
    101 }
    102 
    103 
    104 const G4FissionParameters & G4FissionParameters::operator=(const G4FissionParameters &)
    105 {
    106     throw G4HadronicException(__FILE__, __LINE__, "G4FissionParameters::operator= meant to not be accessable");
    107     return *this;
    108 }
    109 
    110 
    111 G4bool G4FissionParameters::operator==(const G4FissionParameters &) const
    112 {
    113     return false;
    114 }
    115 
    116 G4bool G4FissionParameters::operator!=(const G4FissionParameters &) const
    117 {
    118     return true;
    119 }
    120 
    121 

  • trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionProbability.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4FissionProbability.cc,v 1.9 2009/02/15 17:03:25 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4FissionProbability.cc,v 1.10 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    3838#include "G4PairingCorrection.hh"
    3939
     40G4FissionProbability::G4FissionProbability()
     41{}
     42
     43G4FissionProbability::~G4FissionProbability()
     44{}
    4045
    4146
    42 G4FissionProbability::G4FissionProbability(const G4FissionProbability &) : G4VEmissionProbability()
     47G4double
     48G4FissionProbability::EmissionProbability(const G4Fragment & fragment,
     49                                          G4double MaximalKineticEnergy)
     50  // Compute integrated probability of fission channel
    4351{
    44     throw G4HadronicException(__FILE__, __LINE__, "G4FissionProbability::copy_constructor meant to not be accessable");
     52  if (MaximalKineticEnergy <= 0.0) return 0.0;
     53  G4int A = fragment.GetA_asInt();
     54  G4int Z = fragment.GetZ_asInt();
     55  G4double U = fragment.GetExcitationEnergy();
     56 
     57  G4double Ucompound = U -
     58    G4PairingCorrection::GetInstance()->GetPairingCorrection(A,Z);
     59
     60  G4double Ufission = U -
     61    G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
     62 
     63  G4double SystemEntropy =
     64    2.0*std::sqrt(theEvapLDP.LevelDensityParameter(A,Z,Ucompound)*Ucompound);
     65       
     66  G4double afission = theFissLDP.LevelDensityParameter(A,Z,Ufission);
     67
     68  G4double Cf = 2.0*std::sqrt(afission*MaximalKineticEnergy);
     69
     70  //    G4double Q1 = 1.0 + (Cf - 1.0)*std::exp(Cf);
     71  //    G4double Q2 = 4.0*pi*afission*std::exp(SystemEntropy);
     72 
     73  //    G4double probability = Q1/Q2;
     74   
     75  G4double Exp1 = 0.0;
     76  if (SystemEntropy <= 160.0) { Exp1 = std::exp(-SystemEntropy); }
     77  // @@@@@@@@@@@@@@@@@ hpw changed max to min - cannot notify vicente now
     78  G4double Exp2 = std::exp( std::min(700.0,Cf-SystemEntropy) );
     79
     80  // JMQ 14/02/09 BUG fixed in fission probability (missing parenthesis
     81  // at denominator)
     82  //AH fix from Vincente: G4double probability =
     83  //        (Exp1 + (1.0-Cf)*Exp2) / 4.0*pi*afission;
     84  //    G4double probability = (Exp1 + (Cf-1.0)*Exp2) / 4.0*pi*afission;
     85  G4double probability = (Exp1 + (Cf-1.0)*Exp2) / (4.0*pi*afission);
     86
     87  return probability;
    4588}
    4689
    47 
    48 
    49 
    50 const G4FissionProbability & G4FissionProbability::operator=(const G4FissionProbability &)
    51 {
    52     throw G4HadronicException(__FILE__, __LINE__, "G4FissionProbability::operator= meant to not be accessable");
    53     return *this;
    54 }
    55 
    56 
    57 G4bool G4FissionProbability::operator==(const G4FissionProbability &) const
    58 {
    59     return false;
    60 }
    61 
    62 G4bool G4FissionProbability::operator!=(const G4FissionProbability &) const
    63 {
    64     return true;
    65 }
    66 
    67 
    68 G4double G4FissionProbability::EmissionProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy)
    69     // Compute integrated probability of fission channel
    70 {
    71     if (MaximalKineticEnergy <= 0.0) return 0.0;
    72     G4double A = fragment.GetA();
    73     G4double Z = fragment.GetZ();
    74     G4double U = fragment.GetExcitationEnergy();
    75  
    76     G4double Ucompound = U - G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(A),
    77                                                                                       static_cast<G4int>(Z));
    78     G4double Ufission = U - G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(static_cast<G4int>(A),
    79                                                                                             static_cast<G4int>(Z));
    80  
    81     G4double SystemEntropy = 2.0*std::sqrt(theEvapLDP.LevelDensityParameter(static_cast<G4int>(A),
    82                                                                        static_cast<G4int>(Z),
    83                                                                        Ucompound)*Ucompound);
    84        
    85     G4double afission = theFissLDP.LevelDensityParameter(static_cast<G4int>(A),
    86                                                          static_cast<G4int>(Z),
    87                                                          Ufission);
    88 
    89     G4double Cf = 2.0*std::sqrt(afission*MaximalKineticEnergy);
    90 
    91     //    G4double Q1 = 1.0 + (Cf - 1.0)*std::exp(Cf);
    92     //    G4double Q2 = 4.0*pi*afission*std::exp(SystemEntropy);
    93        
    94     //    G4double probability = Q1/Q2;
    95  
    96    
    97     G4double Exp1 = 0.0;
    98     if (SystemEntropy <= 160.0) Exp1 = std::exp(-SystemEntropy);
    99     // @@@@@@@@@@@@@@@@@ hpw changed max to min - cannot notify vicente now
    100     G4double Exp2 = std::exp( std::min(700.0,Cf-SystemEntropy) );
    101 
    102     // JMQ 14/02/09 BUG fixed in fission probability (missing parenthesis at denominator)
    103     //AH fix from Vincente:    G4double probability = (Exp1 + (1.0-Cf)*Exp2) / 4.0*pi*afission;
    104     //    G4double probability = (Exp1 + (Cf-1.0)*Exp2) / 4.0*pi*afission;
    105     G4double probability = (Exp1 + (Cf-1.0)*Exp2) / (4.0*pi*afission);
    106 
    107 
    108    
    109     return probability;
    110 }
    111 

  • trunk/source/processes/hadronic/models/de_excitation/fission/src/G4VFissionBarrier.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4VFissionBarrier.cc,v 1.4 2006/06/29 20:13:43 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4VFissionBarrier.cc,v 1.5 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    3434#include "G4VFissionBarrier.hh"
    3535
    36 G4VFissionBarrier::G4VFissionBarrier(const G4VFissionBarrier & )
    37 {
    38   throw G4HadronicException(__FILE__, __LINE__, "G4VFissionBarrier::copy_constructor meant to not be accessable.");
    39 }
     36G4VFissionBarrier::G4VFissionBarrier() {}
     37G4VFissionBarrier::~G4VFissionBarrier() {}
    4038
    41 
    42 const G4VFissionBarrier & G4VFissionBarrier::operator=(const G4VFissionBarrier & )
    43 {
    44  throw G4HadronicException(__FILE__, __LINE__, "G4VFissionBarrier::operator= meant to not be accessable.");
    45  return *this;
    46 }
    47 
    48 G4bool G4VFissionBarrier::operator==(const G4VFissionBarrier & ) const
    49 {
    50  return false;
    51 }
    52 
    53 G4bool G4VFissionBarrier::operator!=(const G4VFissionBarrier & ) const
    54 {
    55  return true;
    56 }
    57 
Note: See TracChangeset for help on using the changeset viewer.