00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 #ifndef G4GEMProbability_h
00042 #define G4GEMProbability_h 1
00043
00044 #include <CLHEP/Units/SystemOfUnits.h>
00045
00046 #include "G4VEmissionProbability.hh"
00047 #include "G4VLevelDensityParameter.hh"
00048 #include "G4EvaporationLevelDensityParameter.hh"
00049 #include "G4VCoulombBarrier.hh"
00050 #include "G4PairingCorrection.hh"
00051 #include "G4Pow.hh"
00052
00053 class G4GEMProbability : public G4VEmissionProbability
00054 {
00055 public:
00056
00057 G4GEMProbability(G4int anA, G4int aZ, G4double aSpin);
00058
00059 virtual ~G4GEMProbability();
00060
00061 G4double EmissionProbability(const G4Fragment & fragment, G4double anEnergy);
00062
00063 inline G4int GetZ_asInt(void) const;
00064
00065 inline G4int GetA_asInt(void) const;
00066
00067 inline G4double GetZ(void) const;
00068
00069 inline G4double GetA(void) const;
00070
00071 inline G4double GetSpin(void) const;
00072
00073 inline G4double GetNormalization(void) const;
00074
00075 inline void SetCoulomBarrier(const G4VCoulombBarrier * aCoulombBarrierStrategy);
00076
00077 inline G4double GetCoulombBarrier(const G4Fragment& fragment) const;
00078
00079 inline G4double CalcAlphaParam(const G4Fragment & ) const;
00080
00081 inline G4double CalcBetaParam(const G4Fragment & ) const;
00082
00083 private:
00084
00085 G4double CalcProbability(const G4Fragment & fragment,
00086 G4double MaximalKineticEnergy,
00087 G4double V);
00088
00089 inline G4double CCoeficient(G4int) const;
00090
00091 inline G4double I0(G4double t);
00092 inline G4double I1(G4double t, G4double tx);
00093 inline G4double I2(G4double s0, G4double sx);
00094 G4double I3(G4double s0, G4double sx);
00095
00096
00097 G4GEMProbability();
00098 G4GEMProbability(const G4GEMProbability &right);
00099 const G4GEMProbability & operator=(const G4GEMProbability &right);
00100 G4bool operator==(const G4GEMProbability &right) const;
00101 G4bool operator!=(const G4GEMProbability &right) const;
00102
00103
00104 G4Pow* fG4pow;
00105 G4PairingCorrection* fPairCorr;
00106
00107 G4VLevelDensityParameter * theEvapLDPptr;
00108
00109 G4int theA;
00110 G4int theZ;
00111
00112
00113 G4double Spin;
00114
00115
00116 const G4VCoulombBarrier * theCoulombBarrierPtr;
00117
00118
00119 G4double Normalization;
00120
00121 protected:
00122
00123 G4double fPlanck;
00124
00125
00126 std::vector<G4double> ExcitEnergies;
00127
00128
00129 std::vector<G4double> ExcitSpins;
00130
00131
00132 std::vector<G4double> ExcitLifetimes;
00133
00134 };
00135
00136 inline G4int G4GEMProbability::GetZ_asInt(void) const
00137 {
00138 return theZ;
00139 }
00140
00141 inline G4int G4GEMProbability::GetA_asInt(void) const
00142 {
00143 return theA;
00144 }
00145
00146 inline G4double G4GEMProbability::GetZ(void) const
00147 {
00148 return theZ;
00149 }
00150
00151 inline G4double G4GEMProbability::GetA(void) const
00152 {
00153 return theA;
00154 }
00155
00156 inline G4double G4GEMProbability::GetSpin(void) const
00157 {
00158 return Spin;
00159 }
00160
00161 inline G4double G4GEMProbability::GetNormalization(void) const
00162 {
00163 return Normalization;
00164 }
00165
00166 inline void
00167 G4GEMProbability::SetCoulomBarrier(const G4VCoulombBarrier * aCoulombBarrierStrategy)
00168 {
00169 theCoulombBarrierPtr = aCoulombBarrierStrategy;
00170 }
00171
00172 inline G4double
00173 G4GEMProbability::GetCoulombBarrier(const G4Fragment& fragment) const
00174 {
00175 G4double res = 0.0;
00176 if (theCoulombBarrierPtr) {
00177 G4int Acomp = fragment.GetA_asInt();
00178 G4int Zcomp = fragment.GetZ_asInt();
00179 res = theCoulombBarrierPtr->GetCoulombBarrier(Acomp-theA, Zcomp-theZ,
00180 fragment.GetExcitationEnergy() -
00181 fPairCorr->GetPairingCorrection(Acomp,Zcomp));
00182 }
00183 return res;
00184 }
00185
00186 inline G4double G4GEMProbability::CCoeficient(G4int aZ) const
00187 {
00188
00189
00190
00191 G4double C = 0.0;
00192 if (aZ >= 50){
00193 C=-0.10/G4double(theA);
00194 } else if (aZ > 20) {
00195 C=(0.123482-0.00534691*aZ-0.0000610624*aZ*aZ+5.93719*1e-7*aZ*aZ*aZ+
00196 1.95687*1e-8*aZ*aZ*aZ*aZ)/G4double(theA);
00197 }
00198 return C;
00199 }
00200
00201
00202 inline G4double G4GEMProbability::CalcAlphaParam(const G4Fragment & fragment) const
00203 {
00204
00205
00206 G4double res;
00207 if(GetZ_asInt() == 0) {
00208 res = 0.76+1.93/fG4pow->Z13(fragment.GetA_asInt()-GetA_asInt());
00209 } else {
00210 res = 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ_asInt());
00211 }
00212 return res;
00213 }
00214
00215 inline G4double
00216 G4GEMProbability::CalcBetaParam(const G4Fragment & fragment) const
00217 {
00218
00219
00220 G4double res;
00221 if(GetZ_asInt() == 0) {
00222 res = (1.66/fG4pow->Z23(fragment.GetA_asInt()-GetA_asInt())-0.05)*CLHEP::MeV/
00223 CalcAlphaParam(fragment);
00224 } else {
00225 res = -GetCoulombBarrier(fragment);
00226 }
00227 return res;
00228 }
00229
00230 inline G4double G4GEMProbability::I0(G4double t)
00231 {
00232 return std::exp(t) - 1.0;
00233 }
00234
00235 inline G4double G4GEMProbability::I1(G4double t, G4double tx)
00236 {
00237 return (t - tx + 1.0)*std::exp(tx) - t - 1.0;
00238 }
00239
00240
00241 inline G4double G4GEMProbability::I2(G4double s0, G4double sx)
00242 {
00243 G4double S = 1.0/std::sqrt(s0);
00244 G4double Sx = 1.0/std::sqrt(sx);
00245
00246 G4double p1 = S*S*S*( 1.0 + S*S*( 1.5 + 3.75*S*S) );
00247 G4double p2 = Sx*Sx*Sx*( 1.0 + Sx*Sx*( 1.5 + 3.75*Sx*Sx) )*std::exp(sx-s0);
00248
00249 return p1-p2;
00250 }
00251
00252
00253 #endif