Changeset 1315 for trunk/source/processes/hadronic/models/de_excitation

trunk/source/processes/hadronic/models/de_excitation/History

-                      r1228
+                      r1315
      * Please list in reverse chronological order (last date on top)
      ---------------------------------------------------------------
+June 2010 Vladimir Ivanchenko (hadr-deex-V09-03-18)
+----------------------------------------------------------
+- G4Evaporation - fixed problem of isotope production for high Z
+                  fragments introduced in previous tags
+May 2010 Vladimir Ivanchenko (hadr-deex-V09-03-17)
+----------------------------------------------------------
+- G4E1Probability - added numerical protection to avoid division by zero
+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)
+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
+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
+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
+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
+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
+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
+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
+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
+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
+April 2010 Vladimir Ivanchenko (hadr-deex-V09-03-06)
+----------------------------------------------------------
+- G4GEMProbability - (JQM + VI) fixed numerical problem (division by zero)
+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
+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
+March 2010 Vladimir Ivanchenko (hadr-deex-V09-03-03)
+----------------------------------------------------------
+- G4Evaporation - set as a default variant evaporation combined
+                  standard + GEM probabilities
+February 2009 Vladimir Ivanchenko (hadr-deex-V09-03-02)
+----------------------------------------------------------
+- G4ExcitationHandler - deactivate Multi-Fragmentation model
+February 2009 Vladimir Ivanchenko (hadr-deex-V09-03-01)
+----------------------------------------------------------
+- G4ExcitationHandler - activate FermiBreakUp and Multi-Fragmentation models
+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,
 December 2009 Gunter Folger  (hadr-deex-V09-02-23)

trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4Evaporation.hh

-                      r1228
+                      r1315
 //
 //
 // $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
 …
 //
 //
+// 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
 …
 #include "G4VEvaporationChannel.hh"
 #include "G4Fragment.hh"
+#include "G4UnstableFragmentBreakUp.hh"
 class G4VEvaporationFactory;
+class G4NistManager;
 //#define debug
 …
 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);
 …
 public:
   G4FragmentVector * BreakItUp(const G4Fragment &theNucleus);
 …
 #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;
-  };
 };

trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationFactory.hh

-                      r1228
+                      r1315
 //
 //
 // $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
 …
+{
 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;

trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4VEvaporation.hh

-                      r1228
+                      r1315
 //
 //
 // $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
 …
+{
 public:
   G4VEvaporation() {};
   virtual ~G4VEvaporation() {}; // *
+  G4VEvaporation();
+  virtual ~G4VEvaporation();
 private:
 …
 public:
   virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus) = 0;
+  virtual void Initialise();
   // for inverse cross section choice

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc

-                      r1055
+                      r1315
   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.;

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc

-                      r1055
+                      r1315
   //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;
 …
   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.;

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc

-                      r1228
+                      r1315
 //
 //
 // $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
 …
 // 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"
 …
 #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);
+    }
+}
 …
   if (theChannelFactory != 0) delete theChannelFactory;
   theChannelFactory = new G4EvaporationFactory();
   theChannels = theChannelFactory->GetChannel();
+  InitialiseEvaporation();
+}
 …
   if (theChannelFactory != 0) delete theChannelFactory;
   theChannelFactory = new G4EvaporationGEMFactory();
   theChannels = theChannelFactory->GetChannel();
+  InitialiseEvaporation();
+}
 …
   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
 …
     return theResult;
+}
+*/

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationDefaultGEMFactory.cc

-                      r1228
+                      r1315
 //
 //
 // $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
 …
   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
 …
   theChannel->push_back( new G4Mg28GEMChannel() );     // Mg28
-  theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
-  theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
   return theChannel;

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc

-                      r1228
+                      r1315
 //
 //
 // $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
 …
 #include "G4EvaporationFactory.hh"
 #include "G4NeutronEvaporationChannel.hh"
 …
 #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*> *
 …
   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
 …
   theChannel->push_back( new G4AlphaEvaporationChannel() );    // Alpha
-  theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
-  theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
   return theChannel;

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc

-                      r1055
+                      r1315
 #include "G4EvaporationProbability.hh"
 #include "G4PairingCorrection.hh"
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc

-                      r1055
+                      r1315
   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.;

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc

