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
00046
00047
00048
00049 #include "G4Evaporation.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4EvaporationFactory.hh"
00052 #include "G4EvaporationGEMFactory.hh"
00053 #include "G4EvaporationDefaultGEMFactory.hh"
00054 #include "G4HadronicException.hh"
00055 #include "G4NistManager.hh"
00056 #include "G4FermiFragmentsPool.hh"
00057 #include "G4PhotonEvaporation.hh"
00058
00059 G4Evaporation::G4Evaporation()
00060 : theChannels(0),nChannels(0)
00061 {
00062 SetPhotonEvaporation(new G4PhotonEvaporation());
00063 theChannelFactory = new G4EvaporationDefaultGEMFactory(thePhotonEvaporation);
00064 SetParameters();
00065 InitialiseEvaporation();
00066 }
00067
00068 G4Evaporation::G4Evaporation(G4VEvaporationChannel* photoEvaporation)
00069 : theChannels(0),nChannels(0)
00070 {
00071 if(photoEvaporation) { SetPhotonEvaporation(photoEvaporation); }
00072 else { SetPhotonEvaporation(new G4PhotonEvaporation()); }
00073
00074 theChannelFactory = new G4EvaporationDefaultGEMFactory(thePhotonEvaporation);
00075 SetParameters();
00076 InitialiseEvaporation();
00077 }
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099 G4Evaporation::~G4Evaporation()
00100 {
00101 CleanChannels();
00102 delete thePhotonEvaporation;
00103 delete theChannelFactory;
00104 }
00105
00106 void G4Evaporation::CleanChannels()
00107 {
00108 for (size_t i=1; i<nChannels; ++i) { delete (*theChannels)[i]; }
00109 delete theChannels;
00110 }
00111
00112 void G4Evaporation::SetParameters()
00113 {
00114 nist = G4NistManager::Instance();
00115 minExcitation = CLHEP::keV;
00116 maxZforFBU = G4FermiFragmentsPool::Instance()->GetMaxZ();
00117 maxAforFBU = G4FermiFragmentsPool::Instance()->GetMaxA();
00118 probabilities.reserve(68);
00119 }
00120
00121 void G4Evaporation::InitialiseEvaporation()
00122 {
00123 CleanChannels();
00124 theChannels = theChannelFactory->GetChannel();
00125 nChannels = theChannels->size();
00126 probabilities.resize(nChannels, 0.0);
00127 Initialise();
00128 }
00129
00130 void G4Evaporation::Initialise()
00131 {
00132 for(size_t i=0; i<nChannels; ++i) {
00133 (*theChannels)[i]->SetOPTxs(OPTxs);
00134 (*theChannels)[i]->UseSICB(useSICB);
00135 }
00136 }
00137
00138 void G4Evaporation::SetDefaultChannel()
00139 {
00140 delete theChannelFactory;
00141 theChannelFactory = new G4EvaporationFactory(thePhotonEvaporation);
00142 InitialiseEvaporation();
00143 }
00144
00145 void G4Evaporation::SetGEMChannel()
00146 {
00147 delete theChannelFactory;
00148 theChannelFactory = new G4EvaporationGEMFactory(thePhotonEvaporation);
00149 InitialiseEvaporation();
00150 }
00151
00152 void G4Evaporation::SetCombinedChannel()
00153 {
00154 delete theChannelFactory;
00155 theChannelFactory = new G4EvaporationDefaultGEMFactory(thePhotonEvaporation);
00156 InitialiseEvaporation();
00157 }
00158
00159 void G4Evaporation::SetPhotonEvaporation(G4VEvaporationChannel* ptr)
00160 {
00161 if(ptr) { G4VEvaporation::SetPhotonEvaporation(ptr); }
00162 if(0 < nChannels) { (*theChannels)[0] = ptr; }
00163 }
00164
00165 G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
00166 {
00167 G4FragmentVector * theResult = new G4FragmentVector;
00168 G4FragmentVector * theTempResult;
00169 const G4double Elimit = 3*MeV;
00170
00171
00172 G4Fragment* theResidualNucleus = new G4Fragment(theNucleus);
00173
00174 G4double totprob, prob, oldprob = 0.0;
00175 size_t maxchannel, i;
00176
00177 G4int Amax = theResidualNucleus->GetA_asInt();
00178
00179
00180
00181 for(G4int ia=0; ia<Amax; ++ia) {
00182
00183
00184 G4int Z = theResidualNucleus->GetZ_asInt();
00185 G4int A = theResidualNucleus->GetA_asInt();
00186
00187
00188 if(maxZforFBU > Z && maxAforFBU >= A) {
00189 theResult->push_back(theResidualNucleus);
00190 return theResult;
00191 }
00192
00193
00194 G4double abun = nist->GetIsotopeAbundance(Z, A);
00195
00196
00197
00198
00199
00200
00201
00202 G4double Eex = theResidualNucleus->GetExcitationEnergy();
00203 if(Eex <= minExcitation && abun > 0.0) {
00204 theResult->push_back(theResidualNucleus);
00205 return theResult;
00206 }
00207
00208 totprob = 0.0;
00209 maxchannel = nChannels;
00210
00211
00212
00213
00214
00215 for(i=0; i<nChannels; ++i) {
00216 prob = (*theChannels)[i]->GetEmissionProbability(theResidualNucleus);
00217
00218
00219 totprob += prob;
00220 probabilities[i] = totprob;
00221
00222 if(i>=8) {
00223 if(prob <= totprob*1.e-8 && oldprob <= totprob*1.e-8) {
00224 maxchannel = i+1;
00225 break;
00226 }
00227 }
00228 oldprob = prob;
00229
00230 if(7 == i && Eex > Elimit*A) {
00231 maxchannel = 8;
00232 break;
00233 }
00234 }
00235
00236
00237
00238 if(0.0 < totprob && probabilities[0] == totprob) {
00239
00240 theTempResult = (*theChannels)[0]->BreakUpFragment(theResidualNucleus);
00241 if(theTempResult) {
00242 size_t nsec = theTempResult->size();
00243 for(size_t j=0; j<nsec; ++j) {
00244 theResult->push_back((*theTempResult)[j]);
00245 }
00246 delete theTempResult;
00247 }
00248 totprob = 0.0;
00249 }
00250
00251
00252 if(0.0 == totprob) {
00253
00254
00255 if(0.0 == abun && Z < 20) {
00256
00257 theTempResult = unstableBreakUp.BreakUpFragment(theResidualNucleus);
00258 if(theTempResult) {
00259 size_t nsec = theTempResult->size();
00260 for(size_t j=0; j<nsec; ++j) {
00261 theResult->push_back((*theTempResult)[j]);
00262 }
00263 delete theTempResult;
00264 }
00265 }
00266
00267
00268 theResult->push_back(theResidualNucleus);
00269 return theResult;
00270 }
00271
00272
00273
00274 totprob *= G4UniformRand();
00275
00276 for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
00277
00278
00279 if(i >= nChannels) { i = nChannels - 1; }
00280
00281
00282
00283 if(0 == i) {
00284
00285 G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus);
00286 if(gamma) { theResult->push_back(gamma); }
00287
00288
00289 } else if(1 == i) {
00290
00291 theTempResult = (*theChannels)[1]->BreakUp(*theResidualNucleus);
00292 if(theTempResult) {
00293 size_t nsec = theTempResult->size();
00294 G4bool deletePrimary = true;
00295 for(size_t j=0; j<nsec; ++j) {
00296 if(theResidualNucleus == (*theTempResult)[j]) { deletePrimary = false; }
00297 theResult->push_back((*theTempResult)[j]);
00298 }
00299 if(deletePrimary) { delete theResidualNucleus; }
00300 delete theTempResult;
00301 return theResult;
00302 }
00303
00304
00305 } else {
00306
00307 theTempResult = (*theChannels)[i]->BreakUp(*theResidualNucleus);
00308 if(theTempResult) {
00309 size_t nsec = theTempResult->size();
00310 if(nsec > 0) {
00311 --nsec;
00312 for(size_t j=0; j<nsec; ++j) {
00313 theResult->push_back((*theTempResult)[j]);
00314 }
00315
00316
00317 if(theResidualNucleus != (*theTempResult)[nsec] ) {
00318 delete theResidualNucleus;
00319 theResidualNucleus = (*theTempResult)[nsec];
00320 }
00321 }
00322 delete theTempResult;
00323 }
00324 }
00325 }
00326
00327
00328 theResult->push_back(theResidualNucleus);
00329
00330 return theResult;
00331 }
00332