Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationChannel.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationChannel.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4AlphaEvaporationChannel.cc,v 1.4 2006/06/29 20:10:17 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc	(revision 962)
@@ -24,44 +24,18 @@
 // ********************************************************************
 //
-//
-// $Id: G4AlphaEvaporationProbability.cc,v 1.4 2006/06/29 20:10:19 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
-// by V. Lara (Nov 1999)
-//
-
+// by V. Lara (Oct 1998)
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
 
 #include "G4AlphaEvaporationProbability.hh"
 
 G4AlphaEvaporationProbability::G4AlphaEvaporationProbability() :
-    G4EvaporationProbability(4,2,4) // A,Z,Gamma
-{
-    //  const G4int NumExcitedStates = 31+1;
-    std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
-    std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
-    ExcitEnergies.reserve(NumExcitedStatesEnergy);
-    ExcitSpins.reserve(NumExcitedStatesSpin);
-    ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
-    ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
-//    for (G4int i = 0; i < NumExcitedStates; i++) {
-//      ExcitEnergies(i) = 0.0;
-//      ExcitSpins(i) = 0;
-//    }
-
-    ExcitEnergies[18] = 7.98*MeV;
-    ExcitEnergies[20] = 6.90*MeV;
-    ExcitEnergies[25] = 5.83*MeV;
-    ExcitEnergies[26] = 8.57*MeV;
-    ExcitEnergies[31] = 5.33*MeV;
-
-    ExcitSpins[18] = 4;
-    ExcitSpins[20] = 6;
-    ExcitSpins[25] = 7;
-    ExcitSpins[26] = 4;
-    ExcitSpins[31] = 13;
-	
-    SetExcitationEnergiesPtr(&ExcitEnergies);
-    SetExcitationSpinsPtr(&ExcitSpins);		
+    G4EvaporationProbability(4,2,1,&theCoulombBarrier) // A,Z,Gamma,&theCoumlombBarrier
+{
+	
 }
 
@@ -93,5 +67,12 @@
 }
 
-G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ) const
+
+ G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)  
+  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
+	
+ G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &) 
+  { return 0.0; }
+
+G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ) 
 {
     // Data comes from 
@@ -116,2 +97,151 @@
     return C;
 }
+
+///////////////////////////////////////////////////////////////////////////////////
+//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
+//OPT=0 Dostrovski's parameterization
+//OPT=1,2 Chatterjee's paramaterization 
+//OPT=3,4 Kalbach's parameterization 
+// 
+G4double G4AlphaEvaporationProbability::CrossSection(const  G4Fragment & fragment, 
+						     const  G4double K)
+{
+  theA=GetA();
+  theZ=GetZ();
+  ResidualA=fragment.GetA()-theA;
+  ResidualZ=fragment.GetZ()-theZ; 
+  
+  ResidualAthrd=std::pow(ResidualA,0.33333);
+  FragmentA=fragment.GetA();
+  FragmentAthrd=std::pow(FragmentA,0.33333);
+  
+  
+  if (OPTxs==0) {std::ostringstream errOs;
+    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (Alpha's)!!"  
+	  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;}
+  
+  if( OPTxs==1 || OPTxs==2) return G4AlphaEvaporationProbability::GetOpt12( K);
+  else if (OPTxs==3 || OPTxs==4)  return G4AlphaEvaporationProbability::GetOpt34( K);
+  else{
+    std::ostringstream errOs;
+    errOs << "BAD Alpha CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;
+  }
+}
+
+
+//
+//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+//(fitting to cross section from Bechetti & Greenles OM potential)
+
+G4double G4AlphaEvaporationProbability::GetOpt12(const  G4double K)
+{
+// c     ** alpha from huizenga and igo
+  G4double Kc=K;
+  
+  // JMQ xsec is set constat above limit of validity
+  if (K>50) Kc=50;
+  
+  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
+  
+  G4double     p0 = 10.95;
+  G4double     p1 = -85.2;
+  G4double     p2 = 1146.;
+  G4double     landa0 = 0.0643;
+  G4double     landa1 = -13.96;
+  G4double     mu0 = 781.2;
+  G4double     mu1 = 0.29;
+  G4double     nu0 = -304.7;
+  G4double     nu1 = -470.0;
+  G4double     nu2 = -8.580;   
+  G4double     delta=1.2;          
+
+ 
+  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
+  p = p0 + p1/Ec + p2/(Ec*Ec);
+  landa = landa0*ResidualA + landa1;
+  mu = mu0*std::pow(ResidualA,mu1);
+  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+  q = landa - nu/(Ec*Ec) - 2*p*Ec;
+  r = mu + 2*nu/Ec + p*(Ec*Ec);
+  
+  ji=std::max(Kc,Ec);
+  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
+  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
+                 
+  if (xs <0.0) {xs=0.0;}
+  
+  return xs;
+  
+}
+
+// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
+G4double G4AlphaEvaporationProbability::GetOpt34(const  G4double K)
+// c     ** alpha from huizenga and igo
+{
+  
+  G4double landa, mu, nu, p , signor(1.),sig;
+  G4double ec,ecsq,xnulam,etest(0.),a; 
+  G4double b,ecut,cut,ecut2,geom,elab;
+
+
+  G4double     flow = 1.e-18;
+  G4double     spill= 1.e+18;
+
+  G4double       p0 = 10.95;
+  G4double     p1 = -85.2;
+  G4double     p2 = 1146.;
+  G4double     landa0 = 0.0643;
+  G4double     landa1 = -13.96;
+  G4double     mu0 = 781.2;
+  G4double     mu1 = 0.29;
+  G4double     nu0 = -304.7;
+  G4double     nu1 = -470.0;
+  G4double     nu2 = -8.580;        
+  
+  G4double      ra=1.20;
+  
+  //JMQ 13/02/09 increase of reduced radius to lower the barrier
+  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
+  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
+  ecsq = ec * ec;
+  p = p0 + p1/ec + p2/ecsq;
+  landa = landa0*ResidualA + landa1;
+  a = std::pow(ResidualA,mu1);
+  mu = mu0 * a;
+  nu = a* (nu0+nu1*ec+nu2*ecsq);  
+  xnulam = nu / landa;
+  if (xnulam > spill) xnulam=0.;
+  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
+  
+  a = -2.*p*ec + landa - nu/ecsq;
+  b = p*ecsq + mu + 2.*nu/ec;
+  ecut = 0.;
+  cut = a*a - 4.*p*b;
+  if (cut > 0.) ecut = std::sqrt(cut);
+  ecut = (ecut-a) / (p+p);
+  ecut2 = ecut;
+  if (cut < 0.) ecut2 = ecut - 2.;
+  elab = K * FragmentA / ResidualA;
+  sig = 0.;
+  
+  if (elab <= ec) { //start for E<Ec
+    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+  }           //end for E<Ec
+  else {           //start for E>Ec
+    sig = (landa*elab+mu+nu/elab) * signor;
+    geom = 0.;
+    if (xnulam < flow || elab < etest) return sig;
+    geom = std::sqrt(theA*K);
+    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
+    geom = 31.416 * geom * geom;
+    sig = std::max(geom,sig);
+  }           //end for E>Ec
+  return sig;
+  
+}
+
+//   ************************** end of cross sections ******************************* 
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationChannel.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationChannel.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4DeuteronEvaporationChannel.cc,v 1.4 2006/06/29 20:10:21 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc	(revision 962)
@@ -24,66 +24,23 @@
 // ********************************************************************
 //
-//
-// $Id: G4DeuteronEvaporationProbability.cc,v 1.4 2006/06/29 20:10:23 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
-// by V. Lara (Nov 1999)
-//
+// by V. Lara (Oct 1998)
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
 
 
 #include "G4DeuteronEvaporationProbability.hh"
 
+
 G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability() :
-    G4EvaporationProbability(2,1,6) // A,Z,Gamma
-{
-    std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
-    std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
-    ExcitEnergies.reserve(NumExcitedStatesEnergy);
-    ExcitSpins.reserve(NumExcitedStatesSpin);
-    ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
-    ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
-
-
-    ExcitEnergies[15] = 6.18*MeV;
-    ExcitEnergies[17] = 2.15*MeV;
-    ExcitEnergies[18] = 5.02*MeV;
-    ExcitEnergies[19] = 2.65*MeV;
-    ExcitEnergies[20] = 4.80*MeV;
-    ExcitEnergies[22] = 3.85*MeV;
-    ExcitEnergies[23] = 6.96*MeV;
-    ExcitEnergies[25] = 4.92*MeV;
-    ExcitEnergies[26] = 7.22*MeV;
-    ExcitEnergies[27] = 0.40*MeV;
-    ExcitEnergies[28] = 6.83*MeV;
-    ExcitEnergies[29] = 7.12*MeV;
-    ExcitEnergies[30] = 3.84*MeV;
-    ExcitEnergies[31] = 3.92*MeV;
-
-    ExcitSpins[15] = 1;
-    ExcitSpins[17] = 3;
-    ExcitSpins[18] = 4;
-    ExcitSpins[19] = 4;
-    ExcitSpins[20] = 4;
-    ExcitSpins[22] = 6;
-    ExcitSpins[23] = 6;
-    ExcitSpins[25] = 1;
-    ExcitSpins[26] = 10;
-    ExcitSpins[27] = 3;
-    ExcitSpins[28] = 10;
-    ExcitSpins[29] = 3;
-    ExcitSpins[30] = 6;
-    ExcitSpins[31] = 5;
-	
-    SetExcitationEnergiesPtr(&ExcitEnergies);
-    SetExcitationSpinsPtr(&ExcitSpins);	
-	
-}
+    G4EvaporationProbability(2,1,3,&theCoulombBarrier) // A,Z,Gamma (fixed JMQ)
+{    	}
 G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability(const G4DeuteronEvaporationProbability &) : G4EvaporationProbability()
 {
     throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationProbability::copy_constructor meant to not be accessable");
 }
-
-
 
 
@@ -101,4 +58,5 @@
 }
 