-                      r962
+                      r1315
 ///////////////////////////////////////////////////////////////////////////////////
 //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)
 …
   FragmentAthrd=std::pow(FragmentA,0.33333);
   U=fragment.GetExcitationEnergy();
   if (OPTxs==0) {std::ostringstream errOs;
 …
 //     ** 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;
 …
   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.;
 …
     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.;

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc

-                      r1055
+                      r1315
   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.;

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc

-                      r1228
+                      r1315
 // ********************************************************************
 //
+//
+// $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
 …
 //
+#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()
+{}

trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/include/G4FermiPhaseSpaceDecay.hh

-                      r819
+                      r1315
 #include "G4LorentzVector.hh"
 #include "G4ParticleMomentum.hh"
+#include "Randomize.hh"
+#include "G4Pow.hh"
 #include <CLHEP/Random/RandGamma.h>
 …
+{
 public:
+  inline G4FermiPhaseSpaceDecay();
+  inline ~G4FermiPhaseSpaceDecay();
+  G4FermiPhaseSpaceDecay();
+  ~G4FermiPhaseSpaceDecay();
   inline std::vector<G4LorentzVector*> *
 …
 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;
 …
 };
+#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

trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiConfiguration.cc

-                      r1228
+                      r1315
 //
 //
 // $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"
 …
 // 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

trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiFragmentsPool.cc

-                      r1196
+                      r1315
 // 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"
 …
   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;
+    }

trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/src/G4FermiPhaseSpaceDecay.cc

-                      r819
+                      r1315
 #include <functional>
+G4FermiPhaseSpaceDecay::G4FermiPhaseSpaceDecay()
+{
+}
+G4FermiPhaseSpaceDecay::~G4FermiPhaseSpaceDecay()
+{
+}
 std::vector<G4LorentzVector*> *

trunk/source/processes/hadronic/models/de_excitation/fission/include/G4FissionBarrier.hh

-                      r1228
+                      r1315
 //
 //
 // $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

trunk/source/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc

-                      r1228
+                      r1315
 //
 //
 // $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
 …
 #include "G4CompetitiveFission.hh"
 #include "G4PairingCorrection.hh"
+#include "G4ParticleMomentum.hh"
 G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission")

trunk/source/processes/hadronic/models/de_excitation/fission/src/G4FissionBarrier.cc

-                      r1228
+                      r1315
 //
 //
 // $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

trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4EvaporationGEMFactory.hh

-                      r1228
+                      r1315
 //
 //
 // $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
 …
+{
 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;

trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4GEMProbability.hh

-                      r819
+                      r1315
 // ********************************************************************
 //
+// $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
 //
 //
 …
 // 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
 …
 #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

trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4EvaporationGEMFactory.cc

-                      r1228
+                      r1315
 //
 //
 // $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
 …
 #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()
 …
     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
 …
   theChannel->push_back( new G4Mg28GEMChannel() );     // Mg28
-  theChannel->push_back( new G4CompetitiveFission() ); // Fission Channel
-  theChannel->push_back( new G4PhotonEvaporation() );  // Photon Channel
   return theChannel;
+}

trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMProbability.cc

-                      r1196
+                      r1315
 // ********************************************************************
 //
+// $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)
 //
 //
 …
 //       -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;
+}
 …
     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
 …
         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;
+          }
 …
   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,
 …
   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)

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

-                      r1228
+                      r1315
 // ********************************************************************
 //
+//
+// $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
 …
 #include "G4VEvaporation.hh"
 #include "G4VPhotonEvaporation.hh"
+#include "G4VEvaporationChannel.hh"
 #include "G4Fragment.hh"
 #include "G4DynamicParticle.hh"
 …
 #include "G4IonConstructor.hh"
-//#define debug
 class G4ExcitationHandler
+{
 …
   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;
 …
   void SetFermiModel(G4VFermiBreakUp *const  value);
   void SetPhotonEvaporation(G4VPhotonEvaporation * const value);
+  void SetPhotonEvaporation(G4VEvaporationChannel * const value);
   void SetMaxZForFermiBreakUp(G4int aZ);
 …
   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;
 …
   const G4VFermiBreakUp * GetFermiModel() const;
   const G4VPhotonEvaporation * GetPhotonEvaporation() const;
+  const G4VEvaporationChannel * GetPhotonEvaporation() const;
   G4int GetMaxZ() const;
 …
                          G4FragmentVector * Result) const;
 #endif
 private:
 …
   G4VFermiBreakUp *theFermiModel;
   G4VPhotonEvaporation * thePhotonEvaporation;
