Index: trunk/source/processes/hadronic/models/de_excitation/History
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/History	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/History	(revision 1315)
@@ -14,4 +14,141 @@
      * Please list in reverse chronological order (last date on top)
      ---------------------------------------------------------------
+
+9 June 2010 Vladimir Ivanchenko (hadr-deex-V09-03-18)
+----------------------------------------------------------
+- G4Evaporation - fixed problem of isotope production for high Z
+                  fragments introduced in previous tags
+
+25 May 2010 Vladimir Ivanchenko (hadr-deex-V09-03-17)
+----------------------------------------------------------
+- G4E1Probability - added numerical protection to avoid division by zero
+
+19 May 2010 Vladimir Ivanchenko (hadr-deex-V09-03-16)
+----------------------------------------------------------
+- G4UnstableFragmentBreakUp - fix selection of decay channel by addition of check 
+                              on residual fragment Z and A (addressed
+                              problem reported stt for the tag 03-14)
+
+19 May 2010 Vladimir Ivanchenko (hadr-deex-V09-03-15)
+----------------------------------------------------------
+- G4UnstableFragmentBreakUp - fixed selection of decay channel
+- G4E1Probability - code cleanup and optimisation by usage of G4Pow, integer A 
+                    and introduction of const members
+- G4GEMProbability - code cleanup and optimisation by usage of G4Pow and integer Z,A
+- G4ExcitationHandler - forced  FermiBreakUp for A < 5
+- G4PhotonEvaporation - (F.Lei) added correction of electron state after emission
+- G4FermiFragmentsPool - JMQ fixed fragment 111
+
+11 May 2010 Vladimir Ivanchenko (hadr-deex-V09-03-14)
+----------------------------------------------------------
+- G4UnstableFragmentBreakUp - new class to decay exotic fragmnet (like 2n or 2p)
+- G4Evaporation - added call to G4UnstableFragmentBreakUp if natural abandances
+                  of cold fragment is zero; optimized logic of stopping of 
+                  evaporation loop
+- G4PhotonEvaporation - cleanup new methods EmittedFragment and BreakUpFragment
+- G4ExcitationHandler - use BreakUpFragment method for photon deexcitation
+
+10 May 2010 Vladimir Ivanchenko (hadr-deex-V09-03-13)
+----------------------------------------------------------
+- G4VGammaDeexcitation - take into account bounding energy of electron 
+                         in the case of electron emmision; fixed kinematic
+- G4DiscreteGammaTransition - Removed unphysical corretions of gamma energy; 
+                         fixed default particle as gamma; do not subtract 
+                         bounding energy in case of electron emmision
+- G4ExcitationHandler - allowed emmision e- instead of gamma due to internal conversion
+
+
+03 May 2010 Vladimir Ivanchenko (hadr-deex-V09-03-12)
+----------------------------------------------------------
+- G4Evaporation - improved condition how to stop deexcitation loop
+- G4VGammaDeexcitation - take into account electron mass in the case of
+                         electron emmision
+- G4PhotonEvaporation - improved printout
+- G4ExcitationHandler - disable MFM
+- G4StatMFFragment, G4CompetitiveFission, G4EvaporationProbability,
+  G4GEMProbability - use correct header files
+
+28 April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-11)
+----------------------------------------------------------
+- G4ExcitationHandler - (JQM+VI) add check on stability of primary;
+                                 do evaporation if FermiBreakUp or MFM
+                                 cannot deexcite a fragment; 
+	                         added SetParameters method
+- G4Evaporation - rewrite BreakUp method; added Initialise method, where setup
+                  all options and not at run time; added InitialiseEvaporation
+                  method to setup decay channels; changed order of decay
+                  channels - first photon evaporation, second - fision,
+                  after all other channels as before
+- virtual interfaces - moved constructors and destructors to source files
+
+26 April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-10)
+----------------------------------------------------------
+- G4FermiConfiguration - (JQM) parameter of Coulomb energy Kappa is changed 
+                          from 1 to 6 according to recommendation of original
+                          author of the model A. Botvina
+- G4FermiPhaseSpaceDecay - (JQM) improved model of sampling of kinetic energy;
+                         - (VI) cleanup the code by using G4Pow; moved constructor
+                            and destructor to the source 
+
+25 April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-09)
+----------------------------------------------------------
+- G4ExcitationHandler - (JQM+VI) apply FermiBreakUp to fragments with A>1 
+                        (before was A>4) in order to reduce number of 
+                        fake gamma produced in order deexcite light
+                        fragments; added parameter minExcitation = 1 keV
+- G4VEvaporationChannel, G4PhotonEvaporation - added 2 new virtual methods 
+                        EmittedFragment and BreakUpFragment
+
+23 April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-08)
+----------------------------------------------------------
+- G4VGammaDeexcitation - kinematic of 2-body decay rewritten
+- G4DiscreteGammaTransition, G4DiscreteGammaDeexcitation, 
+  G4ContinuumGammaDeexcitation - set energy tolerance 1 keV;
+                                 set destructors to be virtual
+
+21 April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-07)
+----------------------------------------------------------
+- G4FermiFragmentsPool - (JMQ) extended list of stable fragments
+- G4DiscreteGammaTransition - (JMQ) make transition depended on Z and A
+                                    (before was only Z) and added 
+                                    energy tolerance
+- G4ContinuumGammaDeexcitation - (JQM) more accurate Lorentz computations
+- G4VGammaDeexcitation - (JMQ) improved Lorentz computations
+
+20 April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-06)
+----------------------------------------------------------
+- G4GEMProbability - (JQM + VI) fixed numerical problem (division by zero)
+
+16 April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-05)
+----------------------------------------------------------
+- G4ExcitationHandler - enable Multi-Fragmentation model
+- G4StatMFMacroTemperature - cleanup logic of solving equation for
+           temperature; moved constructors and destructor to source
+
+09 April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-04)
+----------------------------------------------------------
+- G4ProtonEvaporationProbability, G4DeuteronEvaporationProbability,
+  G4TritonEvaporationProbability, G4He3EvaporationProbability,
+  G4AlphaEvaporationProbability - (JMQ) return back to published
+               variant OPT3 (Kalbach) parameterization of inverse 
+	       cross section 
+
+05 March 2010 Vladimir Ivanchenko (hadr-deex-V09-03-03)
+----------------------------------------------------------
+- G4Evaporation - set as a default variant evaporation combined 
+                  standard + GEM probabilities
+
+17 February 2009 Vladimir Ivanchenko (hadr-deex-V09-03-02)
+----------------------------------------------------------
+- G4ExcitationHandler - deactivate Multi-Fragmentation model
+
+05 February 2009 Vladimir Ivanchenko (hadr-deex-V09-03-01)
+----------------------------------------------------------
+- G4ExcitationHandler - activate FermiBreakUp and Multi-Fragmentation models
+
+18 January 2010  Vladimir Ivanchenko (hadr-deex-V09-03-00)
+----------------------------------------------------------
+- Move 9.3 version to the head of cvs for following files: G4FissionBarrier.hh,
+  G4FissionBarrier.cc, G4VGammaDeexcitation.cc, G4VGammaDeexcitation.hh,
 
 09 December 2009 Gunter Folger  (hadr-deex-V09-02-23)
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4Evaporation.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4Evaporation.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4Evaporation.hh	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Evaporation.hh,v 1.7 2009/07/27 10:32:05 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4Evaporation.hh,v 1.12 2010/05/11 11:34:09 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -32,6 +32,8 @@
 //
 //
-// Alex Howard - added protection for negative probabilities in the sum, 14/2/07
-// V.Ivanchenko - added Combined decay channels (default + GEM) 27/07/09
+// 14/02/2007 Alex Howard - added protection for negative probabilities in the sum 
+// 27/07/2009 V.Ivanchenko - added Combined decay channels (default + GEM) 
+// 11/05/2010 V.Ivanchenko - rewrited technical part do not "new" and "delete" 
+//                           of small objects
 
 #ifndef G4Evaporation_h
@@ -43,6 +45,8 @@
 #include "G4VEvaporationChannel.hh"
 #include "G4Fragment.hh"
+#include "G4UnstableFragmentBreakUp.hh"
 
 class G4VEvaporationFactory;
+class G4NistManager;
 
 //#define debug
@@ -52,11 +56,14 @@
 public:
   G4Evaporation();
-  G4Evaporation(std::vector<G4VEvaporationChannel*> * aChannelsVector) :
-    theChannels(aChannelsVector), theChannelFactory(0)
-  {};
+  G4Evaporation(std::vector<G4VEvaporationChannel*> * aChannelsVector);
 	 
   virtual ~G4Evaporation();
 
+  virtual void Initialise();
+
 private:
+
+  void InitialiseEvaporation();
+
   G4Evaporation(const G4Evaporation &right);
 
@@ -66,4 +73,5 @@
 
 public:
+
   G4FragmentVector * BreakItUp(const G4Fragment &theNucleus);
 
@@ -77,27 +85,12 @@
 #endif
 
-
   std::vector<G4VEvaporationChannel*> * theChannels;
+  std::vector<G4double>   probabilities;
   G4VEvaporationFactory * theChannelFactory;
-  
-  
-  class SumProbabilities : public std::binary_function<G4double,G4double,G4double>
-  {
-  public:
-    SumProbabilities() : total(0.0) {}
-    G4double operator() (G4double& /* probSoFar */, G4VEvaporationChannel*& frag)
-    { 
-      G4double temp_prob = frag->GetEmissionProbability();
-      if(temp_prob >= 0.0) total += temp_prob;
-      //      total += frag->GetEmissionProbability();
-      return total;
-    }
+  G4int nChannels;
+  G4double minExcitation;
+  G4NistManager* nist;
+  G4UnstableFragmentBreakUp unstableBreakUp;
     
-    G4double GetTotal() { return total; }
-  public:
-    G4double total;
-    
-  };
-
 };
 
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationFactory.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationFactory.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationFactory.hh	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4EvaporationFactory.hh,v 1.3 2006/06/29 20:09:55 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4EvaporationFactory.hh,v 1.4 2010/04/27 11:43:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -41,9 +41,9 @@
 {
 public:
-  G4EvaporationFactory() {};
-  virtual ~G4EvaporationFactory() {};
+  G4EvaporationFactory();
+  virtual ~G4EvaporationFactory();
 
 private:
-  G4EvaporationFactory(const G4EvaporationFactory & ) : G4VEvaporationFactory() {};
+  G4EvaporationFactory(const G4EvaporationFactory & );
   const G4EvaporationFactory & operator=(const G4EvaporationFactory & val);
   G4bool operator==(const G4EvaporationFactory & val) const;
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4VEvaporation.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4VEvaporation.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4VEvaporation.hh	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4VEvaporation.hh,v 1.4 2008/09/19 13:32:54 ahoward Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4VEvaporation.hh,v 1.6 2010/04/27 14:00:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -46,6 +46,6 @@
 {
 public:
-  G4VEvaporation() {};
-  virtual ~G4VEvaporation() {}; // *
+  G4VEvaporation();
+  virtual ~G4VEvaporation(); 
 
 private:  
@@ -57,5 +57,8 @@
   
 public:
+
   virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus) = 0;
+
+  virtual void Initialise();
 
   // for inverse cross section choice
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc	(revision 1315)
@@ -225,5 +225,8 @@
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-  if (cut < 0.) ecut2 = ecut - 2.;
+//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
+//  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) ecut2 = ecut;
   elab = K * FragmentA / ResidualA;
   sig = 0.;
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc	(revision 1315)
@@ -206,6 +206,7 @@
         
   //JMQ 13/02/09 increase of reduced radius to lower the barrier
-  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
-  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
+  //  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
+   ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
+
   ecsq = ec * ec;
   p = p0 + p1/ec + p2/ecsq;
@@ -225,5 +226,8 @@
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-  if (cut < 0.) ecut2 = ecut - 2.;
+//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
+//  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) ecut2 = ecut;
   elab = K * FragmentA / ResidualA;
   sig = 0.;
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Evaporation.cc,v 1.15 2009/09/16 15:32:25 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4Evaporation.cc,v 1.25 2010/06/09 11:56:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -38,5 +38,10 @@
 // superimposed Coulomb barrier (if useSICBis set true, by default is false) 
 //
-// V.Ivanchenko added G4EvaporationDefaultGEMFactory option, 27/07/09
+// V.Ivanchenko (27 July 2009) added G4EvaporationDefaultGEMFactory option
+// V.Ivanchenko (10 May 2010) rewrited BreakItUp method: do not make new/delete
+//                            photon channel first, fission second,
+//                            added G4UnstableFragmentBreakUp to decay 
+//                            unphysical fragments (like 2n or 2p),
+//                            use Z and A integer
 
 #include "G4Evaporation.hh"
@@ -45,17 +50,46 @@
 #include "G4EvaporationDefaultGEMFactory.hh"
 #include "G4HadronicException.hh"
-#include <numeric>
+#include "G4NistManager.hh"
 
 G4Evaporation::G4Evaporation() 
 {
-  theChannelFactory = new G4EvaporationFactory();
-  //theChannelFactory = new G4EvaporationDefaultGEMFactory();
-  theChannels = theChannelFactory->GetChannel();
+  //theChannelFactory = new G4EvaporationFactory();
+  theChannelFactory = new G4EvaporationDefaultGEMFactory();
+  InitialiseEvaporation();
+}
+
+G4Evaporation::G4Evaporation(std::vector<G4VEvaporationChannel*> * aChannelsVector) 
+  : theChannels(aChannelsVector), theChannelFactory(0), nChannels(0)
+{
+  InitialiseEvaporation();
 }
 
 G4Evaporation::~G4Evaporation()
 {
-  if (theChannels != 0) theChannels = 0;
-  if (theChannelFactory != 0) delete theChannelFactory;
+  if (theChannels != 0) { theChannels = 0; }
+  if (theChannelFactory != 0) { delete theChannelFactory; }
+}
+
+void G4Evaporation::InitialiseEvaporation()
+{
+  nist = G4NistManager::Instance();
+  minExcitation = CLHEP::keV;
+  if(theChannelFactory) { theChannels = theChannelFactory->GetChannel(); }
+  nChannels = theChannels->size();
+  probabilities.resize(nChannels, 0.0);
+  Initialise();
+}
+
+void G4Evaporation::Initialise()
+{
+  // loop over evaporation channels
+  std::vector<G4VEvaporationChannel*>::iterator i;
+  for (i=theChannels->begin(); i != theChannels->end(); ++i) 
+    {
+      // for inverse cross section choice
+      (*i)->SetOPTxs(OPTxs);
+      // for superimposed Coulomb Barrier for inverse cross sections
+      (*i)->UseSICB(useSICB);
+    }
 }
 
@@ -64,5 +98,5 @@
   if (theChannelFactory != 0) delete theChannelFactory;
   theChannelFactory = new G4EvaporationFactory();
-  theChannels = theChannelFactory->GetChannel();
+  InitialiseEvaporation();
 }
 
@@ -71,5 +105,5 @@
   if (theChannelFactory != 0) delete theChannelFactory;
   theChannelFactory = new G4EvaporationGEMFactory();
-  theChannels = theChannelFactory->GetChannel();
+  InitialiseEvaporation();
 }
 
@@ -78,11 +112,176 @@
   if (theChannelFactory != 0) delete theChannelFactory;
   theChannelFactory = new G4EvaporationDefaultGEMFactory();
-  theChannels = theChannelFactory->GetChannel();
-}
-
+  InitialiseEvaporation();
+}
 
 G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
 {
-    G4FragmentVector * theResult = new G4FragmentVector;
+  G4FragmentVector * theResult = new G4FragmentVector;
+  G4FragmentVector * theTempResult;
+
+  // The residual nucleus (after evaporation of each fragment)
+  G4Fragment* theResidualNucleus = new G4Fragment(theNucleus);
+
+  G4double totprob, prob, oldprob = 0.0;
+  G4int maxchannel, i;
+
+  G4int Amax = theResidualNucleus->GetA_asInt();
+
+  // Starts loop over evaporated particles, loop is limited by number
+  // of nucleons
+  for(G4int ia=0; ia<Amax; ++ia) {
+ 
+    // g,n,p - evaporation is finished
+    G4int A = theResidualNucleus->GetA_asInt();
+    if(1 >= A) {
+      theResult->push_back(theResidualNucleus);
+      return theResult;
+    }
+
+    // check if it is stable, then finish evaporation
+    G4int Z = theResidualNucleus->GetZ_asInt();
+    G4double abun = nist->GetIsotopeAbundance(Z, A); 
+    /*
+    G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z
+	   << " A= " << A << " Eex(MeV)= " 
+	   << theResidualNucleus->GetExcitationEnergy()
+	   << " aban= " << abun << G4endl;
+    */
+    if(theResidualNucleus->GetExcitationEnergy() <= minExcitation && 
+       (abun > 0.0)) 
+      {
+	theResult->push_back(theResidualNucleus);
+	return theResult;
+      }
+
+    totprob = 0.0;
+    maxchannel = nChannels;
+
+    //G4cout << "### Evaporation loop #" << ia 
+    //	   << "  Fragment: " << theResidualNucleus << G4endl;
+
+    // loop over evaporation channels
+    for(i=0; i<nChannels; ++i) {
+      (*theChannels)[i]->Initialize(*theResidualNucleus);
+      prob = (*theChannels)[i]->GetEmissionProbability();
+      //G4cout << "  Channel# " << i << "  prob= " << prob << G4endl; 
+
+      //if(0 == i && 0.0 == abun) { prob = 0.0; }
+
+      totprob += prob;
+      probabilities[i] = totprob;
+      // if two recent probabilities are near zero stop computations
+      if(i>=8) {
+	if(prob <= totprob*1.e-8 && oldprob <= totprob*1.e-8) {
+	  maxchannel = i+1; 
+	  break;
+	}
+      }
+      oldprob = prob;
+    }
+
+    // photon evaporation in the case of no other channels available
+    // do evaporation chain and reset total probability
+    if(0.0 < totprob && probabilities[0] == totprob) {
+      theTempResult = (*theChannels)[0]->BreakUpFragment(theResidualNucleus);
+      if(theTempResult) {
+	size_t nsec = theTempResult->size();
+	for(size_t j=0; j<nsec; ++j) {
+	  theResult->push_back((*theTempResult)[j]);
+	}
+        delete theTempResult;
+      }
+      totprob = 0.0;
+    }
+
+    // stable fragnent - evaporation is finished
+    if(0.0 == totprob) {
+
+      // if fragment is exotic, then try to decay it
+      if(0.0 == abun && Z < 20) {
+        //G4cout << "$$$ Decay exotic fragment" << G4endl;
+	theTempResult = unstableBreakUp.BreakUpFragment(theResidualNucleus);
+	if(theTempResult) {
+	  size_t nsec = theTempResult->size();
+	  for(size_t j=0; j<nsec; ++j) {
+	    theResult->push_back((*theTempResult)[j]);
+	  }
+	  delete theTempResult;
+	}
+      }
+
+      // save residual fragment
+      theResult->push_back(theResidualNucleus);
+      return theResult;
+    }
+
+
+    // select channel
+    totprob *= G4UniformRand();
+    // loop over evaporation channels
+    for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
+
+    // this should not happen
+    if(i >= nChannels) { i = nChannels - 1; }
+
+
+    // single photon evaporation, primary pointer is kept
+    if(0 == i) {
+      G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus);
+      if(gamma) { theResult->push_back(gamma); }
+
+      // fission, return results to the main loop if fission is succesful
+    } else if(1 == i) {
+      theTempResult = (*theChannels)[1]->BreakUp(*theResidualNucleus);
+      if(theTempResult) {
+	size_t nsec = theTempResult->size();
+        G4bool deletePrimary = true;
+	for(size_t j=0; j<nsec; ++j) {
+          if(theResidualNucleus == (*theTempResult)[j]) { deletePrimary = false; }
+	  theResult->push_back((*theTempResult)[j]);
+	}
+        if(deletePrimary) { delete theResidualNucleus; }
+        delete theTempResult;
+	return theResult;
+      }
+
+      // other channels
+    } else {
+      theTempResult = (*theChannels)[i]->BreakUp(*theResidualNucleus);
+      if(theTempResult) {
+	size_t nsec = theTempResult->size();
+        if(nsec > 0) {
+          --nsec;
+	  for(size_t j=0; j<nsec; ++j) {
+	    theResult->push_back((*theTempResult)[j]);
+	  }
+	  // if the residual change its pointer 
+	  // then delete previous residual fragment and update to the new
+          if(theResidualNucleus != (*theTempResult)[nsec] ) { 
+	    delete theResidualNucleus; 
+	    theResidualNucleus = (*theTempResult)[nsec];
+	  }
+	}
+        delete theTempResult;
+      }
+    }
+  }
+
+  // loop is stopped, save residual
+  theResult->push_back(theResidualNucleus);
+  
+#ifdef debug
+  G4cout << "======== Evaporation Conservation Test ===========\n"
+	 << "==================================================\n";
+  CheckConservation(theNucleus,theResult);
+  G4cout << "==================================================\n";
+#endif
+  return theResult;
+}
+
+/*
+G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
+{
+  G4FragmentVector * theResult = new G4FragmentVector;
 
     // CHECK that Excitation Energy != 0
@@ -233,5 +432,5 @@
     return theResult;
 }
-
+*/
 
 
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationDefaultGEMFactory.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationDefaultGEMFactory.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationDefaultGEMFactory.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4EvaporationDefaultGEMFactory.cc,v 1.1 2009/07/27 10:20:13 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4EvaporationDefaultGEMFactory.cc,v 1.2 2010/04/27 11:43:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -122,4 +122,7 @@
   theChannel->reserve(68);
 
+  theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
+  theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
+
   // JMQ 220709 standard particle evaporation channels (Z<3,A<5)
   theChannel->push_back( new G4NeutronEvaporationChannel() );  // n
@@ -192,7 +195,4 @@
   theChannel->push_back( new G4Mg28GEMChannel() );     // Mg28
 
-  theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
-  theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
-
   return theChannel;
 
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4EvaporationFactory.cc,v 1.4 2006/06/29 20:10:29 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4EvaporationFactory.cc,v 1.5 2010/04/27 11:43:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -32,5 +32,4 @@
 
 #include "G4EvaporationFactory.hh"
-
 
 #include "G4NeutronEvaporationChannel.hh"
@@ -44,27 +43,9 @@
 #include "G4PhotonEvaporation.hh"
 
+G4EvaporationFactory::G4EvaporationFactory()
+{}
 
