Changeset 1228 for trunk/source/processes/hadronic/models/de_excitation/handler

Timestamp:

Jan 8, 2010, 11:56:51 AM (14 years ago)

Author:

garnier

Message:

update geant4.9.3 tag

Location:

trunk/source/processes/hadronic/models/de_excitation/handler

Files:

• include/G4ExcitationHandler.hh (modified) (1 diff)
• src/G4ExcitationHandler.cc (modified) (4 diffs)

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

@@ -25 +25 @@
 //
 //
-// $Id: G4ExcitationHandler.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $
-// GEANT4 tag $Name: geant4-09-03-cand-01 $
+// $Id: G4ExcitationHandler.hh,v 1.10 2009/11/24 11:18:48 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // Hadronic Process: Nuclear De-excitations

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

@@ -52 +52 @@
 
 #include "G4ExcitationHandler.hh"
+#include "globals.hh"
+#include "G4LorentzVector.hh"
 #include <list>
 
@@ -121 +123 @@
   
   // Pointer which will be used to return the final production vector
-  G4FragmentVector * theResult = 0;
+  //G4FragmentVector * theResult = new G4FragmentVector;
   
   // Variables existing until end of method
-  G4Fragment * theInitialStatePtr = const_cast<G4Fragment*>(&theInitialState);
-  //  G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
-  G4Fragment theExcitedNucleus;              // object to be passed in BreakItUp methods
+  //G4Fragment * theInitialStatePtr = const_cast<G4Fragment*>(&theInitialState);
+  G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
   G4FragmentVector * theTempResult = 0;      // pointer which receives temporal results
   std::list<G4Fragment*> theEvapList;        // list to apply Evaporation, SMF or Fermi Break-Up
   std::list<G4Fragment*> theEvapStableList;  // list to apply PhotonEvaporation
-  std::list<G4Fragment*> theFinalStableList; // list to store final result
+  std::list<G4Fragment*> theResults;         // list to store final result
   std::list<G4Fragment*>::iterator iList;
   //
-  
+  //G4cout << "@@@@@@@@@@ Start G4Exitation Handler @@@@@@@@@@@@@" << G4endl;  
+  //G4cout << theInitialState << G4endl;  
   
   // Variables to describe the excited configuration
@@ -145 +147 @@
   if (A <= 4)
     {
-      // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later
-      
+      // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later      
       theEvapStableList.push_back( theInitialStatePtr );
     }
@@ -172 +173 @@
         }
       
+      G4bool deletePrimary = true;
+      if(theTempResult->size() > 0) 
+        {      
+          // Store original state in theEvapList
+          G4FragmentVector::iterator j;
+          for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
+            {  
+              if((*j) == theInitialStatePtr) { deletePrimary = false; }
+              A = static_cast<G4int>((*j)->GetA()+0.5);  // +0.5 to avoid bad truncation
+
+              if(A <= 1)      { theResults.push_back(*j); }         // gamma, p, n
+              else if(A <= 4) { theEvapStableList.push_back(*j); }  // evaporation is not possible
+              else            { theEvapList.push_back(*j); }        // evaporation is possible
+            }
+        }
+      if( deletePrimary ) { delete theInitialStatePtr; }
+      delete theTempResult;
+    }
+  //
+  // JMQ 150909: Further steps in de-excitation chain follow ..
       