+  G4VEvaporationChannel * thePhotonEvaporation;
   G4int maxZForFermiBreakUp;
   G4int maxAForFermiBreakUp;
   G4double minEForMultiFrag;
+  G4double minExcitation;
   G4ParticleTable *theTableOfParticles;
 …
   G4int OPTxs;
   G4bool useSICB;
   struct DeleteFragment
+  {
 …
+    }
   };
 };
+inline void G4ExcitationHandler::SetOPTxs(G4int opt)
+{
+  OPTxs = opt;
+  SetParameters();
+}
+inline void G4ExcitationHandler::UseSICB()
+{
+  useSICB = true;
+  SetParameters();
+}
 inline const G4VEvaporation * G4ExcitationHandler::GetEvaporation() const
 …
   MyOwnEvaporationClass = false;
   theEvaporation = value;
+  SetParameters();
+}
 …
 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;

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

-                      r1228
+                      r1315
 // ********************************************************************
 //
+// $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
 …
 //
 //
+// 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
 …
 //              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
 //
 …
 #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)
 …
   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()
 …
+}
+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
 …
       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
+                }
+            }
+        }
 …
   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
 …
   // -----------------------
+  // 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
 …
       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();

trunk/source/processes/hadronic/models/de_excitation/management/include/G4VEvaporationChannel.hh

-                      r1228
+                      r1315
 // ********************************************************************
 //
+//
+// $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
 …
+{
 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;}
 …
   // for superimposed Coulomb Barrier for inverse cross sections
   inline void UseSICB(G4bool use) { useSICB = use; }
 protected:
   G4int OPTxs;
   G4bool useSICB;
 private:

trunk/source/processes/hadronic/models/de_excitation/management/include/G4VEvaporationFactory.hh

-                      r1228
+                      r1315
 //
 //
 // $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
 …
+{
 public:
   G4VEvaporationFactory() : _channel(0) {};
+  G4VEvaporationFactory();
   virtual ~G4VEvaporationFactory();

trunk/source/processes/hadronic/models/de_excitation/management/src/G4VEvaporationChannel.cc

-                      r1228
+                      r1315
 // ********************************************************************
 //
+//
+// $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;
+}

trunk/source/processes/hadronic/models/de_excitation/management/src/G4VEvaporationFactory.cc

-                      r1228
+                      r1315
 //
 //
 // $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
 …
 #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()

trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFFragment.hh

-                      r1228
+                      r1315
 //
 //
 // $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
 …
 #include "G4StatMFParameters.hh"
 #include "G4ThreeVector.hh"
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
 #include "G4Fragment.hh"

trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMFMacroTemperature.hh

-                      r1228
+                      r1315
 //
 //
 // $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
 …
                              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)
 …
 private:
     // Default constructor
     G4StatMFMacroTemperature() {};
+    G4StatMFMacroTemperature();
     // copy constructor

trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroTemperature.cc

-                      r1228
+                      r1315
 //
 //
 // $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
 …
 //          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
 …
     G4int iterations = 0;
     while (fTa < 0.0 && iterations++ < 10) {
+    while (fTa < 0.0 && ++iterations < 10) {
         Ta -= 0.5*Ta;
         fTa = this->operator()(Ta);
 …
     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);
+    }
 …
       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.");
+    }
 …
     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();
 …
     // 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;
+}

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4ContinuumGammaDeexcitation.hh

-                      r819
+                      r1315
 // ********************************************************************
 //
+// $Id: G4ContinuumGammaDeexcitation.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
   // Destructor
   ~G4ContinuumGammaDeexcitation();
+  virtual ~G4ContinuumGammaDeexcitation();
   // Functions
 …
   virtual G4VGammaTransition* CreateTransition();
   virtual G4bool CanDoTransition() const;
+  virtual G4bool CanDoTransition();
 private:

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaDeexcitation.hh

-                      r819
+                      r1315
 // ********************************************************************
 //