-const G4EvaporationFactory & 
-G4EvaporationFactory::operator=(const G4EvaporationFactory & )
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationFactory::operator= meant to not be accessable.");
-  return *this;
-}
-
-G4bool 
-G4EvaporationFactory::operator==(const G4EvaporationFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationFactory::operator== meant to not be accessable.");
-  return false;
-}
-
-G4bool 
-G4EvaporationFactory::operator!=(const G4EvaporationFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationFactory::operator!= meant to not be accessable.");
-  return true;
-}
-
-
+G4EvaporationFactory::~G4EvaporationFactory()
+{}
 
 std::vector<G4VEvaporationChannel*> * 
@@ -75,4 +56,7 @@
   theChannel->reserve(8);
 
+  theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
+  theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
+
   theChannel->push_back( new G4NeutronEvaporationChannel() );  // n
   theChannel->push_back( new G4ProtonEvaporationChannel() );   // p
@@ -82,6 +66,4 @@
   theChannel->push_back( new G4AlphaEvaporationChannel() );    // Alpha
 
-  theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
-  theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
 
   return theChannel;
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc	(revision 1315)
@@ -41,5 +41,6 @@
 #include "G4EvaporationProbability.hh"
 #include "G4PairingCorrection.hh"
-
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
 
 
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc	(revision 1315)
@@ -221,5 +221,8 @@
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-  if (cut < 0.) ecut2 = ecut - 2.;
+//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
+//  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) ecut2 = ecut;
   elab = K * FragmentA / ResidualA;
   sig = 0.;
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc	(revision 1315)
@@ -91,8 +91,9 @@
 
 ///////////////////////////////////////////////////////////////////////////////////
-//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
+//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections for protons
 //OPT=0 Dostrovski's parameterization
-//OPT=1,2 Chatterjee's paramaterization 
-//OPT=3,4 Kalbach's parameterization 
+//OPT=1 Chatterjee's parameterization 
+//OPT=2,4 Wellisch's parameterization 
+//OPT=3 Kalbach's parameterization
 // 
 G4double G4ProtonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
@@ -110,5 +111,4 @@
   FragmentAthrd=std::pow(FragmentA,0.33333);
   U=fragment.GetExcitationEnergy();
-
 
   if (OPTxs==0) {std::ostringstream errOs;
@@ -223,6 +223,5 @@
 //     ** p from  becchetti and greenlees (but modified with sub-barrier
 //     ** correction function and xp2 changed from -449)
-// JMQ (june 2008) : refinement of proton cross section for light systems
-//
+
   G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
   G4double p, p0, p1, p2;
@@ -279,5 +278,8 @@
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-  if (cut < 0.) ecut2 = ecut - 2.;
+//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
+//  if (cut < 0.) ecut2 = ecut - 2.;
+ if (cut < 0.) ecut2 = ecut;
   elab = K * FragmentA / ResidualA;
   sig = 0.;
@@ -286,18 +288,8 @@
     signor2 = (ec-elab-c) / w;
     signor2 = 1. + std::exp(signor2);
-    sig = sig / signor2;
-    // signor2 is empirical global corr'n at low elab for protons in PRECO, not enough for light targets
-    //  refinement for proton cross section
-    if (ResidualZ<=26)
-      sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));      
-                       }              //end for E<Ec
+    sig = sig / signor2;     
+                       }       //end for E<=Ec
   else            {           //start for  E>Ec
     sig = (landa*elab+mu+nu/elab) * signor;
-
-    //  refinement for proton cross section
-    if ( ResidualZ<=26 && elab <=(afit*ResidualZ+bfit)*ec) 
-           sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));
-                                                                               
-//
     geom = 0.;
 
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc	(revision 1315)
@@ -219,5 +219,8 @@
   ecut = (ecut-a) / (p+p);
   ecut2 = ecut;
-  if (cut < 0.) ecut2 = ecut - 2.;
+//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
+//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
+//  if (cut < 0.) ecut2 = ecut - 2.;
+  if (cut < 0.) ecut2 = ecut;
   elab = K * FragmentA / ResidualA;
   sig = 0.;
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc	(revision 1315)
@@ -24,7 +24,6 @@
 // ********************************************************************
 //
-//
-// $Id: G4VEvaporation.cc,v 1.5 2006/06/29 20:10:49 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4VEvaporation.cc,v 1.7 2010/04/27 14:00:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -32,32 +31,14 @@
 //
 
+#include "G4VEvaporation.hh"
 
-#include "G4VEvaporation.hh"
-#include "G4HadronicException.hh"
+G4VEvaporation::G4VEvaporation() 
+{}
 
-G4VEvaporation::G4VEvaporation(const G4VEvaporation &)
-{
- throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporation::copy_constructor meant to not be accessable");
-}
+G4VEvaporation::~G4VEvaporation() 
+{}
 
-
-
-
-const G4VEvaporation & G4VEvaporation::operator=(const G4VEvaporation &)
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporation::operator= meant to not be accessable");
-  return *this;
-}
-
-
-G4bool G4VEvaporation::operator==(const G4VEvaporation &) const
-{
-  return false;
-}
-
-G4bool G4VEvaporation::operator!=(const G4VEvaporation &) const
-{
-  return true;
-}
+void G4VEvaporation::Initialise()
+{}
 
 
Index: trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/include/G4FermiPhaseSpaceDecay.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/include/G4FermiPhaseSpaceDecay.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/include/G4FermiPhaseSpaceDecay.hh	(revision 1315)
@@ -33,4 +33,6 @@
 #include "G4LorentzVector.hh"
 #include "G4ParticleMomentum.hh"
+#include "Randomize.hh"
+#include "G4Pow.hh"
 #include <CLHEP/Random/RandGamma.h>
 
@@ -41,6 +43,7 @@
 {
 public:
-  inline G4FermiPhaseSpaceDecay();
-  inline ~G4FermiPhaseSpaceDecay();
+
+  G4FermiPhaseSpaceDecay();
+  ~G4FermiPhaseSpaceDecay();
   
   inline std::vector<G4LorentzVector*> * 
@@ -48,8 +51,9 @@
 
 private:
-  inline G4FermiPhaseSpaceDecay(const G4FermiPhaseSpaceDecay&);
-  inline const G4FermiPhaseSpaceDecay & operator=(const G4FermiPhaseSpaceDecay &); 
-  inline G4bool operator==(const G4FermiPhaseSpaceDecay&);
-  inline G4bool operator!=(const G4FermiPhaseSpaceDecay&);
+
+  G4FermiPhaseSpaceDecay(const G4FermiPhaseSpaceDecay&);
+  const G4FermiPhaseSpaceDecay & operator=(const G4FermiPhaseSpaceDecay &); 
+  G4bool operator==(const G4FermiPhaseSpaceDecay&);
+  G4bool operator!=(const G4FermiPhaseSpaceDecay&);
 
   inline G4double PtwoBody(G4double E, G4double P1, G4double P2) const;
@@ -69,5 +73,47 @@
 };
 
-#include "G4FermiPhaseSpaceDecay.icc"
+inline G4double 
+G4FermiPhaseSpaceDecay::PtwoBody(G4double E, G4double P1, G4double P2) const
+{
+  G4double res = -1.0;
+  G4double P = (E+P1+P2)*(E+P1-P2)*(E-P1+P2)*(E-P1-P2)/(4.0*E*E);
+  if (P>0.0) { res = std::sqrt(P); }
+  return res;
+}
+
+inline std::vector<G4LorentzVector*> * G4FermiPhaseSpaceDecay::
+Decay(const G4double parent_mass, const std::vector<G4double>& fragment_masses) const
+{
+  return KopylovNBodyDecay(parent_mass,fragment_masses);
+}
+
+inline G4double 
+G4FermiPhaseSpaceDecay::BetaKopylov(const G4int K) const
+{
+  //JMQ 250410 old algorithm has been commented
+  // Notice that alpha > beta always
+  // const G4double beta = 1.5;
+  // G4double alpha = 1.5*(K-1);
+  // G4double Y1 = CLHEP::RandGamma::shoot(alpha,1);
+  // G4double Y2 = CLHEP::RandGamma::shoot(beta,1);
+  
+  // return Y1/(Y1+Y2);
+
+  G4Pow* g4pow = G4Pow::GetInstance();
+  G4int N = 3*K - 5;
+  G4double xN = G4double(N);
+  G4double F;
+  //G4double Fmax = std::pow((3.*K-5.)/(3.*K-4.),(3.*K-5.)/2.)*std::sqrt(1-((3.*K-5.)/(3.*K-4.))); 
+  // VI variant
+  G4double Fmax = std::sqrt(g4pow->powZ(N, xN/(xN + 1))/(xN + 1)); 
+  G4double chi;
+  do 
+    {
+      chi = G4UniformRand();
+      F = std::sqrt(g4pow->powZ(N, chi)*(1-chi));      
+    } while ( Fmax*G4UniformRand() > F);  
+  return chi;
+
+}
 
 #endif
Index: trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiConfiguration.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiConfiguration.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiConfiguration.cc	(revision 1315)
@@ -25,12 +25,12 @@
 //
 //
-// $Id: G4FermiConfiguration.cc,v 1.12 2009/12/16 17:51:09 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4FermiConfiguration.cc,v 1.13 2010/04/26 11:14:28 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
 // by V. Lara (Nov 1998)
 //
-// J.M.Quesada (12 October 2009) new implementation of Gamma function in configuration weight 
-
+// J. M. Quesada (12 October 2009) new implementation of Gamma function in configuration weight 
+// J. M. Quesada (09 March 2010) Kappa is set to 6.
 
 #include "G4FermiConfiguration.hh"
@@ -40,5 +40,8 @@
 // Kappa = V/V_0 it is used in calculation of Coulomb energy
 // Kappa is adimensional
-const G4double G4FermiConfiguration::Kappa = 1.0;
+// JMQ 090310 according to model developer (A. Botvina) no theoretical constraint for kappa below 10
+// kappa values  larger than 1 seem to provide better results. 6 is a good choice  
+// const G4double G4FermiConfiguration::Kappa = 1.0;
+const G4double G4FermiConfiguration::Kappa = 6.0;
 
 // r0 is the nuclear radius
Index: trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiFragmentsPool.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiFragmentsPool.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiFragmentsPool.cc	(revision 1315)
@@ -28,10 +28,12 @@
 // by V. Lara
 //
-// J.M.Quesada (JMQ) July 2009, bug fixed in excitation energies: 
+// J.M.Quesada,  July 2009, bug fixed in excitation energies: 
 // ALL of them are in MeV instead of keV (as they were expressed previously)
 // source:  http://www.nndc.bnl.gov/chart
 // Unknown excitation energies in He5  and Li5 have been suppressed
 // Long lived levels (half-lives of the order ps-fs have been included)   
-
+//
+// J. M. Quesada,  April 2010: excitation energies according to tabulated values 
+// in PhotonEvaporatoion2.0. Fake photons eliminated. 
 
 #include "G4FermiFragmentsPool.hh"
@@ -51,140 +53,215 @@
   static const G4StableFermiFragment Fragment04(  3, 2,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment04);
   static const G4StableFermiFragment Fragment05(  4, 2,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment05);
-//JMQ 30/06/09 unknown levels have been supressed
+  //JMQ 30/06/09 unknown levels have been supressed
   static const G4He5FermiFragment    Fragment06(  5, 2,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment06);// He5
   static const G4Li5FermiFragment    Fragment07(  5, 3,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment07);// Li5 
   static const G4StableFermiFragment Fragment08(  6, 2,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment08);
   static const G4StableFermiFragment Fragment09(  6, 3,  3,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment09);
-  static const G4StableFermiFragment Fragment10(  6, 3,  1,  3.56*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment10);
+  //  static const G4StableFermiFragment Fragment10(  6, 3,  1,  3.56*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment10);
+  static const G4StableFermiFragment Fragment10(  6, 3,  1,  3.562880*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment10);
   static const G4StableFermiFragment Fragment11(  7, 3,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment11);
-  static const G4StableFermiFragment Fragment12(  7, 3,  2,  0.48*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment12);
+  //  static const G4StableFermiFragment Fragment12(  7, 3,  2,  0.48*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment12);
+  static const G4StableFermiFragment Fragment12(  7, 3,  2,  0.4776120*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment12);
   static const G4StableFermiFragment Fragment13(  7, 4,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment13);
-  static const G4StableFermiFragment Fragment14(  7, 4,  2,  0.43*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment14);
+  //  static const G4StableFermiFragment Fragment14(  7, 4,  2,  0.43*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment14);
+  static const G4StableFermiFragment Fragment14(  7, 4,  2,  0.4290800*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment14);
   static const G4StableFermiFragment Fragment15(  8, 3,  5,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment15);
-  static const G4StableFermiFragment Fragment16(  8, 3,  3,  0.98*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment16);
+  //  static const G4StableFermiFragment Fragment16(  8, 3,  3,  0.98*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment16);
+  static const G4StableFermiFragment Fragment16(  8, 3,  3,  0.9808000*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment16);
   static const G4Be8FermiFragment    Fragment17(  8, 4,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment17); // Be8
   static const G4StableFermiFragment Fragment18(  9, 4,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment18);
   static const G4B9FermiFragment     Fragment19(  9, 5,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment19); // B9  
   static const G4StableFermiFragment Fragment20( 10, 4,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment20);
-  static const G4StableFermiFragment Fragment21( 10, 4,  5,  3.37*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment21);
-  static const G4StableFermiFragment Fragment22( 10, 4,  8,  5.96*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment22);
-  static const G4StableFermiFragment Fragment23( 10, 4,  1,  6.18*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment23);
-  static const G4StableFermiFragment Fragment24( 10, 4,  5,  6.26*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment24);
+  //  static const G4StableFermiFragment Fragment21( 10, 4,  5,  3.37*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment21);
+  static const G4StableFermiFragment Fragment21( 10, 4,  5,  3.368030*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment21);
+//  static const G4StableFermiFragment Fragment22( 10, 4,  8,  5.96*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment22);
+  static const G4StableFermiFragment Fragment22( 10, 4,  8,  5.958390*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment22);
+  //  static const G4StableFermiFragment Fragment23( 10, 4,  1,  6.18*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment23);
+  static const G4StableFermiFragment Fragment23( 10, 4,  1,  6.179300*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment23);
+  //  static const G4StableFermiFragment Fragment24( 10, 4,  5,  6.26*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment24);
+  static const G4StableFermiFragment Fragment24( 10, 4,  5,  6.263300*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment24);
   static const G4StableFermiFragment Fragment25( 10, 5,  7,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment25);
-  static const G4StableFermiFragment Fragment26( 10, 5,  3,  0.72*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment26);
-  static const G4StableFermiFragment Fragment27( 10, 5,  1,  1.74*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment27);
-  static const G4StableFermiFragment Fragment28( 10, 5,  3,  2.15*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment28);
-  static const G4StableFermiFragment Fragment29( 10, 5,  5,  3.59*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment29);
-  
+  //  static const G4StableFermiFragment Fragment26( 10, 5,  3,  0.72*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment26);
+  static const G4StableFermiFragment Fragment26( 10, 5,  3,  0.7183500*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment26);
+  //  static const G4StableFermiFragment Fragment27( 10, 5,  1,  1.74*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment27);
+  static const G4StableFermiFragment Fragment27( 10, 5,  1,  1.740150*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment27);
+  //  static const G4StableFermiFragment Fragment28( 10, 5,  3,  2.15*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment28);
+  static const G4StableFermiFragment Fragment28( 10, 5,  3,  2.154300*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment28);
+  //  static const G4StableFermiFragment Fragment29( 10, 5,  5,  3.59*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment29);
+  static const G4StableFermiFragment Fragment29( 10, 5,  5,  3.587100*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment29);  
   static const G4StableFermiFragment Fragment30( 10, 6,  3,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment30);
-  static const G4StableFermiFragment Fragment31( 10, 6,  5,  3.35*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment31);
+  //  static const G4StableFermiFragment Fragment31( 10, 6,  5,  3.35*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment31);
+  static const G4StableFermiFragment Fragment31( 10, 6,  5,  3.353600*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment31);
   static const G4StableFermiFragment Fragment32( 11, 5,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment32);
