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
00042
00043
00044
00045 #include "G4PreCompoundEmission.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4Pow.hh"
00049 #include "Randomize.hh"
00050 #include "G4PreCompoundParameters.hh"
00051 #include "G4PreCompoundEmissionFactory.hh"
00052 #include "G4HETCEmissionFactory.hh"
00053 #include "G4HadronicException.hh"
00054
00055 G4PreCompoundEmission::G4PreCompoundEmission()
00056 {
00057 theFragmentsFactory = new G4PreCompoundEmissionFactory();
00058 theFragmentsVector =
00059 new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
00060 g4pow = G4Pow::GetInstance();
00061 theParameters = G4PreCompoundParameters::GetAddress();
00062 }
00063
00064 G4PreCompoundEmission::~G4PreCompoundEmission()
00065 {
00066 if (theFragmentsFactory) { delete theFragmentsFactory; }
00067 if (theFragmentsVector) { delete theFragmentsVector; }
00068 }
00069
00070 void G4PreCompoundEmission::SetDefaultModel()
00071 {
00072 if (theFragmentsFactory) { delete theFragmentsFactory; }
00073 theFragmentsFactory = new G4PreCompoundEmissionFactory();
00074 if (theFragmentsVector)
00075 {
00076 theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector());
00077 }
00078 else
00079 {
00080 theFragmentsVector =
00081 new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
00082 }
00083 return;
00084 }
00085
00086 void G4PreCompoundEmission::SetHETCModel()
00087 {
00088 if (theFragmentsFactory) delete theFragmentsFactory;
00089 theFragmentsFactory = new G4HETCEmissionFactory();
00090 if (theFragmentsVector)
00091 {
00092 theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector());
00093 }
00094 else
00095 {
00096 theFragmentsVector =
00097 new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
00098 }
00099 return;
00100 }
00101
00102 G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)
00103 {
00104
00105 G4VPreCompoundFragment * thePreFragment = theFragmentsVector->ChooseFragment();
00106 if (thePreFragment == 0)
00107 {
00108 G4cout << "G4PreCompoundEmission::PerformEmission : I couldn't choose a fragment\n"
00109 << "while trying to de-excite\n"
00110 << aFragment << G4endl;
00111 throw G4HadronicException(__FILE__, __LINE__, "");
00112 }
00113
00114
00115
00116
00117
00118 G4double kinEnergyOfEmittedFragment = thePreFragment->GetKineticEnergy(aFragment);
00119
00120
00121
00122
00123
00124 if(kinEnergyOfEmittedFragment < 0.0) { kinEnergyOfEmittedFragment = 0.0; }
00125
00126
00127 AngularDistribution(thePreFragment,aFragment,kinEnergyOfEmittedFragment);
00128
00129
00130 G4double EmittedMass = thePreFragment->GetNuclearMass();
00131
00132
00133 G4LorentzVector Emitted4Momentum(theFinalMomentum,
00134 EmittedMass + kinEnergyOfEmittedFragment);
00135
00136
00137 G4LorentzVector Rest4Momentum = aFragment.GetMomentum();
00138 Emitted4Momentum.boost(Rest4Momentum.boostVector());
00139
00140
00141 thePreFragment->SetMomentum(Emitted4Momentum);
00142
00143
00144
00145
00146 Rest4Momentum -= Emitted4Momentum;
00147
00148
00149
00150
00151
00152 aFragment.SetZandA_asInt(thePreFragment->GetRestZ(),
00153 thePreFragment->GetRestA());
00154
00155
00156 aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
00157 thePreFragment->GetA());
00158
00159 aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()-
00160 thePreFragment->GetZ());
00161
00162
00163
00164 aFragment.SetMomentum(Rest4Momentum);
00165
00166
00167 G4ReactionProduct * MyRP = thePreFragment->GetReactionProduct();
00168
00169
00170
00171
00172
00173 return MyRP;
00174 }
00175
00176 void
00177 G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment* thePreFragment,
00178 const G4Fragment& aFragment,
00179 G4double ekin)
00180 {
00181 G4int p = aFragment.GetNumberOfParticles();
00182 G4int h = aFragment.GetNumberOfHoles();
00183 G4double U = aFragment.GetExcitationEnergy();
00184
00185
00186 G4double Bemission = thePreFragment->GetBindingEnergy();
00187
00188
00189 G4double Ef = theParameters->GetFermiEnergy();
00190
00191
00192
00193
00194
00195 G4double gg = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
00196
00197
00198 G4double Eav = 2*p*(p+1)/((p+h)*gg);
00199
00200
00201 G4double Uf = std::max(U - (p - h)*Ef , 0.0);
00202
00203
00204 G4double w_num = rho(p+1, h, gg, Uf, Ef);
00205 G4double w_den = rho(p, h, gg, Uf, Ef);
00206 if (w_num > 0.0 && w_den > 0.0)
00207 {
00208 Eav *= (w_num/w_den);
00209 Eav += - Uf/(p+h) + Ef;
00210 }
00211 else
00212 {
00213 Eav = Ef;
00214 }
00215
00216
00217
00218 G4double an = 0.0;
00219 G4double Eeff = ekin + Bemission + Ef;
00220 if(ekin > DBL_MIN && Eeff > DBL_MIN) {
00221
00222 G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV));
00223
00224
00225
00226
00227
00228 G4double ProjEnergy = aFragment.GetExcitationEnergy();
00229
00230 an = 3*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav);
00231
00232 G4int ne = aFragment.GetNumberOfExcitons() - 1;
00233 if ( ne > 1 ) { an /= (G4double)ne; }
00234
00235
00236 if ( an > 10. ) { an = 10.; }
00237 }
00238
00239
00240 G4double random = G4UniformRand();
00241 G4double cost;
00242
00243 if(an < 0.1) { cost = 1. - 2*random; }
00244 else {
00245 G4double exp2an = std::exp(-2*an);
00246 cost = 1. + std::log(1-random*(1-exp2an))/an;
00247 if(cost > 1.) { cost = 1.; }
00248 else if(cost < -1.) {cost = -1.; }
00249 }
00250
00251 G4double phi = CLHEP::twopi*G4UniformRand();
00252
00253
00254 G4double pmag = std::sqrt(ekin*(ekin + 2.0*thePreFragment->GetNuclearMass()));
00255
00256 G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
00257
00258 theFinalMomentum.set(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost);
00259
00260
00261 G4ThreeVector theIncidentDirection = aFragment.GetMomentum().vect().unit();
00262 theFinalMomentum.rotateUz(theIncidentDirection);
00263 }
00264
00265 G4double G4PreCompoundEmission::rho(G4int p, G4int h, G4double gg,
00266 G4double E, G4double Ef) const
00267 {
00268
00269 G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*gg);
00270
00271
00272 if ( E - Aph < 0.0) { return 0.0; }
00273
00274 G4double logConst = (p+h)*std::log(gg)
00275 - g4pow->logfactorial(p+h-1) - g4pow->logfactorial(p) - g4pow->logfactorial(h);
00276
00277
00278
00279 G4double t1=1;
00280 G4double t2=1;
00281 G4double logt3 = (p+h-1) * std::log(E-Aph) + logConst;
00282 const G4double logmax = 200.;
00283 if(logt3 > logmax) { logt3 = logmax; }
00284 G4double tot = std::exp( logt3 );
00285
00286
00287
00288 G4double Eeff = E - Aph;
00289 for(G4int j=1; j<=h; ++j)
00290 {
00291 Eeff -= Ef;
00292 if(Eeff < 0.0) { break; }
00293 t1 *= -1.;
00294 t2 *= (G4double)(h+1-j)/(G4double)j;
00295 logt3 = (p+h-1) * std::log( Eeff) + logConst;
00296 if(logt3 > logmax) { logt3 = logmax; }
00297 tot += t1*t2*std::exp(logt3);
00298 }
00299
00300 return tot;
00301 }