+// $Id: G4DiscreteGammaDeexcitation.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
 // -------------------------------------------------------------------
 //
 #ifndef G4DiscreteGammaDeexcitation_hh
 #define G4DiscreteGammaDeexcitation_hh
 …
   // Destructor
   ~G4DiscreteGammaDeexcitation();
+  virtual ~G4DiscreteGammaDeexcitation();
   // Functions
 …
   virtual G4VGammaTransition * CreateTransition();
   virtual G4bool CanDoTransition() const;
+  virtual G4bool CanDoTransition();
   void SetICM(G4bool hl) { _icm = hl; };

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaTransition.hh

-                      r819
+                      r1315
 // ********************************************************************
 //
+// $Id: G4DiscreteGammaTransition.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
 //        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)
+//
 // -------------------------------------------------------------------
 …
 #include "G4NuclearLevel.hh"
+//JMQ 180410
+class G4NuclearLevelManager;
 class G4DiscreteGammaTransition : public G4VGammaTransition
+{
 …
   // 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
 …
   G4double _excitation;
   G4double _gammaCreationTime;
+  //JMQ 180410
+  G4int _A;
+  G4int _Z;
+  G4NuclearLevelManager * _levelManager;
+  //JMQ 190410
+  G4double _tolerance;
 };

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4E1Probability.hh

-                      r819
+                      r1315
 // ********************************************************************
 //
+// $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
 //
 //
 …
 #include "G4VEmissionProbability.hh"
 #include "G4Fragment.hh"
+#include "G4VLevelDensityParameter.hh"
+class G4Pow;
 class G4E1Probability : public G4VEmissionProbability
 …
 public:
   G4E1Probability() {};
+  G4E1Probability();
   ~G4E1Probability();
+  virtual ~G4E1Probability();
   G4double EmissionProbability(const G4Fragment& frag, const G4double excite);
 …
 private:
-  // G4E1Probability() {};
   G4E1Probability(const G4E1Probability& right);
 …
   // 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;
 };

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4PhotonEvaporation.hh

-                      r962
+                      r1315
 // ********************************************************************
 //
+// $Id: G4PhotonEvaporation.hh,v 1.7 2010/05/11 11:22:14 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
 //      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
 …
 //            G4ElectronOccupancy _eOccupancy;
 //            G4int _vShellNumber;
+//
+// 11 May 2010, V.Ivanchenko added EmittedFragment and BreakUpFragment
+//                           methods
 //
 // -------------------------------------------------------------------
 …
 #include "globals.hh"
-#include "G4VPhotonEvaporation.hh"
 #include "G4VEvaporationChannel.hh"
 #include "G4VEmissionProbability.hh"
 …
 #include "G4ElectronOccupancy.hh"
-//#define debug
 class G4Fragment;
 class G4PhotonEvaporation : public G4VPhotonEvaporation, public G4VEvaporationChannel {
+class G4PhotonEvaporation : public G4VEvaporationChannel {
 public:
 …
     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);
 …
     void SetEOccupancy( G4ElectronOccupancy  eOccupancy) ;
     G4ElectronOccupancy GetEOccupancy () { return _eOccupancy;} ;
 …
     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

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VGammaDeexcitation.hh

-                      r819
+                      r1315
 // ********************************************************************
 //
+// $Id: G4VGammaDeexcitation.hh,v 1.8 2010/04/28 14:22:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
 //
 //      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
 …
 //          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
 …
 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

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VPhotonEvaporation.hh

r819	r1315
24	24	// ********************************************************************
25	25	//
	26	// $Id: G4VPhotonEvaporation.hh,v 1.3 2010/04/25 18:43:21 vnivanch Exp $
	27	// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
26	28	//
27	29	// -------------------------------------------------------------------

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4ContinuumGammaDeexcitation.cc

-                      r819
+                      r1315
 // ********************************************************************
 //
+// $Id: G4ContinuumGammaDeexcitation.cc,v 1.7 2010/04/30 16:08:03 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
 //
 //      02 May 2003,   Vladimir Ivanchenko change interface to G4NuclearlevelManager
+//
+//      19 April 2010 J. M. Quesada: smaller value of tolerance parameter
 //
 // -------------------------------------------------------------------
 …
 //
+G4ContinuumGammaDeexcitation::G4ContinuumGammaDeexcitation(): _nucleusZ(0), _nucleusA(0), _levelManager(0)
+G4ContinuumGammaDeexcitation::G4ContinuumGammaDeexcitation()
+  : _nucleusZ(0), _nucleusA(0), _levelManager(0)
 { }
 …
 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)
 …
