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