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 #include "G4GEMChannel.hh"
00039 #include "G4PairingCorrection.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4Pow.hh"
00043
00044 G4GEMChannel::G4GEMChannel(G4int theA, G4int theZ, const G4String & aName,
00045 G4GEMProbability * aEmissionStrategy,
00046 G4VCoulombBarrier * aCoulombBarrier) :
00047 G4VEvaporationChannel(aName),
00048 A(theA),
00049 Z(theZ),
00050 theEvaporationProbabilityPtr(aEmissionStrategy),
00051 theCoulombBarrierPtr(aCoulombBarrier),
00052 EmissionProbability(0.0),
00053 MaximalKineticEnergy(-CLHEP::GeV)
00054 {
00055 theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
00056 MyOwnLevelDensity = true;
00057 EvaporatedMass = G4NucleiProperties::GetNuclearMass(A, Z);
00058 ResidualMass = CoulombBarrier = 0.0;
00059 fG4pow = G4Pow::GetInstance();
00060 ResidualZ = ResidualA = 0;
00061 }
00062
00063 G4GEMChannel::~G4GEMChannel()
00064 {
00065 if (MyOwnLevelDensity) { delete theLevelDensityPtr; }
00066 }
00067
00068 G4double G4GEMChannel::GetEmissionProbability(G4Fragment* fragment)
00069 {
00070 G4int anA = fragment->GetA_asInt();
00071 G4int aZ = fragment->GetZ_asInt();
00072 ResidualA = anA - A;
00073 ResidualZ = aZ - Z;
00074
00075
00076
00077
00078 if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
00079 (ResidualA == ResidualZ && ResidualA > 1))
00080 {
00081 CoulombBarrier = 0.0;
00082 MaximalKineticEnergy = -CLHEP::GeV;
00083 EmissionProbability = 0.0;
00084 }
00085 else
00086 {
00087
00088
00089
00090
00091
00092
00093
00094 G4double ExEnergy = fragment->GetExcitationEnergy() -
00095 G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
00096
00097
00098
00099 if( ExEnergy <= 0.0) {
00100 CoulombBarrier = 0.0;
00101 MaximalKineticEnergy = -1000.0*MeV;
00102 EmissionProbability = 0.0;
00103
00104 } else {
00105
00106 ResidualMass = G4NucleiProperties::GetNuclearMass(ResidualA, ResidualZ);
00107
00108
00109 CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
00110
00111
00112
00113 MaximalKineticEnergy =
00114 CalcMaximalKineticEnergy(fragment->GetGroundStateMass()+ExEnergy);
00115
00116
00117
00118 if (MaximalKineticEnergy <= 0.0)
00119 {
00120 EmissionProbability = 0.0;
00121 }
00122 else
00123 {
00124
00125 EmissionProbability =
00126 theEvaporationProbabilityPtr->EmissionProbability(*fragment,
00127 MaximalKineticEnergy);
00128 }
00129 }
00130 }
00131
00132 return EmissionProbability;
00133 }
00134
00135 G4FragmentVector * G4GEMChannel::BreakUp(const G4Fragment & theNucleus)
00136 {
00137 G4double EvaporatedKineticEnergy = CalcKineticEnergy(theNucleus);
00138 G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
00139
00140 G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
00141 (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
00142
00143 momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
00144
00145 G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
00146 EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
00147 G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
00148
00149 G4double theExEnergy = theNucleus.GetExcitationEnergy();
00150 G4double theMass = theNucleus.GetGroundStateMass();
00151 G4double ResidualEnergy =
00152 theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
00153
00154 G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
00155 ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
00156
00157 G4Fragment * ResidualFragment = new G4Fragment( ResidualA, ResidualZ, ResidualMomentum );
00158
00159 G4FragmentVector * theResult = new G4FragmentVector;
00160
00161 theResult->push_back(EvaporatedFragment);
00162 theResult->push_back(ResidualFragment);
00163 return theResult;
00164 }
00165
00166 G4double G4GEMChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
00167
00168
00169 {
00170 G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
00171 (2.0*NucleusTotalE) - EvaporatedMass - CoulombBarrier;
00172
00173 return T;
00174 }
00175
00176 G4double G4GEMChannel::CalcKineticEnergy(const G4Fragment & fragment)
00177
00178 {
00179 G4double U = fragment.GetExcitationEnergy();
00180
00181 G4double Alpha = theEvaporationProbabilityPtr->CalcAlphaParam(fragment);
00182 G4double Beta = theEvaporationProbabilityPtr->CalcBetaParam(fragment);
00183
00184 G4double Normalization = theEvaporationProbabilityPtr->GetNormalization();
00185
00186
00187
00188 G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(ResidualA,ResidualZ);
00189 G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
00190 G4double Ex = Ux + delta0;
00191 G4double InitialLevelDensity;
00192
00193
00194
00195
00196
00197 G4double deltaCN =
00198 G4PairingCorrection::GetInstance()->GetPairingCorrection(fragment.GetA_asInt(),
00199 fragment.GetZ_asInt());
00200 G4double aCN = theLevelDensityPtr->LevelDensityParameter(fragment.GetA_asInt(),
00201 fragment.GetZ_asInt(),U-deltaCN);
00202 G4double UxCN = (2.5 + 150.0/fragment.GetA())*MeV;
00203 G4double ExCN = UxCN + deltaCN;
00204 G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
00205
00206
00207
00208 if ( U < ExCN )
00209 {
00210 G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0
00211 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
00212 InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
00213 }
00214 else
00215 {
00216 G4double x = U-deltaCN;
00217 G4double x1 = std::sqrt(aCN*x);
00218 InitialLevelDensity = (pi/12.0)*std::exp(2*x1)/(x*std::sqrt(x1));
00219
00220
00221 }
00222
00223 const G4double Spin = theEvaporationProbabilityPtr->GetSpin();
00224
00225
00226
00227 G4double gg = (2.0*Spin+1.0)*EvaporatedMass/(pi2* hbarc*hbarc);
00228
00229
00230
00231 G4double Rb = 0.0;
00232 if (A > 4)
00233 {
00234 G4double Ad = fG4pow->Z13(ResidualA);
00235 G4double Aj = fG4pow->Z13(A);
00236
00237 Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
00238 Rb *= fermi;
00239 }
00240 else if (A>1)
00241 {
00242 G4double Ad = fG4pow->Z13(ResidualA);
00243 G4double Aj = fG4pow->Z13(A);
00244 Rb=1.5*(Aj+Ad)*fermi;
00245 }
00246 else
00247 {
00248 G4double Ad = fG4pow->Z13(ResidualA);
00249 Rb = 1.5*Ad*fermi;
00250 }
00251
00252 G4double GeometricalXS = pi*Rb*Rb;
00253
00254 G4double ConstantFactor = gg*GeometricalXS*Alpha/InitialLevelDensity;
00255 ConstantFactor *= pi/12.0;
00256
00257
00258
00259 G4double theEnergy = MaximalKineticEnergy + CoulombBarrier;
00260 G4double KineticEnergy;
00261 G4double Probability;
00262
00263 do
00264 {
00265 KineticEnergy = CoulombBarrier + G4UniformRand()*MaximalKineticEnergy;
00266 Probability = ConstantFactor*(KineticEnergy + Beta);
00267 G4double a =
00268 theLevelDensityPtr->LevelDensityParameter(ResidualA,ResidualZ,theEnergy-KineticEnergy-delta0);
00269 G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
00270
00271
00272 if ( theEnergy-KineticEnergy < Ex)
00273 {
00274 G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0
00275 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
00276 Probability *= std::exp((theEnergy-KineticEnergy-E0)/T)/T;
00277 }
00278 else
00279 {
00280 Probability *= std::exp(2*std::sqrt(a*(theEnergy-KineticEnergy-delta0)))/
00281 std::pow(a*fG4pow->powN(theEnergy-KineticEnergy-delta0,5), 0.25);
00282 }
00283 }
00284 while (Normalization*G4UniformRand() > Probability);
00285
00286 return KineticEnergy;
00287 }
00288
00289
00290 G4ThreeVector G4GEMChannel::IsotropicVector(const G4double Magnitude)
00291
00292
00293 {
00294 G4double CosTheta = 1.0 - 2.0*G4UniformRand();
00295 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
00296 G4double Phi = twopi*G4UniformRand();
00297 G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
00298 Magnitude*std::sin(Phi)*SinTheta,
00299 Magnitude*CosTheta);
00300 return Vector;
00301 }
00302
00303
00304