Changeset 962 for trunk/source/processes/hadronic/models/pre_equilibrium
- Timestamp:
- Apr 6, 2009, 12:30:29 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/pre_equilibrium
- Files:
-
- 35 edited
-
History (modified) (1 diff)
-
exciton_model/include/G4GNASHTransitions.hh (modified) (1 diff)
-
exciton_model/include/G4LowEIonFragmentation.hh (modified) (1 diff)
-
exciton_model/include/G4PreCompoundAlpha.hh (modified) (4 diffs)
-
exciton_model/include/G4PreCompoundDeuteron.hh (modified) (3 diffs)
-
exciton_model/include/G4PreCompoundEmission.hh (modified) (2 diffs)
-
exciton_model/include/G4PreCompoundEmission.icc (modified) (3 diffs)
-
exciton_model/include/G4PreCompoundFragment.hh (modified) (3 diffs)
-
exciton_model/include/G4PreCompoundFragmentVector.hh (modified) (2 diffs)
-
exciton_model/include/G4PreCompoundFragmentVector.icc (modified) (2 diffs)
-
exciton_model/include/G4PreCompoundHe3.hh (modified) (3 diffs)
-
exciton_model/include/G4PreCompoundIon.hh (modified) (2 diffs)
-
exciton_model/include/G4PreCompoundModel.hh (modified) (6 diffs)
-
exciton_model/include/G4PreCompoundNeutron.hh (modified) (4 diffs)
-
exciton_model/include/G4PreCompoundNucleon.hh (modified) (3 diffs)
-
exciton_model/include/G4PreCompoundParameters.hh (modified) (2 diffs)
-
exciton_model/include/G4PreCompoundProton.hh (modified) (3 diffs)
-
exciton_model/include/G4PreCompoundTransitions.hh (modified) (4 diffs)
-
exciton_model/include/G4PreCompoundTriton.hh (modified) (3 diffs)
-
exciton_model/include/G4VPreCompoundFragment.hh (modified) (6 diffs)
-
exciton_model/include/G4VPreCompoundFragment.icc (modified) (2 diffs)
-
exciton_model/include/G4VPreCompoundIon.hh (modified) (2 diffs)
-
exciton_model/include/G4VPreCompoundNucleon.hh (modified) (1 diff)
-
exciton_model/include/G4VPreCompoundTransitions.hh (modified) (2 diffs)
-
exciton_model/src/G4PreCompoundEmission.cc (modified) (7 diffs)
-
exciton_model/src/G4PreCompoundFragment.cc (modified) (5 diffs)
-
exciton_model/src/G4PreCompoundFragmentVector.cc (modified) (3 diffs)
-
exciton_model/src/G4PreCompoundIon.cc (modified) (2 diffs)
-
exciton_model/src/G4PreCompoundModel.cc (modified) (9 diffs)
-
exciton_model/src/G4PreCompoundNucleon.cc (modified) (4 diffs)
-
exciton_model/src/G4PreCompoundParameters.cc (modified) (1 diff)
-
exciton_model/src/G4PreCompoundTransitions.cc (modified) (5 diffs)
-
exciton_model/src/G4VPreCompoundFragment.cc (modified) (4 diffs)
-
exciton_model/src/G4VPreCompoundIon.cc (modified) (2 diffs)
-
exciton_model/src/G4VPreCompoundNucleon.cc (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/pre_equilibrium/History
r819 r962 14 14 * Please list in reverse chronological order (last date on top) 15 15 --------------------------------------------------------------- 16 17 13-February 2009 V.Ivanchenko hadr-pre-V09-02-03 18 --------------------------------------------------- 19 G4PreCompoundXXX - changed the shape of probabilities (return back to 9.2) 20 for d, t, He3, alpha (JMQ) 21 22 12-February 2009 V.Ivanchenko hadr-pre-V09-02-02 23 --------------------------------------------------- 24 G4PreCompoundXXX - changed the shape of probabilities 25 for d, t, He3, alpha near the Coulomb barrier (JMQ) 26 27 11-February 2009 V.Ivanchenko hadr-pre-V09-02-01 28 --------------------------------------------------- 29 G4PreCompoundXXX - set default Opt3 back, add decrease Coulomb barrier 30 for d, t, a, he3 (JMQ) 31 32 10-February 2009 V.Ivanchenko hadr-pre-V09-02-00 33 --------------------------------------------------- 34 Some clean up of comments 35 G4PreCompoundIon - fixed probability of light ion emmision (JMQ) 36 G4PreCompoundXXX - by default Opt1 is used for d, t, a, he3, 37 Opt3 for n, p (JMQ) 38 39 09-December 2008 A.Howard hadr-pre-V09-01-15 40 --------------------------------------------------- 41 Added protection for close to zero excitation energy in 42 G4PreCompoundModel.cc (according to JMQ, MAC) to not try to de-excite. 43 44 09-December 2008 A.Howard hadr-pre-V09-01-14 45 --------------------------------------------------- 46 Added protection for close to zero excitation energy in 47 G4PreCompoundTransitions.cc (returns 0.0), prevents FPE later on. 48 49 27-November 2008 A.Howard hadr-pre-V09-01-13 50 --------------------------------------------------- 51 Added data member initialisation to G4VPreCompoundTransitions.hh 52 53 19-November 2008 A.Howard hadr-pre-V09-01-12 54 --------------------------------------------------- 55 JMQ fix to G4PreCompoundNeutron.cc for Zirconium (if (nu < 0.)nu=-nu). 56 57 30-September 2008 A.Howard hadr-pre-V09-01-11 58 --------------------------------------------------- 59 Trivial protection against negative probabilities for incident protons on 60 targets with A < Carbon. 61 62 22-September 2008 A.Howard hadr-pre-V09-01-10 63 --------------------------------------------------- 64 JMQ's latest developments - which are an extension of hadr-pre-V09-01-08 with 65 cross-section options (equivalent to the de-excitation modifications) and other 66 options (SICB, Never Go Back,...). In addition Gunter's fix to the factorials 67 is included from hadr-pre-V09-01-09. 68 Added new files: G4PreCompoundAlpha.cc G4PreCompoundDeuteron.cc 69 G4PreCompoundHe3.cc G4PreCompoundNeutron.cc G4PreCompoundProton.cc 70 G4PreCompoundTriton.cc (Classes previously existed. Now they have more 71 involved cross-section calculations). 72 73 19-September 2008 A.Howard hadr-pre-V09-01-09 74 --------------------------------------------------- 75 Including Gunter's fix (see 11-August below) on top of ref-07 76 (hadr-pre-V09-00-04), i.e. without JMQ's developments. 77 78 11-August 2008 G.Folger 79 ----------------------------------------------- 80 Rewrite algorithm in G4PreCompoundEmission::rho() to avoid frequent 81 floating point overflow when using Precompund in combination 82 with Binary Cascade. 83 84 85 24 July 2008 J. M. Quesada hadr-pre-V09-01-08 86 --------------------------------------------------- 87 Minor fixes and Coulomb barrier just for Wellisch's proton cross section 88 (OPT=2) 89 90 91 23 July 2008 V.Ivanchenko hadr-pre-V09-01-07 92 --------------------------------------------------- 93 Fixed mistake in tagging 94 95 23 July 2008 V.Ivanchenko hadr-pre-V09-01-06 96 --------------------------------------------------- 97 Return back Coulomb barrier initialisation (J.M.Quesada) 98 99 27 June 2008 V.Ivanchenko hadr-pre-V09-01-05 100 --------------------------------------------------- 101 Fixed G4PreCompoundNucleon and G4PreCompoundIon for 102 Coulomb barrier, added protections for zero cross sections, 103 OPT=2 is used (J.M.Quesada) 104 105 05 June 2008 J. M. Quesada hadr-pre-V09-01-04 106 --------------------------------------------------- 107 Bug fixed in OPT=1 (Chatterjee) charged particle cross sections. 108 Unphysical values at very low emission energies have been corrected (set 109 to 0); OPT=1 ( OPT=2 in previous tag) 110 111 112 15-May 2008 J. M.Quesada hadr-pre-V09-01-03 113 ----------------------------------------------- 114 The retrieval of transition probabilities in G4PreCompoundModel.cc 115 (for protection against unphysically crossing the landa_+=landa_- condition) 116 was misplaced. Now it has been placed inside the preequilibrium transitions 117 loop. 118 119 8-May 2008 J. M.Quesada hadr-pre-V09-01-02 120 ----------------------------------------------- 121 Protection against non physical situation has been set: Equilibrium exciton number 122 Neq (when reached, equilibrium regime starts) should correspond to equal transition 123 probabilities "back" and "forth". Nevertheless, for heavy target (Neq is big even 124 for low incident energies) after first emission this condition (equal trans. prob.) 125 is reached far before reaching corresponding Neq of the residual. Unless this "jump" 126 to equilibrium is forced in this case, preequilibrium will be 127 spuriously prolonged (as a side effect, with huge CPU consumption). 128 G4PreCompoundModel and G4PreCompoundTransitions classes have been modified. 129 130 131 1-May 2008 J. M.Quesada hadr-pre-V09-01-01 132 ----------------------------------------------- 133 - First trial with cvs. 00 skipped by mistake. 134 New inverse cross sections: 135 OPT=1 Chetterjee's parameterization to reaction cross sections from optical potential global fittings. 136 OPT=2 as OPT=1 but for protons the Wellisch's parameterization for protons is used (DEFAULT) 137 OPT=3 Kalbach's modifications of Chatterjee's parameterization of cross sections 138 OPT=4 as OPT=3 but for protons the Wellisch's parameterization for protons is used 139 140 Coulomb barrier has been suppressed as it enters through inverse cross sections. Also methods related to former Dostrovski's cross sections (alpha , beta and C parameters) have been suppressed. In all particle header files: 141 G4PreCompoundNucleon.hh ,G4PreCompoundNeutron.hh,G4PreCompoundProton.hh, G4PreCompoundIon.hh, G4PreCompoundDeuteron.hh, G4PreCompoundTriton.hh ,G4PreCompoundHe3.hh , G4PreCompoundAlpha.hh 142 143 144 Soft cutoff transition from equilibrium to equilibrium has been suppressed (unphysical and not necessary anymore).In G4PreCompoundModel. 145 146 New transition probabilities : 147 OPT=1 Gudima's formulae 148 OPT=2 Machner parameterization of matrix elements for transitions (DEFAULT) 149 150 Level density is set to A/10 (Gudima's prescription). 151 152 Several bugs fixed (average xs in nucleon-nucleon xs in G4PreCompundTransitions.cc and emission probability formula in G4PreCompoundIon.cc ) . 16 153 17 154 25-October 2007 A.Howard hadr-pre-V09-00-04 -
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.