+
 G4bool G4DeuteronEvaporationProbability::operator!=(const G4DeuteronEvaporationProbability &) const
 {
@@ -106,5 +64,18 @@
 }
 
-G4double G4DeuteronEvaporationProbability::CCoeficient(const G4double aZ) const
+
+G4double G4DeuteronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
+{
+    return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));
+}
+
+
+G4double G4DeuteronEvaporationProbability::CalcBetaParam(const G4Fragment & ) 
+{
+    return 0.0;
+}
+
+
+G4double G4DeuteronEvaporationProbability::CCoeficient(const G4double aZ) 
 {
     // Data comes from 
@@ -127,2 +98,150 @@
 	
 }
+
+
+///////////////////////////////////////////////////////////////////////////////////
+//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
+//OPT=0 Dostrovski's parameterization
+//OPT=1,2 Chatterjee's paramaterization 
+//OPT=3,4 Kalbach's parameterization 
+// 
+G4double G4DeuteronEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
+{
+       theA=GetA();
+       theZ=GetZ();
+       ResidualA=fragment.GetA()-theA;
+       ResidualZ=fragment.GetZ()-theZ; 
+ 
+       ResidualAthrd=std::pow(ResidualA,0.33333);
+       FragmentA=fragment.GetA();
+       FragmentAthrd=std::pow(FragmentA,0.33333);
+
+       if (OPTxs==0) {std::ostringstream errOs;
+         errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (deuterons)!!"  <<G4endl;
+         throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+         return 0.;}
+       if( OPTxs==1 || OPTxs==2) return G4DeuteronEvaporationProbability::GetOpt12( K);
+       else if (OPTxs==3 || OPTxs==4)  return G4DeuteronEvaporationProbability::GetOpt34( K);
+       else{
+         std::ostringstream errOs;
+         errOs << "BAD Deuteron CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
+         throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+         return 0.;
+       }
+}
+
+
+//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+//(fitting to cross section from Bechetti & Greenles OM potential)
+
+G4double G4DeuteronEvaporationProbability::GetOpt12(const  G4double K)
+{
+
+  G4double Kc=K;
+
+// JMQ xsec is set constat above limit of validity
+ if (K>50) Kc=50;
+
+ G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
+
+ 
+ G4double    p0 = -38.21;
+ G4double    p1 = 922.6;
+ G4double    p2 = -2804.;
+ G4double    landa0 = -0.0323;
+ G4double    landa1 = -5.48;
+ G4double    mu0 = 336.1;
+ G4double    mu1 = 0.48;
+ G4double    nu0 = 524.3;
+ G4double    nu1 = -371.8;
+ G4double    nu2 = -5.924;  
+ G4double    delta=1.2;            
+ 
+
+ Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
+ p = p0 + p1/Ec + p2/(Ec*Ec);
+ landa = landa0*ResidualA + landa1;
+ mu = mu0*std::pow(ResidualA,mu1);
+ nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+ q = landa - nu/(Ec*Ec) - 2*p*Ec;
+ r = mu + 2*nu/Ec + p*(Ec*Ec);
+ 
+ ji=std::max(Kc,Ec);
+ if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
+ else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
+                 
+ if (xs <0.0) {xs=0.0;}
+ 
+ return xs;
+
+}
+
+
+// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
+G4double G4DeuteronEvaporationProbability::GetOpt34(const  G4double K)
+//     ** d from o.m. of perey and perey
+{
+
+  G4double landa, mu, nu, p , signor(1.),sig;
+  G4double ec,ecsq,xnulam,etest(0.),a; 
+  G4double b,ecut,cut,ecut2,geom,elab;
+
+
+  G4double     flow = 1.e-18;
+  G4double     spill= 1.e+18;
+
+
+  G4double     p0 = 0.798;
+  G4double     p1 = 420.3;
+  G4double     p2 = -1651.;
+  G4double     landa0 = 0.00619;
+  G4double     landa1 = -7.54;
+  G4double     mu0 = 583.5;
+  G4double     mu1 = 0.337;
+  G4double     nu0 = 421.8;
+  G4double     nu1 = -474.5;
+  G4double     nu2 = -3.592;      
+  
+  G4double      ra=0.80;
+        
+  //JMQ 13/02/09 increase of reduced radius to lower the barrier
+  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
+  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
+  ecsq = ec * ec;
+  p = p0 + p1/ec + p2/ecsq;
+  landa = landa0*ResidualA + landa1;
+  a = std::pow(ResidualA,mu1);
+  mu = mu0 * a;
+  nu = a* (nu0+nu1*ec+nu2*ecsq);  
+  xnulam = nu / landa;
+  if (xnulam > spill) xnulam=0.;
+  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
+
+  a = -2.*p*ec + landa - nu/ecsq;
+  b = p*ecsq + mu + 2.*nu/ec;
+  ecut = 0.;
+  cut = a*a - 4.*p*b;
+  if (cut > 0.) ecut = std::sqrt(cut);
+  ecut = (ecut-a) / (p+p);
+  ecut2 = ecut;
+  if (cut < 0.) ecut2 = ecut - 2.;
+  elab = K * FragmentA / ResidualA;
+  sig = 0.;
+
+  if (elab <= ec) { //start for E<Ec
+    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+  }           //end for E<Ec
+  else {           //start for E>Ec
+    sig = (landa*elab+mu+nu/elab) * signor;
+    geom = 0.;
+    if (xnulam < flow || elab < etest) return sig;
+    geom = std::sqrt(theA*K);
+    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
+    geom = 31.416 * geom * geom;
+    sig = std::max(geom,sig);
+  }           //end for E>Ec
+  return sig;
+  
+}
+
+//   ************************** end of cross sections ******************************* 
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc	(revision 962)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Evaporation.cc,v 1.7 2007/02/14 13:37:49 ahoward Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4Evaporation.cc,v 1.12 2008/12/09 17:57:36 ahoward Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -33,4 +33,9 @@
 // Alex Howard - added protection for negative probabilities in the sum, 14/2/07
 //
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
+// JMQ (06 September 2008) Also external choices have been added for 
+// superimposed Coulomb barrier (if useSICBis set true, by default is false) 
+
 #include "G4Evaporation.hh"
 #include "G4EvaporationFactory.hh"
@@ -89,5 +94,4 @@
 }
 
-
 G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
 {
@@ -105,13 +109,18 @@
     // Number of channels
     G4int TotNumberOfChannels = theChannels->size();  
-	
 
     // Starts loop over evaporated particles
     for (;;) 
-      {
+ 
+   {    
 	// loop over evaporation channels
 	std::vector<G4VEvaporationChannel*>::iterator i;
 	for (i=theChannels->begin(); i != theChannels->end(); i++) 
 	  {
+  // for inverse cross section choice
+            (*i)->SetOPTxs(OPTxs);
+  // for superimposed Coulomb Barrier for inverse cross sections
+            (*i)->UseSICB(useSICB);
+
 	    (*i)->Initialize(theResidualNucleus);
 	  }
@@ -172,4 +181,10 @@
 	    if( j >= TotNumberOfChannels ) 
 	      {
+		G4cerr << " Residual A: " << theResidualNucleus.GetA() << " Residual Z: " << theResidualNucleus.GetZ() << " Excitation Energy: " << theResidualNucleus.GetExcitationEnergy() << G4endl;
+		G4cerr << " j has not chosen a channel, j = " << j << " TotNumberOfChannels " << TotNumberOfChannels << " Total Probability: " << TotalProbability << G4endl;
+		for (j=0; j < TotNumberOfChannels; j++) 
+		  {
+		    G4cerr << " j: " << j << " EmissionProbChannel: " << EmissionProbChannel[j] << " and shoot: " << shoot << " (<ProbChannel?) " << G4endl;
+		  }		
 		throw G4HadronicException(__FILE__, __LINE__,  "G4Evaporation::BreakItUp: Can't define emission probability of the channels!" );
 	      } 
@@ -216,4 +231,5 @@
 		  theResidualNucleus.SetCreatorModel(G4String("ResidualNucleus"));
 #endif
+
 		}
 	      }
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc	(revision 962)
@@ -25,10 +25,14 @@
 //
 //
-// $Id: G4EvaporationChannel.cc,v 1.5 2006/06/29 20:10:27 gunter Exp $
-// GEANT4 tag $Name:  $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
 // by V. Lara (Oct 1998)
 //
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
+// JMQ (06 September 2008) Also external choices have been added for 
+// superimposed Coulomb barrier (if useSICB is set true, by default is false) 
+
 
 #include "G4EvaporationChannel.hh"
@@ -36,9 +40,11 @@
 
 
-G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ,
+
+G4EvaporationChannel::G4EvaporationChannel(const G4int anA, const G4int aZ, const G4String & aName,
 					   G4VEmissionProbability * aEmissionStrategy,
-					   G4VCoulombBarrier * aCoulombBarrier):
-    A(theA),
-    Z(theZ),
+                                           G4VCoulombBarrier * aCoulombBarrier):
+    G4VEvaporationChannel(aName),
+    theA(anA),
+    theZ(aZ),
     theEvaporationProbabilityPtr(aEmissionStrategy),
     theCoulombBarrierPtr(aCoulombBarrier),