-      // Store original state in theEvapList
-      G4FragmentVector::iterator j;
-      for (j = theTempResult->begin(); j != theTempResult->end(); j++)
-        {  
-          theEvapList.push_back(*j);
-          //  This memory release works with a list but not with a G4FragmentVector 
-          //     delete (*j);
+  //G4cout << "## After first step " << theEvapList.size() << " for evap;  "
+  //     << theEvapStableList.size() << " for photo-evap; " 
+  //     << theResults.size() << " results. " << G4endl; 
+
+  // ------------------------------
+  // De-excitation loop
+  // ------------------------------
+      
+  for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList)
+    {
+      A = static_cast<G4int>((*iList)->GetA()+0.5);  // +0.5 to avoid bad truncation
+      Z = static_cast<G4int>((*iList)->GetZ()+0.5);
+          
+      // In case A <= 4 the fragment will not perform any nucleon emission
+      if (A <= 4)
+        {
+          // storing G4Fragment* in theEvapStableList to apply thePhotonEvaporation later    
+          theEvapStableList.push_back(*iList );
         }
-      delete theTempResult;
-      //theTempResult->clear();
-      //
-      // JMQ 150909: Further steps in de-excitation chain follow ..
+          
+      else  // If A > 4 we try to apply theFermiModel or theEvaporation
+        {   
+          // stable fragment  
+          if ((*iList)->GetExcitationEnergy() <= 0.1*eV) 
+            { 
+              theResults.push_back(*iList); 
+            }
+          else
+            {
+              if ( A < GetMaxA() && Z < GetMaxZ() ) // if satisfied apply Fermi Break-Up
+                {
+                  theTempResult = theFermiModel->BreakItUp(*(*iList));
+                }
+              else // apply Evaporation in another case
+                {
+                  theTempResult = theEvaporation->BreakItUp(*(*iList));
+                }
+                  
+              // New configuration is stored in theTempResult, so we can free
+              // the memory where the previous configuration is
+                  
+              G4bool deletePrimary = true;
+              G4int nsec = theTempResult->size();
+                  
+              // The number of secondaries tells us if the configuration has changed  
+              if ( nsec > 0 )
+                {
+                  G4FragmentVector::iterator j;
+                  for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
+                    {
+                      if((*j) == (*iList)) { deletePrimary = false; }
+                      A = static_cast<G4int>((*j)->GetA()+0.5);  // +0.5 to avoid bad truncation
+
+                      if(A <= 1)                   { theResults.push_back(*j); }        // gamma, p, n
+                      else if(A <= 4 || 1 == nsec) { theEvapStableList.push_back(*j); } // no evaporation 
+                      else                         { theEvapList.push_back(*j); }      
+                    }
+                }
+              if( deletePrimary ) { delete (*iList); }
+              delete theTempResult;
+              
+            }
+        } // endif (A <=4)
+    } // end of the loop over theEvapList
+
+  //G4cout << "## After 2nd step " << theEvapList.size() << " was evap;  "
+  //     << theEvapStableList.size() << " for photo-evap; " 
+  //     << theResults.size() << " results. " << G4endl; 
       
-      // ------------------------------
-      // De-excitation loop
-      // ------------------------------
-      
-      for (iList = theEvapList.begin(); iList != theEvapList.end(); iList++)
-        {
-          A = static_cast<G4int>((*iList)->GetA()+0.5);  // +0.5 to avoid bad truncation
-          Z = static_cast<G4int>((*iList)->GetZ()+0.5);
-          
-          // In case A <= 4 the fragment will not perform any nucleon emission
-          if (A <= 4)
-            {
-              // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later
-              
-              theEvapStableList.push_back(*iList );
-            }
-          
-          else  // If A > 4 we try to apply theFermiModel or theEvaporation
-            {
-              exEnergy = (*iList)->GetExcitationEnergy();
-              
-              if (exEnergy > 0.0)
-                {
-                  // Check conditions for each model
-                  theExcitedNucleus = *(*iList);
-                  
-                  if ( A < GetMaxA() && Z < GetMaxZ() ) // if satisfied apply Fermi Break-Up
-                    {
-                      theTempResult = theFermiModel->BreakItUp(theExcitedNucleus);
-                    }
-                  else // apply Evaporation in another case
-                    {
-                      theTempResult = theEvaporation->BreakItUp(theExcitedNucleus);
-                    }
-                  
-                  // New configuration is stored in theTempResult, so we can free
-                  // the memory where the previous configuration is
-                  
-                  delete (*iList);
-                  
-                  // And now the theTempResult->size() tells us if the configuration has changed
-                  
-                  if ( theTempResult->size() > 1 )
-                    {
-                      // push_back the result to the end of theEvapList (this same list)
-                      
-                      for (G4FragmentVector::iterator j = theTempResult->begin();
-                           j != theTempResult->end(); j++)
-                        {
-                          theEvapList.push_back(*j);
-                        }
-                    }
-                  else
-                    {
-                      // push_back the result to theEvapStableList, because
-                      // is still excited, but cannot emmit more nucleons
-                      
-                      for (G4FragmentVector::iterator j = theTempResult->begin();
-                           j != theTempResult->end(); j++)
-                        {
-                          theEvapStableList.push_back(*j);
-                        }
-                    }
-                  
-                  // after working with theTempResult, clear and delete it
-                  theTempResult->clear();
-                  delete theTempResult;
-                  
-                }
-              else // exEnergy = 0.0
-                {
-                  // if this fragment is at ground state,
-                  // store it in theFinalStableList
-                  
-                  theFinalStableList.push_back(*iList);
-                  
-                } // endif (exEnergy > 0.0)
-            } // endif (A <=4)
-        } // end of the loop over theEvapList
-      
-      theEvapList.clear();   // clear all the list and do not free memory pointed by
-                             // each element because this have been done before!
-    }
-  
-  // Now we try to deexcite by means of PhotonEvaporation those fragments
-  // which are still excited.
-  
-  
   // -----------------------
   // Photon-Evaporation loop
   // -----------------------
   