-  static const G4StableFermiFragment Fragment33( 11, 5,  2,  2.13*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment33);
-  static const G4StableFermiFragment Fragment34( 11, 5,  6,  4.44*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment34);
-  static const G4StableFermiFragment Fragment35( 11, 5,  4,  5.02*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment35);
-  static const G4StableFermiFragment Fragment36( 11, 5, 10,  6.76*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment36);
-  static const G4StableFermiFragment Fragment37( 11, 5,  6,  7.29*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment37);
-  static const G4StableFermiFragment Fragment38( 11, 5,  4,  7.98*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment38);
-  static const G4StableFermiFragment Fragment39( 11, 5,  6,  8.56*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment39);
-  
-  static const G4StableFermiFragment Fragment40( 11, 6,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment40);
-  static const G4StableFermiFragment Fragment41( 11, 6,  2,  2.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment41);
-  static const G4StableFermiFragment Fragment42( 11, 6,  6,  4.32*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment42);
-  static const G4StableFermiFragment Fragment43( 11, 6,  4,  4.80*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment43);
-  static const G4StableFermiFragment Fragment44( 11, 6,  2,  6.34*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment44);
-  static const G4StableFermiFragment Fragment45( 11, 6,  8,  6.48*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment45);
-  static const G4StableFermiFragment Fragment46( 11, 6,  6,  6.90*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment46);
-  static const G4StableFermiFragment Fragment47( 11, 6,  4,  7.50*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment47);
-  static const G4StableFermiFragment Fragment48( 11, 6,  4,  8.10*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment48);
-  static const G4StableFermiFragment Fragment49( 11, 6,  6,  8.42*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment49);
-
-  static const G4StableFermiFragment Fragment50( 12, 5,  3,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment50);
-  static const G4StableFermiFragment Fragment51( 12, 5,  5,  0.95*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment51);
-  static const G4StableFermiFragment Fragment52( 12, 5,  5,  1.67*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment52);
-  static const G4StableFermiFragment Fragment53( 12, 5,  4,  2.65*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment53);
-  static const G4StableFermiFragment Fragment54( 12, 6,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment54);
-  static const G4StableFermiFragment Fragment55( 12, 6,  5,  4.44*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment55);
-  static const G4StableFermiFragment Fragment56( 13, 6,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment56);
-  static const G4StableFermiFragment Fragment57( 13, 6,  2,  3.09*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment57);
-  static const G4StableFermiFragment Fragment58( 13, 6,  4,  3.68*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment58);
-  
-  static const G4StableFermiFragment Fragment59( 13, 6,  6,  3.85*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment59);
-  static const G4StableFermiFragment Fragment60( 13, 7,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment60);
-  static const G4StableFermiFragment Fragment61( 14, 6,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment61);
-  static const G4StableFermiFragment Fragment62( 14, 6,  3,  6.09*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment62);
-// JMQ 010709 corrected excitation energies for 64-66, according to http://www.nndc.bnl.gov/chart
-  static const G4StableFermiFragment Fragment63( 14, 6,  1,  6.59*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment63);
-  static const G4StableFermiFragment Fragment64( 14, 6,  7,  6.73*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment64);
-  static const G4StableFermiFragment Fragment65( 14, 6,  1,  6.90*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment65);
-  static const G4StableFermiFragment Fragment66( 14, 6,  5,  7.01*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment66);
-  static const G4StableFermiFragment Fragment67( 14, 6,  5,  7.34*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment67);
-
-  static const G4StableFermiFragment Fragment68( 14, 7,  3,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment68);
-  static const G4StableFermiFragment Fragment69( 14, 7,  1,  2.31*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment69);
-  static const G4StableFermiFragment Fragment70( 14, 7,  3,  3.95*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment70);
-  
-  static const G4StableFermiFragment Fragment71( 14, 7,  1,  4.92*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment71);
-  static const G4StableFermiFragment Fragment72( 14, 7,  5,  5.11*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment72);
-  static const G4StableFermiFragment Fragment73( 14, 7,  3,  5.69*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment73);
-  static const G4StableFermiFragment Fragment74( 14, 7,  7,  5.83*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment74);
-  static const G4StableFermiFragment Fragment75( 14, 7,  3,  6.20*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment75);
-  static const G4StableFermiFragment Fragment76( 14, 7,  7,  6.44*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment76);
-  static const G4StableFermiFragment Fragment77( 14, 7,  5,  7.03*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment77);
-  static const G4StableFermiFragment Fragment78( 15, 7,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment78);
-// JMQ 010709 two very close levels instead of only one, with their own spins
-  static const G4StableFermiFragment Fragment79( 15, 7,  6,  5.27*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment79);
-  static const G4StableFermiFragment Fragment80( 15, 7,  2,  5.30*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment80);
-  static const G4StableFermiFragment Fragment81( 15, 7,  4,  6.32*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment81);
-//JMQ 010709 new level and corrected energy and spins
-  static const G4StableFermiFragment Fragment82( 15, 7,  6,  7.15*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment82);
-  static const G4StableFermiFragment Fragment83( 15, 7,  4,  7.30*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment83);
-  static const G4StableFermiFragment Fragment84( 15, 7,  8,  7.57*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment84);
-  static const G4StableFermiFragment Fragment85( 15, 7,  2,  8.31*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment85);
-  static const G4StableFermiFragment Fragment86( 15, 7,  4,  8.57*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment86);
-  static const G4StableFermiFragment Fragment87( 15, 7,  2,  9.05*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment87);
-//JMQ 010709 new levels for N15
-  static const G4StableFermiFragment Fragment88( 15, 7,  4,  9.151*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment88);
-  static const G4StableFermiFragment Fragment89( 15, 7,  6,  9.154*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment89);
-  static const G4StableFermiFragment Fragment90( 15, 7,  2,  9.22*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment90);
-  static const G4StableFermiFragment Fragment91( 15, 7,  6,  9.76*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment91);
-  static const G4StableFermiFragment Fragment92( 15, 7,  8,  9.83*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment92);
-  static const G4StableFermiFragment Fragment93( 15, 7,  4,  9.93*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment93);
-  static const G4StableFermiFragment Fragment94( 15, 7,  4, 10.07*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment94);
-
-
-  static const G4StableFermiFragment Fragment95( 15, 8,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment95);
-//JMQ 010709 new level and spins
-  static const G4StableFermiFragment Fragment96( 15, 8,  2,  5.18*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment96);
-  static const G4StableFermiFragment Fragment97( 15, 8,  6,  5.24*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment97);
-  static const G4StableFermiFragment Fragment98( 15, 8,  4,  6.18*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment98);  
-  static const G4StableFermiFragment Fragment99( 15, 8,  4,  6.79*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment99);
-  static const G4StableFermiFragment Fragment100( 15, 8,  6,  6.86*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment100);
-  static const G4StableFermiFragment Fragment101( 15, 8,  8,  7.28*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment101);
-
-
-  static const G4StableFermiFragment Fragment102( 16, 7,  5,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment102);
-  static const G4StableFermiFragment Fragment103( 16, 7,  1,  0.12*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment103);
-  static const G4StableFermiFragment Fragment104( 16, 7,  7,  0.30*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment104);
-  static const G4StableFermiFragment Fragment105( 16, 7,  3,  0.40*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment105);
-
-//JMQ 010709   some energies and spins have been changed 
-  static const G4StableFermiFragment Fragment106( 16, 8,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment106);
-  static const G4StableFermiFragment Fragment107( 16, 8,  1,  6.05*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment107);
-  static const G4StableFermiFragment Fragment108( 16, 8,  7,  6.13*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment108);
-  static const G4StableFermiFragment Fragment109( 16, 8,  5,  6.92*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment109);
-  static const G4StableFermiFragment Fragment110( 16, 8,  3,  7.12*MeV ); if(MapIsEmpty) fragment_pool.push_back
-(&Fragment110);
+//  static const G4StableFermiFragment Fragment33( 11, 5,  2,  2.13*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment33);
+  static const G4StableFermiFragment Fragment33( 11, 5,  2,  2.124693*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment33);
+  //  static const G4StableFermiFragment Fragment34( 11, 5,  6,  4.44*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment34);
+  static const G4StableFermiFragment Fragment34( 11, 5,  6,  4.444890*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment34);
+  //  static const G4StableFermiFragment Fragment35( 11, 5,  4,  5.02*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment35);
+  static const G4StableFermiFragment Fragment35( 11, 5,  4,  5.020310*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment35);
+//  static const G4StableFermiFragment Fragment36( 11, 5, 10,  6.76*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment36);
+  static const G4StableFermiFragment Fragment36( 11, 5, 8,  6.742900*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment36);
+  //JMQ 190410 new level, fragment numbering shifted accordingly from here onwards
+  static const G4StableFermiFragment Fragment37( 11, 5, 2,  6.791800*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment37);
+  //  static const G4StableFermiFragment Fragment37( 11, 5,  6,  7.29*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment37);
+  static const G4StableFermiFragment Fragment38( 11, 5,  6,  7.285510*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment38);
+  //  static const G4StableFermiFragment Fragment38( 11, 5,  4,  7.98*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment38);
+  static const G4StableFermiFragment Fragment39( 11, 5,  4,  7.977840*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment39);
+  //  static const G4StableFermiFragment Fragment39( 11, 5,  6,  8.56*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment39); 
+  static const G4StableFermiFragment Fragment40( 11, 5,  6,  8.560300*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment40);   
+  static const G4StableFermiFragment Fragment41( 11, 6,  4,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment41);
+  static const G4StableFermiFragment Fragment42( 11, 6,  2,  2.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment42);
+  //  static const G4StableFermiFragment Fragment42( 11, 6,  6,  4.32*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment42);
+  static const G4StableFermiFragment Fragment43( 11, 6,  6,  4.318800*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment43);
+  //  static const G4StableFermiFragment Fragment43( 11, 6,  4,  4.80*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment43);
+  static const G4StableFermiFragment Fragment44( 11, 6,  4,  4.804200*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment44);
+  //  static const G4StableFermiFragment Fragment44( 11, 6,  2,  6.34*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment44);
+  static const G4StableFermiFragment Fragment45( 11, 6,  2,  6.339200*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment45);
+  //  static const G4StableFermiFragment Fragment45( 11, 6,  8,  6.48*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment45);
+  static const G4StableFermiFragment Fragment46( 11, 6,  8,  6.478200*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment46);
+  //  static const G4StableFermiFragment Fragment46( 11, 6,  6,  6.90*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment46);
+  static const G4StableFermiFragment Fragment47( 11, 6,  6,  6.904800*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment47);
+  //  static const G4StableFermiFragment Fragment47( 11, 6,  4,  7.50*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment47);
+  static const G4StableFermiFragment Fragment48( 11, 6,  4,  7.499700*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment48);
+  //  static const G4StableFermiFragment Fragment48( 11, 6,  4,  8.10*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment48);
+  static const G4StableFermiFragment Fragment49( 11, 6,  4,  8.104500*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment49);
+  //  static const G4StableFermiFragment Fragment49( 11, 6,  6,  8.42*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment49);
+  static const G4StableFermiFragment Fragment50( 11, 6,  6,  8.420000*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment50);
+  static const G4StableFermiFragment Fragment51( 12, 5,  3,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment51);
+  //  static const G4StableFermiFragment Fragment51( 12, 5,  5,  0.95*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment51);
+  static const G4StableFermiFragment Fragment52( 12, 5,  5,  0.9531400*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment52);
+  //  static const G4StableFermiFragment Fragment52( 12, 5,  5,  1.67*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment52);
+  static const G4StableFermiFragment Fragment53( 12, 5,  5,  1.673650*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment53);
+  //  static const G4StableFermiFragment Fragment53( 12, 5,  4,  2.65*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment53);
+  static const G4StableFermiFragment Fragment54( 12, 5,  3,  2.620800*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment54);
+  static const G4StableFermiFragment Fragment55( 12, 6,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment55);
+  //  static const G4StableFermiFragment Fragment55( 12, 6,  5,  4.44*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment55);
+  static const G4StableFermiFragment Fragment56( 12, 6,  5,  4.438910*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment56);
+  static const G4StableFermiFragment Fragment57( 13, 6,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment57);
+  //  static const G4StableFermiFragment Fragment57( 13, 6,  2,  3.09*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment57);
+  static const G4StableFermiFragment Fragment58( 13, 6,  2,  3.089443*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment58);
+  //  static const G4StableFermiFragment Fragment58( 13, 6,  4,  3.68*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment58);
+  static const G4StableFermiFragment Fragment59( 13, 6,  4,  3.684507*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment59);  
+  //  static const G4StableFermiFragment Fragment59( 13, 6,  6,  3.85*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment59);
+  static const G4StableFermiFragment Fragment60( 13, 6,  6,  3.853807*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment60);
+  static const G4StableFermiFragment Fragment61( 13, 7,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment61);
+  static const G4StableFermiFragment Fragment62( 14, 6,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment62);
+  //  static const G4StableFermiFragment Fragment62( 14, 6,  3,  6.09*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment62);
+  static const G4StableFermiFragment Fragment63( 14, 6,  3,  6.093800*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment63);
+  // JMQ 010709 corrected excitation energies for 64-66, according to http://www.nndc.bnl.gov/chart
+  //  static const G4StableFermiFragment Fragment63( 14, 6,  1,  6.59*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment63);
+  static const G4StableFermiFragment Fragment64( 14, 6,  1,  6.589400*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment64);
+  //  static const G4StableFermiFragment Fragment64( 14, 6,  7,  6.73*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment64);
+  static const G4StableFermiFragment Fragment65( 14, 6,  7,  6.728200*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment65);
+  //  static const G4StableFermiFragment Fragment65( 14, 6,  1,  6.90*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment65);
+  static const G4StableFermiFragment Fragment66( 14, 6,  1,  6.902600*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment66);
+  //  static const G4StableFermiFragment Fragment66( 14, 6,  5,  7.01*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment66);
+  static const G4StableFermiFragment Fragment67( 14, 6,  5,  7.012000*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment67);
+  //  static const G4StableFermiFragment Fragment67( 14, 6,  5,  7.34*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment67);
+  static const G4StableFermiFragment Fragment68( 14, 6,  5,  7.341000*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment68);
+  
+  static const G4StableFermiFragment Fragment69( 14, 7,  3,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment69);
+  //  static const G4StableFermiFragment Fragment69( 14, 7,  1,  2.31*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment69);
+  static const G4StableFermiFragment Fragment70( 14, 7,  1,  2.312798*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment70);
+  //  static const G4StableFermiFragment Fragment70( 14, 7,  3,  3.95*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment70);
+  static const G4StableFermiFragment Fragment71( 14, 7,  3,  3.948100*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment71);  
+  //  static const G4StableFermiFragment Fragment71( 14, 7,  1,  4.92*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment71);
+  static const G4StableFermiFragment Fragment72( 14, 7,  1,  4.915100*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment72);
+  //  static const G4StableFermiFragment Fragment72( 14, 7,  5,  5.11*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment72);
+  static const G4StableFermiFragment Fragment73( 14, 7,  5,  5.105890*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment73);
+  //  static const G4StableFermiFragment Fragment73( 14, 7,  3,  5.69*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment73);
+  static const G4StableFermiFragment Fragment74( 14, 7,  3,  5.691440*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment74);
+  //  static const G4StableFermiFragment Fragment74( 14, 7,  7,  5.83*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment74);
+  static const G4StableFermiFragment Fragment75( 14, 7,  7,  5.834250*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment75);
+  //  static const G4StableFermiFragment Fragment75( 14, 7,  3,  6.20*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment75);
+  static const G4StableFermiFragment Fragment76( 14, 7,  3,  6.203500*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment76);
+  //  static const G4StableFermiFragment Fragment76( 14, 7,  7,  6.44*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment76);
+  static const G4StableFermiFragment Fragment77( 14, 7,  7,  6.446170*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment77);
+  //  static const G4StableFermiFragment Fragment77( 14, 7,  5,  7.03*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment77);
+  static const G4StableFermiFragment Fragment78( 14, 7,  5,  7.029120*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment78);
+  static const G4StableFermiFragment Fragment79( 15, 7,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment79);
+  // JMQ 010709 two very close levels instead of only one, with their own spins
+  //  static const G4StableFermiFragment Fragment79( 15, 7,  6,  5.27*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment79);
+  static const G4StableFermiFragment Fragment80( 15, 7,  6,  5.270155*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment80);
+  //  static const G4StableFermiFragment Fragment80( 15, 7,  2,  5.30*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment80);
+  static const G4StableFermiFragment Fragment81( 15, 7,  2,  5.298822*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment81);
+  //  static const G4StableFermiFragment Fragment81( 15, 7,  4,  6.32*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment81);
+  static const G4StableFermiFragment Fragment82( 15, 7,  4,  6.323780*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment82);
+  //JMQ 010709 new level and corrected energy and spins
+  //  static const G4StableFermiFragment Fragment82( 15, 7,  6,  7.15*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment82);
+  static const G4StableFermiFragment Fragment83( 15, 7,  6,  7.155050*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment83);
+  //  static const G4StableFermiFragment Fragment83( 15, 7,  4,  7.30*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment83);
+  static const G4StableFermiFragment Fragment84( 15, 7,  4,  7.300830*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment84);
+  //  static const G4StableFermiFragment Fragment84( 15, 7,  8,  7.57*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment84);
+  static const G4StableFermiFragment Fragment85( 15, 7,  8,  7.567100*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment85);
+  //  static const G4StableFermiFragment Fragment85( 15, 7,  2,  8.31*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment85);
+  static const G4StableFermiFragment Fragment86( 15, 7,  2,  8.312620*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment86);
+  //  static const G4StableFermiFragment Fragment86( 15, 7,  4,  8.57*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment86);
+  static const G4StableFermiFragment Fragment87( 15, 7,  4,  8.571400*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment87);
+  //  static const G4StableFermiFragment Fragment87( 15, 7,  2,  9.05*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment87);
+  static const G4StableFermiFragment Fragment88( 15, 7,  2,  9.049710*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment88);
+  //JMQ 010709 new levels for N15
+  //  static const G4StableFermiFragment Fragment88( 15, 7,  4,  9.151*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment88);
+  static const G4StableFermiFragment Fragment89( 15, 7,  4,  9.151900*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment89);
+  //  static const G4StableFermiFragment Fragment89( 15, 7,  6,  9.154*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment89);
+  static const G4StableFermiFragment Fragment90( 15, 7,  6,  9.154900*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment90);
+  //  static const G4StableFermiFragment Fragment90( 15, 7,  2,  9.22*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment90);
+  static const G4StableFermiFragment Fragment91( 15, 7,  2,  9.222100*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment91);
+  //  static const G4StableFermiFragment Fragment91( 15, 7,  6,  9.76*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment91);
+  static const G4StableFermiFragment Fragment92( 15, 7,  6,  9.760000*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment92);
+  //  static const G4StableFermiFragment Fragment92( 15, 7,  8,  9.83*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment92);
+  static const G4StableFermiFragment Fragment93( 15, 7,  8,  9.829000*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment93);
+  //  static const G4StableFermiFragment Fragment93( 15, 7,  4,  9.93*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment93);
+  static const G4StableFermiFragment Fragment94( 15, 7,  4,  9.925000*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment94);
+  //  static const G4StableFermiFragment Fragment94( 15, 7,  4, 10.07*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment94);
+  static const G4StableFermiFragment Fragment95( 15, 7,  4, 10.06600*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment95);
+  
+  static const G4StableFermiFragment Fragment96( 15, 8,  2,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment96);
+  //JMQ 010709 new level and spins
+  //  static const G4StableFermiFragment Fragment96( 15, 8,  2,  5.18*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment96);
+  static const G4StableFermiFragment Fragment97( 15, 8,  2,  5.183000*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment97);
+  //  static const G4StableFermiFragment Fragment97( 15, 8,  6,  5.24*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment97);
+  static const G4StableFermiFragment Fragment98( 15, 8,  6,  5.240900*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment98);
+  //  static const G4StableFermiFragment Fragment98( 15, 8,  4,  6.18*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment98);  
+  static const G4StableFermiFragment Fragment99( 15, 8,  4,  6.176300*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment99);  
+  //  static const G4StableFermiFragment Fragment99( 15, 8,  4,  6.79*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment99);
+  static const G4StableFermiFragment Fragment100( 15, 8,  4,  6.793100*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment100);
+  //  static const G4StableFermiFragment Fragment100( 15, 8,  6,  6.86*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment100);
+  static const G4StableFermiFragment Fragment101( 15, 8,  6,  6.859400*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment101);
+  //  static const G4StableFermiFragment Fragment101( 15, 8,  8,  7.28*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment101);
+  static const G4StableFermiFragment Fragment102( 15, 8,  8,  7.275900*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment102);
+  
+  
+  static const G4StableFermiFragment Fragment103( 16, 7,  5,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment103);
+  //  static const G4StableFermiFragment Fragment103( 16, 7,  1,  0.12*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment103);
+  static const G4StableFermiFragment Fragment104( 16, 7,  1,  0.1204200*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment104);
+  //  static const G4StableFermiFragment Fragment104( 16, 7,  7,  0.30*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment104);
+  static const G4StableFermiFragment Fragment105( 16, 7,  7,  0.2982200*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment105);
+  //  static const G4StableFermiFragment Fragment105( 16, 7,  3,  0.40*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment105);
+  static const G4StableFermiFragment Fragment106( 16, 7,  3,  0.3972700*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment106);
+  
+  //JMQ 010709   some energies and spins have been changed 
+  static const G4StableFermiFragment Fragment107( 16, 8,  1,  0.00*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment107);
+  //  static const G4StableFermiFragment Fragment107( 16, 8,  1,  6.05*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment107);
+  static const G4StableFermiFragment Fragment108( 16, 8,  1,  6.049400*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment108);
+  //  static const G4StableFermiFragment Fragment108( 16, 8,  7,  6.13*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment108);
+  static const G4StableFermiFragment Fragment109( 16, 8,  7,  6.129890*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment109);
+  //  static const G4StableFermiFragment Fragment109( 16, 8,  5,  6.92*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment109);
+  static const G4StableFermiFragment Fragment110( 16, 8,  5,  6.917100*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment110);
+  //  static const G4StableFermiFragment Fragment110( 16, 8,  3,  7.12*MeV ); if(MapIsEmpty) fragment_pool.push_back
+  //JMQ 180510 fixed fragment 111
+  static const G4StableFermiFragment Fragment111( 16, 8,  3,  7.116850*MeV ); if(MapIsEmpty) fragment_pool.push_back(&Fragment111);
 
   static std::multimap<const std::pair<G4int,G4int>, const G4VFermiFragment* , std::less<const std::pair<G4int,G4int> > >  
     theMapOfFragments;
-
+  
   if (MapIsEmpty) 
     {
       for(size_t i=0; i<fragment_pool.size(); i++)
-      {
-        theMapOfFragments.insert(std::pair<const std::pair<G4int,G4int>, 
-				 const G4VFermiFragment* >(std::pair<G4int,G4int>(fragment_pool[i]->GetA(),
+	{
+	  theMapOfFragments.insert(std::pair<const std::pair<G4int,G4int>, 
+				   const G4VFermiFragment* >(std::pair<G4int,G4int>(fragment_pool[i]->GetA(),
 										  fragment_pool[i]->GetZ()),fragment_pool[i]));
-      }
+	}
       MapIsEmpty = false;
     }
Index: trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiPhaseSpaceDecay.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiPhaseSpaceDecay.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiPhaseSpaceDecay.cc	(revision 1315)
@@ -37,4 +37,11 @@
 #include <functional>
 
+G4FermiPhaseSpaceDecay::G4FermiPhaseSpaceDecay()
+{
+}
+
+G4FermiPhaseSpaceDecay::~G4FermiPhaseSpaceDecay()
+{
+}
 
 std::vector<G4LorentzVector*> *
Index: trunk/source/processes/hadronic/models/de_excitation/fission/include/G4FissionBarrier.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/fission/include/G4FissionBarrier.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/fission/include/G4FissionBarrier.hh	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4FissionBarrier.hh,v 1.3 2006/06/29 20:13:21 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4FissionBarrier.hh,v 1.5 2009/12/16 16:50:07 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4CompetitiveFission.cc,v 1.11 2009/03/13 18:57:17 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4CompetitiveFission.cc,v 1.12 2010/05/03 16:49:19 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -38,4 +38,5 @@
 #include "G4CompetitiveFission.hh"
 #include "G4PairingCorrection.hh"
+#include "G4ParticleMomentum.hh"
 
 G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission")
Index: trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionBarrier.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionBarrier.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionBarrier.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4FissionBarrier.cc,v 1.5 2006/06/29 20:13:35 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4FissionBarrier.cc,v 1.7 2009/12/16 16:50:07 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4EvaporationGEMFactory.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4EvaporationGEMFactory.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4EvaporationGEMFactory.hh	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4EvaporationGEMFactory.hh,v 1.7 2009/09/15 12:54:16 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4EvaporationGEMFactory.hh,v 1.8 2010/04/27 11:43:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -41,9 +41,9 @@
 {
 public:
-  G4EvaporationGEMFactory() {}
-  virtual ~G4EvaporationGEMFactory() {}
+  G4EvaporationGEMFactory();
+  virtual ~G4EvaporationGEMFactory(); 
 
 private:
-  G4EvaporationGEMFactory(const G4EvaporationGEMFactory & ) : G4VEvaporationFactory() {}
+  G4EvaporationGEMFactory(const G4EvaporationGEMFactory & );
   const G4EvaporationGEMFactory & operator=(const G4EvaporationGEMFactory & val);
   G4bool operator==(const G4EvaporationGEMFactory & val) const;
Index: trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4GEMProbability.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4GEMProbability.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4GEMProbability.hh	(revision 1315)
@@ -24,4 +24,10 @@
 // ********************************************************************
 //
+// $Id: G4GEMProbability.hh,v 1.4 2010/05/19 10:21:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
+//
+//---------------------------------------------------------------------
+//
+// Geant4 header G4GEMProbability
 //
 //
@@ -29,6 +35,8 @@
 // by V. Lara (Sept 2001) 
 //
-
-
+// 18.05.2010 V.Ivanchenko trying to speedup the most slow method
+//            by usage of G4Pow, integer Z and A; moved constructor, 
+//            destructor and virtual functions to source
+//
 
 #ifndef G4GEMProbability_h
@@ -42,133 +50,150 @@
 #include "G4PairingCorrection.hh"
 
+class G4Pow;
+
 class G4GEMProbability : public G4VEmissionProbability
 {
 public:
-    // Only available constructor
-    G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin) : 
-        theA(anA), theZ(aZ), Spin(aSpin), theCoulombBarrierPtr(0),
-        ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
-        {
-            theEvapLDPptr = new G4EvaporationLevelDensityParameter;
-        }
+
+  // Default constructor - should not be used
+  G4GEMProbability();
+
+  // Only available constructor
+  G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin);
     
-    ~G4GEMProbability()
-        {
-            if (theEvapLDPptr != 0) delete theEvapLDPptr;
-        }
+  virtual ~G4GEMProbability();
 
+  inline G4int GetZ_asInt(void) const { return theZ; }
+	
+  inline G4int GetA_asInt(void) const { return theA;}
+	
+  inline G4double GetZ(void) const { return theZ; }
+	
+  inline G4double GetA(void) const { return theA;}
 
-	
-    G4double GetZ(void) const { return theZ; }
-	
-    G4double GetA(void) const { return theA;}
+  inline G4double GetSpin(void) const { return Spin; }
 
-    G4double GetSpin(void) const { return Spin; }
+  inline G4double GetNormalization(void) const { return Normalization; }
+    
+  inline void SetCoulomBarrier(const G4VCoulombBarrier * aCoulombBarrierStrategy)
+  {
+    theCoulombBarrierPtr = aCoulombBarrierStrategy;
+  }
 
-    G4double GetNormalization(void) const { return Normalization; }
+  inline G4double GetCoulombBarrier(const G4Fragment& fragment) const 
+  {
+    G4double res = 0.0;
+    if (theCoulombBarrierPtr) 
+      {
+	G4int Acomp = fragment.GetA_asInt();
+	G4int Zcomp = fragment.GetZ_asInt();
+	res = theCoulombBarrierPtr->GetCoulombBarrier(Acomp-theA, Zcomp-theZ,
+	  fragment.GetExcitationEnergy()-fPairCorr->GetPairingCorrection(Acomp,Zcomp));
+      }
+    return res;
+  }
     
-    void SetCoulomBarrier(const G4VCoulombBarrier * aCoulombBarrierStrategy)
-        {
-            theCoulombBarrierPtr = aCoulombBarrierStrategy;
-        }
-
-    G4double GetCoulombBarrier(const G4Fragment& fragment) const 
-        {
-            if (theCoulombBarrierPtr) 
-	      {
-		G4int Acompound = static_cast<G4int>(fragment.GetA());
-		G4int Zcompound = static_cast<G4int>(fragment.GetZ());
-		return theCoulombBarrierPtr->GetCoulombBarrier(Acompound-theA, Zcompound-theZ,
-							       fragment.GetExcitationEnergy()-
-							       G4PairingCorrection::GetInstance()->
-							       GetPairingCorrection(Acompound,Zcompound));
-	      }
-            else 
-	      {
-		return 0.0;
-	      }
-        }
-    
-
-    virtual G4double CalcAlphaParam(const G4Fragment & ) const {return 1.0;}
-    virtual G4double CalcBetaParam(const G4Fragment & ) const {return 1.0;}
+  virtual G4double CalcAlphaParam(const G4Fragment & ) const;
+  virtual G4double CalcBetaParam(const G4Fragment & ) const;
     
 protected:
   
-    void SetExcitationEnergiesPtr(std::vector<G4double> * anExcitationEnergiesPtr) 
-        {
-            ExcitationEnergies = anExcitationEnergiesPtr;
-        }
+  inline void SetExcitationEnergiesPtr(std::vector<G4double> * anExcitationEnergiesPtr) 
+  {
+    ExcitationEnergies = anExcitationEnergiesPtr;
+  }
   
-    void SetExcitationSpinsPtr(std::vector<G4double> * anExcitationSpinsPtr)
-        {
-            ExcitationSpins = anExcitationSpinsPtr;
-        }
+  inline void SetExcitationSpinsPtr(std::vector<G4double> * anExcitationSpinsPtr)
+  {
+    ExcitationSpins = anExcitationSpinsPtr;
+  }
 
-    void SetExcitationLifetimesPtr(std::vector<G4double> * anExcitationLifetimesPtr)
-        {
-            ExcitationLifetimes = anExcitationLifetimesPtr;
-        }
+  inline void SetExcitationLifetimesPtr(std::vector<G4double> * anExcitationLifetimesPtr)
+  {
+    ExcitationLifetimes = anExcitationLifetimesPtr;
+  }
 
-    void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier)
-        {
-            theCoulombBarrierPtr = aCoulombBarrier;
-        }
+  inline void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier)
+  {
+    theCoulombBarrierPtr = aCoulombBarrier;
+  }
 
-  
-    // Default constructor
-    G4GEMProbability() {}
 private:
-    // Copy constructor
-    G4GEMProbability(const G4GEMProbability &right);
+
+  // Copy constructor
+  G4GEMProbability(const G4GEMProbability &right);
     
-    const G4GEMProbability & operator=(const G4GEMProbability &right);
-    G4bool operator==(const G4GEMProbability &right) const;
-    G4bool operator!=(const G4GEMProbability &right) const;
+  const G4GEMProbability & operator=(const G4GEMProbability &right);
+  G4bool operator==(const G4GEMProbability &right) const;
+  G4bool operator!=(const G4GEMProbability &right) const;
     
 public:
-    G4double EmissionProbability(const G4Fragment & fragment, const G4double anEnergy);
+
+  G4double EmissionProbability(const G4Fragment & fragment, const G4double anEnergy);
   
 private:
 
-    G4double CalcProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy,
-                             const G4double V);
-    virtual G4double CCoeficient(const G4double ) const {return 0.0;};
+  G4double CalcProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy,
+			   const G4double V);
 