@@ -47,45 +53,11 @@
 { 
     theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
-    MyOwnLevelDensity = true;
-}
-
-G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String & aName,
-					   G4VEmissionProbability * aEmissionStrategy,
-					   G4VCoulombBarrier * aCoulombBarrier):
-    G4VEvaporationChannel(aName),
-    A(theA),
-    Z(theZ),
-    theEvaporationProbabilityPtr(aEmissionStrategy),
-    theCoulombBarrierPtr(aCoulombBarrier),
-    EmissionProbability(0.0),
-    MaximalKineticEnergy(-1000.0)
-{ 
-    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
-    MyOwnLevelDensity = true;
-}
-
-G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String * aName,
-					   G4VEmissionProbability * aEmissionStrategy,
-					   G4VCoulombBarrier * aCoulombBarrier):
-    G4VEvaporationChannel(aName),
-    A(theA),
-    Z(theZ),
-    theEvaporationProbabilityPtr(aEmissionStrategy),
-    theCoulombBarrierPtr(aCoulombBarrier),
-    EmissionProbability(0.0),
-    MaximalKineticEnergy(-1000.0)
-{ 
-    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
-    MyOwnLevelDensity = true;
-}
-
+//    MyOwnLevelDensity = true;  
+}
 
 G4EvaporationChannel::~G4EvaporationChannel()
 {
-    if (MyOwnLevelDensity) delete theLevelDensityPtr;
-}
-
-
-
+  delete theLevelDensityPtr;
+}
 
 G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel()
@@ -112,154 +84,172 @@
 }
 
-
+//void G4EvaporationChannel::SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)
+//  {
+//    if (MyOwnLevelDensity) delete theLevelDensityPtr;
+//    theLevelDensityPtr = aLevelDensity;
+//    MyOwnLevelDensity = false;
+//  }
 
 void G4EvaporationChannel::Initialize(const G4Fragment & fragment)
 {
-
-    G4int anA = static_cast<G4int>(fragment.GetA());
-    G4int aZ = static_cast<G4int>(fragment.GetZ());
-    AResidual = anA - A;
-    ZResidual = aZ - Z;
-
-    // Effective excitation energy
-    G4double ExEnergy = fragment.GetExcitationEnergy() - 
-      G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
-
-    // We only take into account channels which are physically allowed
-    if (AResidual <= 0 || ZResidual <= 0 || AResidual < ZResidual ||
-	(AResidual == ZResidual && AResidual > 1) || ExEnergy <= 0.0) {
-	// 		LevelDensityParameter = 0.0;
-	CoulombBarrier = 0.0;
-	// 		BindingEnergy = 0.0;
-	MaximalKineticEnergy = -1000.0*MeV;
-	EmissionProbability = 0.0;
-    } else {
-	// // Get Level Density
-	// LevelDensityParameter = theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
-
-	// Coulomb Barrier calculation
-	CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(AResidual,ZResidual,ExEnergy);
-	
-	// // Binding Enegy (for separate fragment from nucleus)
-	// BindingEnergy = CalcBindingEnergy(anA,aZ);
-
-	// Maximal Kinetic Energy
-	MaximalKineticEnergy = CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()->
-							GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy);
-		
-	// Emission probability
-	if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0;
-	else { 
-	    // Total emission probability for this channel
-	    EmissionProbability = theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
-
-	}
+  //for inverse cross section choice
+  theEvaporationProbabilityPtr->SetOPTxs(OPTxs);
+  // for superimposed Coulomb Barrier for inverse cross sections
+  theEvaporationProbabilityPtr->UseSICB(useSICB);
+  
+  
+  G4int FragmentA = static_cast<G4int>(fragment.GetA());
+  G4int FragmentZ = static_cast<G4int>(fragment.GetZ());
+  ResidualA = FragmentA - theA;
+  ResidualZ = FragmentZ - theZ;
+  
+  //Effective excitation energy
+  G4double ExEnergy = fragment.GetExcitationEnergy() - 
+    G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ);
+  
+  
+  // Only channels which are physically allowed are taken into account 
+  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
+      (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) {
+    CoulombBarrier=0.0;
+    MaximalKineticEnergy = -1000.0*MeV;
+    EmissionProbability = 0.0;
+  } else {
+    CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
+    // Maximal Kinetic Energy
+    MaximalKineticEnergy = CalcMaximalKineticEnergy
+      (G4ParticleTable::GetParticleTable()->
+       GetIonTable()->GetNucleusMass(FragmentZ,FragmentA)+ExEnergy);
+    
+    // Emission probability
+    // Protection for the case Tmax<V. If not set in this way we could end up in an 
+    // infinite loop in  the method GetKineticEnergy if OPTxs!=0 && useSICB=true. 
+    // Of course for OPTxs=0 we have the Coulomb barrier 
+    
+    G4double limit;
+    if (OPTxs==0 || (OPTxs!=0 && useSICB)) 
+      limit= CoulombBarrier;
+    else limit=0.;
+  
+    // The threshold for charged particle emission must be  set to 0 if Coulomb 
+    //cutoff  is included in the cross sections
+    //if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0;  
+    if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0;
+    else { 
+      // Total emission probability for this channel
+      EmissionProbability = theEvaporationProbabilityPtr->
+        EmissionProbability(fragment, MaximalKineticEnergy);
     }
-	
-    return;
-}
-
+  }
+  
+  return;
+}
 
 G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus)
 {
-
-    G4double EvaporatedKineticEnergy = CalcKineticEnergy();
-    G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
-    G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
-
-    G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
-						(EvaporatedKineticEnergy+2.0*EvaporatedMass))));
-  
-    momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
-
-    G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
-    EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
-
-    G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
+  G4double EvaporatedKineticEnergy=GetKineticEnergy(theNucleus);
+  
+  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
+    GetIonMass(theZ,theA);
+  G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
+  
+  G4ThreeVector momentum(IsotropicVector
+                         (std::sqrt(EvaporatedKineticEnergy*
+                                    (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
+  
+  momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
+  
+  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
+  EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
+  
+  G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum);
 #ifdef PRECOMPOUND_TEST
-    EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
+  EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
 #endif
-    // ** And now the residual nucleus ** 
-    G4double theExEnergy = theNucleus.GetExcitationEnergy();
-    G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
-	GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
-		       static_cast<G4int>(theNucleus.GetA()));
-    G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
-	
-    G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
-    ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
-	
-    G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum );
+  // ** And now the residual nucleus ** 
+  G4double theExEnergy = theNucleus.GetExcitationEnergy();
+  G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
+    GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
+                   static_cast<G4int>(theNucleus.GetA()));
+  G4double ResidualEnergy = theMass + 
+    (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
+  
+  G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
+  ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
+  
+  G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum );
+  
 #ifdef PRECOMPOUND_TEST
-    ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
+  ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
 #endif
-    G4FragmentVector * theResult = new G4FragmentVector;
-
+  G4FragmentVector * theResult = new G4FragmentVector;
+  
 #ifdef debug
-    G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
-    G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
-    if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
-	G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
-	G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :" 
-	       <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV 
-	       << " MeV" << G4endl;
-    }
-    if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
-	std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
-	std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
-	G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
-	G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :" 
-	       <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
-	       << " MeV" << G4endl;   
-
-    }
+  G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
+  G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
+  if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
+    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
+    G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :" 
+           <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV 
+           << " MeV" << G4endl;
+  }
+  if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
+      std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
+      std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
+    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
+    G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :" 
+           <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
+           << " MeV" << G4endl;   
+    
+  }
 #endif
-    theResult->push_back(EvaporatedFragment);
-    theResult->push_back(ResidualFragment);
-    return theResult; 
+  theResult->push_back(EvaporatedFragment);
+  theResult->push_back(ResidualFragment);
+  return theResult; 
 } 
 
-
-
-// G4double G4EvaporationChannel::CalcBindingEnergy(const G4int anA, const G4int aZ)
-// // Calculate Binding Energy for separate fragment from nucleus
-// {
-// 	// Mass Excess for residual nucleus
-// 	G4double ResNucMassExcess = G4NucleiProperties::GetNuclearMass(AResidual,ZResidual);
-// 	// Mass Excess for fragment
-// 	G4double FragmentMassExcess = G4NucleiProperties::GetNuclearMass(A,Z);
-// 	// Mass Excess for Compound Nucleus
-// 	G4double NucleusMassExcess = G4NucleiProperties::GetNuclearMass(anA,aZ);
-// 
-// 	return ResNucMassExcess + FragmentMassExcess - NucleusMassExcess;
-// }
-
-
+/////////////////////////////////////////
+// Calculates the maximal kinetic energy that can be carried by fragment.
 G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
