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 #include "G4EvaporationChannel.hh"
00042 #include "G4PairingCorrection.hh"
00043 #include "G4NucleiProperties.hh"
00044 #include "G4Pow.hh"
00045 #include "G4EvaporationLevelDensityParameter.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "Randomize.hh"
00049 #include "G4Alpha.hh"
00050
00051 G4EvaporationChannel::G4EvaporationChannel(G4int anA, G4int aZ,
00052 const G4String & aName,
00053 G4EvaporationProbability* aEmissionStrategy,
00054 G4VCoulombBarrier* aCoulombBarrier):
00055 G4VEvaporationChannel(aName),
00056 theA(anA),
00057 theZ(aZ),
00058 theEvaporationProbabilityPtr(aEmissionStrategy),
00059 theCoulombBarrierPtr(aCoulombBarrier),
00060 EmissionProbability(0.0),
00061 MaximalKineticEnergy(-1000.0)
00062 {
00063 ResidualA = 0;
00064 ResidualZ = 0;
00065 ResidualMass = CoulombBarrier=0.0;
00066 EvaporatedMass = G4NucleiProperties::GetNuclearMass(theA, theZ);
00067 theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
00068 }
00069
00070 G4EvaporationChannel::G4EvaporationChannel():
00071 G4VEvaporationChannel(""),
00072 theA(0),
00073 theZ(0),
00074 theEvaporationProbabilityPtr(0),
00075 theCoulombBarrierPtr(0),
00076 EmissionProbability(0.0),
00077 MaximalKineticEnergy(-1000.0)
00078 {
00079 ResidualA = 0;
00080 ResidualZ = 0;
00081 EvaporatedMass = ResidualMass = CoulombBarrier = 0.0;
00082 theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
00083 }
00084
00085 G4EvaporationChannel::~G4EvaporationChannel()
00086 {
00087 delete theLevelDensityPtr;
00088 }
00089
00090
00091
00092
00093 G4double G4EvaporationChannel::GetEmissionProbability(G4Fragment* fragment)
00094 {
00095
00096 theEvaporationProbabilityPtr->SetOPTxs(OPTxs);
00097
00098 theEvaporationProbabilityPtr->UseSICB(useSICB);
00099
00100 G4int FragmentA = fragment->GetA_asInt();
00101 G4int FragmentZ = fragment->GetZ_asInt();
00102 ResidualA = FragmentA - theA;
00103 ResidualZ = FragmentZ - theZ;
00104
00105
00106
00107
00108 G4double ExEnergy = fragment->GetExcitationEnergy() -
00109 G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ);
00110
00111
00112 if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
00113 (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) {
00114 CoulombBarrier = ResidualMass = 0.0;
00115 MaximalKineticEnergy = -1000.0*MeV;
00116 EmissionProbability = 0.0;
00117 } else {
00118 ResidualMass = G4NucleiProperties::GetNuclearMass(ResidualA, ResidualZ);
00119 G4double FragmentMass = fragment->GetGroundStateMass();
00120 CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
00121
00122 MaximalKineticEnergy = CalcMaximalKineticEnergy(FragmentMass + ExEnergy);
00123
00124
00125
00126
00127
00128
00129
00130
00131 G4double limit;
00132 if (OPTxs==0 || (OPTxs!=0 && useSICB))
00133 limit= CoulombBarrier;
00134 else limit=0.;
00135 limit = std::max(limit, FragmentMass - ResidualMass - EvaporatedMass);
00136
00137
00138
00139 if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0;
00140 else {
00141
00142 EmissionProbability = theEvaporationProbabilityPtr->
00143 EmissionProbability(*fragment, MaximalKineticEnergy);
00144 }
00145 }
00146
00147
00148 return EmissionProbability;
00149 }
00150
00151 G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus)
00152 {
00153
00154
00155
00156
00157
00158
00159 G4double EvaporatedEnergy = GetKineticEnergy(theNucleus) + EvaporatedMass;
00160
00161 G4ThreeVector momentum(IsotropicVector
00162 (std::sqrt((EvaporatedEnergy - EvaporatedMass)*
00163 (EvaporatedEnergy + EvaporatedMass))));
00164
00165 G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
00166 G4LorentzVector ResidualMomentum = theNucleus.GetMomentum();
00167 EvaporatedMomentum.boost(ResidualMomentum.boostVector());
00168
00169 G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum);
00170 ResidualMomentum -= EvaporatedMomentum;
00171
00172 G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum);
00173
00174 G4FragmentVector * theResult = new G4FragmentVector;
00175
00176 #ifdef debug
00177 G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
00178 G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
00179 if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
00180 G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
00181 G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV Final :"
00182 <<Efinal/MeV << " MeV Delta: " << (Efinal-theNucleus.GetMomentum().e())/MeV
00183 << " MeV" << G4endl;
00184 }
00185 if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
00186 std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
00187 std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
00188 G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
00189 G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV Final :"
00190 <<Pfinal/MeV << " MeV Delta: " << Pfinal-theNucleus.GetMomentum().vect()
00191 << " MeV" << G4endl;
00192
00193 }
00194 #endif
00195 theResult->push_back(EvaporatedFragment);
00196 theResult->push_back(ResidualFragment);
00197 return theResult;
00198 }
00199
00201
00202 G4double G4EvaporationChannel::CalcMaximalKineticEnergy(G4double NucleusTotalE)
00203 {
00204
00205 G4double Tmax =
00206 ((NucleusTotalE-ResidualMass)*(NucleusTotalE+ResidualMass) + EvaporatedMass*EvaporatedMass)
00207 /(2.0*NucleusTotalE) - EvaporatedMass;
00208
00209
00210
00211
00212
00213
00214 if(OPTxs==0)
00215 Tmax=Tmax- CoulombBarrier;
00216
00217 return Tmax;
00218 }
00219
00221
00222 G4double G4EvaporationChannel::GetKineticEnergy(const G4Fragment & aFragment)
00223 {
00224
00225 if (OPTxs==0) {
00226
00227
00228
00229
00230
00231
00232 if (MaximalKineticEnergy < 0.0) {
00233 throw G4HadronicException(__FILE__, __LINE__,
00234 "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0");
00235 }
00236 G4double Rb = 4.0*theLevelDensityPtr->
00237 LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)*
00238 MaximalKineticEnergy;
00239 G4double RbSqrt = std::sqrt(Rb);
00240 G4double PEX1 = 0.0;
00241 if (RbSqrt < 160.0) PEX1 = std::exp(-RbSqrt);
00242 G4double Rk = 0.0;
00243 G4double FRk = 0.0;
00244 do {
00245 G4double RandNumber = G4UniformRand();
00246 Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
00247 G4double Q1 = 1.0;
00248 G4double Q2 = 1.0;
00249 if (theZ == 0) {
00250 G4double Beta = (2.12/G4Pow::GetInstance()->Z23(ResidualA) - 0.05)*MeV/
00251 (0.76 + 2.2/G4Pow::GetInstance()->Z13(ResidualA));
00252 Q1 = 1.0 + Beta/(MaximalKineticEnergy);
00253 Q2 = Q1*std::sqrt(Q1);
00254 }
00255
00256 FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
00257
00258 } while (FRk < G4UniformRand());
00259
00260 G4double result = MaximalKineticEnergy * (1.0-Rk*Rk) + CoulombBarrier;
00261 return result;
00262 } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) {
00263
00264
00265 G4double V = 0;
00266 if(useSICB) { V= CoulombBarrier; }
00267
00268 V = std::max(V, aFragment.GetGroundStateMass()-EvaporatedMass-ResidualMass);
00269
00270 G4double Tmax=MaximalKineticEnergy;
00271 G4double T(0.0);
00272 G4double NormalizedProbability(1.0);
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 do
00298 {
00299 T=V+G4UniformRand()*(Tmax-V);
00300 NormalizedProbability =
00301 theEvaporationProbabilityPtr->ProbabilityDistributionFunction(aFragment,T)/
00302 GetEmissionProbability(const_cast<G4Fragment*>(&aFragment));
00303
00304 }
00305 while (G4UniformRand() > NormalizedProbability);
00306
00307 return T;
00308 } else{
00309 std::ostringstream errOs;
00310 errOs << "Bad option for energy sampling in evaporation" <<G4endl;
00311 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00312 }
00313 }
00314
00315 G4ThreeVector G4EvaporationChannel::IsotropicVector(G4double Magnitude)
00316
00317
00318 {
00319 G4double CosTheta = 1.0 - 2.0*G4UniformRand();
00320 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
00321 G4double Phi = twopi*G4UniformRand();
00322 G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
00323 Magnitude*std::sin(Phi)*SinTheta,
00324 Magnitude*CosTheta);
00325 return Vector;
00326 }
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336