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 #include "G4PreCompoundFragment.hh"
00038
00039 G4PreCompoundFragment::
00040 G4PreCompoundFragment(const G4ParticleDefinition* part,
00041 G4VCoulombBarrier* aCoulombBarrier)
00042 : G4VPreCompoundFragment(part,aCoulombBarrier)
00043 {}
00044
00045 G4PreCompoundFragment::~G4PreCompoundFragment()
00046 {}
00047
00048 G4double G4PreCompoundFragment::
00049 CalcEmissionProbability(const G4Fragment & aFragment)
00050 {
00051
00052
00053
00054 G4double limit = 0.0;
00055 if(OPTxs==0 || useSICB) { limit = theCoulombBarrier; }
00056 if (GetMaximalKineticEnergy() <= limit)
00057 {
00058 theEmissionProbability = 0.0;
00059 return 0.0;
00060 }
00061
00062
00063
00064
00065 G4double LowerLimit = limit;
00066
00067
00068
00069 G4double UpperLimit = GetMaximalKineticEnergy();
00070
00071 theEmissionProbability =
00072 IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment);
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 return theEmissionProbability;
00083 }
00084
00085 G4double G4PreCompoundFragment::
00086 IntegrateEmissionProbability(G4double Low, G4double Up,
00087 const G4Fragment & aFragment)
00088 {
00089 static const G4int N = 10;
00090
00091 static const G4double w[N] = {
00092 0.0666713443086881,
00093 0.149451349150581,
00094 0.219086362515982,
00095 0.269266719309996,
00096 0.295524224714753,
00097 0.295524224714753,
00098 0.269266719309996,
00099 0.219086362515982,
00100 0.149451349150581,
00101 0.0666713443086881
00102 };
00103 static const G4double x[N] = {
00104 -0.973906528517172,
00105 -0.865063366688985,
00106 -0.679409568299024,
00107 -0.433395394129247,
00108 -0.148874338981631,
00109 0.148874338981631,
00110 0.433395394129247,
00111 0.679409568299024,
00112 0.865063366688985,
00113 0.973906528517172
00114 };
00115
00116 G4double Total = 0.0;
00117
00118 for (G4int i = 0; i < N; ++i)
00119 {
00120 G4double KineticE = 0.5*((Up-Low)*x[i]+(Up+Low));
00121 Total += w[i]*ProbabilityDistributionFunction(KineticE, aFragment);
00122 }
00123 Total *= 0.5*(Up-Low);
00124 return Total;
00125 }
00126
00127 G4double G4PreCompoundFragment::
00128 GetKineticEnergy(const G4Fragment & aFragment)
00129 {
00130
00131 G4double V = 0.0;
00132 if(OPTxs==0 || useSICB) { V = theCoulombBarrier; }
00133
00134 G4double Tmax = GetMaximalKineticEnergy();
00135 if(Tmax < V) { return 0.0; }
00136 G4double T(0.0);
00137 G4double Probability(1.0);
00138 G4double maxProbability = GetEmissionProbability();
00139 do
00140 {
00141 T = V + G4UniformRand()*(Tmax-V);
00142 Probability = ProbabilityDistributionFunction(T,aFragment);
00143 } while (maxProbability*G4UniformRand() > Probability);
00144 return T;
00145 }