-    // Calculate maximal kinetic energy that can be carried by fragment.
-{
-    G4double ResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( ZResidual, AResidual );
-    G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( Z, A );
-	
-    G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
-	(2.0*NucleusTotalE) -
-	EvaporatedMass - CoulombBarrier;
-	
-    return T;
-}
-
-
-
-
-G4double G4EvaporationChannel::CalcKineticEnergy(void)
-    // Samples fragment kinetic energy.
+{
+  G4double ResidualMass = G4ParticleTable::GetParticleTable()->
+    GetIonTable()->GetNucleusMass( ResidualZ, ResidualA );
+  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
+    GetIonTable()->GetNucleusMass( theZ, theA );
+
+  // This is the "true" assimptotic kinetic energy (from energy conservation)	
+  G4double Tmax = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - 		
+                   ResidualMass*ResidualMass)/(2.0*NucleusTotalE) - EvaporatedMass;
+  
+  //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated
+  //at the Coulomb barrier
+  //IMPORTANT: meaning of Tmax differs in OPTxs=0 and OPTxs!=0
+  //When OPTxs!=0 Tmax is the TRUE (assimptotic) maximal kinetic energy
+  
+  if(OPTxs==0) 
+    Tmax=Tmax- CoulombBarrier;
+  
+  return Tmax;
+}
+
+///////////////////////////////////////////
+//JMQ: New method for MC sampling of kinetic energy. Substitutes old CalcKineticEnergy
+G4double G4EvaporationChannel::GetKineticEnergy(const G4Fragment & aFragment)
+{
+  
+  if (OPTxs==0) {
     // It uses Dostrovsky's approximation for the inverse reaction cross
-    // in the probability for fragment emisson
-{
-    if (MaximalKineticEnergy < 0.0) 
-	throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic energy is less than 0");
-
-    G4double Rb = 4.0*theLevelDensityPtr->LevelDensityParameter(AResidual+A,ZResidual+Z,MaximalKineticEnergy)*
-	MaximalKineticEnergy;
+    // in the probability for fragment emission
+    // MaximalKineticEnergy energy in the original version (V.Lara) was calculated at 
+    //the Coulomb barrier.
+    
+    
+    if (MaximalKineticEnergy < 0.0)   
+      throw G4HadronicException(__FILE__, __LINE__, 
+                                "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0");
+    
+    G4double Rb = 4.0*theLevelDensityPtr->
+      LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)*
+      MaximalKineticEnergy;
     G4double RbSqrt = std::sqrt(Rb);
     G4double PEX1 = 0.0;
@@ -268,24 +258,68 @@
     G4double FRk = 0.0;
     do {
-	G4double RandNumber = G4UniformRand();
-	Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
-	G4double Q1 = 1.0;
-	G4double Q2 = 1.0;
-	if (Z == 0) { // for emitted neutron
-	    G4double Beta = (2.12/std::pow(G4double(AResidual),2./3.) - 0.05)*MeV/
-		(0.76 + 2.2/std::pow(G4double(AResidual),1./3.));
-	    Q1 = 1.0 + Beta/(MaximalKineticEnergy);
-	    Q2 = Q1*std::sqrt(Q1);
-	} 
-    
-	FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
-    
+      G4double RandNumber = G4UniformRand();
+      Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
+      G4double Q1 = 1.0;
+      G4double Q2 = 1.0;
+      if (theZ == 0) { // for emitted neutron
+        G4double Beta = (2.12/std::pow(G4double(ResidualA),2./3.) - 0.05)*MeV/
+          (0.76 + 2.2/std::pow(G4double(ResidualA),1./3.));
+        Q1 = 1.0 + Beta/(MaximalKineticEnergy);
+        Q2 = Q1*std::sqrt(Q1);
+      } 
+      
+      FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
+      
     } while (FRk < G4UniformRand());
-
+    
     G4double result =  MaximalKineticEnergy * (1.0-Rk*Rk) + CoulombBarrier;
-
     return result;
-} 
- 
+  } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) {
+    
+    G4double V;
+    if(useSICB) V= CoulombBarrier;
+    else V=0.;
+    //If Coulomb barrier is just included  in the cross sections
+    // 	G4double V=0.;
+
+    G4double Tmax=MaximalKineticEnergy;
+    G4double T(0.0);
+    G4double NormalizedProbability(1.0);
+    
+    // A pointer is created in order to access the distribution function.
+    G4EvaporationProbability * G4EPtemp;
+    
+    if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability();
+    else if (theA==1 && theZ==1) G4EPtemp=new G4ProtonEvaporationProbability();
+    else if (theA==2 && theZ==1 ) G4EPtemp=new G4DeuteronEvaporationProbability();
+    else if (theA==3 && theZ==1 ) G4EPtemp=new G4TritonEvaporationProbability();
+    else if (theA==3 && theZ==2 ) G4EPtemp=new G4He3EvaporationProbability();
+    else if (theA==4 && theZ==2) G4EPtemp=new G4AlphaEvaporationProbability(); 
+    else {
+      std::ostringstream errOs;
+      errOs << "ejected particle out of range in G4EvaporationChannel"  << G4endl;
+      throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    }
+
+      //for cross section selection and superimposed Coulom Barrier for xs
+      G4EPtemp->SetOPTxs(OPTxs);
+      G4EPtemp->UseSICB(useSICB);
+
+    do
+      {  
+        T=V+G4UniformRand()*(Tmax-V);
+        NormalizedProbability=G4EPtemp->ProbabilityDistributionFunction(aFragment,T)/
+          (this->GetEmissionProbability());
+        
+      }
+    while (G4UniformRand() > NormalizedProbability);
+   delete G4EPtemp;
+   return T;
+  } else{
+    std::ostringstream errOs;
+    errOs << "Bad option for energy sampling in evaporation"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+  }
+}
 
 G4ThreeVector G4EvaporationChannel::IsotropicVector(const G4double Magnitude)
@@ -300,6 +334,13 @@
 			 Magnitude*CosTheta);
     return Vector;
-}
-
-
-
+            }
+
+
+
+
+
+   
+
+
+  
+
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4EvaporationFactory.cc,v 1.4 2006/06/29 20:10:29 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc	(revision 962)
@@ -24,12 +24,18 @@
 // ********************************************************************
 //
-//
-// $Id: G4EvaporationProbability.cc,v 1.5 2006/06/29 20:10:31 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
 // by V. Lara (Oct 1998)
 //
-
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
+// JMQ (06 September 2008) Also external choices have been added for 
+// superimposed Coulomb barrier (if useSICB is set true, by default is false) 
+//
+// JMQ (14 february 2009) bug fixed in emission width: hbarc instead of hbar_Planck in the denominator
+//
+#include <iostream>
+using namespace std;
 
 #include "G4EvaporationProbability.hh"
@@ -63,37 +69,28 @@
     return true;
 }
-
+  
 G4double G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, const G4double anEnergy)
 {
     G4double EmissionProbability = 0.0;
     G4double MaximalKineticEnergy = anEnergy;
+
     if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
-	EmissionProbability = CalcProbability(fragment,MaximalKineticEnergy);
-	// 		// Next there is a loop over excited states for this channel summing probabilities
-	// 		G4double SavedGamma = Gamma;
-	// 		for (G4int i = 0; i < ExcitationEnergies->length(); i++) {
-	// 			if (ExcitationSpins->operator()(i) < 0.1) continue;
-	// 			Gamma = ExcitationSpins->operator()(i)*A; // A is the channel mass number
-	// 			// substract excitation energies
-	// 			MaximalKineticEnergy -= ExcitationEnergies->operator()(i);
-	// 			// update probability
-	// 			EmissionProbability += CalcProbability(fragment,MaximalKineticEnergy);
-	// 				EmissionProbability += tmp;
-	// 			}
-	// 		// restore Gamma and MaximalKineticEnergy
-	// 		MaximalKineticEnergy = SavedMaximalKineticEnergy;
-	// 		Gamma = SavedGamma;
-	// 		}
+	EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy);
+
     }
     return EmissionProbability;
 }
 
-G4double G4EvaporationProbability::CalcProbability(const G4Fragment & fragment, 
-						   const G4double MaximalKineticEnergy)
-    // Calculate integrated probability (width) for rvaporation channel
-{	
-    G4double ResidualA = static_cast<G4double>(fragment.GetA() - theA);
-    G4double ResidualZ = static_cast<G4double>(fragment.GetZ() - theZ);
+////////////////////////////////////
+
+// Computes the integrated probability of evaporation channel
+G4double G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy)
+{
+    G4double ResidualA = fragment.GetA() - theA;
+    G4double ResidualZ = fragment.GetZ() - theZ;
     G4double U = fragment.GetExcitationEnergy();
+   
+ if (OPTxs==0) {
+
 	
     G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA);
@@ -107,5 +104,5 @@
 				      (U-delta0));
 								  
-    // compute the integrated probability of evaporation channel
+
     G4double RN = 1.5*fermi;
 
@@ -133,7 +130,130 @@
 	
     return Width;