+  virtual G4double CCoeficient(const G4double ) const;
+
+  inline G4double I0(const G4double t);
+  inline G4double I1(const G4double t, const G4double tx);
+  inline G4double I2(const G4double s, const G4double sx);
+  G4double I3(const G4double s, const G4double sx);
     
-    G4double I0(const G4double t);
-    G4double I1(const G4double t, const G4double tx);
-    G4double I2(const G4double s, const G4double sx);
-    G4double I3(const G4double s, const G4double sx);
+  // Data Members
+
+  G4Pow*   fG4pow;
+  G4PairingCorrection* fPairCorr;
     
-    // Data Members
+  G4VLevelDensityParameter * theEvapLDPptr;
+	
+  G4int theA;
+  G4int theZ;
     
-    G4VLevelDensityParameter * theEvapLDPptr;
-	
-    G4int theA;
-    G4int theZ;
+  // Spin is fragment spin
+  G4double Spin;
+
+  // Coulomb Barrier
+  const G4VCoulombBarrier * theCoulombBarrierPtr;
+  
+  // Resonances Energy
+  std::vector<G4double> * ExcitationEnergies;
     
-    // Spin is fragment spin
-    G4double Spin;
+  // Resonances Spin 
+  std::vector<G4double> * ExcitationSpins;
 
-    // Coulomb Barrier
-    const G4VCoulombBarrier * theCoulombBarrierPtr;
+  // Resonances half lifetime
+  std::vector<G4double> * ExcitationLifetimes;
 
-    
-    // Resonances Energy
-    std::vector<G4double> * ExcitationEnergies;
-    
-    // Resonances Spin 
-    std::vector<G4double> * ExcitationSpins;
-
-    // Resonances half lifetime
-    std::vector<G4double> * ExcitationLifetimes;
-
-
-    // Normalization
-    G4double Normalization;
+  // Normalization
+  G4double Normalization;
     
 };
 
+inline G4double G4GEMProbability::I0(const G4double t)
+{
+  return std::exp(t) - 1.0;
+}
+
+inline G4double G4GEMProbability::I1(const G4double t, const G4double tx)
+{
+  return (t - tx + 1.0)*std::exp(tx) - t - 1.0;
+}
+
+
+inline G4double G4GEMProbability::I2(const G4double s, const G4double sx)
+{
+  G4double S = 1.0/std::sqrt(s);
+  G4double Sx = 1.0/std::sqrt(sx);
+  
+  G4double p1 = S*S*S*( 1.0 + S*S*( 1.5 + 3.75*S*S) );
+  G4double p2 = Sx*Sx*Sx*( 1.0 + Sx*Sx*( 1.5 + 3.75*Sx*Sx) )*std::exp(sx-s);
+  
+  return p1-p2;
+}
+
 
 #endif
Index: trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4EvaporationGEMFactory.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4EvaporationGEMFactory.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4EvaporationGEMFactory.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4EvaporationGEMFactory.cc,v 1.9 2009/09/15 12:54:16 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4EvaporationGEMFactory.cc,v 1.10 2010/04/27 11:43:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -103,24 +103,9 @@
 #include "G4PhotonEvaporation.hh"
 
-const G4EvaporationGEMFactory & 
-G4EvaporationGEMFactory::operator=(const G4EvaporationGEMFactory & )
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationGEMFactory::operator= meant to not be accessable.");
-  return *this;
-}
-
-G4bool 
-G4EvaporationGEMFactory::operator==(const G4EvaporationGEMFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationGEMFactory::operator== meant to not be accessable.");
-  return false;
-}
-
-G4bool 
-G4EvaporationGEMFactory::operator!=(const G4EvaporationGEMFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationGEMFactory::operator!= meant to not be accessable.");
-  return true;
-}
+G4EvaporationGEMFactory::G4EvaporationGEMFactory()
+{} 
+  
+G4EvaporationGEMFactory::~G4EvaporationGEMFactory() 
+{}
                  
 std::vector<G4VEvaporationChannel*> * G4EvaporationGEMFactory::CreateChannel()
@@ -129,4 +114,7 @@
     new std::vector<G4VEvaporationChannel*>;
   theChannel->reserve(68);
+
+  theChannel->push_back( new G4PhotonEvaporation() );  // Photon Channel
+  theChannel->push_back( new G4CompetitiveFission() ); // Fission Channel
 
   theChannel->push_back( new G4NeutronGEMChannel() );  // n
@@ -197,7 +185,4 @@
   theChannel->push_back( new G4Mg28GEMChannel() );     // Mg28
 
-  theChannel->push_back( new G4CompetitiveFission() ); // Fission Channel
-  theChannel->push_back( new G4PhotonEvaporation() );  // Photon Channel
-
   return theChannel;
 }
Index: trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMProbability.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMProbability.cc	(revision 1315)
@@ -24,4 +24,14 @@
 // ********************************************************************
 //
+// $Id: G4GEMProbability.cc,v 1.13 2010/05/19 10:21:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
+//
+//---------------------------------------------------------------------
+//
+// Geant4 class G4GEMProbability
+//
+//
+// Hadronic Process: Nuclear De-excitations
+// by V. Lara (Sept 2001) 
 //
 //
@@ -35,34 +45,51 @@
 //       -InitialLeveldensity calculated according to the right conditions of the 
 //        initial excited nucleus.
+// J. M. Quesada 19/04/2010 fix in  emission probability calculation.
+// V.Ivanchenko  20/04/2010 added usage of G4Pow and use more safe computation
+// V.Ivanchenko 18/05/2010 trying to speedup the most slow method
+//            by usage of G4Pow, integer Z and A; moved constructor, 
+//            destructor and virtual functions to source
 
 #include "G4GEMProbability.hh"
 #include "G4PairingCorrection.hh"
-
-
-
-G4GEMProbability::G4GEMProbability(const G4GEMProbability &) : G4VEmissionProbability()
-{
-    throw G4HadronicException(__FILE__, __LINE__, "G4GEMProbability::copy_constructor meant to not be accessable");
-}
-
-
-
-
-const G4GEMProbability & G4GEMProbability::
-operator=(const G4GEMProbability &)
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4GEMProbability::operator= meant to not be accessable");
-  return *this;
-}
-
-
-G4bool G4GEMProbability::operator==(const G4GEMProbability &) const
-{
-  return false;
-}
-
-G4bool G4GEMProbability::operator!=(const G4GEMProbability &) const
-{
-  return true;
+#include "G4Pow.hh"
+#include "G4IonTable.hh"
+
+G4GEMProbability:: G4GEMProbability() : 
+  theA(0), theZ(0), Spin(0.0), theCoulombBarrierPtr(0),
+  ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
+{
+  theEvapLDPptr = new G4EvaporationLevelDensityParameter;
+  fG4pow = G4Pow::GetInstance(); 
+  fPairCorr = G4PairingCorrection::GetInstance();
+}
+
+G4GEMProbability:: G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin) : 
+  theA(anA), theZ(aZ), Spin(aSpin), theCoulombBarrierPtr(0),
+  ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
+{
+  theEvapLDPptr = new G4EvaporationLevelDensityParameter;
+  fG4pow = G4Pow::GetInstance(); 
+  fPairCorr = G4PairingCorrection::GetInstance();
+}
+    
+G4GEMProbability::~G4GEMProbability()
+{
+  delete theEvapLDPptr;
+}
+
+G4double G4GEMProbability::CalcAlphaParam(const G4Fragment & ) const 
+{
+  return 1.0;
+}
+
+G4double G4GEMProbability::CalcBetaParam(const G4Fragment & ) const 
+{
+  return 1.0;
+}
+
+G4double G4GEMProbability::CCoeficient(const G4double ) const 
+{
+  return 0.0;
 }
 
@@ -80,5 +107,5 @@
     if (ExcitationEnergies  &&  ExcitationSpins && ExcitationLifetimes) {
       G4double SavedSpin = Spin;
-      for (unsigned int i = 0; i < ExcitationEnergies->size(); i++) {
+      for (size_t i = 0; i < ExcitationEnergies->size(); ++i) {
 	Spin = ExcitationSpins->operator[](i);
 	// substract excitation energies
@@ -86,6 +113,8 @@
 	if (Tmax > 0.0) {
 	  G4double width = CalcProbability(fragment,Tmax,CoulombBarrier);
+	  //JMQ April 2010 added condition to prevent reported crash
 	  // update probability
-	  if (hbar_Planck*std::log(2.0)/width < ExcitationLifetimes->operator[](i)) {
+	  if (width > 0. && 
+	      hbar_Planck*fG4pow->logZ(2) < width*ExcitationLifetimes->operator[](i)) {
 	    EmissionProbability += width;
 	  }
@@ -98,4 +127,139 @@
   return EmissionProbability;
 }
+
+
+G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment, 
+                                           const G4double MaximalKineticEnergy,
+                                           const G4double V)
+
+// Calculate integrated probability (width) for evaporation channel
+{
+  G4int A = fragment.GetA_asInt();
+  G4int Z = fragment.GetZ_asInt();
+
+  G4int ResidualA = A - theA;
+  G4int ResidualZ = Z - theZ;
+  G4double U = fragment.GetExcitationEnergy();
+  
+  G4double NuclearMass = fragment.ComputeGroundStateMass(theZ, theA);
+  
+  G4double Alpha = CalcAlphaParam(fragment);
+  G4double Beta = CalcBetaParam(fragment);
+  
+  
+  //                       ***RESIDUAL***
+  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
+  
+  G4double delta0 = fPairCorr->GetPairingCorrection(ResidualA, ResidualZ);  
+  
+  G4double a = theEvapLDPptr->LevelDensityParameter(ResidualA,
+  						    ResidualZ,MaximalKineticEnergy+V-delta0);
+  G4double Ux = (2.5 + 150.0/G4double(ResidualA))*MeV;
+  G4double Ex = Ux + delta0;
+  G4double T  = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
+  //JMQ fixed bug in units
+  G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
+  //                      ***end RESIDUAL ***
+  
+  //                       ***PARENT***
+  //JMQ (September 2009) the following quantities refer to the PARENT:
+     
+  G4double deltaCN=fPairCorr->GetPairingCorrection(A, Z);				       
+  G4double aCN = theEvapLDPptr->LevelDensityParameter(A, Z, U-deltaCN);
+  G4double UxCN = (2.5 + 150.0/G4double(A))*MeV;
+  G4double ExCN = UxCN + deltaCN;
+  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
+  //JMQ fixed bug in units
+  G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
+//                       ***end PARENT***
+
+  G4double Width;
+  G4double InitialLevelDensity;
+  G4double t = MaximalKineticEnergy/T;
+  if ( MaximalKineticEnergy < Ex ) {
+    //JMQ 190709 bug in I1 fixed (T was  missing)
+    Width = (I1(t,t)*T + (Beta+V)*I0(t))/std::exp(E0/T);
+    //JMQ 160909 fix:  InitialLevelDensity has been taken away 
+    //(different conditions for initial CN..) 
+  } else {
+
+    //VI minor speedup
+    G4double expE0T = std::exp(E0/T);
+    const G4double sqrt2 = std::sqrt(2.0);
+
+    G4double tx = Ex/T;
+    G4double s = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
+    G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
+    Width = I1(t,tx)*T/expE0T + I3(s,sx)*std::exp(s)/(sqrt2*a);
+    // For charged particles (Beta+V) = 0 beacuse Beta = -V
+    if (theZ == 0) {
+      Width += (Beta+V)*(I0(tx)/expE0T + 2.0*sqrt2*I2(s,sx)*std::exp(s));
+    }
+  }
+  
+  //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of hbar_planck must be used
+  //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
+  G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
+  
+  //JMQ 190709 fix on Rb and  geometrical cross sections according to Furihata's paper 
+  //                      (JAERI-Data/Code 2001-105, p6)
+  //    G4double RN = 0.0;
+  G4double Rb = 0.0;
+  if (theA > 4) 
+    {
+      //	G4double R1 = std::pow(ResidualA,1.0/3.0);
+      //	G4double R2 = std::pow(G4double(theA),1.0/3.0);
+      G4double Ad = fG4pow->Z13(ResidualA);
+      G4double Aj = fG4pow->Z13(theA);
+      //	RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
+      Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
+      Rb *= fermi;
+    }
+  else if (theA>1)
+    {
+      G4double Ad = fG4pow->Z13(ResidualA);
+      G4double Aj = fG4pow->Z13(theA);
+      Rb=1.5*(Aj+Ad)*fermi;
+    }
+  else
+    {
+      G4double Ad = fG4pow->Z13(ResidualA);
+      Rb = 1.5*Ad*fermi;
+    }
+  //   G4double GeometricalXS = pi*RN*RN*std::pow(ResidualA,2./3.); 
+  G4double GeometricalXS = pi*Rb*Rb; 
+  //end of JMQ fix on Rb by 190709
+  
+
+
+//JMQ 160909 fix: initial level density must be calculated according to the 
+// conditions at the initial compound nucleus 
+// (it has been removed from previous "if" for the residual) 
+
+   if ( U < ExCN ) 
+      {
+        InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
+      } 
+    else 
+      {
+	//        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
+	//VI speedup
+        G4double x  = U-deltaCN;
+        G4double x2 = x*x;
+        G4double x3 = x2*x;
+        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*x))/std::sqrt(std::sqrt(aCN*x2*x3));
+      }
+ //
+
+
+  //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
+  //    Width *= std::sqrt(pi)*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity); 
+  Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity); 
+   
+
+  return Width;
+}
+
+/*
 
 G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment, 
@@ -219,33 +383,5 @@
   return Width;
 }
-
-
-
-
-G4double G4GEMProbability::I0(const G4double t)
-{
-  G4double result = (std::exp(t) - 1.0);
-  return result;
-}
-
-G4double G4GEMProbability::I1(const G4double t, const G4double tx)
-{
-  G4double result = t - tx + 1.0;
-  result *= std::exp(tx);
-  result -= (t + 1.0);
-  return result;
-}
-
-
-G4double G4GEMProbability::I2(const G4double s, const G4double sx)
-{
-  G4double S = 1.0/std::sqrt(s);
-  G4double Sx = 1.0/std::sqrt(sx);
-  
-  G4double p1 = S*S*S*( 1.0 + S*S*( 1.5 + 3.75*S*S) );
-  G4double p2 = Sx*Sx*Sx*( 1.0 + Sx*Sx*( 1.5 + 3.75*Sx*Sx) )*std::exp(sx-s);
-  
-  return p1-p2;
-}
+*/
 
 G4double G4GEMProbability::I3(const G4double s, const G4double sx)
Index: trunk/source/processes/hadronic/models/de_excitation/handler/include/G4ExcitationHandler.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/handler/include/G4ExcitationHandler.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/handler/include/G4ExcitationHandler.hh	(revision 1315)
@@ -24,7 +24,6 @@
 // ********************************************************************
 //
-//
-// $Id: G4ExcitationHandler.hh,v 1.10 2009/11/24 11:18:48 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4ExcitationHandler.hh,v 1.12 2010/04/27 14:00:23 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -55,4 +54,5 @@
 #include "G4VEvaporation.hh"
 #include "G4VPhotonEvaporation.hh"
+#include "G4VEvaporationChannel.hh"
 #include "G4Fragment.hh"
 #include "G4DynamicParticle.hh"
@@ -68,6 +68,4 @@
 #include "G4IonConstructor.hh" 
 
-//#define debug
-
 class G4ExcitationHandler 
 {
@@ -75,13 +73,14 @@
   G4ExcitationHandler(); 
   ~G4ExcitationHandler();
+
 private:
+
   G4ExcitationHandler(const G4ExcitationHandler &right);
-
   const G4ExcitationHandler & operator=(const G4ExcitationHandler &right);
   G4bool operator==(const G4ExcitationHandler &right) const;
   G4bool operator!=(const G4ExcitationHandler &right) const;
   
-
 public:
+
   G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState) const;
 
@@ -92,5 +91,5 @@
   void SetFermiModel(G4VFermiBreakUp *const  value);
 
-  void SetPhotonEvaporation(G4VPhotonEvaporation * const value);
+  void SetPhotonEvaporation(G4VEvaporationChannel * const value);
 
   void SetMaxZForFermiBreakUp(G4int aZ);
@@ -99,11 +98,13 @@
   void SetMinEForMultiFrag(G4double anE);
 
-// for inverse cross section choice
-  inline void SetOPTxs(G4int opt) { OPTxs = opt;}
-// for superimposed Coulomb Barrir for inverse cross sections
-  inline void UseSICB(){useSICB=true;}
+  // for inverse cross section choice
+  inline void SetOPTxs(G4int opt);
+  // for superimposed Coulomb Barrir for inverse cross sections
+  inline void UseSICB();
 
 private:
 
+  void SetParameters();
+
   G4ReactionProductVector * Transform(G4FragmentVector * theFragmentVector) const;
 
@@ -114,5 +115,5 @@
   const G4VFermiBreakUp * GetFermiModel() const;
 
-  const G4VPhotonEvaporation * GetPhotonEvaporation() const;
+  const G4VEvaporationChannel * GetPhotonEvaporation() const;
 
   G4int GetMaxZ() const;
@@ -124,4 +125,5 @@
 			 G4FragmentVector * Result) const;
 #endif
+
 private:
   
@@ -132,9 +134,10 @@
   G4VFermiBreakUp *theFermiModel;
  
-  G4VPhotonEvaporation * thePhotonEvaporation;
+  G4VEvaporationChannel * thePhotonEvaporation;
 
   G4int maxZForFermiBreakUp;
   G4int maxAForFermiBreakUp;
   G4double minEForMultiFrag;
+  G4double minExcitation;
 
   G4ParticleTable *theTableOfParticles;
@@ -147,5 +150,5 @@
   G4int OPTxs;
   G4bool useSICB;
-
+  
   struct DeleteFragment 
   {
@@ -156,9 +159,18 @@
     }
   };
-
-
+  
 };
 
-
+inline void G4ExcitationHandler::SetOPTxs(G4int opt) 
+{ 
+  OPTxs = opt; 
+  SetParameters();
+}
+
+inline void G4ExcitationHandler::UseSICB()
+{ 
+  useSICB = true; 
+  SetParameters();
+}
 
 inline const G4VEvaporation * G4ExcitationHandler::GetEvaporation() const
@@ -172,4 +184,5 @@
   MyOwnEvaporationClass = false;
   theEvaporation = value;
+  SetParameters();
 }
 
@@ -199,10 +212,10 @@
 
 
-inline const G4VPhotonEvaporation * G4ExcitationHandler::GetPhotonEvaporation() const
+inline const G4VEvaporationChannel * G4ExcitationHandler::GetPhotonEvaporation() const
 {
   return thePhotonEvaporation;
 }
 
