Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4GNASHTransitions.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4GNASHTransitions.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4GNASHTransitions.hh	(revision 962)
@@ -52,4 +52,17 @@
   virtual G4Fragment PerformTransition(const G4Fragment & aFragment);
 
+//JMQ 03/01/08
+ public:
+// inline G4double GetTransitionProb1() const
+G4double GetTransitionProb1() 
+{return 0;}
+// inline G4double GetTransitionProb2() const
+G4double GetTransitionProb2() 
+{return 0;}
+// inline G4double GetTransitionProb3() const
+G4double GetTransitionProb3() 
+{return 0;}
+//
+
 
 };
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4LowEIonFragmentation.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4LowEIonFragmentation.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4LowEIonFragmentation.hh	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4LowEIonFragmentation.hh,v 1.3 2006/06/29 20:58:04 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by H.P. Wellisch
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundAlpha.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundAlpha.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundAlpha.hh	(revision 962)
@@ -24,9 +24,8 @@
 // ********************************************************************
 //
+// by V. Lara
 //
-// $Id: G4PreCompoundAlpha.hh,v 1.6 2007/10/01 10:41:59 ahoward Exp $
-// GEANT4 tag $Name:  $
-//
-// by V. Lara
+//J. M. Quesada (July 08) 
+
 
 #ifndef G4PreCompoundAlpha_h
@@ -36,7 +35,6 @@
 #include "G4ReactionProduct.hh"
 #include "G4Alpha.hh"
-
 #include "G4AlphaCoulombBarrier.hh"
-
+#include "G4PreCompoundParameters.hh"
 
 class G4PreCompoundAlpha : public G4PreCompoundIon
@@ -47,5 +45,5 @@
 
   // copy constructor
-  G4PreCompoundAlpha(const G4PreCompoundAlpha &right): G4PreCompoundIon(right) {}
+  G4PreCompoundAlpha(const G4PreCompoundAlpha &right): G4PreCompoundIon(right){}
 
   // destructor
@@ -66,80 +64,44 @@
 
 
-  G4ReactionProduct * GetReactionProduct() const
-  {
-    G4ReactionProduct * theReactionProduct =
-      new G4ReactionProduct(G4Alpha::AlphaDefinition());
-    theReactionProduct->SetMomentum(GetMomentum().vect());
-    theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-    return theReactionProduct;
-  }   
-    
+  G4ReactionProduct * GetReactionProduct() const;
+
 private:
 
-// added Rj method according to literature and JMQ
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged)
-  {
-    G4double rj = 1.0;
-    G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2)*(NumberParticles-3);
-    if(denominator !=0) rj = 6.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); //JMQ 23/8/07
+  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged);
 
-    return rj;
-  }
+  virtual G4double CrossSection(const  G4double K) ; 
+
+  virtual G4double FactorialFactor(const G4double N, const G4double P);
+
+  virtual G4double CoalescenceFactor(const G4double A);
+
+  G4double GetOpt0(const G4double K);
+  G4double GetOpt12(const G4double K);
+  G4double GetOpt34(const G4double K);
+
+  G4double GetAlpha();
+  
+  G4double GetBeta();
+
+//data members
+
+      G4AlphaCoulombBarrier theAlphaCoulombBarrier;
+      G4double ResidualA;
+      G4double ResidualZ; 
+      G4double theA;
+      G4double theZ;
+      G4double ResidualAthrd;
+      G4double FragmentA;
+      G4double FragmentAthrd;
+
+};
+#endif
 
 
-  virtual G4double GetAlpha()
-  {
-    G4double C = 0.0;
-    G4double aZ = GetZ() + GetRestZ();
-    if (aZ <= 30) 
-      {
-	C = 0.10;
-      } 
-    else if (aZ <= 50) 
-      {
-	C = 0.1 + -((aZ-50.)/20.)*0.02;
-      } 
-    else if (aZ < 70) 
-      {
-	C = 0.08 + -((aZ-70.)/20.)*0.02;
-      }
-    else 
-      {
-	C = 0.06;
-      }
-    return 1.0+C;
-  }
 
-  virtual G4double GetBeta()
-  {
-    return -GetCoulombBarrier();
-  }
 
-  virtual G4double FactorialFactor(const G4double N, const G4double P)
-  {
-    return 
-      (N-4.0)*(P-3.0)*(
-		       (((N-3.0)*(P-2.0))/2.0) *(
-						 (((N-2.0)*(P-1.0))/3.0) *(
-									   (((N-1.0)*P)/2.0)
-									   )
-						 )
-		       );
-  }
 
-  virtual G4double CoalescenceFactor(const G4double A)
-  {
-    return 4096.0/(A*A*A);
-  }    
-private:
-
-  G4AlphaCoulombBarrier theAlphaCoulombBarrier;
-
-};
-
-#endif
  
 
+
+
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundDeuteron.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundDeuteron.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundDeuteron.hh	(revision 962)
@@ -24,9 +24,7 @@
 // ********************************************************************
 //
+// by V. Lara
 //
-// $Id: G4PreCompoundDeuteron.hh,v 1.7 2007/10/01 10:41:59 ahoward Exp $
-// GEANT4 tag $Name:  $
-//
-// by V. Lara
+//J. M. Quesada (July 08) 
 
 #ifndef G4PreCompoundDeuteron_h
@@ -36,6 +34,6 @@
 #include "G4ReactionProduct.hh"
 #include "G4Deuteron.hh"
-
 #include "G4DeuteronCoulombBarrier.hh"
+#include "G4PreCompoundParameters.hh"
 
 
@@ -66,67 +64,37 @@
 
 
-  G4ReactionProduct * GetReactionProduct() const
-  {
-    G4ReactionProduct * theReactionProduct =
-      new G4ReactionProduct(G4Deuteron::DeuteronDefinition());
-    theReactionProduct->SetMomentum(GetMomentum().vect());
-    theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-    return theReactionProduct;
-  }   
+  G4ReactionProduct * GetReactionProduct() const;
+
   
 private:
 
-// added Rj method according to literature and JMQ - formula from Jose Quesada
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged)
-  {
-    G4double rj = 1.0;
-    G4double denominator = NumberParticles*(NumberParticles-1);
-    if(denominator != 0) rj = 2.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged))/static_cast<G4double>(denominator); //JMQ re-corrected 23/8/07
-    //    return static_cast<G4double>(NumberCharged*(NumberCharged-1))/static_cast<G4double>(NumberParticles*(NumberParticles-1));
+  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged);
 
-    return rj;
-  }
+  virtual G4double CrossSection(const  G4double K) ; 
 
+  virtual G4double FactorialFactor(const G4double N, const G4double P);
 
+  virtual G4double CoalescenceFactor(const G4double A);
 
-  virtual G4double GetAlpha()
-  {
-    G4double C = 0.0;
-    G4double aZ = GetZ() + GetRestZ();
-    if (aZ >= 70) 
-      {
-	C = 0.10;
-      } 
-    else 
-      {
-	C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 
-      }
-    return 1.0 + C/2.0;
-  }
-    
-  virtual G4double GetBeta()
-  {
-    return -GetCoulombBarrier();
-  }
+  G4double GetOpt0(const G4double K);
+  G4double GetOpt12(const G4double K);
+  G4double GetOpt34(const G4double K);
 
-  virtual G4double FactorialFactor(const G4double N, const G4double P)
-  {
-    return (N-1.0)*(N-2.0)*(P-1.0)*P/2.0;
-  }
+  G4double GetAlpha();
   
-  virtual G4double CoalescenceFactor(const G4double A)
-  {
-    return 16.0/A;
-  }    
-private:
-  
-  G4DeuteronCoulombBarrier theDeuteronCoulombBarrier;
+  G4double GetBeta();
+
+//data members
+
+      G4DeuteronCoulombBarrier theDeuteronCoulombBarrier;
+      G4double ResidualA;
+      G4double ResidualZ; 
+      G4double theA;
+      G4double theZ;
+      G4double ResidualAthrd;
+      G4double FragmentA;
+      G4double FragmentAthrd;
 
 };
+#endif
 
-#endif
- 
-
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.hh	(revision 962)
@@ -24,9 +24,14 @@
 // ********************************************************************
 //
-// $Id: G4PreCompoundEmission.hh,v 1.3 2006/06/29 20:58:10 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4PreCompoundEmission.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear Preequilibrium
 // by V. Lara 
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
+// JMQ (06 September 2008) Also external choice has been added for:
+//                      - superimposed Coulomb barrier (if useSICB=true) 
 
 #ifndef G4PreCompoundEmission_h
@@ -78,4 +83,13 @@
 	       const G4double E, const G4double Ef) const;
 
+  G4double factorial(G4double a) const;
+
+  //for inverse cross section choice
+public:
+  inline void SetOPTxs(G4int);
+  //for superimposed CoulomBarrier for inverse cross sections
+  inline void UseSICB(G4bool);
+
+
   //==============
   // Data Members
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.icc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.icc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.icc	(revision 962)
@@ -24,4 +24,16 @@
 // ********************************************************************
 //
+// $Id: G4PreCompoundEmission.icc,v 1.5 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
+//
+//
+// Author:         V.Lara
+//
+// Modified:  
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
+// JMQ (06 September 2008) Also external choice has been added for:
+//                      - superimposed Coulomb barrier (if useSICB=true) 
+
 inline void G4PreCompoundEmission::Initialize(const G4Fragment & aFragment) 
 {
@@ -29,5 +41,4 @@
   return;
 }
-
 
 inline G4double G4PreCompoundEmission::GetTotalProbability(const G4Fragment & aFragment) 
@@ -46,2 +57,16 @@
   return;
 }
+
+//for inverse cross section choice
+inline void G4PreCompoundEmission::SetOPTxs(G4int opt)
+{
+  theFragmentsVector->SetOPTxs(opt);
+  return;
+}
+//for superimposed Coumlomb Barrier for inverse cross sections 
+inline void G4PreCompoundEmission::UseSICB(G4bool use)
+{
+  theFragmentsVector->UseSICB(use);
+  return;
+}
+
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragment.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragment.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragment.hh	(revision 962)
@@ -24,5 +24,12 @@
 // ********************************************************************
 //
-// by V. Lara
+//J. M. Quesada (August 2008).  
+//Based  on previous work by V. Lara
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option (default OPTxs=2)
+// JMQ (06 September 2008) Also external choice has been added for
+// superimposed Coulomb barrier (if useSICB=true, default false) 
+
 
 #ifndef G4PreCompoundFragment_h
@@ -62,5 +69,5 @@
   
   // Initialization method
-  void Initialize(const G4Fragment & aFragment);
+//  void Initialize(const G4Fragment & aFragment);
     
   // ================================================