-}
-
-
-
-
+             
+ } else if (OPTxs==1 || OPTxs==2 ||OPTxs==3 || OPTxs==4) {
+
+   G4double limit;
+   if (useSICB) limit=theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U);
+   else limit=0.;
+
+   if (MaximalKineticEnergy <= limit)  return 0.0;
+
+
+   // if Coulomb barrier cutoff is superimposed for all cross sections the limit is the Coulomb Barrier
+   G4double LowerLimit= limit;
+
+   //  MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel)     
+
+   G4double UpperLimit = MaximalKineticEnergy;
+
+
+   G4double Width = IntegrateEmissionProbability(fragment,LowerLimit,UpperLimit);
+
+   return Width;
+ } else{
+   std::ostringstream errOs;
+   errOs << "Bad option for cross sections at evaporation"  <<G4endl;
+   throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+ }
+  
+}
+
+/////////////////////////////////////////////////////////////////////
+
+G4double G4EvaporationProbability::
+IntegrateEmissionProbability(const G4Fragment & fragment, const G4double & Low, const G4double & Up )
+{
+
+  static const G4int N = 10;
+  // 10-Points Gauss-Legendre abcisas and weights
+  static const G4double w[N] = {
+    0.0666713443086881,
+    0.149451349150581,
+    0.219086362515982,
+    0.269266719309996,
+    0.295524224714753,
+    0.295524224714753,
+    0.269266719309996,
+    0.219086362515982,
+    0.149451349150581,
+    0.0666713443086881
+  };
+  static const G4double x[N] = {
+    -0.973906528517172,
+    -0.865063366688985,
+    -0.679409568299024,
+    -0.433395394129247,
+    -0.148874338981631,
+    0.148874338981631,
+    0.433395394129247,
+    0.679409568299024,
+    0.865063366688985,
+    0.973906528517172
+  };
+
+  G4double Total = 0.0;
+
+
+  for (G4int i = 0; i < N; i++) 
+    {
+
+      G4double KineticE = ((Up-Low)*x[i]+(Up+Low))/2.0;
+
+      Total += w[i]*ProbabilityDistributionFunction(fragment, KineticE);
+
+    }
+  Total *= (Up-Low)/2.0;
+  return Total;
+}
+
+
+/////////////////////////////////////////////////////////
+//New method (OPT=1,2,3,4)
+
+G4double G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment, const G4double K)
+{ 
+
+ 
+
+
+  G4double ResidualA = fragment.GetA() - theA;
+  G4double ResidualZ = fragment.GetZ() - theZ;  
+  G4double U = fragment.GetExcitationEnergy();
+
+  //        if(K <= theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return 0.0;   
+
+  G4double delta1 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ));
+
+ 
+  G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));
+
+  
+  G4double ParticleMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA);
+
+  G4double theSeparationEnergy= G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),static_cast<G4int>(theZ)) +
+    G4NucleiProperties::GetMassExcess(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ)) -
+    G4NucleiProperties::GetMassExcess(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));
+
+  G4double a0 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()),U - delta0);
+
+  G4double a1 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ),U - theSeparationEnergy - delta1);
+  
+  
+  G4double E0=U-delta0;
+
+  G4double E1=U-theSeparationEnergy-delta1-K;
+
+  if (E1<0.) return 0.;
+
+  //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead of hbar_Planck 
+  //Without 1/hbar_Panck remains as a width
+  //  G4double  Prob=Gamma*ParticleMass/((pi*hbar_Planck)*(pi*hbar_Planck)*
+  //std::exp(2*std::sqrt(a0*E0)))*K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn;
+
+  G4double Prob=Gamma*ParticleMass/((pi*hbarc)*(pi*hbarc)*std::exp(2*std::sqrt(a0*E0)))
+    *K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn;
+
+  return Prob;
+}
+
+
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationChannel.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationChannel.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4He3EvaporationChannel.cc,v 1.4 2006/06/29 20:10:33 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc	(revision 962)
@@ -24,41 +24,18 @@
 // ********************************************************************
 //
-//
-// $Id: G4He3EvaporationProbability.cc,v 1.4 2006/06/29 20:10:35 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
-// by V. Lara (Nov 1999)
-//
-
+// by V. Lara (Oct 1998)
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
 
 #include "G4He3EvaporationProbability.hh"
 
 G4He3EvaporationProbability::G4He3EvaporationProbability() :
-    G4EvaporationProbability(3,2,6) // A,Z,Gamma
-{
-    std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
-    std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
-    ExcitEnergies.reserve(NumExcitedStatesEnergy);
-    ExcitSpins.reserve(NumExcitedStatesSpin);
-    ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
-    ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
-
-
-	
-    ExcitEnergies[18] = 7.29*MeV;
-    ExcitEnergies[20] = 6.48*MeV;
-    ExcitEnergies[25] = 5.69*MeV;
-    ExcitEnergies[26] = 8.31*MeV;
-    ExcitEnergies[31] = 5.10*MeV;
-
-    ExcitSpins[18] = 6;
-    ExcitSpins[20] = 8;
-    ExcitSpins[25] = 3;
-    ExcitSpins[26] = 2;
-    ExcitSpins[31] = 7;	
-	
-    SetExcitationEnergiesPtr(&ExcitEnergies);
-    SetExcitationSpinsPtr(&ExcitSpins);	
+   G4EvaporationProbability(3,2,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
+{
+
 }
 
@@ -89,5 +66,12 @@
 }
 
-G4double G4He3EvaporationProbability::CCoeficient(const G4double aZ) const
+  G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)  
+  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
+	
+  G4double G4He3EvaporationProbability::CalcBetaParam(const G4Fragment & )  
+  { return 0.0; }
+
+
+G4double G4He3EvaporationProbability::CCoeficient(const G4double aZ) 
 {
     // Data comes from 
@@ -113,2 +97,148 @@
     return C*(4.0/3.0);
 }