-inline void G4ExcitationHandler::SetPhotonEvaporation(G4VPhotonEvaporation *const  value)
+inline void G4ExcitationHandler::SetPhotonEvaporation(G4VEvaporationChannel *const  value)
 {
   if (thePhotonEvaporation != 0 && MyOwnPhotonEvaporationClass) delete thePhotonEvaporation;
Index: trunk/source/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4ExcitationHandler.cc,v 1.39 2010/05/18 15:46:11 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -29,16 +31,6 @@
 //
 //
-// Modif (September 2009) by J. M. Quesada: 
-// according to Igor Pshenichnov, SMM will be applied (just in case) only once .
-//
-// Modif (September 2008) by J. M. Quesada. External choices have been added for :
-//                   -inverse cross section option (default OPTxs=3)
-//                   -superimposed Coulomb barrier (if useSICB is set true, by default it is false) 
-//
-// Modif (24 Jul 2008) by M. A. Cortes Giraldo:
-//      -Max Z,A for Fermi Break-Up turns to 9,17 by default
-//      -BreakItUp() reorganised and bug in Evaporation loop fixed
-//      -Transform() optimised
-// Modif (30 June 1998) by V. Lara:
+// Modified:
+// (30 June 1998) by V. Lara:
 //      -Modified the Transform method for use G4ParticleTable and 
 //       therefore G4IonTable. It makes possible to convert all kind 
@@ -49,4 +41,19 @@
 //              MultiFragmentation: G4StatMF 
 //              Fermi Breakup model: G4FermiBreakUp
+// (24 Jul 2008) by M. A. Cortes Giraldo:
+//      -Max Z,A for Fermi Break-Up turns to 9,17 by default
+//      -BreakItUp() reorganised and bug in Evaporation loop fixed
+//      -Transform() optimised
+// (September 2008) by J. M. Quesada. External choices have been added for :
+//      -inverse cross section option (default OPTxs=3)
+//      -superimposed Coulomb barrier (if useSICB is set true, by default it is false) 
+// (September 2009) by J. M. Quesada: 
+//      -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
+// (27 Nov 2009) by V.Ivanchenko: 
+//      -cleanup the logic, reduce number internal vectors, fixed memory leak.
+// (11 May 2010) by V.Ivanchenko: 
+//      -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for 
+//       final photon deexcitation; used check on adundance of a fragment, decay 
+//       unstable fragments with A <5
 //
 
@@ -54,15 +61,15 @@
 #include "globals.hh"
 #include "G4LorentzVector.hh"
+#include "G4NistManager.hh"
 #include <list>
-
-//#define debugphoton
-
 
 G4ExcitationHandler::G4ExcitationHandler():
   // JMQ 160909 Fermi BreakUp & MultiFrag are on by default 
   // This is needed for activation of such models when G4BinaryLightIonReaction is used
-  //  since no interface (for external activation via macro input file) is still available in this case.
-  //maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(3.0*MeV),
-  maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV),
+  // since no interface (for external activation via macro input file) is still available.
+  maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4.0*GeV),
+  //  maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(3.0*MeV),
+  //maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV),
+  minExcitation(CLHEP::keV),
   MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true),
   MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false)
@@ -74,11 +81,6 @@
   theFermiModel = new G4FermiBreakUp;
   thePhotonEvaporation = new G4PhotonEvaporation;
+  SetParameters();
 }
-
-G4ExcitationHandler::G4ExcitationHandler(const G4ExcitationHandler &)
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::copy_constructor: is meant to not be accessable! ");
-}
-
 
 G4ExcitationHandler::~G4ExcitationHandler()
@@ -90,81 +92,77 @@
 }
 
-
-const G4ExcitationHandler & G4ExcitationHandler::operator=(const G4ExcitationHandler &)
+void G4ExcitationHandler::SetParameters()
 {
-  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator=: is meant to not be accessable! ");
-  
-  return *this;
-}
-
-
-G4bool G4ExcitationHandler::operator==(const G4ExcitationHandler &) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator==: is meant to not be accessable! ");
-  return false;
-} 
-
-G4bool G4ExcitationHandler::operator!=(const G4ExcitationHandler &) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator!=: is meant to not be accessable! ");
-  return true;
-}
-
-////////////////////////////////////////////////////////////////////////////////////////////////
-/// 25/07/08 16:45  Proposed by MAC ////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////
-
-G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
-{	
   //for inverse cross section choice
   theEvaporation->SetOPTxs(OPTxs);
   //for the choice of superimposed Coulomb Barrier for inverse cross sections
   theEvaporation->UseSICB(useSICB);
-  
-  // Pointer which will be used to return the final production vector
-  //G4FragmentVector * theResult = new G4FragmentVector;
+  theEvaporation->Initialise();
+}
+
+G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
+{	
   
   // Variables existing until end of method
-  //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*> theEvapList;        // list to apply Evaporation or Fermi Break-Up
   std::list<G4Fragment*> theEvapStableList;  // list to apply PhotonEvaporation
   std::list<G4Fragment*> theResults;         // list to store final result
   std::list<G4Fragment*>::iterator iList;
   //
-  //G4cout << "@@@@@@@@@@ Start G4Exitation Handler @@@@@@@@@@@@@" << G4endl;  
+  //G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@" << G4endl;  
   //G4cout << theInitialState << G4endl;  
   
   // Variables to describe the excited configuration
   G4double exEnergy = theInitialState.GetExcitationEnergy();
-  G4int A = static_cast<G4int>( theInitialState.GetA() +0.5 );
-  G4int Z = static_cast<G4int>( theInitialState.GetZ() +0.5 );
+  G4int A = theInitialState.GetA_asInt();
+  G4int Z = theInitialState.GetZ_asInt();
+
+  G4NistManager* nist = G4NistManager::Instance();
   
   // JMQ 150909:  first step in de-excitation chain (SMM will be used only here)
   
-  // In case A <= 4 the fragment will not perform any nucleon emission
-  if (A <= 4)
+  // In case A <= 1 the fragment will not perform any nucleon emission
+  if (A <= 1)
     {
-      // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later      
-      theEvapStableList.push_back( theInitialStatePtr );
+      theResults.push_back( theInitialStatePtr );
     }
-  
-  else  // If A > 4 we try to apply theFermiModel, theMultiFragmentation or theEvaporation
+  // check if a fragment is stable
+  else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0)
     {
-      
+      theResults.push_back( theInitialStatePtr );
+    }  
+  else  // In all cases apply once theFermiModel, theMultiFragmentation or theEvaporation
+    {      
       // JMQ 150909: first step in de-excitation is treated separately 
       // Fragments after the first step are stored in theEvapList 
-      // Statistical Multifragmentation will take place (just in case) only here
+      // Statistical Multifragmentation will take place only once
       //
-      // Test applicability
-      // Initial State De-Excitation 
-      if(A<GetMaxA()&&Z<GetMaxZ()) 
+      // Fermi Break-Up always first
+      if((A < 5) || (A<GetMaxA() && Z<GetMaxZ())) 
         {
           theTempResult = theFermiModel->BreakItUp(theInitialState);
+	  // fragment was not decaied try to evaporate
+	  if(1 == theTempResult->size()) {
+            if((*theTempResult)[0] !=  theInitialStatePtr) { 
+	      delete (*theTempResult)[0]; 
+	    }
+            delete theTempResult;
+	    theTempResult = theEvaporation->BreakItUp(theInitialState);
+	  }
         }
       else   if (exEnergy>GetMinE()*A) 
         {
           theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
+	  // fragment was not decaied try to evaporate
+	  if(1 == theTempResult->size()) {
+            if((*theTempResult)[0] !=  theInitialStatePtr) { 
+	      delete (*theTempResult)[0]; 
+	    }
+            delete theTempResult;
+	    theTempResult = theEvaporation->BreakItUp(theInitialState);
+	  }
         }
       else 
@@ -174,16 +172,28 @@
       
       G4bool deletePrimary = true;
-      if(theTempResult->size() > 0) 
+      G4int nsec=theTempResult->size();
+      if(nsec > 0) 
 	{      
-	  // Store original state in theEvapList
+	  // Sort out secondary fragments
 	  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
+	      A = (*j)->GetA_asInt();  
+
+	      if(A <= 1) { theResults.push_back(*j); }   // gamma, p, n
+	      else if(1 == nsec) { theEvapStableList.push_back(*j); } // evaporation is not possible  
+	      else  // Analyse fragment
+		{
+		  G4double exEnergy = (*j)->GetExcitationEnergy();
+		  if(exEnergy < minExcitation) {
+		    Z = (*j)->GetZ_asInt(); 
+                    if(nist->GetIsotopeAbundance(Z, A) > 0.0) { 
+		      theResults.push_back(*j); // stable fragment 
+		    } else {   
+		      theEvapList.push_back(*j); // unstable cold fragment 
+		    }
+		  } else { theEvapList.push_back(*j); } // hot fragment
+		}
 	    }
 	}
@@ -204,57 +214,53 @@
   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);
+      A = (*iList)->GetA_asInt(); 
+      Z = (*iList)->GetZ_asInt();
 	  
-      // In case A <= 4 the fragment will not perform any nucleon emission
-      if (A <= 4)
+      // Fermi Break-Up 
+      if ((A < 5) || (A < GetMaxA() && Z < GetMaxZ())) 
 	{
-	  // storing G4Fragment* in theEvapStableList to apply thePhotonEvaporation later    
-	  theEvapStableList.push_back(*iList );
+	  theTempResult = theFermiModel->BreakItUp(*(*iList));
+	  // fragment was not decaied try to evaporate
+	  if(1 == theTempResult->size()) {
+            if((*theTempResult)[0] !=  theInitialStatePtr) { delete (*theTempResult)[0]; }
+            delete theTempResult;
+	    theTempResult = theEvaporation->BreakItUp(*(*iList));
+	  }
 	}
-	  
-      else  // If A > 4 we try to apply theFermiModel or theEvaporation
-	{   
-	  // stable fragment  
-	  if ((*iList)->GetExcitationEnergy() <= 0.1*eV) 
-	    { 
-	      theResults.push_back(*iList); 
+      else // apply Evaporation in another case
+	{
+	  theTempResult = theEvaporation->BreakItUp(*(*iList));
+	}
+      
+      G4bool deletePrimary = true;
+      G4int nsec = theTempResult->size();
+		  
+      // Sort out secondary fragments
+      if ( nsec > 0 )
+	{
+	  G4FragmentVector::iterator j;
+	  for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
+	    {
+	      if((*j) == (*iList)) { deletePrimary = false; }
+	      A = (*j)->GetA_asInt();
+
+	      if(A <= 1) { theResults.push_back(*j); }                // gamma, p, n
+	      else if(1 == nsec) { theEvapStableList.push_back(*j); } // no evaporation 
+	      else  // Analyse fragment
+		{
+		  G4double exEnergy = (*j)->GetExcitationEnergy();
+		  if(exEnergy < minExcitation) {
+		    Z = (*j)->GetZ_asInt();
+                    if(nist->GetIsotopeAbundance(Z, A) > 0.0) { 
+		      theResults.push_back(*j); // stable fragment 
+		    } else {   
+		      theEvapList.push_back(*j); // unstable cold fragment 
+		    }
+		  } else { theEvapList.push_back(*j); } // hot fragment
+		}
 	    }
-	  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)
+	}
+      if( deletePrimary ) { delete (*iList); }
+      delete theTempResult;
     } // end of the loop over theEvapList
 
@@ -267,67 +273,25 @@
   // -----------------------
   
+  // normally should not reach this point - it is kind of work around	      
   for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); ++iList)
     {
-      // take out stable particles and fragments
-      A = static_cast<G4int>((*iList)->GetA()+0.5);
-      if ( A <= 1 )                                       { theResults.push_back(*iList); }
-      else if ((*iList)->GetExcitationEnergy() <= 0.1*eV) { theResults.push_back(*iList); }
-
-      else
+      // photon-evaporation is applied
+      theTempResult = thePhotonEvaporation->BreakUpFragment(*iList);	  
+      G4int nsec = theTempResult->size();
+	  
+      // if there is a gamma emission then
+      if (nsec > 0)
 	{
-	  // photon-evaporation is applied
-	  theTempResult = thePhotonEvaporation->BreakItUp(*(*iList));
-	  
-	  G4bool deletePrimary = true;
-	  G4int nsec = theTempResult->size();
-	  
-	  // if there is a gamma emission then
-	  if (nsec > 1)
+	  G4FragmentVector::iterator j;
+	  for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
 	    {
-	      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); }      
-		}
+	      theResults.push_back(*j); 
 	    }
-	      
-	  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) 
-		{
-		  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;
-	}  
+	}
+      delete theTempResult;
+
+      // priamry fragment is kept
+      theResults.push_back(*iList); 
+
     } // end of photon-evaporation loop
 
@@ -356,6 +320,6 @@
       theFragmentMomentum = (*i)->GetMomentum();
       G4ParticleDefinition* theKindOfFragment = 0;
-      if (theFragmentA == 0 && theFragmentZ == 0) {       // photon
-	theKindOfFragment = G4Gamma::GammaDefinition();   
+      if (theFragmentA == 0) {       // photon or e-
+	theKindOfFragment = (*i)->GetParticleDefinition();   
       } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
 	theKindOfFragment = G4Neutron::NeutronDefinition();
Index: trunk/source/processes/hadronic/models/de_excitation/management/include/G4VEvaporationChannel.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/management/include/G4VEvaporationChannel.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/management/include/G4VEvaporationChannel.hh	(revision 1315)
@@ -24,17 +24,18 @@
 // ********************************************************************
 //
-//
-// $Id: G4VEvaporationChannel.hh,v 1.4 2008/09/19 13:32:54 ahoward Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4VEvaporationChannel.hh,v 1.6 2010/05/11 11:16:04 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
 // by V. Lara (Oct 1998)
 //
-// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
-// cross section option
-// JMQ (06 September 2008) Also external choices have been added for 
-// superimposed Coulomb barrier (if useSICB is set true, by default is false) 
-
-
+// Modified:
+// 03.09.2008 (J.M.Quesada) for external choice of inverse cross section option
+// 06.09.2008 (J.M.Quesada) external choices have been added for superimposed 
+//                          Coulomb barrier (if useSICB is set true, by default is false) 
+// 24.04.2010 (V.Ivanchenko) moved constructor and destructor to source; added two
+//                          new virtual methods EmittedFragment(s) to allow more optimal
+//                          work with G4Fragment objects
+//                          
 
 #ifndef G4VEvaporationChannel_h
@@ -47,24 +48,33 @@
 {
 public:
-  G4VEvaporationChannel() : Name("Anonymous") {};
-  G4VEvaporationChannel(const G4String & aName) : Name(aName) {};
-  G4VEvaporationChannel(const G4String * aName) : Name(*aName) {};
-  virtual ~G4VEvaporationChannel() {};
+
+  G4VEvaporationChannel(const G4String & aName = "Anonymous");
+  virtual ~G4VEvaporationChannel();
 
 private:
+
   G4VEvaporationChannel(const G4VEvaporationChannel & right);
+  const G4VEvaporationChannel & operator=(const G4VEvaporationChannel & right);
 
-  const G4VEvaporationChannel & operator=(const G4VEvaporationChannel & right);
 public:
+
   G4bool operator==(const G4VEvaporationChannel & right) const;
   G4bool operator!=(const G4VEvaporationChannel & right) const;
 
-public:
   virtual void Initialize(const G4Fragment & fragment) = 0;
 
+  // return emitted fragment, initial fragment is modified
+  // and not deleted
+  virtual G4Fragment* EmittedFragment(G4Fragment* theNucleus);
+
+  // return vector of emitted fragments, initial fragment is modified
+  // but not included in this vector
+  virtual G4FragmentVector* BreakUpFragment(G4Fragment* theNucleus);
+
+  // old method initial fragment is not modified, its copy included 
+  // in the list of emitted fragments
   virtual G4FragmentVector * BreakUp(const G4Fragment & theNucleus) = 0;
 
   virtual G4double GetEmissionProbability(void) const = 0;
-
 
   G4String GetName() const {return Name;}
@@ -75,8 +85,9 @@
   // for superimposed Coulomb Barrier for inverse cross sections 	
   inline void UseSICB(G4bool use) { useSICB = use; }	
+
 protected:
+
   G4int OPTxs;
   G4bool useSICB;
-
 
 private:
Index: trunk/source/processes/hadronic/models/de_excitation/management/include/G4VEvaporationFactory.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/management/include/G4VEvaporationFactory.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/management/include/G4VEvaporationFactory.hh	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4VEvaporationFactory.hh,v 1.4 2008/09/19 13:32:54 ahoward Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4VEvaporationFactory.hh,v 1.5 2010/04/27 11:43:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -42,5 +42,5 @@
 {
 public:
-  G4VEvaporationFactory() : _channel(0) {};
+  G4VEvaporationFactory();
   virtual ~G4VEvaporationFactory();
 
Index: trunk/source/processes/hadronic/models/de_excitation/management/src/G4VEvaporationChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/management/src/G4VEvaporationChannel.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/management/src/G4VEvaporationChannel.cc	(revision 1315)
@@ -24,41 +24,52 @@
 // ********************************************************************
 //
-//
-// $Id: G4VEvaporationChannel.cc,v 1.5 2006/06/29 20:23:55 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4VEvaporationChannel.cc,v 1.6 2010/04/25 18:43:08 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
 // by V. Lara (Oct 1998)
 //
+// Modified:
+// 24.04.2010 (V.Ivanchenko) moved constructor and destructor to source; added two
+//                          new virtual methods EmittedFragment(s) to allow more optimal
+//                          work with G4Fragment objects; removed unnecesary exceptions
 
 #include "G4VEvaporationChannel.hh"
 #include "G4HadronicException.hh"
 
-G4VEvaporationChannel::G4VEvaporationChannel(const G4VEvaporationChannel &)
-{
- throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporationChannel::copy_constructor meant to not be accessable");
-}
+G4VEvaporationChannel::G4VEvaporationChannel(const G4String & aName) 
+  : Name(aName) 
+{}
 
+G4VEvaporationChannel::~G4VEvaporationChannel() 
+{}
 
-
-
-const G4VEvaporationChannel & G4VEvaporationChannel::operator=(const G4VEvaporationChannel &)
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporationChannel::operator= meant to not be accessable");
-  return *this;
-}
-
+//G4VEvaporationChannel::G4VEvaporationChannel(const G4VEvaporationChannel &)
+//{
+// throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporationChannel::copy_constructor meant to not be accessable");
+//}
+//const G4VEvaporationChannel & G4VEvaporationChannel::operator=(const G4VEvaporationChannel &)
+//{
+//  throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporationChannel::operator= meant to not be accessable");
+//  return *this;
+//}
 
 G4bool G4VEvaporationChannel::operator==(const G4VEvaporationChannel &right) const
 {
-    return (this == (G4VEvaporationChannel *) &right);
-    //  return false;
+  return (this == (G4VEvaporationChannel *) &right);
 }
 
 G4bool G4VEvaporationChannel::operator!=(const G4VEvaporationChannel &right) const
 {
-    return (this != (G4VEvaporationChannel *) &right);
-    //  return true;
+  return (this != (G4VEvaporationChannel *) &right);
 }
 
+G4Fragment* G4VEvaporationChannel::EmittedFragment(G4Fragment*)
+{
+  return 0;
+}
 
+G4FragmentVector* G4VEvaporationChannel::BreakUpFragment(G4Fragment*)
+{
+  return 0;
+}
Index: trunk/source/processes/hadronic/models/de_excitation/management/src/G4VEvaporationFactory.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/management/src/G4VEvaporationFactory.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/management/src/G4VEvaporationFactory.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4VEvaporationFactory.cc,v 1.6 2008/09/19 13:32:54 ahoward Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4VEvaporationFactory.cc,v 1.7 2010/04/27 11:43:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -32,29 +32,7 @@
 
 #include "G4VEvaporationFactory.hh"
-#include "G4HadronicException.hh"
 
-const G4VEvaporationFactory & 
-G4VEvaporationFactory::operator=(const G4VEvaporationFactory & )
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporationFactory::operator= meant to not be accessable.");
-  return *this;
-}
-
-G4bool 
-G4VEvaporationFactory::operator==(const G4VEvaporationFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporationFactory::operator== meant to not be accessable.");
-  return false;
-}
-
-G4bool 
-G4VEvaporationFactory::operator!=(const G4VEvaporationFactory & ) const
-{
-  throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporationFactory::operator!= meant to not be accessable.");
-  return true;
-}
-
-
-
+G4VEvaporationFactory::G4VEvaporationFactory() : _channel(0) 
+{}
 
 G4VEvaporationFactory::~G4VEvaporationFactory()
Index: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFFragment.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFFragment.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFFragment.hh	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4StatMFFragment.hh,v 1.3 2006/06/29 20:24:07 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4StatMFFragment.hh,v 1.5 2010/05/03 16:49:19 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -36,4 +36,6 @@
 #include "G4StatMFParameters.hh"
 #include "G4ThreeVector.hh"
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
 #include "G4Fragment.hh"
 
Index: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFMacroTemperature.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFMacroTemperature.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFMacroTemperature.hh	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4StatMFMacroTemperature.hh,v 1.3 2006/06/29 20:24:21 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4StatMFMacroTemperature.hh,v 1.4 2010/04/16 17:04:08 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -48,18 +48,7 @@
 			     const G4double ExEnergy, const G4double FreeE0, 
 			     const G4double kappa, 
-			     std::vector<G4VStatMFMacroCluster*> * ClusterVector) :
-	theA(anA),
-	theZ(aZ),
-	_ExEnergy(ExEnergy),
-	_FreeInternalE0(FreeE0),
-	_Kappa(kappa),
-	_MeanMultiplicity(0.0),
-	_MeanTemperature(0.0),
-	_ChemPotentialMu(0.0),
-	_ChemPotentialNu(0.0),
-	_theClusters(ClusterVector) 
-	{};
+			     std::vector<G4VStatMFMacroCluster*> * ClusterVector);
 	
-    ~G4StatMFMacroTemperature() {};
+    ~G4StatMFMacroTemperature();
    
     G4double operator()(const G4double T)
@@ -67,6 +56,7 @@
 
 private:
+
     // Default constructor
-    G4StatMFMacroTemperature() {};
+    G4StatMFMacroTemperature();
 
     // copy constructor
Index: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroTemperature.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroTemperature.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroTemperature.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4StatMFMacroTemperature.cc,v 1.7 2008/11/19 14:33:31 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4StatMFMacroTemperature.cc,v 1.8 2010/04/16 17:04:08 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -36,6 +36,30 @@
 //          Moscow, pshenich@fias.uni-frankfurt.de) make algorithm closer to
 //          original MF model
+// 16.04.10 V.Ivanchenko improved logic of solving equation for tempetature
+//          to protect code from rare unwanted exception; moved constructor 
+//          and destructor to source  
 
 #include "G4StatMFMacroTemperature.hh"
+
+G4StatMFMacroTemperature::G4StatMFMacroTemperature(const G4double anA, const G4double aZ, 
+  const G4double ExEnergy, const G4double FreeE0, const G4double kappa, 
+  std::vector<G4VStatMFMacroCluster*> * ClusterVector) :
+  theA(anA),
+  theZ(aZ),
+  _ExEnergy(ExEnergy),
+  _FreeInternalE0(FreeE0),
+  _Kappa(kappa),
+  _MeanMultiplicity(0.0),
+  _MeanTemperature(0.0),
+  _ChemPotentialMu(0.0),
+  _ChemPotentialNu(0.0),
+  _theClusters(ClusterVector) 
+{}
+
+G4StatMFMacroTemperature::G4StatMFMacroTemperature() 
+{}
+	
+G4StatMFMacroTemperature::~G4StatMFMacroTemperature() 
+{}
 
 // operators definitions
@@ -82,5 +106,5 @@
 
     G4int iterations = 0;  
-    while (fTa < 0.0 && iterations++ < 10) {
+    while (fTa < 0.0 && ++iterations < 10) {
 	Ta -= 0.5*Ta;
 	fTa = this->operator()(Ta);
@@ -89,5 +113,5 @@
     iterations = 0;  
     while (fTa*fTb > 0.0 && iterations++ < 10) {
-	Tb += 2.*std::abs(Tb-Ta);
+	Tb += 2.*std::fabs(Tb-Ta);
 	fTb = this->operator()(Tb);
     }
@@ -96,5 +120,5 @@
       G4cerr <<"G4StatMFMacroTemperature:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
       G4cerr <<"G4StatMFMacroTemperature:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
-	throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution.");
+      throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution.");
     }
 
@@ -102,6 +126,6 @@
     theSolver->SetIntervalLimits(Ta,Tb);
     if (!theSolver->Crenshaw(*this)){ 
-      G4cerr <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
-      G4cerr <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
+      G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
+      G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
     }
     _MeanTemperature = theSolver->GetRoot();
@@ -111,24 +135,28 @@
     // Verify if the root is found and it is indeed within the physical domain, 
     // say, between 1 and 50 MeV, otherwise try Brent method:
-    if (_MeanTemperature < 1. || _MeanTemperature > 50. || std::abs(FunctionValureAtRoot) > 5.e-2) {
-    G4cout << "Crenshaw method failed; function = " << FunctionValureAtRoot << " solution? = " << _MeanTemperature << " MeV " << G4endl;
-    G4Solver<G4StatMFMacroTemperature> * theSolverBrent = new G4Solver<G4StatMFMacroTemperature>(200,1.e-3);
-    theSolverBrent->SetIntervalLimits(Ta,Tb);
-    if (!theSolverBrent->Brent(*this)){
-      G4cerr <<"G4StatMFMacroTemperature, Brent method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
-      G4cerr <<"G4StatMFMacroTemperature, Brent method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl; 
+    if (std::fabs(FunctionValureAtRoot) > 5.e-2) {
+      if (_MeanTemperature < 1. || _MeanTemperature > 50.) {
+	G4cout << "Crenshaw method failed; function = " << FunctionValureAtRoot 
+	       << " solution? = " << _MeanTemperature << " MeV " << G4endl;
+	G4Solver<G4StatMFMacroTemperature> * theSolverBrent = new G4Solver<G4StatMFMacroTemperature>(200,1.e-3);
+	theSolverBrent->SetIntervalLimits(Ta,Tb);
+	if (!theSolverBrent->Brent(*this)){
+	  G4cout <<"G4StatMFMacroTemperature, Brent method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
+	  G4cout <<"G4StatMFMacroTemperature, Brent method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl; 
+	  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
+	}
+
+	_MeanTemperature = theSolverBrent->GetRoot();
+	FunctionValureAtRoot =  this->operator()(_MeanTemperature);
+	delete theSolverBrent;
+      }
+      if (std::abs(FunctionValureAtRoot) > 5.e-2) {
+	//if (_MeanTemperature < 1. || _MeanTemperature > 50. || std::abs(FunctionValureAtRoot) > 5.e-2) {
+	G4cout << "Brent method failed; function = " << FunctionValureAtRoot << " solution? = " << _MeanTemperature << " MeV " << G4endl;
 	throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
-    }
-
-    _MeanTemperature = theSolverBrent->GetRoot();
-    FunctionValureAtRoot =  this->operator()(_MeanTemperature);
-    delete theSolverBrent;
-
-     if (_MeanTemperature < 1. || _MeanTemperature > 50. || std::abs(FunctionValureAtRoot) > 5.e-2) {
-    G4cout << "Brent method failed; function = " << FunctionValureAtRoot << " solution? = " << _MeanTemperature << " MeV " << G4endl;
-	throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
-     }
-    }
-
+      }
+    }
+    //G4cout << "G4StatMFMacroTemperature::CalcTemperature: function = " << FunctionValureAtRoot 
+    //	   << " T(MeV)= " << _MeanTemperature << G4endl;
     return _MeanTemperature;
 }
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4ContinuumGammaDeexcitation.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4ContinuumGammaDeexcitation.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4ContinuumGammaDeexcitation.hh	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4ContinuumGammaDeexcitation.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -64,5 +66,5 @@
 
   // Destructor
-  ~G4ContinuumGammaDeexcitation();
+  virtual ~G4ContinuumGammaDeexcitation();
 
   // Functions
@@ -72,5 +74,5 @@
   virtual G4VGammaTransition* CreateTransition();
 
-  virtual G4bool CanDoTransition() const;
+  virtual G4bool CanDoTransition();
 
 private:
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaDeexcitation.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaDeexcitation.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaDeexcitation.hh	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4DiscreteGammaDeexcitation.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -51,4 +53,5 @@
 // -------------------------------------------------------------------
 //
+
 #ifndef G4DiscreteGammaDeexcitation_hh
 #define G4DiscreteGammaDeexcitation_hh
@@ -71,5 +74,5 @@
 
   // Destructor
-  ~G4DiscreteGammaDeexcitation();
+  virtual ~G4DiscreteGammaDeexcitation();
 
   // Functions
@@ -79,5 +82,5 @@
   virtual G4VGammaTransition * CreateTransition();
 
-  virtual G4bool CanDoTransition() const;
+  virtual G4bool CanDoTransition();
 
   void SetICM(G4bool hl) { _icm = hl; };
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaTransition.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaTransition.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaTransition.hh	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4DiscreteGammaTransition.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -52,5 +54,10 @@
 //        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
 //              Added creation time evaluation for products of evaporation
+//
 //      
+//        19 April 2010, J. M. Quesada. 
+//              Corrections added for taking into account mismatch between tabulated 
+//              gamma energies and level energy differences (fake photons eliminated) 
+//   
 // -------------------------------------------------------------------
 
@@ -62,4 +69,7 @@
 #include "G4NuclearLevel.hh"
 
+//JMQ 180410
+class G4NuclearLevelManager;
+
 class G4DiscreteGammaTransition : public G4VGammaTransition
 {
@@ -68,5 +78,7 @@
   // Constructor
   G4DiscreteGammaTransition(const G4NuclearLevel& level);
-  G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z);
+//JMQ 180410
+//  G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z);
+  G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z, G4int A);
 
   // Destructor
@@ -99,4 +111,10 @@
   G4double _excitation;
   G4double _gammaCreationTime;
+  //JMQ 180410
+  G4int _A;
+  G4int _Z;
+  G4NuclearLevelManager * _levelManager;
+  //JMQ 190410
+  G4double _tolerance;
 };
 
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4E1Probability.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4E1Probability.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4E1Probability.hh	(revision 1315)
@@ -24,4 +24,16 @@
 // ********************************************************************
 //
+// $Id: G4E1Probability.hh,v 1.5 2010/05/19 10:21:44 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
+//
+//---------------------------------------------------------------------
+//
+// Geant4 header G4E1Probability
+//
+// by V. Lara (May 2003)
+//
+// Modifications:
+// 18.05.2010 V.Ivanchenko trying to speedup the most slow method
+//            by usage of G4Pow, integer A and introduction of const members
 //
 //
@@ -33,5 +45,6 @@
 #include "G4VEmissionProbability.hh"
 #include "G4Fragment.hh"
-#include "G4VLevelDensityParameter.hh"
+
+class G4Pow;
 
 class G4E1Probability : public G4VEmissionProbability
@@ -40,7 +53,7 @@
 public:
 
-  G4E1Probability() {};
+  G4E1Probability();
 
-  ~G4E1Probability();
+  virtual ~G4E1Probability();
 
   G4double EmissionProbability(const G4Fragment& frag, const G4double excite);
@@ -48,6 +61,4 @@
 
 private:
-
-  // G4E1Probability() {};
 
   G4E1Probability(const G4E1Probability& right);
@@ -58,10 +69,12 @@
 
   // Integrator (simple Gaussian quadrature)
-
   G4double EmissionIntegration(const G4Fragment& frag, const G4double excite,
                                const G4double lowLim, const G4double upLim,
                                const G4int numIters);
 
-  // G4VLevelDensityParameter* _levelDensity; // Don't need this
+  // members
+  G4Pow*   fG4pow;
+  G4double theLevelDensityParameter;
+  G4double normC;
 
 };
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4PhotonEvaporation.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4PhotonEvaporation.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4PhotonEvaporation.hh	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4PhotonEvaporation.hh,v 1.7 2010/05/11 11:22:14 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -36,7 +38,7 @@
 //      Creation date: 23 October 1998
 //
-//      Modifications:
+//Modifications:
 //
-//        18 October 2002, Fan Lei (flei@space.qinetiq.com)
+// 18 October 2002, Fan Lei (flei@space.qinetiq.com)
 //   
 //        Implementation of Internal Convertion process in discrete deexcitation
@@ -54,4 +56,7 @@
 //            G4ElectronOccupancy _eOccupancy;
 //            G4int _vShellNumber;
+//
+// 11 May 2010, V.Ivanchenko added EmittedFragment and BreakUpFragment
+//                           methods
 // 
 // -------------------------------------------------------------------
@@ -61,5 +66,4 @@
 
 #include "globals.hh"
-#include "G4VPhotonEvaporation.hh"
 #include "G4VEvaporationChannel.hh"
 #include "G4VEmissionProbability.hh"
@@ -67,9 +71,7 @@
 #include "G4ElectronOccupancy.hh"
 
-//#define debug
-
 class G4Fragment;
 
-class G4PhotonEvaporation : public G4VPhotonEvaporation, public G4VEvaporationChannel {
+class G4PhotonEvaporation : public G4VEvaporationChannel {
 
 public:
@@ -79,7 +81,11 @@
     virtual ~G4PhotonEvaporation();
 
+    virtual void Initialize(const G4Fragment & fragment);
+
+    virtual G4Fragment* EmittedFragment(G4Fragment* theNucleus);
+
+    virtual G4FragmentVector* BreakUpFragment(G4Fragment* theNucleus);
+
     virtual G4FragmentVector * BreakItUp(const G4Fragment & nucleus);
-
-    virtual void Initialize(const G4Fragment & fragment);
 
     virtual G4FragmentVector * BreakUp(const G4Fragment & nucleus);
@@ -99,5 +105,4 @@
     void SetEOccupancy( G4ElectronOccupancy  eOccupancy) ;
 
-
     G4ElectronOccupancy GetEOccupancy () { return _eOccupancy;} ;
    
@@ -111,23 +116,20 @@
     G4VGammaDeexcitation * _discrDeexcitation;
     G4VGammaDeexcitation * _contDeexcitation;
-  //    G4VGammaDeexcitation * _cdDeexcitation;
 
     G4ElectronOccupancy _eOccupancy;
     G4int _vShellNumber;
 
-    G4Fragment _nucleus;
+    G4Fragment* _nucleus;
     G4double _gammaE;
 
     G4PhotonEvaporation(const G4PhotonEvaporation & right);
-
     const G4PhotonEvaporation & operator = (const G4PhotonEvaporation & right);
 
-    // MGP - Check == and != multiple inheritance... must be a mess!
     G4bool operator == (const G4PhotonEvaporation & right) const;
     G4bool operator != (const G4PhotonEvaporation & right) const;
 
-#ifdef debug
-    void CheckConservation(const G4Fragment & theInitialState, G4FragmentVector * Result) const;
-#endif
+  //#ifdef debug
+  //  void CheckConservation(const G4Fragment & theInitialState, G4FragmentVector * Result) const;
+  //#endif
 
 
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VGammaDeexcitation.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VGammaDeexcitation.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VGammaDeexcitation.hh	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4VGammaDeexcitation.hh,v 1.8 2010/04/28 14:22:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -37,4 +39,15 @@
 //
 //      Modifications:
+//          
+//        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
+//              Added creation time evaluation for products of evaporation
+//
+//        21 Nov 2001, Fan Lei (flei@space.qinetiq.com)
+//           Modified GenerateGamma() and UpdateUncleus() for implementation
+//           of Internal Conversion processs
+// 
+//        8 March 2002, Fan Lei (flei@space.qinetiq.com)
+//          Added  SetEO () , GetEO(), UpdateElectrons() to allow the assignment 
+//          and modification of electron configuration.
 //
 //        18 October 2002, F. Lei
@@ -43,17 +56,9 @@
 //          Added SetVaccantSN(). It is need to to re-set _vSN after each 
 //          IC happened.
-// 
-//        8 March 2002, Fan Lei (flei@space.qinetiq.com)
-//          Added  SetEO () , GetEO(), UpdateElectrons() to allow the assignment 
-//          and modification of electron configuration.
-//          
-//        21 Nov 2001, Fan Lei (flei@space.qinetiq.com)
-//           Modified GenerateGamma() and UpdateUncleus() for implementation
-//           of Internal Conversion processs
 //
-//        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
-//              Added creation time evaluation for products of evaporation
+//        28 April 2010, V.Ivanchenko cleanup methods
 //
 // -------------------------------------------------------------------
+//
 
 #ifndef G4VGAMMADEEXCITATION_HH
@@ -70,53 +75,76 @@
 public:
 
-    G4VGammaDeexcitation();
+  G4VGammaDeexcitation();
 
-    virtual ~G4VGammaDeexcitation();
+  virtual ~G4VGammaDeexcitation();
 
-    virtual G4VGammaTransition * CreateTransition() = 0;
-    virtual G4bool CanDoTransition() const = 0;
+  virtual G4VGammaTransition * CreateTransition() = 0;
+  virtual G4bool CanDoTransition() = 0;
 
-    // Single gamma transition
-    virtual G4FragmentVector * DoTransition();
+  // Single gamma transition
+  G4FragmentVector * DoTransition();
 
-    // Chain of gamma transitions
-    virtual G4FragmentVector * DoChain();
+  // Chain of gamma transitions
+  G4FragmentVector * DoChain();
 
-    virtual G4Fragment * GenerateGamma();
+  G4Fragment * GenerateGamma();
 
-    virtual const G4Fragment & GetNucleus() const;
+  inline G4Fragment* GetNucleus();
 
-    virtual void SetNucleus(const G4Fragment & nucleus);
+  inline void SetNucleus(G4Fragment* nucleus);
 
-    virtual void SetVerboseLevel(G4int verbose);
+  inline void SetVerboseLevel(G4int verbose);
 
-    void SetEO(G4ElectronOccupancy eo) { _electronO = eo; };
-    void SetVaccantSN( G4int val ) { _vSN = val;};
+  inline void Initialize();
+
+  void SetEO(G4ElectronOccupancy eo) { _electronO = eo; };
+  void SetVaccantSN( G4int val ) { _vSN = val;};
   
-    G4ElectronOccupancy GetEO() { return _electronO; };    
-    G4int GetVacantSN() {return _vSN;};
+  G4ElectronOccupancy GetEO() { return _electronO; };    
+  G4int GetVacantSN() {return _vSN;};
+
 protected:
 
-    void Initialize();
-    void UpdateNucleus(const G4Fragment * gamma);
-    void UpdateElectrons ();
-    void Update();
+  void Update();
 
-    G4VGammaTransition * _transition; // Owned pointer
-    G4int _verbose;
+  G4VGammaTransition* _transition; // Owned pointer
+  G4int _verbose;
 
 private:
 
-    G4Fragment _nucleus;
-    G4ElectronOccupancy _electronO;
-    G4int _vSN;
+  G4Fragment* _nucleus;
+  G4ElectronOccupancy _electronO;
+  G4int _vSN;
 
-    G4VGammaDeexcitation(const G4VGammaDeexcitation & right);
-
-    const G4VGammaDeexcitation & operator = (const G4VGammaDeexcitation & right);
-    G4bool operator == (const G4VGammaDeexcitation & right) const;
-    G4bool operator != (const G4VGammaDeexcitation & right) const;
+  G4VGammaDeexcitation(const G4VGammaDeexcitation & right);
+  const G4VGammaDeexcitation & operator = (const G4VGammaDeexcitation & right);
+  G4bool operator == (const G4VGammaDeexcitation & right) const;
+  G4bool operator != (const G4VGammaDeexcitation & right) const;
 
 };
+
+inline G4Fragment* G4VGammaDeexcitation::GetNucleus() 
+{
+  return _nucleus; 
+}
+
+inline void G4VGammaDeexcitation::SetNucleus(G4Fragment* nucleus)
+{
+  _nucleus = nucleus;
+}
+
+inline void G4VGammaDeexcitation::SetVerboseLevel(G4int verbose)
+{
+  _verbose = verbose;
+}
+
+inline void G4VGammaDeexcitation::Initialize()
+{
+  if (_transition != 0) { delete _transition; }
+  _transition = CreateTransition();
+  if (_transition != 0) {
+    _transition->SetEnergyFrom(_nucleus->GetExcitationEnergy());
+  }
+}
 
 #endif
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VPhotonEvaporation.hh
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VPhotonEvaporation.hh	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VPhotonEvaporation.hh	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4VPhotonEvaporation.hh,v 1.3 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4ContinuumGammaDeexcitation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4ContinuumGammaDeexcitation.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4ContinuumGammaDeexcitation.cc	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4ContinuumGammaDeexcitation.cc,v 1.7 2010/04/30 16:08:03 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -40,4 +42,6 @@
 //
 //      02 May 2003,   Vladimir Ivanchenko change interface to G4NuclearlevelManager
+//
+//      19 April 2010 J. M. Quesada: smaller value of tolerance parameter
 //      
 // -------------------------------------------------------------------
@@ -61,5 +65,6 @@
 //
 
-G4ContinuumGammaDeexcitation::G4ContinuumGammaDeexcitation(): _nucleusZ(0), _nucleusA(0), _levelManager(0)
+G4ContinuumGammaDeexcitation::G4ContinuumGammaDeexcitation()
+  : _nucleusZ(0), _nucleusA(0), _levelManager(0)
 { }
 
@@ -71,8 +76,8 @@
 G4VGammaTransition* G4ContinuumGammaDeexcitation::CreateTransition()
 {
-  G4Fragment nucleus = GetNucleus();
-  G4int Z = static_cast<G4int>(nucleus.GetZ());
-  G4int A = static_cast<G4int>(nucleus.GetA());
-  G4double excitation = nucleus.GetExcitationEnergy();
+  G4Fragment* nucleus = GetNucleus();
+  G4int Z = static_cast<G4int>(nucleus->GetZ());
+  G4int A = static_cast<G4int>(nucleus->GetA());
+  G4double excitation = nucleus->GetExcitationEnergy();
 
   if (_nucleusA != A || _nucleusZ != Z)
@@ -83,7 +88,7 @@
     }
 
-  if (_verbose > 1)
+  if (_verbose > 1) {
     G4cout << "G4ContinuumGammaDeexcitation::CreateTransition - Created" << G4endl;
-
+  }
   G4VGammaTransition* gt =  new G4ContinuumGammaTransition(_levelManager,Z,A,excitation,_verbose );
 
@@ -92,60 +97,61 @@
     
 
-G4bool G4ContinuumGammaDeexcitation::CanDoTransition() const
+G4bool G4ContinuumGammaDeexcitation::CanDoTransition() 
 {
- G4bool canDo = true;
+  //JMQ: far too small, creating sometimes continuum gammas instead of the right discrete ones
+  // (when excitation energy is slightly over maximum discrete  energy): changed
+  //  G4double tolerance = 10*eV;
+  const G4double tolerance = CLHEP::keV;
 
   if (_transition == 0) 
     {
-      canDo = false;
-
-      if (_verbose > 0)
-	G4cout 
-	  << "G4ContinuumGammaDeexcitation::CanDoTransition - Null transition "
-	  << G4endl;
+      if (_verbose > 0) {
+	G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition - Null transition "
+	       << G4endl;
+      }
+      return false;
     }
 
-  G4Fragment nucleus = GetNucleus();
-  G4double excitation = nucleus.GetExcitationEnergy();
+  G4Fragment* nucleus = GetNucleus();
+  G4double excitation = nucleus->GetExcitationEnergy();
 
-  G4double A = nucleus.GetA();
-  G4double Z = nucleus.GetZ();
-  if (A <2 || Z<3)
+  //  G4int A = (G4int)nucleus->GetA();
+  // G4int Z = (G4int)nucleus->GetZ();
+  if (_nucleusA<2 || _nucleusZ<3)
     {
-      canDo = false;
-      if (_verbose > 0) 
-	G4cout 
-	  << "G4ContinuumGammaDeexcitation::CanDoTransition - n/p/H"
-	  << G4endl;
+      if (_verbose > 1) { 
+	G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition - n/p/H"
+	       << G4endl;
+      }
+      return false;
     }
 
-
-
-  if (excitation <= 0.) 
+  if (excitation <= tolerance) 
     {
-      canDo = false;
-      if (_verbose > 0) 
-	G4cout 
-	  << "G4ContinuumGammaDeexcitation::CanDoTransition -  Excitation <= 0"
-	  << G4endl;
+      if (_verbose > 1) { 
+	G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition -  Excitation "
+	       << excitation/CLHEP::keV << " keV is too small"
+	       << G4endl;
+      }
+      return false;
     }
-  G4double tolerance = 10*eV;
-  if (excitation <= (_levelManager->MaxLevelEnergy()+ tolerance)) 
-    {
-      canDo = false;  
-      if (_verbose > 0)
-      G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition -  Excitation " 
-	     << excitation << " below max discrete level " 
-	     << _levelManager->MaxLevelEnergy() << G4endl;
+  if (excitation <= (_levelManager->MaxLevelEnergy() + tolerance)) 
+    {  
+      if (_verbose > 0) {
+	G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition -  Excitation " 
+	       << excitation << " below max discrete level " 
+	       << _levelManager->MaxLevelEnergy() << G4endl;
+      }
+      return false;
     }
   
-  if (canDo)
-    { if (_verbose > 1) 
-      G4cout <<"G4ContinuumGammaDeexcitation::CanDoTransition - CanDo" 
-	     << G4endl; 
-    }
-
-  return canDo;
-
+  if (_verbose > 1) {
+    G4cout <<"G4ContinuumGammaDeexcitation::CanDoTransition - CanDo" 
+	   << " Eex(keV)= " << excitation/CLHEP::keV 
+	   << " Emax(keV)= " << _levelManager->MaxLevelEnergy()/CLHEP::keV 
+	   << " Z= " << _nucleusZ << " A= " << _nucleusA
+	   << G4endl;
+  }
+  return true;
 }
 
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaDeexcitation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaDeexcitation.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaDeexcitation.cc	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4DiscreteGammaDeexcitation.cc,v 1.15 2010/05/10 07:20:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -65,17 +67,17 @@
   _rdm(false), _levelManager(0)
 {
-  _tolerance = 0.1 * MeV;
+  _tolerance = CLHEP::keV;
 }
 
-
-G4DiscreteGammaDeexcitation::~G4DiscreteGammaDeexcitation() {}
-
+G4DiscreteGammaDeexcitation::~G4DiscreteGammaDeexcitation() 
+{}
 
 G4VGammaTransition* G4DiscreteGammaDeexcitation::CreateTransition()
 {
-  G4Fragment nucleus = GetNucleus();
-  G4int A = static_cast<G4int>(nucleus.GetA());
-  G4int Z = static_cast<G4int>(nucleus.GetZ());
-
+  G4Fragment* nucleus = GetNucleus();
+  G4int A = static_cast<G4int>(nucleus->GetA());
+  G4int Z = static_cast<G4int>(nucleus->GetZ());
+  //  _verbose =2;
+  //  G4cout << "G4DiscreteGammaDeexcitation::CreateTransition: " << nucleus << G4endl;
   if (_nucleusA != A || _nucleusZ != Z) 
     {
@@ -84,6 +86,4 @@
       _levelManager = G4NuclearLevelStore::GetInstance()->GetManager(Z,A);
     }
-
-
 
   if (_levelManager->IsValid()) 
@@ -96,16 +96,16 @@
 	}
 	
-      G4double excitation = nucleus.GetExcitationEnergy();
-      //      const G4NuclearLevel* level =_levelManager.NearestLevel(excitation, _tolerance);
+      G4double excitation = nucleus->GetExcitationEnergy();
       const G4NuclearLevel* level =_levelManager->NearestLevel(excitation);
 	
       if (level != 0)  
 	{ 
-	  if (_verbose > 0)
+	  if (_verbose > 0) {
 	    G4cout 
 	      << "G4DiscreteGammaDeexcitation::CreateTransition - Created from level energy " 
 	      << level->Energy() << ", excitation is " 
 	      << excitation << G4endl;
-	  G4DiscreteGammaTransition* dtransit = new G4DiscreteGammaTransition(*level,Z);
+	  }
+	  G4DiscreteGammaTransition* dtransit = new G4DiscreteGammaTransition(*level,Z,A);
 	  dtransit->SetICM(_icm);  
 	  return dtransit;
@@ -113,17 +113,17 @@
       else 
 	{ 
-	  if (_verbose > 0)
+	  if (_verbose > 0) {
 	    G4cout 
 	      << "G4DiscreteGammaDeexcitation::CreateTransition - No transition created from "
 	      << excitation << " within tolerance " << _tolerance << G4endl;
-	    
+	  }
 	  return 0; 
 	}
     }
-  else return 0;
+  return 0;
 }
 
 
-G4bool G4DiscreteGammaDeexcitation::CanDoTransition() const
+G4bool G4DiscreteGammaDeexcitation::CanDoTransition() 
 {
 
@@ -138,9 +138,6 @@
 	<< G4endl;
   } 
-  G4Fragment nucleus = GetNucleus();
   if (canDo)  {
-    G4double A = nucleus.GetA();
-    G4double Z = nucleus.GetZ();
-    if (Z<2 || A<3 || Z>98)
+    if (_nucleusZ<2 || _nucleusA<3 || _nucleusZ>98)
       {
 	canDo = false;
@@ -152,19 +149,23 @@
   }
 
-  G4double excitation = nucleus.GetExcitationEnergy();
+  G4Fragment* nucleus = GetNucleus();
+  G4double excitation = nucleus->GetExcitationEnergy();
+  //G4cout << "G4DiscreteGammaDeexcitation::CanDoTransition: " << nucleus << G4endl;
 
   if (canDo) {
-    if (excitation <= 0.) {
+    if (excitation <= _tolerance) {
       canDo = false;
-      if (_verbose > 0) 
+      if (_verbose > 0) {
 	G4cout 
 	  << "G4DiscreteGammaDeexcitation::CanDoTransition -  Excitation <= 0" 
+	  << excitation << "  " << excitation - _tolerance
 	  << G4endl;
+      }
     } else { 
       if (excitation > _levelManager->MaxLevelEnergy() + _tolerance) canDo = false;
-      if (excitation < _levelManager->MinLevelEnergy() - _tolerance) canDo = false;  
+      //if (excitation < _levelManager->MinLevelEnergy() - _tolerance) canDo = false;  
       // The following is a protection to avoid looping in case of elements with very low
       // ensdf levels
-      if (excitation < _levelManager->MinLevelEnergy() * 0.9) canDo = false;  
+      //if (excitation < _levelManager->MinLevelEnergy() * 0.9) canDo = false;  
   
       if (_verbose > 0) {
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaTransition.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaTransition.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaTransition.cc	(revision 1315)
@@ -37,15 +37,23 @@
 //
 //      Modifications:
+//        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
+//              Added creation time evaluation for products of evaporation
+//
+//        21 Nov. 2001, Fan Lei (flei@space.qinetiq.com)
+//              i) added G4int _nucleusZ initialise it through the constructor
+//              ii) modified SelectGamma() to allow the generation of conversion electrons    
+//              iii) added #include G4AtomicShells.hh
+//      
 //        09 Sep. 2002, Fan Lei  (flei@space.qinetiq.com)
 //              Added renormalization to determine whether transition leads to
 //              electron or gamma in SelectGamma()
 //
-//        21 Nov. 2001, Fan Lei (flei@space.qinetiq.com)
-//              i) added G4int _nucleusZ initialise it through the constructor
-//              ii) modified SelectGamma() to allow the generation of conversion electrons      
-//              iii) added #include G4AtomicShells.hh
-//      
-//        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
-//              Added creation time evaluation for products of evaporation
+//        19 April 2010, J. M. Quesada. 
+//              Corrections added for taking into account mismatch between tabulated 
+//              gamma energies and level energy differences (fake photons eliminated) 
+//
+//        9 May 2010, V.Ivanchenko
+//              Removed unphysical corretions of gamma energy; fixed default particle 
+//              as gamma; do not subtract bounding energy in case of electron emmision
 //      
 // -------------------------------------------------------------------
@@ -55,4 +63,7 @@
 #include "G4RandGeneralTmp.hh"
 #include "G4AtomicShells.hh"
+//JMQ: 
+#include "G4NuclearLevelStore.hh"
+#include "G4Pow.hh"
 
 G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level): 
@@ -60,9 +71,13 @@
 { }
 
-G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z): 
+//JMQ: now A is also needed in the constructor
+//G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z): 
+G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z, G4int A): 
   _nucleusZ(Z), _orbitE(-1), _bondE(0.), _aGamma(true), _icm(false), _gammaEnergy(0.), 
-  _level(level), _excitation(0.),  _gammaCreationTime(0.)
+  _level(level), _excitation(0.),  _gammaCreationTime(0.),_A(A),_Z(Z)
 {
   _verbose = 0;
+  //JMQ: added tolerence in the mismatch
+  _tolerance = CLHEP::keV;
 }
 
@@ -74,27 +89,44 @@
 void G4DiscreteGammaTransition::SelectGamma()
 {
-
+  // default gamma 
+  _aGamma = true;    
   _gammaEnergy = 0.;
-
+  
   G4int nGammas = _level.NumberOfGammas();
   if (nGammas > 0)
     {
       G4double random = G4UniformRand();
-
-      G4int iGamma = 0;
-      for(iGamma=0;iGamma < nGammas;iGamma++)
+      
+      G4int iGamma;
+      for(iGamma=0; iGamma<nGammas; ++iGamma)
 	{
 	  if(random <= (_level.GammaCumulativeProbabilities())[iGamma])
-	    break;
+	    { break; }
 	}
-
-
+      
       // Small correction due to the fact that there are mismatches between 
       // nominal level energies and emitted gamma energies
-
-      G4double eCorrection = _level.Energy() - _excitation;
-
+      
+      // 09.05.2010 VI : it is an error ?
+      G4double eCorrection = _level.Energy() - _excitation;      
       _gammaEnergy = (_level.GammaEnergies())[iGamma] - eCorrection;
-
+            
+      //JMQ: 
+      //1)If chosen gamma energy is close enough to excitation energy, the later
+      //  is used instead for gamma dacey to gs (it guarantees energy conservation)
+      //2)For energy conservation, level energy differences instead of  tabulated 
+      //  gamma energies must be used (origin of final fake photons)
+      
+      if(_excitation - _gammaEnergy < _tolerance)      
+	{ 
+	  _gammaEnergy =_excitation;
+	}
+      else
+	{		 
+	  _levelManager = G4NuclearLevelStore::GetInstance()->GetManager(_Z,_A);
+	  _gammaEnergy = _excitation - 
+	    _levelManager->NearestLevel(_excitation - _gammaEnergy)->Energy();
+	}
+      
       //  Warning: the following check is needed to avoid loops:
       //  Due essentially to missing nuclear levels in data files, it is
@@ -106,6 +138,10 @@
       //          but this change needs a more complex revision of actual design.
       //          I leave this for a later revision.
-
-      if (_gammaEnergy < _level.Energy()*10e-5) _gammaEnergy = _excitation;
+      
+      if (_gammaEnergy < _level.Energy()*10.e-5) _gammaEnergy = _excitation;
+
+      //G4cout << "G4DiscreteGammaTransition::SelectGamma: " << _gammaEnergy 
+      //	     << " _icm: " << _icm << G4endl;
+
       // now decide whether Internal Coversion electron should be emitted instead
       if (_icm) {
@@ -144,20 +180,31 @@
 	      }
 	    }
-	    if (_verbose > 0)
+	    _bondE = G4AtomicShells::GetBindingEnergy(_nucleusZ, iShell);
+	    if (_verbose > 0) {
 	      G4cout << "G4DiscreteGammaTransition: _nucleusZ = " <<_nucleusZ 
 		     << " , iShell = " << iShell  
-		     << " , Shell binding energy = " << G4AtomicShells::GetBindingEnergy(_nucleusZ, iShell) / keV
+		     << " , Shell binding energy = " << _bondE/keV
 		     << " keV " << G4endl;
-	    _bondE = G4AtomicShells::GetBindingEnergy(_nucleusZ, iShell);
-	    _gammaEnergy = _gammaEnergy - _bondE;
+	    }
+
+	    // 09.05.2010 VI : it is an error - cannot subtract bond energy from 
+	    //                 transition energy here
+	    //_gammaEnergy = _gammaEnergy - _bondE; 
+	    //G4cout << "_gammaEnergy = " << _gammaEnergy << G4endl;
+
 	    _orbitE = iShell;	  
 	    _aGamma = false ;   // emitted is not a gamma now 
 	  }
       }
-    
-      G4double tau = _level.HalfLife() / std::log(2.0);
-
-      G4double tMin = 0;
-      G4double tMax = 10.0 * tau;
+      
+      G4double tau = _level.HalfLife() / G4Pow::GetInstance()->logZ(2);
+
+      //09.05.2010 VI rewrite samling of decay time 
+      //              assuming ordinary exponential low
+      _gammaCreationTime = 0.;      
+      if(tau > 0.0) {  _gammaCreationTime = -tau*log(G4UniformRand()); }
+
+      //G4double tMin = 0;
+      //G4double tMax = 10.0 * tau;
       //  Original code, not very efficent
       //  G4int nBins = 200;
@@ -172,9 +219,10 @@
       //  G4RandGeneralTmp randGeneral(sampleArray, nBins);
       //G4double random = randGeneral.shoot();
-  
+      
       //_gammaCreationTime = tMin + (tMax - tMin) * random;
 
       // new code by Fan Lei
       //
+      /*
       if (tau != 0 ) 
       {
@@ -186,15 +234,9 @@
 	  //	   << _gammaCreationTime/second << G4endl;
        } else { _gammaCreationTime=0.; }
+      */
     }
   return;
 }
 
-
-//G4bool G4DiscreteGammaTransition::IsAGamma()
-//{
-//  return _aGamma;
-//}
-
-
 G4double G4DiscreteGammaTransition::GetGammaEnergy()
 {
@@ -210,10 +252,9 @@
 {
   _excitation = energy;
-  return;
-}
-
-
-
-
-
-
+}
+
+
+
+
+
+
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4E1Probability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4E1Probability.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4E1Probability.cc	(revision 1315)
@@ -25,42 +25,37 @@
 //
 //
-//  Class G4E1Probability.cc
+// $Id: G4E1Probability.cc,v 1.9 2010/05/25 10:47:24 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
+//
+//---------------------------------------------------------------------
+//
+// Geant4 class G4E1Probability
+//
+// by V. Lara (May 2003)
+//
+// Modifications:
+// 18.05.2010 V.Ivanchenko trying to speedup the most slow method
+//            by usage of G4Pow, integer A and introduction of const members
+//
 //
 
 #include "G4E1Probability.hh"
-#include "G4ConstantLevelDensityParameter.hh"
+//#include "G4ConstantLevelDensityParameter.hh"
 #include "Randomize.hh"
+#include "G4Pow.hh"
 
 // Constructors and operators
 //
 
-G4E1Probability::G4E1Probability(const G4E1Probability& ) : G4VEmissionProbability()
-{
-
-  throw G4HadronicException(__FILE__, __LINE__, "G4E1Probability::copy_constructor meant to not be accessible");
-
-}
-
-const G4E1Probability& G4E1Probability::
-operator=(const G4E1Probability& ) 
-{
-
-  throw G4HadronicException(__FILE__, __LINE__, "G4E1Probability::operator= meant to not be accessible");
-  return *this;
-}
-
-G4bool G4E1Probability::operator==(const G4E1Probability& ) const
-{
-
-  return false;
-
-}
-
-G4bool G4E1Probability::operator!=(const G4E1Probability& ) const
-{
-
-  return true;
-
-}
+G4E1Probability::G4E1Probability():G4VEmissionProbability()
+{
+  G4double x = CLHEP::pi*CLHEP::hbarc;
+  normC = 1.0 / (x*x);
+  theLevelDensityParameter = 0.125/MeV;
+  fG4pow = G4Pow::GetInstance(); 
+}
+
+G4E1Probability::~G4E1Probability()
+{}
 
 // Calculate the emission probability
@@ -68,5 +63,5 @@
 
 G4double G4E1Probability::EmissionProbDensity(const G4Fragment& frag, 
-                                                 const G4double gammaE)
+					      const G4double gammaE)
 {
 
@@ -76,14 +71,12 @@
   // the probability density for photon evaporation from U to U - gammaE
   // (U = nucleus excitation energy, gammaE = total evaporated photon
-  // energy).
-  // fragment = nuclear fragment BEFORE de-excitation
+  // energy). Fragment = nuclear fragment BEFORE de-excitation
 
   G4double theProb = 0.0;
 
-  const G4double Afrag = frag.GetA();
-  const G4double Zfrag = frag.GetZ();
-  const G4double Uexcite = frag.GetExcitationEnergy();
-
-  if( (Uexcite-gammaE) < 0.0 || gammaE < 0 || Uexcite <= 0) return theProb;
+  G4int Afrag = frag.GetA_asInt();
+  G4double Uexcite = frag.GetExcitationEnergy();
+
+  if( (Uexcite-gammaE) < 0.0 || gammaE < 0) { return theProb; }
 
   // Need a level density parameter.
@@ -91,13 +84,14 @@
   // nuclei).
 
-  static G4ConstantLevelDensityParameter a;
-
-  G4double aLevelDensityParam = a.LevelDensityParameter(static_cast<G4int>(Afrag),
-							static_cast<G4int>(Zfrag),
-							Uexcite);
-
-  G4double levelDensBef = std::exp(2.0*std::sqrt(aLevelDensityParam*Uexcite));
-  G4double levelDensAft = std::exp(2.0*std::sqrt(aLevelDensityParam*(Uexcite-gammaE)));
-
+  G4double aLevelDensityParam = Afrag*theLevelDensityParameter;
+
+  //  G4double levelDensBef = std::exp(2*std::sqrt(aLevelDensityParam*Uexcite));
+  //  G4double levelDensAft = std::exp(2*std::sqrt(aLevelDensityParam*(Uexcite-gammaE)));
+  // VI reduce number of calls to exp 
+  G4double levelDens = 1.0;
+  if( aLevelDensityParam > 0.0 ) {
+    levelDens = std::exp(2*(std::sqrt(aLevelDensityParam*(Uexcite-gammaE)) 
+      - std::sqrt(aLevelDensityParam*Uexcite))); 
+  }
   // Now form the probability density
 
@@ -107,9 +101,7 @@
   G4double sigma0 = 2.5 * Afrag * millibarn;  // millibarns
 
-  G4double Egdp = (40.3 / std::pow(Afrag,0.2) )*MeV;
+  G4double Egdp   = (40.3 / fG4pow->powZ(Afrag,0.2) )*MeV;
   G4double GammaR = 0.30 * Egdp;
  
-  static G4double normC = 1.0 / ((pi * hbarc)*(pi * hbarc));
-
   // CD
   //cout<<"  PROB TESTS "<<G4endl;
@@ -124,12 +116,18 @@
   //cout<<" normC = "<<normC<<G4endl;
 
-  G4double numerator = sigma0 * gammaE*gammaE * GammaR*GammaR;
-  G4double denominator = (gammaE*gammaE - Egdp*Egdp)*
-           (gammaE*gammaE - Egdp*Egdp) + GammaR*GammaR*gammaE*gammaE;
-
-  G4double sigmaAbs = numerator/denominator; 
-
-  theProb = normC * sigmaAbs * gammaE*gammaE *
-            levelDensAft/levelDensBef;
+  // VI implementation 18.05.2010
+  G4double gammaE2 = gammaE*gammaE;
+  G4double gammaR2 = gammaE2*GammaR*GammaR;
+  G4double egdp2   = gammaE2 - Egdp*Egdp;
+  G4double sigmaAbs = sigma0*gammaR2/(egdp2*egdp2 + gammaR2); 
+  theProb = normC * sigmaAbs * gammaE2 * levelDens;
+
+  // old implementation
+  //  G4double numerator = sigma0 * gammaE*gammaE * GammaR*GammaR;
+  // G4double denominator = (gammaE*gammaE - Egdp*Egdp)*
+  //         (gammaE*gammaE - Egdp*Egdp) + GammaR*GammaR*gammaE*gammaE;
+
+  //G4double sigmaAbs = numerator/denominator; 
+  //theProb = normC * sigmaAbs * gammaE2 * levelDensAft/levelDensBef;
 
   // CD
@@ -211,5 +209,3 @@
 }
 
-G4E1Probability::~G4E1Probability() {}
-
-
+
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PhotonEvaporation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PhotonEvaporation.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PhotonEvaporation.cc	(revision 1315)
@@ -24,5 +24,6 @@
 // ********************************************************************
 //
-
+// $Id: G4PhotonEvaporation.cc,v 1.13 2010/05/17 11:47:43 flei Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -37,7 +38,7 @@
 //      Creation date: 23 October 1998
 //
-//      Modifications: 
+// Modifications: 
 //      
-//        8 March 2002, Fan Lei (flei@space.qinetiq.com)
+// 8 March 2002, Fan Lei (flei@space.qinetiq.com)
 //   
 //        Implementation of Internal Convertion process in discrete deexcitation
@@ -50,5 +51,9 @@
 //       G4ElectronOccupancy GetEOccupancy () ;
 //
+// 11 May 2010, V.Ivanchenko add implementation of EmittedFragment and 
+//                           BreakUpFragment methods; cleanup logic
+//
 // -------------------------------------------------------------------
+//
 
 #include "G4PhotonEvaporation.hh"
