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 <set>
00038
00039 #include "G4FermiConfigurationList.hh"
00040 #include "G4FermiFragmentsPool.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "Randomize.hh"
00043 #include "G4Pow.hh"
00044
00045 const G4double G4FermiConfigurationList::Kappa = 6.0;
00046 const G4double G4FermiConfigurationList::r0 = 1.3*CLHEP::fermi;
00047
00048 G4FermiConfigurationList::G4FermiConfigurationList()
00049 {
00050 thePool = G4FermiFragmentsPool::Instance();
00051 g4pow = G4Pow::GetInstance();
00052 Coef = 0.6*(CLHEP::elm_coupling/r0)/g4pow->Z13(1+G4int(Kappa));
00053 ConstCoeff = g4pow->powN(r0/hbarc,3)*Kappa*std::sqrt(2.0/pi)/3.0;
00054
00055
00056 nmax = 50;
00057 NormalizedWeights.resize(nmax,0.0);
00058 }
00059
00060 G4FermiConfigurationList::~G4FermiConfigurationList()
00061 {}
00062
00063 G4double
00064 G4FermiConfigurationList::CoulombBarrier(
00065 const std::vector<const G4VFermiFragment*>& conf)
00066 {
00067
00068 G4int SumA = 0;
00069 G4int SumZ = 0;
00070 G4double CoulombEnergy = 0.;
00071 size_t nn = conf.size();
00072 for (size_t i=0; i<nn; ++i) {
00073 G4int z = conf[i]->GetZ();
00074 G4int a = conf[i]->GetA();
00075 CoulombEnergy += G4double(z*z)/g4pow->Z13(a);
00076 SumA += a;
00077 SumZ += z;
00078 }
00079 CoulombEnergy -= SumZ*SumZ/g4pow->Z13(SumA);
00080 return -Coef * CoulombEnergy;
00081 }
00082
00083 G4double
00084 G4FermiConfigurationList::DecayProbability(G4int A, G4double TotalE,
00085 G4FermiConfiguration* conf)
00086
00087 {
00088
00089
00090
00091 G4double KineticEnergy = TotalE;
00092 G4double ProdMass = 1.0;
00093 G4double SumMass = 0.0;
00094 G4double S_n = 1.0;
00095 std::set<G4int> combSet;
00096 std::multiset<G4int> combmSet;
00097
00098 const std::vector<const G4VFermiFragment*> flist =
00099 conf->GetFragmentList();
00100
00101
00102 G4int K = flist.size();
00103
00104 for (G4int i=0; i<K; ++i) {
00105 G4int a = flist[i]->GetA();
00106 combSet.insert(a);
00107 combmSet.insert(a);
00108 G4double mass = flist[i]->GetFragmentMass();
00109 ProdMass *= mass;
00110 SumMass += mass;
00111
00112 S_n *= flist[i]->GetPolarization();
00113 KineticEnergy -= mass + flist[i]->GetExcitationEnergy();
00114 }
00115
00116
00117 KineticEnergy -= CoulombBarrier(flist);
00118 if (KineticEnergy <= 0.0) { return 0.0; }
00119
00120 G4double MassFactor = ProdMass/SumMass;
00121 MassFactor *= std::sqrt(MassFactor);
00122
00123
00124 G4double Coeff = g4pow->powN(ConstCoeff*A, K-1);
00125
00126
00127
00128
00129 G4double Gamma = 1.0;
00130
00131
00132
00133
00134
00135
00136
00137
00138 if ((K-1)%2 != 1)
00139
00140 {
00141 G4double arg = 3.0*(K-1)/2.0 - 1.0;
00142 while (arg > 1.1)
00143 {
00144 Gamma *= arg;
00145 arg--;
00146 }
00147 }
00148 else {
00149 G4double n= 3.0*K/2.0-2.0;
00150 G4double arg2=2*n-1;
00151 while (arg2>1.1)
00152 {
00153 Gamma *= arg2;
00154 arg2 -= 2;
00155 }
00156 Gamma *= std::sqrt(pi)/g4pow->powZ(2,n);
00157 }
00158
00159
00160
00161
00162 G4double G_n = 1.0;
00163 for (std::set<G4int>::iterator itr = combSet.begin();
00164 itr != combSet.end(); ++itr)
00165 {
00166 for (G4int ni = combmSet.count(*itr); ni > 1; ni--) { G_n *= ni; }
00167 }
00168
00169 G4double Weight = Coeff * MassFactor * (S_n / G_n) / Gamma;
00170 Weight *= std::sqrt(g4pow->powN(KineticEnergy,3*(K-1)))/KineticEnergy;
00171
00172 return Weight;
00173 }
00174
00175 G4FragmentVector*
00176 G4FermiConfigurationList::GetFragments(const G4Fragment& theNucleus)
00177 {
00178
00179 G4double M = theNucleus.GetMomentum().m();
00180 const std::vector<const G4VFermiFragment*>* conf =
00181 SelectConfiguration(theNucleus.GetZ_asInt(),
00182 theNucleus.GetA_asInt(), M);
00183
00184
00185 G4FragmentVector* theResult = new G4FragmentVector();
00186 size_t nn = conf->size();
00187 if(1 >= nn) {
00188 theResult->push_back(new G4Fragment(theNucleus));
00189 delete conf;
00190 return theResult;
00191 }
00192
00193 G4ThreeVector boostVector = theNucleus.GetMomentum().boostVector();
00194 std::vector<G4double> mr;
00195 mr.reserve(nn);
00196 for(size_t i=0; i<nn; ++i) {
00197 mr.push_back( (*conf)[i]->GetTotalEnergy() );
00198 }
00199 std::vector<G4LorentzVector*>* mom = thePhaseSpace.Decay(M,mr);
00200 if(!mom) {
00201 delete conf;
00202 return theResult;
00203 }
00204
00205 size_t nmom = mom->size();
00206
00207
00208 if(0 < nmom) {
00209 for (size_t j=0; j<nmom; ++j) {
00210 G4LorentzVector* FourMomentum = (*mom)[j];
00211
00212
00213 FourMomentum->boost(boostVector);
00214
00215 G4FragmentVector* fragment = (*conf)[j]->GetFragment(*FourMomentum);
00216
00217 size_t nfrag = fragment->size();
00218 for (size_t k=0; k<nfrag; ++k) { theResult->push_back((*fragment)[k]); }
00219 delete fragment;
00220 delete (*mom)[j];
00221 }
00222 }
00223
00224 delete mom;
00225 delete conf;
00226 return theResult;
00227 }
00228
00229 const std::vector<const G4VFermiFragment*>*
00230 G4FermiConfigurationList::SelectConfiguration(G4int Z, G4int A, G4double mass)
00231 {
00232 std::vector<const G4VFermiFragment*>* res =
00233 new std::vector<const G4VFermiFragment*>;
00234 const std::vector<G4FermiConfiguration*>* conflist =
00235 thePool->GetConfigurationList(Z, A, mass);
00236 if(!conflist) { return res; }
00237 size_t nn = conflist->size();
00238 if(0 < nn) {
00239 size_t idx = 0;
00240 if(1 < nn) {
00241 if(nn > nmax) {
00242 nmax = nn;
00243 NormalizedWeights.resize(nmax,0.0);
00244 }
00245 G4double prob = 0.0;
00246 for(size_t i=0; i<nn; ++i) {
00247 prob += DecayProbability(A, mass, (*conflist)[i]);
00248 NormalizedWeights[i] = prob;
00249 }
00250 prob *= G4UniformRand();
00251 for(idx=0; idx<nn; ++idx) {
00252 if(NormalizedWeights[idx] >= prob) { break; }
00253 }
00254 }
00255 const std::vector<const G4VFermiFragment*> flist =
00256 (*conflist)[idx]->GetFragmentList();
00257 size_t nf = flist.size();
00258 for(size_t i=0; i<nf; ++i) { res->push_back(flist[i]); }
00259
00260
00261 }
00262 delete conflist;
00263 return res;
00264 }