+    }
   if (_verbose > 1)
+  if (_verbose > 1) {
     G4cout << "G4ContinuumGammaDeexcitation::CreateTransition - Created" << G4endl;
+  }
   G4VGammaTransition* gt =  new G4ContinuumGammaTransition(_levelManager,Z,A,excitation,_verbose );
 …
 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;
+}

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaDeexcitation.cc

-                      r962
+                      r1315
 // ********************************************************************
 //
+// $Id: G4DiscreteGammaDeexcitation.cc,v 1.15 2010/05/10 07:20:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
   _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)
+    {
 …
       _levelManager = G4NuclearLevelStore::GetInstance()->GetManager(Z,A);
+    }
   if (_levelManager->IsValid())
 …
+        }
+      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;
 …
       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()
+{
 …
         << 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;
 …
+  }
+  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) {

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaTransition.cc

-                      r819
+                      r1315
 //
 //      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
 //
 // -------------------------------------------------------------------
 …
 #include "G4RandGeneralTmp.hh"
 #include "G4AtomicShells.hh"
+//JMQ:
+#include "G4NuclearLevelStore.hh"
+#include "G4Pow.hh"
 G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level):
 …
 { }
+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;
+}
 …
 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
 …
       //          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) {
 …
+              }
+            }
+            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;
 …
       //  G4RandGeneralTmp randGeneral(sampleArray, nBins);
       //G4double random = randGeneral.shoot();
       //_gammaCreationTime = tMin + (tMax - tMin) * random;
       // new code by Fan Lei
       //
+      /*
       if (tau != 0 )
+      {
 …
           //       << _gammaCreationTime/second << G4endl;
        } else { _gammaCreationTime=0.; }
+      */
+    }
   return;
+}
-//G4bool G4DiscreteGammaTransition::IsAGamma()
-//{
-//  return _aGamma;
-//}
 G4double G4DiscreteGammaTransition::GetGammaEnergy()
+{
 …
+{
   _excitation = energy;
+  return;
+}
+}

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4E1Probability.cc

-                      r819
+                      r1315
 //
 //
+//  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
 …
 G4double G4E1Probability::EmissionProbDensity(const G4Fragment& frag,
                                                  const G4double gammaE)
+                                              const G4double gammaE)
+{
 …
   // 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.
 …
   // 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
 …
   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;
 …
   //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
 …
+}
+G4E1Probability::~G4E1Probability() {}

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PhotonEvaporation.cc

-                      r962
+                      r1315
 // ********************************************************************
 //
+// $Id: G4PhotonEvaporation.cc,v 1.13 2010/05/17 11:47:43 flei Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
 //      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
 …
 //       G4ElectronOccupancy GetEOccupancy () ;
 //
+// 11 May 2010, V.Ivanchenko add implementation of EmittedFragment and
+//                           BreakUpFragment methods; cleanup logic
+//
 // -------------------------------------------------------------------
+//
 #include "G4PhotonEvaporation.hh"
 …
 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
 …
   // 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;
+}
 …
+{
   G4double prob = 0.;
+  if (_probAlgorithm != 0) prob = _probAlgorithm->EmissionProbability(_nucleus,_gammaE);
+  if (_probAlgorithm != 0) {
+    prob = _probAlgorithm->EmissionProbability(*_nucleus,_gammaE);
+  }
   return prob;
+}

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VGammaDeexcitation.cc

-                      r819
+                      r1315
 // ********************************************************************
 //
+// $Id: G4VGammaDeexcitation.cc,v 1.17 2010/05/09 17:31:23 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
 …
 //              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
 //
 // -------------------------------------------------------------------
 …
 #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;
+}
 …
 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;
+}
 …
       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
 …
       // 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;
+}

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VPhotonEvaporation.cc

-                      r819
+                      r1315
 //
 // -------------------------------------------------------------------
+//
+// $Id: G4VPhotonEvaporation.cc,v 1.3 2010/04/25 18:43:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $

Context Navigation

Legend:

Download in other formats: