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 #include "G4HETCFragment.hh"
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4PreCompoundParameters.hh"
00037
00038 G4HETCFragment::
00039 G4HETCFragment(const G4ParticleDefinition* part,
00040 G4VCoulombBarrier* aCoulombBarrier)
00041 : G4VPreCompoundFragment(part, aCoulombBarrier)
00042 {
00043 G4double r0 = theParameters->Getr0();
00044 r2norm = r0*r0/(CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc*CLHEP::hbarc);
00045 }
00046
00047 G4HETCFragment::~G4HETCFragment()
00048 {}
00049
00050 G4double G4HETCFragment::
00051 CalcEmissionProbability(const G4Fragment & aFragment)
00052 {
00053 if (GetEnergyThreshold() <= 0.0)
00054 {
00055 theEmissionProbability = 0.0;
00056 return 0.0;
00057 }
00058
00059
00060 G4double LowerLimit = theCoulombBarrier;
00061
00062
00063
00064 G4double UpperLimit = GetMaximalKineticEnergy();
00065
00066 theEmissionProbability =
00067 IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment);
00068
00069 return theEmissionProbability;
00070 }
00071
00072 G4double G4HETCFragment::
00073 IntegrateEmissionProbability(const G4double & Low, const G4double & Up,
00074 const G4Fragment & aFragment)
00075 {
00076
00077 if ( !IsItPossible(aFragment) ) { return 0.0; }
00078
00079 G4double U = aFragment.GetExcitationEnergy();
00080
00081 G4int P = aFragment.GetNumberOfParticles();
00082 G4int H = aFragment.GetNumberOfHoles();
00083 G4int N = P + H;
00084 G4int Pb = P - GetA();
00085 G4int Nb = Pb + H;
00086 if (Nb <= 0.0) { return 0.0; }
00087 G4double ga = (6.0/pi2)*aFragment.GetA()*theParameters->GetLevelDensity();
00088 G4double gb = (6.0/pi2)*GetRestA()*theParameters->GetLevelDensity();
00089
00090 G4double A = G4double(P*P+H*H+P-3*H)/(4.0*ga);
00091 G4double Ab = G4double(Pb*Pb+H*H+Pb-3*H)/(4.0*gb);
00092 U = std::max(U-A,0.0);
00093 if (U <= 0.0) { return 0.0; }
00094
00095 G4int Pf = P;
00096 G4int Hf = H;
00097 G4int Nf = N-1;
00098 for (G4int i = 1; i < GetA(); ++i)
00099 {
00100 Pf *= (P-i);
00101 Hf *= (H-i);
00102 Nf *= (N-1-i);
00103 }
00104
00105 G4double X = std::max(Up - Ab + GetBeta(),0.0);
00106 G4double Y = std::max(Up - Ab - Low, 0.0);
00107
00108 G4double Probability = r2norm*GetSpinFactor()*GetReducedMass()*GetAlpha()
00109 *g4pow->Z23(GetRestA())*Pf*Hf*Nf*K(aFragment)*(X/Nb - Y/(Nb+1))
00110 *U*g4pow->powN(gb*Y,Nb)/g4pow->powN(ga*U,N);
00111
00112
00113
00114
00115
00116
00117
00118 return Probability;
00119 }