Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/hadronic/models/de_excitation/handler
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/handler/include/G4ExcitationHandler.hh

    r819 r962  
    2525//
    2626//
    27 // $Id: G4ExcitationHandler.hh,v 1.7 2006/08/19 19:55:59 dennis Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4ExcitationHandler.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (May 1998)
     32//
     33// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     34// cross section option
     35//
    3236// Modif (30 June 1998) by V. Lara:
    3337//      -Using G4ParticleTable and therefore G4IonTable
    … …  
    3842//              MultiFragmentation: G4DummyMF (a dummy one)
    3943//              Fermi Breakup model: G4StatFermiBreakUp
    40 
     44//
     45// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     46// cross section option
     47// JMQ (06 September 2008) Also external choices have been added for
     48// superimposed Coulomb barrier (if useSICBis set true, by default is false)
    4149
    4250#ifndef G4ExcitationHandler_h
    … …  
    9199  void SetMinEForMultiFrag(G4double anE);
    92100
     101// for inverse cross section choice
     102  inline void SetOPTxs(G4int opt) { OPTxs = opt;}
     103// for superimposed Coulomb Barrir for inverse cross sections
     104  inline void UseSICB(){useSICB=true;}
     105
    93106private:
    94107
    … …  
    131144  G4bool MyOwnFermiBreakUpClass;
    132145  G4bool MyOwnPhotonEvaporationClass;
     146
     147  G4int OPTxs;
     148  G4bool useSICB;
    133149
    134150  struct DeleteFragment

  • trunk/source/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc

    r819 r962  
    2727// Hadronic Process: Nuclear De-excitations
    2828// by V. Lara (May 1998)
     29//
     30// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     31// cross section option
     32//
     33// Modif (24 Jul 2008) by M. A. Cortes Giraldo:
     34//      -Max Z,A for Fermi Break-Up turns to 9,17 by default
     35//      -BreakItUp() reorganised and bug in Evaporation loop fixed
     36//      -Transform() optimised
    2937// Modif (30 June 1998) by V. Lara:
    3038//      -Modified the Transform method for use G4ParticleTable and
    … …  
    3341//       G4DynamicParticle
    3442//      -It uses default algorithms for:
    35 //              Evaporation: G4StatEvaporation
    36 //              MultiFragmentation: G4DummyMF (a dummy one)
    37 //              Fermi Breakup model: G4StatFermiBreakUp
    38 
     43//              Evaporation: G4Evaporation
     44//              MultiFragmentation: G4StatMF
     45//              Fermi Breakup model: G4FermiBreakUp
     46//
     47// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
     48// cross section option (default OPTxs=3)
     49// JMQ (06 September 2008) Also external choices have been added for
     50// superimposed Coulomb barrier (if useSICBis set true, by default is false)
    3951
    4052#include "G4ExcitationHandler.hh"
    … …  
    4557
    4658G4ExcitationHandler::G4ExcitationHandler():
     59  // Fermi BreakUp is on and MultiFrag is off by default
     60  //  maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4.0*GeV),
    4761  maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV),
    4862  MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true),
    49   MyOwnPhotonEvaporationClass(true)
     63  MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false)
    5064{                                                                         
    5165  theTableOfParticles = G4ParticleTable::GetParticleTable();
    … …  
    92106}
    93107
    94 
    95 G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment &theInitialState) const
    96 {
    97 
    98   G4FragmentVector * theResult = 0;
     108////////////////////////////////////////////////////////////////////////////////////////////////
     109/// 25/07/08 16:45  Proposed by MAC ////////////////////////////////////////////////////////////
     110////////////////////////////////////////////////////////////////////////////////////////////////
     111
     112G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
     113{       
     114  //for inverse cross section choice
     115  theEvaporation->SetOPTxs(OPTxs);
     116  //for the choice of superimposed Coulomb Barrier for inverse cross sections
     117  theEvaporation->UseSICB(useSICB);
     118
     119  // Pointer which will be used to return the final production vector
     120  G4FragmentVector * theResult = 0;
     121
     122  // Variables existing until end of method
     123  G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
     124  G4Fragment theExcitedNucleus;              // object to be passed in BreakItUp methods
     125  G4FragmentVector * theTempResult = 0;      // pointer which receives temporal results
     126  std::list<G4Fragment*> theEvapList;        // list to apply Evaporation, SMF or Fermi Break-Up
     127  std::list<G4Fragment*> theEvapStableList;  // list to apply PhotonEvaporation
     128  std::list<G4Fragment*> theFinalStableList; // list to store final result
     129  std::list<G4Fragment*>::iterator iList;
     130
     131  // Variables to describe the excited configuration
    99132  G4double exEnergy = theInitialState.GetExcitationEnergy();
    100   //  G4cout << " first exEnergy in MeV: " << exEnergy/MeV << G4endl;
    101   G4double A = theInitialState.GetA();
    102   G4int Z = static_cast<G4int>(theInitialState.GetZ());
    103   G4FragmentVector* theTempResult = 0;
    104   G4Fragment theExcitedNucleus;
    105 
    106   // Test applicability
    107   if (A > 4)
     133  G4int A = static_cast<G4int>( theInitialState.GetA() +0.5 );
     134  G4int Z = static_cast<G4int>( theInitialState.GetZ() +0.5 );
     135
     136  // In case A <= 4 the fragment will not perform any nucleon emission
     137  if (A <= 4)
    108138    {
    109       // Initial State De-Excitation
    110       if(A<GetMaxA()&&Z<GetMaxZ())
    111         //     && exEnergy>G4NucleiPropertiesTable::GetBindingEnergy(Z,A)) {
    112         {
    113           theResult = theFermiModel->BreakItUp(theInitialState);
    114         }
    115       else   if (exEnergy>GetMinE()*A)
    116         {
    117           theResult = theMultiFragmentation->BreakItUp(theInitialState);
    118         }
    119       else
    120         {
    121           theResult = theEvaporation->BreakItUp(theInitialState);
    122         }
    123    
    124 
    125 
    126 
    127       // De-Excitation loop
    128       // ------------------
    129       // Check if there are excited fragments
    130       std::list<G4Fragment*> theResultList;
    131       G4FragmentVector::iterator j;
    132       std::list<G4Fragment*>::iterator i;
    133       for (j = theResult->begin(); j != theResult->end();j++)
    134         {
    135           theResultList.push_back(*j);
    136         }
    137       theResult->clear();
    138       for (i = theResultList.begin(); i != theResultList.end(); i++)
    139         {
    140           exEnergy = (*i)->GetExcitationEnergy();
    141           //      G4cout << " exEnergy in MeV: " << exEnergy/MeV << G4endl;
    142           if (exEnergy > 0.0)
    143             {
    144               A = (*i)->GetA();
    145               Z = static_cast<G4int>((*i)->GetZ());
    146               theExcitedNucleus = *(*i);
    147               // try to de-excite this fragment
    148               if( A < GetMaxA() && Z < GetMaxZ() )
    149                 // && exEnergy>G4NucleiPropertiesTable::GetBindingEnergy(Z,A))
    150                 {
    151                   // Fermi Breakup not now called for for exotic fragments for good reasons...
    152                   // theTempResult = theFermiModel->BreakItUp(theExcitedNucleus);
    153                   //if (theTempResult->size() == 1)
    154                   //  {
    155                   //    std::for_each(theTempResult->begin(),theTempResult->end(), G4Delete());
    156                   //    delete theTempResult;
    157                   //  }
    158                   theTempResult = theEvaporation->BreakItUp(theExcitedNucleus);
    159                 }
    160               else
    161                 {
    162                   // Evaporation
    163                   theTempResult = theEvaporation->BreakItUp(theExcitedNucleus);
    164                 }
    165               // The Nucleus has been fragmented?
    166               if (theTempResult->size() > 1)
    167                 // If so :
    168                 {
    169                   // Remove excited fragment from the result
    170                   //    delete theResult->removeAt(i--);
    171                   delete (*i);
    172                   i = theResultList.erase(i);
    173                   // and add theTempResult elements to theResult
    174                   for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin();
    175                        ri != theTempResult->rend(); ++ri)
    176                     {
    177                       theResultList.push_back(*ri);
    178                     }
    179                   delete theTempResult;
    180                 }
    181               else
    182                 // If not :
    183                 {
    184                   // it doesn't matter, we Follow with the next fragment but
    185                   // I have to clean up
    186                   std::for_each(theTempResult->begin(),theTempResult->end(), G4Delete());
    187                   delete theTempResult;
    188                 }
    189             }
    190         }
    191       for (i = theResultList.begin(); i != theResultList.end(); i++)
    192         {
    193           theResult->push_back(*i);
    194         }
    195       theResultList.clear();
     139      // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later
     140
     141      theEvapStableList.push_back( theInitialStatePtr );
    196142    }
    197   else   // if A > 4
     143
     144  else  // If A > 4 we try to apply theFermiModel, theMultiFragmentation or theEvaporation
    198145    {
    199       theResult = new G4FragmentVector();
    200       theResult->push_back(new G4Fragment(theInitialState));
     146
     147      // Store original state in theEvapList
     148      theEvapList.push_back( theInitialStatePtr );
     149
     150      // ------------------------------
     151      // De-excitation loop
     152      // ------------------------------
     153
     154      for (iList = theEvapList.begin(); iList != theEvapList.end(); iList++)
     155        {
     156          exEnergy = (*iList)->GetExcitationEnergy();
     157
     158          if (exEnergy > 0.0)
     159            {
     160              // Check conditions for each model
     161              theExcitedNucleus = *(*iList);
     162              A = static_cast<G4int>((*iList)->GetA()+0.5);  // +0.5 to avoid bad truncation
     163              Z = static_cast<G4int>((*iList)->GetZ()+0.5);
     164
     165              if ( A < GetMaxA() && Z < GetMaxZ() ) // if satisfied apply Fermi Break-Up
     166                {
     167                  theTempResult = theFermiModel->BreakItUp(theExcitedNucleus);
     168                }
     169              else if (exEnergy > GetMinE()*A)   // if satisfied apply SMF
     170                {
     171                  theTempResult = theMultiFragmentation->BreakItUp(theExcitedNucleus);
     172                }
     173              else // apply Evaporation in another case
     174                {
     175                  theTempResult = theEvaporation->BreakItUp(theExcitedNucleus);
     176                }
     177
     178              // New configuration is stored in theTempResult, so we can free
     179              // the memory where the previous configuration is
     180
     181              delete (*iList);
     182
     183              // And now the theTempResult->size() tells us if the configuration has changed
     184
     185              if ( theTempResult->size() > 1 )
     186                {
     187                  // push_back the result to the end of theEvapList (this same list)
     188
     189                  for (G4FragmentVector::iterator j = theTempResult->begin();
     190                       j != theTempResult->end(); j++)
     191                    {
     192                      theEvapList.push_back(*j);
     193                    }
     194                }
     195              else
     196                {
     197                  // push_back the result to theEvapStableList, because
     198                  // is still excited, but cannot emmit more nucleons
     199
     200                  for (G4FragmentVector::iterator j = theTempResult->begin();
     201                       j != theTempResult->end(); j++)
     202                    {
     203                      theEvapStableList.push_back(*j);
     204                    }
     205                }
     206
     207              // after working with theTempResult, clear and delete it
     208              theTempResult->clear();
     209              delete theTempResult;
     210
     211            }
     212          else // exEnergy = 0.0
     213            {
     214              // if this fragment is at ground state,
     215              // store it in theFinalStableList
     216
     217              theFinalStableList.push_back(*iList);
     218
     219            } // endif (exEnergy > 0.0)
     220
     221        } // end of the loop over theEvapList
     222
     223      theEvapList.clear();   // clear all the list and do not free memory pointed by
     224                             // each element because this have been done before!
    201225    }
    202226
    203227  // Now we try to deexcite by means of PhotonEvaporation those fragments
    204   // which are excited.
    205  
    206   theTempResult = 0;
    207   std::list<G4Fragment*> theFinalResultList;
    208   //AHtest  std::list<G4Fragment*> theFinalPhotonResultList;
    209   std::list<G4Fragment*> theResultList;
    210   std::list<G4Fragment*>::iterator j;
    211   G4FragmentVector::iterator i;
    212   for (i = theResult->begin(); i != theResult->end();i++)
     228  // which are still excited.
     229
     230
     231  // -----------------------
     232  // Photon-Evaporation loop
     233  // -----------------------
     234
     235  for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); iList++)
    213236    {
    214       theResultList.push_back(*i);
    215       //      G4cout << " Before loop list energy in MeV: " << ((*i)->GetExcitationEnergy())/MeV << G4endl;
    216     }
    217   theResult->clear();
    218 
    219   for (j = theResultList.begin(); j != theResultList.end(); j++) {
    220     //    G4cout << " Test loop list: " << (*j)->GetExcitationEnergy() << " size: " << theResultList.size() << G4endl;
    221   }
    222 
    223   //  for (j = theResultList.begin(); j != theResultList.end(); j++)
    224   j = theResultList.begin();  //AH
    225   while (j != theResultList.end()) //AH
    226     {
    227       if ((*j)->GetA() > 1 && (*j)->GetExcitationEnergy() > 0.1*eV)
    228         {
    229           theExcitedNucleus = *(*j);
    230           theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus);
    231           // If Gamma Evaporation has succeed then
    232           if (theTempResult->size() > 1)
    233             {
    234               // Remove excited fragment from the result
    235               //              delete (*j);
    236               //              theResultList.erase(j--);
    237               //              theResultList.erase(j); don't delete as there's no push back...
    238               // and add theTempResult elements to theResult
    239               for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin();
     237      A = static_cast<G4int>((*iList)->GetA()+0.5);
     238      exEnergy = (*iList)->GetExcitationEnergy();
     239
     240      if ( A > 1 && exEnergy > 0.1*eV ) // if so, photon-evaporation is applied
     241        {
     242          theExcitedNucleus = *(*iList);
     243          theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus);
     244
     245          // if there is a gamma emission then
     246          if (theTempResult->size() > 1)
     247            {
     248              // first free the memory occupied by the previous state
     249              delete (*iList);
     250
     251              // and now add the final state from gamma emission to the end of
     252              // theEvapStableList
     253              for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin();
    240254                   ri != theTempResult->rend(); ++ri)
    241                 {
     255                // reversed is applied in order to have residual nucleus in first position
     256                {
     257
    242258#ifdef PRECOMPOUND_TEST
    243                   if ((*ri)->GetA() == 0)
    244                     (*ri)->SetCreatorModel(G4String("G4PhotonEvaporation"));
    245                   else
    246                     (*ri)->SetCreatorModel(G4String("ResidualNucleus"));
    247 #endif
    248                   theResultList.push_back(*ri);
    249                   //AHtest                theFinalPhotonResultList.push_back(*ri);
    250                   //              theFinalResultList.push_back(*ri); don't add to final result as they'll go through the loop
    251                 }
    252               delete theTempResult;
    253             }
    254           // In other case, just clean theTempResult and continue
    255           else
    256             {
    257               std::for_each(theTempResult->begin(), theTempResult->end(), DeleteFragment());
    258               delete theTempResult;
     259                  if ((*ri)->GetA() == 0)
     260                    (*ri)->SetCreatorModel(G4String("G4PhotonEvaporation"));
     261                  else
     262                    (*ri)->SetCreatorModel(G4String("ResidualNucleus"));
     263#endif
     264
     265                  theEvapStableList.push_back(*ri);
     266                }
     267             
     268              // now we clean and remove the temporal vector
     269              theTempResult->clear();
     270              delete theTempResult;
     271
     272            }
     273          else  // if theTempResult->size() = 1
     274            {
     275              // if there is not any gamma emission from this excited fragment
     276              // we have to emmit a gamma which forces the deexcitation
     277
     278              // First I clean completely theTempResult
     279
     280              for (G4FragmentVector::iterator j = theTempResult->begin();
     281                   j != theTempResult->end(); j++)
     282                {
     283                  delete (*j);
     284                }
     285
     286              theTempResult->clear();
     287              delete theTempResult;
     288
    259289#ifdef debugphoton
    260290              G4cout << "G4ExcitationHandler: Gamma Evaporation could not deexcite the nucleus: \n"
    … …  
    263293                     << "-----------------------------------------------------------------------\n";
    264294#endif
    265               G4double GammaEnergy = (*j)->GetExcitationEnergy();
    266               G4double cosTheta = 1. - 2. * G4UniformRand();
    267               G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
     295
     296              // Let's create a G4Fragment pointer representing the gamma emmited
     297              G4double GammaEnergy = (*iList)->GetExcitationEnergy();
     298              G4double cosTheta = 1. - 2. * G4UniformRand();
     299              G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
    268300              G4double phi = twopi * G4UniformRand();
    269301              G4ThreeVector GammaP(GammaEnergy * sinTheta * std::cos(phi),
    … …  
    272304              G4LorentzVector Gamma4P(GammaP,GammaEnergy);
    273305              G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition());
    274              
    275              
    276              
    277               G4double Mass = (*j)->GetGroundStateMass();
    278               G4ThreeVector ResidualP((*j)->GetMomentum().vect() - GammaP);
     306
     307              // And now we update momentum and energy for the nucleus
     308              G4double Mass = (*iList)->GetGroundStateMass();
     309              G4ThreeVector ResidualP((*iList)->GetMomentum().vect() - GammaP);
    279310              G4double ResidualE = std::sqrt(ResidualP*ResidualP + Mass*Mass);
    280311              G4LorentzVector Residual4P(ResidualP,ResidualE);
    281               (*j)->SetMomentum(Residual4P);
    282              
    283              
    284        
     312              (*iList)->SetMomentum(Residual4P); // Now this fragment has been deexcited!
     313
     314              // we store the deexcited fragment in theFinalStableList
     315              theFinalStableList.push_back(*iList);
     316
    285317#ifdef PRECOMPOUND_TEST
    286               theHandlerPhoton->SetCreatorModel("G4ExcitationHandler");
    287 #endif
    288 //            theFinalPhotonResultList.push_back( theHandlerPhoton );
    289 //            G4cout << " adding photon fragment " << G4endl;
    290               theResultList.push_back( theHandlerPhoton );
    291               //              theFinalResultList.push_back( theHandlerPhoton );
    292               theFinalResultList.push_back(*j);
     318              theHandlerPhoton->SetCreatorModel("G4ExcitationHandler");
     319#endif
     320
     321              // Finally, we add theHandlerPhoton to theFinalStableList
     322              theFinalStableList.push_back(theHandlerPhoton);
     323
    293324#ifdef debugphoton
    294325              G4cout << "Emmited photon:\n"
    295                      << theResultList.back() << '\n'
     326                     << theFinalStableList.back() << '\n'
    296327                     << "Residual nucleus after photon emission:\n"
    297                      << *(*j) << '\n'
     328                     << *(*iList) << '\n'
    298329                     << "-----------------------------------------------------------------------\n";
    299330#endif
    300               //test          j++; // AH only increment if not erased:
    301             }   
    302         } else {
    303           //test          j++; // AH increment iterator if a proton or excitation energy small
    304           theFinalResultList.push_back(*j);
     331
     332            }
    305333        }
    306 //      G4cout << " Inside loop list: " << (*j)->GetExcitationEnergy() << " size: " << theFinalResultList.size() << G4endl;
    307       j++;
     334      else // case of a nucleon, gamma or very small excitation energy
     335        {
     336          // we don't have to do anything, just store the fragment in theFinalStableList
     337
     338          theFinalStableList.push_back(*iList);
     339       
     340        } // A > 1 && exEnergy > 0.1*eV
     341
     342    } // end of photon-evaporation loop
     343
     344
     345  // The deexcitation from fragments inside theEvapStableList has been finished, so...
     346  theEvapStableList.clear();
     347
     348  // Now the final state is in theFinalStableList, and we have to send it to theResult vector
     349
     350  theResult = new G4FragmentVector;
     351  theResult->reserve( theFinalStableList.size() * sizeof(G4Fragment*) );
     352  // We reserve enough memory to optimise the storing speed
     353
     354  for (iList = theFinalStableList.begin(); iList != theFinalStableList.end(); iList++)
     355    {
     356      theResult->push_back(*iList);
    308357    }
    309   //  for (j = theResultList.begin(); j != theResultList.end(); j++)
    310   for (j = theFinalResultList.begin(); j != theFinalResultList.end(); j++)
    311     {
    312       theResult->push_back(*j);
    313     }
    314 
    315 //AHtest   for (j = theFinalPhotonResultList.begin(); j != theFinalPhotonResultList.end(); j++)
    316 //AHtest     {
    317 //AHtest       theResult->push_back(*j);
    318 //AHtest       number_results++;
    319 //AHtest     }
    320 
    321 
    322   theResultList.clear();
    323   theFinalResultList.clear();
    324   //AHtest  theFinalPhotonResultList.clear();
    325  
    326  
     358
     359  // After storing the final state , we can clear theFinalStableList
     360  theFinalStableList.clear();
     361
    327362#ifdef debug
    328363  CheckConservation(theInitialState,theResult);
    329364#endif
    330   // Change G4FragmentVector by G4DynamicParticle
     365
     366  // Change G4FragmentVector* to G4ReactionProductVector*
    331367  return Transform(theResult);
    332368}
     369
     370
    333371
    334372G4ReactionProductVector *
    … …  
    348386  theNeutron->SetVerboseLevel(2);
    349387  G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
     388
     389  // MAC (24/07/08)
     390  // To optimise the storing speed, we reserve space in memory for the vector
     391  theReactionProductVector->reserve( theFragmentVector->size() * sizeof(G4ReactionProduct*) );
     392
    350393  G4int theFragmentA, theFragmentZ;
    351394  G4LorentzVector theFragmentMomentum;
Note: See TracChangeset for help on using the changeset viewer.