G4Evaporation.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 //
00027 // $Id$
00028 //
00029 // Hadronic Process: Nuclear De-excitations
00030 // by V. Lara (Oct 1998)
00031 //
00032 // Alex Howard - added protection for negative probabilities in the sum, 14/2/07
00033 //
00034 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
00035 // cross section option
00036 // JMQ (06 September 2008) Also external choices have been added for 
00037 // superimposed Coulomb barrier (if useSICBis set true, by default is false) 
00038 //
00039 // V.Ivanchenko (27 July 2009)  added G4EvaporationDefaultGEMFactory option
00040 // V.Ivanchenko (10 May  2010)  rewrited BreakItUp method: do not make new/delete
00041 //                              photon channel first, fission second,
00042 //                              added G4UnstableFragmentBreakUp to decay 
00043 //                              unphysical fragments (like 2n or 2p),
00044 //                              use Z and A integer
00045 // V.Ivanchenko (22 April 2011) added check if a fragment can be deexcited 
00046 //                              by the FermiBreakUp model
00047 // V.Ivanchenko (23 January 2012) added pointer of G4VPhotonEvaporation 
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 G4Evaporation::G4Evaporation(std::vector<G4VEvaporationChannel*>* channels) 
00081   : theChannels(channels),theChannelFactory(0),nChannels(0)
00082 {
00083   // are input relible?
00084   G4bool accepted = true;
00085   if(!theChannels) { accepted = false; }
00086   else if(0 == theChannels->size()) { accepted = false; }
00087 
00088   if(accepted) {
00089     SetPhotonEvaporation((*channels)[0]); 
00090   } else {
00091     SetPhotonEvaporation(new G4PhotonEvaporation());
00092     theChannelFactory = new G4EvaporationDefaultGEMFactory(thePhotonEvaporation); 
00093   }
00094   SetParameters();
00095   InitialiseEvaporation();
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   // The residual nucleus (after evaporation of each fragment)
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   // Starts loop over evaporated particles, loop is limited by number
00180   // of nucleons
00181   for(G4int ia=0; ia<Amax; ++ia) {
00182  
00183     // g,n,p and light fragments - evaporation is finished
00184     G4int Z = theResidualNucleus->GetZ_asInt();
00185     G4int A = theResidualNucleus->GetA_asInt();
00186 
00187     // stop deecitation loop if can be deexcited by FBU
00188     if(maxZforFBU > Z && maxAforFBU >= A) {
00189       theResult->push_back(theResidualNucleus);
00190       return theResult;
00191     }
00192 
00193     // check if it is stable, then finish evaporation
00194     G4double abun = nist->GetIsotopeAbundance(Z, A); 
00195     /*
00196     G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z
00197            << " A= " << A << " Eex(MeV)= " 
00198            << theResidualNucleus->GetExcitationEnergy()
00199            << " aban= " << abun << G4endl;
00200     */
00201     // stop deecitation loop in the case of the cold stable fragment 
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     G4cout << "### Evaporation loop #" << ia 
00212            << "  Fragment: " << theResidualNucleus << G4endl;
00213     */
00214     // loop over evaporation channels
00215     for(i=0; i<nChannels; ++i) {
00216       prob = (*theChannels)[i]->GetEmissionProbability(theResidualNucleus);
00217       //G4cout << "  Channel# " << i << "  prob= " << prob << G4endl; 
00218 
00219       totprob += prob;
00220       probabilities[i] = totprob;
00221       // if two recent probabilities are near zero stop computations
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       // protection for very excited fragment - avoid GEM
00230       if(7 == i && Eex > Elimit*A) {
00231         maxchannel = 8;
00232         break;
00233       }
00234     }
00235 
00236     // photon evaporation in the case of no other channels available
00237     // do evaporation chain and reset total probability
00238     if(0.0 < totprob && probabilities[0] == totprob) {
00239       //G4cout << "Start gamma evaporation" << G4endl;
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     // stable fragnent - evaporation is finished
00252     if(0.0 == totprob) {
00253 
00254       // if fragment is exotic, then try to decay it
00255       if(0.0 == abun && Z < 20) {
00256         //G4cout << "$$$ Decay exotic fragment" << G4endl;
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       // save residual fragment
00268       theResult->push_back(theResidualNucleus);
00269       return theResult;
00270     }
00271 
00272 
00273     // select channel
00274     totprob *= G4UniformRand();
00275     // loop over evaporation channels
00276     for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
00277 
00278     // this should not happen
00279     if(i >= nChannels) { i = nChannels - 1; }
00280 
00281 
00282     // single photon evaporation, primary pointer is kept
00283     if(0 == i) {
00284       //G4cout << "Single gamma" << G4endl;
00285       G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus);
00286       if(gamma) { theResult->push_back(gamma); }
00287 
00288       // fission, return results to the main loop if fission is succesful
00289     } else if(1 == i) {
00290       //G4cout << "Fission" << G4endl;
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       // other channels
00305     } else {
00306       //G4cout << "Channel # " << i << G4endl;
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           // if the residual change its pointer 
00316           // then delete previous residual fragment and update to the new
00317           if(theResidualNucleus != (*theTempResult)[nsec] ) { 
00318             delete theResidualNucleus; 
00319             theResidualNucleus = (*theTempResult)[nsec];
00320           }
00321         }
00322         delete theTempResult;
00323       }
00324     }
00325   }
00326 
00327   // loop is stopped, save residual
00328   theResult->push_back(theResidualNucleus);
00329   
00330   return theResult;
00331 }
00332 

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