+
+///////////////////////////////////////////////////////////////////////////////////
+//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
+//OPT=0 Dostrovski's parameterization
+//OPT=1,2 Chatterjee's paramaterization 
+//OPT=3,4 Kalbach's parameterization 
+// 
+G4double G4He3EvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
+{
+       theA=GetA();
+       theZ=GetZ();
+       ResidualA=fragment.GetA()-theA;
+       ResidualZ=fragment.GetZ()-theZ; 
+ 
+       ResidualAthrd=std::pow(ResidualA,0.33333);
+       FragmentA=fragment.GetA();
+       FragmentAthrd=std::pow(FragmentA,0.33333);
+
+
+       if (OPTxs==0) {std::ostringstream errOs;
+         errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (He3's)!!"  
+	       <<G4endl;
+         throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+         return 0.;}
+       if( OPTxs==1 || OPTxs==2) return G4He3EvaporationProbability::GetOpt12( K);
+       else if (OPTxs==3 || OPTxs==4)  return G4He3EvaporationProbability::GetOpt34( K);
+       else{
+         std::ostringstream errOs;
+         errOs << "BAD He3's CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
+         throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+         return 0.;
+       }
+}
+//
+//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+//(fitting to cross section from Bechetti & Greenles OM potential)
+
+G4double G4He3EvaporationProbability::GetOpt12(const  G4double K)
+{
+
+G4double Kc=K;
+
+// JMQ xsec is set constat above limit of validity
+ if (K>50) Kc=50;
+
+ G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
+
+ G4double     p0 = -3.06;
+ G4double     p1 = 278.5;
+ G4double     p2 = -1389.;
+ G4double     landa0 = -0.00535;
+ G4double     landa1 = -11.16;
+ G4double     mu0 = 555.5;
+ G4double     mu1 = 0.40;
+ G4double     nu0 = 687.4;
+ G4double     nu1 = -476.3;
+ G4double     nu2 = 0.509;    
+ G4double     delta=1.2;              
+
+      Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
+      p = p0 + p1/Ec + p2/(Ec*Ec);
+      landa = landa0*ResidualA + landa1;
+      mu = mu0*std::pow(ResidualA,mu1);
+      nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+      q = landa - nu/(Ec*Ec) - 2*p*Ec;
+      r = mu + 2*nu/Ec + p*(Ec*Ec);
+
+   ji=std::max(Kc,Ec);
+   if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
+   else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
+   
+   if (xs <0.0) {xs=0.0;}
+              
+   return xs;
+
+}
+
+// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
+G4double G4He3EvaporationProbability::GetOpt34(const  G4double K)
+//c     ** 3he from o.m. of gibson et al
+{
+  
+  G4double landa, mu, nu, p , signor(1.),sig;
+  G4double ec,ecsq,xnulam,etest(0.),a; 
+  G4double b,ecut,cut,ecut2,geom,elab;
+  
+  
+  G4double     flow = 1.e-18;
+  G4double     spill= 1.e+18;
+
+
+  G4double     p0 = -2.88;
+  G4double     p1 = 205.6;
+  G4double     p2 = -1487.;
+  G4double     landa0 = 0.00459;
+  G4double     landa1 = -8.93;
+  G4double     mu0 = 611.2;
+  G4double     mu1 = 0.35;
+  G4double     nu0 = 473.8;
+  G4double     nu1 = -468.2;
+  G4double     nu2 = -2.225;      
+  
+  G4double      ra=0.80;
+  
+  //JMQ 13/02/09 increase of reduced radius to lower the barrier
+  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
+  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
+  ecsq = ec * ec;
+  p = p0 + p1/ec + p2/ecsq;
+  landa = landa0*ResidualA + landa1;
+  a = std::pow(ResidualA,mu1);
+  mu = mu0 * a;
+  nu = a* (nu0+nu1*ec+nu2*ecsq);  
+  xnulam = nu / landa;
+  if (xnulam > spill) xnulam=0.;
+  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
+  
+  a = -2.*p*ec + landa - nu/ecsq;
+  b = p*ecsq + mu + 2.*nu/ec;
+  ecut = 0.;
+  cut = a*a - 4.*p*b;
+  if (cut > 0.) ecut = std::sqrt(cut);
+  ecut = (ecut-a) / (p+p);
+  ecut2 = ecut;
+  if (cut < 0.) ecut2 = ecut - 2.;
+  elab = K * FragmentA / ResidualA;
+  sig = 0.;
+
+  if (elab <= ec) { //start for E<Ec
+    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+  }           //end for E<Ec
+  else {           //start for E>Ec
+    sig = (landa*elab+mu+nu/elab) * signor;
+    geom = 0.;
+    if (xnulam < flow || elab < etest) return sig;
+    geom = std::sqrt(theA*K);
+    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
+    geom = 31.416 * geom * geom;
+    sig = std::max(geom,sig);
+  }           //end for E>Ec
+  return sig;
+  
+}
+
+//   ************************** end of cross sections ******************************* 
+
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationChannel.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationChannel.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4NeutronEvaporationChannel.cc,v 1.4 2006/06/29 20:10:37 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationProbability.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationProbability.cc	(revision 962)
@@ -24,69 +24,18 @@
 // ********************************************************************
 //
-//
-// $Id: G4NeutronEvaporationProbability.cc,v 1.4 2006/06/29 20:10:39 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
-// by V. Lara (Nov 1999)
-//
-
+// by V. Lara (Oct 1998)
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
 
 #include "G4NeutronEvaporationProbability.hh"
 
 G4NeutronEvaporationProbability::G4NeutronEvaporationProbability() :
-    G4EvaporationProbability(1,0,2) // A,Z,Gamma
-{
-    std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
-    std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
-    ExcitEnergies.reserve(NumExcitedStatesEnergy);
-    ExcitSpins.reserve(NumExcitedStatesSpin);
-    ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
-    ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
-
-	
-    ExcitEnergies[ 9] = 3.56*MeV;
-    ExcitEnergies[10] = 0.48*MeV;
-    ExcitEnergies[11] = 0.98*MeV;
-    ExcitEnergies[12] = 0.43*MeV;
-    ExcitEnergies[15] = 3.37*MeV;
-    ExcitEnergies[17] = 0.72*MeV;
-    ExcitEnergies[18] = 2.13*MeV;
-    ExcitEnergies[19] = 0.95*MeV;
-    ExcitEnergies[20] = 2.00*MeV;
-    ExcitEnergies[21] = 4.44*MeV;
-    ExcitEnergies[22] = 3.09*MeV;
-    ExcitEnergies[23] = 6.09*MeV;
-    ExcitEnergies[25] = 2.31*MeV;
-    ExcitEnergies[26] = 5.28*MeV;
-    ExcitEnergies[27] = 0.12*MeV;
-    ExcitEnergies[28] = 5.22*MeV;
-    ExcitEnergies[29] = 6.10*MeV;
-    ExcitEnergies[30] = 0.87*MeV;
-    ExcitEnergies[31] = 1.98*MeV;
-
-
-    ExcitSpins[ 9] = 1;
-    ExcitSpins[10] = 2;
-    ExcitSpins[11] = 3;
-    ExcitSpins[12] = 2;
-    ExcitSpins[15] = 5;
-    ExcitSpins[17] = 3;
-    ExcitSpins[18] = 2;
-    ExcitSpins[19] = 5;
-    ExcitSpins[20] = 2;
-    ExcitSpins[21] = 5;
-    ExcitSpins[22] = 2;
-    ExcitSpins[23] = 3;
-    ExcitSpins[25] = 1;
-    ExcitSpins[26] = 8;
-    ExcitSpins[27] = 1;
-    ExcitSpins[28] = 8;
-    ExcitSpins[29] = 8;
-    ExcitSpins[30] = 2;
-    ExcitSpins[31] = 5;
-
-    SetExcitationEnergiesPtr(&ExcitEnergies);
-    SetExcitationSpinsPtr(&ExcitSpins);
+    G4EvaporationProbability(1,0,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
+{
+  
 }
 
@@ -117,2 +66,170 @@
     return true;
 }
+
+ G4double G4NeutronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
+  { return 0.76+2.2/std::pow(static_cast<G4double>(fragment.GetA()-GetA()),1.0/3.0);}
+
+	
+  G4double G4NeutronEvaporationProbability::CalcBetaParam(const G4Fragment &  fragment) 
+  { return (2.12/std::pow(static_cast<G4double>(fragment.GetA()-GetA()),2.0/3.0) - 0.05)*MeV/
+      CalcAlphaParam(fragment); }
+
+
+////////////////////////////////////////////////////////////////////////////////////
+//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
+//OPT=0 Dostrovski's parameterization
+//OPT=1,2 Chatterjee's paramaterization 
+//OPT=3,4 Kalbach's parameterization 
+// 
+ G4double G4NeutronEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
+{ 
+  theA=GetA();
+  theZ=GetZ();
+  ResidualA=fragment.GetA()-theA;
+  ResidualZ=fragment.GetZ()-theZ; 
+ 
+  ResidualAthrd=std::pow(ResidualA,0.33333);
+  FragmentA=fragment.GetA();
+  FragmentAthrd=std::pow(FragmentA,0.33333);
+
+
+  if (OPTxs==0) {std::ostringstream errOs;
+    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (neutrons)!!"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;}
+  else if( OPTxs==1 ||OPTxs==2) return GetOpt12( K);
+  else if (OPTxs==3 || OPTxs==4)  return GetOpt34( K);
+  else{
+    std::ostringstream errOs;
+    errOs << "BAD NEUTRON CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;
+  }
+}
+
+
+
+
+
+
+//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+//(fitting to cross section from Bechetti & Greenles OM potential)
+
+G4double G4NeutronEvaporationProbability::GetOpt12(const  G4double K)
+{
+
+  G4double Kc=K;
+
+// JMQ  xsec is set constat above limit of validity
+  if (K>50) Kc=50;
+
+  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
+
+  landa0 = 18.57;
+  landa1 = -22.93;
+  mu0 = 381.7;
+  mu1 = 24.31;
+  nu0 = 0.172;
+  nu1 = -15.39;
+  nu2 = 804.8;
+  landa = landa0/ResidualAthrd + landa1;
+  mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd;
+  nu = nu0*ResidualAthrd*ResidualA + nu1*ResidualAthrd*ResidualAthrd + nu2 ;
+  xs=landa*Kc + mu + nu/Kc;
+  if (xs <= 0.0 ){
+    std::ostringstream errOs;
+    G4cout<<"WARNING:  NEGATIVE OPT=1 neutron cross section "<<G4endl;     
+    errOs << "RESIDUAL: Ar=" << ResidualA << " Zr=" << ResidualZ <<G4endl;
+    errOs <<"  xsec("<<Kc<<" MeV) ="<<xs <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+  }
+  return xs;
+}
+
+
+// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
+G4double G4NeutronEvaporationProbability::GetOpt34(const  G4double K)
+{
+
+  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
+  G4double p, p0, p1, p2;
+  G4double flow,spill,ec,ecsq,xnulam,etest(0.),ra(0.),a,signor(1.),sig; 
+  G4double b,ecut,cut,ecut2,geom,elab;
+
+//safety initialization
+  landa0=0;
+  landa1=0;
+  mu0=0.;
+  mu1=0.;
+  nu0=0.;
+  nu1=0.;
+  nu2=0.;
+  p0=0.;
+  p1=0.;
+  p2=0.;
+
+
+  flow = 1.e-18;
+  spill= 1.0e+18; 
+
+  
+
+// PRECO xs for neutrons is choosen
+
+  p0 = -312.;
+  p1= 0.;
+  p2 = 0.;
+  landa0 = 12.10;
+  landa1=  -11.27;
+  mu0 = 234.1;
+  mu1 = 38.26;
+  nu0 = 1.55;
+  nu1 = -106.1;
+  nu2 = 1280.8; 
+
+  if (ResidualA < 40.) signor=0.7+ResidualA*0.0075;
+  if (ResidualA > 210.) signor = 1. + (ResidualA-210.)/250.;
+  landa = landa0/ResidualAthrd + landa1;
+  mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd;
+  nu = nu0*ResidualAthrd*ResidualA + nu1*ResidualAthrd*ResidualAthrd + nu2;
+
+  // JMQ very low energy behaviour corrected (problem  for A (apprx.)>60)
+  if (nu < 0.)nu=-nu;
+
+  ec = 0.5;
+  ecsq = 0.25;
+  p = p0;
+  xnulam = 1.;
+  etest = 32.;
+//          ** etest is the energy above which the rxn cross section is
+//          ** compared with the geometrical limit and the max taken.
+//          ** xnulam here is a dummy value to be used later.
+
+
+
+  a = -2.*p*ec + landa - nu/ecsq;
+  b = p*ecsq + mu + 2.*nu/ec;
+  ecut = 0.;
+  cut = a*a - 4.*p*b;
+  if (cut > 0.) ecut = std::sqrt(cut);
+  ecut = (ecut-a) / (p+p);
+  ecut2 = ecut;
+  if (cut < 0.) ecut2 = ecut - 2.;
+  elab = K * FragmentA / ResidualA;
+  sig = 0.;
+  if (elab <= ec) { //start for E<Ec 
+    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;    
+  }              //end for E<Ec
+  else {           //start for  E>Ec
+    sig = (landa*elab+mu+nu/elab) * signor;
+    geom = 0.;
+    if (xnulam < flow || elab < etest) return sig;
+    geom = std::sqrt(theA*K);
+    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
+    geom = 31.416 * geom * geom;
+    sig = std::max(geom,sig);
+  }
+  return sig;}
+
+//   ************************** end of cross sections ******************************* 
+
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationChannel.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationChannel.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4ProtonEvaporationChannel.cc,v 1.4 2006/06/29 20:10:41 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc	(revision 962)
@@ -24,61 +24,18 @@
 // ********************************************************************
 //
-//
-// $Id: G4ProtonEvaporationProbability.cc,v 1.4 2006/06/29 20:10:43 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
-// by V. Lara (Nov 1999)
-//
-
+// by V. Lara (Oct 1998)
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
 
 #include "G4ProtonEvaporationProbability.hh"
 
 G4ProtonEvaporationProbability::G4ProtonEvaporationProbability() :
-    G4EvaporationProbability(1,1,2) // A,Z,Gamma
-{
-    std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
-    std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
-    ExcitEnergies.reserve(NumExcitedStatesEnergy);
-    ExcitSpins.reserve(NumExcitedStatesSpin);
-    ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
-    ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
-
-
-
-    ExcitEnergies[15] = 5.96*MeV;
-    ExcitEnergies[17] = 1.74*MeV;
-    ExcitEnergies[18] = 4.44*MeV;
-    ExcitEnergies[19] = 1.67*MeV;
-    ExcitEnergies[20] = 4.32*MeV;
-    ExcitEnergies[22] = 3.68*MeV;
-    ExcitEnergies[23] = 6.69*MeV;
-    ExcitEnergies[25] = 3.95*MeV;
-    ExcitEnergies[26] = 6.32*MeV;
-    ExcitEnergies[27] = 0.30*MeV;
-    ExcitEnergies[28] = 6.18*MeV;
-    ExcitEnergies[29] = 6.92*MeV;
-    ExcitEnergies[30] = 3.06*MeV;
-    ExcitEnergies[31] = 3.57*MeV;
-
-
-    ExcitSpins[15] = 8;
-    ExcitSpins[17] = 1;
-    ExcitSpins[18] = 6;
-    ExcitSpins[19] = 5;
-    ExcitSpins[20] = 6;
-    ExcitSpins[22] = 4;
-    ExcitSpins[23] = 8;
-    ExcitSpins[25] = 3;
-    ExcitSpins[26] = 4;
-    ExcitSpins[27] = 7;
-    ExcitSpins[28] = 4;
-    ExcitSpins[29] = 5;
-    ExcitSpins[30] = 2;
-    ExcitSpins[31] = 10;
-	
-    SetExcitationEnergiesPtr(&ExcitEnergies);
-    SetExcitationSpinsPtr(&ExcitSpins);	
-	
+    G4EvaporationProbability(1,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
+{
+   
 }
 
@@ -106,5 +63,11 @@
 }
 
-G4double G4ProtonEvaporationProbability::CCoeficient(const G4double aZ) const
+  G4double G4ProtonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
+  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
+	
+  G4double G4ProtonEvaporationProbability::CalcBetaParam(const G4Fragment & )  
+  { return 0.0; }
+
+  G4double G4ProtonEvaporationProbability::CCoeficient(const G4double aZ) 
 {
     // Data comes from 
@@ -126,2 +89,232 @@
 	
 }
+
+///////////////////////////////////////////////////////////////////////////////////
+//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
+//OPT=0 Dostrovski's parameterization
+//OPT=1,2 Chatterjee's paramaterization 
+//OPT=3,4 Kalbach's parameterization 
+// 
+G4double G4ProtonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
+{
+//  G4cout<<" In G4ProtonEVaporationProbability OPTxs="<<OPTxs<<G4endl;
+//  G4cout<<" In G4ProtonEVaporationProbability useSICB="<<useSICB<<G4endl;
+
+  theA=GetA();
+  theZ=GetZ();
+  ResidualA=fragment.GetA()-theA;
+  ResidualZ=fragment.GetZ()-theZ; 
+ 
+  ResidualAthrd=std::pow(ResidualA,0.33333);
+  FragmentA=fragment.GetA();
+  FragmentAthrd=std::pow(FragmentA,0.33333);
+  U=fragment.GetExcitationEnergy();
+
+
+  if (OPTxs==0) {std::ostringstream errOs;
+    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (protons)!!"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;}
+  else if( OPTxs==1 ) return GetOpt1( K);
+  else if( OPTxs==2 ||OPTxs==4) return GetOpt2( K);
+  else if (OPTxs==3 )  return GetOpt3( K);
+  else{
+    std::ostringstream errOs;
+    errOs << "BAD PROTON CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;
+  }
+}
+//********************* OPT=1 : Chatterjee's cross section ************************ 
+//(fitting to cross section from Bechetti & Greenles OM potential)
+
+G4double G4ProtonEvaporationProbability::GetOpt1(const  G4double K)
+{
+  G4double Kc=K; 
+
+// JMQ  xsec is set constat above limit of validity
+  if (K>50)  Kc=50;
+
+  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
+  G4double p, p0, p1, p2,Ec,delta,q,r,ji;
+  
+  p0 = 15.72;
+  p1 = 9.65;
+  p2 = -449.0;
+  landa0 = 0.00437;
+  landa1 = -16.58;
+  mu0 = 244.7;
+  mu1 = 0.503;
+  nu0 = 273.1;
+  nu1 = -182.4;
+  nu2 = -1.872;  
+  delta=0.;  
+
+  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
+  p = p0 + p1/Ec + p2/(Ec*Ec);
+  landa = landa0*ResidualA + landa1;
+  mu = mu0*std::pow(ResidualA,mu1);
+  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+  q = landa - nu/(Ec*Ec) - 2*p*Ec;
+  r = mu + 2*nu/Ec + p*(Ec*Ec);
+
+  ji=std::max(Kc,Ec);
+  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
+  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
+   if (xs <0.0) {xs=0.0;}
+
+   return xs; 
+
+}
+
+
+
+//************* OPT=2 : Wellisch's proton reaction cross section ************************ 
+
+G4double G4ProtonEvaporationProbability::GetOpt2(const  G4double K)
+{
+
+  G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
+ 
+//This is redundant when the Coulomb  barrier is overimposed to all cross sections 
+//It should be kept when Coulomb barrier only imposed at OPTxs=2, this is why ..
+
+         if(!useSICB && K <= theCoulombBarrier.GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return xine_th=0.0;
+
+  eekin=K;
+  rnpro=ResidualZ;
+  rnneu=ResidualA-ResidualZ;
+  ekin=eekin/1000;
+
+  r0=1.36*1.e-15;
+  fac=pi*r0*r0;
+  b0=2.247-0.915*(1.-1./ResidualAthrd);
+  fac1=b0*(1.-1./ResidualAthrd);
+  fac2=1.;
+  if(rnneu > 1.5) fac2=std::log(rnneu);
+  xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1);
+  xine_th=(1.-0.15*std::exp(-ekin))*xine_th/(1.00-0.0007*ResidualA);	
+  ff1=0.70-0.0020*ResidualA ;
+  ff2=1.00+1/ResidualA;
+  ff3=0.8+18/ResidualA-0.002*ResidualA;
+  fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2))));
+  xine_th=xine_th*(1.+ff3*fac);
+  ff1=1.-1/ResidualA-0.001*ResidualA;
+  ff2=1.17-2.7/ResidualA-0.0014*ResidualA;
+  fac=-8.*ff1*(std::log10(ekin)+2.0*ff2);
+  fac=1./(1.+std::exp(fac));
+  xine_th=xine_th*fac;               
+  if (xine_th < 0.0){
+    std::ostringstream errOs;
+    G4cout<<"WARNING:  negative Wellisch cross section "<<G4endl; 
+    errOs << "RESIDUAL: A=" << ResidualA << " Z=" << ResidualZ <<G4endl;
+    errOs <<"  xsec("<<ekin<<" MeV) ="<<xine_th <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+  }
+
+  return xine_th;
+            
+}
+
+
+// *********** OPT=3 : Kalbach's cross sections (from PRECO code)*************
+G4double G4ProtonEvaporationProbability::GetOpt3(const  G4double K)
+{
+//     ** p from  becchetti and greenlees (but modified with sub-barrier
+//     ** correction function and xp2 changed from -449)
+// JMQ (june 2008) : refinement of proton cross section for light systems
+//
+  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
+  G4double p, p0, p1, p2;
+  p0 = 15.72;
+  p1 = 9.65;
+  p2 = -300.;
+  landa0 = 0.00437;
+  landa1 = -16.58;
+  mu0 = 244.7;
+  mu1 = 0.503;
+  nu0 = 273.1;
+  nu1 = -182.4;
+  nu2 = -1.872;
+
+// parameters for  proton cross section refinement 
+  G4double afit,bfit,a2,b2;
+  afit=-0.0785656;
+  bfit=5.10789;
+  a2= -0.00089076;
+  b2= 0.0231597;  
+
+  G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig; 
+  G4double b,ecut,cut,ecut2,geom,elab;
+
+
+  G4double	flow = 1.e-18;
+  G4double       spill= 1.e+18; 
+
+
+  if (ResidualA <= 60.)  signor = 0.92;
+  else if (ResidualA < 100.) signor = 0.8 + ResidualA*0.002;
+
+
+  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
+  ecsq = ec * ec;
+  p = p0 + p1/ec + p2/ecsq;
+  landa = landa0*ResidualA + landa1;
+  a = std::pow(ResidualA,mu1);
+  mu = mu0 * a;
+  nu = a* (nu0+nu1*ec+nu2*ecsq);
+ 
+  c =std::min(3.15,ec*0.5);
+  w = 0.7 * c / 3.15; 
+
+  xnulam = nu / landa;
+  if (xnulam > spill) xnulam=0.;
+  if (xnulam >= flow) etest =std::sqrt(xnulam) + 7.;
+
+  a = -2.*p*ec + landa - nu/ecsq;
+  b = p*ecsq + mu + 2.*nu/ec;
+  ecut = 0.;
+  cut = a*a - 4.*p*b;
+  if (cut > 0.) ecut = std::sqrt(cut);
+  ecut = (ecut-a) / (p+p);
+  ecut2 = ecut;
+  if (cut < 0.) ecut2 = ecut - 2.;
+  elab = K * FragmentA / ResidualA;
+  sig = 0.;
+  if (elab <= ec) { //start for E<Ec 
+    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+    signor2 = (ec-elab-c) / w;
+    signor2 = 1. + std::exp(signor2);
+    sig = sig / signor2;
+    // signor2 is empirical global corr'n at low elab for protons in PRECO, not enough for light targets
+    //  refinement for proton cross section
+    if (ResidualZ<=26)
+      sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));      
+                       }              //end for E<Ec
+  else            {           //start for  E>Ec
+    sig = (landa*elab+mu+nu/elab) * signor;
+
+    //  refinement for proton cross section
+    if ( ResidualZ<=26 && elab <=(afit*ResidualZ+bfit)*ec) 
+           sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));
+                                                                               
+//
+    geom = 0.;
+
+    if (xnulam < flow || elab < etest)
+     {
+        if (sig <0.0) {sig=0.0;}
+        return sig;
+      }
+    geom = std::sqrt(theA*K);
+    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
+    geom = 31.416 * geom * geom;
+    sig = std::max(geom,sig);
+
+  }   //end for E>Ec
+ return sig;}
+
+
+
+//   ************************** end of cross sections ******************************* 
+
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationChannel.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationChannel.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4TritonEvaporationChannel.cc,v 1.4 2006/06/29 20:10:45 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc	(revision 962)
@@ -24,48 +24,18 @@
 // ********************************************************************
 //