@@ -84,5 +91,6 @@
   virtual G4double 
   ProbabilityDistributionFunction(const G4double K, 
-				  const G4Fragment & aFragment) = 0;  
+				  const G4Fragment & aFragment) = 0; 
+
 };
 
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.hh	(revision 962)
@@ -25,17 +25,19 @@
 //
 //
-// $Id: G4PreCompoundFragmentVector.hh,v 1.3 2006/06/29 20:58:18 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4PreCompoundFragmentVector.hh,v 1.6 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear Preequilibrium
-// by V. Lara 
+// by V. Lara
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option 
+// JMQ (06 September 2008) Also external choice has been added for:
+//                      - superimposed Coulomb barrier (if useSICB=true) 
 
 #ifndef G4PreCompoundFragmentVector_h
 #define G4PreCompoundFragmentVector_h 1
 
-
 #include "G4VPreCompoundFragment.hh"
-
-
 
 class G4PreCompoundFragmentVector 
@@ -69,4 +71,11 @@
   G4double TotalEmissionProbability;
 
+//for inverse cross section choice
+public:
+  inline void SetOPTxs(G4int);
+  //for superimposed CoulomBarrier for inverse cross sections
+  inline void UseSICB(G4bool);
+
+
 };
 
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.icc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.icc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundFragmentVector.icc	(revision 962)
@@ -25,10 +25,14 @@
 //
 //
-// $Id: G4PreCompoundFragmentVector.icc,v 1.3 2006/06/29 20:58:20 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4PreCompoundFragmentVector.icc,v 1.4 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear Preequilibrium
 // by V. Lara 
-
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
+// JMQ (06 September 2008) Also external choice has been added for:
+//                      - superimposed Coulomb barrier (if useSICB=true) 
 
 inline G4PreCompoundFragmentVector::G4PreCompoundFragmentVector(pcfvector * avector) :
@@ -64,2 +68,20 @@
   return;
 }
+
+//for inverse cross section choice
+inline void G4PreCompoundFragmentVector::
+SetOPTxs(G4int opt)
+{    
+  for (pcfvector::iterator i = theChannels->begin(); i != theChannels->end(); ++i) 
+    (*i)->SetOPTxs(opt);
+  return;
+}
+
+//for superimposed Coulomb Barrier for inverse  cross sections 
+inline void G4PreCompoundFragmentVector::
+UseSICB(G4bool use)
+{    
+  for (pcfvector::iterator i = theChannels->begin(); i != theChannels->end(); ++i) 
+    (*i)->UseSICB(use);
+  return;
+}
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundHe3.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundHe3.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundHe3.hh	(revision 962)
@@ -24,9 +24,7 @@
 // ********************************************************************
 //
+// by V. Lara
 //
-// $Id: G4PreCompoundHe3.hh,v 1.6 2007/10/01 10:42:00 ahoward Exp $
-// GEANT4 tag $Name:  $
-//
-// by V. Lara
+//J. M. Quesada (July 08) 
 
 #ifndef G4PreCompoundHe3_h
@@ -36,7 +34,6 @@
 #include "G4ReactionProduct.hh"
 #include "G4He3.hh"
-
 #include "G4He3CoulombBarrier.hh"
-
+#include "G4PreCompoundParameters.hh"
 
 class G4PreCompoundHe3 : public G4PreCompoundIon
@@ -66,79 +63,43 @@
 
 
-  G4ReactionProduct * GetReactionProduct() const
-  {
-    G4ReactionProduct * theReactionProduct =
-      new G4ReactionProduct(G4He3::He3Definition());
-    theReactionProduct->SetMomentum(GetMomentum().vect());
-    theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-    return theReactionProduct;
-  }   
-    
+  G4ReactionProduct * GetReactionProduct() const;
+ 
+
 private:
 
-// added Rj method according to literature and JMQ
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged)
-  {
-    G4double rj = 1.0;
-    G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2);
-    if(denominator != 0) rj = 3.0*static_cast<G4double>(NumberCharged*(NumberCharged-1)*(NumberParticles-NumberCharged))/static_cast<G4double>(denominator); //JMQ 23/8/07
+  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged);
 
-    return rj;
-  }
+  virtual G4double CrossSection(const  G4double K) ; 
+
+  virtual G4double FactorialFactor(const G4double N, const G4double P);
+
+  virtual G4double CoalescenceFactor(const G4double A);
+
+  G4double GetOpt0(const G4double K);
+  G4double GetOpt12(const G4double K);
+  G4double GetOpt34(const G4double K);
+
+  G4double GetAlpha();
+  
+  G4double GetBeta();
+
+//data members
+
+      G4He3CoulombBarrier theHe3CoulombBarrier;
+        G4double ResidualA;
+      G4double ResidualZ; 
+      G4double theA;
+      G4double theZ;
+      G4double ResidualAthrd;
+      G4double FragmentA;
+      G4double FragmentAthrd;
+
+
+};
+#endif
 
 
 
-  virtual G4double GetAlpha()
-  {
-    G4double C = 0.0;
-    G4double aZ = GetZ() + GetRestZ();
-    if (aZ <= 30) 
-      {
-	C = 0.10;
-      }
-    else if (aZ <= 50) 
-      {
-	C = 0.1 + -((aZ-50.)/20.)*0.02;
-      } 
-    else if (aZ < 70) 
-      {
-	C = 0.08 + -((aZ-70.)/20.)*0.02;
-      }
-    else 
-      {
-	C = 0.06;
-      }
-    return 1.0 + C*(4.0/3.0);
-  }
 
-  virtual G4double GetBeta()
-  {
-    return -GetCoulombBarrier();
-  }
 
-  virtual G4double FactorialFactor(const G4double N, const G4double P)
-  {
-    return 
-      (N-3.0)*(P-2.0)*(
-		       (((N-2.0)*(P-1.0))/2.0) *(
-						 (((N-1.0)*P)/3.0) 
-						 )
-		       );
-  }
-
-  virtual G4double CoalescenceFactor(const G4double A)
-  {
-    return 243.0/(A*A);
-  }    
-private:
-
-  G4He3CoulombBarrier theHe3CoulombBarrier;
-
-};
-
-#endif
  
-
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundIon.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundIon.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundIon.hh	(revision 962)
@@ -24,5 +24,7 @@
 // ********************************************************************
 //
-// by V. Lara
+//J. M. Quesada (August 2008).  
+//Based  on previous work by V. Lara
+//
 
 #ifndef G4PreCompoundIon_h
@@ -70,22 +72,22 @@
   }
     
-  virtual G4double ProbabilityDistributionFunction(const G4double eKin,
-						   const G4Fragment& aFragment);
+  virtual G4double ProbabilityDistributionFunction(const G4double eKin, 
+                                                   const G4Fragment& aFragment);
 
-protected:
-  G4bool IsItPossible(const G4Fragment& aFragment) 
-  {
-    G4int pplus = aFragment.GetNumberOfCharged();   
-    G4int pneut = aFragment.GetNumberOfParticles()-pplus;
-    return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
-  }
+  private:
+
+  G4bool IsItPossible(const G4Fragment& aFragment) ;
   
-  virtual G4double GetAlpha() = 0;
-  virtual G4double GetBeta() = 0;
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0;
+  protected:
+
+  virtual G4double CrossSection(const G4double ekin)=0; 
+
+  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0; 
+
   virtual G4double FactorialFactor(const G4double N, const G4double P) = 0;
+
   virtual G4double CoalescenceFactor(const G4double A) = 0; 
-    
-};
+
+   };
 
 #endif
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundModel.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundModel.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundModel.hh	(revision 962)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4PreCompoundModel.hh,v 1.3 2006/06/29 20:58:26 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4PreCompoundModel.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
@@ -36,4 +36,12 @@
 // transport, or any of the string-parton models.
 // Class Description - End
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option.(default OPTxs=3)
+// JMQ (06 September 2008) Also external choices have been added for:
+//                - superimposed Coulomb barrier (if useSICB=true, default false) 
+//                - "never go back"  hipothesis (if useNGB=true, default false) 
+//                - soft cutoff from preeq. to equlibrium (if useSCO=true, default false)
+//                - CEM transition probabilities (if useCEMtr=true, default false)  
 
 #ifndef G4PreCompoundModel_h
@@ -50,4 +58,11 @@
 #include "Randomize.hh"
 
+//#include "G4PreCompoundEmission.hh"
+
+#include "G4DynamicParticle.hh"
+#include "G4ReactionProductVector.hh"
+#include "G4ReactionProduct.hh"
+#include "G4ParticleTypes.hh"
+#include "G4ParticleTable.hh"
 
 //#define debug
@@ -56,8 +71,9 @@
 class G4PreCompoundModel : public G4VPreCompoundModel
 {
+ 
 public:
-    
   G4PreCompoundModel(G4ExcitationHandler * const value) : 
-    G4VPreCompoundModel(value), useHETCEmission(false), useGNASHTransition(false) {}
+    G4VPreCompoundModel(value), useHETCEmission(false), useGNASHTransition(false), 
+    OPTxs(3), useSICB(false), useNGB(false), useSCO(false), useCEMtr(false) {}
 
   ~G4PreCompoundModel() {}
@@ -89,8 +105,17 @@
   inline void UseDefaultTransition() { useGNASHTransition = false; }
 
+ //for cross section selection
+  inline void SetOPTxs(G4int opt) { OPTxs = opt; }
+//for the rest of external choices
+  inline void UseSICB() { useSICB = true; }
+  inline void UseNGB()  { useNGB = true; }
+  inline void UseSCO()  { useSCO = true; }
+  inline void UseCEMtr() { useCEMtr = true; }
 private:  
 
   void PerformEquilibriumEmission(const G4Fragment & aFragment, 
 				  G4ReactionProductVector * theResult) const;
+
+private:
 
 #ifdef debug				  
@@ -104,6 +129,18 @@
   //==============
 
+
+
   G4bool           useHETCEmission;
   G4bool           useGNASHTransition;
+
+//for cross section options
+  G4int OPTxs;
+//for the rest of external choices
+  G4bool useSICB;
+  G4bool useNGB;
+  G4bool useSCO;
+  G4bool useCEMtr;
+
+
     G4HadFinalState theResult;
 
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundNeutron.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundNeutron.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundNeutron.hh	(revision 962)
@@ -25,8 +25,5 @@
 //
 //
-// $Id: G4PreCompoundNeutron.hh,v 1.7 2007/10/01 10:42:00 ahoward Exp $
-// GEANT4 tag $Name:  $
-//
-// by V. Lara
+//J. M. Quesada (July  08) 
 
 
@@ -39,5 +36,4 @@
 #include "G4PreCompoundParameters.hh"
 #include "Randomize.hh"
-
 #include "G4NeutronCoulombBarrier.hh"
 
@@ -47,5 +43,5 @@
 public:
   // default constructor
-  G4PreCompoundNeutron() : G4PreCompoundNucleon(1,0,&theNeutronCoulomBarrier,"Neutron") {}
+  G4PreCompoundNeutron() : G4PreCompoundNucleon(1,0,&theNeutronCoulombBarrier,"Neutron") {}
 
   // copy constructor
@@ -68,49 +64,37 @@
 
 
-  G4ReactionProduct * GetReactionProduct() const
-  {
-    G4ReactionProduct * theReactionProduct = 
-      new G4ReactionProduct(G4Neutron::NeutronDefinition());
-    theReactionProduct->SetMomentum(GetMomentum().vect());
-    theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-    return theReactionProduct;
-  }
+  G4ReactionProduct * GetReactionProduct() const;
   
 private:
 
-// added Rj method according to literature and JMQ - formula from Jose Quesada
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged)
-  {
-    G4double rj = 1.0;
-    if(NumberParticles != 0) rj = static_cast<G4double>(NumberParticles - NumberCharged)/static_cast<G4double>(NumberParticles);
-    return rj;
-  }
+  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged); 
+
+  virtual G4double CrossSection(const  G4double K); 
 
 
-  virtual G4double GetAlpha()
-  {
-    return 0.76+2.2/std::pow(GetRestA(),1.0/3.0);
-  }
+  G4double GetOpt0(const G4double K);
+  G4double GetOpt12(const G4double K);
+  G4double GetOpt34(const G4double K);
+
   
-  virtual G4double GetBeta() 
-  {
-    return (2.12/std::pow(GetRestA(),2.0/3.0)-0.05)*MeV/GetAlpha();
-  }
+  G4double GetAlpha();
   
-  virtual G4bool IsItPossible(const G4Fragment& aFragment)
-  {
-    return ((aFragment.GetNumberOfParticles()-aFragment.GetNumberOfCharged()) >= 1);  
-  }
+  G4double GetBeta();
   
+//data members
 
-private:
-  
-  G4NeutronCoulombBarrier theNeutronCoulomBarrier;
+      G4NeutronCoulombBarrier theNeutronCoulombBarrier;
+      G4double ResidualA;
+      G4double ResidualZ; 
+      G4double theA;
+      G4double theZ;
+      G4double ResidualAthrd;
+      G4double FragmentA;
+      G4double FragmentAthrd;
 
+ 
 };
 
+ 
 #endif
 
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundNucleon.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundNucleon.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundNucleon.hh	(revision 962)
@@ -24,5 +24,8 @@
 // ********************************************************************
 //
-// by V. Lara
+//J. M. Quesada (August 2008).  
+//Based  on previous work by V. Lara
+//
+
 
 #ifndef G4PreCompoundNucleon_h
@@ -46,8 +49,9 @@
   G4PreCompoundNucleon(const G4double anA, 
 		       const G4double aZ, 
-		       G4VCoulombBarrier* aCoulombBarrier,
-		       const G4String & aName): 
-    G4PreCompoundFragment(anA,aZ,aCoulombBarrier,aName) {}
-  
+		       G4VCoulombBarrier* aCoulombBarrier,		       
+		       const G4String & aName) :
+    G4PreCompoundFragment(anA,aZ,aCoulombBarrier,aName) {} 
+
+
   virtual ~G4PreCompoundNucleon() {}
 
@@ -72,14 +76,16 @@
   virtual G4double ProbabilityDistributionFunction(const G4double eKin,
 						   const G4Fragment& aFragment);
+  
+  private:
 
+  G4bool IsItPossible(const G4Fragment&) ;   
     
-protected:
+ protected:
 
-// added Rj method according to literature and JMQ
+  virtual G4double CrossSection(const G4double ekin)=0; 
+
   virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0; 
-  virtual G4double GetAlpha() = 0;
-  virtual G4double GetBeta() = 0;
-  virtual G4bool IsItPossible(const G4Fragment&) = 0;     
-};
+
+ };
 
 #endif
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundParameters.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundParameters.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundParameters.hh	(revision 962)
@@ -25,8 +25,10 @@
 //
 //
-// $Id: G4PreCompoundParameters.hh,v 1.3 2006/06/29 20:58:32 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4PreCompoundParameters.hh,v 1.5 2008/05/08 10:34:25 quesada Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
+//
+//J. M. Quesada (Apr. 2008) Level density set to A/10 at preequilibrium
 
 
@@ -42,5 +44,7 @@
 
     // default constructor
-    G4PreCompoundParameters() : theLevelDensity(0.125/MeV),
+//    G4PreCompoundParameters() : theLevelDensity(0.125/MeV),
+//JMQ level density parameter  set to  A/10 at preequilibrium
+    G4PreCompoundParameters() : theLevelDensity(0.10/MeV),
       r0(1.5*fermi),Transitions_r0(0.6*fermi),FermiEnergy(35.0*MeV) 
 	{}
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundProton.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundProton.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundProton.hh	(revision 962)
@@ -25,8 +25,5 @@
 //
 //
-// $Id: G4PreCompoundProton.hh,v 1.7 2007/10/01 10:42:00 ahoward Exp $
-// GEANT4 tag $Name:  $
-//
-// by V. Lara
+//J. M. Quesada (July 08) 
 
 #ifndef G4PreCompoundProton_h
@@ -38,7 +35,5 @@
 #include "G4PreCompoundParameters.hh"
 #include "Randomize.hh"
-
 #include "G4ProtonCoulombBarrier.hh"
-
 
 class G4PreCompoundProton : public G4PreCompoundNucleon
@@ -68,59 +63,38 @@
 
 
-  G4ReactionProduct * GetReactionProduct() const
-  {
-    G4ReactionProduct * theReactionProduct = 
-      new G4ReactionProduct(G4Proton::ProtonDefinition());
-    theReactionProduct->SetMomentum(GetMomentum().vect());
-    theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-    return theReactionProduct;
-  }
+  G4ReactionProduct * GetReactionProduct() const;
+ 
   
 private:
 
-// added Rj method according to literature and JMQ - formula from Jose Quesada
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged)
-  {
-    G4double rj = 1.0;
-    if(NumberParticles != 0) rj = static_cast<G4double>(NumberCharged)/static_cast<G4double>(NumberParticles);
-    return rj;
-  }
+//JMQ (Sep. 07) combinatorial factor Rj
+  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged);
+
+  virtual G4double CrossSection(const  G4double K) ; 
 
 
-  virtual G4double GetAlpha()
-  {
-    G4double aZ = static_cast<G4double>(GetRestZ());
-    G4double C = 0.0;
-    if (aZ >= 70) 
-      {
-	C = 0.10;
-      } 
-    else 
-      {
-	C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
-      }
-    return 1.0 + C;
-  }
+  G4double GetOpt0(const G4double K);
+  G4double GetOpt1(const G4double K);
+  G4double GetOpt2(const G4double K);
+  G4double GetOpt3(const G4double K);
 
-  virtual G4double GetBeta()
-  {
-    return -GetCoulombBarrier();
-  }
+  G4double GetAlpha();
   
-  virtual G4bool IsItPossible(const G4Fragment& aFragment)
-  {
-    return (aFragment.GetNumberOfCharged() >= 1);
-  }
-  
-  
-private:
-  
-  G4ProtonCoulombBarrier theProtonCoulombBarrier;
-  
+  G4double GetBeta();
+
+//data members
+
+      G4ProtonCoulombBarrier theProtonCoulombBarrier;
+       G4double ResidualA;
+      G4double ResidualZ; 
+      G4double theA;
+      G4double theZ;
+      G4double ResidualAthrd;
+      G4double FragmentA;
+      G4double FragmentAthrd;
+
+
+ 
 };
-
 #endif
  
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTransitions.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTransitions.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTransitions.hh	(revision 962)
@@ -25,8 +25,11 @@
 //
 //
-// $Id: G4PreCompoundTransitions.hh,v 1.3 2006/06/29 20:58:36 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4PreCompoundTransitions.hh,v 1.6 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
+// J. M. Quesada . New methods for accessing to individual transition probabilities (landa+, landa-, landa0) from
+// G4PreCompoundModel.cc  
+
 
 #ifndef G4PreCompoundTransitions_h
@@ -50,4 +53,5 @@
 
 class G4PreCompoundTransitions : public G4VPreCompoundTransitions
+
 {
 public:
@@ -55,5 +59,5 @@
   // Calculates transition probabilities with Delta N = +2 (Trans1) -2 (Trans2) and 0 (Trans3)
 
-  G4PreCompoundTransitions() : TransitionProb1(0.0), TransitionProb2(0.0), TransitionProb3(0.0) {}
+  G4PreCompoundTransitions() : TransitionProb1(0.0), TransitionProb2(0.0), TransitionProb3(0.0){}
 
   virtual ~G4PreCompoundTransitions() {}
@@ -81,4 +85,23 @@
   G4double TransitionProb3;
 
+
+  //J. M.Quesada (May. 08)
+public:
+  // inline G4double GetTransitionProb1() const
+  G4double GetTransitionProb1() 
+  {
+    return TransitionProb1;
+  }
+  // inline G4double GetTransitionProb2() const
+  G4double GetTransitionProb2() 
+  {
+    return TransitionProb2;
+  }
+  // inline G4double GetTransitionProb3() const
+  G4double GetTransitionProb3() 
+  {
+    return TransitionProb3;
+  }
+
 };
 
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTriton.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTriton.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundTriton.hh	(revision 962)
@@ -25,8 +25,7 @@
 //
 //
-// $Id: G4PreCompoundTriton.hh,v 1.6 2007/10/01 10:42:00 ahoward Exp $
-// GEANT4 tag $Name:  $
+// by V. Lara
 //
-// by V. Lara
+//J. M. Quesada (July 08) 
 
 #ifndef G4PreCompoundTriton_h
@@ -36,6 +35,6 @@
 #include "G4ReactionProduct.hh"
 #include "G4Triton.hh"
-
 #include "G4TritonCoulombBarrier.hh"
+#include "G4PreCompoundParameters.hh"
 
 
@@ -53,81 +52,54 @@
 
   // operators  
-  const G4PreCompoundTriton & operator=(const G4PreCompoundTriton &right) {
-    if (&right != this) this->G4PreCompoundIon::operator=(right);
-    return *this;
-  }
+//  const G4PreCompoundTriton & operator=(const G4PreCompoundTriton &right) {
+//    if (&right != this) this->G4PreCompoundIon::operator=(right);
+//    return *this;
+//  }
 
-  G4bool operator==(const G4PreCompoundTriton &right) const
-  { return G4PreCompoundIon::operator==(right);}
+//  G4bool operator==(const G4PreCompoundTriton &right) const
+//  { return G4PreCompoundIon::operator==(right);}
 
   
-  G4bool operator!=(const G4PreCompoundTriton &right) const
-  { return G4PreCompoundIon::operator!=(right);}
+//  G4bool operator!=(const G4PreCompoundTriton &right) const
+//  { return G4PreCompoundIon::operator!=(right);}
 
 
-  G4ReactionProduct * GetReactionProduct() const
-  {
-    G4ReactionProduct * theReactionProduct =
-      new G4ReactionProduct(G4Triton::TritonDefinition());
-    theReactionProduct->SetMomentum(GetMomentum().vect());
-    theReactionProduct->SetTotalEnergy(GetMomentum().e());
-#ifdef PRECOMPOUND_TEST
-    theReactionProduct->SetCreatorModel("G4PrecompoundModel");
-#endif
-    return theReactionProduct;
-  }   
-    
+  G4ReactionProduct * GetReactionProduct() const;
+ 
+
 private:
 
-// added Rj method according to literature and JMQ
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged)
-  {
-    G4double rj = 1.0;
-    G4double denominator = NumberParticles*(NumberParticles-1)*(NumberParticles-2);
-    if(denominator != 0) rj = 3.0*static_cast<G4double>(NumberCharged*(NumberParticles-NumberCharged)*(NumberParticles-NumberCharged-1))/static_cast<G4double>(denominator); //JMQ 23/8/07
+  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged);
 
-    return rj;
-  }
+  virtual G4double CrossSection(const  G4double K) ; 
 
-  virtual G4double GetAlpha()
-  {
-    G4double C = 0.0;
-    G4double aZ = GetZ() + GetRestZ();
-    if (aZ >= 70) 
-      {
-	C = 0.10;
-      } 
-    else 
-      {
-	C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375; 
-      }
- 
-    return 1.0 + C/3.0;
-  }
+  virtual G4double FactorialFactor(const G4double N, const G4double P);
 
-  virtual G4double GetBeta()
-  {
-    return -GetCoulombBarrier();
-  }
+  virtual G4double CoalescenceFactor(const G4double A);
 
-  virtual G4double FactorialFactor(const G4double N, const G4double P)
-  {
-    return 
-      (N-3.0)*(P-2.0)*(
-		       (((N-2.0)*(P-1.0))/2.0) *(
-						 (((N-1.0)*P)/3.0) 
-						 )
-		       );
-  }
+  G4double GetOpt0(const G4double K);
+  G4double GetOpt12(const G4double K);
+  G4double GetOpt34(const G4double K);
 
-  virtual G4double CoalescenceFactor(const G4double A)
-  {
-    return 243.0/(A*A);
-  }    
-private:
+  G4double GetAlpha();
+  
+  G4double GetBeta();
 
-  G4TritonCoulombBarrier theTritonCoulombBarrier;
+//data members
 
+      G4TritonCoulombBarrier theTritonCoulombBarrier;
+       G4double ResidualA;
+      G4double ResidualZ; 
+      G4double theA;
+      G4double theZ;
+      G4double ResidualAthrd;
+      G4double FragmentA;
+      G4double FragmentAthrd;
 };
 
 #endif
+
+
+
+    
+
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.hh	(revision 962)
@@ -25,8 +25,14 @@
 //
 //
-// $Id: G4VPreCompoundFragment.hh,v 1.4 2006/08/20 01:07:28 dennis Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4VPreCompoundFragment.hh,v 1.10 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
-// by V. Lara
+// J. M. Quesada (August 2008).  
+// Based  on previous work by V. Lara
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option 
+// JMQ (06 September 2008) Also external choice has been added for:
+//                      - superimposed Coulomb barrier (if useSICB=true) 
 
 #ifndef G4VPreCompoundFragment_h
@@ -60,4 +66,5 @@
 			 G4VCoulombBarrier * aCoulombBarrier,
 			 const G4String &  aName);
+
   
   virtual ~G4VPreCompoundFragment();
@@ -135,8 +142,15 @@
   inline void IncrementStage();
 
+  //for inverse cross section choice
+  inline void SetOPTxs(G4int);
+  //for superimposed Coulomb Barrier for inverse cross sections
+  inline void UseSICB(G4bool);
+
+
 
   // =============
   // Data members
   // =============
+
 
 private:
@@ -145,4 +159,5 @@
   
   G4double theZ;
+private:
   
   G4double theRestNucleusA;
@@ -155,5 +170,5 @@
   
   G4double theBindingEnergy;
-  
+
   G4double theMaximalKineticEnergy;
   
@@ -165,5 +180,11 @@
   G4String theFragmentName;
 
-  G4int theStage;    
+  G4int theStage; 
+
+protected:
+//for inverse cross section choice
+  G4int OPTxs;
+//for superimposed Coulomb Barrier for inverse cross sections
+  G4bool useSICB;
 };
 
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.icc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.icc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundFragment.icc	(revision 962)
@@ -25,9 +25,13 @@
 //
 //
-// $Id: G4VPreCompoundFragment.icc,v 1.4 2006/08/20 01:07:40 dennis Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4VPreCompoundFragment.icc,v 1.7 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
-
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option 
+// JMQ (06 September 2008) Also external choice has been added for:
+//                      - superimposed Coulomb barrier (if useSICB=true) 
 
 inline G4double G4VPreCompoundFragment::GetA() const 
@@ -137,2 +141,15 @@
 }
 
+//for inverse cross section choice
+inline void G4VPreCompoundFragment::SetOPTxs(G4int opt)
+{
+  OPTxs=opt;
+  return;
+}
+
+//for superimposed Coulomb Barrier for inverse cross sections
+inline void G4VPreCompoundFragment::UseSICB(G4bool use)
+{
+  useSICB=use;
+  return;
+}
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundIon.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundIon.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundIon.hh	(revision 962)
@@ -25,70 +25,17 @@
 //
 //
-// $Id: G4VPreCompoundIon.hh,v 1.5 2007/08/23 16:29:01 ahoward Exp $
-// GEANT4 tag $Name:  $
+
+// $Id: G4VPreCompoundIon.hh,v 1.9 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
-
-#ifndef G4VPreCompoundIon_h
-#define G4VPreCompoundIon_h 1
-
-#include "G4VPreCompoundFragment.hh"
-#include "G4VCoulombBarrier.hh"
+//
+//J.M. Quesada (Apr. 2008) DUMMY abstract base class for ions.
+// Not ihherited by anything
+// Suppressed
 
 
-class G4VPreCompoundIon : public G4VPreCompoundFragment
-{
-protected:
-  // default constructor
-  G4VPreCompoundIon() {}
 
-public:
-
-  // copy constructor
-  G4VPreCompoundIon(const G4VPreCompoundIon &right): 
-    G4VPreCompoundFragment(right) {}
-    
-  // constructor  
-  G4VPreCompoundIon(const G4double anA, 
-		    const G4double aZ, 
-		    G4VCoulombBarrier* aCoulombBarrier,
-		    const G4String & aName): 
-    G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName) {}
-    
-  virtual ~G4VPreCompoundIon() {}
-    
-  // operators  
-  const G4VPreCompoundIon & 
-  operator=(const G4VPreCompoundIon &right) {
-    if (&right != this) this->G4VPreCompoundFragment::operator=(right);
-    return *this;
-  }
-    
-  G4bool operator==(const G4VPreCompoundIon &right) const 
-  { return G4VPreCompoundFragment::operator==(right);}
-    
-  G4bool operator!=(const G4VPreCompoundIon &right) const 
-  { return G4VPreCompoundFragment::operator!=(right);}
-    
-  virtual G4double ProbabilityDistributionFunction(const G4double eKin,
-						   const G4Fragment& aFragment);
-
-protected:
-  G4bool IsItPossible(const G4Fragment& aFragment) 
-  {
-    G4int pplus = aFragment.GetNumberOfCharged();   
-    G4int pneut = aFragment.GetNumberOfParticles()-pplus;
-    return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
-  }
-
-  virtual G4double GetAlpha() = 0;
-  virtual G4double GetBeta() = 0;
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0;
-  virtual G4double FactorialFactor(const G4double N, const G4double P) = 0;
-  virtual G4double CoalescenceFactor(const G4double A) = 0; 
-    
-};
-
-#endif
+ 
 
 
@@ -98,3 +45,2 @@
 
 
-
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundNucleon.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundNucleon.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundNucleon.hh	(revision 962)
@@ -25,64 +25,10 @@
 //
 //
-// $Id: G4VPreCompoundNucleon.hh,v 1.4 2007/07/23 09:56:40 ahoward Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4VPreCompoundNucleon.hh,v 1.8 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
+//J.M. Quesada (Apr. 2008) DUMMY abstract base class for nucleons. 
+// Not ihherited by anything
+// Suppressed
 
-#ifndef G4VPreCompoundNucleon_h
-#define G4VPreCompoundNucleon_h 1
-
-#include "G4VPreCompoundFragment.hh"
-#include "G4VCoulombBarrier.hh"
-
-
-class G4VPreCompoundNucleon : public G4VPreCompoundFragment
-{
-protected:
-  // default constructor
-  G4VPreCompoundNucleon() {}
-
-public:
-
-  // copy constructor
-  G4VPreCompoundNucleon(const G4VPreCompoundNucleon &right): 
-    G4VPreCompoundFragment(right) {}
-
-  // constructor  
-  G4VPreCompoundNucleon(const G4double anA, 
-			const G4double aZ, 
-			G4VCoulombBarrier* aCoulombBarrier,
-			const G4String & aName): 
-    G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName) {}
-
-  virtual ~G4VPreCompoundNucleon() {}
-
-  // operators  
-  const G4VPreCompoundNucleon & 
-  operator=(const G4VPreCompoundNucleon &right) {
-    if (&right != this) this->G4VPreCompoundFragment::operator=(right);
-    return *this;
-  }
-
-  G4bool operator==(const G4VPreCompoundNucleon &right) const 
-  { return G4VPreCompoundFragment::operator==(right);}
-    
-  G4bool operator!=(const G4VPreCompoundNucleon &right) const 
-  { return G4VPreCompoundFragment::operator!=(right);}
-    
-  virtual G4double ProbabilityDistributionFunction(const G4double eKin,
-						   const G4Fragment& aFragment);
-
-    
-protected:
-
-// added Rj method according to literature and JMQ
-  virtual G4double GetRj(const G4int NumberParticles, const G4int NumberCharged) = 0;
-  virtual G4double GetAlpha() = 0;
-  virtual G4double GetBeta() = 0;
-  virtual G4bool IsItPossible(const G4Fragment&) = 0; 
-    
-    
-};
-
-#endif
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundTransitions.hh
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundTransitions.hh	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4VPreCompoundTransitions.hh	(revision 962)
@@ -24,4 +24,9 @@
 // ********************************************************************
 //
+//J. M. Quesada (May 08). New virtual classes have been added 
+// JMQ (06 September 2008) Also external choices have been added for:
+//                      - "never go back"  hipothesis (useNGB=true) 
+//                      - CEM transition probabilities (useCEMtr=true)  
+
 #ifndef G4VPreCompoundTransitions_hh
 #define G4VPreCompoundTransitions_hh 1
@@ -33,10 +38,22 @@
 public:
 
-  G4VPreCompoundTransitions() {}
+  G4VPreCompoundTransitions():useNGB(false),useCEMtr(false) {}
   virtual ~G4VPreCompoundTransitions() {}
 
   virtual G4double CalculateProbability(const G4Fragment& aFragment) = 0;
   virtual G4Fragment PerformTransition(const G4Fragment&  aFragment) = 0;
+//J. M. Quesada (May.08) New virtual classes
+  virtual G4double GetTransitionProb1()=0;
+  virtual G4double GetTransitionProb2()=0;
+  virtual G4double GetTransitionProb3()=0;
 
+  // for never go back hypothesis (if useNGB=true, default=false)
+  inline void UseNGB(G4bool use){useNGB=use;}
+  //for use of CEM transition probabilities (if useCEMtr=true, defaut false)
+  inline void UseCEMtr(G4bool use){useCEMtr=use;}
+
+protected:
+  G4bool useNGB;
+  G4bool useCEMtr;
 };
 
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc	(revision 962)
@@ -25,7 +25,18 @@
 //
 //
-// Hadronic Process: Nuclear Preequilibrium
-// by V. Lara 
-
+// $Id: G4PreCompoundEmission.cc,v 1.19 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     G4PreCompoundEmission
+//
+// Author:         V.Lara
+//
+// Modified:  
+//
 
 #include "G4PreCompoundEmission.hh"
@@ -42,10 +53,10 @@
 
 
-  G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
+G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
 {
   return false;
 }
 
-    G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
+G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
 {
   return true;
@@ -95,6 +106,4 @@
   return;
 }
-
-
 
 G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)
@@ -163,4 +172,5 @@
   // Update nucleus parameters:
   // --------------------------
+
   // Number of excitons
   aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
@@ -175,5 +185,5 @@
   // Charge
   aFragment.SetZ(theFragment->GetRestZ());
-    
+
     
   // Perform Lorentz boosts
@@ -275,7 +285,46 @@
 }
 
-
 G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g, 
 				    const G4double E, const G4double Ef) const
+{
+	
+  G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g);
+  G4double alpha = (p*p+h*h)/(2.0*g);
+  
+  if ( (E-alpha) < 0 ) return 0;
+
+  G4double factp=factorial(p);
+
+  G4double facth=factorial(h);
+
+  G4double factph=factorial(p+h-1);
+  
+  G4double logConst =  (p+h)*std::log(g) - std::log (factph) - std::log(factp) - std::log(facth);
+
+// initialise values using j=0
+
+  G4double t1=1;
+  G4double t2=1;
+  G4double logt3=(p+h-1) * std::log(E-Aph);
+  G4double tot = std::exp( logt3 + logConst );
+
+// and now sum rest of terms 
+  G4int j(1);  
+  while ( (j <= h) && ((E - alpha - j*Ef) > 0.0) ) 
+    {
+	  t1 *= -1.;
+	  t2 *= (h+1-j)/j;
+	  logt3 = (p+h-1) * std::log( E - j*Ef - Aph) + logConst;
+	  G4double t3 = std::exp(logt3);
+	  tot += t1*t2*t3;
+	  j++;
+    }
+        
+  return tot;
+}
+
+
+
+G4double G4PreCompoundEmission::factorial(G4double a) const
 {
   // Values of factorial function from 0 to 60
@@ -349,77 +398,21 @@
   //      fact[n] = fact[n-1]*static_cast<G4double>(n); 
   //    }
-	
-  G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g);
-  G4double alpha = (p*p+h*h)/(2.0*g);
-
-  G4double tot = 0.0;
-  for (G4int j = 0; j <= h; j++) 
-    {
-      if (E-alpha-static_cast<G4double>(j)*Ef > 0.0)
-	{
-	  G4double t1 = std::pow(-1.0, static_cast<G4double>(j));
-	  G4double t2 = fact[static_cast<G4int>(h)]/ (fact[static_cast<G4int>(h)-j]*fact[j]);
-	  G4double t3 = E - static_cast<G4double>(j)*Ef - Aph;
-	  t3 = std::pow(t3,p+h-1);
-	  tot += t1*t2*t3;
-	}
-      else 
+  G4double result(0.0);
+  G4int ia = static_cast<G4int>(a);
+  if (ia < factablesize) 
+    {
+      result = fact[ia];
+    }
+  else
+    {
+      result = fact[factablesize-1];
+      for (G4int n = factablesize; n < ia+1; ++n)
         {
-          break;
+          result *= static_cast<G4double>(n);
         }
     }
     
-
-  G4double factp(0.0);
-  G4int ph = static_cast<G4int>(p);
-  if (ph < factablesize) 
-    {
-      factp = fact[ph];
-    }
-  else
-    {
-      factp = fact[factablesize-1];
-      for (G4int n = factablesize; n < ph+1; ++n)
-        {
-          factp *= static_cast<G4double>(n);
-        }
-    }
-
-  G4double facth(0.0);
-  ph = static_cast<G4int>(h);
-  if (ph < factablesize) 
-    {
-      facth = fact[ph];
-    }
-  else
-    {
-      facth = fact[factablesize-1];
-      for (G4int n = factablesize; n < ph+1; ++n)
-        {
-          facth *= static_cast<G4double>(n);
-        }
-    }
-  G4double factph(0.0);
-  ph = static_cast<G4int>(p+h)-1;
-  if (ph < factablesize)
-    {
-      factph = fact[ph];
-    }
-  else
-    {
-      factph = fact[factablesize-1];
-      for (G4int n = factablesize; n < ph+1; ++n)
-        {
-          factph *= static_cast<G4double>(n);
-        }
-    }
-  
-  tot *= std::pow(g,p+h)/factph;
-  tot /= factp;
-  tot /= facth;
-    
-  return tot;
-}
-
+    return result;
+}
 
 
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragment.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragment.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragment.cc	(revision 962)
@@ -24,7 +24,12 @@
 // ********************************************************************
 //
+// $Id: G4PreCompoundFragment.cc,v 1.8 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
-// by V. Lara
- 
+// J. M. Quesada (August 2008).  
+// Based  on previous work by V. Lara
+// JMQ (06 September 2008) Also external choice has been added for:
+//                      - superimposed Coulomb barrier (if useSICB=true) 
+//
 #include "G4PreCompoundFragment.hh"
 
@@ -41,5 +46,6 @@
 		      const G4String & aName):
   G4VPreCompoundFragment(anA,aZ,aCoulombBarrier,aName)
-{}
+{
+}
 
 
@@ -71,20 +77,25 @@
 CalcEmissionProbability(const G4Fragment & aFragment)
 {
-  if (GetEnergyThreshold() <= 0.0) 
+// If  theCoulombBarrier effect is included in the emission probabilities
+//if (GetMaximalKineticEnergy() <= 0.0) 
+  G4double limit;
+  if(OPTxs==0 ||  useSICB) limit= theCoulombBarrier;
+  else limit=0.;
+  if (GetMaximalKineticEnergy() <= limit) 
     {
       theEmissionProbability = 0.0;
       return 0.0;
   }    
-  // Coulomb barrier is the lower limit 
-  // of integration over kinetic energy
-  G4double LowerLimit = theCoulombBarrier;
-  
-  // Excitation energy of nucleus after fragment emission is the upper limit
-  // of integration over kinetic energy
-  G4double UpperLimit = this->GetMaximalKineticEnergy();
+// If  theCoulombBarrier effect is included in the emission probabilities
+//  G4double LowerLimit = 0.;
+// Coulomb barrier is the lower limit 
+// of integration over kinetic energy
+  G4double LowerLimit = limit;
+
+// Excitation energy of nucleus after fragment emission is the upper limit of integration over kinetic energy
+  G4double UpperLimit = GetMaximalKineticEnergy();
   
   theEmissionProbability = 
     IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment);
-  
   return theEmissionProbability;
 }
@@ -130,5 +141,4 @@
     }
   Total  *= (Up-Low)/2.0;
-
   return Total;
 }
@@ -140,18 +150,19 @@
 GetKineticEnergy(const G4Fragment & aFragment) 
 {
-  G4double V = this->GetCoulombBarrier();
-  G4double Tmax =  this->GetMaximalKineticEnergy() - V;
-  
+
+//	G4double V = this->GetCoulombBarrier();// alternative way for accessing the Coulomb barrier
+//                                             //should be equivalent (in fact it is)
+  G4double V;
+  if(OPTxs==0 || useSICB) V= theCoulombBarrier;//let's keep this way for consistency with CalcEmissionProbability method
+  else V=0.;
+
+  G4double Tmax =  GetMaximalKineticEnergy() ;  
   G4double T(0.0);
   G4double NormalizedProbability(1.0);
   do 
     {
-      T = V + G4UniformRand()*Tmax;
-      NormalizedProbability = this->ProbabilityDistributionFunction(T,aFragment)/
-	this->GetEmissionProbability();
-      
-    } 
-  while (G4UniformRand() > NormalizedProbability);
-  
+      T =V+ G4UniformRand()*(Tmax-V);
+      NormalizedProbability = ProbabilityDistributionFunction(T,aFragment)/GetEmissionProbability();      
+    }   while (G4UniformRand() > NormalizedProbability);  
   return T;
 }
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragmentVector.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragmentVector.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundFragmentVector.cc	(revision 962)
@@ -24,7 +24,6 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundFragmentVector.cc,v 1.5 2006/06/29 20:59:21 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4PreCompoundFragmentVector.cc,v 1.11 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear Preequilibrium
@@ -86,4 +85,5 @@
   for (i = theChannels->begin(); i != theChannels->end(); ++i) {
     accumulation += (*i)->GetEmissionProbability();
+
     running.push_back(accumulation);
   }
@@ -126,5 +126,5 @@
 	  (*i)->IncrementStage();
 	}
-    } 
+    }
 
   return theChannels->operator[](ChosenChannel);
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundIon.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundIon.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundIon.cc	(revision 962)
@@ -24,10 +24,29 @@
 // ********************************************************************
 //
-// by V. Lara
+// $Id: G4PreCompoundIon.cc,v 1.16 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     G4PreCompoundIon
+//
+// Author:         V.Lara
+//
+// Modified:  
+// 10.02.2009 J. M. Quesada fixed bug in density level of light fragments  
+//
 
 #include "G4PreCompoundIon.hh"
 #include "G4PreCompoundParameters.hh"
-//#include "G4EvaporationLevelDensityParameter.hh"
 
+G4bool G4PreCompoundIon::IsItPossible(const G4Fragment& aFragment) 
+{
+  G4int pplus = aFragment.GetNumberOfCharged();   
+  G4int pneut = aFragment.GetNumberOfParticles()-pplus;
+  return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
+}
 
 G4double G4PreCompoundIon::
@@ -43,38 +62,43 @@
   G4double H = aFragment.GetNumberOfHoles();
   G4double N = P + H;
-  
-  // g = (6.0/pi2)*a*A
-  //  G4EvaporationLevelDensityParameter theLDP;
+
   G4double g0 = (6.0/pi2)*aFragment.GetA() * 
     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);
+ 
   G4double g1 = (6.0/pi2)*GetRestA() * 
     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(GetRestA(),GetRestZ(),U);
-  //no longer used:  G4double gj = (6.0/pi2)*GetA() *
-  //  -----"-----    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(GetA(),GetZ(),U);
 
-  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
-  G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1);
-  G4double Aj = GetA()*(GetA()+1.0)/4.0;
+  //JMQ 06/02/209  This is  THE BUG that was killing cluster emission
+  // G4double gj = (6.0/pi2)*GetA() *
+  //   G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+
+  G4double gj = g1;
+
+  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0; 
+
+  G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1); 
+
+  G4double Aj = GetA()*(GetA()+1.0)/4.0/gj; 
+
 
   G4double E0 = std::max(0.0,U - A0);
   if (E0 == 0.0) return 0.0;
-  //  G4double E1 = std::max(0.0,GetMaximalKineticEnergy() + GetCoulombBarrier() - eKin - A1);
-  G4double E1 = (std::max(0.0,GetMaximalKineticEnergy() - eKin - A1))/g1; //JMQ fix
-  //  G4double Ej = std::max(0.0,U + GetBindingEnergy() - Aj);
-  G4double Ej = std::max(0.0,eKin + GetBindingEnergy() - Aj); //JMQ fix
 
-  G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*(eKin+GetBindingEnergy()))))*
-      GetAlpha()*(eKin+GetBeta())/(r0*std::pow(GetRestA(),1.0/3.0)) *
-      CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P);
+  G4double E1 = (std::max(0.0,GetMaximalKineticEnergy() - eKin - A1)); 
+
+  G4double Ej = std::max(0.0,eKin + GetBindingEnergy() -Aj); 
+
+  // JMQ 10/02/09 reshaping of the formula (unnecessary std::pow elimitated)
+  G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*
+		(eKin+GetBindingEnergy()))))/(pi * r0 * r0 *r0* GetRestA())* 
+                eKin*CrossSection(eKin)*millibarn* 
+                CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P)* 
+                GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged());
 
   G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0);
+ 
+  G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0; 
 
-  //  G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0;
-  G4double pC = std::pow((g1*Ej)/(g0*E0),GetA()-1.0)*(g1/g0)/E0; //JMQ fix
-
-  G4double Probability = pA * pB * pC * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged());
+  G4double Probability = pA * pB * pC;
 
   return Probability;
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundModel.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundModel.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundModel.cc	(revision 962)
@@ -25,8 +25,21 @@
 //
 //
-// $Id: G4PreCompoundModel.cc,v 1.11 2007/10/11 14:19:36 ahoward Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4PreCompoundModel.cc,v 1.17 2008/12/09 14:09:59 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
+//
+//J. M. Quesada (Apr.08). Several changes. Soft cut-off switched off. 
+//(May. 08). Protection against non-physical preeq. transitional regime has 
+// been set
+//
+// 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 (useSICB=true) 
+//                      - "never go back"  hipothesis (useNGB=true) 
+//                      - soft cutoff from preeq. to equlibrium (useSCO=true)
+//                      - CEM transition probabilities (useCEMtr=true)  
+
 
 #include "G4PreCompoundModel.hh"
@@ -34,4 +47,5 @@
 #include "G4PreCompoundTransitions.hh"
 #include "G4GNASHTransitions.hh"
+#include "G4ParticleDefinition.hh"
 
 
@@ -160,4 +174,11 @@
   // Copy of the initial state 
   G4Fragment aFragment(theInitialState);
+
+  if (aFragment.GetExcitationEnergy() < 10*eV)
+    {
+      // Perform Equilibrium Emission
+      PerformEquilibriumEmission(aFragment,Result);
+      return Result;
+    }
   
   if (aFragment.GetA() < 5) {
@@ -175,5 +196,20 @@
   if (useHETCEmission) aEmission.SetHETCModel();
   aEmission.SetUp(theInitialState);
-
+  
+  //for cross section options
+  
+  if (OPTxs!= 0 && OPTxs!=1 && OPTxs !=2 && OPTxs !=3 && OPTxs !=4  ) {
+    std::ostringstream errOs;
+    errOs << "BAD CROSS SECTION OPTION in G4PreCompoundModel.cc !!"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());}
+  else aEmission.SetOPTxs(OPTxs); 
+  
+  //for the choice of superimposed Coulomb Barrier for inverse cross sections
+  
+   aEmission.UseSICB(useSICB);
+ 
+  
+  //----------
+  
   G4VPreCompoundTransitions * aTransition = 0;
   if (useGNASHTransition) 
@@ -184,48 +220,93 @@
     {
       aTransition = new G4PreCompoundTransitions;
-    }
-
-
+      // for the choice of "never go back" hypothesis and CEM transition probabilities  
+      if (useNGB) aTransition->UseNGB(useNGB);
+      if (useCEMtr) aTransition->UseCEMtr(useCEMtr);
+    }
+  
   // Main loop. It is performed until equilibrium deexcitation.
+  //G4int fragment=0;
+  
   for (;;) {
+    
+    //fragment++;
+    //G4cout<<"-------------------------------------------------------------------"<<G4endl;
+    //G4cout<<"Fragment number .. "<<fragment<<G4endl;
+    
     // Initialize fragment according with the nucleus parameters
     aEmission.Initialize(aFragment);
-
-    // Equilibrium exciton number
-    //    G4EvaporationLevelDensityParameter theLDP;
-    //    G4double g = (6.0/pi2)*aFragment.GetA()*
-    //      theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),
-    //				   aFragment.GetExcitationEnergy());
+    
+    
     
     G4double g = (6.0/pi2)*aFragment.GetA()*
       G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-    G4int EquilibriumExcitonNumber = static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy())
-                                                        + 0.5);
-    
-    // Loop for transitions, it is performed while there are preequilibrium transitions.
+    
+    
+    
+    
+    G4int EquilibriumExcitonNumber = static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy())+ 0.5);
+//   
+//    G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
+//
+// J. M. Quesada (Jan. 08)  equilibrium hole number could be used as preeq.- evap. delimiter (IAEA report)
+//    G4int EquilibriumHoleNumber = static_cast<G4int>(0.2*std::sqrt(g*aFragment.GetExcitationEnergy())+ 0.5);
+    
+// Loop for transitions, it is performed while there are preequilibrium transitions.
     G4bool ThereIsTransition = false;
+    
+    //        G4cout<<"----------------------------------------"<<G4endl;
+    //        G4double NP=aFragment.GetNumberOfParticles();
+    //        G4double NH=aFragment.GetNumberOfHoles();
+    //        G4double NE=aFragment.GetNumberOfExcitons();
+    //        G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
+    //        G4cout<<"N. excitons="<<NE<<"  N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
+        
+    
+    //G4int transition=0;
     do 
       {
+        //transition++;
+        //G4cout<<"transition number .."<<transition<<G4endl;
+        //G4cout<<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
         G4bool go_ahead = false;
+        // soft cutoff criterium as an "ad-hoc" solution to force increase in  evaporation  
+        //       G4double test = static_cast<G4double>(aFragment.GetNumberOfHoles());
         G4double test = static_cast<G4double>(aFragment.GetNumberOfExcitons());
-        if (test < EquilibriumExcitonNumber)
-          {
-            test /= static_cast<G4double>(EquilibriumExcitonNumber);
-            test -= 1.0;
-            test = test*test;
-            test /= 0.32;
-            test = 1.0 - std::exp(-test);
-            go_ahead = (G4UniformRand() < test);
-          }
-      
+        
+        
+        if (test < EquilibriumExcitonNumber) go_ahead=true;
+        //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
+        if (useSCO) {
+          if (test < EquilibriumExcitonNumber)
+            //  if (test < EquilibriumHoleNumber)
+            {
+              test /= static_cast<G4double>(EquilibriumExcitonNumber); 
+              //     test /= static_cast<G4double>(EquilibriumHoleNumber);
+              test -= 1.0;
+              test = test*test;
+              test /= 0.32;
+              test = 1.0 - std::exp(-test);
+              go_ahead = (G4UniformRand() < test);
+              
+            }
+        } 
+        
+        //JMQ: WARNING:  CalculateProbability MUST be called prior to Get methods !! (O values would be returned otherwise)
+        G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);
+        G4double P1=aTransition->GetTransitionProb1();
+        G4double P2=aTransition->GetTransitionProb2();
+        G4double P3=aTransition->GetTransitionProb3();
+        //       G4cout<<"P1="<<P1<<" P2="<<P2<<"  P3="<<P3<<G4endl;
+        
+        
+        //J.M. Quesada (May. 08). Physical criterium (lamdas)  PREVAILS over approximation (critical exciton number)
+        if(P1<=(P2+P3)) go_ahead=false;
+        
         if (go_ahead &&  aFragment.GetA() > 4)
-	  {
-	    //  	if (aFragment.GetNumberOfParticles() < 1) {
-	    //  	  aFragment.SetNumberOfHoles(aFragment.GetNumberOfHoles()+1);
-	    //  	  aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()+1);       
-	    //  	}
-				
+	  {				
 	    G4double TotalEmissionProbability = aEmission.GetTotalProbability(aFragment);
-      	
+            //
+            //  G4cout<<"TotalEmissionProbability="<<TotalEmissionProbability<<G4endl;
+            //
 	    // Check if number of excitons is greater than 0
 	    // else perform equilibrium emission
@@ -242,17 +323,17 @@
 	    
 	    //	    G4PreCompoundTransitions aTransition(aFragment);
-	
+            
+            //J.M.Quesada (May 08) this has already been done in order to decide what to do (preeq-eq) 
 	    // Sum of transition probabilities
-	    G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);
-	
+            //	    G4double TotalTransitionProbability = aTransition->CalculateProbability(aFragment);
+            
 	    // Sum of all probabilities
 	    G4double TotalProbability = TotalEmissionProbability + TotalTransitionProbability;
-	
+            
 	    // Select subprocess
 	    if (G4UniformRand() > TotalEmissionProbability/TotalProbability) 
 	      {
 		// It will be transition to state with a new number of excitons
-		ThereIsTransition = true;
-		
+		ThereIsTransition = true;		
 		// Perform the transition
 		aFragment = aTransition->PerformTransition(aFragment);
@@ -262,6 +343,4 @@
 		// It will be fragment emission
 		ThereIsTransition = false;
-
-		// Perform the emission and Add emitted fragment to Result
 		Result->push_back(aEmission.PerformEmission(aFragment));
 	      }
@@ -288,4 +367,5 @@
 {
   G4ReactionProductVector * theEquilibriumResult;
+
   theEquilibriumResult = GetExcitationHandler()->BreakItUp(aFragment);
   
@@ -353,2 +433,5 @@
 
 #endif
+
+
+
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNucleon.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNucleon.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundNucleon.cc	(revision 962)
@@ -24,9 +24,30 @@
 // ********************************************************************
 //
-// by V. Lara
+//
+// $Id: G4PreCompoundNucleon.cc,v 1.13 2009/02/11 18:06:00 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     G4PreCompoundNucleon
+//
+// Author:         V.Lara
+//
+// Modified:  
+// 10.02.2009 J. M. Quesada cleanup
+//
 
 #include "G4PreCompoundNucleon.hh"
 #include "G4PreCompoundParameters.hh"
 
+G4bool G4PreCompoundNucleon::IsItPossible(const G4Fragment& aFragment) 
+{
+  G4int pplus = aFragment.GetNumberOfCharged();   
+  G4int pneut = aFragment.GetNumberOfParticles()-pplus;
+  return (pneut >= (GetA()-GetZ()) && pplus >= GetZ());
+}
 
 G4double G4PreCompoundNucleon::
@@ -36,6 +57,4 @@
   if ( !IsItPossible(aFragment) ) return 0.0;
 
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-
   G4double U = aFragment.GetExcitationEnergy();
   G4double P = aFragment.GetNumberOfParticles();
@@ -43,15 +62,13 @@
   G4double N = P + H;
   
-  // g = (6.0/pi2)*a*A
-  //  G4EvaporationLevelDensityParameter theLDP;
   G4double g0 = (6.0/pi2)*aFragment.GetA() * 
     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);
+ 
   G4double g1 = (6.0/pi2)*GetRestA() * 
     G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-    //    theLDP.LevelDensityParameter(static_cast<G4int>(GetRestA()),static_cast<G4int>(GetRestZ()),U);
 
   G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
-  G4double A1 = (A0*g0 - P/2.0)/g1;
+
+  G4double A1 = (A0 - P/2.0)/g1;
   
   G4double E0 = std::max(0.0,U - A0);
@@ -61,9 +78,35 @@
 
 
-  G4double Probability = 2.0/(hbarc*hbarc*hbarc) * GetReducedMass() * 
-      GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) * r0 * r0 * std::pow(GetRestA(),2.0/3.0) * GetAlpha() * (eKin + GetBeta()) *
-      P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 *
-      g1/(g0*g0);
+  G4double Probability = 1.0/pi2*2.0/(hbarc*hbarc*hbarc) * GetReducedMass()
+    * GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())  
+    * eKin*CrossSection(eKin)*millibarn * P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 * g1/(g0*g0);
+
+  if (GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())<0.0 
+      || CrossSection(eKin) <0) {  
+    std::ostringstream errOs;
+    G4cout<<"WARNING:  NEGATIVE VALUES "<<G4endl;     
+    errOs << "Rj=" << GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged())
+	  <<G4endl;
+    errOs <<"  xsec("<<eKin<<" MeV) ="<<CrossSection(eKin)<<G4endl;
+    errOs <<"  A="<<GetA()<<"  Z="<<GetZ()<<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+  }
 
   return Probability;
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundParameters.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundParameters.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundParameters.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4PreCompoundParameters.cc,v 1.3 2006/06/29 20:59:29 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTransitions.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTransitions.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTransitions.cc	(revision 962)
@@ -24,9 +24,21 @@
 // ********************************************************************
 //
-//
-// $Id: G4PreCompoundTransitions.cc,v 1.13 2007/07/23 12:48:54 ahoward Exp $
-// GEANT4 tag $Name:  $
-//
-// by V. Lara
+// $Id: G4PreCompoundTransitions.cc,v 1.20 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     G4PreCompoundIon
+//
+// Author:         V.Lara
+//
+// Modified:  
+// 16.02.2008 J. M. Quesada fixed bugs 
+// 06.09.2008 J. M. Quesada added external choices for:
+//                      - "never go back"  hipothesis (useNGB=true) 
+//                      -  CEM transition probabilities (useCEMtr=true)
 
 #include "G4PreCompoundTransitions.hh"
@@ -55,4 +67,7 @@
 CalculateProbability(const G4Fragment & aFragment)
 {
+  //G4cout<<"In G4PreCompoundTransitions.cc  useNGB="<<useNGB<<G4endl;
+  //G4cout<<"In G4PreCompoundTransitions.cc  useCEMtr="<<useCEMtr<<G4endl;
+
   // Fermi energy
   const G4double FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
@@ -74,95 +89,143 @@
   G4double Z = static_cast<G4double>(aFragment.GetZ());
   G4double U = aFragment.GetExcitationEnergy();
-
-  // Relative Energy (T_{rel})
-  G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N;
-  
-  // Sample kind of nucleon-projectile 
-  G4bool ChargedNucleon(false);
-  G4double chtest = 0.5;
-  if (P > 0) chtest = aFragment.GetNumberOfCharged()/P;
-  if (G4UniformRand() < chtest) ChargedNucleon = true;
-
-  // Relative Velocity: 
-  // <V_{rel}>^2
-  G4double RelativeVelocitySqr(0.0);
-  if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2;
-  else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2;
-
-  // <V_{rel}>
-  G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
-
-  // Proton-Proton Cross Section
-  G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn;
-  // Proton-Neutron Cross Section
-  G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn;
-  
-  // Averaged Cross Section: \sigma(V_{rel})
-  //  G4double AveragedXSection = (ppXSection+npXSection)/2.0;
-  G4double AveragedXSection(0.0);
-  if (ChargedNucleon)
-    {
-      AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0);
-    }
-  else 
-    {
-      AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0);
-    }
-
-  // Fermi relative energy ratio
-  G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
-
-  // This factor is introduced to take into account the Pauli principle
-  G4double PauliFactor = 1.0 - (7.0/5.0)*FermiRelRatio;
-  if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
-
-  // Interaction volume 
-  G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
-
-  // Transition probability for \Delta n = +2
-  //  TransitionProb1 = 0.00332*AveragedXSection*PauliFactor*std::sqrt(RelativeEnergy)/
-  //    std::pow(1.2 + 1.0/(4.7*RelativeVelocity), 3.0);
-  TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint;
-  if (TransitionProb1 < 0.0) TransitionProb1 = 0.0; 
-
-  // g = (6.0/pi2)*aA; 
-  //  G4double a = theLDP.LevelDensityParameter(A,Z,U-G4PairingCorrection::GetInstance()->GetPairingCorrection(A,Z));
-  G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  // GE = g*E where E is Excitation Energy
-  G4double GE = (6.0/pi2)*a*A*U;
-
-
-  // F(p,h) = 0.25*(p^2 + h^2 + p - h) - 0.5*h
-  G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
-  G4bool NeverGoBack(false);
-  //AH  if (U-Fph < 0.0) NeverGoBack = true;
-  if (GE-Fph < 0.0) NeverGoBack = true;
-  // F(p+1,h+1)
-  G4double Fph1 = Fph + N/2.0;
-  // (n+1)/n ((g*E - F(p,h))/(g*E - F(p+1,h+1)))^(n+1)
-  G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0);
-
-
-  if (NeverGoBack)
-    {
+  
+  if(U<10*eV) return 0.0;
+  
+  //J. M. Quesada (Feb. 08) new physics
+  // OPT=1 Transitions are calculated according to Gudima's paper (original in G4PreCompound from VL) 
+  // OPT=2 Transitions are calculated according to Gupta's formulae
+  //
+  
+  
+  
+  if (useCEMtr){
+
+    
+    // Relative Energy (T_{rel})
+    G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N;
+    
+    // Sample kind of nucleon-projectile 
+    G4bool ChargedNucleon(false);
+    G4double chtest = 0.5;
+    if (P > 0) chtest = aFragment.GetNumberOfCharged()/P;
+    if (G4UniformRand() < chtest) ChargedNucleon = true;
+    
+    // Relative Velocity: 
+    // <V_{rel}>^2
+    G4double RelativeVelocitySqr(0.0);
+    if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2;
+    else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2;
+    
+    // <V_{rel}>
+    G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
+    
+    // Proton-Proton Cross Section
+    G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn;
+    // Proton-Neutron Cross Section
+    G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn;
+    
+    // Averaged Cross Section: \sigma(V_{rel})
+    //  G4double AveragedXSection = (ppXSection+npXSection)/2.0;
+    G4double AveragedXSection(0.0);
+    if (ChargedNucleon)
+      {
+        //JMQ: small bug fixed
+        //      AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0);
+        AveragedXSection = ((Z-1.0) * ppXSection + (A-Z) * npXSection) / (A-1.0);
+      }
+    else 
+      {
+        AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0);
+      }
+    
+    // Fermi relative energy ratio
+    G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
+    
+    // This factor is introduced to take into account the Pauli principle
+    G4double PauliFactor = 1.0 - (7.0/5.0)*FermiRelRatio;
+    if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
+    
+    // Interaction volume 
+    //  G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
+    G4double xx=2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity);
+    G4double Vint = (4.0/3.0)*pi*xx*xx*xx;
+    
+    // Transition probability for \Delta n = +2
+    
+    TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint;
+    if (TransitionProb1 < 0.0) TransitionProb1 = 0.0; 
+    
+    G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+    // GE = g*E where E is Excitation Energy
+    G4double GE = (6.0/pi2)*a*A*U;
+    
+    G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
+    
+    //G4bool NeverGoBack(false);
+    G4bool NeverGoBack;
+    if(useNGB)  NeverGoBack=true;
+    else NeverGoBack=false;
+    
+    
+    //JMQ/AH  bug fixed: if (U-Fph < 0.0) NeverGoBack = true;
+    if (GE-Fph < 0.0) NeverGoBack = true;
+    
+    // F(p+1,h+1)
+    G4double Fph1 = Fph + N/2.0;
+    
+    G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0);
+    
+    
+    if (NeverGoBack)
+      {
       TransitionProb2 = 0.0;
       TransitionProb3 = 0.0;
-    }
-  else 
-    {
-      // Transition probability for \Delta n = -2 (at F(p,h) = 0)
-      //  TransitionProb2 = max(0, (TransitionProb1*P*H*(P+H+1.0)*(P+H-2.0))/(GE*GE));
-      //  TransitionProb2 = (TransitionProb1*P*H*(P+H+1.0)*(P+H-2.0))/(GE*GE);
-      TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph));
+      }
+    else 
+      {
+        // Transition probability for \Delta n = -2 (at F(p,h) = 0)
+        TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph));
+        if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 
+        
+        // Transition probability for \Delta n = 0 (at F(p,h) = 0)
+        TransitionProb3 = TransitionProb1* ((N+1.0)/N) * ProbFactor  * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph);
+        if (TransitionProb3 < 0.0) TransitionProb3 = 0.0; 
+      }
+    
+    //  G4cout<<"U = "<<U<<G4endl;
+    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
+    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl; 
+    return TransitionProb1 + TransitionProb2 + TransitionProb3;
+  }
+  
+  else  {
+    //JMQ: Transition probabilities from Gupta's work
+    
+    G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
+    // GE = g*E where E is Excitation Energy
+    G4double GE = (6.0/pi2)*a*A*U;
+    
+    G4double Kmfp=2.;
+     
+    
+    TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
+    if (TransitionProb1 < 0.0) TransitionProb1 = 0.0;
+    
+    if (useNGB){
+      TransitionProb2=0.;
+      TransitionProb3=0.;
+    }
+    else{        
+      if (N<=1) TransitionProb2=0. ;    
+      else  TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);      
       if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 
-      
-      
-      // Transition probability for \Delta n = 0 (at F(p,h) = 0)
-      //  TransitionProb3 = TransitionProb1*(P+H+1.0)*(P*(P-1.0)+4.0*P*H+H*(H-1.0))/((P+H)*GE);
-      TransitionProb3 = TransitionProb1 * ProbFactor * ((N+1.0)/N) * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph);
-      if (TransitionProb3 < 0.0) TransitionProb3 = 0.0; 
-    }
-  
-  return TransitionProb1 + TransitionProb2 + TransitionProb3;
+      TransitionProb3=0.;
+    }
+    
+      //  G4cout<<"U = "<<U<<G4endl;
+    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
+    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl; 
+    return TransitionProb1 + TransitionProb2 + TransitionProb3;
+  }
 }
 
@@ -185,19 +248,21 @@
     }
 
-  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and in proportion to the number charges w.r.t. number of particles
-  if(deltaN < 0 && G4UniformRand() <= static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) && (result.GetNumberOfCharged() >= 1))
+  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and in proportion 
+  // to the number charges w.r.t. number of particles,  PROVIDED that there are charged particles
+  if(deltaN < 0 && G4UniformRand() <= 
+     static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) 
+     && (result.GetNumberOfCharged() >= 1)) {
     result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative!
-
-  // result.SetNumberOfParticles was here
-  // result.SetNumberOfHoles was here
-  // the following lines have to be before SetNumberOfCharged, otherwise the check on number of charged vs. number of particles fails
+  }
+
+  // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on 
+  // number of charged vs. number of particles fails
   result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
   result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2); 
 
-  // With weight Z/A, number of charged particles is decreased on +1
-  //  if ((deltaN > 0 || result.GetNumberOfCharged() > 0) && // AH/JMQ check is now in initialize within G4VPreCompoundFragment
+  // With weight Z/A, number of charged particles is increased with +1
   if ( ( deltaN > 0 ) &&
       (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/
-		  std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
+       std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
     {
       result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2);
@@ -210,8 +275,4 @@
     }
   
-  // moved from above to make code more readable
-  //  result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
-  //  result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2); 
-
   return result;
 }
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundFragment.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundFragment.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundFragment.cc	(revision 962)
@@ -24,9 +24,10 @@
 // ********************************************************************
 //
+// $Id: G4VPreCompoundFragment.cc,v 1.12 2009/02/10 16:01:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
-// $Id: G4VPreCompoundFragment.cc,v 1.4 2007/07/23 09:56:40 ahoward Exp $
-// GEANT4 tag $Name:  $
+// J. M. Quesada (August 2008).  
+// Based  on previous work by V. Lara
 //
-// by V. Lara
  
 #include "G4VPreCompoundFragment.hh"
@@ -142,6 +143,5 @@
   if ((theRestNucleusA < theRestNucleusZ) ||
       (theRestNucleusA < theA) ||
-      (theRestNucleusZ < theZ) ||
-      (aFragment.GetNumberOfCharged() < theZ)) // AH last argument from JMQ
+      (theRestNucleusZ < theZ)) 
     {
       // In order to be sure that emission probability will be 0.
@@ -155,5 +155,6 @@
     GetCoulombBarrier(static_cast<G4int>(theRestNucleusA),static_cast<G4int>(theRestNucleusZ),
 		      aFragment.GetExcitationEnergy());
-  
+
+
   // Compute Binding Energies for fragments 
   // (needed to separate a fragment from the nucleus)
@@ -164,9 +165,10 @@
   
   // Compute Maximal Kinetic Energy which can be carried by fragments after separation
+  // This is the true (assimptotic) maximal kinetic energy
   G4double m = aFragment.GetMomentum().m();
   G4double rm = GetRestNuclearMass();
   G4double em = GetNuclearMass();
   theMaximalKineticEnergy = ((m - rm)*(m + rm) + em*em)/(2.0*m) - em;
-  // - theCoulombBarrier;
+ 
   
   return;
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundIon.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundIon.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundIon.cc	(revision 962)
@@ -25,60 +25,12 @@
 //
 //
-// $Id: G4VPreCompoundIon.cc,v 1.5 2007/07/23 09:56:40 ahoward Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4VPreCompoundIon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // by V. Lara
-
-#include "G4VPreCompoundIon.hh"
-#include "G4PreCompoundParameters.hh"
-//#include "G4EvaporationLevelDensityParameter.hh"
-
-
-G4double G4VPreCompoundIon::
-ProbabilityDistributionFunction(const G4double eKin, 
-				const G4Fragment& aFragment)
-{
-  if ( !IsItPossible(aFragment) ) return 0.0;
-
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-
-  G4double U = aFragment.GetExcitationEnergy();
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
-  G4double N = P + H;
-
-  // g = (6.0/pi2)*a*A
-  //  G4EvaporationLevelDensityParameter theLDP;
-  G4double g0 = (6.0/pi2)*aFragment.GetA() * 
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(static_cast<G4int>(aFragment.GetA()),static_cast<G4int>(aFragment.GetZ()),U);
-  G4double g1 = (6.0/pi2)*GetRestA() * 
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(GetRestA(),GetRestZ(),U);
-  G4double gj = (6.0/pi2)*GetA() *
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(GetA(),GetZ(),U);
-
-  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
-  G4double A1 = std::max(0.0,(A0*g0 + GetA()*(GetA()-2.0*P-1.0)/4.0)/g1);
-  G4double Aj = GetA()*(GetA()+1.0)/4.0;
-
-  G4double E0 = std::max(0.0,U - A0);
-  if (E0 == 0.0) return 0.0;
-  G4double E1 = std::max(0.0,GetMaximalKineticEnergy() + GetCoulombBarrier() - eKin - A1);
-  G4double Ej = std::max(0.0,U + GetBindingEnergy() - Aj);
-
-  G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()*(eKin+GetBindingEnergy()))))*
-      GetAlpha()*(eKin+GetBeta())/(r0*std::pow(GetRestA(),1.0/3.0)) *
-      CoalescenceFactor(aFragment.GetA()) * FactorialFactor(N,P);
-
-  G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0);
-
-  G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0;
-
-  G4double Probability = pA * pB * pC;
-
-  return Probability;
-}
+//
+//J.M. Quesada (Apr. 2008) DUMMY abstract base class for ions
+// it is not ihherited by anything
+// Suppressed
 
 
@@ -88,2 +40,7 @@
 
 
+
+
+
+
+
Index: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundNucleon.cc
===================================================================
--- trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundNucleon.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4VPreCompoundNucleon.cc	(revision 962)
@@ -25,51 +25,14 @@
 //
 //
-// $Id: G4VPreCompoundNucleon.cc,v 1.5 2007/07/23 09:56:40 ahoward Exp $
-// GEANT4 tag $Name:  $
+
+// $Id: G4VPreCompoundNucleon.cc,v 1.9 2008/09/22 10:18:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
+
 //
 // by V. Lara
-
-#include "G4VPreCompoundNucleon.hh"
-#include "G4PreCompoundParameters.hh"
-//#include "G4EvaporationLevelDensityParameter.hh"
+//
+//J.M. Quesada (Apr. 2008) DUMMY abstract base class for nucleons.
+// it is not ihherited by anything
+// Suppressed
 
 
-G4double G4VPreCompoundNucleon::
-ProbabilityDistributionFunction(const G4double eKin, 
-				const G4Fragment& aFragment)
-{
-  if ( !IsItPossible(aFragment) ) return 0.0;
-
-  const G4double r0 = G4PreCompoundParameters::GetAddress()->Getr0();
-
-  G4double U = aFragment.GetExcitationEnergy();
-  G4double P = aFragment.GetNumberOfParticles();
-  G4double H = aFragment.GetNumberOfHoles();
-  G4double N = P + H;
-
-  // g = (6.0/pi2)*a*A
-  //  G4EvaporationLevelDensityParameter theLDP;
-  G4double g0 = (6.0/pi2)*aFragment.GetA() * 
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-  //    theLDP.LevelDensityParameter(G4int(aFragment.GetA()),G4int(aFragment.GetZ()),U);
-  G4double g1 = (6.0/pi2)*GetRestA() * 
-    G4PreCompoundParameters::GetAddress()->GetLevelDensity();
-    //    theLDP.LevelDensityParameter(G4int(GetRestA()),G4int(GetRestZ()),U);
-
-
-  G4double A0 = ((P*P+H*H+P-H)/4.0 - H/2.0)/g0;
-  G4double A1 = (A0*g0 - P/2.0)/g1;
-
-  G4double E0 = std::max(0.0,U - A0);
-  if (E0 == 0.0) return 0.0;
-  G4double E1 = std::max(0.0,U - eKin - GetBindingEnergy() - A1);
-
-  G4double Probability = 2.0/(hbarc*hbarc*hbarc) * GetReducedMass() * 
-      GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()) * r0 * r0 * std::pow(GetRestA(),2.0/3.0) * GetAlpha() * (eKin + GetBeta()) *
-      P*(N-1.0) * std::pow(g1*E1/(g0*E0),N-2.0)/E0 *
-      g1/(g0*g0);
-
-  return Probability;
-}
-
-
