- Timestamp:
- Apr 6, 2009, 12:30:29 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model
- Files:
-
- 34 edited
-
include/G4GNASHTransitions.hh (modified) (1 diff)
-
include/G4LowEIonFragmentation.hh (modified) (1 diff)
-
include/G4PreCompoundAlpha.hh (modified) (4 diffs)
-
include/G4PreCompoundDeuteron.hh (modified) (3 diffs)
-
include/G4PreCompoundEmission.hh (modified) (2 diffs)
-
include/G4PreCompoundEmission.icc (modified) (3 diffs)
-
include/G4PreCompoundFragment.hh (modified) (3 diffs)
-
include/G4PreCompoundFragmentVector.hh (modified) (2 diffs)
-
include/G4PreCompoundFragmentVector.icc (modified) (2 diffs)
-
include/G4PreCompoundHe3.hh (modified) (3 diffs)
-
include/G4PreCompoundIon.hh (modified) (2 diffs)
-
include/G4PreCompoundModel.hh (modified) (6 diffs)
-
include/G4PreCompoundNeutron.hh (modified) (4 diffs)
-
include/G4PreCompoundNucleon.hh (modified) (3 diffs)
-
include/G4PreCompoundParameters.hh (modified) (2 diffs)
-
include/G4PreCompoundProton.hh (modified) (3 diffs)
-
include/G4PreCompoundTransitions.hh (modified) (4 diffs)
-
include/G4PreCompoundTriton.hh (modified) (3 diffs)
-
include/G4VPreCompoundFragment.hh (modified) (6 diffs)
-
include/G4VPreCompoundFragment.icc (modified) (2 diffs)
-
include/G4VPreCompoundIon.hh (modified) (2 diffs)
-
include/G4VPreCompoundNucleon.hh (modified) (1 diff)
-
include/G4VPreCompoundTransitions.hh (modified) (2 diffs)
-
src/G4PreCompoundEmission.cc (modified) (7 diffs)
-
src/G4PreCompoundFragment.cc (modified) (5 diffs)
-
src/G4PreCompoundFragmentVector.cc (modified) (3 diffs)
-
src/G4PreCompoundIon.cc (modified) (2 diffs)
-
src/G4PreCompoundModel.cc (modified) (9 diffs)
-
src/G4PreCompoundNucleon.cc (modified) (4 diffs)
-
src/G4PreCompoundParameters.cc (modified) (1 diff)
-
src/G4PreCompoundTransitions.cc (modified) (5 diffs)
-
src/G4VPreCompoundFragment.cc (modified) (4 diffs)
-
src/G4VPreCompoundIon.cc (modified) (2 diffs)
-
src/G4VPreCompoundNucleon.cc (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4GNASHTransitions.hh
r819 r962 52 52 virtual G4Fragment PerformTransition(const G4Fragment & aFragment); 53 53 54 //JMQ 03/01/08 55 public: 56 // inline G4double GetTransitionProb1() const 57 G4double GetTransitionProb1() 58 {return 0;} 59 // inline G4double GetTransitionProb2() const 60 G4double GetTransitionProb2() 61 {return 0;} 62 // inline G4double GetTransitionProb3() const 63 G4double GetTransitionProb3() 64 {return 0;} 65 // 66 54 67 55 68 }; -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4LowEIonFragmentation.hh
r819 r962 26 26 // 27 27 // $Id: G4LowEIonFragmentation.hh,v 1.3 2006/06/29 20:58:04 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by H.P. Wellisch -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundAlpha.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // by V. Lara 26 27 // 27 // $Id: G4PreCompoundAlpha.hh,v 1.6 2007/10/01 10:41:59 ahoward Exp $ 28 // GEANT4 tag $Name: $ 29 // 30 // by V. Lara 28 //J. M. Quesada (July 08) 29 31 30 32 31 #ifndef G4PreCompoundAlpha_h … … 36 35 #include "G4ReactionProduct.hh" 37 36 #include "G4Alpha.hh" 38 39 37 #include "G4AlphaCoulombBarrier.hh" 40 38 #include "G4PreCompoundParameters.hh" 41 39 42 40 class G4PreCompoundAlpha : public G4PreCompoundIon … … 47 45 48 46 // copy constructor 49 G4PreCompoundAlpha(const G4PreCompoundAlpha &right): G4PreCompoundIon(right) 47 G4PreCompoundAlpha(const G4PreCompoundAlpha &right): G4PreCompoundIon(right){} 50 48 51 49 // destructor … … 66 64 67 65 68 G4ReactionProduct * GetReactionProduct() const 69 { 70 G4ReactionProduct * theReactionProduct = 71 new G4ReactionProduct(G4Alpha::AlphaDefinition()); 72 theReactionProduct->SetMomentum(GetMomentum().vect()); 73 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 74 #ifdef PRECOMPOUND_TEST 75 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 76 #endif 77 return theReactionProduct; 78 } 79 66 G4ReactionProduct * GetReactionProduct() const; 67 80 68 private: 81 69 82 // added Rj method according to literature and JMQ 83 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) 84 { 85 G4double rj = 1.0; 86 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2)*(NumberParticles-3); 87 if(denominator !=0) rj = 6.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); //JMQ 23/8/07 70 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged); 88 71 89 return rj; 90 } 72 virtual G4double CrossSection(const G4double K) ; 73 74 virtual G4double FactorialFactor(const G4double N, const G4double P); 75 76 virtual G4double CoalescenceFactor(const G4double A); 77 78 G4double GetOpt0(const G4double K); 79 G4double GetOpt12(const G4double K); 80 G4double GetOpt34(const G4double K); 81 82 G4double GetAlpha(); 83 84 G4double GetBeta(); 85 86 //data members 87 88 G4AlphaCoulombBarrier theAlphaCoulombBarrier; 89 G4double ResidualA; 90 G4double ResidualZ; 91 G4double theA; 92 G4double theZ; 93 G4double ResidualAthrd; 94 G4double FragmentA; 95 G4double FragmentAthrd; 96 97 }; 98 #endif 91 99 92 100 93 virtual G4double GetAlpha()94 {95 G4double C = 0.0;96 G4double aZ = GetZ() + GetRestZ();97 if (aZ <= 30)98 {99 C = 0.10;100 }101 else if (aZ <= 50)102 {103 C = 0.1 + -((aZ-50.)/20.)*0.02;104 }105 else if (aZ < 70)106 {107 C = 0.08 + -((aZ-70.)/20.)*0.02;108 }109 else110 {111 C = 0.06;112 }113 return 1.0+C;114 }115 101 116 virtual G4double GetBeta()117 {118 return -GetCoulombBarrier();119 }120 102 121 virtual G4double FactorialFactor(const G4double N, const G4double P)122 {123 return124 (N-4.0)*(P-3.0)*(125 (((N-3.0)*(P-2.0))/2.0) *(126 (((N-2.0)*(P-1.0))/3.0) *(127 (((N-1.0)*P)/2.0)128 )129 )130 );131 }132 103 133 virtual G4double CoalescenceFactor(const G4double A)134 {135 return 4096.0/(A*A*A);136 }137 private:138 139 G4AlphaCoulombBarrier theAlphaCoulombBarrier;140 141 };142 143 #endif144 104 145 105 106 107 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundDeuteron.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // by V. Lara 26 27 // 27 // $Id: G4PreCompoundDeuteron.hh,v 1.7 2007/10/01 10:41:59 ahoward Exp $ 28 // GEANT4 tag $Name: $ 29 // 30 // by V. Lara 28 //J. M. Quesada (July 08) 31 29 32 30 #ifndef G4PreCompoundDeuteron_h … … 36 34 #include "G4ReactionProduct.hh" 37 35 #include "G4Deuteron.hh" 38 39 36 #include "G4DeuteronCoulombBarrier.hh" 37 #include "G4PreCompoundParameters.hh" 40 38 41 39 … … 66 64 67 65 68 G4ReactionProduct * GetReactionProduct() const 69 { 70 G4ReactionProduct * theReactionProduct = 71 new G4ReactionProduct(G4Deuteron::DeuteronDefinition()); 72 theReactionProduct->SetMomentum(GetMomentum().vect()); 73 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 74 #ifdef PRECOMPOUND_TEST 75 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 76 #endif 77 return theReactionProduct; 78 } 66 G4ReactionProduct * GetReactionProduct() const; 67 79 68 80 69 private: 81 70 82 // added Rj method according to literature and JMQ - formula from Jose Quesada 83 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) 84 { 85 G4double rj = 1.0; 86 G4double denominator = NumberParticles*(NumberParticles-1); 87 if(denominator != 0) rj = 2.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged))/static_cast<G4double>(denominator); //JMQ re-corrected 23/8/07 88 // return static_cast<G4double>(NumberCharged*(NumberCharged-1))/static_cast<G4double>(NumberParticles*(NumberParticles-1)); 71 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged); 89 72 90 return rj; 91 } 73 virtual G4double CrossSection(const G4double K) ; 92 74 75 virtual G4double FactorialFactor(const G4double N, const G4double P); 93 76 77 virtual G4double CoalescenceFactor(const G4double A); 94 78 95 virtual G4double GetAlpha() 96 { 97 G4double C = 0.0; 98 G4double aZ = GetZ() + GetRestZ(); 99 if (aZ >= 70) 100 { 101 C = 0.10; 102 } 103 else 104 { 105 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 106 } 107 return 1.0 + C/2.0; 108 } 109 110 virtual G4double GetBeta() 111 { 112 return -GetCoulombBarrier(); 113 } 79 G4double GetOpt0(const G4double K); 80 G4double GetOpt12(const G4double K); 81 G4double GetOpt34(const G4double K); 114 82 115 virtual G4double FactorialFactor(const G4double N, const G4double P) 116 { 117 return (N-1.0)*(N-2.0)*(P-1.0)*P/2.0; 118 } 83 G4double GetAlpha(); 119 84 120 virtual G4double CoalescenceFactor(const G4double A) 121 { 122 return 16.0/A; 123 } 124 private: 125 126 G4DeuteronCoulombBarrier theDeuteronCoulombBarrier; 85 G4double GetBeta(); 86 87 //data members 88 89 G4DeuteronCoulombBarrier theDeuteronCoulombBarrier; 90 G4double ResidualA; 91 G4double ResidualZ; 92 G4double theA; 93 G4double theZ; 94 G4double ResidualAthrd; 95 G4double FragmentA; 96 G4double FragmentAthrd; 127 97 128 98 }; 99 #endif 129 100 130 #endif131 132 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PreCompoundEmission.hh,v 1. 3 2006/06/29 20:58:10 gunterExp $27 // GEANT4 tag $Name: $26 // $Id: G4PreCompoundEmission.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 28 28 // 29 29 // Hadronic Process: Nuclear Preequilibrium 30 30 // by V. Lara 31 // 32 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 33 // cross section option 34 // JMQ (06 September 2008) Also external choice has been added for: 35 // - superimposed Coulomb barrier (if useSICB=true) 31 36 32 37 #ifndef G4PreCompoundEmission_h … … 78 83 const G4double E, const G4double Ef) const; 79 84 85 G4double factorial(G4double a) const; 86 87 //for inverse cross section choice 88 public: 89 inline void SetOPTxs(G4int); 90 //for superimposed CoulomBarrier for inverse cross sections 91 inline void UseSICB(G4bool); 92 93 80 94 //============== 81 95 // Data Members -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.icc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PreCompoundEmission.icc,v 1.5 2009/02/10 16:01:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 28 // 29 // 30 // Author: V.Lara 31 // 32 // Modified: 33 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 34 // cross section option 35 // JMQ (06 September 2008) Also external choice has been added for: 36 // - superimposed Coulomb barrier (if useSICB=true) 37 26 38 inline void G4PreCompoundEmission::Initialize(const G4Fragment & aFragment) 27 39 { … … 29 41 return; 30 42 } 31 32 43 33 44 inline G4double G4PreCompoundEmission::GetTotalProbability(const G4Fragment & aFragment) … … 46 57 return; 47 58 } 59 60 //for inverse cross section choice 61 inline void G4PreCompoundEmission::SetOPTxs(G4int opt) 62 { 63 theFragmentsVector->SetOPTxs(opt); 64 return; 65 } 66 //for superimposed Coumlomb Barrier for inverse cross sections 67 inline void G4PreCompoundEmission::UseSICB(G4bool use) 68 { 69 theFragmentsVector->UseSICB(use); 70 return; 71 } 72 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragment.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // by V. Lara 26 //J. M. Quesada (August 2008). 27 //Based on previous work by V. Lara 28 // 29 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 30 // cross section option (default OPTxs=2) 31 // JMQ (06 September 2008) Also external choice has been added for 32 // superimposed Coulomb barrier (if useSICB=true, default false) 33 27 34 28 35 #ifndef G4PreCompoundFragment_h … … 62 69 63 70 // Initialization method 64 void Initialize(const G4Fragment & aFragment);71 // void Initialize(const G4Fragment & aFragment); 65 72 66 73 // ================================================ … … 84 91 virtual G4double 85 92 ProbabilityDistributionFunction(const G4double K, 86 const G4Fragment & aFragment) = 0; 93 const G4Fragment & aFragment) = 0; 94 87 95 }; 88 96 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundFragmentVector.hh,v 1. 3 2006/06/29 20:58:18 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4PreCompoundFragmentVector.hh,v 1.6 2009/02/10 16:01:37 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear Preequilibrium 31 // by V. Lara 31 // by V. Lara 32 // 33 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 34 // cross section option 35 // JMQ (06 September 2008) Also external choice has been added for: 36 // - superimposed Coulomb barrier (if useSICB=true) 32 37 33 38 #ifndef G4PreCompoundFragmentVector_h 34 39 #define G4PreCompoundFragmentVector_h 1 35 40 36 37 41 #include "G4VPreCompoundFragment.hh" 38 39 40 42 41 43 class G4PreCompoundFragmentVector … … 69 71 G4double TotalEmissionProbability; 70 72 73 //for inverse cross section choice 74 public: 75 inline void SetOPTxs(G4int); 76 //for superimposed CoulomBarrier for inverse cross sections 77 inline void UseSICB(G4bool); 78 79 71 80 }; 72 81 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.icc
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundFragmentVector.icc,v 1. 3 2006/06/29 20:58:20 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4PreCompoundFragmentVector.icc,v 1.4 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear Preequilibrium 31 31 // by V. Lara 32 32 // 33 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 34 // cross section option 35 // JMQ (06 September 2008) Also external choice has been added for: 36 // - superimposed Coulomb barrier (if useSICB=true) 33 37 34 38 inline G4PreCompoundFragmentVector::G4PreCompoundFragmentVector(pcfvector * avector) : … … 64 68 return; 65 69 } 70 71 //for inverse cross section choice 72 inline void G4PreCompoundFragmentVector:: 73 SetOPTxs(G4int opt) 74 { 75 for (pcfvector::iterator i = theChannels->begin(); i != theChannels->end(); ++i) 76 (*i)->SetOPTxs(opt); 77 return; 78 } 79 80 //for superimposed Coulomb Barrier for inverse cross sections 81 inline void G4PreCompoundFragmentVector:: 82 UseSICB(G4bool use) 83 { 84 for (pcfvector::iterator i = theChannels->begin(); i != theChannels->end(); ++i) 85 (*i)->UseSICB(use); 86 return; 87 } -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundHe3.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // by V. Lara 26 27 // 27 // $Id: G4PreCompoundHe3.hh,v 1.6 2007/10/01 10:42:00 ahoward Exp $ 28 // GEANT4 tag $Name: $ 29 // 30 // by V. Lara 28 //J. M. Quesada (July 08) 31 29 32 30 #ifndef G4PreCompoundHe3_h … … 36 34 #include "G4ReactionProduct.hh" 37 35 #include "G4He3.hh" 38 39 36 #include "G4He3CoulombBarrier.hh" 40 37 #include "G4PreCompoundParameters.hh" 41 38 42 39 class G4PreCompoundHe3 : public G4PreCompoundIon … … 66 63 67 64 68 G4ReactionProduct * GetReactionProduct() const 69 { 70 G4ReactionProduct * theReactionProduct = 71 new G4ReactionProduct(G4He3::He3Definition()); 72 theReactionProduct->SetMomentum(GetMomentum().vect()); 73 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 74 #ifdef PRECOMPOUND_TEST 75 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 76 #endif 77 return theReactionProduct; 78 } 79 65 G4ReactionProduct * GetReactionProduct() const; 66 67 80 68 private: 81 69 82 // added Rj method according to literature and JMQ 83 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) 84 { 85 G4double rj = 1.0; 86 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2); 87 if(denominator != 0) rj = 3.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged))/static_cast<G4double>(denominator); //JMQ 23/8/07 70 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged); 88 71 89 return rj; 90 } 72 virtual G4double CrossSection(const G4double K) ; 73 74 virtual G4double FactorialFactor(const G4double N, const G4double P); 75 76 virtual G4double CoalescenceFactor(const G4double A); 77 78 G4double GetOpt0(const G4double K); 79 G4double GetOpt12(const G4double K); 80 G4double GetOpt34(const G4double K); 81 82 G4double GetAlpha(); 83 84 G4double GetBeta(); 85 86 //data members 87 88 G4He3CoulombBarrier theHe3CoulombBarrier; 89 G4double ResidualA; 90 G4double ResidualZ; 91 G4double theA; 92 G4double theZ; 93 G4double ResidualAthrd; 94 G4double FragmentA; 95 G4double FragmentAthrd; 96 97 98 }; 99 #endif 91 100 92 101 93 102 94 virtual G4double GetAlpha()95 {96 G4double C = 0.0;97 G4double aZ = GetZ() + GetRestZ();98 if (aZ <= 30)99 {100 C = 0.10;101 }102 else if (aZ <= 50)103 {104 C = 0.1 + -((aZ-50.)/20.)*0.02;105 }106 else if (aZ < 70)107 {108 C = 0.08 + -((aZ-70.)/20.)*0.02;109 }110 else111 {112 C = 0.06;113 }114 return 1.0 + C*(4.0/3.0);115 }116 103 117 virtual G4double GetBeta()118 {119 return -GetCoulombBarrier();120 }121 104 122 virtual G4double FactorialFactor(const G4double N, const G4double P)123 {124 return125 (N-3.0)*(P-2.0)*(126 (((N-2.0)*(P-1.0))/2.0) *(127 (((N-1.0)*P)/3.0)128 )129 );130 }131 132 virtual G4double CoalescenceFactor(const G4double A)133 {134 return 243.0/(A*A);135 }136 private:137 138 G4He3CoulombBarrier theHe3CoulombBarrier;139 140 };141 142 #endif143 105 144 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundIon.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // by V. Lara 26 //J. M. Quesada (August 2008). 27 //Based on previous work by V. Lara 28 // 27 29 28 30 #ifndef G4PreCompoundIon_h … … 70 72 } 71 73 72 virtual G4double ProbabilityDistributionFunction(const G4double eKin, 73 74 virtual G4double ProbabilityDistributionFunction(const G4double eKin, 75 const G4Fragment& aFragment); 74 76 75 protected: 76 G4bool IsItPossible(const G4Fragment& aFragment) 77 { 78 G4int pplus = aFragment.GetNumberOfCharged(); 79 G4int pneut = aFragment.GetNumberOfParticles()-pplus; 80 return (pneut >= (GetA()-GetZ()) && pplus >= GetZ()); 81 } 77 private: 78 79 G4bool IsItPossible(const G4Fragment& aFragment) ; 82 80 83 virtual G4double GetAlpha() = 0; 84 virtual G4double GetBeta() = 0; 85 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0; 81 protected: 82 83 virtual G4double CrossSection(const G4double ekin)=0; 84 85 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0; 86 86 87 virtual G4double FactorialFactor(const G4double N, const G4double P) = 0; 88 87 89 virtual G4double CoalescenceFactor(const G4double A) = 0; 88 89 };90 91 }; 90 92 91 93 #endif -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundModel.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundModel.hh,v 1. 3 2006/06/29 20:58:26 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4PreCompoundModel.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by V. Lara … … 36 36 // transport, or any of the string-parton models. 37 37 // Class Description - End 38 // 39 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 40 // cross section option.(default OPTxs=3) 41 // JMQ (06 September 2008) Also external choices have been added for: 42 // - superimposed Coulomb barrier (if useSICB=true, default false) 43 // - "never go back" hipothesis (if useNGB=true, default false) 44 // - soft cutoff from preeq. to equlibrium (if useSCO=true, default false) 45 // - CEM transition probabilities (if useCEMtr=true, default false) 38 46 39 47 #ifndef G4PreCompoundModel_h … … 50 58 #include "Randomize.hh" 51 59 60 //#include "G4PreCompoundEmission.hh" 61 62 #include "G4DynamicParticle.hh" 63 #include "G4ReactionProductVector.hh" 64 #include "G4ReactionProduct.hh" 65 #include "G4ParticleTypes.hh" 66 #include "G4ParticleTable.hh" 52 67 53 68 //#define debug … … 56 71 class G4PreCompoundModel : public G4VPreCompoundModel 57 72 { 73 58 74 public: 59 60 75 G4PreCompoundModel(G4ExcitationHandler * const value) : 61 G4VPreCompoundModel(value), useHETCEmission(false), useGNASHTransition(false) {} 76 G4VPreCompoundModel(value), useHETCEmission(false), useGNASHTransition(false), 77 OPTxs(3), useSICB(false), useNGB(false), useSCO(false), useCEMtr(false) {} 62 78 63 79 ~G4PreCompoundModel() {} … … 89 105 inline void UseDefaultTransition() { useGNASHTransition = false; } 90 106 107 //for cross section selection 108 inline void SetOPTxs(G4int opt) { OPTxs = opt; } 109 //for the rest of external choices 110 inline void UseSICB() { useSICB = true; } 111 inline void UseNGB() { useNGB = true; } 112 inline void UseSCO() { useSCO = true; } 113 inline void UseCEMtr() { useCEMtr = true; } 91 114 private: 92 115 93 116 void PerformEquilibriumEmission(const G4Fragment & aFragment, 94 117 G4ReactionProductVector * theResult) const; 118 119 private: 95 120 96 121 #ifdef debug … … 104 129 //============== 105 130 131 132 106 133 G4bool useHETCEmission; 107 134 G4bool useGNASHTransition; 135 136 //for cross section options 137 G4int OPTxs; 138 //for the rest of external choices 139 G4bool useSICB; 140 G4bool useNGB; 141 G4bool useSCO; 142 G4bool useCEMtr; 143 144 108 145 G4HadFinalState theResult; 109 146 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundNeutron.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundNeutron.hh,v 1.7 2007/10/01 10:42:00 ahoward Exp $ 28 // GEANT4 tag $Name: $ 29 // 30 // by V. Lara 27 //J. M. Quesada (July 08) 31 28 32 29 … … 39 36 #include "G4PreCompoundParameters.hh" 40 37 #include "Randomize.hh" 41 42 38 #include "G4NeutronCoulombBarrier.hh" 43 39 … … 47 43 public: 48 44 // default constructor 49 G4PreCompoundNeutron() : G4PreCompoundNucleon(1,0,&theNeutronCoulom Barrier,"Neutron") {}45 G4PreCompoundNeutron() : G4PreCompoundNucleon(1,0,&theNeutronCoulombBarrier,"Neutron") {} 50 46 51 47 // copy constructor … … 68 64 69 65 70 G4ReactionProduct * GetReactionProduct() const 71 { 72 G4ReactionProduct * theReactionProduct = 73 new G4ReactionProduct(G4Neutron::NeutronDefinition()); 74 theReactionProduct->SetMomentum(GetMomentum().vect()); 75 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 76 #ifdef PRECOMPOUND_TEST 77 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 78 #endif 79 return theReactionProduct; 80 } 66 G4ReactionProduct * GetReactionProduct() const; 81 67 82 68 private: 83 69 84 // added Rj method according to literature and JMQ - formula from Jose Quesada 85 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) 86 { 87 G4double rj = 1.0; 88 if(NumberParticles != 0) rj = static_cast<G4double>(NumberParticles - NumberCharged)/static_cast<G4double>(NumberParticles); 89 return rj; 90 } 70 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged); 71 72 virtual G4double CrossSection(const G4double K); 91 73 92 74 93 virtual G4double GetAlpha()94 {95 return 0.76+2.2/std::pow(GetRestA(),1.0/3.0);96 } 75 G4double GetOpt0(const G4double K); 76 G4double GetOpt12(const G4double K); 77 G4double GetOpt34(const G4double K); 78 97 79 98 virtual G4double GetBeta() 99 { 100 return (2.12/std::pow(GetRestA(),2.0/3.0)-0.05)*MeV/GetAlpha(); 101 } 80 G4double GetAlpha(); 102 81 103 virtual G4bool IsItPossible(const G4Fragment& aFragment) 104 { 105 return ((aFragment.GetNumberOfParticles()-aFragment.GetNumberOfCharged()) >= 1); 106 } 82 G4double GetBeta(); 107 83 84 //data members 108 85 109 private: 110 111 G4NeutronCoulombBarrier theNeutronCoulomBarrier; 86 G4NeutronCoulombBarrier theNeutronCoulombBarrier; 87 G4double ResidualA; 88 G4double ResidualZ; 89 G4double theA; 90 G4double theZ; 91 G4double ResidualAthrd; 92 G4double FragmentA; 93 G4double FragmentAthrd; 112 94 95 113 96 }; 114 97 98 115 99 #endif 116 100 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundNucleon.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // by V. Lara 26 //J. M. Quesada (August 2008). 27 //Based on previous work by V. Lara 28 // 29 27 30 28 31 #ifndef G4PreCompoundNucleon_h … … 46 49 G4PreCompoundNucleon(const G4double anA, 47 50 const G4double aZ, 48 G4VCoulombBarrier* aCoulombBarrier, 49 const G4String & aName): 50 G4PreCompoundFragment(anA,aZ,aCoulombBarrier,aName) {} 51 51 G4VCoulombBarrier* aCoulombBarrier, 52 const G4String & aName) : 53 G4PreCompoundFragment(anA,aZ,aCoulombBarrier,aName) {} 54 55 52 56 virtual ~G4PreCompoundNucleon() {} 53 57 … … 72 76 virtual G4double ProbabilityDistributionFunction(const G4double eKin, 73 77 const G4Fragment& aFragment); 78 79 private: 74 80 81 G4bool IsItPossible(const G4Fragment&) ; 75 82 76 protected:83 protected: 77 84 78 // added Rj method according to literature and JMQ 85 virtual G4double CrossSection(const G4double ekin)=0; 86 79 87 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0; 80 virtual G4double GetAlpha() = 0; 81 virtual G4double GetBeta() = 0; 82 virtual G4bool IsItPossible(const G4Fragment&) = 0; 83 }; 88 89 }; 84 90 85 91 #endif -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundParameters.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundParameters.hh,v 1. 3 2006/06/29 20:58:32 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4PreCompoundParameters.hh,v 1.5 2008/05/08 10:34:25 quesada Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by V. Lara 31 // 32 //J. M. Quesada (Apr. 2008) Level density set to A/10 at preequilibrium 31 33 32 34 … … 42 44 43 45 // default constructor 44 G4PreCompoundParameters() : theLevelDensity(0.125/MeV), 46 // G4PreCompoundParameters() : theLevelDensity(0.125/MeV), 47 //JMQ level density parameter set to A/10 at preequilibrium 48 G4PreCompoundParameters() : theLevelDensity(0.10/MeV), 45 49 r0(1.5*fermi),Transitions_r0(0.6*fermi),FermiEnergy(35.0*MeV) 46 50 {} -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundProton.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundProton.hh,v 1.7 2007/10/01 10:42:00 ahoward Exp $ 28 // GEANT4 tag $Name: $ 29 // 30 // by V. Lara 27 //J. M. Quesada (July 08) 31 28 32 29 #ifndef G4PreCompoundProton_h … … 38 35 #include "G4PreCompoundParameters.hh" 39 36 #include "Randomize.hh" 40 41 37 #include "G4ProtonCoulombBarrier.hh" 42 43 38 44 39 class G4PreCompoundProton : public G4PreCompoundNucleon … … 68 63 69 64 70 G4ReactionProduct * GetReactionProduct() const 71 { 72 G4ReactionProduct * theReactionProduct = 73 new G4ReactionProduct(G4Proton::ProtonDefinition()); 74 theReactionProduct->SetMomentum(GetMomentum().vect()); 75 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 76 #ifdef PRECOMPOUND_TEST 77 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 78 #endif 79 return theReactionProduct; 80 } 65 G4ReactionProduct * GetReactionProduct() const; 66 81 67 82 68 private: 83 69 84 // added Rj method according to literature and JMQ - formula from Jose Quesada 85 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) 86 { 87 G4double rj = 1.0; 88 if(NumberParticles != 0) rj = static_cast<G4double>(NumberCharged)/static_cast<G4double>(NumberParticles); 89 return rj; 90 } 70 //JMQ (Sep. 07) combinatorial factor Rj 71 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged); 72 73 virtual G4double CrossSection(const G4double K) ; 91 74 92 75 93 virtual G4double GetAlpha() 94 { 95 G4double aZ = static_cast<G4double>(GetRestZ()); 96 G4double C = 0.0; 97 if (aZ >= 70) 98 { 99 C = 0.10; 100 } 101 else 102 { 103 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 104 } 105 return 1.0 + C; 106 } 76 G4double GetOpt0(const G4double K); 77 G4double GetOpt1(const G4double K); 78 G4double GetOpt2(const G4double K); 79 G4double GetOpt3(const G4double K); 107 80 108 virtual G4double GetBeta() 109 { 110 return -GetCoulombBarrier(); 111 } 81 G4double GetAlpha(); 112 82 113 virtual G4bool IsItPossible(const G4Fragment& aFragment) 114 { 115 return (aFragment.GetNumberOfCharged() >= 1); 116 } 117 118 119 private: 120 121 G4ProtonCoulombBarrier theProtonCoulombBarrier; 122 83 G4double GetBeta(); 84 85 //data members 86 87 G4ProtonCoulombBarrier theProtonCoulombBarrier; 88 G4double ResidualA; 89 G4double ResidualZ; 90 G4double theA; 91 G4double theZ; 92 G4double ResidualAthrd; 93 G4double FragmentA; 94 G4double FragmentAthrd; 95 96 97 123 98 }; 124 125 99 #endif 126 100 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTransitions.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundTransitions.hh,v 1. 3 2006/06/29 20:58:36 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4PreCompoundTransitions.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by V. Lara 31 // J. M. Quesada . New methods for accessing to individual transition probabilities (landa+, landa-, landa0) from 32 // G4PreCompoundModel.cc 33 31 34 32 35 #ifndef G4PreCompoundTransitions_h … … 50 53 51 54 class G4PreCompoundTransitions : public G4VPreCompoundTransitions 55 52 56 { 53 57 public: … … 55 59 // Calculates transition probabilities with Delta N = +2 (Trans1) -2 (Trans2) and 0 (Trans3) 56 60 57 G4PreCompoundTransitions() : TransitionProb1(0.0), TransitionProb2(0.0), TransitionProb3(0.0) 61 G4PreCompoundTransitions() : TransitionProb1(0.0), TransitionProb2(0.0), TransitionProb3(0.0){} 58 62 59 63 virtual ~G4PreCompoundTransitions() {} … … 81 85 G4double TransitionProb3; 82 86 87 88 //J. M.Quesada (May. 08) 89 public: 90 // inline G4double GetTransitionProb1() const 91 G4double GetTransitionProb1() 92 { 93 return TransitionProb1; 94 } 95 // inline G4double GetTransitionProb2() const 96 G4double GetTransitionProb2() 97 { 98 return TransitionProb2; 99 } 100 // inline G4double GetTransitionProb3() const 101 G4double GetTransitionProb3() 102 { 103 return TransitionProb3; 104 } 105 83 106 }; 84 107 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTriton.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundTriton.hh,v 1.6 2007/10/01 10:42:00 ahoward Exp $ 28 // GEANT4 tag $Name: $ 27 // by V. Lara 29 28 // 30 // by V. Lara29 //J. M. Quesada (July 08) 31 30 32 31 #ifndef G4PreCompoundTriton_h … … 36 35 #include "G4ReactionProduct.hh" 37 36 #include "G4Triton.hh" 38 39 37 #include "G4TritonCoulombBarrier.hh" 38 #include "G4PreCompoundParameters.hh" 40 39 41 40 … … 53 52 54 53 // operators 55 const G4PreCompoundTriton & operator=(const G4PreCompoundTriton &right) {56 if (&right != this) this->G4PreCompoundIon::operator=(right);57 return *this;58 }54 // const G4PreCompoundTriton & operator=(const G4PreCompoundTriton &right) { 55 // if (&right != this) this->G4PreCompoundIon::operator=(right); 56 // return *this; 57 // } 59 58 60 G4bool operator==(const G4PreCompoundTriton &right) const61 { return G4PreCompoundIon::operator==(right);}59 // G4bool operator==(const G4PreCompoundTriton &right) const 60 // { return G4PreCompoundIon::operator==(right);} 62 61 63 62 64 G4bool operator!=(const G4PreCompoundTriton &right) const65 { return G4PreCompoundIon::operator!=(right);}63 // G4bool operator!=(const G4PreCompoundTriton &right) const 64 // { return G4PreCompoundIon::operator!=(right);} 66 65 67 66 68 G4ReactionProduct * GetReactionProduct() const 69 { 70 G4ReactionProduct * theReactionProduct = 71 new G4ReactionProduct(G4Triton::TritonDefinition()); 72 theReactionProduct->SetMomentum(GetMomentum().vect()); 73 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 74 #ifdef PRECOMPOUND_TEST 75 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 76 #endif 77 return theReactionProduct; 78 } 79 67 G4ReactionProduct * GetReactionProduct() const; 68 69 80 70 private: 81 71 82 // added Rj method according to literature and JMQ 83 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) 84 { 85 G4double rj = 1.0; 86 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2); 87 if(denominator != 0) rj = 3.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); //JMQ 23/8/07 72 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged); 88 73 89 return rj; 90 } 74 virtual G4double CrossSection(const G4double K) ; 91 75 92 virtual G4double GetAlpha() 93 { 94 G4double C = 0.0; 95 G4double aZ = GetZ() + GetRestZ(); 96 if (aZ >= 70) 97 { 98 C = 0.10; 99 } 100 else 101 { 102 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 103 } 104 105 return 1.0 + C/3.0; 106 } 76 virtual G4double FactorialFactor(const G4double N, const G4double P); 107 77 108 virtual G4double GetBeta() 109 { 110 return -GetCoulombBarrier(); 111 } 78 virtual G4double CoalescenceFactor(const G4double A); 112 79 113 virtual G4double FactorialFactor(const G4double N, const G4double P) 114 { 115 return 116 (N-3.0)*(P-2.0)*( 117 (((N-2.0)*(P-1.0))/2.0) *( 118 (((N-1.0)*P)/3.0) 119 ) 120 ); 121 } 80 G4double GetOpt0(const G4double K); 81 G4double GetOpt12(const G4double K); 82 G4double GetOpt34(const G4double K); 122 83 123 virtual G4double CoalescenceFactor(const G4double A) 124 { 125 return 243.0/(A*A); 126 } 127 private: 84 G4double GetAlpha(); 85 86 G4double GetBeta(); 128 87 129 G4TritonCoulombBarrier theTritonCoulombBarrier; 88 //data members 130 89 90 G4TritonCoulombBarrier theTritonCoulombBarrier; 91 G4double ResidualA; 92 G4double ResidualZ; 93 G4double theA; 94 G4double theZ; 95 G4double ResidualAthrd; 96 G4double FragmentA; 97 G4double FragmentAthrd; 131 98 }; 132 99 133 100 #endif 101 102 103 104 105 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4VPreCompoundFragment.hh,v 1. 4 2006/08/20 01:07:28 dennisExp $28 // GEANT4 tag $Name: $27 // $Id: G4VPreCompoundFragment.hh,v 1.10 2009/02/10 16:01:37 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 // by V. Lara 30 // J. M. Quesada (August 2008). 31 // Based on previous work by V. Lara 32 // 33 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 34 // cross section option 35 // JMQ (06 September 2008) Also external choice has been added for: 36 // - superimposed Coulomb barrier (if useSICB=true) 31 37 32 38 #ifndef G4VPreCompoundFragment_h … … 60 66 G4VCoulombBarrier * aCoulombBarrier, 61 67 const G4String & aName); 68 62 69 63 70 virtual ~G4VPreCompoundFragment(); … … 135 142 inline void IncrementStage(); 136 143 144 //for inverse cross section choice 145 inline void SetOPTxs(G4int); 146 //for superimposed Coulomb Barrier for inverse cross sections 147 inline void UseSICB(G4bool); 148 149 137 150 138 151 // ============= 139 152 // Data members 140 153 // ============= 154 141 155 142 156 private: … … 145 159 146 160 G4double theZ; 161 private: 147 162 148 163 G4double theRestNucleusA; … … 155 170 156 171 G4double theBindingEnergy; 157 172 158 173 G4double theMaximalKineticEnergy; 159 174 … … 165 180 G4String theFragmentName; 166 181 167 G4int theStage; 182 G4int theStage; 183 184 protected: 185 //for inverse cross section choice 186 G4int OPTxs; 187 //for superimposed Coulomb Barrier for inverse cross sections 188 G4bool useSICB; 168 189 }; 169 190 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.icc
r819 r962 25 25 // 26 26 // 27 // $Id: G4VPreCompoundFragment.icc,v 1. 4 2006/08/20 01:07:40 dennisExp $28 // GEANT4 tag $Name: $27 // $Id: G4VPreCompoundFragment.icc,v 1.7 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by V. Lara 31 31 // 32 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 33 // cross section option 34 // JMQ (06 September 2008) Also external choice has been added for: 35 // - superimposed Coulomb barrier (if useSICB=true) 32 36 33 37 inline G4double G4VPreCompoundFragment::GetA() const … … 137 141 } 138 142 143 //for inverse cross section choice 144 inline void G4VPreCompoundFragment::SetOPTxs(G4int opt) 145 { 146 OPTxs=opt; 147 return; 148 } 149 150 //for superimposed Coulomb Barrier for inverse cross sections 151 inline void G4VPreCompoundFragment::UseSICB(G4bool use) 152 { 153 useSICB=use; 154 return; 155 } -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundIon.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4VPreCompoundIon.hh,v 1.5 2007/08/23 16:29:01 ahoward Exp $ 28 // GEANT4 tag $Name: $ 27 28 // $Id: G4VPreCompoundIon.hh,v 1.9 2008/09/22 10:18:36 ahoward Exp $ 29 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 30 // 30 31 // by V. Lara 31 32 #ifndef G4VPreCompoundIon_h 33 #define G4VPreCompoundIon_h 1 34 35 #include "G4VPreCompoundFragment.hh" 36 #include "G4VCoulombBarrier.hh" 32 // 33 //J.M. Quesada (Apr. 2008) DUMMY abstract base class for ions. 34 // Not ihherited by anything 35 // Suppressed 37 36 38 37 39 class G4VPreCompoundIon : public G4VPreCompoundFragment40 {41 protected:42 // default constructor43 G4VPreCompoundIon() {}44 38 45 public: 46 47 // copy constructor 48 G4VPreCompoundIon(const G4VPreCompoundIon &right): 49 G4VPreCompoundFragment(right) {} 50 51 // constructor 52 G4VPreCompoundIon(const G4double anA, 53 const G4double aZ, 54 G4VCoulombBarrier* aCoulombBarrier, 55 const G4String & aName): 56 G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName) {} 57 58 virtual ~G4VPreCompoundIon() {} 59 60 // operators 61 const G4VPreCompoundIon & 62 operator=(const G4VPreCompoundIon &right) { 63 if (&right != this) this->G4VPreCompoundFragment::operator=(right); 64 return *this; 65 } 66 67 G4bool operator==(const G4VPreCompoundIon &right) const 68 { return G4VPreCompoundFragment::operator==(right);} 69 70 G4bool operator!=(const G4VPreCompoundIon &right) const 71 { return G4VPreCompoundFragment::operator!=(right);} 72 73 virtual G4double ProbabilityDistributionFunction(const G4double eKin, 74 const G4Fragment& aFragment); 75 76 protected: 77 G4bool IsItPossible(const G4Fragment& aFragment) 78 { 79 G4int pplus = aFragment.GetNumberOfCharged(); 80 G4int pneut = aFragment.GetNumberOfParticles()-pplus; 81 return (pneut >= (GetA()-GetZ()) && pplus >= GetZ()); 82 } 83 84 virtual G4double GetAlpha() = 0; 85 virtual G4double GetBeta() = 0; 86 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0; 87 virtual G4double FactorialFactor(const G4double N, const G4double P) = 0; 88 virtual G4double CoalescenceFactor(const G4double A) = 0; 89 90 }; 91 92 #endif 39 93 40 94 41 … … 98 45 99 46 100 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundNucleon.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4VPreCompoundNucleon.hh,v 1. 4 2007/07/23 09:56:40ahoward Exp $28 // GEANT4 tag $Name: $27 // $Id: G4VPreCompoundNucleon.hh,v 1.8 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by V. Lara 31 //J.M. Quesada (Apr. 2008) DUMMY abstract base class for nucleons. 32 // Not ihherited by anything 33 // Suppressed 31 34 32 #ifndef G4VPreCompoundNucleon_h33 #define G4VPreCompoundNucleon_h 134 35 #include "G4VPreCompoundFragment.hh"36 #include "G4VCoulombBarrier.hh"37 38 39 class G4VPreCompoundNucleon : public G4VPreCompoundFragment40 {41 protected:42 // default constructor43 G4VPreCompoundNucleon() {}44 45 public:46 47 // copy constructor48 G4VPreCompoundNucleon(const G4VPreCompoundNucleon &right):49 G4VPreCompoundFragment(right) {}50 51 // constructor52 G4VPreCompoundNucleon(const G4double anA,53 const G4double aZ,54 G4VCoulombBarrier* aCoulombBarrier,55 const G4String & aName):56 G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName) {}57 58 virtual ~G4VPreCompoundNucleon() {}59 60 // operators61 const G4VPreCompoundNucleon &62 operator=(const G4VPreCompoundNucleon &right) {63 if (&right != this) this->G4VPreCompoundFragment::operator=(right);64 return *this;65 }66 67 G4bool operator==(const G4VPreCompoundNucleon &right) const68 { return G4VPreCompoundFragment::operator==(right);}69 70 G4bool operator!=(const G4VPreCompoundNucleon &right) const71 { return G4VPreCompoundFragment::operator!=(right);}72 73 virtual G4double ProbabilityDistributionFunction(const G4double eKin,74 const G4Fragment& aFragment);75 76 77 protected:78 79 // added Rj method according to literature and JMQ80 virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0;81 virtual G4double GetAlpha() = 0;82 virtual G4double GetBeta() = 0;83 virtual G4bool IsItPossible(const G4Fragment&) = 0;84 85 86 };87 88 #endif -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundTransitions.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 //J. M. Quesada (May 08). New virtual classes have been added 27 // JMQ (06 September 2008) Also external choices have been added for: 28 // - "never go back" hipothesis (useNGB=true) 29 // - CEM transition probabilities (useCEMtr=true) 30 26 31 #ifndef G4VPreCompoundTransitions_hh 27 32 #define G4VPreCompoundTransitions_hh 1 … … 33 38 public: 34 39 35 G4VPreCompoundTransitions() {}40 G4VPreCompoundTransitions():useNGB(false),useCEMtr(false) {} 36 41 virtual ~G4VPreCompoundTransitions() {} 37 42 38 43 virtual G4double CalculateProbability(const G4Fragment& aFragment) = 0; 39 44 virtual G4Fragment PerformTransition(const G4Fragment& aFragment) = 0; 45 //J. M. Quesada (May.08) New virtual classes 46 virtual G4double GetTransitionProb1()=0; 47 virtual G4double GetTransitionProb2()=0; 48 virtual G4double GetTransitionProb3()=0; 40 49 50 // for never go back hypothesis (if useNGB=true, default=false) 51 inline void UseNGB(G4bool use){useNGB=use;} 52 //for use of CEM transition probabilities (if useCEMtr=true, defaut false) 53 inline void UseCEMtr(G4bool use){useCEMtr=use;} 54 55 protected: 56 G4bool useNGB; 57 G4bool useCEMtr; 41 58 }; 42 59 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc
r819 r962 25 25 // 26 26 // 27 // Hadronic Process: Nuclear Preequilibrium 28 // by V. Lara 29 27 // $Id: G4PreCompoundEmission.cc,v 1.19 2009/02/10 16:01:37 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 // 30 // ------------------------------------------------------------------- 31 // 32 // GEANT4 Class file 33 // 34 // 35 // File name: G4PreCompoundEmission 36 // 37 // Author: V.Lara 38 // 39 // Modified: 40 // 30 41 31 42 #include "G4PreCompoundEmission.hh" … … 42 53 43 54 44 55 G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const 45 56 { 46 57 return false; 47 58 } 48 59 49 60 G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const 50 61 { 51 62 return true; … … 95 106 return; 96 107 } 97 98 99 108 100 109 G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment) … … 163 172 // Update nucleus parameters: 164 173 // -------------------------- 174 165 175 // Number of excitons 166 176 aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()- … … 175 185 // Charge 176 186 aFragment.SetZ(theFragment->GetRestZ()); 177 187 178 188 179 189 // Perform Lorentz boosts … … 275 285 } 276 286 277 278 287 G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g, 279 288 const G4double E, const G4double Ef) const 289 { 290 291 G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g); 292 G4double alpha = (p*p+h*h)/(2.0*g); 293 294 if ( (E-alpha) < 0 ) return 0; 295 296 G4double factp=factorial(p); 297 298 G4double facth=factorial(h); 299 300 G4double factph=factorial(p+h-1); 301 302 G4double logConst = (p+h)*std::log(g) - std::log (factph) - std::log(factp) - std::log(facth); 303 304 // initialise values using j=0 305 306 G4double t1=1; 307 G4double t2=1; 308 G4double logt3=(p+h-1) * std::log(E-Aph); 309 G4double tot = std::exp( logt3 + logConst ); 310 311 // and now sum rest of terms 312 G4int j(1); 313 while ( (j <= h) && ((E - alpha - j*Ef) > 0.0) ) 314 { 315 t1 *= -1.; 316 t2 *= (h+1-j)/j; 317 logt3 = (p+h-1) * std::log( E - j*Ef - Aph) + logConst; 318 G4double t3 = std::exp(logt3); 319 tot += t1*t2*t3; 320 j++; 321 } 322 323 return tot; 324 } 325 326 327 328 G4double G4PreCompoundEmission::factorial(G4double a) const 280 329 { 281 330 // Values of factorial function from 0 to 60 … … 349 398 // fact[n] = fact[n-1]*static_cast<G4double>(n); 350 399 // } 351 352 G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g); 353 G4double alpha = (p*p+h*h)/(2.0*g); 354 355 G4double tot = 0.0; 356 for (G4int j = 0; j <= h; j++) 357 { 358 if (E-alpha-static_cast<G4double>(j)*Ef > 0.0) 359 { 360 G4double t1 = std::pow(-1.0, static_cast<G4double>(j)); 361 G4double t2 = fact[static_cast<G4int>(h)]/ (fact[static_cast<G4int>(h)-j]*fact[j]); 362 G4double t3 = E - static_cast<G4double>(j)*Ef - Aph; 363 t3 = std::pow(t3,p+h-1); 364 tot += t1*t2*t3; 365 } 366 else 400 G4double result(0.0); 401 G4int ia = static_cast<G4int>(a); 402 if (ia < factablesize) 403 { 404 result = fact[ia]; 405 } 406 else 407 { 408 result = fact[factablesize-1]; 409 for (G4int n = factablesize; n < ia+1; ++n) 367 410 { 368 break;411 result *= static_cast<G4double>(n); 369 412 } 370 413 } 371 414 372 373 G4double factp(0.0); 374 G4int ph = static_cast<G4int>(p); 375 if (ph < factablesize) 376 { 377 factp = fact[ph]; 378 } 379 else 380 { 381 factp = fact[factablesize-1]; 382 for (G4int n = factablesize; n < ph+1; ++n) 383 { 384 factp *= static_cast<G4double>(n); 385 } 386 } 387 388 G4double facth(0.0); 389 ph = static_cast<G4int>(h); 390 if (ph < factablesize) 391 { 392 facth = fact[ph]; 393 } 394 else 395 { 396 facth = fact[factablesize-1]; 397 for (G4int n = factablesize; n < ph+1; ++n) 398 { 399 facth *= static_cast<G4double>(n); 400 } 401 } 402 G4double factph(0.0); 403 ph = static_cast<G4int>(p+h)-1; 404 if (ph < factablesize) 405 { 406 factph = fact[ph]; 407 } 408 else 409 { 410 factph = fact[factablesize-1]; 411 for (G4int n = factablesize; n < ph+1; ++n) 412 { 413 factph *= static_cast<G4double>(n); 414 } 415 } 416 417 tot *= std::pow(g,p+h)/factph; 418 tot /= factp; 419 tot /= facth; 420 421 return tot; 422 } 423 415 return result; 416 } 424 417 425 418 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragment.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PreCompoundFragment.cc,v 1.8 2009/02/10 16:01:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 26 28 // 27 // by V. Lara 28 29 // J. M. Quesada (August 2008). 30 // Based on previous work by V. Lara 31 // JMQ (06 September 2008) Also external choice has been added for: 32 // - superimposed Coulomb barrier (if useSICB=true) 33 // 29 34 #include "G4PreCompoundFragment.hh" 30 35 … … 41 46 const G4String & aName): 42 47 G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName) 43 {} 48 { 49 } 44 50 45 51 … … 71 77 CalcEmissionProbability(const G4Fragment & aFragment) 72 78 { 73 if (GetEnergyThreshold() <= 0.0) 79 // If theCoulombBarrier effect is included in the emission probabilities 80 //if (GetMaximalKineticEnergy() <= 0.0) 81 G4double limit; 82 if(OPTxs==0 || useSICB) limit= theCoulombBarrier; 83 else limit=0.; 84 if (GetMaximalKineticEnergy() <= limit) 74 85 { 75 86 theEmissionProbability = 0.0; 76 87 return 0.0; 77 88 } 78 // Coulomb barrier is the lower limit 79 // of integration over kinetic energy 80 G4double LowerLimit = theCoulombBarrier; 81 82 // Excitation energy of nucleus after fragment emission is the upper limit 83 // of integration over kinetic energy 84 G4double UpperLimit = this->GetMaximalKineticEnergy(); 89 // If theCoulombBarrier effect is included in the emission probabilities 90 // G4double LowerLimit = 0.; 91 // Coulomb barrier is the lower limit 92 // of integration over kinetic energy 93 G4double LowerLimit = limit; 94 95 // Excitation energy of nucleus after fragment emission is the upper limit of integration over kinetic energy 96 G4double UpperLimit = GetMaximalKineticEnergy(); 85 97 86 98 theEmissionProbability = 87 99 IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment); 88 89 100 return theEmissionProbability; 90 101 } … … 130 141 } 131 142 Total *= (Up-Low)/2.0; 132 133 143 return Total; 134 144 } … … 140 150 GetKineticEnergy(const G4Fragment & aFragment) 141 151 { 142 G4double V = this->GetCoulombBarrier(); 143 G4double Tmax = this->GetMaximalKineticEnergy() - V; 144 152 153 // G4double V = this->GetCoulombBarrier();// alternative way for accessing the Coulomb barrier 154 // //should be equivalent (in fact it is) 155 G4double V; 156 if(OPTxs==0 || useSICB) V= theCoulombBarrier;//let's keep this way for consistency with CalcEmissionProbability method 157 else V=0.; 158 159 G4double Tmax = GetMaximalKineticEnergy() ; 145 160 G4double T(0.0); 146 161 G4double NormalizedProbability(1.0); 147 162 do 148 163 { 149 T = V + G4UniformRand()*Tmax; 150 NormalizedProbability = this->ProbabilityDistributionFunction(T,aFragment)/ 151 this->GetEmissionProbability(); 152 153 } 154 while (G4UniformRand() > NormalizedProbability); 155 164 T =V+ G4UniformRand()*(Tmax-V); 165 NormalizedProbability = ProbabilityDistributionFunction(T,aFragment)/GetEmissionProbability(); 166 } while (G4UniformRand() > NormalizedProbability); 156 167 return T; 157 168 } -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragmentVector.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4PreCompoundFragmentVector.cc,v 1.5 2006/06/29 20:59:21 gunter Exp $ 28 // GEANT4 tag $Name: $ 26 // $Id: G4PreCompoundFragmentVector.cc,v 1.11 2009/02/10 16:01:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 28 // 30 29 // Hadronic Process: Nuclear Preequilibrium … … 86 85 for (i = theChannels->begin(); i != theChannels->end(); ++i) { 87 86 accumulation += (*i)->GetEmissionProbability(); 87 88 88 running.push_back(accumulation); 89 89 } … … 126 126 (*i)->IncrementStage(); 127 127 } 128 } 128 } 129 129 130 130 return theChannels->operator[](ChosenChannel); -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundIon.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // by V. Lara 26 // $Id: G4PreCompoundIon.cc,v 1.16 2009/02/10 16:01:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 28 // 29 // ------------------------------------------------------------------- 30 // 31 // GEANT4 Class file 32 // 33 // 34 // File name: G4PreCompoundIon 35 // 36 // Author: V.Lara 37 // 38 // Modified: 39 // 10.02.2009 J. M. Quesada fixed bug in density level of light fragments 40 // 27 41 28 42 #include "G4PreCompoundIon.hh" 29 43 #include "G4PreCompoundParameters.hh" 30 //#include "G4EvaporationLevelDensityParameter.hh"31 44 45 G4bool G4PreCompoundIon::IsItPossible(const G4Fragment& aFragment) 46 { 47 G4int pplus = aFragment.GetNumberOfCharged(); 48 G4int pneut = aFragment.GetNumberOfParticles()-pplus; 49 return (pneut >= (GetA()-GetZ()) && pplus >= GetZ()); 50 } 32 51 33 52 G4double G4PreCompoundIon:: … … 43 62 G4double H = aFragment.GetNumberOfHoles(); 44 63 G4double N = P + H; 45 46 // g = (6.0/pi2)*a*A 47 // G4EvaporationLevelDensityParameter theLDP; 64 48 65 G4double g0 = (6.0/pi2)*aFragment.GetA() * 49 66 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 50 // theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);67 51 68 G4double g1 = (6.0/pi2)*GetRestA() * 52 69 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 53 // theLDP.LevelDensityParameter(GetRestA(),GetRestZ(),U);54 //no longer used: G4double gj = (6.0/pi2)*GetA() *55 // -----"----- G4PreCompoundParameters::GetAddress()->GetLevelDensity();56 // theLDP.LevelDensityParameter(GetA(),GetZ(),U);57 70 58 G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0; 59 G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1); 60 G4double Aj = GetA()*(GetA()+1.0)/4.0; 71 //JMQ 06/02/209 This is THE BUG that was killing cluster emission 72 // G4double gj = (6.0/pi2)*GetA() * 73 // G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 74 75 G4double gj = g1; 76 77 G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0; 78 79 G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1); 80 81 G4double Aj = GetA()*(GetA()+1.0)/4.0/gj; 82 61 83 62 84 G4double E0 = std::max(0.0,U - A0); 63 85 if (E0 == 0.0) return 0.0; 64 // G4double E1 = std::max(0.0,GetMaximalKineticEnergy() + GetCoulombBarrier() - eKin - A1);65 G4double E1 = (std::max(0.0,GetMaximalKineticEnergy() - eKin - A1))/g1; //JMQ fix66 // G4double Ej = std::max(0.0,U + GetBindingEnergy() - Aj);67 G4double Ej = std::max(0.0,eKin + GetBindingEnergy() - Aj); //JMQ fix68 86 69 G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*(eKin+GetBindingEnergy()))))* 70 GetAlpha()*(eKin+GetBeta())/(r0*std::pow(GetRestA(),1.0/3.0)) * 71 CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P); 87 G4double E1 = (std::max(0.0,GetMaximalKineticEnergy() - eKin - A1)); 88 89 G4double Ej = std::max(0.0,eKin + GetBindingEnergy() -Aj); 90 91 // JMQ 10/02/09 reshaping of the formula (unnecessary std::pow elimitated) 92 G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()* 93 (eKin+GetBindingEnergy()))))/(pi * r0 * r0 *r0* GetRestA())* 94 eKin*CrossSection(eKin)*millibarn* 95 CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P)* 96 GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()); 72 97 73 98 G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0); 99 100 G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0; 74 101 75 // G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0; 76 G4double pC = std::pow((g1*Ej)/(g0*E0),GetA()-1.0)*(g1/g0)/E0; //JMQ fix 77 78 G4double Probability = pA * pB * pC * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()); 102 G4double Probability = pA * pB * pC; 79 103 80 104 return Probability; -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundModel.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4PreCompoundModel.cc,v 1.1 1 2007/10/11 14:19:36ahoward Exp $28 // GEANT4 tag $Name: $27 // $Id: G4PreCompoundModel.cc,v 1.17 2008/12/09 14:09:59 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by V. Lara 31 // 32 //J. M. Quesada (Apr.08). Several changes. Soft cut-off switched off. 33 //(May. 08). Protection against non-physical preeq. transitional regime has 34 // been set 35 // 36 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 37 // cross section option 38 // JMQ (06 September 2008) Also external choices have been added for: 39 // - superimposed Coulomb barrier (useSICB=true) 40 // - "never go back" hipothesis (useNGB=true) 41 // - soft cutoff from preeq. to equlibrium (useSCO=true) 42 // - CEM transition probabilities (useCEMtr=true) 43 31 44 32 45 #include "G4PreCompoundModel.hh" … … 34 47 #include "G4PreCompoundTransitions.hh" 35 48 #include "G4GNASHTransitions.hh" 49 #include "G4ParticleDefinition.hh" 36 50 37 51 … … 160 174 // Copy of the initial state 161 175 G4Fragment aFragment(theInitialState); 176 177 if (aFragment.GetExcitationEnergy() < 10*eV) 178 { 179 // Perform Equilibrium Emission 180 PerformEquilibriumEmission(aFragment,Result); 181 return Result; 182 } 162 183 163 184 if (aFragment.GetA() < 5) { … … 175 196 if (useHETCEmission) aEmission.SetHETCModel(); 176 197 aEmission.SetUp(theInitialState); 177 198 199 //for cross section options 200 201 if (OPTxs!= 0 && OPTxs!=1 && OPTxs !=2 && OPTxs !=3 && OPTxs !=4 ) { 202 std::ostringstream errOs; 203 errOs << "BAD CROSS SECTION OPTION in G4PreCompoundModel.cc !!" <<G4endl; 204 throw G4HadronicException(__FILE__, __LINE__, errOs.str());} 205 else aEmission.SetOPTxs(OPTxs); 206 207 //for the choice of superimposed Coulomb Barrier for inverse cross sections 208 209 aEmission.UseSICB(useSICB); 210 211 212 //---------- 213 178 214 G4VPreCompoundTransitions * aTransition = 0; 179 215 if (useGNASHTransition) … … 184 220 { 185 221 aTransition = new G4PreCompoundTransitions; 186 } 187 188 222 // for the choice of "never go back" hypothesis and CEM transition probabilities 223 if (useNGB) aTransition->UseNGB(useNGB); 224 if (useCEMtr) aTransition->UseCEMtr(useCEMtr); 225 } 226 189 227 // Main loop. It is performed until equilibrium deexcitation. 228 //G4int fragment=0; 229 190 230 for (;;) { 231 232 //fragment++; 233 //G4cout<<"-------------------------------------------------------------------"<<G4endl; 234 //G4cout<<"Fragment number .. "<<fragment<<G4endl; 235 191 236 // Initialize fragment according with the nucleus parameters 192 237 aEmission.Initialize(aFragment); 193 194 // Equilibrium exciton number 195 // G4EvaporationLevelDensityParameter theLDP; 196 // G4double g = (6.0/pi2)*aFragment.GetA()* 197 // theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()), 198 // aFragment.GetExcitationEnergy()); 238 239 199 240 200 241 G4double g = (6.0/pi2)*aFragment.GetA()* 201 242 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 202 G4int EquilibriumExcitonNumber = static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy()) 203 + 0.5); 204 205 // Loop for transitions, it is performed while there are preequilibrium transitions. 243 244 245 246 247 G4int EquilibriumExcitonNumber = static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy())+ 0.5); 248 // 249 // G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl; 250 // 251 // J. M. Quesada (Jan. 08) equilibrium hole number could be used as preeq.- evap. delimiter (IAEA report) 252 // G4int EquilibriumHoleNumber = static_cast<G4int>(0.2*std::sqrt(g*aFragment.GetExcitationEnergy())+ 0.5); 253 254 // Loop for transitions, it is performed while there are preequilibrium transitions. 206 255 G4bool ThereIsTransition = false; 256 257 // G4cout<<"----------------------------------------"<<G4endl; 258 // G4double NP=aFragment.GetNumberOfParticles(); 259 // G4double NH=aFragment.GetNumberOfHoles(); 260 // G4double NE=aFragment.GetNumberOfExcitons(); 261 // G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl; 262 // G4cout<<"N. excitons="<<NE<<" N. Part="<<NP<<"N. Holes ="<<NH<<G4endl; 263 264 265 //G4int transition=0; 207 266 do 208 267 { 268 //transition++; 269 //G4cout<<"transition number .."<<transition<<G4endl; 270 //G4cout<<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl; 209 271 G4bool go_ahead = false; 272 // soft cutoff criterium as an "ad-hoc" solution to force increase in evaporation 273 // G4double test = static_cast<G4double>(aFragment.GetNumberOfHoles()); 210 274 G4double test = static_cast<G4double>(aFragment.GetNumberOfExcitons()); 211 if (test < EquilibriumExcitonNumber) 212 { 213 test /= static_cast<G4double>(EquilibriumExcitonNumber); 214 test -= 1.0; 215 test = test*test; 216 test /= 0.32; 217 test = 1.0 - std::exp(-test); 218 go_ahead = (G4UniformRand() < test); 219 } 220 275 276 277 if (test < EquilibriumExcitonNumber) go_ahead=true; 278 //J. M. Quesada (Apr. 08): soft-cutoff switched off by default 279 if (useSCO) { 280 if (test < EquilibriumExcitonNumber) 281 // if (test < EquilibriumHoleNumber) 282 { 283 test /= static_cast<G4double>(EquilibriumExcitonNumber); 284 // test /= static_cast<G4double>(EquilibriumHoleNumber); 285 test -= 1.0; 286 test = test*test; 287 test /= 0.32; 288 test = 1.0 - std::exp(-test); 289 go_ahead = (G4UniformRand() < test); 290 291 } 292 } 293 294 //JMQ: WARNING: CalculateProbability MUST be called prior to Get methods !! (O values would be returned otherwise) 295 G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment); 296 G4double P1=aTransition->GetTransitionProb1(); 297 G4double P2=aTransition->GetTransitionProb2(); 298 G4double P3=aTransition->GetTransitionProb3(); 299 // G4cout<<"P1="<<P1<<" P2="<<P2<<" P3="<<P3<<G4endl; 300 301 302 //J.M. Quesada (May. 08). Physical criterium (lamdas) PREVAILS over approximation (critical exciton number) 303 if(P1<=(P2+P3)) go_ahead=false; 304 221 305 if (go_ahead && aFragment.GetA() > 4) 222 { 223 // if (aFragment.GetNumberOfParticles() < 1) { 224 // aFragment.SetNumberOfHoles(aFragment.GetNumberOfHoles()+1); 225 // aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()+1); 226 // } 227 306 { 228 307 G4double TotalEmissionProbability = aEmission.GetTotalProbability(aFragment); 229 308 // 309 // G4cout<<"TotalEmissionProbability="<<TotalEmissionProbability<<G4endl; 310 // 230 311 // Check if number of excitons is greater than 0 231 312 // else perform equilibrium emission … … 242 323 243 324 // G4PreCompoundTransitions aTransition(aFragment); 244 325 326 //J.M.Quesada (May 08) this has already been done in order to decide what to do (preeq-eq) 245 327 // Sum of transition probabilities 246 G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);247 328 // G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment); 329 248 330 // Sum of all probabilities 249 331 G4double TotalProbability = TotalEmissionProbability + TotalTransitionProbability; 250 332 251 333 // Select subprocess 252 334 if (G4UniformRand() > TotalEmissionProbability/TotalProbability) 253 335 { 254 336 // It will be transition to state with a new number of excitons 255 ThereIsTransition = true; 256 337 ThereIsTransition = true; 257 338 // Perform the transition 258 339 aFragment = aTransition->PerformTransition(aFragment); … … 262 343 // It will be fragment emission 263 344 ThereIsTransition = false; 264 265 // Perform the emission and Add emitted fragment to Result266 345 Result->push_back(aEmission.PerformEmission(aFragment)); 267 346 } … … 288 367 { 289 368 G4ReactionProductVector * theEquilibriumResult; 369 290 370 theEquilibriumResult = GetExcitationHandler()->BreakItUp(aFragment); 291 371 … … 353 433 354 434 #endif 435 436 437 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNucleon.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // by V. Lara 26 // 27 // $Id: G4PreCompoundNucleon.cc,v 1.13 2009/02/11 18:06:00 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 // 30 // ------------------------------------------------------------------- 31 // 32 // GEANT4 Class file 33 // 34 // 35 // File name: G4PreCompoundNucleon 36 // 37 // Author: V.Lara 38 // 39 // Modified: 40 // 10.02.2009 J. M. Quesada cleanup 41 // 27 42 28 43 #include "G4PreCompoundNucleon.hh" 29 44 #include "G4PreCompoundParameters.hh" 30 45 46 G4bool G4PreCompoundNucleon::IsItPossible(const G4Fragment& aFragment) 47 { 48 G4int pplus = aFragment.GetNumberOfCharged(); 49 G4int pneut = aFragment.GetNumberOfParticles()-pplus; 50 return (pneut >= (GetA()-GetZ()) && pplus >= GetZ()); 51 } 31 52 32 53 G4double G4PreCompoundNucleon:: … … 36 57 if ( !IsItPossible(aFragment) ) return 0.0; 37 58 38 const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();39 40 59 G4double U = aFragment.GetExcitationEnergy(); 41 60 G4double P = aFragment.GetNumberOfParticles(); … … 43 62 G4double N = P + H; 44 63 45 // g = (6.0/pi2)*a*A46 // G4EvaporationLevelDensityParameter theLDP;47 64 G4double g0 = (6.0/pi2)*aFragment.GetA() * 48 65 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 49 // theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);66 50 67 G4double g1 = (6.0/pi2)*GetRestA() * 51 68 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 52 // theLDP.LevelDensityParameter(static_cast<G4int>(GetRestA()),static_cast<G4int>(GetRestZ()),U);53 69 54 70 G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0; 55 G4double A1 = (A0*g0 - P/2.0)/g1; 71 72 G4double A1 = (A0 - P/2.0)/g1; 56 73 57 74 G4double E0 = std::max(0.0,U - A0); … … 61 78 62 79 63 G4double Probability = 2.0/(hbarc*hbarc*hbarc) * GetReducedMass() * 64 GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) * r0 * r0 * std::pow(GetRestA(),2.0/3.0) * GetAlpha() * (eKin + GetBeta()) * 65 P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 * 66 g1/(g0*g0); 80 G4double Probability = 1.0/pi2*2.0/(hbarc*hbarc*hbarc) * GetReducedMass() 81 * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) 82 * eKin*CrossSection(eKin)*millibarn * P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 * g1/(g0*g0); 83 84 if (GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())<0.0 85 || CrossSection(eKin) <0) { 86 std::ostringstream errOs; 87 G4cout<<"WARNING: NEGATIVE VALUES "<<G4endl; 88 errOs << "Rj=" << GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) 89 <<G4endl; 90 errOs <<" xsec("<<eKin<<" MeV) ="<<CrossSection(eKin)<<G4endl; 91 errOs <<" A="<<GetA()<<" Z="<<GetZ()<<G4endl; 92 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 93 } 67 94 68 95 return Probability; 69 96 } 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundParameters.cc
r819 r962 26 26 // 27 27 // $Id: G4PreCompoundParameters.cc,v 1.3 2006/06/29 20:59:29 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTransitions.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4PreCompoundTransitions.cc,v 1.13 2007/07/23 12:48:54 ahoward Exp $ 28 // GEANT4 tag $Name: $ 29 // 30 // by V. Lara 26 // $Id: G4PreCompoundTransitions.cc,v 1.20 2009/02/10 16:01:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 28 // 29 // ------------------------------------------------------------------- 30 // 31 // GEANT4 Class file 32 // 33 // 34 // File name: G4PreCompoundIon 35 // 36 // Author: V.Lara 37 // 38 // Modified: 39 // 16.02.2008 J. M. Quesada fixed bugs 40 // 06.09.2008 J. M. Quesada added external choices for: 41 // - "never go back" hipothesis (useNGB=true) 42 // - CEM transition probabilities (useCEMtr=true) 31 43 32 44 #include "G4PreCompoundTransitions.hh" … … 55 67 CalculateProbability(const G4Fragment & aFragment) 56 68 { 69 //G4cout<<"In G4PreCompoundTransitions.cc useNGB="<<useNGB<<G4endl; 70 //G4cout<<"In G4PreCompoundTransitions.cc useCEMtr="<<useCEMtr<<G4endl; 71 57 72 // Fermi energy 58 73 const G4double FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy(); … … 74 89 G4double Z = static_cast<G4double>(aFragment.GetZ()); 75 90 G4double U = aFragment.GetExcitationEnergy(); 76 77 // Relative Energy (T_{rel}) 78 G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N; 79 80 // Sample kind of nucleon-projectile 81 G4bool ChargedNucleon(false); 82 G4double chtest = 0.5; 83 if (P > 0) chtest = aFragment.GetNumberOfCharged()/P; 84 if (G4UniformRand() < chtest) ChargedNucleon = true; 85 86 // Relative Velocity: 87 // <V_{rel}>^2 88 G4double RelativeVelocitySqr(0.0); 89 if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2; 90 else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2; 91 92 // <V_{rel}> 93 G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr); 94 95 // Proton-Proton Cross Section 96 G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn; 97 // Proton-Neutron Cross Section 98 G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn; 99 100 // Averaged Cross Section: \sigma(V_{rel}) 101 // G4double AveragedXSection = (ppXSection+npXSection)/2.0; 102 G4double AveragedXSection(0.0); 103 if (ChargedNucleon) 104 { 105 AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0); 106 } 107 else 108 { 109 AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0); 110 } 111 112 // Fermi relative energy ratio 113 G4double FermiRelRatio = FermiEnergy/RelativeEnergy; 114 115 // This factor is introduced to take into account the Pauli principle 116 G4double PauliFactor = 1.0 - (7.0/5.0)*FermiRelRatio; 117 if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0); 118 119 // Interaction volume 120 G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0); 121 122 // Transition probability for \Delta n = +2 123 // TransitionProb1 = 0.00332*AveragedXSection*PauliFactor*std::sqrt(RelativeEnergy)/ 124 // std::pow(1.2 + 1.0/(4.7*RelativeVelocity), 3.0); 125 TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint; 126 if (TransitionProb1 < 0.0) TransitionProb1 = 0.0; 127 128 // g = (6.0/pi2)*aA; 129 // G4double a = theLDP.LevelDensityParameter(A,Z,U-G4PairingCorrection::GetInstance()->GetPairingCorrection(A,Z)); 130 G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 131 // GE = g*E where E is Excitation Energy 132 G4double GE = (6.0/pi2)*a*A*U; 133 134 135 // F(p,h) = 0.25*(p^2 + h^2 + p - h) - 0.5*h 136 G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0); 137 G4bool NeverGoBack(false); 138 //AH if (U-Fph < 0.0) NeverGoBack = true; 139 if (GE-Fph < 0.0) NeverGoBack = true; 140 // F(p+1,h+1) 141 G4double Fph1 = Fph + N/2.0; 142 // (n+1)/n ((g*E - F(p,h))/(g*E - F(p+1,h+1)))^(n+1) 143 G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0); 144 145 146 if (NeverGoBack) 147 { 91 92 if(U<10*eV) return 0.0; 93 94 //J. M. Quesada (Feb. 08) new physics 95 // OPT=1 Transitions are calculated according to Gudima's paper (original in G4PreCompound from VL) 96 // OPT=2 Transitions are calculated according to Gupta's formulae 97 // 98 99 100 101 if (useCEMtr){ 102 103 104 // Relative Energy (T_{rel}) 105 G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N; 106 107 // Sample kind of nucleon-projectile 108 G4bool ChargedNucleon(false); 109 G4double chtest = 0.5; 110 if (P > 0) chtest = aFragment.GetNumberOfCharged()/P; 111 if (G4UniformRand() < chtest) ChargedNucleon = true; 112 113 // Relative Velocity: 114 // <V_{rel}>^2 115 G4double RelativeVelocitySqr(0.0); 116 if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2; 117 else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2; 118 119 // <V_{rel}> 120 G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr); 121 122 // Proton-Proton Cross Section 123 G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn; 124 // Proton-Neutron Cross Section 125 G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn; 126 127 // Averaged Cross Section: \sigma(V_{rel}) 128 // G4double AveragedXSection = (ppXSection+npXSection)/2.0; 129 G4double AveragedXSection(0.0); 130 if (ChargedNucleon) 131 { 132 //JMQ: small bug fixed 133 // AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0); 134 AveragedXSection = ((Z-1.0) * ppXSection + (A-Z) * npXSection) / (A-1.0); 135 } 136 else 137 { 138 AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0); 139 } 140 141 // Fermi relative energy ratio 142 G4double FermiRelRatio = FermiEnergy/RelativeEnergy; 143 144 // This factor is introduced to take into account the Pauli principle 145 G4double PauliFactor = 1.0 - (7.0/5.0)*FermiRelRatio; 146 if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0); 147 148 // Interaction volume 149 // G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0); 150 G4double xx=2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity); 151 G4double Vint = (4.0/3.0)*pi*xx*xx*xx; 152 153 // Transition probability for \Delta n = +2 154 155 TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint; 156 if (TransitionProb1 < 0.0) TransitionProb1 = 0.0; 157 158 G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 159 // GE = g*E where E is Excitation Energy 160 G4double GE = (6.0/pi2)*a*A*U; 161 162 G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0); 163 164 //G4bool NeverGoBack(false); 165 G4bool NeverGoBack; 166 if(useNGB) NeverGoBack=true; 167 else NeverGoBack=false; 168 169 170 //JMQ/AH bug fixed: if (U-Fph < 0.0) NeverGoBack = true; 171 if (GE-Fph < 0.0) NeverGoBack = true; 172 173 // F(p+1,h+1) 174 G4double Fph1 = Fph + N/2.0; 175 176 G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0); 177 178 179 if (NeverGoBack) 180 { 148 181 TransitionProb2 = 0.0; 149 182 TransitionProb3 = 0.0; 150 } 151 else 152 { 153 // Transition probability for \Delta n = -2 (at F(p,h) = 0) 154 // TransitionProb2 = max(0, (TransitionProb1*P*H*(P+H+1.0)*(P+H-2.0))/(GE*GE)); 155 // TransitionProb2 = (TransitionProb1*P*H*(P+H+1.0)*(P+H-2.0))/(GE*GE); 156 TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph)); 183 } 184 else 185 { 186 // Transition probability for \Delta n = -2 (at F(p,h) = 0) 187 TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph)); 188 if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 189 190 // Transition probability for \Delta n = 0 (at F(p,h) = 0) 191 TransitionProb3 = TransitionProb1* ((N+1.0)/N) * ProbFactor * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph); 192 if (TransitionProb3 < 0.0) TransitionProb3 = 0.0; 193 } 194 195 // G4cout<<"U = "<<U<<G4endl; 196 // G4cout<<"N="<<N<<" P="<<P<<" H="<<H<<G4endl; 197 // G4cout<<"l+ ="<<TransitionProb1<<" l- ="<< TransitionProb2<<" l0 ="<< TransitionProb3<<G4endl; 198 return TransitionProb1 + TransitionProb2 + TransitionProb3; 199 } 200 201 else { 202 //JMQ: Transition probabilities from Gupta's work 203 204 G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 205 // GE = g*E where E is Excitation Energy 206 G4double GE = (6.0/pi2)*a*A*U; 207 208 G4double Kmfp=2.; 209 210 211 TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U); 212 if (TransitionProb1 < 0.0) TransitionProb1 = 0.0; 213 214 if (useNGB){ 215 TransitionProb2=0.; 216 TransitionProb3=0.; 217 } 218 else{ 219 if (N<=1) TransitionProb2=0. ; 220 else TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U); 157 221 if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 158 159 160 // Transition probability for \Delta n = 0 (at F(p,h) = 0) 161 // TransitionProb3 = TransitionProb1*(P+H+1.0)*(P*(P-1.0)+4.0*P*H+H*(H-1.0))/((P+H)*GE); 162 TransitionProb3 = TransitionProb1 * ProbFactor * ((N+1.0)/N) * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph); 163 if (TransitionProb3 < 0.0) TransitionProb3 = 0.0; 164 } 165 166 return TransitionProb1 + TransitionProb2 + TransitionProb3; 222 TransitionProb3=0.; 223 } 224 225 // G4cout<<"U = "<<U<<G4endl; 226 // G4cout<<"N="<<N<<" P="<<P<<" H="<<H<<G4endl; 227 // G4cout<<"l+ ="<<TransitionProb1<<" l- ="<< TransitionProb2<<" l0 ="<< TransitionProb3<<G4endl; 228 return TransitionProb1 + TransitionProb2 + TransitionProb3; 229 } 167 230 } 168 231 … … 185 248 } 186 249 187 // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and in proportion to the number charges w.r.t. number of particles 188 if(deltaN < 0 && G4UniformRand() <= static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) && (result.GetNumberOfCharged() >= 1)) 250 // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and in proportion 251 // to the number charges w.r.t. number of particles, PROVIDED that there are charged particles 252 if(deltaN < 0 && G4UniformRand() <= 253 static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) 254 && (result.GetNumberOfCharged() >= 1)) { 189 255 result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative! 190 191 // result.SetNumberOfParticles was here 192 // result.SetNumberOfHoles was here193 // the following lines have to be before SetNumberOfCharged, otherwise the check onnumber of charged vs. number of particles fails256 } 257 258 // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on 259 // number of charged vs. number of particles fails 194 260 result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2); 195 261 result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2); 196 262 197 // With weight Z/A, number of charged particles is decreased on +1 198 // if ((deltaN > 0 || result.GetNumberOfCharged() > 0) && // AH/JMQ check is now in initialize within G4VPreCompoundFragment 263 // With weight Z/A, number of charged particles is increased with +1 199 264 if ( ( deltaN > 0 ) && 200 265 (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/ 201 266 std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.))) 202 267 { 203 268 result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); … … 210 275 } 211 276 212 // moved from above to make code more readable213 // result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);214 // result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2);215 216 277 return result; 217 278 } -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundFragment.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VPreCompoundFragment.cc,v 1.12 2009/02/10 16:01:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 26 28 // 27 // $Id: G4VPreCompoundFragment.cc,v 1.4 2007/07/23 09:56:40 ahoward Exp $28 // GEANT4 tag $Name: $29 // J. M. Quesada (August 2008). 30 // Based on previous work by V. Lara 29 31 // 30 // by V. Lara31 32 32 33 #include "G4VPreCompoundFragment.hh" … … 142 143 if ((theRestNucleusA < theRestNucleusZ) || 143 144 (theRestNucleusA < theA) || 144 (theRestNucleusZ < theZ) || 145 (aFragment.GetNumberOfCharged() < theZ)) // AH last argument from JMQ 145 (theRestNucleusZ < theZ)) 146 146 { 147 147 // In order to be sure that emission probability will be 0. … … 155 155 GetCoulombBarrier(static_cast<G4int>(theRestNucleusA),static_cast<G4int>(theRestNucleusZ), 156 156 aFragment.GetExcitationEnergy()); 157 157 158 158 159 // Compute Binding Energies for fragments 159 160 // (needed to separate a fragment from the nucleus) … … 164 165 165 166 // Compute Maximal Kinetic Energy which can be carried by fragments after separation 167 // This is the true (assimptotic) maximal kinetic energy 166 168 G4double m = aFragment.GetMomentum().m(); 167 169 G4double rm = GetRestNuclearMass(); 168 170 G4double em = GetNuclearMass(); 169 171 theMaximalKineticEnergy = ((m - rm)*(m + rm) + em*em)/(2.0*m) - em; 170 // - theCoulombBarrier;172 171 173 172 174 return; -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundIon.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4VPreCompoundIon.cc,v 1. 5 2007/07/23 09:56:40ahoward Exp $28 // GEANT4 tag $Name: $27 // $Id: G4VPreCompoundIon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // by V. Lara 31 32 #include "G4VPreCompoundIon.hh" 33 #include "G4PreCompoundParameters.hh" 34 //#include "G4EvaporationLevelDensityParameter.hh" 35 36 37 G4double G4VPreCompoundIon:: 38 ProbabilityDistributionFunction(const G4double eKin, 39 const G4Fragment& aFragment) 40 { 41 if ( !IsItPossible(aFragment) ) return 0.0; 42 43 const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0(); 44 45 G4double U = aFragment.GetExcitationEnergy(); 46 G4double P = aFragment.GetNumberOfParticles(); 47 G4double H = aFragment.GetNumberOfHoles(); 48 G4double N = P + H; 49 50 // g = (6.0/pi2)*a*A 51 // G4EvaporationLevelDensityParameter theLDP; 52 G4double g0 = (6.0/pi2)*aFragment.GetA() * 53 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 54 // theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U); 55 G4double g1 = (6.0/pi2)*GetRestA() * 56 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 57 // theLDP.LevelDensityParameter(GetRestA(),GetRestZ(),U); 58 G4double gj = (6.0/pi2)*GetA() * 59 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 60 // theLDP.LevelDensityParameter(GetA(),GetZ(),U); 61 62 G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0; 63 G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1); 64 G4double Aj = GetA()*(GetA()+1.0)/4.0; 65 66 G4double E0 = std::max(0.0,U - A0); 67 if (E0 == 0.0) return 0.0; 68 G4double E1 = std::max(0.0,GetMaximalKineticEnergy() + GetCoulombBarrier() - eKin - A1); 69 G4double Ej = std::max(0.0,U + GetBindingEnergy() - Aj); 70 71 G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*(eKin+GetBindingEnergy()))))* 72 GetAlpha()*(eKin+GetBeta())/(r0*std::pow(GetRestA(),1.0/3.0)) * 73 CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P); 74 75 G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0); 76 77 G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0; 78 79 G4double Probability = pA * pB * pC; 80 81 return Probability; 82 } 31 // 32 //J.M. Quesada (Apr. 2008) DUMMY abstract base class for ions 33 // it is not ihherited by anything 34 // Suppressed 83 35 84 36 … … 88 40 89 41 42 43 44 45 46 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundNucleon.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4VPreCompoundNucleon.cc,v 1.5 2007/07/23 09:56:40 ahoward Exp $ 28 // GEANT4 tag $Name: $ 27 28 // $Id: G4VPreCompoundNucleon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $ 29 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 30 29 31 // 30 32 // by V. Lara 31 32 #include "G4VPreCompoundNucleon.hh" 33 #include "G4PreCompoundParameters.hh" 34 // #include "G4EvaporationLevelDensityParameter.hh"33 // 34 //J.M. Quesada (Apr. 2008) DUMMY abstract base class for nucleons. 35 // it is not ihherited by anything 36 // Suppressed 35 37 36 38 37 G4double G4VPreCompoundNucleon::38 ProbabilityDistributionFunction(const G4double eKin,39 const G4Fragment& aFragment)40 {41 if ( !IsItPossible(aFragment) ) return 0.0;42 43 const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();44 45 G4double U = aFragment.GetExcitationEnergy();46 G4double P = aFragment.GetNumberOfParticles();47 G4double H = aFragment.GetNumberOfHoles();48 G4double N = P + H;49 50 // g = (6.0/pi2)*a*A51 // G4EvaporationLevelDensityParameter theLDP;52 G4double g0 = (6.0/pi2)*aFragment.GetA() *53 G4PreCompoundParameters::GetAddress()->GetLevelDensity();54 // theLDP.LevelDensityParameter(G4int(aFragment.GetA()),G4int(aFragment.GetZ()),U);55 G4double g1 = (6.0/pi2)*GetRestA() *56 G4PreCompoundParameters::GetAddress()->GetLevelDensity();57 // theLDP.LevelDensityParameter(G4int(GetRestA()),G4int(GetRestZ()),U);58 59 60 G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;61 G4double A1 = (A0*g0 - P/2.0)/g1;62 63 G4double E0 = std::max(0.0,U - A0);64 if (E0 == 0.0) return 0.0;65 G4double E1 = std::max(0.0,U - eKin - GetBindingEnergy() - A1);66 67 G4double Probability = 2.0/(hbarc*hbarc*hbarc) * GetReducedMass() *68 GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) * r0 * r0 * std::pow(GetRestA(),2.0/3.0) * GetAlpha() * (eKin + GetBeta()) *69 P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 *70 g1/(g0*g0);71 72 return Probability;73 }74 75
Note: See TracChangeset
for help on using the changeset viewer.