-//
-// $Id: G4TritonEvaporationProbability.cc,v 1.4 2006/06/29 20:10:47 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+//J.M. Quesada (August2008). Based on:
 //
 // Hadronic Process: Nuclear De-excitations
-// by V. Lara (Nov 1999)
-//
-
+// by V. Lara (Oct 1998)
+//
+// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
+// cross section option
 
 #include "G4TritonEvaporationProbability.hh"
 
 G4TritonEvaporationProbability::G4TritonEvaporationProbability() :
-    G4EvaporationProbability(3,1,6) // A,Z,Gamma
-{
-    std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1;
-    std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1;
-    ExcitEnergies.reserve(NumExcitedStatesEnergy);
-    ExcitSpins.reserve(NumExcitedStatesSpin);
-    ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0);
-    ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0);
-
-	
-    ExcitEnergies[15] = 6.26*MeV;
-    ExcitEnergies[17] = 3.59*MeV;
-    ExcitEnergies[18] = 6.76*MeV;
-    ExcitEnergies[20] = 6.34*MeV;
-    ExcitEnergies[23] = 7.34*MeV;
-    ExcitEnergies[25] = 5.11*MeV;
-    ExcitEnergies[26] = 7.57*MeV;
-    ExcitEnergies[28] = 7.28*MeV;
-    ExcitEnergies[31] = 4.46*MeV;
-
-    ExcitSpins[15] = 5;
-    ExcitSpins[17] = 5;
-    ExcitSpins[18] = 10;
-    ExcitSpins[20] = 2;
-    ExcitSpins[23] = 5;
-    ExcitSpins[25] = 5;
-    ExcitSpins[26] = 8;
-    ExcitSpins[28] = 8;
-    ExcitSpins[31] = 3;
-
-    SetExcitationEnergiesPtr(&ExcitEnergies);
-    SetExcitationSpinsPtr(&ExcitSpins);	
+    G4EvaporationProbability(3,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
+{
+
 }
 
@@ -96,5 +66,11 @@
 }
 