@@ -86,17 +91,103 @@
 void G4PhotonEvaporation::Initialize(const G4Fragment& fragment)
 {
-  _nucleus = fragment;
-  return;
-}
-
+  _nucleus = const_cast<G4Fragment*>(&fragment);
+}
+
+G4Fragment* G4PhotonEvaporation::EmittedFragment(G4Fragment* nucleus)
+{
+  //G4cout << "G4PhotonEvaporation::EmittedFragment" << G4endl;
+  _nucleus = nucleus;
+  
+  // Do one photon emission by the continues deexcitation  
+  _contDeexcitation->SetNucleus(_nucleus);
+  _contDeexcitation->Initialize();
+
+  if(_contDeexcitation->CanDoTransition()) {  
+    G4Fragment* gamma = _contDeexcitation->GenerateGamma();
+    if(gamma) { 
+      if (_verbose > 0) {
+	G4cout << "G4PhotonEvaporation::EmittedFragment continium deex: "   
+	       << gamma << G4endl;
+        G4cout << "   Residual: " << nucleus << G4endl;
+      }
+      return gamma; 
+    }
+  }
+
+  // Do one photon emission by the discrete deexcitation 
+  _discrDeexcitation->SetNucleus(_nucleus);
+  _discrDeexcitation->Initialize();
+
+  if(_discrDeexcitation->CanDoTransition()) {  
+    G4Fragment* gamma = _discrDeexcitation->GenerateGamma();
+    if(gamma) { 
+      if (_verbose > 0) {
+	G4cout << "G4PhotonEvaporation::EmittedFragment discrete deex: "   
+	       << gamma << G4endl;
+        G4cout << "   Residual: " << nucleus << G4endl;
+      }
+      return gamma; 
+    }
+  }
+
+  if (_verbose > 0) {
+    G4cout << "G4PhotonEvaporation unable emit gamma: " 
+	   << nucleus << G4endl;
+  }
+  return 0;
+}
+
+G4FragmentVector* G4PhotonEvaporation::BreakUpFragment(G4Fragment* nucleus)
+{
+  //G4cout << "G4PhotonEvaporation::BreakUpFragment" << G4endl;
+  // The same pointer of primary nucleus
+  _nucleus = nucleus;
+  _contDeexcitation->SetNucleus(_nucleus);
+  _discrDeexcitation->SetNucleus(_nucleus);
+
+  // Do the whole gamma chain 
+  G4FragmentVector* products = _contDeexcitation->DoChain();  
+  if( !products ) { products = new G4FragmentVector(); }
+
+  if (_verbose > 0) {
+    G4cout << "G4PhotonEvaporation::BreakUpFragment " << products->size() 
+	   << " gammas from ContinuumDeexcitation " << G4endl;
+    G4cout << "   Residual: " << nucleus << G4endl;
+  }
+  // Products from discrete gamma transitions
+  G4FragmentVector* discrProducts = _discrDeexcitation->DoChain();
+  if(discrProducts) {
+    _eOccupancy = _discrDeexcitation->GetEO();
+    _vShellNumber = _discrDeexcitation->GetVacantSN();
+
+    // not sure if the following line is needed!
+    _discrDeexcitation->SetVaccantSN(-1);
+
+    if (_verbose > 0) {
+      G4cout << "G4PhotonEvaporation::BreakUpFragment " << discrProducts->size() 
+	     << " gammas from DiscreteDeexcitation " << G4endl;
+      G4cout << "   Residual: " << nucleus << G4endl;
+    }
+    G4FragmentVector::iterator i;
+    for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
+      {
+	products->push_back(*i);
+      }
+    delete discrProducts;
+  }
+
+  if (_verbose > 0) {
+    G4cout << "*-*-* Photon evaporation: " << products->size() << G4endl;
+  }
+  return products;
+}
 
 G4FragmentVector* G4PhotonEvaporation::BreakUp(const G4Fragment& nucleus)
 {
-  _nucleus = nucleus;
-
-  G4FragmentVector* products = new G4FragmentVector;
- 
-  _contDeexcitation->SetNucleus(nucleus);
-  _discrDeexcitation->SetNucleus(nucleus);
+  //G4cout << "G4PhotonEvaporation::BreakUp" << G4endl;
+  _nucleus = new G4Fragment(nucleus);
+
+  _contDeexcitation->SetNucleus(_nucleus);
+  _discrDeexcitation->SetNucleus(_nucleus);
   
   // Do one photon emission
@@ -104,130 +195,94 @@
   // Products from continuum gamma transitions
 
-  G4FragmentVector* contProducts = _contDeexcitation->DoTransition();  
-
-  G4int nCont = 0;
-  if (contProducts != 0) nCont = contProducts->size();
-
-  G4FragmentVector::iterator i;
-  if (nCont > 0)
-    {
-      G4Fragment modifiedNucleus = _contDeexcitation->GetNucleus();
-      _discrDeexcitation->SetNucleus(modifiedNucleus);
-      for (i = contProducts->begin(); i != contProducts->end(); i++)
-	{
-	  products->push_back(*i);
-	}
-    }
-  else
+  G4FragmentVector* products = _contDeexcitation->DoTransition();  
+  if( !products ) { products = new G4FragmentVector(); }
+  else if(_verbose > 0) {
+    G4cout << "G4PhotonEvaporation::BreakUp " << products->size() 
+	   << " gammas from ContinuesDeexcitation " << G4endl;
+    G4cout << "   Residual: " << nucleus << G4endl;
+  }
+
+  if (0 == products->size())
     {
       // Products from discrete gamma transitions
       G4FragmentVector* discrProducts = _discrDeexcitation->DoTransition();
-      
-      G4int nDiscr = 0;
-      if (discrProducts != 0) nDiscr = discrProducts->size();
-      
-      if (_verbose > 0)
-	G4cout << " = BreakUp = " << nDiscr 
-	       << " gammas from DiscreteDeexcitation " 
-	       << G4endl;
-      
-      for (i = discrProducts->begin(); i != discrProducts->end(); i++)
-	{
-	  products->push_back(*i);
+
+      if (discrProducts) {
+	_eOccupancy = _discrDeexcitation->GetEO();
+	_vShellNumber = _discrDeexcitation->GetVacantSN();
+
+	// not sure if the following line is needed!
+	_discrDeexcitation->SetVaccantSN(-1);
+	//
+	if (_verbose > 0) {
+	  G4cout << " = BreakUp = " << discrProducts->size() 
+		 << " gammas from DiscreteDeexcitation " 
+		 << G4endl;
+	  G4cout << "   Residual: " << nucleus << G4endl;
 	}
-      discrProducts->clear();
-      delete discrProducts;
+	G4FragmentVector::iterator i;
+	for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
+	  {
+	    products->push_back(*i);
+	  }
+	delete discrProducts;
+      }
     }
-
-
-  _gammaE = 0.;
-  if (products->size() > 0)
-    {
-      _gammaE = (*(products->begin()))->GetMomentum().e();
-    }
-
-  contProducts->clear();
-  delete contProducts;  // delete vector, not fragments 
-
+  
   // Add deexcited nucleus to products
-  G4Fragment* finalNucleus = new G4Fragment(_discrDeexcitation->GetNucleus());
-  products->push_back(finalNucleus);
-
-
-  if (_verbose > 0)
+  products->push_back(_nucleus);
+
+  if (_verbose > 0) {
     G4cout << "*-*-*-* Photon evaporation: " << products->size() << G4endl;
+  }
 
   return products;
 }
 
-
 G4FragmentVector* G4PhotonEvaporation::BreakItUp(const G4Fragment& nucleus)
 {
-  _nucleus = nucleus;
-
-  G4FragmentVector* products = new G4FragmentVector;
-
-  _contDeexcitation->SetNucleus(nucleus);
-  _discrDeexcitation->SetNucleus(nucleus);
+  // The same pointer of primary nucleus
+  _nucleus = new G4Fragment(nucleus);
+  _contDeexcitation->SetNucleus(_nucleus);
+  _discrDeexcitation->SetNucleus(_nucleus);
+
+  //G4cout << "G4PhotonEvaporation::BreakItUp:  " << nucleus << G4endl;
 
   // Do the whole gamma chain 
-
-  G4FragmentVector* contProducts = _contDeexcitation->DoChain();  
+  G4FragmentVector* products = _contDeexcitation->DoChain();  
+  if( !products ) { products = new G4FragmentVector; }
 
   // Products from continuum gamma transitions
-  G4int nCont = 0;
-  if (contProducts != 0) nCont = contProducts->size();
-
-  if (_verbose > 0)
-    G4cout << " = BreakItUp = " << nCont 
+  if (_verbose > 0) {
+    G4cout << " = BreakItUp = " << products->size()
 	   << " gammas from ContinuumDeexcitation " << G4endl;
-
-  G4FragmentVector::iterator i;
-  if (nCont > 0)
-    {
-      G4Fragment modifiedNucleus = _contDeexcitation->GetNucleus();
-      _discrDeexcitation->SetNucleus(modifiedNucleus);
-      for (i = contProducts->begin(); i != contProducts->end(); i++)
-	{
-	  products->push_back(*i);
-	}
-    }
+  }
 
   // Products from discrete gamma transitions
   G4FragmentVector* discrProducts = _discrDeexcitation->DoChain();
-  _eOccupancy = _discrDeexcitation->GetEO();
-  _vShellNumber = _discrDeexcitation->GetVacantSN();
-
-  // not sure if the following line is needed!
-  _discrDeexcitation->SetVaccantSN(-1);
-
-  G4int nDiscr = 0;
-  if (discrProducts != 0) nDiscr = discrProducts->size();
-
-  if (_verbose > 0)
-    G4cout << " = BreakItUp = " << nDiscr 
-	   << " gammas from DiscreteDeexcitation " << G4endl;
-
-  for (i = discrProducts->begin(); i != discrProducts->end(); i++)
-    {
-      products->push_back(*i);
+  if(discrProducts) {
+    _eOccupancy = _discrDeexcitation->GetEO();
+    _vShellNumber = _discrDeexcitation->GetVacantSN();
+
+    // not sure if the following line is needed!
+    _discrDeexcitation->SetVaccantSN(-1);
+
+    if (_verbose > 0) {
+      G4cout << " = BreakItUp = " << discrProducts->size() 
+	     << " gammas from DiscreteDeexcitation " << G4endl;
     }
-
+    G4FragmentVector::iterator i;
+    for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
+      {
+	products->push_back(*i);
+      }
+    delete discrProducts;
+  }
   // Add deexcited nucleus to products
-  G4Fragment* finalNucleus = new G4Fragment(_discrDeexcitation->GetNucleus());
-  products->push_back(finalNucleus);
-  if (_verbose > 0)
-    G4cout << " = BreakItUp = Nucleus added to products" << G4endl;
-
-  if (_verbose > 0)
+  products->push_back(_nucleus);
+
+  if (_verbose > 0) {
     G4cout << "*-*-* Photon evaporation: " << products->size() << G4endl;
-
-#ifdef debug
-  CheckConservation(nucleus,products);
-#endif
-  contProducts->clear();
-  discrProducts->clear();
-  delete contProducts;  // delete vector, not fragments 
-  delete discrProducts;
+  }
   return products;
 }
@@ -236,5 +291,7 @@
 {
   G4double prob = 0.;
-  if (_probAlgorithm != 0) prob = _probAlgorithm->EmissionProbability(_nucleus,_gammaE);
+  if (_probAlgorithm != 0) {
+    prob = _probAlgorithm->EmissionProbability(*_nucleus,_gammaE);
+  }
   return prob;
 }
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VGammaDeexcitation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VGammaDeexcitation.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VGammaDeexcitation.cc	(revision 1315)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4VGammaDeexcitation.cc,v 1.17 2010/05/09 17:31:23 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -45,4 +47,9 @@
 //              Added creation time evaluation for products of evaporation
 //
+//        19 April 2010, J. M. Quesada calculations in CM system
+//              pending final boost to lab system  (not critical)
+//
+//        23 April 2010, V.Ivanchenko rewite kinematic part using PDG formula
+//                                    for 2-body decay
 //
 // -------------------------------------------------------------------
@@ -64,38 +71,29 @@
 #include "G4DiscreteGammaTransition.hh"
 
+
 G4VGammaDeexcitation::G4VGammaDeexcitation(): _transition(0), _verbose(0),
 					      _electronO (0), _vSN(-1)
 { }
 
-
 G4VGammaDeexcitation::~G4VGammaDeexcitation()
 { 
-  //  if (_transition != 0) delete _transition;
-}
-
+  if (_transition != 0) { delete _transition; }
+}
 
 G4FragmentVector* G4VGammaDeexcitation::DoTransition()
 {
-  // Template method
-
- Initialize();
- G4FragmentVector* products = new G4FragmentVector;
-  
+  Initialize();
+  G4FragmentVector* products = new G4FragmentVector();
+ 
   if (CanDoTransition())
     {
-     G4Fragment* gamma = GenerateGamma();
-     if (gamma != 0)
-       {
-	 products->push_back(gamma);
-	 UpdateNucleus(gamma);
-	 Update();
-       }
-    }
-
-  if (_verbose > 1)
+      G4Fragment* gamma = GenerateGamma();
+      if (gamma != 0) { products->push_back(gamma); }
+    }
+ 
+  if (_verbose > 1) {
     G4cout << "G4VGammaDeexcitation::DoTransition - Transition deleted " << G4endl;
-
-  if (_transition != 0) delete _transition;
-
+  }
+ 
   return products;
 }
@@ -103,162 +101,95 @@
 G4FragmentVector* G4VGammaDeexcitation::DoChain()
 {
+  if (_verbose > 1) { G4cout << "G4VGammaDeexcitation::DoChain" << G4endl; }
+  const G4double tolerance = CLHEP::keV;
+
   Initialize();
-  G4FragmentVector* products = new G4FragmentVector;
-
+  G4FragmentVector* products = new G4FragmentVector();
+  
   while (CanDoTransition())
-    {
-      if (_verbose > 5) G4cout << "G4VGammaDeexcitation::DoChain -  Looping" << G4endl;
-
+    {      
+      _transition->SetEnergyFrom(_nucleus->GetExcitationEnergy());
       G4Fragment* gamma = GenerateGamma();
       if (gamma != 0) 
 	{
 	  products->push_back(gamma);
-	  UpdateNucleus(gamma);
-	  UpdateElectrons ();
+	  //G4cout << "Eex(keV)= " << _nucleus->GetExcitationEnergy()/keV << G4endl;
+	  if(_nucleus->GetExcitationEnergy() <= tolerance) { break; }
+	  Update();
 	}
-     Update();
     } 
-
-  if (_verbose > 1)
-      G4cout << "G4VGammaDeexcitation::DoChain - Transition deleted, end of chain " << G4endl;
-
-  if (_transition != 0) delete _transition;
+  
+  if (_verbose > 1) {
+    G4cout << "G4VGammaDeexcitation::DoChain - Transition deleted, end of chain " << G4endl;
+  }
   
   return products;
 }
 
-
-const G4Fragment& G4VGammaDeexcitation::GetNucleus() const
-{
-  return _nucleus; 
-}
-
-
-void G4VGammaDeexcitation::SetNucleus(const G4Fragment& nucleus)
-{
-  _nucleus = G4Fragment(nucleus);
-}
-
-
 G4Fragment* G4VGammaDeexcitation::GenerateGamma()
 {
+  // 23/04/10 V.Ivanchenko rewrite complitely
   G4double eGamma = 0.;
-
+  
   if (_transition != 0) {
     _transition->SelectGamma();  // it can be conversion electron too
     eGamma = _transition->GetGammaEnergy(); 
-  }
+    if(eGamma <= 0.0) { return 0; }
+  }
+  G4double excitation = _nucleus->GetExcitationEnergy() - eGamma;
+  if(excitation < 0.0) { excitation = 0.0; } 
   if (_verbose > 1 && _transition != 0 ) 
     {
-      G4cout << "G4VGammaDeexcitation::GenerateGamma - Gamma energy " << eGamma 
-	     << " ** New excitation " << _nucleus.GetExcitationEnergy() - eGamma
+      G4cout << "G4VGammaDeexcitation::GenerateGamma - Edeexc(MeV)= " << eGamma 
+	     << " ** left Eexc(MeV)= " << excitation
 	     << G4endl;
     }
   
-  // Photon momentum isotropically generated
-  // the same for electron 
-  
-  if (eGamma > 0.) 
-    {
-      G4double cosTheta = 1. - 2. * G4UniformRand();
-      G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
-      G4double phi = twopi * G4UniformRand();
-      G4double pM = eGamma;
-      G4DiscreteGammaTransition* dtransition = 0; 
-      dtransition = dynamic_cast <G4DiscreteGammaTransition*> (_transition);
-      if ( dtransition && !( dtransition->IsAGamma()) ) {
-	G4double eMass = G4Electron::ElectronDefinition()->GetPDGMass();
-	pM =std::sqrt(eGamma*(eGamma + 2.0*eMass));
-	eGamma = eGamma + eMass;
-      }
-      G4ThreeVector pGamma( pM * sinTheta * std::cos(phi),
-			    pM * sinTheta * std::sin(phi),
-			    pM * cosTheta );
-      G4LorentzVector gamma(pGamma, eGamma);
-      //	gamma.boost(_nucleus.GetMomentum().boostVector() );
-      G4Fragment* gammaFragment ;
-      if ( dtransition && !(dtransition->IsAGamma()) ){
-	gammaFragment = new G4Fragment(gamma,G4Electron::ElectronDefinition());
-      } else {
-	gammaFragment = new G4Fragment(gamma,G4Gamma::GammaDefinition()); 
-      }
-      G4double gammaTime = _transition->GetGammaCreationTime();
-      gammaTime += _nucleus.GetCreationTime();
-      gammaFragment->SetCreationTime(gammaTime);
-      
-      if (_verbose > 1 && dtransition )
-	G4cout << "G4VGammaDeexcitation::GenerateGamma -  Gamma fragment generated:  " 
-	       << (dtransition->IsAGamma() ? " Gamma" : " Electron" ) << G4endl;
-      return gammaFragment;
-    }
-  else
-    {
-      return 0;
-    }
-}
-
-void G4VGammaDeexcitation::UpdateNucleus(const G4Fragment*  gamma)
-{
-  G4LorentzVector p4Gamma = gamma->GetMomentum();
-  G4ThreeVector pGamma(p4Gamma.vect());
-  
-  G4double eGamma = 0.;
-  if (_transition != 0)
-    eGamma = _transition->GetGammaEnergy();
+  // Do complete Lorentz computation 
+
+  G4LorentzVector lv = _nucleus->GetMomentum();
+  G4double Mass = _nucleus->GetGroundStateMass() + excitation;
+
+  // select secondary
+  G4ParticleDefinition* gamma = G4Gamma::Gamma();
+
   G4DiscreteGammaTransition* dtransition = 0; 
   dtransition = dynamic_cast <G4DiscreteGammaTransition*> (_transition);
-  if (dtransition && !(dtransition->IsAGamma()) )      
-    eGamma += dtransition->GetBondEnergy(); 
-
-  G4LorentzVector p4Nucleus(_nucleus.GetMomentum() );
-
-// New tetravector calculation:
-
-//  G4double Mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(_nucleus.GetZ(),_nucleus.GetA());
-  G4double m1 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(static_cast<G4int>(_nucleus.GetZ()),
-									       static_cast<G4int>(_nucleus.GetA()));
-  G4double m2 = _nucleus.GetZ() *  G4Proton::Proton()->GetPDGMass() + 
-    (_nucleus.GetA()- _nucleus.GetZ())*G4Neutron::Neutron()->GetPDGMass();
-
-  G4double Mass = std::min(m1,m2);
-
-
-  G4double newExcitation = p4Nucleus.mag() - Mass - eGamma;
-  if(newExcitation < 0)
-    newExcitation = 0;
-  
-  G4ThreeVector p3Residual(p4Nucleus.vect() - pGamma);
-  G4double newEnergy = std::sqrt(p3Residual * p3Residual +
-			    (Mass + newExcitation) * (Mass + newExcitation));
-  G4LorentzVector p4Residual(p3Residual, newEnergy);
-  
-  //  G4LorentzVector p4Residual(-pGamma, p4Nucleus.e() - eGamma);
-  //  p4Residual.boost( _nucleus.GetMomentum().boostVector() );
-  
-  // Update excited nucleus parameters
-
-  _nucleus.SetMomentum(p4Residual);
-  _nucleus.SetCreationTime(gamma->GetCreationTime());
-
-//  if (_transition != 0) 
-//  {
-//    G4double excitation =_transition->GetEnergyTo();
-//    if (excitation < 0.) excitation = 0.0;
-//    _nucleus.SetExcitationEnergy(excitation);
-//  }
-
-  return;
-}
-
-void G4VGammaDeexcitation::UpdateElectrons( )
-{
-  G4DiscreteGammaTransition* dtransition = 0; 
-  dtransition = dynamic_cast <G4DiscreteGammaTransition*> (_transition);
-  if (dtransition && !(dtransition->IsAGamma()) ) {    
-    
+  if ( dtransition && !( dtransition->IsAGamma()) ) { 
+    gamma = G4Electron::Electron(); 
     _vSN = dtransition->GetOrbitNumber();   
     _electronO.RemoveElectron(_vSN);
-  }
-  return;
+    lv += G4LorentzVector(0.0,0.0,0.0,CLHEP::electron_mass_c2 - dtransition->GetBondEnergy());
+  }
+
+  // check consistency  
+  G4double eMass = gamma->GetPDGMass();
+
+  G4double Ecm       = lv.mag();
+  G4ThreeVector bst  = lv.boostVector();
+
+  G4double GammaEnergy = 0.5*((Ecm - Mass)*(Ecm + Mass) + eMass*eMass)/Ecm;
+  if(GammaEnergy <= eMass) { return 0; }
+
+  G4double cosTheta = 1. - 2. * G4UniformRand(); 
+  G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
+  G4double phi = twopi * G4UniformRand();
+  G4double mom = sqrt((GammaEnergy - eMass)*(GammaEnergy + eMass));
+  G4LorentzVector Gamma4P(mom * sinTheta * std::cos(phi),
+			  mom * sinTheta * std::sin(phi),
+			  mom * cosTheta,
+			  GammaEnergy);
+  Gamma4P.boost(bst);  
+  G4Fragment * thePhoton = new G4Fragment(Gamma4P,gamma);
+
+  G4double gammaTime = _nucleus->GetCreationTime() + _transition->GetGammaCreationTime();
+  thePhoton->SetCreationTime(gammaTime);
+
+  lv -= Gamma4P;
+  _nucleus->SetMomentum(lv);
+  _nucleus->SetCreationTime(gammaTime);
+
+  //G4cout << "G4VGammaDeexcitation::GenerateGamma left nucleus: " << _nucleus << G4endl;
+  return thePhoton;
 }
 
@@ -269,13 +200,14 @@
       delete _transition;
       _transition = 0;
-      if (_verbose > 1)
+      if (_verbose > 1) {
 	G4cout << "G4VGammaDeexcitation::Update - Transition deleted " << G4endl;
-    }
-
+      }
+    }
+  
   _transition = CreateTransition();
   if (_transition != 0) 
     {
-      _transition->SetEnergyFrom(_nucleus.GetExcitationEnergy());
-      //      if ( _vSN != -1) (dynamic_cast <G4DiscreteGammaTransition*> (_transition))->SetICM(false);
+      _transition->SetEnergyFrom(_nucleus->GetExcitationEnergy());
+      // if ( _vSN != -1) (dynamic_cast <G4DiscreteGammaTransition*> (_transition))->SetICM(false);
       // the above line is commented out for bug fix #952. It was intruduced for reason that
       // the k-shell electron is most likely one to be kicked out and there is no time for 
@@ -283,28 +215,5 @@
       // reported in #952
     }
-
+  
   return;
 }
-
-
-void G4VGammaDeexcitation::Initialize()
-{
-  _transition = CreateTransition();
-  if (_transition != 0) {
-    _transition->SetEnergyFrom(_nucleus.GetExcitationEnergy());
-  }
-  return;
-}
-
-
-void G4VGammaDeexcitation::SetVerboseLevel(G4int verbose)
-{
-  _verbose = verbose;
-  return;
-}
-
-
-
-
-
-
Index: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VPhotonEvaporation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VPhotonEvaporation.cc	(revision 1228)
+++ trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VPhotonEvaporation.cc	(revision 1315)
@@ -39,4 +39,7 @@
 //      
 // -------------------------------------------------------------------
+//
+// $Id: G4VPhotonEvaporation.cc,v 1.3 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 
 