-  for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); iList++)
+  for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); ++iList)
     {
+      // take out stable particles and fragments
       A = static_cast<G4int>((*iList)->GetA()+0.5);
-      exEnergy = (*iList)->GetExcitationEnergy();
-      
-      if ( A > 1 && exEnergy > 0.1*eV ) // if so, photon-evaporation is applied
+      if ( A <= 1 )                                       { theResults.push_back(*iList); }
+      else if ((*iList)->GetExcitationEnergy() <= 0.1*eV) { theResults.push_back(*iList); }
+
+      else
         {
-          theExcitedNucleus = *(*iList);
-          theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus);
+          // photon-evaporation is applied
+          theTempResult = thePhotonEvaporation->BreakItUp(*(*iList));
           
+          G4bool deletePrimary = true;
+          G4int nsec = theTempResult->size();
           
           // if there is a gamma emission then
-          if (theTempResult->size() > 1)
+          if (nsec > 1)
             {
-              // first free the memory occupied by the previous state
-              delete (*iList);
+              G4FragmentVector::iterator j;
+              for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
+                {
+                  if((*j) == (*iList)) { deletePrimary = false; }
+                  A = static_cast<G4int>((*j)->GetA()+0.5);  // +0.5 to avoid bad truncation
+
+                  if(A <= 1)                                     { theResults.push_back(*j); }  // gamma, p, n
+                  else if((*j)->GetExcitationEnergy() <= 0.1*eV) { theResults.push_back(*j); }  // stable fragment
+                  else                                           { theEvapStableList.push_back(*j); }      
+                }
+            }
               
-              // and now add the final state from gamma emission to the end of
-              // theEvapStableList
-              for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin(); 
-                   ri != theTempResult->rend(); ++ri)
-                // reversed is applied in order to have residual nucleus in first position
+          else if(1 == nsec) 
+            {
+              G4FragmentVector::iterator j = theTempResult->begin();
+              if((*j) == (*iList)) { deletePrimary = false; }
+              // Let's create a G4Fragment pointer representing the gamma emmited
+              G4LorentzVector lv = (*j)->GetMomentum();
+              G4double Mass = (*j)->GetGroundStateMass();
+              G4double Ecm = lv.m();
+              if(Ecm - Mass > 0.1*eV) 
                 {
-                  
-#ifdef PRECOMPOUND_TEST
-                  if ((*ri)->GetA() == 0)
-                    (*ri)->SetCreatorModel(G4String("G4PhotonEvaporation"));
-                  else
-                    (*ri)->SetCreatorModel(G4String("ResidualNucleus"));
+                  G4ThreeVector bst = lv.boostVector();
+                  G4double GammaEnergy = 0.5*(Ecm - Mass)*(Ecm + Mass)/Ecm;
+                  G4double cosTheta = 1. - 2. * G4UniformRand(); 
+                  G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
+                  G4double phi = twopi * G4UniformRand();
+                  G4LorentzVector Gamma4P(GammaEnergy * sinTheta * std::cos(phi),
+                                          GammaEnergy * sinTheta * std::sin(phi),
+                                          GammaEnergy * cosTheta,
+                                          GammaEnergy);
+                  Gamma4P.boost(bst);  
+                  G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition());
+                  theResults.push_back(theHandlerPhoton); 
+              
+                  // And now we update momentum and energy for the nucleus
+                  lv -= Gamma4P;
+                  (*j)->SetMomentum(lv); // Now this fragment has been deexcited!
+                }
+              // we store the deexcited fragment 
+              theResults.push_back(*j);
+            }
+          if( deletePrimary ) { delete (*iList); }
+          delete theTempResult;
+        }  
+    } // end of photon-evaporation loop
+
+  //G4cout << "## After 3d step " << theEvapList.size() << " was evap;  "
+  //     << theEvapStableList.size() << " was photo-evap; " 
+  //     << theResults.size() << " results. " << G4endl; 
+    
+#ifdef debug
+  CheckConservation(theInitialState,*theResults);
 #endif
