- Timestamp:
- Apr 20, 2009, 5:54:05 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model
- Files:
-
- 27 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4LowEIonFragmentation.hh
r962 r1007 26 26 // 27 27 // $Id: G4LowEIonFragmentation.hh,v 1.3 2006/06/29 20:58:04 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-02 -ref-02$28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // by H.P. Wellisch -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.hh
r962 r1007 25 25 // 26 26 // $Id: G4PreCompoundEmission.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $ 27 // GEANT4 tag $Name: geant4-09-02 -ref-02$27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // Hadronic Process: Nuclear Preequilibrium -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.icc
r962 r1007 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.Lara31 //32 // Modified:33 26 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 34 27 // cross section option … … 41 34 return; 42 35 } 36 43 37 44 38 inline G4double G4PreCompoundEmission::GetTotalProbability(const G4Fragment & aFragment) -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.hh
r962 r1007 25 25 // 26 26 // 27 // $Id: G4PreCompoundFragmentVector.hh,v 1. 6 2009/02/10 16:01:37 vnivanchExp $28 // GEANT4 tag $Name: geant4-09-02 -ref-02$27 // $Id: G4PreCompoundFragmentVector.hh,v 1.4 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear Preequilibrium … … 39 39 #define G4PreCompoundFragmentVector_h 1 40 40 41 41 42 #include "G4VPreCompoundFragment.hh" 43 44 42 45 43 46 class G4PreCompoundFragmentVector -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundModel.hh
r962 r1007 26 26 // 27 27 // $Id: G4PreCompoundModel.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02 -ref-02$28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundParameters.hh
r962 r1007 26 26 // 27 27 // $Id: G4PreCompoundParameters.hh,v 1.5 2008/05/08 10:34:25 quesada Exp $ 28 // GEANT4 tag $Name: geant4-09-02 -ref-02$28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTransitions.hh
r962 r1007 26 26 // 27 27 // $Id: G4PreCompoundTransitions.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02 -ref-02$28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.hh
r962 r1007 24 24 // ******************************************************************** 25 25 // 26 // 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 // 30 // J. M. Quesada (August 2008). 31 // Based on previous work by V. Lara 26 //J. M. Quesada (August 2008). 27 //Based on previous work by V. Lara 32 28 // 33 29 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse … … 35 31 // JMQ (06 September 2008) Also external choice has been added for: 36 32 // - superimposed Coulomb barrier (if useSICB=true) 33 37 34 38 35 #ifndef G4VPreCompoundFragment_h -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundIon.hh
r962 r1007 27 27 28 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 // GEANT4 tag $Name: geant4-09-02 $ 30 30 // 31 31 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundNucleon.hh
r962 r1007 26 26 // 27 27 // $Id: G4VPreCompoundNucleon.hh,v 1.8 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02 -ref-02$28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundAlpha.cc
r968 r1007 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4PreCompoundAlpha.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 // 30 // ------------------------------------------------------------------- 31 // 32 // GEANT4 Class file 33 // 34 // 35 // File name: G4PreCompoundAlpha 36 // 37 // Author: V.Lara 38 // 39 // Modified: 40 // 21.08.2008 J. M. Quesada add choice of options 41 // 10.02.2009 J. M. Quesada set default opt1 42 // 26 //J.M.Quesada (August 08). New source file 27 // 28 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse 29 // cross section option 43 30 44 31 #include "G4PreCompoundAlpha.hh" 45 32 46 G4ReactionProduct * G4PreCompoundAlpha::GetReactionProduct() const47 {48 G4ReactionProduct * theReactionProduct =49 new G4ReactionProduct(G4Alpha::AlphaDefinition());50 theReactionProduct->SetMomentum(GetMomentum().vect());51 theReactionProduct->SetTotalEnergy(GetMomentum().e());33 G4ReactionProduct * G4PreCompoundAlpha::GetReactionProduct() const 34 { 35 G4ReactionProduct * theReactionProduct = 36 new G4ReactionProduct(G4Alpha::AlphaDefinition()); 37 theReactionProduct->SetMomentum(GetMomentum().vect()); 38 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 52 39 #ifdef PRECOMPOUND_TEST 53 theReactionProduct->SetCreatorModel("G4PrecompoundModel");40 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 54 41 #endif 55 return theReactionProduct; 56 } 57 58 G4double G4PreCompoundAlpha::FactorialFactor(const G4double N, const G4double P) 59 { 60 return 42 return theReactionProduct; 43 } 44 45 46 G4double G4PreCompoundAlpha::FactorialFactor(const G4double N, const G4double P) 47 { 48 return 61 49 (N-4.0)*(P-3.0)*( 62 50 (((N-3.0)*(P-2.0))/2.0) *( … … 66 54 ) 67 55 ); 68 } 69 70 G4double G4PreCompoundAlpha::CoalescenceFactor(const G4double A) 71 { 72 return 4096.0/(A*A*A); 73 } 74 75 G4double G4PreCompoundAlpha::GetRj(const G4int NumberParticles, const G4int NumberCharged) 76 { 77 G4double rj = 0.0; 78 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2)*(NumberParticles-3); 79 if(NumberCharged >=2 && (NumberParticles-NumberCharged) >=2 ) { 80 rj = 6.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged)* 81 (NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); 82 } 83 return rj; 84 } 56 } 57 58 G4double G4PreCompoundAlpha::CoalescenceFactor(const G4double A) 59 { 60 return 4096.0/(A*A*A); 61 } 62 63 64 65 G4double G4PreCompoundAlpha::GetRj(const G4int NumberParticles, const G4int NumberCharged) 66 { 67 G4double rj = 0.0; 68 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2)*(NumberParticles-3); 69 if(NumberCharged >=2 && (NumberParticles-NumberCharged) >=2 ) rj = 6.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); 70 return rj; 71 } 72 73 74 85 75 86 76 //////////////////////////////////////////////////////////////////////////////////// … … 90 80 //OPT=3,4 Kalbach's parameterization 91 81 // 92 G4double G4PreCompoundAlpha::CrossSection(const G4double K)82 G4double G4PreCompoundAlpha::CrossSection(const G4double K) 93 83 { 94 84 … … 125 115 //---------------- 126 116 // 127 G4double G4PreCompoundAlpha::GetAlpha()128 {129 130 G4double aZ = GetZ() + GetRestZ();131 if (aZ <= 30)132 {133 134 }135 else if (aZ <= 50)136 {137 138 }139 else if (aZ < 70)140 {141 142 }143 else144 {145 146 }147 return 1.0+C;148 }117 G4double G4PreCompoundAlpha::GetAlpha() 118 { 119 G4double C = 0.0; 120 G4double aZ = GetZ() + GetRestZ(); 121 if (aZ <= 30) 122 { 123 C = 0.10; 124 } 125 else if (aZ <= 50) 126 { 127 C = 0.1 + -((aZ-50.)/20.)*0.02; 128 } 129 else if (aZ < 70) 130 { 131 C = 0.08 + -((aZ-70.)/20.)*0.02; 132 } 133 else 134 { 135 C = 0.06; 136 } 137 return 1.0+C; 138 } 149 139 // 150 140 //-------------------- 151 141 // 152 G4double G4PreCompoundAlpha::GetBeta()153 {154 return -GetCoulombBarrier();155 }142 G4double G4PreCompoundAlpha::GetBeta() 143 { 144 return -GetCoulombBarrier(); 145 } 156 146 // 157 147 //********************* OPT=1,2 : Chatterjee's cross section ************************ … … 163 153 G4double Kc=K; 164 154 165 // JMQ xsec is set consta nt above limit of validity155 // JMQ xsec is set constat above limit of validity 166 156 if (K>50) Kc=50; 167 157 … … 223 213 G4double ra=1.20; 224 214 225 //JMQ 13/02/09 increase of reduced radius to lower the barrier 226 // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 227 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra); 215 ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 228 216 ecsq = ec * ec; 229 217 p = p0 + p1/ec + p2/ecsq; … … 248 236 249 237 if (elab <= ec) { //start for E<Ec 250 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 238 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 251 239 } //end for E<Ec 252 240 else { //start for E>Ec … … 264 252 265 253 // ************************** end of cross sections ******************************* 254 255 256 257 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundDeuteron.cc
r968 r1007 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PreCompoundDeuteron.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 28 // 29 // ------------------------------------------------------------------- 30 // 31 // GEANT4 Class file 32 // 33 // 34 // File name: G4PreCompoundDeuteron 35 // 36 // Author: V.Lara 37 // 38 // Modified: 39 // 21.08.2008 J. M. Quesada add choice of options 40 // 10.02.2009 J. M. Quesada set default opt1 41 // 26 //J.M.Quesada (August 08). New source file 27 // 28 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse 29 // cross section option 42 30 43 31 #include "G4PreCompoundDeuteron.hh" 44 32 45 G4ReactionProduct * G4PreCompoundDeuteron::GetReactionProduct() const 46 { 47 G4ReactionProduct * theReactionProduct = 48 new G4ReactionProduct(G4Deuteron::DeuteronDefinition()); 49 theReactionProduct->SetMomentum(GetMomentum().vect()); 50 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 33 34 35 G4ReactionProduct * G4PreCompoundDeuteron::GetReactionProduct() const 36 { 37 G4ReactionProduct * theReactionProduct = 38 new G4ReactionProduct(G4Deuteron::DeuteronDefinition()); 39 theReactionProduct->SetMomentum(GetMomentum().vect()); 40 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 51 41 #ifdef PRECOMPOUND_TEST 52 theReactionProduct->SetCreatorModel("G4PrecompoundModel");42 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 53 43 #endif 54 return theReactionProduct;55 }56 57 58 G4double G4PreCompoundDeuteron::FactorialFactor(const G4double N, const G4double P)59 {60 return (N-1.0)*(N-2.0)*(P-1.0)*P/2.0;61 }44 return theReactionProduct; 45 } 46 47 48 G4double G4PreCompoundDeuteron::FactorialFactor(const G4double N, const G4double P) 49 { 50 return (N-1.0)*(N-2.0)*(P-1.0)*P/2.0; 51 } 62 52 63 G4double G4PreCompoundDeuteron::CoalescenceFactor(const G4double A) 64 { 65 return 16.0/A; 66 } 67 68 G4double G4PreCompoundDeuteron::GetRj(const G4int NumberParticles, const G4int NumberCharged) 69 { 70 G4double rj = 0.0; 71 G4double denominator = NumberParticles*(NumberParticles-1); 72 if(NumberCharged >=1 && (NumberParticles-NumberCharged) >=1) { 73 rj = 2.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged)) 74 / static_cast<G4double>(denominator); 75 } 76 return rj; 77 } 53 G4double G4PreCompoundDeuteron::CoalescenceFactor(const G4double A) 54 { 55 return 16.0/A; 56 } 57 58 G4double G4PreCompoundDeuteron::GetRj(const G4int NumberParticles, const G4int NumberCharged) 59 { 60 G4double rj = 0.0; 61 G4double denominator = NumberParticles*(NumberParticles-1); 62 if(NumberCharged >=1 && (NumberParticles-NumberCharged) >=1) rj = 2.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged))/static_cast<G4double>(denominator); 63 64 return rj; 65 } 78 66 79 67 //////////////////////////////////////////////////////////////////////////////////// … … 83 71 //OPT=3,4 Kalbach's parameterization 84 72 // 85 G4double G4PreCompoundDeuteron::CrossSection(const G4double K)73 G4double G4PreCompoundDeuteron::CrossSection(const G4double K) 86 74 { 87 75 … … 118 106 //--------- 119 107 // 120 G4double G4PreCompoundDeuteron::GetAlpha()121 {122 123 G4double aZ = GetZ() + GetRestZ();124 if (aZ >= 70)125 {126 127 }128 else129 {130 131 }132 return 1.0 + C/2.0;133 }108 G4double G4PreCompoundDeuteron::GetAlpha() 109 { 110 G4double C = 0.0; 111 G4double aZ = GetZ() + GetRestZ(); 112 if (aZ >= 70) 113 { 114 C = 0.10; 115 } 116 else 117 { 118 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 119 } 120 return 1.0 + C/2.0; 121 } 134 122 // 135 123 //--------- 136 124 // 137 G4double G4PreCompoundDeuteron::GetBeta()138 {139 return -GetCoulombBarrier();140 }125 G4double G4PreCompoundDeuteron::GetBeta() 126 { 127 return -GetCoulombBarrier(); 128 } 141 129 // 142 130 //********************* OPT=1,2 : Chatterjee's cross section ************************ … … 152 140 153 141 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs; 154 142 //G4double Eo(0),epsilon1(0),epsilon2(0),discri(0); 155 143 156 144 … … 186 174 } 187 175 176 177 178 188 179 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 189 180 G4double G4PreCompoundDeuteron::GetOpt34(const G4double K) … … 214 205 G4double ra=0.80; 215 206 216 //JMQ 13/02/09 increase of reduced radius to lower the barrier 217 // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 218 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra); 207 ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 219 208 ecsq = ec * ec; 220 209 p = p0 + p1/ec + p2/ecsq; … … 237 226 elab = K * FragmentA / ResidualA; 238 227 sig = 0.; 239 228 240 229 if (elab <= ec) { //start for E<Ec 241 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 230 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 242 231 } //end for E<Ec 243 232 else { //start for E>Ec -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc
r962 r1007 25 25 // 26 26 // 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 // 27 // Hadronic Process: Nuclear Preequilibrium 28 // by V. Lara 29 41 30 42 31 #include "G4PreCompoundEmission.hh" … … 53 42 54 43 55 G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const44 G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const 56 45 { 57 46 return false; 58 47 } 59 48 60 G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const49 G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const 61 50 { 62 51 return true; … … 106 95 return; 107 96 } 97 98 108 99 109 100 G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment) -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragment.cc
r962 r1007 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 $ 28 // 29 // J. M. Quesada (August 2008). 30 // Based on previous work by V. Lara 26 //J. M. Quesada (August 2008). 27 //Based on previous work by V. Lara 31 28 // JMQ (06 September 2008) Also external choice has been added for: 32 29 // - superimposed Coulomb barrier (if useSICB=true) -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragmentVector.cc
r962 r1007 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PreCompoundFragmentVector.cc,v 1.11 2009/02/10 16:01:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 26 // 27 // $Id: G4PreCompoundFragmentVector.cc,v 1.7 2008/05/08 10:36:03 quesada Exp $ 28 // GEANT4 tag $Name: geant4-09-02 $ 28 29 // 29 30 // Hadronic Process: Nuclear Preequilibrium -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundHe3.cc
r968 r1007 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4PreCompoundHe3.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 // 30 // ------------------------------------------------------------------- 31 // 32 // GEANT4 Class file 33 // 34 // 35 // File name: G4PreCompoundHe3 36 // 37 // Author: V.Lara 38 // 39 // Modified: 40 // 21.08.2008 J. M. Quesada add choice of options 41 // 10.02.2009 J. M. Quesada set default opt1 42 // 26 //J.M.Quesada (August 08). New source file 27 // 28 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse 29 // cross section option 43 30 44 31 #include "G4PreCompoundHe3.hh" 45 32 46 G4ReactionProduct * G4PreCompoundHe3::GetReactionProduct() const 47 { 48 G4ReactionProduct * theReactionProduct = 49 new G4ReactionProduct(G4He3::He3Definition()); 50 theReactionProduct->SetMomentum(GetMomentum().vect()); 51 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 33 34 G4ReactionProduct * G4PreCompoundHe3::GetReactionProduct() const 35 { 36 G4ReactionProduct * theReactionProduct = 37 new G4ReactionProduct(G4He3::He3Definition()); 38 theReactionProduct->SetMomentum(GetMomentum().vect()); 39 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 52 40 #ifdef PRECOMPOUND_TEST 53 theReactionProduct->SetCreatorModel("G4PrecompoundModel");41 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 54 42 #endif 55 return theReactionProduct; 56 } 57 58 G4double G4PreCompoundHe3::FactorialFactor(const G4double N, const G4double P) 59 { 60 return 43 return theReactionProduct; 44 } 45 46 47 G4double G4PreCompoundHe3::FactorialFactor(const G4double N, const G4double P) 48 { 49 return 61 50 (N-3.0)*(P-2.0)*( 62 51 (((N-2.0)*(P-1.0))/2.0) *( … … 64 53 ) 65 54 ); 66 } 67 68 G4double G4PreCompoundHe3::CoalescenceFactor(const G4double A) 69 { 70 return 243.0/(A*A); 71 } 72 73 G4double G4PreCompoundHe3::GetRj(const G4int NumberParticles, const G4int NumberCharged) 74 { 75 G4double rj = 0.0; 76 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2); 77 if(NumberCharged >=2 && (NumberParticles-NumberCharged) >= 1) { 78 rj = 3.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged)) 79 / static_cast<G4double>(denominator); 80 } 81 return rj; 82 } 55 } 56 57 G4double G4PreCompoundHe3::CoalescenceFactor(const G4double A) 58 { 59 return 243.0/(A*A); 60 } 61 62 63 64 G4double G4PreCompoundHe3::GetRj(const G4int NumberParticles, const G4int NumberCharged) 65 { 66 G4double rj = 0.0; 67 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2); 68 if(NumberCharged >=2 && (NumberParticles-NumberCharged) >= 1) rj = 3.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged))/static_cast<G4double>(denominator); 69 70 return rj; 71 } 72 73 74 83 75 84 76 //////////////////////////////////////////////////////////////////////////////////// … … 88 80 //OPT=3,4 Kalbach's parameterization 89 81 // 90 G4double G4PreCompoundHe3::CrossSection(const G4double K)82 G4double G4PreCompoundHe3::CrossSection(const G4double K) 91 83 { 92 84 ResidualA=GetRestA(); … … 122 114 //---------------- 123 115 // 124 G4double G4PreCompoundHe3::GetAlpha()125 {126 G4double C = 0.0;127 G4double aZ = GetZ() + GetRestZ();128 if (aZ <= 30)129 {130 131 }132 else if (aZ <= 50)133 {134 135 }136 else if (aZ < 70)137 {138 139 }140 else141 {142 143 }144 return 1.0 + C*(4.0/3.0);145 }116 G4double G4PreCompoundHe3::GetAlpha() 117 { 118 G4double C = 0.0; 119 G4double aZ = GetZ() + GetRestZ(); 120 if (aZ <= 30) 121 { 122 C = 0.10; 123 } 124 else if (aZ <= 50) 125 { 126 C = 0.1 + -((aZ-50.)/20.)*0.02; 127 } 128 else if (aZ < 70) 129 { 130 C = 0.08 + -((aZ-70.)/20.)*0.02; 131 } 132 else 133 { 134 C = 0.06; 135 } 136 return 1.0 + C*(4.0/3.0); 137 } 146 138 // 147 139 //-------------------- 148 140 // 149 G4double G4PreCompoundHe3::GetBeta()150 {151 return -GetCoulombBarrier();152 }141 G4double G4PreCompoundHe3::GetBeta() 142 { 143 return -GetCoulombBarrier(); 144 } 153 145 // 154 146 //********************* OPT=1,2 : Chatterjee's cross section ************************ … … 222 214 G4double ra=0.80; 223 215 224 //JMQ 13/02/09 increase of reduced radius to lower the barrier 225 // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 226 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra); 216 ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 227 217 ecsq = ec * ec; 228 218 p = p0 + p1/ec + p2/ecsq; … … 247 237 248 238 if (elab <= ec) { //start for E<Ec 249 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 239 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 250 240 } //end for E<Ec 251 241 else { //start for E>Ec -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundIon.cc
r962 r1007 24 24 // ******************************************************************** 25 25 // 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 26 //J. M. Quesada (August 2008). 27 //Based on previous work by V. Lara 40 28 // 41 29 … … 43 31 #include "G4PreCompoundParameters.hh" 44 32 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 }33 G4bool G4PreCompoundIon::IsItPossible(const G4Fragment& aFragment) 34 { 35 G4int pplus = aFragment.GetNumberOfCharged(); 36 G4int pneut = aFragment.GetNumberOfParticles()-pplus; 37 return (pneut >= (GetA()-GetZ()) && pplus >= GetZ()); 38 } 51 39 52 40 G4double G4PreCompoundIon:: … … 69 57 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 70 58 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; 59 G4double gj = (6.0/pi2)*GetA() * 60 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 76 61 77 62 G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0; … … 89 74 G4double Ej = std::max(0.0,eKin + GetBindingEnergy() -Aj); 90 75 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()); 76 77 G4double pA = 1.e-25*(3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()* 78 (eKin+GetBindingEnergy()))))/(pi * r0 * r0 * std::pow(GetRestA(),2.0/3.0) )* eKin*CrossSection(eKin) /(r0*std::pow(GetRestA(),1.0/3.0)) * CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P)* GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) ; 97 79 98 80 G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0); … … 100 82 G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0; 101 83 84 102 85 G4double Probability = pA * pB * pC; 103 86 104 87 return Probability; 105 88 } 89 90 91 92 93 94 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundModel.cc
r962 r1007 26 26 // 27 27 // $Id: G4PreCompoundModel.cc,v 1.17 2008/12/09 14:09:59 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02 -ref-02$28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNeutron.cc
r968 r1007 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4PreCompoundNeutron.cc,v 1.4 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: G4PreCompoundNeutron 36 // 37 // Author: V.Lara 38 // 39 // Modified: 40 // 21.08.2008 J. M. Quesada add choice of options 41 // 10.02.2009 J. M. Quesada set default opt3 26 //J. M. Quesada (August 2008). New source file 27 // 28 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse 29 // cross section option 42 30 // 43 44 31 #include "G4PreCompoundNeutron.hh" 45 32 46 G4ReactionProduct * G4PreCompoundNeutron::GetReactionProduct() const 47 { 48 G4ReactionProduct * theReactionProduct = 49 new G4ReactionProduct(G4Neutron::NeutronDefinition()); 50 theReactionProduct->SetMomentum(GetMomentum().vect()); 51 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 33 34 G4ReactionProduct * G4PreCompoundNeutron::GetReactionProduct() const 35 { 36 G4ReactionProduct * theReactionProduct = 37 new G4ReactionProduct(G4Neutron::NeutronDefinition()); 38 theReactionProduct->SetMomentum(GetMomentum().vect()); 39 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 52 40 #ifdef PRECOMPOUND_TEST 53 theReactionProduct->SetCreatorModel("G4PrecompoundModel");41 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 54 42 #endif 55 return theReactionProduct;56 }57 58 G4double G4PreCompoundNeutron::GetRj(const G4int NumberParticles, const G4int NumberCharged) 59 { 60 G4double rj = 0.0;61 if(NumberParticles > 0) rj = static_cast<G4double>(NumberParticles - NumberCharged)/62 63 return rj;64 }43 return theReactionProduct; 44 } 45 46 47 G4double G4PreCompoundNeutron::GetRj(const G4int NumberParticles, const G4int NumberCharged) 48 { 49 G4double rj = 0.0; 50 if(NumberParticles > 0) rj = static_cast<G4double>(NumberParticles - NumberCharged)/static_cast<G4double>(NumberParticles); 51 return rj; 52 } 65 53 66 54 //////////////////////////////////////////////////////////////////////////////////// … … 70 58 //OPT=3,4 Kalbach's parameterization 71 59 // 72 G4double G4PreCompoundNeutron::CrossSection(const G4double K)60 G4double G4PreCompoundNeutron::CrossSection(const G4double K) 73 61 { 74 62 ResidualA=GetRestA(); … … 104 92 //------- 105 93 // 106 G4double G4PreCompoundNeutron::GetAlpha()107 {108 94 G4double G4PreCompoundNeutron::GetAlpha() 95 { 96 // return 0.76+2.2/std::pow(GetRestA(),1.0/3.0); 109 97 return 0.76+2.2/ResidualAthrd; 110 }98 } 111 99 // 112 100 //------------ 113 101 // 114 G4double G4PreCompoundNeutron::GetBeta() 115 { 116 // return (2.12/std::pow(GetRestA(),2.0/3.0)-0.05)*MeV/GetAlpha(); 117 return (2.12/(ResidualAthrd*ResidualAthrd)-0.05)*MeV/GetAlpha(); 118 } 119 // 102 G4double G4PreCompoundNeutron::GetBeta() 103 { 104 // return (2.12/std::pow(GetRestA(),2.0/3.0)-0.05)*MeV/GetAlpha(); 105 return (2.12/(ResidualAthrd*ResidualAthrd)-0.05)*MeV/GetAlpha(); 106 } 107 // 108 109 110 120 111 121 112 //********************* OPT=1,2 : Chatterjee's cross section ************************ … … 165 156 G4double b,ecut,cut,ecut2,geom,elab; 166 157 167 158 //safety initialization 168 159 landa0=0; 169 160 landa1=0; … … 181 172 spill= 1.0e+18; 182 173 183 // PRECO xs for neutrons is choosen 174 175 176 // PRECO xs for neutrons is choosen 184 177 185 178 p0 = -312.; … … 208 201 xnulam = 1.; 209 202 etest = 32.; 210 // ** etest is the energy above which the rxn cross section is 211 // ** compared with the geometrical limit and the max taken. 212 // ** xnulam here is a dummy value to be used later. 203 // ** etest is the energy above which the rxn cross section is 204 // ** compared with the geometrical limit and the max taken. 205 // ** xnulam here is a dummy value to be used later. 206 207 213 208 214 209 a = -2.*p*ec + landa - nu/ecsq; -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNucleon.cc
r962 r1007 24 24 // ******************************************************************** 25 25 // 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 26 //J. M. Quesada (August 2008). 27 //Based on previous work by V. Lara 41 28 // 42 29 … … 45 32 46 33 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 } 34 { 35 G4int pplus = aFragment.GetNumberOfCharged(); 36 G4int pneut = aFragment.GetNumberOfParticles()-pplus; 37 return (pneut >= (GetA()-GetZ()) && pplus >= GetZ()); 38 } 39 52 40 53 41 G4double G4PreCompoundNucleon:: … … 57 45 if ( !IsItPossible(aFragment) ) return 0.0; 58 46 47 59 48 G4double U = aFragment.GetExcitationEnergy(); 60 49 G4double P = aFragment.GetNumberOfParticles(); … … 62 51 G4double N = P + H; 63 52 53 54 64 55 G4double g0 = (6.0/pi2)*aFragment.GetA() * 65 56 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); … … 67 58 G4double g1 = (6.0/pi2)*GetRestA() * 68 59 G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 60 61 69 62 70 63 G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0; … … 78 71 79 72 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); 73 G4double Probability =1.e-25* 1.0/pi2*2.0/(hbarc*hbarc*hbarc) * GetReducedMass() 74 * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) * eKin*CrossSection(eKin) * P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 * g1/(g0*g0); 83 75 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 } 76 77 if (GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())<0.0 || CrossSection(eKin) <0) 78 { std::ostringstream errOs; 79 G4cout<<"WARNING: NEGATIVE VALUES "<<G4endl; 80 errOs << "Rj=" << GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())<<G4endl; 81 errOs <<" xsec("<<eKin<<" MeV) ="<<CrossSection(eKin)<<G4endl; 82 errOs <<" A="<<GetA()<<" Z="<<GetZ()<<G4endl; 83 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 84 } 94 85 95 86 return Probability; -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundParameters.cc
r962 r1007 26 26 // 27 27 // $Id: G4PreCompoundParameters.cc,v 1.3 2006/06/29 20:59:29 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-02 -ref-02$28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundProton.cc
r968 r1007 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4PreCompoundProton.cc,v 1.4 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: G4PreCompoundProton 36 // 37 // Author: V.Lara 38 // 39 // Modified: 40 // 21.08.2008 J. M. Quesada added external choice of inverse cross section option 41 // 21.08.2008 J. M. Quesada added external choice for superimposed Coulomb barrier 42 // (if useSICB=true) 43 // 26 //J.M.Quesada (August 08). New source file 27 // 28 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse 29 // cross section option 30 // 31 // JMQ (06 September 2008) Also external choice has been added for: 32 // - superimposed Coulomb barrier (if useSICB=true) 44 33 45 34 #include "G4PreCompoundProton.hh" 46 35 47 G4ReactionProduct * G4PreCompoundProton::GetReactionProduct() const48 {49 G4ReactionProduct * theReactionProduct =50 new G4ReactionProduct(G4Proton::ProtonDefinition());51 theReactionProduct->SetMomentum(GetMomentum().vect());52 theReactionProduct->SetTotalEnergy(GetMomentum().e());36 G4ReactionProduct * G4PreCompoundProton::GetReactionProduct() const 37 { 38 G4ReactionProduct * theReactionProduct = 39 new G4ReactionProduct(G4Proton::ProtonDefinition()); 40 theReactionProduct->SetMomentum(GetMomentum().vect()); 41 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 53 42 #ifdef PRECOMPOUND_TEST 54 theReactionProduct->SetCreatorModel("G4PrecompoundModel");43 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 55 44 #endif 56 return theReactionProduct; 57 } 58 59 G4double G4PreCompoundProton::GetRj(const G4int NumberParticles, const G4int NumberCharged) 60 { 61 G4double rj = 0.0; 62 if(NumberParticles > 0) rj = static_cast<G4double>(NumberCharged)/static_cast<G4double>(NumberParticles); 63 return rj; 64 } 45 return theReactionProduct; 46 } 47 48 49 G4double G4PreCompoundProton::GetRj(const G4int NumberParticles, const G4int NumberCharged) 50 { 51 G4double rj = 0.0; 52 if(NumberParticles > 0) rj = static_cast<G4double>(NumberCharged)/static_cast<G4double>(NumberParticles); 53 return rj; 54 } 55 56 65 57 66 58 //////////////////////////////////////////////////////////////////////////////////// … … 71 63 //OPT=3 Kalbach's parameterization 72 64 // 73 G4double G4PreCompoundProton::CrossSection(const G4double K)65 G4double G4PreCompoundProton::CrossSection(const G4double K) 74 66 { 75 67 //G4cout<<" In G4PreCompoundProton OPTxs="<<OPTxs<<G4endl; … … 83 75 FragmentA=GetA()+GetRestA(); 84 76 FragmentAthrd=std::pow(FragmentA,0.33333); 77 78 85 79 86 80 if (OPTxs==0) return GetOpt0(K); … … 100 94 G4double G4PreCompoundProton::GetOpt0(const G4double K) 101 95 { 102 96 const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0(); 103 97 // cross section is now given in mb (r0 is in mm) for the sake of consistency 104 98 //with the rest of the options 105 99 return 1.e+25*pi*(r0*ResidualAthrd)*(r0*ResidualAthrd)*GetAlpha()*(1.+GetBeta()/K); 106 100 } 107 101 // 108 102 //------------ 109 103 // 110 G4double G4PreCompoundProton::GetAlpha()111 {112 113 G4double C = 0.0;114 if (aZ >= 70)115 {116 117 }118 else119 {120 121 }122 return 1.0 + C;123 }104 G4double G4PreCompoundProton::GetAlpha() 105 { 106 G4double aZ = static_cast<G4double>(GetRestZ()); 107 G4double C = 0.0; 108 if (aZ >= 70) 109 { 110 C = 0.10; 111 } 112 else 113 { 114 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 115 } 116 return 1.0 + C; 117 } 124 118 // 125 119 //------------------- 126 120 // 127 G4double G4PreCompoundProton::GetBeta()128 {121 G4double G4PreCompoundProton::GetBeta() 122 { 129 123 return -GetCoulombBarrier(); 130 } 131 // 132 124 } 125 // 126 127 128 133 129 //********************* OPT=1 : Chatterjee's cross section ************************ 134 130 //(fitting to cross section from Bechetti & Greenles OM potential) … … 136 132 G4double G4PreCompoundProton::GetOpt1(const G4double K) 137 133 { 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 134 G4double Kc=K; 135 136 // JMQ xsec is set constat above limit of validity 137 if (K>50) Kc=50; 138 139 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs; 140 G4double p, p0, p1, p2,Ec,delta,q,r,ji; 141 142 p0 = 15.72; 143 p1 = 9.65; 144 p2 = -449.0; 145 landa0 = 0.00437; 146 landa1 = -16.58; 147 mu0 = 244.7; 148 mu1 = 0.503; 149 nu0 = 273.1; 150 nu1 = -182.4; 151 nu2 = -1.872; 152 delta=0.; 153 154 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta); 155 p = p0 + p1/Ec + p2/(Ec*Ec); 156 landa = landa0*ResidualA + landa1; 157 mu = mu0*std::pow(ResidualA,mu1); 158 nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 159 q = landa - nu/(Ec*Ec) - 2*p*Ec; 160 r = mu + 2*nu/Ec + p*(Ec*Ec); 161 162 ji=std::max(Kc,Ec); 163 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;} 164 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;} 165 if (xs <0.0) {xs=0.0;} 166 167 return xs; 172 168 173 169 } 170 171 174 172 175 173 //************* OPT=2 : Welisch's proton reaction cross section ************************ … … 180 178 G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0); 181 179 182 183 180 //This is redundant when the Coulomb barrier is overimposed to all cross sections 181 //It should be kept when Coulomb barrier only imposed at OPTxs=2 184 182 185 183 if(!useSICB && K<=theCoulombBarrier) return xine_th=0.0; … … 318 316 } //end for E>Ec 319 317 320 return sig; 321 } 318 return sig;} 319 320 322 321 323 322 // ************************** end of cross sections ******************************* -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTransitions.cc
r962 r1007 24 24 // ******************************************************************** 25 25 // 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 26 // 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: 27 //J. M. Quesada (Feb. 08). Base on previous work by V. Lara. 28 // New transition probabilities. Several bugs fixed. 29 // JMQ (06 September 2008) Also external choices have been added for: 41 30 // - "never go back" hipothesis (useNGB=true) 42 31 // - CEM transition probabilities (useCEMtr=true) 43 44 32 #include "G4PreCompoundTransitions.hh" 45 33 #include "G4HadronicException.hh" … … 196 184 // G4cout<<"N="<<N<<" P="<<P<<" H="<<H<<G4endl; 197 185 // G4cout<<"l+ ="<<TransitionProb1<<" l- ="<< TransitionProb2<<" l0 ="<< TransitionProb3<<G4endl; 198 return TransitionProb1 + TransitionProb2 + TransitionProb3; 199 } 186 return TransitionProb1 + TransitionProb2 + TransitionProb3;} 200 187 201 188 else { … … 228 215 return TransitionProb1 + TransitionProb2 + TransitionProb3; 229 216 } 217 218 230 219 } 231 220 … … 248 237 } 249 238 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)) { 239 // 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, PROVIDED that there are charged particles 240 if(deltaN < 0 && G4UniformRand() <= static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) && (result.GetNumberOfCharged() >= 1)) 255 241 result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative! 256 } 257 258 // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on 259 // number of charged vs. number of particles fails 242 243 244 //JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on number of charged vs. number of particles fails 260 245 result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2); 261 246 result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2); 262 247 263 248 // With weight Z/A, number of charged particles is increased with +1 264 249 if ( ( deltaN > 0 ) && 265 250 (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/ 266 251 std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.))) 267 252 { 268 253 result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTriton.cc
r968 r1007 25 25 // 26 26 // 27 // $Id: G4PreCompoundTriton.cc,v 1.5 2009/02/13 18:57:32 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 // 30 // ------------------------------------------------------------------- 31 // 32 // GEANT4 Class file 33 // 34 // 35 // File name: G4PreCompoundTriton 36 // 37 // Author: V.Lara 38 // 39 // Modified: 40 // 21.08.2008 J. M. Quesada add choice of options 41 // 10.02.2009 J. M. Quesada set default opt1 42 // 27 //J.M.Quesada (August 08). New source file 28 // 29 // Modif (21 August 2008) by J. M. Quesada for external choice of inverse 30 // cross section option 43 31 44 32 #include "G4PreCompoundTriton.hh" 45 33 46 34 47 G4ReactionProduct * G4PreCompoundTriton::GetReactionProduct() const48 {49 G4ReactionProduct * theReactionProduct =50 new G4ReactionProduct(G4Triton::TritonDefinition());51 theReactionProduct->SetMomentum(GetMomentum().vect());52 theReactionProduct->SetTotalEnergy(GetMomentum().e());35 G4ReactionProduct * G4PreCompoundTriton::GetReactionProduct() const 36 { 37 G4ReactionProduct * theReactionProduct = 38 new G4ReactionProduct(G4Triton::TritonDefinition()); 39 theReactionProduct->SetMomentum(GetMomentum().vect()); 40 theReactionProduct->SetTotalEnergy(GetMomentum().e()); 53 41 #ifdef PRECOMPOUND_TEST 54 theReactionProduct->SetCreatorModel("G4PrecompoundModel");42 theReactionProduct->SetCreatorModel("G4PrecompoundModel"); 55 43 #endif 56 return theReactionProduct;57 }58 59 G4double G4PreCompoundTriton::FactorialFactor(const G4double N, const G4double P)60 {61 return44 return theReactionProduct; 45 } 46 47 G4double G4PreCompoundTriton::FactorialFactor(const G4double N, const G4double P) 48 { 49 return 62 50 (N-3.0)*(P-2.0)*( 63 51 (((N-2.0)*(P-1.0))/2.0) *( … … 65 53 ) 66 54 ); 67 }55 } 68 56 69 G4double G4PreCompoundTriton::CoalescenceFactor(const G4double A) 70 { 71 return 243.0/(A*A); 72 } 73 74 G4double G4PreCompoundTriton::GetRj(const G4int NumberParticles, const G4int NumberCharged) 75 { 76 G4double rj = 0.0; 77 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2); 78 if(NumberCharged >= 1 && (NumberParticles-NumberCharged) >= 2) { 79 rj = 3.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1)) 80 /static_cast<G4double>(denominator); 81 } 82 return rj; 83 } 57 G4double G4PreCompoundTriton::CoalescenceFactor(const G4double A) 58 { 59 return 243.0/(A*A); 60 } 61 62 63 G4double G4PreCompoundTriton::GetRj(const G4int NumberParticles, const G4int NumberCharged) 64 { 65 G4double rj = 0.0; 66 G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2); 67 if(NumberCharged >= 1 && (NumberParticles-NumberCharged) >= 2) rj = 3.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); 68 return rj; 69 } 70 71 72 84 73 85 74 //////////////////////////////////////////////////////////////////////////////////// … … 89 78 //OPT=3,4 Kalbach's parameterization 90 79 // 91 G4double G4PreCompoundTriton::CrossSection(const G4double K)80 G4double G4PreCompoundTriton::CrossSection(const G4double K) 92 81 { 93 82 ResidualA=GetRestA(); … … 122 111 //--------- 123 112 // 124 G4double G4PreCompoundTriton::GetAlpha()125 {126 G4double C = 0.0;127 G4double aZ = GetZ() + GetRestZ();128 if (aZ >= 70)129 {130 131 }132 else133 {134 135 }136 137 return 1.0 + C/3.0;138 }113 G4double G4PreCompoundTriton::GetAlpha() 114 { 115 G4double C = 0.0; 116 G4double aZ = GetZ() + GetRestZ(); 117 if (aZ >= 70) 118 { 119 C = 0.10; 120 } 121 else 122 { 123 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 124 } 125 126 return 1.0 + C/3.0; 127 } 139 128 // 140 129 //------------- 141 130 // 142 G4double G4PreCompoundTriton::GetBeta()143 {144 return -GetCoulombBarrier();145 }131 G4double G4PreCompoundTriton::GetBeta() 132 { 133 return -GetCoulombBarrier(); 134 } 146 135 // 147 136 //********************* OPT=1,2 : Chatterjee's cross section ************************ … … 194 183 195 184 G4double landa, mu, nu, p , signor(1.),sig; 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 185 G4double ec,ecsq,xnulam,etest(0.),a; 186 G4double b,ecut,cut,ecut2,geom,elab; 187 188 189 G4double flow = 1.e-18; 190 G4double spill= 1.e+18; 191 192 193 G4double p0 = -21.45; 194 G4double p1 = 484.7; 195 G4double p2 = -1608.; 196 G4double landa0 = 0.0186; 197 G4double landa1 = -8.90; 198 G4double mu0 = 686.3; 199 G4double mu1 = 0.325; 200 G4double nu0 = 368.9; 201 G4double nu1 = -522.2; 202 G4double nu2 = -4.998; 214 203 215 204 G4double ra=0.80; 216 205 217 //JMQ 13/02/09 increase of reduced radius to lower the barrier 218 // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 219 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra); 220 ecsq = ec * ec; 221 p = p0 + p1/ec + p2/ecsq; 222 landa = landa0*ResidualA + landa1; 223 a = std::pow(ResidualA,mu1); 224 mu = mu0 * a; 225 nu = a* (nu0+nu1*ec+nu2*ecsq); 226 xnulam = nu / landa; 227 if (xnulam > spill) xnulam=0.; 228 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam); 229 230 a = -2.*p*ec + landa - nu/ecsq; 231 b = p*ecsq + mu + 2.*nu/ec; 232 ecut = 0.; 233 cut = a*a - 4.*p*b; 234 if (cut > 0.) ecut = std::sqrt(cut); 235 ecut = (ecut-a) / (p+p); 236 ecut2 = ecut; 237 if (cut < 0.) ecut2 = ecut - 2.; 238 elab = K * FragmentA / ResidualA; 239 sig = 0.; 240 241 if (elab <= ec) { //start for E<Ec 242 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 243 } //end for E<Ec 244 else { //start for E>Ec 245 sig = (landa*elab+mu+nu/elab) * signor; 246 geom = 0.; 247 if (xnulam < flow || elab < etest) return sig; 248 geom = std::sqrt(theA*K); 249 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 250 geom = 31.416 * geom * geom; 251 sig = std::max(geom,sig); 252 } //end for E>Ec 253 return sig; 206 ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 207 ecsq = ec * ec; 208 p = p0 + p1/ec + p2/ecsq; 209 landa = landa0*ResidualA + landa1; 210 a = std::pow(ResidualA,mu1); 211 mu = mu0 * a; 212 nu = a* (nu0+nu1*ec+nu2*ecsq); 213 xnulam = nu / landa; 214 if (xnulam > spill) xnulam=0.; 215 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam); 216 217 a = -2.*p*ec + landa - nu/ecsq; 218 b = p*ecsq + mu + 2.*nu/ec; 219 ecut = 0.; 220 cut = a*a - 4.*p*b; 221 if (cut > 0.) ecut = std::sqrt(cut); 222 ecut = (ecut-a) / (p+p); 223 ecut2 = ecut; 224 if (cut < 0.) ecut2 = ecut - 2.; 225 elab = K * FragmentA / ResidualA; 226 sig = 0.; 227 228 if (elab <= ec) { //start for E<Ec 229 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 230 } //end for E<Ec 231 else { //start for E>Ec 232 sig = (landa*elab+mu+nu/elab) * signor; 233 geom = 0.; 234 if (xnulam < flow || elab < etest) return sig; 235 geom = std::sqrt(theA*K); 236 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 237 geom = 31.416 * geom * geom; 238 sig = std::max(geom,sig); 239 } //end for E>Ec 240 return sig; 254 241 255 242 } 256 243 257 244 // ************************** end of cross sections ******************************* 245 246 247 248 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundFragment.cc
r962 r1007 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 $28 26 // 29 // 30 // 27 //J. M. Quesada (August 2008). 28 //Based on previous work by V. Lara 31 29 // 32 30 -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundIon.cc
r962 r1007 26 26 // 27 27 // $Id: G4VPreCompoundIon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02 -ref-02$28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // by V. Lara -
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundNucleon.cc
r962 r1007 27 27 28 28 // $Id: G4VPreCompoundNucleon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $ 29 // GEANT4 tag $Name: geant4-09-02 -ref-02$29 // GEANT4 tag $Name: geant4-09-02 $ 30 30 31 31 //
Note: See TracChangeset
for help on using the changeset viewer.