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
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/fission/include/G4CompetitiveFission.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4CompetitiveFission.hh,v 1.3 2006/06/29 20:13:19 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4CompetitiveFission.hh,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
    … …  
    5858private:
    5959  G4CompetitiveFission(const G4CompetitiveFission &right);
    60 
    6160  const G4CompetitiveFission & operator=(const G4CompetitiveFission &right);
    62 public:
    6361  G4bool operator==(const G4CompetitiveFission &right) const;
    6462  G4bool operator!=(const G4CompetitiveFission &right) const;
    … …  
    121119  G4double LevelDensityParameter;
    122120
    123 
    124 
    125 
    126121  //  --------------------
    127122
    128 
    129123  // Sample AtomicNumber of Fission products
    130   G4int FissionAtomicNumber(const G4int A, const G4FissionParameters & theParam);
    131   G4double MassDistribution(const G4double x, const G4double A, const G4FissionParameters & theParam);
     124  G4int FissionAtomicNumber(G4int A, const G4FissionParameters & theParam);
     125  G4double MassDistribution(G4double x, G4double A, const G4FissionParameters & theParam);
    132126
    133127
    134128  // Sample Charge of fission products
    135   G4int FissionCharge(const G4double A, const G4double Z, const G4double Af);
     129  G4int FissionCharge(G4double A, G4double Z, G4double Af);
    136130
    137131
    138132  // Sample Kinetic energy of fission products
    139   G4double FissionKineticEnergy(const G4double A, const G4double Z,
    140                                 const G4double Af1, const G4double Zf1,
    141                                 const G4double Af2, const G4double Zf2,
    142                                 const G4double U, const G4double Tmax,
     133  G4double FissionKineticEnergy(G4int A, G4int Z,
     134                                G4double Af1, G4double Zf1,
     135                                G4double Af2, G4double Zf2,
     136                                G4double U, G4double Tmax,
    143137                                const G4FissionParameters & theParam);
    144138   
     139  G4double Ratio(G4double A, G4double A11, G4double B1, G4double A00);
     140  G4double SymmetricRatio(G4int A, G4double A11);
     141  G4double AsymmetricRatio(G4int A, G4double A11);
    145142
    146 
    147   G4double Ratio(const G4double A,const G4double A11,const G4double B1,const G4double A00);
    148   G4double SymmetricRatio(const G4double A,const G4double A11);
    149   G4double AsymmetricRatio(const G4double A,const G4double A11);
    150 
    151 
    152 
    153   G4ThreeVector IsotropicVector(const G4double Magnitude = 1.0);
    154 
     143  G4ThreeVector IsotropicVector(G4double Magnitude = 1.0);
    155144
    156145#ifdef debug
    … …  
    159148#endif
    160149
    161 
    162150};
    163 
    164 
    165151
    166152#endif

  • trunk/source/processes/hadronic/models/de_excitation/fission/include/G4FissionBarrier.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4FissionBarrier.hh,v 1.5 2009/12/16 16:50:07 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4FissionBarrier.hh,v 1.6 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    4141{
    4242public:
    43   G4FissionBarrier() {};
    44   ~G4FissionBarrier() {};
     43  G4FissionBarrier();
     44  ~G4FissionBarrier();
    4545
    4646private:
    … …  
    5252 
    5353public:
    54   G4double FissionBarrier(const G4int A, const G4int Z, const G4double U);
     54  G4double FissionBarrier(G4int A, G4int Z, G4double U);
    5555
    5656
    5757private:
    5858
    59   G4double BarashenkovFissionBarrier(const G4int A, const G4int Z);
     59  G4double BarashenkovFissionBarrier(G4int A, G4int Z);
    6060 
    61   G4double SellPlusPairingCorrection(const G4int Z, const G4int N)
     61  inline G4double SellPlusPairingCorrection(G4int Z, G4int N)
    6262  {
    6363    G4CameronShellPlusPairingCorrections* SPtr = G4CameronShellPlusPairingCorrections::GetInstance();

  • trunk/source/processes/hadronic/models/de_excitation/fission/include/G4FissionLevelDensityParameter.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4FissionLevelDensityParameter.hh,v 1.3 2006/06/29 20:13:23 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4FissionLevelDensityParameter.hh,v 1.4 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    4545{
    4646public:
    47   G4FissionLevelDensityParameter() {};
    48   virtual ~G4FissionLevelDensityParameter() {};
     47  G4FissionLevelDensityParameter();
     48  virtual ~G4FissionLevelDensityParameter();
    4949
    5050private: 
    … …  
    5656 
    5757public:
    58   G4double LevelDensityParameter(const G4int A,const G4int Z,const G4double U) const;
     58  G4double LevelDensityParameter(G4int A, G4int Z, G4double U) const;
    5959
    6060 

  • trunk/source/processes/hadronic/models/de_excitation/fission/include/G4FissionParameters.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4FissionParameters.hh,v 1.3 2006/06/29 20:13:25 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4FissionParameters.hh,v 1.4 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    4444public:
    4545  // Only available constructor
    46   G4FissionParameters(const G4int A, const G4int Z, const G4double ExEnergy, const G4double FissionBarrier);
     46  G4FissionParameters(G4int A, G4int Z, G4double ExEnergy, G4double FissionBarrier);
    4747
    48   ~G4FissionParameters() {}; 
     48  ~G4FissionParameters(); 
    4949
    5050private: 
    5151  // Default constructor
    52   G4FissionParameters() {};
     52  G4FissionParameters();
    5353
    5454  // Copy constructor

  • trunk/source/processes/hadronic/models/de_excitation/fission/include/G4FissionProbability.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4FissionProbability.hh,v 1.3 2006/06/29 20:13:27 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4FissionProbability.hh,v 1.4 2010/11/17 20:22:46 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    4646public:
    4747  // Default constructor
    48   G4FissionProbability() {};
     48  G4FissionProbability();
    4949
    50   ~G4FissionProbability() {}; 
     50  ~G4FissionProbability(); 
    5151
    5252private: 
    … …  
    5959 
    6060public:
    61         G4double EmissionProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy);
     61  G4double EmissionProbability(const G4Fragment & fragment, G4double MaximalKineticEnergy);
    6262
    6363private:
    64         G4EvaporationLevelDensityParameter theEvapLDP;
    65         G4FissionLevelDensityParameter theFissLDP;
     64  G4EvaporationLevelDensityParameter theEvapLDP;
     65  G4FissionLevelDensityParameter theFissLDP;
    6666
    6767

  • trunk/source/processes/hadronic/models/de_excitation/fission/include/G4ParaFissionModel.hh

    r819 r1347  
    2525//
    2626#ifndef G4ParaFissionModel_h
    27 #define G4ParaFissionModel_h
     27#define G4ParaFissionModel_h 1
    2828
    2929#include "G4CompetitiveFission.hh"
    … …  
    4949    SetMaxEnergy( 60.*MeV );
    5050  }
     51
     52  ~G4ParaFissionModel() {};
    5153 
    5254  virtual G4HadFinalState* ApplyYourself(const G4HadProjectile& aTrack,
    53                                               G4Nucleus& theNucleus)
     55                                        G4Nucleus& theNucleus)
    5456  {
    5557    theParticleChange.Clear();
    … …  
    5961    // prepare the fragment
    6062
    61     G4Fragment anInitialState;
    62     G4double anA = theNucleus.GetN();
    63     G4double aZ = theNucleus.GetZ();
    64     G4double nucMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(G4int(aZ) ,G4int(anA));
    65      
    66     anA += aTrack.GetDefinition()->GetBaryonNumber();
    67     aZ += aTrack.GetDefinition()->GetPDGCharge();
     63    G4int A = theNucleus.GetA_asInt();
     64    G4int Z = theNucleus.GetZ_asInt();
     65    G4double nucMass =
     66      G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
    6867     
    6968    G4int numberOfEx = aTrack.GetDefinition()->GetBaryonNumber();
    70     G4int numberOfCh = G4int(std::abs(aTrack.GetDefinition()->GetPDGCharge()));
     69    G4int numberOfCh = G4int(aTrack.GetDefinition()->GetPDGCharge() + 0.5);
    7170    G4int numberOfHoles = 0;
     71
     72    A += numberOfEx;
     73    Z += numberOfCh;
    7274     
    73     G4ThreeVector exciton3Momentum = aTrack.Get4Momentum().vect();
    74     G4double compoundMass = aTrack.GetTotalEnergy();
    75     compoundMass += nucMass;
    76     compoundMass = std::sqrt(compoundMass*compoundMass - exciton3Momentum*exciton3Momentum);
    77     G4LorentzVector fragment4Momentum(exciton3Momentum,
    78                                std::sqrt(exciton3Momentum.mag2()+compoundMass*compoundMass));
    79    
    80     anInitialState.SetA(anA);
    81     anInitialState.SetZ(aZ);
    82     anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
    83     anInitialState.SetNumberOfCharged(numberOfCh);
    84     anInitialState.SetNumberOfHoles(numberOfHoles);
    85     anInitialState.SetMomentum(fragment4Momentum);
     75    G4LorentzVector v = aTrack.Get4Momentum() + G4LorentzVector(0.0,0.0,0.0,nucMass);
     76    G4Fragment anInitialState(Z,A,v);
     77    anInitialState.SetNumberOfExcitedParticle(numberOfEx,numberOfCh);
     78    anInitialState.SetNumberOfHoles(0,0);
    8679
    8780    // do the fission
    … …  
    10396        {
    10497          G4ReactionProduct* rp0 = (*theExcitationResult)[j];
    105           G4DynamicParticle* p0 = new G4DynamicParticle;
    106           p0->SetDefinition(rp0->GetDefinition() );
    107           p0->SetMomentum(rp0->GetMomentum() );
     98          G4DynamicParticle* p0 =
     99            new G4DynamicParticle(rp0->GetDefinition(),rp0->GetMomentum());
    108100          theParticleChange.AddSecondary(p0);
    109101          delete rp0;
    … …  
    114106      {
    115107        // add secondary
    116         G4DynamicParticle* p0 = new G4DynamicParticle;
    117         p0->SetDefinition(aFragment->GetParticleDefinition());
    118         p0->SetMomentum(aFragment->GetMomentum().vect());
     108        G4DynamicParticle* p0 =
     109          new G4DynamicParticle(aFragment->GetParticleDefinition(),
     110                                aFragment->GetMomentum());
    119111        theParticleChange.AddSecondary(p0);
    120112      }

  • trunk/source/processes/hadronic/models/de_excitation/fission/include/G4VFissionBarrier.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4VFissionBarrier.hh,v 1.4 2006/06/29 20:13:31 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4VFissionBarrier.hh,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
    … …  
    4141{
    4242public:
    43   G4VFissionBarrier() {};
    44   virtual ~G4VFissionBarrier() {};
     43  G4VFissionBarrier();
     44  virtual ~G4VFissionBarrier();
    4545
    4646private:
    … …  
    5252 
    5353public:
    54   virtual G4double FissionBarrier(const G4int A, const G4int Z,const G4double U) = 0;
     54  virtual G4double FissionBarrier(G4int A, G4int Z,const G4double U) = 0;
    5555 
    5656

  • 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.