-                  
-                  theEvapStableList.push_back(*ri);
-                }
-              
-              // now we clean and remove the temporal vector
-              theTempResult->clear();
-              delete theTempResult;
-              
-            }
-          else  // if theTempResult->size() = 1
-            {
-              // if there is not any gamma emission from this excited fragment
-              // we have to emmit a gamma which forces the deexcitation
-              
-              // First I clean completely theTempResult
-              
-              for (G4FragmentVector::iterator j = theTempResult->begin();
-                   j != theTempResult->end(); j++)
-                {
-                  delete (*j);
-                }
-              
-              theTempResult->clear();
-              delete theTempResult;
-              
-#ifdef debugphoton
-              G4cout << "G4ExcitationHandler: Gamma Evaporation could not deexcite the nucleus: \n"
-                     << "-----------------------------------------------------------------------\n"
-                     << theExcitedNucleus << '\n'
-                     << "-----------------------------------------------------------------------\n";
-#endif
-              
-              // Let's create a G4Fragment pointer representing the gamma emmited
-              G4double GammaEnergy = (*iList)->GetExcitationEnergy();
-              G4double cosTheta = 1. - 2. * G4UniformRand();
-              G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
-              G4double phi = twopi * G4UniformRand();
-              G4ThreeVector GammaP(GammaEnergy * sinTheta * std::cos(phi),
-                                   GammaEnergy * sinTheta * std::sin(phi),
-                                   GammaEnergy * cosTheta );
-              G4LorentzVector Gamma4P(GammaP,GammaEnergy);
-              G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition());
-              
-              // And now we update momentum and energy for the nucleus
-              G4double Mass = (*iList)->GetGroundStateMass();
-              G4ThreeVector ResidualP((*iList)->GetMomentum().vect() - GammaP);
-              G4double ResidualE = std::sqrt(ResidualP*ResidualP + Mass*Mass);
-              G4LorentzVector Residual4P(ResidualP,ResidualE);
-              (*iList)->SetMomentum(Residual4P); // Now this fragment has been deexcited!
-              
-              // we store the deexcited fragment in theFinalStableList
-              theFinalStableList.push_back(*iList);
-              
-#ifdef PRECOMPOUND_TEST
-              theHandlerPhoton->SetCreatorModel("G4ExcitationHandler");
-#endif
-              
-              // Finally, we add theHandlerPhoton to theFinalStableList
-              theFinalStableList.push_back(theHandlerPhoton); 
-              
-#ifdef debugphoton
-              G4cout << "Emmited photon:\n"
-                     << theFinalStableList.back() << '\n'
-                     << "Residual nucleus after photon emission:\n"
-                     << *(*iList) << '\n'
-                     << "-----------------------------------------------------------------------\n";
-#endif
-              
-            }
+
+  G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
+
+  // MAC (24/07/08)
+  // To optimise the storing speed, we reserve space in memory for the vector
+  theReactionProductVector->reserve( theResults.size() );
+
+  G4int theFragmentA, theFragmentZ;
+  G4LorentzVector theFragmentMomentum;
+
+  std::list<G4Fragment*>::iterator i;
+  for (i = theResults.begin(); i != theResults.end(); ++i) 
+    {
+      theFragmentA = static_cast<G4int>((*i)->GetA());
+      theFragmentZ = static_cast<G4int>((*i)->GetZ());
+      theFragmentMomentum = (*i)->GetMomentum();
+      G4ParticleDefinition* theKindOfFragment = 0;
+      if (theFragmentA == 0 && theFragmentZ == 0) {       // photon
+        theKindOfFragment = G4Gamma::GammaDefinition();   
+      } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
+        theKindOfFragment = G4Neutron::NeutronDefinition();
+      } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
+        theKindOfFragment = G4Proton::ProtonDefinition();
+      } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
+        theKindOfFragment = G4Deuteron::DeuteronDefinition();
+      } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
+        theKindOfFragment = G4Triton::TritonDefinition();
+      } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
+        theKindOfFragment = G4He3::He3Definition();
+      } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
+        theKindOfFragment = G4Alpha::AlphaDefinition();;
+      } else {
+        theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ);
+      }
+      if (theKindOfFragment != 0) 
+        {
+          G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
+          theNew->SetMomentum(theFragmentMomentum.vect());
+          theNew->SetTotalEnergy(theFragmentMomentum.e());
+          theNew->SetFormationTime((*i)->GetCreationTime());
+          theReactionProductVector->push_back(theNew);
         }
-      else // case of a nucleon, gamma or very small excitation energy
-        {
-          // we don't have to do anything, just store the fragment in theFinalStableList
-          
-          theFinalStableList.push_back(*iList);
-          
-        } // A > 1 && exEnergy > 0.1*eV
-      
-    } // end of photon-evaporation loop
-  
-  
-  // The deexcitation from fragments inside theEvapStableList has been finished, so...
-  theEvapStableList.clear();
-  
-  // Now the final state is in theFinalStableList, and we have to send it to theResult vector
-  
-  theResult = new G4FragmentVector;
-  theResult->reserve( theFinalStableList.size() * sizeof(G4Fragment*) );
-  // We reserve enough memory to optimise the storing speed
-  
-  for (iList = theFinalStableList.begin(); iList != theFinalStableList.end(); iList++)
-    {
-      theResult->push_back(*iList);
+      delete (*i);
     }
-  
-  // After storing the final state , we can clear theFinalStableList
-  theFinalStableList.clear();
-  
-#ifdef debug
-  CheckConservation(theInitialState,theResult);
-#endif
-  
-  
-  // Change G4FragmentVector* to G4ReactionProductVector*
-  return Transform(theResult);
+
+  return theReactionProductVector;
 }