-G4double G4TritonEvaporationProbability::CCoeficient(const G4double aZ) const
+   G4double G4TritonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
+  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
+	
+  G4double G4TritonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 
+  { return 0.0; }
+
+G4double G4TritonEvaporationProbability::CCoeficient(const G4double aZ) 
 {
     // Data comes from 
@@ -117,2 +93,150 @@
 	
 }
+
+///////////////////////////////////////////////////////////////////////////////////
+//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 
+//OPT=0 Dostrovski's parameterization
+//OPT=1,2 Chatterjee's paramaterization 
+//OPT=3,4 Kalbach's parameterization 
+// 
+G4double G4TritonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
+{
+  theA=GetA();
+  theZ=GetZ();
+  ResidualA=fragment.GetA()-theA;
+  ResidualZ=fragment.GetZ()-theZ; 
+  
+  ResidualAthrd=std::pow(ResidualA,0.33333);
+  FragmentA=fragment.GetA();
+  FragmentAthrd=std::pow(FragmentA,0.33333);
+
+  
+  if (OPTxs==0) {std::ostringstream errOs;
+    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (tritons)!!"  
+	  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;}
+  if( OPTxs==1 || OPTxs==2) return G4TritonEvaporationProbability::GetOpt12( K);
+  else if (OPTxs==3 || OPTxs==4)  return G4TritonEvaporationProbability::GetOpt34( K);
+  else{
+    std::ostringstream errOs;
+    errOs << "BAD Triton CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
+    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
+    return 0.;
+  }
+}
+
+//
+//********************* OPT=1,2 : Chatterjee's cross section ************************ 
+//(fitting to cross section from Bechetti & Greenles OM potential)
+
+G4double G4TritonEvaporationProbability::GetOpt12(const  G4double K)
+{
+
+  G4double Kc=K;
+
+// JMQ xsec is set constat above limit of validity
+  if (K>50) Kc=50;
+
+ G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
+//G4double Eo(0),epsilon1(0),epsilon2(0),discri(0);
+
+ 
+ G4double    p0 = -11.04;
+ G4double    p1 = 619.1;
+ G4double    p2 = -2147.;
+ G4double    landa0 = -0.0426;
+ G4double    landa1 = -10.33;
+ G4double    mu0 = 601.9;
+ G4double    mu1 = 0.37;
+ G4double    nu0 = 583.0;
+ G4double    nu1 = -546.2;
+ G4double    nu2 = 1.718;  
+ G4double    delta=1.2;            
+
+ Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
+ p = p0 + p1/Ec + p2/(Ec*Ec);
+ landa = landa0*ResidualA + landa1;
+ mu = mu0*std::pow(ResidualA,mu1);
+ nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
+ q = landa - nu/(Ec*Ec) - 2*p*Ec;
+ r = mu + 2*nu/Ec + p*(Ec*Ec);
+
+ ji=std::max(Kc,Ec);
+ if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
+ else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
+ 
+ if (xs <0.0) {xs=0.0;}
+ 
+ return xs;
+
+}
+
+// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
+G4double G4TritonEvaporationProbability::GetOpt34(const  G4double K)
+//     ** t from o.m. of hafele, flynn et al
+{
+  
+  G4double landa, mu, nu, p , signor(1.),sig;
+  G4double ec,ecsq,xnulam,etest(0.),a; 
+  G4double b,ecut,cut,ecut2,geom,elab;
+  
+  
+  G4double     flow = 1.e-18;
+  G4double     spill= 1.e+18;
+
+  G4double     p0 = -21.45;
+  G4double     p1 = 484.7;
+  G4double     p2 = -1608.;
+  G4double     landa0 = 0.0186;
+  G4double     landa1 = -8.90;
+  G4double     mu0 = 686.3;
+  G4double     mu1 = 0.325;
+  G4double     nu0 = 368.9;
+  G4double     nu1 = -522.2;
+  G4double     nu2 = -4.998;  
+  
+  G4double      ra=0.80;
+        
+  //JMQ 13/02/09 increase of reduced radius to lower the barrier
+  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
+  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
+  ecsq = ec * ec;
+  p = p0 + p1/ec + p2/ecsq;
+  landa = landa0*ResidualA + landa1;
+  a = std::pow(ResidualA,mu1);
+  mu = mu0 * a;
+  nu = a* (nu0+nu1*ec+nu2*ecsq);  
+  xnulam = nu / landa;
+  if (xnulam > spill) xnulam=0.;
+  if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
+  
+  a = -2.*p*ec + landa - nu/ecsq;
+  b = p*ecsq + mu + 2.*nu/ec;
+  ecut = 0.;
+  cut = a*a - 4.*p*b;
+  if (cut > 0.) ecut = std::sqrt(cut);
+  ecut = (ecut-a) / (p+p);
+  ecut2 = ecut;
+  if (cut < 0.) ecut2 = ecut - 2.;
+  elab = K * FragmentA / ResidualA;
+  sig = 0.;
+
+  if (elab <= ec) { //start for E<Ec
+    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
+  }           //end for E<Ec
+  else {           //start for E>Ec  
+    sig = (landa*elab+mu+nu/elab) * signor;
+    geom = 0.;
+    if (xnulam < flow || elab < etest) return sig;
+    geom = std::sqrt(theA*K);
+    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
+    geom = 31.416 * geom * geom;
+    sig = std::max(geom,sig);
+  }           //end for E>Ec
+  return sig;
+  
+}
+
+//   ************************** end of cross sections ******************************* 
+
Index: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc
===================================================================
--- trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc	(revision 819)
+++ trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc	(revision 962)
@@ -26,5 +26,5 @@
 //
 // $Id: G4VEvaporation.cc,v 1.5 2006/06/29 20:10:49 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Hadronic Process: Nuclear De-excitations
