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
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 #include <list>
00066
00067 #include "G4ExcitationHandler.hh"
00068 #include "G4SystemOfUnits.hh"
00069 #include "G4LorentzVector.hh"
00070 #include "G4NistManager.hh"
00071 #include "G4ParticleTable.hh"
00072 #include "G4ParticleTypes.hh"
00073
00074 #include "G4VMultiFragmentation.hh"
00075 #include "G4VFermiBreakUp.hh"
00076 #include "G4VFermiFragment.hh"
00077
00078 #include "G4VEvaporation.hh"
00079 #include "G4VEvaporationChannel.hh"
00080 #include "G4VPhotonEvaporation.hh"
00081 #include "G4Evaporation.hh"
00082 #include "G4StatMF.hh"
00083 #include "G4PhotonEvaporation.hh"
00084 #include "G4FermiBreakUp.hh"
00085 #include "G4FermiFragmentsPool.hh"
00086
00087 G4ExcitationHandler::G4ExcitationHandler():
00088 maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4*GeV),
00089 minExcitation(keV),OPTxs(3),useSICB(false),isEvapLocal(true)
00090 {
00091 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
00092
00093 theMultiFragmentation = new G4StatMF;
00094 theFermiModel = new G4FermiBreakUp;
00095 thePhotonEvaporation = new G4PhotonEvaporation;
00096 theEvaporation = new G4Evaporation(thePhotonEvaporation);
00097 thePool = G4FermiFragmentsPool::Instance();
00098 SetParameters();
00099 }
00100
00101 G4ExcitationHandler::~G4ExcitationHandler()
00102 {
00103 if(isEvapLocal) { delete theEvaporation; }
00104 delete theMultiFragmentation;
00105 delete theFermiModel;
00106 }
00107
00108 void G4ExcitationHandler::SetParameters()
00109 {
00110
00111 theEvaporation->SetOPTxs(OPTxs);
00112
00113 theEvaporation->UseSICB(useSICB);
00114 theEvaporation->Initialise();
00115 }
00116
00117 G4ReactionProductVector *
00118 G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
00119 {
00120
00121
00122
00123 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
00124
00125 G4FragmentVector * theTempResult = 0;
00126 std::list<G4Fragment*> theEvapList;
00127 std::list<G4Fragment*> thePhotoEvapList;
00128 std::list<G4Fragment*> theResults;
00129
00130
00131
00132
00133 G4double exEnergy = theInitialState.GetExcitationEnergy();
00134 G4int A = theInitialState.GetA_asInt();
00135 G4int Z = theInitialState.GetZ_asInt();
00136
00137 G4NistManager* nist = G4NistManager::Instance();
00138
00139
00140 if (A <= 1)
00141 {
00142 theResults.push_back( theInitialStatePtr );
00143 }
00144
00145 else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0)
00146 {
00147 theResults.push_back( theInitialStatePtr );
00148 }
00149 else
00150 {
00151
00152
00153
00154
00155
00156 if((A<maxAForFermiBreakUp && Z<maxZForFermiBreakUp)
00157 || exEnergy <= minEForMultiFrag*A)
00158 {
00159 theEvapList.push_back(theInitialStatePtr);
00160 }
00161 else
00162 {
00163 theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
00164 if(!theTempResult) { theEvapList.push_back(theInitialStatePtr); }
00165 else {
00166 size_t nsec = theTempResult->size();
00167 if(0 == nsec) { theEvapList.push_back(theInitialStatePtr); }
00168 else {
00169
00170
00171 G4bool deletePrimary = true;
00172 G4FragmentVector::iterator j;
00173 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
00174 if((*j) == theInitialStatePtr) { deletePrimary = false; }
00175 A = (*j)->GetA_asInt();
00176
00177
00178 if(A <= 1) { theResults.push_back(*j); }
00179
00180
00181 else {
00182 G4double exEnergy1 = (*j)->GetExcitationEnergy();
00183
00184
00185 if(exEnergy1 < minExcitation) {
00186 Z = (*j)->GetZ_asInt();
00187 if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
00188 theResults.push_back(*j);
00189
00190 } else {
00191
00192
00193 const G4VFermiFragment* ffrag = thePool->GetFragment(Z, A);
00194 if(ffrag) {
00195 if(ffrag->IsStable()) { theResults.push_back(*j); }
00196 else { theEvapList.push_back(*j); }
00197
00198
00199 } else {
00200 theEvapList.push_back(*j);
00201 }
00202 }
00203
00204
00205 } else { theEvapList.push_back(*j); }
00206 }
00207 }
00208 if( deletePrimary ) { delete theInitialStatePtr; }
00209 }
00210 delete theTempResult;
00211 }
00212 }
00213 }
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223 std::list<G4Fragment*>::iterator iList;
00224 for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList)
00225 {
00226
00227
00228 A = (*iList)->GetA_asInt();
00229 Z = (*iList)->GetZ_asInt();
00230
00231
00232 G4bool wasFBU = false;
00233 if (A < maxAForFermiBreakUp && Z < maxZForFermiBreakUp)
00234 {
00235 theTempResult = theFermiModel->BreakItUp(*(*iList));
00236 wasFBU = true;
00237
00238 if(1 == theTempResult->size()) {
00239 delete theTempResult;
00240 theTempResult = theEvaporation->BreakItUp(*(*iList));
00241 }
00242 }
00243 else
00244 {
00245 theTempResult = theEvaporation->BreakItUp(*(*iList));
00246 }
00247
00248 G4bool deletePrimary = true;
00249 size_t nsec = theTempResult->size();
00250
00251
00252
00253 if ( nsec > 0 ) {
00254 G4FragmentVector::iterator j;
00255 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
00256 if((*j) == (*iList)) { deletePrimary = false; }
00257
00258
00259 A = (*j)->GetA_asInt();
00260 exEnergy = (*j)->GetExcitationEnergy();
00261
00262 if(A <= 1) { theResults.push_back(*j); }
00263
00264
00265 else if(1 == nsec) {
00266 if(exEnergy < minExcitation) { theResults.push_back(*j); }
00267 else { thePhotoEvapList.push_back(*j); }
00268
00269 } else {
00270
00271
00272 if(exEnergy < minExcitation) {
00273 Z = (*j)->GetZ_asInt();
00274
00275
00276 if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
00277 theResults.push_back(*j);
00278
00279 } else {
00280 const G4VFermiFragment* ffrag = thePool->GetFragment(Z, A);
00281
00282
00283 if(ffrag) {
00284 if(ffrag->IsStable()) { theResults.push_back(*j); }
00285 else { theEvapList.push_back(*j); }
00286
00287
00288 } else {
00289 theEvapList.push_back(*j);
00290 }
00291 }
00292
00293
00294 } else if (wasFBU) {
00295 thePhotoEvapList.push_back(*j);
00296 } else {
00297 theEvapList.push_back(*j);
00298 }
00299 }
00300 }
00301 }
00302 if( deletePrimary ) { delete (*iList); }
00303 delete theTempResult;
00304 }
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315 for(iList = thePhotoEvapList.begin(); iList != thePhotoEvapList.end(); ++iList)
00316 {
00317
00318
00319 exEnergy = (*iList)->GetExcitationEnergy();
00320
00321
00322 if(exEnergy >= minExcitation) {
00323 theTempResult = thePhotonEvaporation->BreakUpFragment(*iList);
00324 size_t nsec = theTempResult->size();
00325
00326
00327
00328 if (nsec > 0)
00329 {
00330 G4FragmentVector::iterator j;
00331 for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
00332 {
00333 theResults.push_back(*j);
00334 }
00335 }
00336 delete theTempResult;
00337 }
00338
00339
00340 theResults.push_back(*iList);
00341
00342 }
00343
00344
00345
00346
00347
00348 G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
00349
00350
00351
00352 theReactionProductVector->reserve( theResults.size() );
00353
00354 G4int theFragmentA, theFragmentZ;
00355
00356 std::list<G4Fragment*>::iterator i;
00357 for (i = theResults.begin(); i != theResults.end(); ++i)
00358 {
00359 theFragmentA = (*i)->GetA_asInt();
00360 theFragmentZ = (*i)->GetZ_asInt();
00361 G4ParticleDefinition* theKindOfFragment = 0;
00362 if (theFragmentA == 0) {
00363 theKindOfFragment = (*i)->GetParticleDefinition();
00364 } else if (theFragmentA == 1 && theFragmentZ == 0) {
00365 theKindOfFragment = G4Neutron::NeutronDefinition();
00366 } else if (theFragmentA == 1 && theFragmentZ == 1) {
00367 theKindOfFragment = G4Proton::ProtonDefinition();
00368 } else if (theFragmentA == 2 && theFragmentZ == 1) {
00369 theKindOfFragment = G4Deuteron::DeuteronDefinition();
00370 } else if (theFragmentA == 3 && theFragmentZ == 1) {
00371 theKindOfFragment = G4Triton::TritonDefinition();
00372 } else if (theFragmentA == 3 && theFragmentZ == 2) {
00373 theKindOfFragment = G4He3::He3Definition();
00374 } else if (theFragmentA == 4 && theFragmentZ == 2) {
00375 theKindOfFragment = G4Alpha::AlphaDefinition();;
00376 } else {
00377 theKindOfFragment =
00378 theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0);
00379 }
00380 if (theKindOfFragment != 0)
00381 {
00382 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
00383 theNew->SetMomentum((*i)->GetMomentum().vect());
00384 theNew->SetTotalEnergy((*i)->GetMomentum().e());
00385 theNew->SetFormationTime((*i)->GetCreationTime());
00386 theReactionProductVector->push_back(theNew);
00387 }
00388 delete (*i);
00389 }
00390
00391 return theReactionProductVector;
00392 }
00393
00394 void G4ExcitationHandler::SetEvaporation(G4VEvaporation* ptr)
00395 {
00396 if(ptr && ptr != theEvaporation) {
00397 delete theEvaporation;
00398 theEvaporation = ptr;
00399 thePhotonEvaporation = ptr->GetPhotonEvaporation();
00400 SetParameters();
00401 isEvapLocal = false;
00402 }
00403 }
00404
00405 void
00406 G4ExcitationHandler::SetMultiFragmentation(G4VMultiFragmentation* ptr)
00407 {
00408 if(ptr && ptr != theMultiFragmentation) {
00409 delete theMultiFragmentation;
00410 theMultiFragmentation = ptr;
00411 }
00412 }
00413
00414 void G4ExcitationHandler::SetFermiModel(G4VFermiBreakUp* ptr)
00415 {
00416 if(ptr && ptr != theFermiModel) {
00417 delete theFermiModel;
00418 theFermiModel = ptr;
00419 }
00420 }
00421
00422 void
00423 G4ExcitationHandler::SetPhotonEvaporation(G4VEvaporationChannel* ptr)
00424 {
00425 if(ptr && ptr != thePhotonEvaporation) {
00426 thePhotonEvaporation = ptr;
00427 theEvaporation->SetPhotonEvaporation(ptr);
00428 }
00429 }
00430
00431 void G4ExcitationHandler::SetMaxZForFermiBreakUp(G4int aZ)
00432 {
00433 maxZForFermiBreakUp = aZ;
00434 }
00435
00436 void G4ExcitationHandler::SetMaxAForFermiBreakUp(G4int anA)
00437 {
00438 maxAForFermiBreakUp = std::min(5,anA);
00439 }
00440
00441 void G4ExcitationHandler::SetMaxAandZForFermiBreakUp(G4int anA, G4int aZ)
00442 {
00443 SetMaxAForFermiBreakUp(anA);
00444 SetMaxZForFermiBreakUp(aZ);
00445 }
00446
00447 void G4ExcitationHandler::SetMinEForMultiFrag(G4double anE)
00448 {
00449 minEForMultiFrag = anE;
00450 }