G4ExcitationHandler.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4ExcitationHandler.cc 66934 2013-01-21 13:18:35Z vnivanch $
00027 //
00028 // Hadronic Process: Nuclear De-excitations
00029 // by V. Lara (May 1998)
00030 //
00031 //
00032 // Modified:
00033 // 30 June 1998 by V. Lara:
00034 //      -Modified the Transform method for use G4ParticleTable and 
00035 //       therefore G4IonTable. It makes possible to convert all kind 
00036 //       of fragments (G4Fragment) produced in deexcitation to 
00037 //       G4DynamicParticle
00038 //      -It uses default algorithms for:
00039 //              Evaporation: G4Evaporation
00040 //              MultiFragmentation: G4StatMF 
00041 //              Fermi Breakup model: G4FermiBreakUp
00042 // 24 Jul 2008 by M. A. Cortes Giraldo:
00043 //      -Max Z,A for Fermi Break-Up turns to 9,17 by default
00044 //      -BreakItUp() reorganised and bug in Evaporation loop fixed
00045 //      -Transform() optimised
00046 // (September 2008) by J. M. Quesada. External choices have been added for :
00047 //      -inverse cross section option (default OPTxs=3)
00048 //      -superimposed Coulomb barrier (if useSICB is set true, by default it is false) 
00049 // September 2009 by J. M. Quesada: 
00050 //      -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
00051 // 27 Nov 2009 by V.Ivanchenko: 
00052 //      -cleanup the logic, reduce number internal vectors, fixed memory leak.
00053 // 11 May 2010 by V.Ivanchenko: 
00054 //      -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for 
00055 //       final photon deexcitation; used check on adundance of a fragment, decay 
00056 //       unstable fragments with A <5
00057 // 22 March 2011 by V.Ivanchenko: general cleanup and addition of a condition: 
00058 //       products of Fermi Break Up cannot be further deexcited by this model 
00059 // 30 March 2011 by V.Ivanchenko removed private inline methods, moved Set methods 
00060 //       to the source
00061 // 23 January 2012 by V.Ivanchenko general cleanup including destruction of 
00062 //    objects, propagate G4PhotonEvaporation pointer to G4Evaporation class and 
00063 //    not delete it here 
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   //for inverse cross section choice
00111   theEvaporation->SetOPTxs(OPTxs);
00112   //for the choice of superimposed Coulomb Barrier for inverse cross sections
00113   theEvaporation->UseSICB(useSICB);
00114   theEvaporation->Initialise();
00115 }
00116 
00117 G4ReactionProductVector * 
00118 G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
00119 {       
00120   //G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@" << G4endl;
00121   
00122   // Variables existing until end of method
00123   G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
00124 
00125   G4FragmentVector * theTempResult = 0;      // pointer which receives temporal results
00126   std::list<G4Fragment*> theEvapList;        // list to apply Evaporation or Fermi Break-Up
00127   std::list<G4Fragment*> thePhotoEvapList;   // list to apply PhotonEvaporation
00128   std::list<G4Fragment*> theResults;         // list to store final result
00129   //
00130   //G4cout << theInitialState << G4endl;  
00131   
00132   // Variables to describe the excited configuration
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   // In case A <= 1 the fragment will not perform any nucleon emission
00140   if (A <= 1)
00141     {
00142       theResults.push_back( theInitialStatePtr );
00143     }
00144   // check if a fragment is stable
00145   else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0)
00146     {
00147       theResults.push_back( theInitialStatePtr );
00148     }  
00149   else  
00150     {      
00151       // JMQ 150909: first step in de-excitation is treated separately 
00152       // Fragments after the first step are stored in theEvapList 
00153       // Statistical Multifragmentation will take place only once
00154       //
00155       // move to evaporation loop
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               // secondary are produced
00170               // Sort out secondary fragments
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                 // gamma, p, n
00178                 if(A <= 1) { theResults.push_back(*j); }
00179 
00180                 // Analyse fragment A > 1
00181                 else {
00182                   G4double exEnergy1 = (*j)->GetExcitationEnergy();
00183 
00184                   // cold fragments
00185                   if(exEnergy1 < minExcitation) {
00186                     Z = (*j)->GetZ_asInt(); 
00187                     if(nist->GetIsotopeAbundance(Z, A) > 0.0) { 
00188                       theResults.push_back(*j); // stable fragment 
00189 
00190                     } else {
00191 
00192                       // check if the cold fragment is from FBU pool
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                         // cold fragment may be unstable
00199                       } else {
00200                         theEvapList.push_back(*j); 
00201                       }
00202                     }
00203 
00204                     // hot fragments are unstable
00205                   } else { theEvapList.push_back(*j); } 
00206                 }
00207               }
00208               if( deletePrimary ) { delete theInitialStatePtr; }
00209             }
00210             delete theTempResult;
00211           }
00212         }
00213     }
00214   /*
00215   G4cout << "## After first step " << theEvapList.size() << " for evap;  "
00216    << thePhotoEvapList.size() << " for photo-evap; " 
00217    << theResults.size() << " results. " << G4endl; 
00218   */
00219   // -----------------------------------
00220   // FermiBreakUp and De-excitation loop
00221   // -----------------------------------
00222       
00223   std::list<G4Fragment*>::iterator iList;
00224   for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList)
00225     {
00226       //G4cout << "Next evaporate: " << G4endl;  
00227       //G4cout << *iList << G4endl;  
00228       A = (*iList)->GetA_asInt(); 
00229       Z = (*iList)->GetZ_asInt();
00230           
00231       // Fermi Break-Up 
00232       G4bool wasFBU = false;
00233       if (A < maxAForFermiBreakUp && Z < maxZForFermiBreakUp) 
00234         {
00235           theTempResult = theFermiModel->BreakItUp(*(*iList));
00236           wasFBU = true; 
00237           // if initial fragment returned unchanged try to evaporate it
00238           if(1 == theTempResult->size()) {
00239             delete theTempResult;
00240             theTempResult = theEvaporation->BreakItUp(*(*iList)); 
00241           }
00242         }
00243       else // apply Evaporation in another case
00244         {
00245           theTempResult = theEvaporation->BreakItUp(*(*iList));
00246         }
00247       
00248       G4bool deletePrimary = true;
00249       size_t nsec = theTempResult->size();
00250       //G4cout << "Nproducts= " << nsec << G4endl;  
00251                   
00252       // Sort out secondary fragments
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           //G4cout << *j << G4endl;  
00259           A = (*j)->GetA_asInt();
00260           exEnergy = (*j)->GetExcitationEnergy();
00261 
00262           if(A <= 1) { theResults.push_back(*j); }    // gamma, p, n
00263 
00264           // evaporation is not possible
00265           else if(1 == nsec) { 
00266             if(exEnergy < minExcitation) { theResults.push_back(*j); }
00267             else                         { thePhotoEvapList.push_back(*j); }
00268 
00269           } else { // Analyse fragment
00270 
00271             // cold fragment
00272             if(exEnergy < minExcitation) {
00273               Z = (*j)->GetZ_asInt();
00274 
00275               // natural isotope
00276               if(nist->GetIsotopeAbundance(Z, A) > 0.0) { 
00277                 theResults.push_back(*j); // stable fragment 
00278 
00279               } else {
00280                 const G4VFermiFragment* ffrag = thePool->GetFragment(Z, A);
00281 
00282                 // isotope from FBU pool
00283                 if(ffrag) {
00284                   if(ffrag->IsStable()) { theResults.push_back(*j); }
00285                   else                  { theEvapList.push_back(*j); }
00286 
00287                   // isotope may be unstable
00288                 } else {
00289                   theEvapList.push_back(*j);
00290                 }   
00291               }
00292 
00293               // hot fragment
00294             } else if (wasFBU) { 
00295               thePhotoEvapList.push_back(*j); // FBU applied only once 
00296             } else {  
00297               theEvapList.push_back(*j);        
00298             }
00299           }
00300         }
00301       }
00302       if( deletePrimary ) { delete (*iList); }
00303       delete theTempResult;
00304     } // end of the loop over theEvapList
00305 
00306   //G4cout << "## After 2nd step " << theEvapList.size() << " was evap;  "
00307   // << thePhotoEvapList.size() << " for photo-evap; " 
00308   // << theResults.size() << " results. " << G4endl; 
00309       
00310   // -----------------------
00311   // Photon-Evaporation loop
00312   // -----------------------
00313   
00314   // at this point only photon evaporation is possible
00315   for(iList = thePhotoEvapList.begin(); iList != thePhotoEvapList.end(); ++iList)
00316     {
00317       //G4cout << "Next photon evaporate: " << thePhotonEvaporation << G4endl;  
00318       //G4cout << *iList << G4endl;
00319       exEnergy = (*iList)->GetExcitationEnergy();
00320 
00321       // only hot fragments
00322       if(exEnergy >= minExcitation) {  
00323         theTempResult = thePhotonEvaporation->BreakUpFragment(*iList);    
00324         size_t nsec = theTempResult->size();
00325         //G4cout << "Nproducts= " << nsec << G4endl;  
00326           
00327         // if there is a gamma emission then
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       // priamry fragment is kept
00340       theResults.push_back(*iList); 
00341 
00342     } // end of photon-evaporation loop
00343 
00344   //G4cout << "## After 3d step " << theEvapList.size() << " was evap;  "
00345   //     << thePhotoEvapList.size() << " was photo-evap; " 
00346   //     << theResults.size() << " results. " << G4endl; 
00347     
00348   G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
00349 
00350   // MAC (24/07/08)
00351   // To optimise the storing speed, we reserve space in memory for the vector
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) {       // photon or e-
00363         theKindOfFragment = (*i)->GetParticleDefinition();   
00364       } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
00365         theKindOfFragment = G4Neutron::NeutronDefinition();
00366       } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
00367         theKindOfFragment = G4Proton::ProtonDefinition();
00368       } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
00369         theKindOfFragment = G4Deuteron::DeuteronDefinition();
00370       } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
00371         theKindOfFragment = G4Triton::TritonDefinition();
00372       } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
00373         theKindOfFragment = G4He3::He3Definition();
00374       } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
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 }

Generated on Mon May 27 17:48:13 2013 for Geant4 by  doxygen 1.4.7