Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4BoundMuDecInElIntegration.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4BoundMuDecInElIntegration.cc	(revision 968)
+++ 	(revision )
@@ -1,122 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS(Q2) fixed step integration
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-2005
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-//#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-// Calculation of electron spectra integrals for randomization of bound mu->e+nu+nu decay
-
-double espec(double Z , double T) // T=E_kin in MeV
-{
-  static double mZ=0.,mp,mc,mep,mem,me,mpp;
-  if(Z!=mZ)
-  {
-    double Z2=Z*Z;
-    double Z4=Z2*Z2;
-    double az=std::log(Z);
-    mZ=Z;
-    mp=(11.65+330./Z)/(1.+7.6e-17*Z4*Z4);
-    mc=52.83-.2257*az*az;
-    mep=.674;
-    mem=.1431+Z2/15080.;
-    me=5.54*std::exp(-.0327*std::pow(Z,1.241));
-    mpp=2.253*std::exp(-Z2/8922.7);
-  		G4cout<<"z="<<Z<<":p="<<mp<<",c="<<mc<<",ep="<<mep<<",em="<<mem<<",e="<<me<<",pp="<<mpp
-          <<G4endl;
-		}
-  double p=mp;
-  double c=mc;
-  double ep=mep;
-  double em=mem;
-  double e=me;
-  double pp=mpp;
-  double res=std::pow(T,pp)/(1.+e*std::exp(em*std::pow(T,ep)))/(1.+std::pow(T/c,p));
-  //G4cout<<"T="<<T<<", R="<<.0001389*res<<G4endl;
-  return res;
-}
-int main()
-{
-  const int niX=201;             // number of points for randomization table
-  const int liX=niX-1;           // the last index for randomization table (?)
-  //const int niE=1000000;         // number of points for integration table
-  const int niE=1000;           // number of points for integration table
-  // **********************************************************************
-  double Z=20.;                  // Z of the Element **********************
-  // **********************************************************************
-  double X[niX];                 // randomization table
-  double in[niE];                // integration table
-  double Emax=80.;               // Max kin E in MeV
-  double dE=Emax/niE;            // Step in the integration table
-  double spP=0.;                 // Previous spectrum value [sp(0)=0]
-  in[0]=0.;                      // Not yet normed collected integrals
-  double e=0.;                   // running value of energy
-  double inv=0.;                 // Integral cumulative value
-  for(int i=1; i<niE; i++)
-  {
-    e+=dE;
-    double sp=espec(Z,e);
-    inv+=spP+sp;                 // integrate with the Gauss method
-    in[i]=inv;
-  		//G4cout<<"sp="<<sp<<",in[="<<i<<"]="<<inv<<G4endl;
-    spP=sp;
-  } // End of the big loop over log(E)
-  for(int j=1; j<niE; j++) in[j]/=inv; // Normalize the integrals by a unit
-  //for(int n=1; n<niE; n++) G4cout<<"1-in["<<n<<"]="<<1.-in[n]<<G4endl;
-  double dx=1./liX;              // The last in the table is MAXE by definition
-  double d=dx;                   // Start searching with 1, because x(0)=0
-  int m=1;                       // index in the randomization table
-  double oi=0.;                  // previouse integral value
-  X[0]=0.;                       // Interpolated randomization table (energy)
-  e=0.;
-  for(int k=1; k<niE; k++)
-		{
-    e+=dE;
-    double ci=in[k];
-    if(ci>=d)
-				{
-      if(ci==d) X[m]=e;
-      else X[m]=e-(ci-d)*dE/(ci-oi);
-      m++;
-      d+=dx;                     // start searching for another point
-    }
-				oi=ci;
-  }
-  for(int k=0; k<liX; k++)
-		{
-    G4cout<<X[k]<<",";
-    if(!((k+1)%10))G4cout<<G4endl;
-  }
-  return 0;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4EffectivePhotonsTest.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4EffectivePhotonsTest.cc	(revision 968)
+++ 	(revision )
@@ -1,92 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-#include "G4ios.hh"
-#include "G4Electron.hh"
-#include "G4Element.hh"
-#include "G4ElectroNuclearCrossSection.hh"
-
-int main()
-{
-  G4ElectroNuclearCrossSection eACrossSection;
-  const G4int nN=6;                          // A#of plots (nuclei in the test)
-  ofstream file0("h1_epl.out", ios::out);
-  file0.setf( ios::scientific, ios::floatfield );
-  ofstream file1("h1_eph.out", ios::out);
-  file1.setf( ios::scientific, ios::floatfield );
-  ofstream file2("cu_epl.out", ios::out);
-  file2.setf( ios::scientific, ios::floatfield );
-  ofstream file3("cu_eph.out", ios::out);
-  file3.setf( ios::scientific, ios::floatfield );
-  G4Element* theElement[nN]={new G4Element("Hydrogen", "H", 1, 1.*g/mole),
-							 new G4Element("Deuterium", "D", 1, 2.*g/mole),
-							 new G4Element("Carbon", "C", 6, 12.*g/mole),
-							 new G4Element("Aluminum", "Al", 13, 27.*g/mole),
-							 new G4Element("Copper", "Cu", 29, 64.*g/mole),
-							 new G4Element("Lead", "Pb", 82, 208.*g/mole)};
-  G4ParticleDefinition* theParticleDefinition = G4Electron::ElectronDefinition();
-  G4DynamicParticle* theDynamicParticle;
-  for(G4int n=0; n<4; n++)
-  {
-    G4double ekin=1000.;
-    if(n==3) ekin=10000.;
-    else if(n==1) ekin=10000000.;
-    G4int   nel=0;
-    if(n>1) nel=3;
-    theDynamicParticle = new G4DynamicParticle(theParticleDefinition,
-											   G4ParticleMomentum(1.,0.,0.), ekin*MeV);
-	G4double sig = eACrossSection.GetCrossSection(theDynamicParticle,theElement[nel])/millibarn;
-	G4cout<<"n="<<n<<", Cross Section="<<sig<<G4endl;
-	delete theDynamicParticle;
-	G4int prob[100];
-	for(G4int k=0; k<100; k++) prob[k]=0;
-	for(G4int i=0; i<10000000; i++)
-	{
-	  G4double phe=eACrossSection.GetEffectivePhotonEnergy();
-      //if(phe<100.)G4Exception("Too low photon energy");
-	  G4int ind=static_cast<int>(phe/10.);
-	  if(n==3) ind/=10;
-	  else if(n==1) ind/=10000;
-	  prob[ind]++;
-	  if(!(i%1000000)) G4cout<<i<<" equivalent photons are generated"<<G4endl;
-	}
-	//G4cout<<"Now printing"<<G4endl;
-	for(G4int j=0; j<100; j++)
-	{
-	  //G4cout<<"j="<<j<<G4endl;
-	  if     (!n)   file0<<j*10    <<" "<<prob[j]<<G4endl;
-	  else if(n==1) file1<<j*100000<<" "<<prob[j]<<G4endl;
-	  else if(n==2) file2<<j*10    <<" "<<prob[j]<<G4endl;
-	  else if(n==3) file3<<j*100   <<" "<<prob[j]<<G4endl;
-	}
-  }
-  return EXIT_SUCCESS;
-}
-
-
-
-
-
-
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4ElectroNuclearCrossSectionTest.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4ElectroNuclearCrossSectionTest.cc	(revision 968)
+++ 	(revision )
@@ -1,108 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-#define debug
-#define qedebug
-
-//#include "G4UIterminal.hh"
-#include "G4ios.hh"
-//#include "G4ProcessManager.hh"
-//#include "G4DynamicParticle.hh"
-#include "G4Electron.hh"
-#include "G4Element.hh"
-#include "G4ElectroNuclearCrossSection.hh"
-
-int main()
-{
-  G4ElectroNuclearCrossSection eACrossSection;
-  const G4int nN=6;                          // A#of plots (nuclei in the test)
-  const G4int nC=227;                        // A#of points in the plots (energy values)
-#ifdef qedebug
-  const G4int nP=10;                         // A#of points in nu/Q2 for each eE (energy values)
-#endif
-  std::ofstream fileH1("h1_e.out", std::ios::out);
-  fileH1.setf( std::ios::scientific, std::ios::floatfield );
-  std::ofstream fileH2("h2_e.out", std::ios::out);
-  fileH2.setf( std::ios::scientific, std::ios::floatfield );
-  std::ofstream fileC12("c12_e.out", std::ios::out);
-  fileC12.setf( std::ios::scientific, std::ios::floatfield );
-  std::ofstream fileAl27("al27_e.out", std::ios::out);
-  fileAl27.setf( std::ios::scientific, std::ios::floatfield );
-  std::ofstream fileCu("cu_e.out", std::ios::out);
-  fileCu.setf( std::ios::scientific, std::ios::floatfield );
-  std::ofstream filePb("pb_e.out", std::ios::out);
-  filePb.setf( std::ios::scientific, std::ios::floatfield );
-  G4double low[nN]={137. , 2.22, 12.5 , 8.3 , 7.  , 7.35};
-  G4double high[nN]={2.e7, 2.e7, 2.e7, 2.e7, 2.e7, 2.e7};
-  G4int zA[nN]  ={1, 1,  6, 13, 29,  82};
-  G4int aA[nN]  ={1, 2, 12, 27, 64, 208};
-  G4Element* theElement[nN]={new G4Element("Hydrogen",  "H" , zA[0], aA[0]*g/mole),
-							 new G4Element("Deuterium", "D" , zA[1], aA[1]*g/mole),
-							 new G4Element("Carbon",    "C" , zA[2], aA[2]*g/mole),
-							 new G4Element("Aluminum",  "Al", zA[3], aA[3]*g/mole),
-							 new G4Element("Copper",    "Cu", zA[4], aA[4]*g/mole),
-							 new G4Element("Lead",      "Pb", zA[5], aA[5]*g/mole)};
-  G4ParticleDefinition* theParticleDefinition = G4Electron::ElectronDefinition();
-  G4DynamicParticle* theDynamicParticle;
-  G4double sig=0.;
-  for(G4int n=0; n<nN; n++)                 // >>>> LOOP over Elements
-  {
-	G4double lekin = std::log(low[n]);
-	G4double dlekin= std::exp((std::log(std::log(high[n]))-std::log(lekin))/(nC-1));
-	lekin /= dlekin;
-#ifdef debug
-	G4cout<<"G4ElNucCSTest:>>>>n="<<n<<",l="<<low[n]<<",h="<<high[n]<<",n="<<nC<<",d="<<dlekin<<G4endl;
-#endif
-	for(G4int ll=0; ll<nC; ll++)            // >>>> LOOP over electron energies with the log-step
-	{
-	  lekin*=dlekin;
-      G4double ekin=std::exp(lekin);
-	  theDynamicParticle = new G4DynamicParticle(theParticleDefinition,
-												 G4ParticleMomentum(1.,0.,0.), ekin*MeV);
-	  sig = eACrossSection.GetCrossSection(theDynamicParticle,theElement[n])/millibarn;
-	  delete theDynamicParticle;
-#ifdef debug
-      G4cout<<"G4ElNucCSTest:****>>> Z="<<zA[n]<<", A="<<aA[n]<<", eE="<<ekin<<", sigma="<<sig<<G4endl;
-#endif
-#ifdef qedebug
-	  if(sig>0.)for(G4int qe=0; qe<nP; qe++) // >>>> LOOP over equiv. photon energies for particular eE
-	  {
-	    G4double phE = eACrossSection.GetEquivalentPhotonEnergy();
-	    G4double phQ2= eACrossSection.GetEquivalentPhotonQ2(phE);
-	    G4double x = phQ2/1878./phE;
-        G4cout<<"G4ElNucCSTest:.......eE="<<ekin<<", phE="<<phE<<", phQ2="<<phQ2<<", x="<<x<<G4endl;
-	  } // End of the equivalent photon LOOP
-#endif
-      // Fill the result for the electro-nuclear cross section plot
-      if(!n) fileH1<<ekin<<" "<<sig<<G4endl;
-      else if(n==1) fileH2<<ekin<<" "<<sig<<G4endl;
-      else if(n==2) fileC12<<ekin<<" "<<sig<<G4endl;
-      else if(n==3) fileAl27<<ekin<<" "<<sig<<G4endl;
-      else if(n==4) fileCu<<ekin<<" "<<sig<<G4endl;
-      else if(n==5) filePb<<ekin<<" "<<sig<<G4endl;
-    } // End of the LOOP over electron energies
-  } // End of the LOOP over elements
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4GlauberCrossSections.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4GlauberCrossSections.cc	(revision 968)
+++ 	(revision )
@@ -1,1121 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: Glauber Hadron-Nuclear Cross Sections
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Jan-2005
-// 
-//=====================================================================
-
-//#define debug
-//#define bestest
-//#define integrcs
-#define eldiffcs
-
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-#include "G4QBesIKJY.hh"
-#include "G4QPDGCode.hh"
-#include "G4NucleiPropertiesTable.hh"
-#include "G4NucleiProperties.hh"
-//#include <CLHEP/GenericFunctions/Bessel.hh>
-
-// Integrated cross-sections
-G4double TotalCrSec, InelCrSec, ProdCrSec, ElasticCrSec, QuasyElasticCrSec; 
-// Parameters of hadron-nucleon interaction
-G4double HadrTot, HadrSlope, HadrReIm, DDSect2, DDSect3, Tot00;
-// CHIPS parameters of hadron-proton & hadron-neutron interaction (Re/Im is common)
-G4double HadrTElN, HadrTElP, HadrSlpN, HadrSlpP, HadrPexN, HadrPexP;
-
-// CHIPS Parameters (only for protons, nn is necessary for neutrons np+nn)
-void CHIPSParameters(G4int iPDG, G4double HadrMoment) // HadrMoment is in MeV
-{
-  G4double p=HadrMoment/GeV;                   // Momentum in GeV
-  G4double p2=p*p;                             // Squared momentum in GeV^2
-  G4double sp=std::sqrt(p);
-  G4double p2s=p2*sp;
-  G4double p3=p2*p;
-  G4double p4=p2*p2;
-  G4double p5=p4*p;
-  G4double lp=std::log(p);
-  G4double d2=(lp-5.)*(lp-5.);
-  G4double p8=p4*p4;
-  if (iPDG==2212)
-  {
-    HadrTElN=12./(p2s+.05*p+.0001/sp)+.35/p+(6.75+.14*d2+19./p)/(1.+.6/p4); // pn (mb)
-    HadrTElP=2.91/(p2s+.0024)+(6.75+.14*d2+23./p)/(1.+.4/p4);               // pp (mb)
-    HadrSlpN=(7.2+4.32/(p8+.012*p3))/(1.+2.5/p4);                           // pn (GeV^-2) 
-    HadrSlpP=8.*std::pow(p,.055)/(1.+3.64/p4);                              // pp (GeV^-2) 
-    HadrPexN=(6.75+.14*d2+13./p)/(1.+.14/p4)+.6/(p4+.00013);                // pn (mb/sr)
-    HadrPexP=(74.+3.*d2)/(1.+3.4/p5)+(.2/p2+17.*p)/(p4+.001*sp);            // pn (mb/sr)
-    HadrReIm=-.54+.0954*lp+.975/(p2+.09865/p);           // no unit (the same for pp & pn)
-    //G4cout<<"CHIPSPar:p="<<p<<",N="<<HadrTotN<<",P="<<HadrTotN<<",RI="<<HadrReIm<<G4endl;
-  }     
-  else G4cout<<"CHIPSPar: Only proton projectile (2212) is implemented PDG="<<iPDG<<G4endl;
-}
-
-// Dubna Parameters
-G4int CalculateParameters(G4int iPDG, G4double HadrMoment) // HadrMoment is in MeV
-{
-  static const G4double mN=.938;   // Dubna's mass of the nucleon in nuclei in GeV
-  //static const G4double mN=.9315;  // AtomicMassUnit (mass of the nucleon in nuclei) in GeV
-  static const G4double dmN=mN+mN; // Doubled mass of the nucleon in nuclei in GeV
-  static const G4double mN2=mN*mN; // Squared mass of the nucleon in nuclei in GeV^2
-  G4int iHadron=-1;
-  if(iPDG==2212||iPDG==2112||iPDG==3122||iPDG==3222||iPDG==3112||
-     iPDG==3212||iPDG==3312||iPDG==3322||iPDG==3334) iHadron = 0;
-  else if(iPDG==-2212||iPDG==-2112||iPDG==-3122||iPDG==-3222||
-          iPDG==-3112||iPDG==-3212||iPDG==-3312||iPDG==-3322||
-          iPDG==-3334) iHadron = 1;       // Anti-nucleons, Anti-hyperons
-  else if(iPDG== 211)  iHadron = 2;       // Pi+
-  else if(iPDG==-211)  iHadron = 3;       // Pi-
-  else if(iPDG== 321)  iHadron = 4;       // K+
-  else if(iPDG==-321)  iHadron = 5;       // K- @@ What about K0?
-  else
-  {
-    G4cout<<"CalculateParameters: PDG= "<<iPDG<<" is not supported"<<G4endl;
-    return iHadron;
-  }
-  G4double mHadr      = G4QPDGCode(iPDG).GetMass()/1000.; // Hadron Mass in GeV
-  G4double mHadr2     = mHadr*mHadr;               // Squared Hadron Mass in GeV
-           HadrMoment/= 1000.;                     // Momentum in GeV (input parameter)
-  G4double HadrMoment2= HadrMoment*HadrMoment;     // Squared Momentum in GeV^2
-  G4double HadrEnergy2= mHadr2+HadrMoment2;        // Tot energy in GeV
-  G4double HadrEnergy = std::sqrt(HadrEnergy2);    // Tot energy in GeV (global parameter)
-  G4double sHadr      = HadrEnergy*dmN+mN2+mHadr2; // s in GeV^2
-  G4double sqrS       = std::sqrt(sHadr);          // W in GeV
-  //G4cout<<"GHAD: E="<<HadrEnergy<<",dM="<<dmN<<",sN="<<mN2<<",sH="<<mHadr2<<G4endl;
-  switch (iHadron)
-  {
-  case 0:                                      // =====>>> proton or hyperons (the same?)
-    //G4double Delta=1;                                                   // LowEnergy corr
-    //if(HadrEnergy<40.) Delta = 0.916+0.0021*HadrEnergy;
-    HadrTot = 5.2+5.2*std::log(HadrEnergy)+51*std::pow(HadrEnergy,-0.35); // mb
-    HadrSlope = 6.44+0.88*std::log(sHadr)-1;                              // GeV^-2 
-    HadrReIm  = 0.13*std::log(sHadr/350)*std::pow(sHadr,-0.18);           // no unit
-    //G4cout<<"GHAD:S="<<sHadr<<",T="<<HadrTot<<",R="<<HadrReIm<<",B="<<HadrSlope<<G4endl;
-    DDSect2   = 11.;                                                      // mb*GeV^-2
-    DDSect3   = 3.;                                                       // mb*GeV^-2
-    // Hyperon corrections
-    if(iPDG==3122||iPDG==3222||iPDG==3112||iPDG==3212) //Lambda,Sig0+-
-    {
-      HadrTot   *= 0.80;
-      HadrSlope *= 0.85;
-    }
-    else if(iPDG==3312||iPDG==3322)            // ---> Xi0-
-    {
-      HadrTot   *= 0.70;
-      HadrSlope *= 0.75;
-    }           
-    else if(iPDG==3334)                        // ---> Omega-
-    {
-      HadrTot   *= 0.60;
-      HadrSlope *= 0.65;
-    }
-    break;
-  case 1:                                      // =====>>> antiproton and antihyperons
-    HadrTot = 5.2+5.2*std::log(HadrEnergy)+123.2*std::pow(HadrEnergy,-0.5); //  mb
-    HadrSlope = 8.32+0.57*std::log(sHadr);     // GeV^-2 
-    if(HadrEnergy<1000.) HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8);
-    else HadrReIm  = 0.6*std::log(sHadr/350)*std::pow(sHadr,-0.25);
-    DDSect2 = 11.;                             //mb*GeV-2 @@ the same as in case 0
-    DDSect3 = 3.;                              //mb*GeV-2 @@ the same as in case 0
-    // Hyperon corrections
-    if(iPDG==-3122||iPDG==-3222||iPDG==-3112||iPDG==-3212)//AntiLam,Sig
-    {
-      HadrTot   *=0.75;
-      HadrSlope *=0.85;
-    }
-    else if(iPDG==-3312||iPDG==-3322)          // =====>>> AntiXi0-
-    {
-      HadrTot   *=0.65;
-      HadrSlope *=0.75;
-    }
-    if(iPDG==-3334)                            // ======>>>AntiOmega-
-    {
-      HadrTot   *=0.55;
-      HadrSlope *=0.65;
-    }
-    break;
-  case 2:                                      // =====>>>  pi plus
-    if(HadrMoment>2.) HadrTot = 10.6+2.*std::log(HadrEnergy)+25*std::pow(HadrEnergy,-0.43);// mb
-    else HadrTot = 40-50*(HadrMoment-1.5)*(HadrMoment-1.7);
-    HadrSlope = 7.28+0.245*std::log(sHadr);                                // GeV^-2 
-    HadrReIm  = 0.2*std::log(sHadr/100)*std::pow(sHadr,-0.15);             // no dim
-    DDSect2   = 4.6;                                                       // mb*GeV^-2
-    DDSect3   = 1.33;                                                      // mb*GeV^-2
-    break;
-  case 3:                                      // =====>>>  pi minus
-    HadrTot = 10.6+2*std::log(HadrEnergy)+30*std::pow(HadrEnergy,-0.43);   // mb @@ E->inf?
-    if(HadrMoment<1.399) HadrTot = HadrTot+21.0/0.4*(1.4-HadrMoment);      // A huge jump ?
-    HadrSlope = 7.28+0.245*std::log(sHadr);                                // GeV-2 
-    HadrReIm  = 0.2*std::log(sHadr/100)*std::pow(sHadr,-0.15);
-    DDSect2   = 4.6;                           //mb*GeV-2 @@ the same as in case 2
-    DDSect3   = 1.33;                          //mb*GeV-2
-    break;
-  case 4:                                      // =====>>> K plus
-    HadrTot = 10.6+1.8*std::log(HadrEnergy)+9.*std::pow(HadrEnergy,-0.55); // mb 
-    if(HadrEnergy>100) HadrSlope = 15.;              // Jump at 100 GeV (?)
-    else HadrSlope = 1.0+1.76*std::log(sHadr)-2.84*std::pow(sHadr,-0.5);     // GeV-2
-    //else HadrSlope = 5.28+1.76*std::log(sHadr)-2.84*std::pow(sHadr,-0.5);    // GeV-2
-    HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1);
-    DDSect2   = 3.5;                           //mb*GeV-2
-    DDSect3   = 1.03;                          //mb*GeV-2
-    break;
-  case 5:                                      // ======>>> K minus
-    HadrTot = 10+1.8*std::log(HadrEnergy)+25*std::pow(HadrEnergy,-0.5);      // mb 
-    HadrSlope = 6.98+0.127*std::log(sHadr);           // GeV-2 
-    //if(HadrEnergy<8) HadrReIm = 0.7;
-    HadrReIm  = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1);
-    DDSect2   = 3.5;                           //mb*GeV-2
-    DDSect3   = 1.03;                          //mb*GeV-2
-    break;
-  }     
-  return iHadron;
-}
-
-void  CalculateIntegralCS(G4int  Anuc, G4double HadrEnergy)
-{
-  static const G4double eps = 0.0001;          // Accuracy of calculations
-  static const G4double mb2G2 = 2.568;         // Transformation from mb to GeV^-2
-  static const G4double incrf = 10;            // Increase factor of radius for intergation
-  static const G4double m2G10 = incrf*mb2G2;   // Big conversion factor
-  static const G4double dR0   = 1.06+1.06;     // Doubled nominal R0
-  G4double An      = Anuc;                     // A float
-  G4double Stot    = HadrTot*mb2G2;            // Total hadron-nucleon CS in GeV-2
-  G4double Bhad    = HadrSlope;                // Diffraction cone nH In GeV-2
-  G4double Asq     = 1+HadrReIm*HadrReIm;      // |M|^2/(ImM)^2
-
-  G4double R0      = std::sqrt(0.99);          // in fermi (strange value? 1.07?)
-
-  if (Anuc == 58) R0 = std::sqrt(0.6);         // Personal correction for Fe
-  if (Anuc >20)   R0 = std::sqrt((35.34+0.5*Anuc)/(40.97+Anuc));
-  if (Anuc == 16) R0 = std::sqrt(0.75);        // Personal correction for O16
-  if (Anuc >10)   R0 = std::sqrt(0.84);        // The same R0 for A=11-20 (except O16)
-  if (Anuc ==  4) R0 *= 1.2;                   // Personal correction for Helium
-
-  G4double Rnucl   = R0*std::pow(static_cast<double>(Anuc),.3333333); // in fermi
-  G4double Rnuc2   = Rnucl*Rnucl*m2G10;        // in GeV^-2 (10 -> to make it >> R)
-  G4double RB      = Rnuc2+Bhad;               // Effective increase of R2 (Convolution?)
-  G4double R2B     = RB+Bhad;                  // Over-convolution (?)
-  G4double Delta   = Stot/R2B/twopi;           // Integration step for Total & Inelastic CS
-  G4double Delt2   = Delta+Delta;              // Two delta
-  G4double Delt3   = Stot*Asq*R2B/RB/Bhad/16/pi; // Integration step for Production CS
-
-  G4double Tot0    = 0.;                       // =SUM[(-Delta)^(i-1)*(C^A_i)/i]
-  G4double Inel0   = 0.;                       // =SUM[(-Delta)^(i-1)*(2A)!/(2A-i)!/i!/i]
-  G4double N       = -1./Delta; 
-  G4double N1      = N;
-  G4double N3      = -1./Delt2;
-  G4double Prod0   = 0.;
-  G4int    Ai      = Anuc+1;                   // Working value to avoid multiplication
-  G4double RBi     = 0.;                       // Working value to avoid multiplication
-  for (G4int i=1; i<= Anuc; i++)               // ==> Loop over nucleons
-  {
-    Ai--;                                      // Working value to avoid multiplication
-    RBi   += RB;                               // Working value to avoid multiplication
-    G4double Di=Delta/i;                       // Working value to avoid repetition
-    N     *= -Di*Ai;                           // (-Delta)^(i-1)*A!/(A-i)!/i! (C^A_i)
-    N1    *= -Di*(Anuc+Ai);                    // (-Delta)^(i-1)*(2A)!/(2A-i)!/i!
-    N3    *= -Delt2*Ai/i;                      // (-2*Delt)^(i-1)*A!/(A-i)!/i! (C^A_i)
-    Tot0  += N/i;                              // =SUM[(-Delta)^(i-1)*(C^A_i)/i]
-    Inel0 += N1/i;                             // =SUM[(-Delta)^(i-1)*(2A)!/(2A-i)!/i!/i]
-    //G4double N2    = 1.;
-    G4double N4    = 1.;
-    //G4double Inel1 = 0.;
-    //G4double Inel2 = 0.;
-    G4double Prod1 = 0.;
-    G4double Bhl   = 0.;
-    for (G4int l=0; l<= i; l++)                // ==> Loop over nucleon-nucleon pairs (0?)
-    {
-      //Inel1 += N2/(i+l);                 // SUM[(-Delta)^(l-1)*i!(i-1)!/(l-1)!/l!/(i+l)!]
-      Prod1 += N4*RB/(RBi+Bhl);                //
-      //N2    *= -Delt*(i-l)/(l+1);              // (-Delta)^l*i!/l!/(l+1)!
-      N4    *= -Delt3*(i-l)/(l+1);
-      Bhl   += Bhad;
-    }
-    //Inel2 += Inel1*N3;
-    Prod0   += Prod1*N3;
-    if(std::fabs(N1/i/Inel0) < eps)  break;
-  }
-
-  Tot0             *= HadrTot;                 // Multiply by the hN cross section
-  Inel0            *= HadrTot/2;               // Multiply by a half of hN cross section
-  //Inel2            *= HadrTot;                 // Multiply by the hN cross section
-  Prod0            *= HadrTot;                 // Multiply by the hN cross section
-  Tot00             = Tot0;                    // Remember because it will be changed
-  G4double ak       = Rnuc2*twopi/Stot;        // cross-section ratio for integration
-  G4double Ank      = An/ak;                   // Working value to avoid repetition
-  G4double rak      = 1./ak;                   // Working value to avoid repetition
-  G4double ak1      = 1.-rak/4;                // Working value to avoid repetition
-  G4double Ank1     = Ank*ak1;                 // Working value to avoid repetition
-  G4double DDSect1  = (DDSect2+DDSect3*std::log(dR0*HadrEnergy/Rnucl/std::sqrt(m2G10)/4));
-  G4double Dtot     = 8*pi*ak/HadrTot*(1.-(1.+Ank)*std::exp(-Ank))*DDSect1/mb2G2;//TotCSCor
-  G4double bk       = (1.-rak)/Stot/ak1;       // Working value to avoid power(,2)
-  G4double bd       = bk*bk*DDSect1*(1.-(1.+Ank1)*std::exp(-Ank1))*Rnuc2; // Working
-  G4double Dprod    = bd*4*pi2*mb2G2;          // Correction for the Production & Inel CS
-  TotalCrSec        = Tot0-Dtot;               // Total Cross Section         (primary)
-  InelCrSec         = Inel0-Dprod;             // Inelastic Cross Section     (primary)
-  //InelCrSec1        = Inel2;                 // Direct calculation without correction
-  ProdCrSec         = Prod0-Dprod;             // Production Cross Section    (primary)
-  ElasticCrSec      = TotalCrSec-InelCrSec;    // Elastic Cross Section       (secondary)
-  QuasyElasticCrSec = InelCrSec-ProdCrSec;     // Quasy-Elastic Cross Section (secondary)
-}
-
-#ifdef eldiffcs
-// Functions for calculation of differential elastic cross-sections
-// Static globals
-static const G4double  rAmax = 2.5;      // Factor of maximum radius
-static const G4int     NptB  = 500;      // A#of steps in the impact parameter integration
-static const G4int     nFact = 253;      // Maximum A (for Factorials and Factors)
-
-// Globals
-G4double r0 = 1.1;                       // Wood-Saxon's factor
-G4double Mnoj[nFact], Factorials[nFact]; // Factors and factorials
-G4double R1, R2, Pnucl, Aeff;            // Nuclear Parameters
-G4double InCoh;                          // Incoherent quasi-elastic cross section
-
-G4double binom(G4int N, G4int M)
-{
-  if(N>170 && N>M && M)
-		{
-    G4double fN=N;
-    G4double fM=M;
-    G4double fD=N-M;
-    G4double man=N*std::log(fN)-M*std::log(fM)-fD*std::log(fD)-1.;
-    //G4cout<<"G4GCS::binom: N="<<N<<", M="<<M<<", C="<<man<<", E="<<std::exp(man)<<G4endl;
-    return std::exp(man);
-		}
-  else if(N>1 && N>M && M) return Factorials[N]/Factorials[M]/Factorials[N-M];
-  return 1.;
-}
-
-// Initialize factorials and combinatoric coefficients
-void Initialize()
-{
-  //static const G4double  sat = 3.03;      // Stirling saturation of Dubna
-  G4double Val = Factorials[0] = Mnoj[0] = 1.;
-  for (G4int i=1; i<nFact; i++)
-  {
-    G4double Si = 0.;
-             Val *= i;                    // Calculation of the factorial values
-    Factorials[i] = Val;
-    for (G4int l = 0; l<=i; l++)
-    {
-      G4double Sm = 1.;                   // sum of reversed binom coefficients (?)
-      for (G4int m=1; m<=i-l; m++) Sm += 1./binom(i-l,m);
-      //if(i>169)G4cout<<"G4GCS::Init:S="<<Sm<<", b("<<i<<","<<l<<")="<<binom(i,l)<<G4endl;
-      Si += Sm/binom(i,l);
-    } //  End of l LOOP
-    Mnoj[i] = Si;
-    //G4double d=Si-sat;
-    //G4double r=d/sat;
-    //if(i>100)G4cout<<"G4GCS::Init:i="<<i<<":"<<Si<<"-"<<sat<<" = "<<d<<", r="<<r<<G4endl;
-  } // End of i LOOP
-} // End of Initialization
-
-// Dubna paranmeters of nuclei
-void CalculateNuclearParameters(G4int Anuc)  // R1,R2,Pnuc,Aeff
-{//  ==============================================
-  G4double A=Anuc;
-  // Personal corrections fore some nuclei
-  if(Anuc == 4)                                      // Personal correction for He (!)
-  {
-    R1    = 5.5;
-    R2    = 3.7;
-    Pnucl = 0.4;   
-    Aeff  = 0.87;
-  }
-  else if(Anuc == 9)                                 // Personal correction for Be (!)
-  {
-    R1    = 9.0;
-    R2    = 7.0;
-    Pnucl = 0.190;
-    Aeff  = 0.9;
-  }
-  //else if(Anuc == 11)                              // Personal correction for B (!)
-  //{
-  //  R1    = 10.8;
-  //  R2    = 7.5;
-  //  Pnucl = 0.85;
-  //  Aeff  = 1.2;
-  //}
-		//else if(Anuc == 12)                              // Personal correction for C (!)
-  //{
-  //   R1    = 9.336;
-  //   R2    = 5.63;
-  //   Pnucl = 0.197;
-  //   Aeff  = 1.0;
-  //}
-  else if(Anuc == 16)                                // Personal correction for O (!)
-  {
-    R1    = 10.50;
-    R2    = 5.5;
-    Pnucl = 0.7;
-    Aeff  = 0.98;
-    //R1    = 11.3;
-    //R2    = 2.5;
-    //Pnucl = 0.75;
-    //Aeff  = 0.9;
-  }
-  else if(Anuc == 58)                                // Personal correction for Ni (!)
-  {
-    R1    = 15.0;
-    R2    = 9.9;
-    Pnucl = 0.45;
-    Aeff  = 0.85;
-  }
-  else if(Anuc == 90)                                // Personal correction for Zr (!)
-  {
-    //R1    = 16.5;
-    //R2    = 11.62;
-    //Pnucl = 0.4;
-    //Aeff  = 0.9;
-    R1    = 16.5;
-    R2    = 11.62;
-    Pnucl = 0.4;
-    Aeff  = 0.7;
-  }
-  else if(Anuc == 208)                               // Personal correction for Pb (!)
-  {  
-    //      R1 = 20.73; R2 = 15.74.
-    //R1    = 4.1408*std::pow(A,0.3018);
-    //R2    = 3.806*std::pow(A-10.068,0.2685);
-    //Pnucl = 0.9;
-    //Aeff  = 1.1;
-    R1    = 19.5;
-    R2    = 15.74;
-    Pnucl = 0.4;
-    Aeff  = 0.7;
-  }
-  else
-  {
-    R1 = 4.45*std::pow(A-1.,.309);                  // First diffractional radius
-    if(Anuc == 28) R1 *= .955;                       // Personal correction for Si (?!)
-    R2    = 2.3*std::pow(A,.36);                    // Second diffraction radius
-    Pnucl = 0.176+(.00167+.00000869*Anuc)*Anuc;      // Momentum of mean Fermi motion
-    Aeff  = 0.9;                                     // Effective screaning
-  }
-}
-
-// Independent transition method calculation. Has problems for A>95 (short cut).
-G4double DiffElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) // All in MeV
-{//      =================================================================
-  static const G4double eps = 0.000001;                     // Accuracy of calculations
-  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
-  static const G4double piMG  = pi/mb2G2;                   // PiTrans from GeV^-2 to mb
-  static const G4double dR0   = 1.06+1.06;                  // Doubled nominal R0
-  G4double hMom       = HadrMom/1000.;                      // Momentum (GeV, inputParam)
-		G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
-  G4double MassH2     = MassH*MassH;
-  G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
-		//G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl;
-  if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3");
-  if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!");
-  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
-  CalculateParameters(hPDG, HadrMom);                   
-  CalculateIntegralCS(A, HadrEnergy);
-  CalculateNuclearParameters(A);
-  G4double Q2 = aQ2/1000/1000;                              // in GeV^2
-		//MassN           = A*0.938;                              // @@ This is bad!
-		G4double MassN    = A*0.9315;                             // @@ Atomic Unit is bad too
-  G4double MassN2   = MassN*MassN;
-  G4double S        = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s
-  G4double EcmH     = (S-MassN2+MassH2)/2/std::sqrt(S);     // CM energy of a hadron
-  G4double CMMom    = std::sqrt(EcmH*EcmH-MassH2);          // CM momentum
-
-  G4double Stot     = HadrTot*mb2G2;                        // in GeV^-2
-  G4double Bhad     = HadrSlope;                            // in GeV^-2
-  G4double Asq      = 1+HadrReIm*HadrReIm;                  // |M|^2/(ImM)^2
-  G4double Rho2     = std::sqrt(Asq);                       // M/ImM
-  G4double R12      = R1*R1;
-  G4double R22      = R2*R2;
-  G4double R12B     = R12+Bhad+Bhad;
-  G4double R22B     = R22+Bhad+Bhad;
-  G4double R12Ap    = R12+20;
-  G4double R22Ap    = R22+20;
-  G4double R13Ap    = R12*R1/R12Ap;
-  G4double R23Ap    = R22*R2/R22Ap*Pnucl;
-  G4double R23dR13  = R23Ap/R13Ap;
-  G4double R12Apd   = 2./R12Ap;
-  G4double R22Apd   = 2./R22Ap;
-  G4double R122Apd  = (R12Apd+R22Apd)/2;
-  G4double dR0H     = dR0*HadrEnergy/4;
-  G4double DDSec1p  = DDSect2+DDSect3*std::log(dR0H/R1);
-  G4double DDSec2p  = DDSect2+DDSect3*std::log(dR0H/std::sqrt((R12+R22)/2));
-  G4double DDSec3p  = DDSect2+DDSect3*std::log(dR0H/R2);
-  G4double Norm     = (R12*R1-Pnucl*R22*R2)*Aeff;    // Some questionable norming
-  G4double R13      = R12*R1/R12B;
-  G4double R23      = Pnucl*R22*R2/R22B;
-  G4double norFac   = Stot/twopi/Norm;               // totCS (in GeV^-2) is used
-  G4double Unucl    = norFac*R13;
-  G4double UnuclScr = norFac*R13Ap;
-  G4double SinFi    = HadrReIm/Rho2;
-  G4double FiH      = std::asin(SinFi);
-  G4double N        = -1.;
-  G4double N2       = R23/R13;
-  G4double ImElA0   = 0.;
-  G4double ReElA0   = 0.;
-  G4double Tot1     = 0.;
-  for(G4int i=1; i<=A; i++)                          // @@ Make separately for n and p
-  {
-    N              *= -Unucl*(A-i+1)*Rho2/i;         // Includes total cross-section
-    G4double N4     = 1.;
-    G4double Prod1  = std::exp(-Q2*R12B/i/4)*R12B/i; // Includes the slope
-    G4double medTot = R12B/i;
-    for(G4int l=1; l<=i; l++)
-    {
-      G4double exp1 = l/R22B+(i-l)/R12B;
-      N4           *= -(i-l+1)*N2/l;
-      G4double Nexp = N4/exp1;
-      Prod1        += Nexp*std::exp(-Q2/exp1/4);
-      medTot       += Nexp;
-    }
-    G4double FiHi   = FiH*i;
-    G4double nCos   = N*std::cos(FiHi);
-    ReElA0 += Prod1*N*std::sin(FiHi);
-    ImElA0 += Prod1*nCos;
-    Tot1           += medTot*nCos;
-    if(std::abs(Prod1*N/ImElA0) < eps) break;
-  }
-  ImElA0   *= piMG;                                  // The amplitude in mb
-  ReElA0   *= piMG;                                  // The amplitude in mb
-  Tot1             *= piMG+piMG;
-  G4double N1p      = 1.;
-  G4double R13Ap1   = DDSec1p*R13Ap*R13Ap/2;
-  G4double R22Ap1   = DDSec2p*R13Ap*R23Ap;
-  G4double R23Ap1   = DDSec3p*R23Ap*R23Ap*R22Ap/2;
-  G4double R13Ap2   = R13Ap1*R12Ap/4;
-  G4double R22Ap2   = R22Ap1/R122Apd/2;
-  G4double R23Ap2   = R23Ap1*R22Ap/4;
-  G4double Q28      = Q2/8;
-  G4double Q24      = Q28+Q28;
-  G4double Din1     = R13Ap2*std::exp(-Q28*R12Ap)-R22Ap2*std::exp(-(Q28+Q28)/R122Apd)+
-				                  R23Ap2*std::exp(-Q28*R22Ap);   // i=0 start value
-  G4double DTot1    = R13Ap2-R22Ap2+R23Ap2;          // i=0 start value
-  if (A<96)
-		{
-    for(G4int i=1; i<A-1; i++)
-    {
-      N1p              *= -UnuclScr*(A-i-1)/i*Rho2;
-      G4double N2p      = 1.;
-      G4double Din2     = 0.;
-      G4double DmedTot  = 0.;
-      G4double BinCoeff = 1.;                        // Start for Binomial Coefficient
-      for(G4int l = 0; l<=i; l++) 
-      {
-        if(l) BinCoeff *= (i-l+1)/l;
-        G4double exp1   = l/R22B+(i-l)/R12B;
-        G4double exp1p  = exp1+R12Apd;
-        G4double exp2p  = exp1+R122Apd;
-        G4double exp3p  = exp1+R22Apd;
-        G4double NBinC  = N2p*BinCoeff;
-        Din2           += NBinC*(R13Ap1/exp1p*std::exp(-Q24/exp1p)-
-                                 R22Ap1/exp2p*std::exp(-Q24/exp2p)+
-                                 R23Ap1/exp3p*std::exp(-Q24/exp3p));
-        DmedTot        += NBinC * (R13Ap1/exp1p - R22Ap1/exp2p + R23Ap1/exp3p);
-        N2p            *= -R23dR13;
-      }
-      G4double comFac   = N1p*Mnoj[i]/(i+2)/(i+1)*std::cos(FiH*i);
-				  Din1             += Din2*comFac;
-      DTot1            += DmedTot*comFac;
-      if(std::abs(Din2*N1p/Din1) < eps) break;
-    }
-    G4double cFac       = A*(A-1)*4/Norm/Norm;
-    Din1               *= -cFac;
-    DTot1              *= cFac;
-  }
-  //else Din1=0.;
-
-  G4double Corr0      = Tot00/Tot1;
-  ImElA0             *= Corr0;
-                        
-  //return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*CMMom*CMMom*mb2G2/4/pi2; // ds/do
-  return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*mb2G2/(twopi+twopi);
-} // End of DiffElasticCS
-
-G4double CHIPSDiffElCS(G4int hPDG, G4double HadrMom, G4int Z, G4int A, G4double aQ2) // MeV
-{//      =================================================================
-  static const G4double eps = 0.000001;                     // Accuracy of calculations
-  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
-  static const G4double piMG  = pi/mb2G2;                   // PiTrans from GeV^-2 to mb
-  G4double p          = HadrMom/1000.;                      // Momentum (GeV, inputParam)
-  G4double p2         = p*p;
-		G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
-  G4double MassH2     = MassH*MassH;
-  G4double HadrEnergy = std::sqrt(p2+MassH2);               // Tot energy in GeV
-		//G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl;
-  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
-  CalculateParameters(hPDG, HadrMom);                       // ?                   
-  CalculateIntegralCS(A, HadrEnergy);                       // ?
-  CalculateNuclearParameters(A);                            // ?
-  G4double Q2 = aQ2/1000/1000;                              // in GeV^2
-		//G4doubleMassN     = A*0.938;                              // @@ This is bad!
-		G4double MassN    = A*0.9315;                             // @@ Atomic Unit is bad too
-  if(G4NucleiPropertiesTable::IsInTable(Z,A))
-				MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.;    // Geant4 NuclearMass in GeV
-  G4double MassN2   = MassN*MassN;
-  G4double S        = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s
-  G4double EcmH     = (S-MassN2+MassH2)/2/std::sqrt(S);     // CM energy of a hadron
-  G4double CMMom    = std::sqrt(EcmH*EcmH-MassH2);          // CM momentum
-  //G4cout<<"P="<<CMMom<<",E="<<HadrEnergy<<",N="<<MassN2<<",h="<<MassH2<<",p="<<p<<G4endl;
-  G4double p3       = p2*p;
-  G4double p4       = p2*p2;
-  G4double sp       = std::sqrt(p);
-  G4double p2s      = p2*sp;
-  G4double ap       = std::log(p);
-  G4double dl       = ap-3.;
-  G4double dl2      = dl*dl;
-  G4double shPPTot  = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb
-  G4double shPNTot  = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p
-                      +(38.+.3*dl2+8./p)/(1.+1.2/p3);       // SigPP in mb
-  G4double hNTot    = (Z*shPPTot+(A-Z)*shPNTot)/A;          // SighN in mb
-  G4double Stot     = hNTot*mb2G2;                          // in GeV^-2
-  G4double shPPSl   = 8.*std::pow(p,.055)/(1.+3.64/p4);     // PPslope in GeV^-2
-  G4double shPNSl   = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2
-  G4double Bhad     = (Z*shPPSl+(A-Z)*shPNSl)/A;            // B-slope in GeV^-2
-  G4double hNReIm   = -.55+ap*(.12+ap*.0045);               // Re/Im_hN in no unit
-  G4double Asq      = 1+hNReIm*hNReIm;                      // |M|^2/(ImM)^2
-  G4double Rho2     = std::sqrt(Asq);                       // M/ImM
-  G4double r1       = 3.9*std::pow(A-1.,.309);             // Positive diffractionalRadius
-  G4double r2       = 2.*std::pow(A,.36);                  // Negative diffraction radius
-  G4double pN       = Pnucl;                                // Dubna value
-  //G4double pN       = .4;                                   // Screaning factor
-  G4double Ae       = Aeff;                                 // Dubna value
-  //G4double Ae       = .75;                                   // Normalization
-  G4double R12      = r1*r1;
-  G4double R22      = r2*r2;
-  G4double R1C      = R12*r1;
-  G4double R2C      = R22*r2;
-  G4double R12B     = R12+Bhad+Bhad;                        // Slope is used
-  G4double R22B     = R22+Bhad+Bhad;                        // Slope is used
-  G4double Norm     = (R1C-pN*R2C)*Ae;                      // ScreanFac & NormFac are used
-  G4double R13      = R1C/R12B;                             // Slope is used
-  G4double R23      = pN*R2C/R22B;                          // Slope & ScreanFact are used
-  G4double norFac   = Stot/twopi/Norm;                      // totCS (in GeV^-2) is used
-  G4double Unucl    = norFac*R13;
-  G4double SinFi    = hNReIm/Rho2;                          // Real part
-  G4double FiH      = std::asin(SinFi);
-  G4double N        = -1.;
-  G4double N2       = R23/R13;                              // Slope is used
-  G4double ImElA0   = 0.;
-  G4double ReElA0   = 0.;
-  G4double Tot1     = 0.;
-  for(G4int i=1; i<=A; i++)
-  {
-    N              *= -Unucl*(A-i+1)*Rho2/i;                // TotCS is used
-    G4double N4     = 1.;
-    G4double Prod1  = std::exp(-Q2*R12B/i/4)*R12B/i;        // Slope is used
-    G4double medTot = R12B/i;                               // Slope is used
-    for(G4int l=1; l<=i; l++)
-    {
-      G4double exp1 = l/R22B+(i-l)/R12B;                    // Slope is used
-      N4           *= -(i-l+1)*N2/l;                        // Slope is used
-      G4double Nexp = N4/exp1;
-      Prod1        += Nexp*std::exp(-Q2/exp1/4);
-      medTot       += Nexp;
-    }
-    G4double FiHi   = FiH*i;
-    G4double nCos   = N*std::cos(FiHi);
-    ReElA0 += Prod1*N*std::sin(FiHi);
-    ImElA0 += Prod1*nCos;
-    Tot1           += medTot*nCos;
-    if(std::abs(Prod1*N/ImElA0) < eps) break;
-  }
-  ImElA0         *= piMG;                          // The amplitude in mb
-  ReElA0         *= piMG;                          // The amplitude in mb
-  Tot1           *= piMG+piMG;
-  G4double Corr0  = Tot00/Tot1;
-  ImElA0         *= Corr0;
-                        
-  //return (ReElA0*ReElA0+ImElA0*ImElA0)*CMMom*CMMom*mb2G2/4/pi2; // ds/do
-  return (ReElA0*ReElA0+ImElA0*ImElA0)*mb2G2/(twopi+twopi);
-} // End of DiffElasticCS
-
-G4double CHIPS_Tb(G4int  A, G4double b)              // T(b) in fm-2
-{
-  static G4int Am=0;                                 // Associative memory value A
-  static G4double B=0.;                              // Effective edge
-  static G4double C=0.;                              // Normalization factor
-  static G4double D=0.;                              // Slope parameter
-  if(A!=Am)
-		{
-    B=.0008*A*A;                                     // no units
-    D=.42*std::pow(A,-.26);                               // fm^-2
-    C=A*D/pi/std::log(1.+B);                         // fm^-2
-  }
-  G4double E=B*std::exp(-D*b*b);
-  return C*E/(1.+E);                                 // T(b) in fm^-2
-}
-
-G4double Thickness(G4int A, G4double b, G4double R)  // T(b) in fm^-2
-{//      ==========================================  // rAmax=2.5 is a GlobStatConstPar
-  static const G4double bTh  = .545;                 // Surface thicknes
-  static const G4double bTh2 = bTh*bTh;              // SquaredSurface thicknes
-  G4double b2    = b*b;                              // Squared impact parameter
-  G4double dr    = rAmax*R/(NptB-1);                 // Step of integration in z
-  G4double R2    = R*R;                              // Working parameter
-  G4double Norm  = .75/pi/R2/R/(1.+bTh2*pi2/R2);     // Corrected Volume (?)
-  G4double r     = -dr;                              // Pre value for the radius
-  G4double SumZ  = 0.;                               // Prototype of the result
-  //G4double SumN  = 0.;                             // is not used
-  for(G4int k=0; k<NptB; k++)
-  {
-    r    += dr;
-    SumZ += 1./(1.+std::exp((std::sqrt(b2+r*r)-R)/bTh)); // Woods-Saxon density
-    //SumN += r*r/(1.+std::exp((r-rAfm)/bTh));       // is not used
-  }
-  //SumN *= (twopi+twopi)*dr*Norm;                   // is not used
-  return SumZ*(A+A)*dr*Norm;                         // T(b) in fm^-2
-} // End of Thickness
-
-G4double CoherentDifElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) //in MeV
-{//      =========================================================================
-  //static const G4double eps = 0.000001;                   // CalculationAccuracy@@NotUsed
-  //static const G4double mN = .9315;                      // Atomic Unit GeV
-  static const G4double mN = .938;                      // Atomic Unit GeV
-  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
-  static const G4double incrf = 10;            // Increase factor of radius for intergation
-  static const G4double m2G10 = incrf*mb2G2;                // Big conversion factor
-  static const G4double s2G10 = std::sqrt(m2G10);           // sqrt Big conversion factor
-  G4double ReIntegrand[NptB],ImIntegrand[NptB],Thick[NptB]; // Calculated arrays
-  G4QBesIKJY QI0(BessI0);                                   // I0 Bessel function
-  G4QBesIKJY QJ0(BessJ0);                                   // J0 Bessel function
-  G4double hMom       = HadrMom/1000.;                      // Momentum (GeV, inputParam)
-		G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
-  G4double MassH2     = MassH*MassH;
-  G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
-  if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3");
-  if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!");
-  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
-  G4double Re         = std::pow(A, 0.33333333);
-  G4double Re2        = Re*Re;
-  CalculateParameters(hPDG, HadrMom);                   
-  G4double Q2 = aQ2/1000/1000;                              //  GeV
-  // Individual corrections for nuclei which had data
-  if    (A==208) r0   = 1.125;
-  else if(A==90) r0   = 1.12;
-  else if(A==64) r0   = 1.1;
-  else if(A==58) r0   = 1.09;
-  else if(A==48) r0   = 1.07;
-  else if(A==40) r0   = 1.15;
-  else if(A==28) r0   = 0.93;
-  else if(A==16) r0   = 0.92;
-  else if(A==12) r0   = 0.80;
-		else           r0   = 1.16*(1.-1.16/Re2);  // For other nuclei which have not been tested
-		G4double MassN      = A*0.9315;                             // @@ Atomic Unit is bad too
-  G4double MassN2     = MassN*MassN;
-  G4double S          = (MassN+MassN)*HadrEnergy+MassN2+MassH2; // Mondelstam S
-  G4double EcmH       = (S-MassN2+MassH2)/2/std::sqrt(S);   // Hadron CM Energy
-  G4double CMMom      = std::sqrt(EcmH*EcmH-MassH2);        // CM momentum
-  G4double rAfm       = r0*Re;                              // Nuclear radius
-  G4double rAGeV      = rAfm*s2G10;                         // Big integration radius
-  G4double stepB      = rAmax*rAGeV/(NptB-1);               // dr step of integration
-  G4double hTotG2     = HadrTot*mb2G2;
-  G4double ReSum      = 0.;
-  G4double ImSum      = 0.;
-  G4double ValB       = -stepB;
-  for(G4int i=0; i<NptB; i++)
-  {
-    ValB             += stepB;                              // An incident parameter
-    G4double ValB2    = ValB*ValB;                          // A working value
-    G4double IPH      = ValB/HadrSlope;                     // A working value
-    G4double dHS      = HadrSlope+HadrSlope;                // A working value
-    G4double Integ    = 0.;
-    G4double ValS     = 0.;
-    for(G4int j=1; j<NptB; j++)
-    {
-      ValS           += stepB;                              // Impact parameter GeV^-1
-      if(!i) Thick[j] = Thickness(A,ValS/s2G10,rAfm)/m2G10; // Calculate only once
-      G4double FunS   = IPH*ValS;
-      if(FunS > 320.) break;
-      Integ          += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Thick[j];
-    } 
-    G4double InExp    = -hTotG2*Integ*stepB/dHS;
-    G4double expB     = std::exp(InExp);
-    G4double HRIE     = HadrReIm*InExp;
-    ReIntegrand[i]    = (1.-expB*std::cos(HRIE));           // Real part of the amplitude
-    ImIntegrand[i]    = expB*std::sin(HRIE);                // Imaginary part of the amplit
-  } 
-  InCoh     = 0.;                                           // incohirent (quasi-elastic)
-  ValB       = -stepB;
-  for(G4int k=0; k<NptB; k++)                               // Third integration (?)
-  {
-    ValB              += stepB;
-    InCoh             += Thick[k]*ValB*std::exp(-hTotG2*Thick[k]);
-    G4double J0qb      = QJ0(std::sqrt(Q2)*ValB)*ValB;      // Bessel0(Q*b)
-    ReSum             += J0qb*ReIntegrand[k];
-    ImSum             += J0qb*ImIntegrand[k];
-  }
-  //G4cout<<"GHAD:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2="
-  //      <<Q2<<",m="<<mb2G2<<G4endl;
-  InCoh *= stepB*hTotG2*hTotG2*(1.+HadrReIm*HadrReIm)*std::exp(-HadrSlope*Q2)/8/mb2G2;
-  //G4cout<<"GHAD:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl;
-  //return (ReSum*ReSum+ImSum*ImSum)*m2G10*CMMom*CMMom*stepB*stepB/twopi; // ds/do
-  return (ReSum*ReSum+ImSum*ImSum)*m2G10*stepB*stepB/12; // ds/dt
-}  
-
-G4double CHIPSDifElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2)
-{//      ============================================================================
-  static const G4int Npb      = 500;                        // A#of intergation points
-  //static const G4double mN = .9315;                         // Atomic Unit GeV
-  static const G4double mN = .938;                          // Atomic Unit GeV
-  static const G4double hc2   = .3893793;                   // Transform from GeV^-2 to mb
-  static const G4double mb2G2 = 1./hc2;                     // Transform from mb to GeV^-2
-  static const G4double f22mb = 10;                         // Transform from fermi^2 to mb
-  static const G4double f22G2 = f22mb*mb2G2;                // Transform from fm2 to GeV^-2
-  static const G4double f2Gm1 = std::sqrt(f22G2);           // Transform from fm to GeV^-1
-  G4double Re         = std::pow(A,.33333333);              // A-dep coefficient
-  G4double Lim        = 100*Re;                             // Integration accuracy limit
-  G4double Tb[Npb];                                         // Calculated T(b) array
-  G4QBesIKJY QI0(BessI0);                                   // I0 Bessel function
-  G4QBesIKJY QJ0(BessJ0);                                   // J0 Bessel function
-  G4double p          = hMom/1000.;                         // Momentum (GeV, inputParam)
-  G4double p2         = p*p;
-		G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
-  G4double MassH2     = MassH*MassH;                        // Squared mass of the hadron
-  G4double hEnergy    = std::sqrt(p2+MassH2);               // Tot energy in GeV
-  G4double Q2         = aQ2/1000000.;                       // -t in GeV
-		G4double MassN    = A*0.9315;                             // @@ Atomic Unit is bad too
-  if(G4NucleiPropertiesTable::IsInTable(Z,A))
-				MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.;    // Geant4 NuclearMass in GeV
-  G4double MassN2     = MassN*MassN;                        // Squared mass of the target
-  G4double S          = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s
-  G4double EcmH       = (S-MassN2+MassH2)/2/std::sqrt(S);   // Hadron CM Energy
-  G4double CMMom      = std::sqrt(EcmH*EcmH-MassH2);        // CM momentum (to norm CS)
-  //G4cout<<"2: P="<<CMMom<<",E="<<hEnergy<<",N="<<MassN2<<",h="<<MassH2<<",p="<<p<<G4endl;
-  // The mean value of the total can be used
-  G4double p3         = p2*p;
-  G4double p4         = p2*p2;
-  G4double sp         = std::sqrt(p);
-  G4double p2s        = p2*sp;
-  G4double ap         = std::log(p);
-  G4double dl         = ap-3.;
-  G4double dl2        = dl*dl;
-  G4double shPPTot    = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb
-  G4double shPNTot    = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p
-                        +(38.+.3*dl2+8./p)/(1.+1.2/p3);     // SigPP in mb
-  G4double shNTot     = (Z*shPPTot+(A-Z)*shPNTot)/A;        // SighN in mb
-#ifdef debug
-  G4cout<<"CHIPS:SI,p="<<p<<",n="<<shNTot<<",P="<<shPPTot<<",N="<<shPNTot
-        <<",Z="<<Z<<",A="<<A<<G4endl;
-#endif
-  G4double shPPSl     = 8.*std::pow(p,.055)/(1.+3.64/p4);   // PPslope in GeV^-2
-  G4double shPNSl     = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2
-  G4double shNSl      = (Z*shPPSl+(A-Z)*shPNSl)/A;          // B-slope in GeV^-2
-#ifdef debug
-  G4cout<<"CHIPS:SL,n="<<shNSl<<",P="<<shPPSl<<",N="<<shPNSl<<G4endl;
-#endif
-  G4double shNReIm    = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit
-  //G4cout<<"CHPS: s="<<S<<",T="<<shNTot<<",R="<<shNReIm<<",B="<<shNSl<<G4endl;
-  // @@ End of temporary ^^^^^^^
-  G4double dHS        = shNSl+shNSl;                        // Working: doubled B-slope
-  G4double stepB      = (Re+Re+2.7)*f2Gm1/(Npb-1);          // in GeV^-1, step of integral
-  G4double hTotG2     = shNTot*mb2G2;                       // sigma_hN in GeV^-2
-  G4double ReSum      = 0.;                                 // Integration of RePart of Amp
-  G4double ImSum      = 0.;                                 // Integration of ImPart of Amp
-  G4double ValB       = -stepB;
-  for(G4int i=0; i<Npb; i++)                                // First integration over b
-  {
-    ValB             += stepB;                              // An incident parameter
-    G4double ValB2    = ValB*ValB;                          // A working value b^2
-    G4double IPH      = ValB/shNSl;                         // A working value slope
-    G4double Integ    = 0.;                                 // Integral over ImpactParam.
-    G4double ValS     = 0.;                                 // Prototype of ImpactParameter
-    for(G4int j=1; j<Npb; j++)                              // Second integration over b
-    {
-      ValS           += stepB; //  back to fm               // Impact parameter GeV^-1
-      if(!i) Tb[j]    = CHIPS_Tb(A,ValS/f2Gm1)/f22G2;       // GeV^2, calculate only once
-      //if(!i) Tb[j]    = Thickness(A,ValS/f2Gm1,rAfm)/f22G2; // Calculate T(b) only once
-      G4double FunS   = IPH*ValS;                           // b1*b2/slope
-      if(FunS > Lim) break;                                 // To avoid NAN
-      Integ          += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Tb[j]; // BessI0
-    } 
-    G4double InExp    = -hTotG2*Integ*stepB/dHS;            // Integrated absorption 
-    G4double expB     = std::exp(InExp);                    // Exponential absorption
-    G4double HRIE     = shNReIm*InExp;                      // Phase shift
-    G4double J0qb      = QJ0(std::sqrt(Q2)*ValB)*ValB;      // Bessel0(Q*b)
-    ReSum             += J0qb*(1.-expB*std::cos(HRIE));
-    ImSum             += J0qb*expB*std::sin(HRIE);
-  } 
-  //G4cout<<"CHPS:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl;
-  //return (ReSum*ReSum+ImSum*ImSum)*mb2G2*CMMom*CMMom*stepB*stepB/twopi;
-  //return (ReSum*ReSum+ImSum*ImSum)*f22G2*CMMom*CMMom*stepB*stepB/twopi; // ds/do
-  return (ReSum*ReSum+ImSum*ImSum)*f22G2*stepB*stepB/12; // ds/dt
-}
-
-// Separate quasielastic calculation
-G4double CHIPSDifQuasiElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2)
-{//      =================================================================================
-  static const G4int Npb      = 500;                        // A#of intergation points
-  static const G4double hc2   = .3893793;                   // Transform from GeV^-2 to mb
-  static const G4double mb2G2 = 1./hc2;                     // Transform from mb to GeV^-2
-  //static const G4double mb2G2 = 2.568;                     // Transform from mb to GeV^-2
-  static const G4double f22mb = 10;                         // Transform from fermi^2 to mb
-  static const G4double f22G2 = f22mb*mb2G2;                // Transform from fm2 to GeV^-2
-  static const G4double f2Gm1 = std::sqrt(f22G2);           // Transform from fm to GeV^-1
-  G4double p          = hMom/1000.;                         // Momentum (GeV, inputParam)
-  G4double p2         = p*p;
-  G4double Q2         = aQ2/1000000.;                       // -t in GeV
-  // @@ Temporary only for nucleons
-  G4double p3         = p2*p;
-  G4double p4         = p2*p2;
-  G4double sp         = std::sqrt(p);
-  G4double p2s        = p2*sp;
-  G4double ap         = std::log(p);
-  G4double dl         = ap-3.;
-  G4double dl2        = dl*dl;
-  G4double shPPTot    = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb
-  G4double shPNTot    = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p
-                        +(38.+.3*dl2+8./p)/(1.+1.2/p3);     // SigPP in mb
-  G4double shNTot     = (Z*shPPTot+(A-Z)*shPNTot)/A;        // SighN in mb
-  G4double shPPSl     = 8.*std::pow(p,.055)/(1.+3.64/p4);   // PPslope in GeV^-2
-  G4double shPNSl     = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2
-  G4double shNSl      = (Z*shPPSl+(A-Z)*shPNSl)/A;          // B-slope in GeV^-2
-  G4double shNReIm    = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit
-  // @@ End of temporary ^^^^^^^
-  G4double Re         = std::pow(A,.33333333);              // A-dep coefficient
-  G4double stepB      = (Re+Re+2.7)*f2Gm1/(Npb-1);          // in GeV^-1, step of integral
-  G4double InQE       = 0.;                                 // Quasielastic integral
-  G4double ValB       = -stepB;                             // Impact parameter
-  G4double hTotG2     = shNTot*mb2G2;                       // sigma_hN in GeV^-2
-  for(G4int k=0; k<Npb; k++)
-  {
-    ValB        += stepB;
-    //G4double Tb  = Thickness(A,ValB/f2Gm1,rAfm)/f22G2;      // Calculate T(b) only once
-    G4double Tb  = CHIPS_Tb(A,ValB/f2Gm1)/f22G2;            // GeV^2, calculate only once
-    InQE        += Tb*ValB*std::exp(-hTotG2*Tb);
-  }
-  //G4cout<<"CHPS:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2="
-  //      <<Q2<<",m="<<mb2G2<<G4endl;
-  return InQE*stepB*hTotG2*hTotG2*(1.+shNReIm*shNReIm)*std::exp(-shNSl*Q2)/8/mb2G2;
-}
-#endif
-
-int main()
-{
-#ifdef bestest
-  // *** Test of CHIPS implementation of Bessel functions (accuracy 1.E-8 and better) ***
-  const G4int nj=21;
-  const G4int ny=16;
-  const G4int ni=20;
-  const G4int nk=20;
-  const G4double jX[nj]=
-          				{-5.,-4.,-3.,-2.,-1.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.};
-  const G4double yX[ny]={.1,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.};
-  const G4double iX[ni]=
-                    {0.,.2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,6.,8.,10.};
-  const G4double kX[nk]=
-                    {.1,.2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,6.,8.,10.};
-  const G4double vBJ0[nj]=
-				              {-.1775968,-.3971498,-.2600520,0.2238908,0.7651976,1.0000000,0.7651977,
-                   0.2238908,-.2600520,-.3971498,-.1775968,0.1506453,0.3000793,0.1716508,
-                   -.0903336,-.2459358,-.1711903,0.0476893,0.2069261,0.1710735,-.0142245};
-  const G4double vBJ1[nj]=
-              				{0.3275791,0.0660433,-.3390590,-.5767248,-.4400506,0.0000000,0.4400506,
-                   0.5767248,0.3390590,-.0660433,-.3275791,-.2766839,-.0046828,0.2346364,
-                   0.2453118,0.0434728,-.1767853,-.2234471,-.0703181,0.1333752,0.2051040};
-  const G4double vBY0[ny]=
-        {-1.534239,0.0882570,.51037567,.37685001,-.0169407,-.3085176,-.2881947,-.0259497,
-         0.2235215,0.2499367,0.0556712,-.1688473,-.2252373,-.0782079,0.1271926,0.2054743};
-  const G4double vBY1[ny]=
-    				{-6.458951,-.7812128,-.1070324,0.3246744,0.3979257,0.1478631,-.1750103,-.3026672,
-         -.1580605,0.1043146,0.2490154,0.1637055,-.0570992,-.2100814,-.1666448,0.0210736};
-  const G4double vBI0[ni]={1.0000000,1.0100250,1.0404018,1.0920453,1.1665149,
-                           1.2660658,1.3937256,1.5533951,1.7499807,1.9895593,
-                           2.2795852,3.2898391,4.8807925,7.3782035,11.301922,
-                           17.481172,27.239871,67.234406,427.56411,2815.7167};
-  const G4double vBI1[ni]={.00000000,.10050083,.20402675,.31370403,.43286480,
-                           .56515912,.71467794,.88609197,1.0848107,1.3171674,
-                           1.5906369,2.5167163,3.9533700,6.2058350,9.7594652,
-                           15.389221,24.335643,61.341937,399.87313,2670.9883};
-  const G4double vBK0[nk]={2.4270690,1.7527038,1.1145291,.77752208,.56534710,
-                           .42102445,.31850821,.24365506,.18795475,.14593140,
-                           .11389387,.062347553,.0347395,.019598897,.011159676,
-                           .0063998572,.0036910983,.0012439943,.00014647071,.000017780062};
-  const G4double vBK1[nk]={9.8538451,4.7759725,2.1843544,1.3028349,.86178163,
-                           .60190724,.43459241,.32083589,.24063392,.18262309,
-																											.13986588,.073890816,.040156431,.022239393,.012483499,
-                           .0070780949,.0040446134,.0013439197,.00015536921,.000018648773};
-  //Bessel J0('J',0); // CLHEP J/Y is possible insead of 'J'
-  //Bessel J1('J',1);
-  G4QBesIKJY QI0(BessI0); // CHIPS I/K/J/Y 0/1 Bessel functions
-  G4QBesIKJY QI1(BessI1);
-  G4QBesIKJY QJ0(BessJ0);
-  G4QBesIKJY QJ1(BessJ1);
-  G4QBesIKJY QK0(BessK0);
-  G4QBesIKJY QK1(BessK1);
-  G4QBesIKJY QY0(BessY0);
-  G4QBesIKJY QY1(BessY1);
-  G4cout<<"G4GCS: ***J0*** Test: ================================================"<<G4endl;
-  for(G4int j0=0; j0<nj; j0++)
-		{
-    G4double F=QJ0(jX[j0]);
-    G4double d=F-vBJ0[j0];
-    G4double r=std::abs(d/F);
-    G4cout<<"G4GCS: x="<<jX[j0]<<", J0="<<F<<" - "<<vBJ0[j0]<<" = "<<d<<", r="<<r<<G4endl;
-  }
-  G4cout<<"G4GCS: ***J1*** Test: ================================================"<<G4endl;
-  for(G4int j1=0; j1<nj; j1++)
-		{
-    G4double F=QJ1(jX[j1]);
-    G4double d=F-vBJ1[j1];
-    G4double r=std::abs(d/F);
-    G4cout<<"G4GCS: x="<<jX[j1]<<", J1="<<F<<" - "<<vBJ1[j1]<<" = "<<d<<", r="<<r<<G4endl;
-  }
-  G4cout<<"G4GCS: ***Y0*** Test: ================================================"<<G4endl;
-  for(G4int y0=0; y0<ny; y0++)
-		{
-    G4double F=QY0(yX[y0]);
-    G4double d=F-vBY0[y0];
-    G4double r=std::abs(d/F);
-    G4cout<<"G4GCS: x="<<yX[y0]<<", Y0="<<F<<" - "<<vBY0[y0]<<" = "<<d<<", r="<<r<<G4endl;
-  }
-  G4cout<<"G4GCS: ***Y1*** Test: ================================================"<<G4endl;
-  for(G4int y1=0; y1<ny; y1++)
-		{
-    G4double F=QY1(yX[y1]);
-    G4double d=F-vBY1[y1];
-    G4double r=std::abs(d/F);
-    G4cout<<"G4GCS: x="<<yX[y1]<<", Y1="<<F<<" - "<<vBY1[y1]<<" = "<<d<<", r="<<r<<G4endl;
-  }
-  G4cout<<"G4GCS: ***I0*** Test: ================================================"<<G4endl;
-  for(G4int i0=0; i0<ni; i0++)
-		{
-    G4double F=QI0(iX[i0]);
-    G4double d=F-vBI0[i0];
-    G4double r=std::abs(d/F);
-    G4cout<<"G4GCS: x="<<iX[i0]<<", I0="<<F<<" - "<<vBI0[i0]<<" = "<<d<<", r="<<r<<G4endl;
-  }
-  G4cout<<"G4GCS: ***I1*** Test: ================================================"<<G4endl;
-  for(G4int i1=0; i1<ni; i1++)
-		{
-    G4double F=QI1(iX[i1]);
-    G4double d=F-vBI1[i1];
-    G4double r=std::abs(d/F);
-    G4cout<<"G4GCS: x="<<iX[i1]<<", I1="<<F<<" - "<<vBI1[i1]<<" = "<<d<<", r="<<r<<G4endl;
-  }
-  G4cout<<"G4GCS: ***K0*** Test: ================================================"<<G4endl;
-  for(G4int k0=0; k0<nk; k0++)
-		{
-    G4double F=QK0(kX[k0]);
-    G4double d=F-vBK0[k0];
-    G4double r=std::abs(d/F);
-    G4cout<<"G4GCS: x="<<kX[k0]<<", I1="<<F<<" - "<<vBK0[k0]<<" = "<<d<<", r="<<r<<G4endl;
-  }
-  G4cout<<"G4GCS: ***K1*** Test: ================================================"<<G4endl;
-  for(G4int k1=0; k1<nk; k1++)
-		{
-    G4double F=QK1(kX[k1]);
-    G4double d=F-vBK1[k1];
-    G4double r=std::abs(d/F);
-    G4cout<<"G4GCS: x="<<kX[k1]<<", I1="<<F<<" - "<<vBK1[k1]<<" = "<<d<<", r="<<r<<G4endl;
-  }
-  G4cout<<"End .................................................................."<<G4endl;
-#endif
-
-  // Test of integrated cross sections
-  //const G4int na=12;
-  const G4int na=1;
-  //
-  const G4int np=7;
-  const G4int nm=1;
-  ////                He Be C O  Al Ti Ni Cu  Sn Ta  Pb   U
-  //const G4int A[na]={4,9,12,16,27,48,58,64,120,181,207,238}; // A's of target nuclei
-  //                 He Al Pb
-  const G4int A[na]={119}; // A's of target nuclei
-  //                     p    n  pi+  pi-  K+  K-  antip
-  const G4int pdg[np]={2212,2112,211,-211,321,-321,-2212}; // projectiles
-  const G4double mom[nm]={120.}; // momentum in MeV/c
-#ifdef integrc
-  for(G4int ip=0; ip<np; ip++)
-		{
-    G4int PDG=pdg[ip];
-    for(G4int im=0; im<nm; im++)
-		  {
-      CalculateParameters(PDG, mom[im]);
-      G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<":A,Tot,Inel,Prod,El,QEl"<<G4endl;
-      for(G4int ia=0; ia<na; ia++)
-		    {
-        CalculateIntegralCrossSections(A[ia]);
-        G4cout<<A[ia]<<" "<<TotalCrSec<<" "<<InelCrSec<<" "<<ProdCrSec
-              <<" "<<ElasticCrSec<<" "<<QuasyElasticCrSec<<G4endl;
-      }
-    }
-  }
-#endif
-
-#ifdef eldiffcs
-  const G4double ms=.000001;
-  const G4int nt=200;
-  G4double t[nt]; // Q2=-t in Mev^2
-  const G4double tMin=200.;
-  const G4double tMax=2000000.;
-  const G4double ltMi=std::log(tMin);
-  const G4double ltMa=std::log(tMax);
-  const G4double dlt=(ltMa-ltMi)/(nt-1);
-  G4double lt=ltMi-dlt;
-  for(G4int it=0; it<nt; it++)
-		{
-    lt+=dlt;
-    t[it]=std::exp(lt);
-  }
-  ////                He Be C O  Al Ti Ni Cu  Sn Ta  Pb   U
-  //const G4int Z[na]={2,4, 6, 8,13,22,28,29, 50, 73, 82, 92}; // Z's of target nuclei
-  //                He Al Pb
-  const G4int Z[na]={50}; // Z's of target nuclei
-  // Test of differential ellastic cross sections
-		Initialize();
-  //for(G4int ip=0; ip<np; ip++)
-  for(G4int ip=0; ip<1; ip++)
-		{
-    G4int PDG=pdg[ip];
-    for(G4int im=0; im<nm; im++)
-		  {
-      CalculateParameters(PDG, mom[im]);
-      G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<G4endl;
-      for(G4int ia=0; ia<na; ia++)
-		    {
-        G4cout<<"G4GCS:-------------------A="<<A[ia]<<G4endl;
-        for(G4int it=0; it<nt; it++)
-		      {
-          G4double Sig1 = CHIPSDiffElCS(PDG, mom[im], Z[ia], A[ia], t[it]);
-          //G4double Sig1 = DiffElasticCS(PDG, mom[im], A[ia], t[it]); //Doesn't work
-          G4double Sig2 = CoherentDifElasticCS(PDG, mom[im], A[ia], t[it]);
-          G4double Sig3 = CHIPSDifElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]);
-          G4double CQEl = CHIPSDifQuasiElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]);
-								  G4cout<<std::setw(12)<<ms*t[it]<<std::setw(12)<<Sig1<<std::setw(12)<<Sig2
-                <<std::setw(12)<<InCoh<<std::setw(12)<<Sig3<<std::setw(12)<<CQEl<<G4endl;
-        }
-      }
-    }
-  }
-#endif
-  return 0;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4MultyFragmentationTest.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4MultyFragmentationTest.cc	(revision 968)
+++ 	(revision )
@@ -1,354 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//#define debug
-//#define pdebug
-
-//#include "G4UIterminal.hh"
-#include "G4ios.hh"
-#include "G4QEnvironment.hh"
-
-int main()
-{
-  //G4double theEnergyLossPerFermi = 1.*GeV;
-  G4int nop = 152;                               // clusters (A<6)
-  G4double fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40)
-  G4double fractionOfPairedQuasiFreeNucleons = 0.05;
-  G4double clusteringCoefficient = 5.;
-  G4double temperature = 180.;
-  G4double halfTheStrangenessOfSee = 0.3; // = s/d = s/u
-  G4double etaToEtaPrime = 0.3;
-  G4double fusionToExchange = 100.;
-  //G4double theInnerCoreDensityCut = 50.;
-#ifdef debug
-	 G4cout<<"G4MulFragTst: nop="<<nop<<",fN="<<fractionOfSingleQuasiFreeNucleons
-        <<",fD="<<fractionOfPairedQuasiFreeNucleons<<",wC="<<clusteringCoefficient
-        <<",vn="<<fusionToExchange<<",T="<<temperature
-        <<",ss="<<halfTheStrangenessOfSee<<",se="<<etaToEtaPrime<<G4endl;
-#endif
-  //-------- Initialize CHIPS
-  G4QCHIPSWorld* theW=G4QCHIPSWorld::Get();
-  theW->GetParticles(nop);           // Create CHIPS World of nop particles
-  G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
-                            fractionOfPairedQuasiFreeNucleons,
-			                         clusteringCoefficient, 
-					                       fusionToExchange);
-  G4Quasmon::SetParameters(temperature, halfTheStrangenessOfSee, etaToEtaPrime);
-
-  G4int tPDG = 90006005;    // PDG of the Target Nucleus
-  G4int nEvents= 100000;
-  //#ifdef debug
-	 G4cout<<"G4MulFragTst: targetPDG="<<tPDG<<", nEv="<<nEvents<<G4endl;
-  //#endif
-  G4QHadronVector projHV;
-
-  G4int PDG1=211;
-  G4LorentzVector LV1(-485.221,-216.252,267.887,611.103);
-  //#ifdef debug
-	 G4cout<<"G4MulFragTst: #1 projPDG="<<PDG1<<", 4M="<<LV1<<G4endl;
-  //#endif
-
-  G4int PDG2=-211;
-  G4LorentzVector LV2(-201.944,475.785,952.224,1092.41);
-  //#ifdef debug
-	 G4cout<<"G4MulFragTst: #2 projPDG="<<PDG2<<", 4M="<<LV2<<G4endl;
-  //#endif
-
-  G4int PDG3=111;
-  G4LorentzVector LV3(213.293,183.683,996.566,1044.32);
-  //#ifdef debug
-	 G4cout<<"G4MulFragTst: #3 projPDG="<<PDG3<<", 4M="<<LV3<<G4endl;
-  //#endif
-
-  //std::ifstream inFile("chipstest.in", std::ios::in);
-  //G4double temperature;
-  //G4double ssin2g;
-  //G4double eteps;
-  //G4double momb;
-  //G4double enb;
-  //G4double cP;
-  //G4double fN;
-  //G4double fD;
-  //G4double rM;
-  //G4double sA;
-  //G4int    nop;
-  //G4int    pPDG;
-  //G4int    tPDG;
-  //G4int    nEvt;
-  ////G4int    nofdecays;
-  ////G4int    decmask=0;
-  //inFile>>temperature>>ssin2g>>eteps>>nop>>momb>>enb>>pPDG>>tPDG>>nEvt>>fN>>fD>>cP>>rM>>sA;
-  //G4cout<<"CHIPStest Par's: Temp="<<temperature<<",SSs="<<ssin2g<<",P/S="<<eteps<<",nop="
-  //      <<nop<<",p="<<momb<<",e="<<enb<<",projCode="<<pPDG<<",targCode="<<tPDG<<",nEv="
-  //      <<nEvt<<",fN="<<fN<<",fD="<<fD<<",cP="<<cP<<",rM="<<rM<<",sA="<<sA<<G4endl;
-  //-------- Initialize CHIPS
-  //G4QCHIPSWorld* theW=G4QCHIPSWorld::Get();
-  //theW->GetParticles(nop);           // Create CHIPS World of nop particles
-  //G4Exception("***CHIPStest: TMP");
-  //G4QNucleus::SetParameters(fN,fD,cP,rM);
-  //G4Quasmon::SetParameters(temperature,ssin2g,eteps);
-  //G4QEnvironment::SetParameters(sA);                 // SolAngle (pbar-A secondary capture)
-  //-------- write results onto a file --------
-  //std::ofstream outFile( "chipstest.out", std::ios::out);
-  //outFile.setf( std::ios::scientific, std::ios::floatfield );
-
-  //--------- Example how to write in the outFile -----------
-  //outFile << " " << endl;
-
-  // *********** Now momb is a momentum of the incident particle *************
-  //G4double mp=G4QPDGCode(pPDG).GetMass();             // @@ just for the check
-  //G4QContent pQC=G4QPDGCode(pPDG).GetQuarkContent();
-  //G4int    cp=pQC.GetCharge();
-  //G4int    bnp=pQC.GetBaryonNumber();
-  //G4double mp1=G4QPDGCode(PDG1).GetMass();             // @@ just for the check
-  G4QContent pQC1=G4QPDGCode(PDG1).GetQuarkContent();
-  G4int    cp1=pQC1.GetCharge();
-  G4int    bnp1=pQC1.GetBaryonNumber();
-  //G4double mp2=G4QPDGCode(PDG2).GetMass();             // @@ just for the check
-  G4QContent pQC2=G4QPDGCode(PDG2).GetQuarkContent();
-  G4int    cp2=pQC2.GetCharge();
-  G4int    bnp2=pQC2.GetBaryonNumber();
-  //G4double mp3=G4QPDGCode(PDG3).GetMass();             // @@ just for the check
-  G4QContent pQC3=G4QPDGCode(PDG3).GetQuarkContent();
-  G4int    cp3=pQC3.GetCharge();
-  G4int    bnp3=pQC3.GetBaryonNumber();
-  //momb=momb;
-  //G4double ep=sqrt(mp*mp+momb*momb);                  // @@ just for the check
-  //if(enb>0.) ep=enb;
-  G4double mt=G4QPDGCode(tPDG).GetMass();             // @@ just for the check
-  G4QContent tQC=G4QPDGCode(tPDG).GetQuarkContent();
-  G4int    ct=tQC.GetCharge();
-  G4int    bnt=tQC.GetBaryonNumber();
-  G4int    totC=cp1+cp2+cp3+ct;
-  G4int    totBN=bnp1+bnp2+bnp3+bnt;
-  G4LorentzVector preSumLV(0.,0.,0.,mt);
-  preSumLV+=LV1+LV2+LV3;
-  G4double fEvt=nEvents;
-  G4double sumE=0.;
-  G4double sumK=0.;
-  G4double sumG=0.;
-  G4double sumT=0.;
-  G4double sumN=0.;
-  G4double sum0=0.;
-  G4double sumP=0.;
-  G4double sum1N=0.;
-  G4double sumNN=0.;
-  G4double sumPP=0.;
-  G4double sumAL=0.;
-  G4double time=clock()/CLOCKS_PER_SEC;
-  // Main LOOP over events ======================================
-  for (G4int ir=0; ir<nEvents; ir++)
-  {
-    // Randomization loop: cycle random generator, using 2 lower digits in nEvents
-	   G4int    iRandCount = nEvents%100;
-	   G4double vRandCount = 0.;
-	   while (iRandCount>0)
-	   {
-	     vRandCount = G4UniformRand();
- 	    iRandCount--;
-	   }
-    if(!(ir%1000) && ir) G4cout<<"G4MultFragTst: "<<ir<<" events are simulated"<<G4endl;
-    //G4cout<<"G4MultFragTst: "<<ir<<" events are simulated"<<G4endl;
-    G4LorentzVector totSum = preSumLV;
-    G4QHadron* H1 = new G4QHadron(PDG1,LV1);
-    G4QHadron* H2 = new G4QHadron(PDG2,LV2);
-    G4QHadron* H3 = new G4QHadron(PDG3,LV3);
-    G4int           totCharge = totC;
-    G4int           totBaryN = totBN;
-    G4QHadronVector projHV;
-    projHV.push_back(H1);                                 // DESTROYED over 3 line
-    projHV.push_back(H2);                                 // DESTROYED over 2 line
-    projHV.push_back(H3);                                 // DESTROYED over 1 line
-    G4QEnvironment* pan= new G4QEnvironment(projHV,tPDG); // DELETED over 8 lines
-    std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
-    projHV.clear();
-#ifdef debug
-    G4cout<<"CHIPStest:===>>> Now call Fragment (HadronizeQuasmon) function" << G4endl;
-#endif
-    G4QHadronVector* output; // Prototype of the output
-    try
-    {
-      output = pan->Fragment();// DESTROYED in the end of the LOOP work space
-    }
-    catch (G4QException& error)
-    {
-#ifdef pdebug
-      G4cout<<"***CHIPStest: Exception is catched"<<G4endl;
-#endif
-      G4cerr<<"***CHIPStest Abort: "<<error.GetMessage()<<G4endl;
-      abort();
-    }
-#ifdef debug
-    G4cout<<"CHIPStest:--->>>Now come out of Fragment (HadronizeQuasmon) function"<<G4endl;
-#endif
-    delete pan; // Destruct theHadronVector (& theCandidateVector) of the Quasmon
-#ifdef debug
-    G4cout << "CHIPStest: >>> Here the histograms are filled" << G4endl;
-#endif
-    G4int tNH = output->size();
-    G4int npt=0;
-    G4int nGamma=0;
-    G4double EGamma=0;
-    G4int nP0=0;
-    G4int nPP=0;
-    G4int nPN=0;
-    G4int nKaons=0;
-    G4int nEta=0;
-    G4int nAlphas=0;
-    G4int nPhotons=0;
-    G4int nProtons=0;
-    G4int nNeutrons=0;
-    G4int nSpNeut=0;
-    G4int nSpAlph=0;
-    G4int nOmega=0;
-    G4int nDec=0;
-    G4int dirN=0;
-#ifdef pdebug
-    G4cout<<"----------DONE^^^^^^^************^^^^^^^^^^^:ir="<<ir<<": #ofH="<<tNH<<G4endl;
-    if(!(ir%100)) G4cerr<<"#"<<ir<<G4endl;
-#endif
-    G4bool alarm=false;
-    G4bool rad=false;
-    G4bool hyp=false;
-    G4bool badPDG=false;
-    for (G4int ind=0; ind<tNH; ind++)
-    {
-      G4QHadron* curH=output->operator[](ind);
-      G4double      m = curH->GetMass();            // Mass of the particle
-      G4LorentzVector lorV = curH->Get4Momentum();  // 4-momentum of the particle
-      if(std::fabs(m-lorV.m())>.005)
-	     {
-		      G4cerr<<"***CHIPStest: m="<<lorV.m()<<" # "<<m<<", d="<<lorV.m()-m<<G4endl;
-        alarm=true;
-	     }
-      if(!(lorV.e()>=0||lorV.e()<0)   || !(lorV.px()>=0||lorV.px()<0) ||
-         !(lorV.py()>=0||lorV.py()<0) || !(lorV.pz()>=0||lorV.pz()<0))
-	     {
-		      G4cerr<<"***CHIPStest: NAN in LorentzVector="<<lorV<<G4endl;
-        alarm=true;
-	     }
-      G4int         d=curH->GetNFragments(); // In how many particles this particle decayed
-      G4ThreeVector p = lorV.vect();         // 3-momentum of the particle
-      G4double      e = lorV.e();            // Energy of the particle
-      G4int         c=curH->GetPDGCode();    // PDG Code of the particle
-      //if(!d&&(c==90000002||c==90002000||c==92000000||c==221||c==331))
-      if(!d&&(c==90000002||c==90002000||c==92000000))
-      {
-        //G4cout<<"***CHIPStest:***Dibaryon or Eta*** ind="<<ind<<", PDG="<<c<<G4endl;
-        G4cout<<"***CHIPStest:***Dibaryon *** ind="<<ind<<", PDG="<<c<<G4endl;
-        alarm=true;
-      }
-      if(!d&&(c==90000003||c==90003000||c==93000000))
-      {
-        G4cout<<"***CHIPStest:***Tribaryon *** ind="<<ind<<", PDG="<<c<<G4endl;
-        alarm=true;
-      }
-      if(!d) npt++;
-      if(d) nDec+=d;
-      if(c==223) nOmega++;
-      if(c==22) nPhotons++;
-      if(c==311||c==321||c==-311||c==-321) nKaons++; // kaons
-      if(c==221) nEta++;                             // etas
-      if(c==90002002) nAlphas++;                     // Alphas
-      if(c==2212 || c==90001000) nProtons++;         // Protons
-      if(c==90000001 || c==90001000) dirN++;         // Ditrect nucleons
-      if(c==2112 || c==90000001) nNeutrons++;        // Neutrons
-      if((c==2112 || c==90000001) && std::fabs(e-1005.)<3.) nSpNeut++;// Dibar-Neutrons
-      if(!d && c==90002002 && e-m<7.) nSpAlph++;     // Special Alphas
-      if(c==111) nP0++;                              // Neutral  pions
-      if(c==-211) nPN++;                             // Negative pions
-      if(c==211) nPP++;                              // Positive pions
-      if(c==22) nGamma++;                            // Gammas
-      if(c==22) EGamma+=e;                           // Energy of gammas
-      if(!d) totCharge-=curH->GetCharge();
-      if(!d) totBaryN-=curH->GetBaryonNumber();
-      if(!d) totSum    -= lorV;
-      if(c>80000000 && (c<90000000 || c%1000>500 || c%1000000>500000) ||
-         !(c>=0 || c<0))
-      {
-        G4cout<<"***OUTPUT ERROR*** CHIPStest: bad PDG is found. It is "<<c<<G4endl;
-        badPDG=true;
-	     }
-      rad=rad||c==90002000||c==90003000||c==90004000||c==90000002||c==90000003||c==90000004
-             ||c==90002003||c==90003002||c==90004002||c==90002005||c==90005002||c==90004004
-             ||c==90006002;
-      hyp = hyp || c>90999999;
-#ifdef pdebug
-      G4cout<<"#"<<ind<<"(d="<<d<<"), PDG="<<c<<",4M="<<lorV<<m<<",T="<<lorV.e()-m<<G4endl;
-#endif
-    }
-#ifdef pdebug
-    G4cout<<"CHECK: 4M="<<totSum<<", Charge="<<totCharge<<", BaryN="<<totBaryN<<G4endl;
-#endif
-    G4double ss=std::fabs(totSum.t())+std::fabs(totSum.x())+std::fabs(totSum.y())+std::fabs(totSum.z());    
-	   if (totCharge ||totBaryN || !(ss<.01) || alarm || nGamma&&!EGamma || badPDG)
-	   //if (totCharge || ss>.01 || alarm || nSpNeut)
-	   //if (totCharge || ss>.01 || alarm || nSpAlph)
-    {
-      G4cerr<<"***CHIPStest:#"<<ir<<":n="<<tNH<<",4M="<<totSum<<",Charge="<<totCharge
-            <<",BaryN="<<totBaryN<<G4endl;
-      if(nGamma&&!EGamma)G4cerr<<"***CHIPStest: Egamma=0"<<G4endl;
-      totSum    = preSumLV;
-      for (int indx=0; indx<tNH; indx++)
-      {
-        G4QHadron* curH=output->operator[](indx);
-        G4double m = curH->GetMass();
-        G4LorentzVector lorV = curH->Get4Momentum();
-        G4int d=curH->GetNFragments();
-        G4int c=curH->GetPDGCode();
-        if(!d) totSum    -= lorV;
-        G4cerr<<"#"<<indx<<"("<<d<<"), PDG="<<c<<", m/LV="<<m<<lorV<<", T="<<lorV.e()-m
-              <<", d4M="<<totSum<<G4endl;
-      }
-      G4Exception("***CHIPStest: ALARM or charge/energy/momentum is not conserved");
-    }
-    sum0+=nP0;
-    sumP+=nPP;
-    sumN+=nPN;
-    G4int nPions=nP0+nPP+nPN;
-    sumG+=nGamma;
-    sumT+=EGamma;
-    if(nAlphas)sumAL++;
-    if(nProtons)sumPP++;
-    if(nNeutrons)sumNN++;
-    if(nNeutrons&&!nProtons&&!nAlphas) sum1N++;
-    if(nKaons)sumK++;
-    if(nEta)sumE++;
-    if (nPhotons) nPions=11;
-    if (nKaons==2&&npt==2) nPions=1;
-    else if (nKaons) nPions=10;
-    //histPi.fill(nPions);
-    //histNeut.fill(nNeutrons);
-    std::for_each(output->begin(), output->end(), DeleteQHadron());
-    output->clear();
-    delete output;
-  }
-  time=(clock()/CLOCKS_PER_SEC-time)/fEvt;
-  G4cerr<<"CHIPStest::t="<<time<<",Yields:pi-="<<sumN/fEvt<<",pi+="<<sumP/fEvt<<",pi0="
-        <<sum0/fEvt<<",K="<<sumK/fEvt<<",eta="<<sumE/fEvt<<",gamma="<<sumG/fEvt<<"(<E>="
-        <<sumT/sumG<<"),n="<<sumNN/fEvt<<",p="<<sumPP/fEvt<<",alpha="<<sumAL/fEvt
-        <<",onlyN="<<sum1N/fEvt<<G4endl;
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuElXQ2Integration.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuElXQ2Integration.cc	(revision 968)
+++ 	(revision )
@@ -1,290 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuE DIS (x,Q2) approximation is integrated over x
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 30-Sept-05
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-void strucf(int A, double x, double Q2, double& f2, double& xf3, double& fL)
-{
-  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN2=MN*MN;       // M_N^2 in GeV^2
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double Wt=MN+mpi;       // Delta threshold
-  //const double W2t=Wt*Wt;       // Squared Delta threshold
-  const Genfun::LogGamma lGam;
-  static int mA=0;
-  static double mQ2=0., mN, mD, mDel, mU2, mU3, mV;
-  //static double mUU;
-  double N=3., D=0., Del=0., U2=0., U3=0., V=0.;
-  //double UU=0.;
-  if(A==mA && Q2==mQ2)         // Associative memory for acceleration
-		{
-    N  =mN;
-    D  =mD;
-    Del=mDel;
-    U2 =mU2;
-    U3 =mU3;
-    //UU =mUU;
-    V  =mV;
-  }
-  else
-  {
-		  //double r=0.;
-    double min=.077;             // Delta(0)
-				double max=.47;              // Hard Pomeron
-    //double H=1.22;
-    if(Q2<=0.) G4cout<<"strucf: Q2="<<Q2<<" <= 0"<< G4endl;
-    double Q=std::sqrt(Q2);
-    //if(A==1)                   // Proton
-		  //{
-    //  r=std::sqrt(Q2/1.66);
-    //  max=.5;
-    //}
-    //else if(A<13)              // Light nuclei
-		  //{
-    //  double f=Q2/4.62;
-    //  r=f*f;
-    //  max=.3;
-    //  if(A>2) H=1.;
-    //}
-    //else if(A>0)               // Heavy nuclei
-		  //{
-    //  double f=Q2/3.4;
-    //  double ff=f*f;
-    //  r=ff*ff;
-    //  max=.5;
-    //  H=1.;
-    //}
-    //else G4cout<<"strucf: A="<<A<<" <= 0"<< G4endl;
-    //
-    N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
-    Del=min+(max-min)/(1.+5./Q);
-    D=0.68*std::pow(1.+.145/Q2,-1.-Del); // 0.68=1-0.34, m2=.145 GeV2
-    V=3*(1.-D)*(N-1.);
-    double c3=.75;                       // 3*0.25
-    double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
-    U2=(c3+N-3.)*uu;
-    U3=c3*uu;
-    mA  = A;
-    mQ2 = Q2;
-    mN  = N;
-    mD  = D;
-    mDel=Del;
-    mU2 =U2;
-    mU3 =U3;
-    //mUU =UU; // @@
-    mV  =V;
-  }  
-  // From here the Q2 coefficients are used
-  double x1=std::pow(1.-x,N-2.);
-  double pp=D*std::pow(x,-Del)*x1;
-  double dir=V*x*x1;
-  double per=U2*pp;
-		f2 = per + dir;
-  //double W2=MN2-MN2*x+Q2/x-Q2;
-  //if(W2<W2t)
-		//{
-  //  per=UU*pp;
-  //  xf3= per+dir;
-  //}
-  //else
-  xf3= U3*pp+dir;
-  fL = per/5.;                             // 20%
-  return;
-}
-
-void getFun(int A, double lx, double Q2, double* f)
-{
-  double f2=0., xf3=0., fL=0.;
-  if (lx>0.5) G4cerr<<"***getFun: ln(x)="<<lx<<">.5"<<G4endl;
-  double x=std::exp(lx);
-  double x2=x*x;
-  strucf(A, x, Q2, f2, xf3, fL);
-  f[0]=f2;             // direct part
-  f[1]=(-f2+xf3)/x;    // *y (neutrino) part
-  f[2]=(-f2-xf3)/x;    // *y (anti-neutrino) part
-  f[3]=(f2-fL-xf3)/x2; // *y2 (neutrino) part
-  f[4]=(f2-fL+xf3)/x2; // *y2 (anti-neutrino) part
-}
-
-int main()
-{
-  const double reps=.001;        // relative accuracy of the total Q2 integral calculation
-  const double xeps=.0001;      // relative accuracy of the total X integral calculation
-  //           =========
-  const double GF=1.16637e-5;   // Fermi constant in GeV^-2
-  const double GF2=GF*GF;       // Squared Fermi constant in GeV^-4
-  const double MW=80.425;       // Mass of W-boson in GeV
-  const double MW2=MW*MW;       // Squared mass of W-boson in GeV^2
-  const double MW4=MW2*MW2;     // Quadro mass of W-boson in GeV^4
-  const double hc2=38937932300.;// (hc)^2 in GeV^2*10^-38cm2 to convert GeV^-2 to 10^-38cm2
-  const double pif=3.14159265*4;// 4pi
-  const double sik=GF2*hc2/pif; // precalculated coefficient
-  //const double mpi=.1349766;    // pi0 meson mass in GeV
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
-  const double me=.00051099892; // electron mass in GeV
-  const double me2=me*me;       // m_e^2 in GeV^2
-  const double hme2=me2/2;      // .5*m_e^2 in GeV^2
-  //const double mmu=.105658369;  // mu meson mass in GeV
-  //const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
-  //const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
-  //const double mtau=1.777;      // tau meson mass in GeV
-  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
-  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
-  const double mp=.93827203;    // proton mass in GeV
-  const double mn=.93956536;    // neutron mass in GeV
-  //const double md=1.87561282;   // deuteron mass in GeV
-  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
-  const double MD=1.232;        // proton mass in GeV
-  //const double mp2=mp*mp;       // m_p^2 in GeV^2
-  const double MN2=MN*MN;       // M_N^2 in GeV^2
-  const double dMN=MN+MN;       // 2*M_N in GeV
-  const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
-  const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
-  const double EminE=me+me2/dMN;// Threshold for muon production
-  //const double EminMu=mmu+mmu2/dMN;  // Threshold for muon production
-  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
-  //
-  //const double mc=.3;           // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  const double mc=std::sqrt((MN*MN+MD*MD)/2)-MN; // Squared mean between N and \Delta
-  //const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
-  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
-  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
-  // _____ Begin of Test Area
-  //Genfun::LogGamma logGamma;
-  //double n=4.9;
-		//double g=exp(logGamma(n));
-  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
-  // ^^^^^ End of Test Area
-  //
-  double f[5];                    // A working array
-  int    A=12;                    // Neucleus for which calculations should be done
-  double lEnuMin=0;               // LogLog of Minimum energy of neutrino
-  double lEnuMax=std::log(1.+std::log(300./EminE)); // LogLog of MaximumEnergy of neutrino
-  int    nE=49;                   // Number of points
-  double dlE=(lEnuMax-lEnuMin)/nE;
-  lEnuMin+=dlE/10;
-  lEnuMax+=dlE/5;
-  G4cout<<"Emin="<<EminE<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
-  for(double lEnu=lEnuMin; lEnu<lEnuMax; lEnu+=dlE)
-		{
-    double Enu=std::exp(std::exp(lEnu)-1.)*EminE; // Energy of neutrino/anti-neutrino
-    double dEnu=Enu+Enu;          // doubled energy of nu/anu
-    double Enu2=Enu*Enu;          // squared energy of nu/anu
-    double Ee=Enu-me;           // Free Energy of neutrino/anti-neutrino
-    double Ee2=Ee*Ee;          // squared energy of nu/anu
-    double ME=Enu*MN;             // M*E
-    double dME=ME+ME;             // 2*M*E
-    double DIStsig=1.;            // Total curent DIS cross-section to be integrated
-    double DISmsig=1.e20;         // Total remembered DIS cross-section
-    double dEMN=(dEnu+MN)*ME;
-    double MEm=ME-hme2;
-    double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
-    double E2M=MN*Enu2-(Enu+MN)*hme2;
-    double ymax=(E2M+sqE)/dEMN;
-    double ymin=(E2M-sqE)/dEMN;
-    double rmin=1.-ymin;
-    double rhm2E=hme2/Enu2;
-    double Q2min=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
-    double Q2max=dME*ymax;
-    int    nQ2=8;
-    //G4cout<<"*** E="<<Enu<<", Q2i="<<Q2min<<" < Q2a="<<Q2max<<", yi="<<ymin<<" < ya="
-    //      <<ymax<<G4endl;
-    while(std::fabs(DIStsig-DISmsig)/DIStsig>reps)
-    {
-      DISmsig=DIStsig;
-      DIStsig=0.;
-      nQ2*=2;
-      double dQ2=(Q2max-Q2min)/nQ2;
-      for(double Q2=Q2min+dQ2/2; Q2<Q2max; Q2+=dQ2)
-				  {
-        double DISxint=1.;        // Curent DIS x-integral
-        double DISmint=1.e20;     // Remembered DIS x-integral
-        double Q2M=Q2+MW2;
-        double dik=MW4/Q2M/Q2M;
-        double qmc=Q2+mcV;
-        double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2);//QuasielasticBoundary
-        //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c)
-        //double xN=Q2/dME;
-        double xN=Q2/MN/(Ee+std::sqrt(Ee2+Q2));
-        //double lXmin=log(xN/ymax);
-        double lXmin=std::log(xN);
-        // ****** QE ********
-        if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES >>>>>>>>>>>>>>>>
-        // *** End of QE^^^^^
-        double lXmax=0.;               // QES is in DIS
-        //double lXmax=lXQES;          // Cut off quasielastic
-        int    nX=8;
-        while(std::fabs(DISxint-DISmint)/DISxint>xeps)
-        {
-          DISmint=DISxint;
-          DISxint=0.;
-          nX*=2;
-          double dlX=(lXmax-lXmin)/nX;
-          for(double lX=lXmin+dlX/2; lX<lXmax; lX+=dlX)
-				      {
-            getFun(A, lX, Q2, f);
-            // ***** Neutrino/Antineutrino switch ******>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-            DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
-            //DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
-            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
-          }
-          DISxint*=dlX;
-          //G4cout<<"--- E="<<Enu<<" --- Q2="<<Q2<<" --- nX="<<nX<<", iX="<<DISxint
-          //      <<", mX="<<DISmint<<", rX="<<(DISxint-DISmint)/DISxint<<G4endl;
-        }
-        //G4cout<<"(E="<<Enu<<"), Q2="<<Q2<<", I="<<DISxint/dik/dik<<G4endl;
-        DIStsig+=DISxint*dik;
-      }
-      DIStsig*=dQ2;
-      //G4cout<<"=== E="<<Enu<<" ===> nQ="<<nQ2<<", iQ="<<DIStsig<<", mQ="<<DISmsig
-						//      <<", rQ="<<(DIStsig-DISmsig)/DIStsig<<G4endl;
-    }
-    //===== tot/qe choice ====
-    DIStsig*=sik/Enu;
-    //G4cout<<"***total-neutrino*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
-    //G4cout<<"***total-antineutrino*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
-    G4cout<<"***quasiel-nu*** E="<<Enu<<" ,sig/E="<<DIStsig*Enu<<G4endl;
-    //G4cout<<"***quasiel-antinu*** E="<<Enu<<" ,sig/E="<<DIStsig*Enu<<G4endl;
-    //...................
-    //DIStsig*=sik;
-    //G4cout<<"***qelas*** E="<<Enu<<",sig="<<DIStsig<<G4endl;
-    //===== End of the choice
-		} // End of the Enery LOOP
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuElXQ2Integration_N.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuElXQ2Integration_N.cc	(revision 968)
+++ 	(revision )
@@ -1,294 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuEl DIS (x,Q2) approximation is integrated over x & Q2
-// .....................................................[ for a <nucleon> ]
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Oct-07
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-void strucf(int A, double x, double Q2, double& f2, double& xf3, double& fL)
-{
-  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN2=MN*MN;       // M_N^2 in GeV^2
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double Wt=MN+mpi;       // Delta threshold
-  //const double W2t=Wt*Wt;       // Squared Delta threshold
-  const Genfun::LogGamma lGam;
-  static int mA=0;
-  static double mQ2=0., mN, mD, mDel, mU2, mU3, mV;
-  //static double mUU;
-  double N=3., D=0., Del=0., U2=0., U3=0., V=0.;
-  //double UU=0.;
-  if(A==mA && Q2==mQ2)         // Associative memory for acceleration
-		{
-    N  =mN;
-    D  =mD;
-    Del=mDel;
-    U2 =mU2;
-    U3 =mU3;
-    //UU =mUU;
-    V  =mV;
-  }
-  else
-  {
-    double Q=std::sqrt(Q2);
-    N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
-    Del=.077+.393/(1.+5./Q);
-    D=.68*std::pow(1.+.145/Q2,-1.-Del);
-    V=3*(N-1.);
-    double c3=.75;
-    double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
-    U2=(c3+N-3.)*uu;
-    U3=c3*uu;
-    //UU=uu+uu+uu; // @@
-    mA  = A;
-    mQ2 = Q2;
-    mN  = N;
-    mD  = D;
-    mDel=Del;
-    mU2 =U2;
-    mU3 =U3;
-    //mUU =UU; // @@
-    mV  =V;
-  }  
-  // From here the Q2 coefficients are used
-  double x1=std::pow(1.-x,N-2.);
-  double pp=D*std::pow(x,-Del)*x1;
-  double dir=(1-D)*V*x*x1;
-  double per=U2*pp;
-		f2 = per + dir;
-  //double W2=MN2-MN2*x+Q2/x-Q2;
-  //if(W2<W2t)
-		//{
-  //  per=UU*pp;
-  //  xf3= per+dir;
-  //}
-  //else
-  xf3= U3*pp+dir;
-  fL = per/5.;
-  return;
-}
-
-void getFun(int A, double lx, double Q2, double* f)
-{
-  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  const double MN2=MN*MN;   // Squared Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  double f2=0., xf3=0., fL=0.;
-  if (lx>0.5) G4cerr<<"***getFun: ln(x)="<<lx<<">.5"<<G4endl;
-  double x=std::exp(lx);
-  double x2=x*x;
-  double c=0.;
-  if(Q2>.0001) c=x2*MN2/Q2;
-  double c2=1+c+c;     // Q^2/nu^2 correction
-  strucf(A, x, Q2, f2, xf3, fL);
-  f[0]=f2;             // direct part
-  f[1]=(-f2+xf3)/x;    // *y (neutrino) part
-  f[2]=(-f2-xf3)/x;    // *y (anti-neutrino) part
-  f[3]=(f2*c2-fL-xf3)/x2; // *y2 (neutrino) part
-  f[4]=(f2*c2-fL+xf3)/x2; // *y2 (anti-neutrino) part
-}
-
-int main()
-{
-  const double reps=.001;        // relative accuracy of the total Q2 integral calculation
-  const double xeps=.0001;      // relative accuracy of the total X integral calculation
-  //           =========
-  const double GF=1.16637e-5;   // Fermi constant in GeV^-2
-  const double GF2=GF*GF;       // Squared Fermi constant in GeV^-4
-  const double MW=80.425;       // Mass of W-boson in GeV
-  const double MW2=MW*MW;       // Squared mass of W-boson in GeV^2
-  const double MW4=MW2*MW2;     // Quadro mass of W-boson in GeV^4
-  const double hc2=38937932300.;// (hc)^2 in GeV^2*10^-38cm2 to convert GeV^-2 to 10^-38cm2
-  const double pif=3.14159265*4;// 4pi
-  const double sik=GF2*hc2/pif; // precalculated coefficient
-  //const double mpi=.1349766;    // pi0 meson mass in GeV
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
-  const double me=.00051099892; // electron mass in GeV
-  const double me2=me*me;       // m_e^2 in GeV^2
-  const double hme2=me2/2;      // .5*m_e^2 in GeV^2
-  //const double mmu=.105658369;  // mu meson mass in GeV
-  //const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
-  //const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
-  //const double mtau=1.777;      // tau meson mass in GeV
-  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
-  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
-  //const double mp=.93827203;    // proton mass in GeV
-  //const double mn=.93956536;    // neutron mass in GeV
-  //const double md=1.87561282;   // deuteron mass in GeV
-  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
-  //const double MD=1.232;        // proton mass in GeV
-  //const double mp2=mp*mp;       // m_p^2 in GeV^2
-  const double MN2=MN*MN;       // M_N^2 in GeV^2
-  const double dMN=MN+MN;       // 2*M_N in GeV
-  const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
-  const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
-  const double EminE=me+me2/dMN;// Threshold for muon production
-  //const double EminMu=mmu+mmu2/dMN;  // Threshold for muon production
-  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
-  //
-  const double mc=.261;           // parameter of W2>(M_N^2+M_D^2)/2 cut for QuasiEl/Delta
-  //const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
-  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
-  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
-  // _____ Begin of Test Area
-  //Genfun::LogGamma logGamma;
-  //double n=4.9;
-		//double g=exp(logGamma(n));
-  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
-  // ^^^^^ End of Test Area
-  //
-  double f[5];                    // A working array
-  int    A=12;                    // Neucleus for which calculations should be done
-  double lEnuMin=0;               // LogLog of Minimum energy of neutrino
-  double lEnuMax=std::log(1.+std::log(300./EminE)); // LogLog of MaximumEnergy of neutrino
-  //double lEnuMax=std::log(1.+std::log(300./EminMu));// LogLog of MaximumEnergy of neutrino
-  int    nE=63;
-  double dlE=(lEnuMax-lEnuMin)/nE;
-  lEnuMin+=dlE/10;
-  lEnuMax+=dlE/5;
-  G4cout<<"Emin="<<EminE<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
-  //G4cout<<"Emin="<<EminMu<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
-  for(double lEnu=lEnuMin; lEnu<lEnuMax; lEnu+=dlE)
-		{
-    double Enu=std::exp(std::exp(lEnu)-1.)*EminE; // Energy of neutrino/anti-neutrino
-    //double Enu=std::exp(std::exp(lEnu)-1.)*EminMu; // Energy of neutrino/anti-neutrino
-    double dEnu=Enu+Enu;          // doubled energy of nu/anu
-    double Enu2=Enu*Enu;          // squared energy of nu/anu
-    double Ee=Enu-me;           // Free Energy of neutrino/anti-neutrino
-    double Ee2=Ee*Ee;          // squared energy of nu/anu
-    //double Emu=Enu-mmu;           // Free Energy of neutrino/anti-neutrino
-    //double Emu2=Emu*Emu;          // squared energy of nu/anu
-    double ME=Enu*MN;             // M*E
-    double dME=ME+ME;             // 2*M*E
-    double DIStsig=1.;            // Total curent DIS cross-section to be integrated
-    double DISmsig=1.e20;         // Total remembered DIS cross-section
-    double dEMN=(dEnu+MN)*ME;
-    double MEm=ME-hme2;
-    //double MEm=ME-hmmu2;
-    double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
-    //double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
-    double E2M=MN*Enu2-(Enu+MN)*hme2;
-    //double E2M=MN*Enu2-(Enu+MN)*hmmu2;
-    double ymax=(E2M+sqE)/dEMN;
-    double ymin=(E2M-sqE)/dEMN;
-    double rmin=1.-ymin;
-    double rhm2E=hme2/Enu2;
-    //double rhm2E=hmmu2/Enu2;
-    double Q2min=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
-    double Q2max=dME*ymax;
-    int    nQ2=8;
-    //G4cout<<"*** E="<<Enu<<", Q2i="<<Q2min<<" < Q2a="<<Q2max<<", yi="<<ymin<<" < ya="
-    //      <<ymax<<G4endl;
-    while(std::fabs(DIStsig-DISmsig)/DIStsig>reps)
-    {
-      DISmsig=DIStsig;
-      DIStsig=0.;
-      nQ2*=2;
-      double dQ2=(Q2max-Q2min)/nQ2;
-      for(double Q2=Q2min+dQ2/2; Q2<Q2max; Q2+=dQ2)
-				  {
-        double DISxint=1.;        // Curent DIS x-integral
-        double DISmint=1.e20;     // Remembered DIS x-integral
-        double Q2M=Q2+MW2;
-        double dik=MW4/Q2M/Q2M;
-        double qmc=Q2+mcV;
-        double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2); // QuasielastBoundary
-        //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c)
-        //double xN=Q2/dME;
-        double xN=Q2/MN/(Ee+std::sqrt(Ee2+Q2));
-        //double xN=Q2/MN/(Emu+std::sqrt(Emu2+Q2));
-        //double lXmin=log(xN/ymax);
-        double lXmin=std::log(xN);
-        // ****** QE ******** tot/qe
-        //if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES
-        // *** End of QE^^^^^
-        double lXmax=0.;               // QES is in DIS
-        //double lXmax=lXQES;          // Cut off quasielastic
-        int    nX=8;
-        while(std::fabs(DISxint-DISmint)/DISxint>xeps)
-        {
-          DISmint=DISxint;
-          DISxint=0.;
-          nX*=2;
-          double dlX=(lXmax-lXmin)/nX;
-          for(double lX=lXmin+dlX/2; lX<lXmax; lX+=dlX)
-				      {
-            getFun(A, lX, Q2, f);
-            DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
-            //DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
-            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
-          }
-          DISxint*=dlX;
-          //G4cout<<"--- E="<<Enu<<" --- Q2="<<Q2<<" --- nX="<<nX<<", iX="<<DISxint
-          //      <<", mX="<<DISmint<<", rX="<<(DISxint-DISmint)/DISxint<<G4endl;
-        }
-        //G4cout<<"(E="<<Enu<<"), Q2="<<Q2<<", I="<<DISxint/dik/dik<<G4endl;
-        DIStsig+=DISxint*dik;
-      }
-      DIStsig*=dQ2;
-      //G4cout<<"=== E="<<Enu<<" ===> nQ="<<nQ2<<", iQ="<<DIStsig<<", mQ="<<DISmsig
-						//      <<", rQ="<<(DIStsig-DISmsig)/DIStsig<<G4endl;
-    }
-    //===== tot/qe choice ====
-    DIStsig*=sik/Enu;
-    G4cout<<"***total*** E="<<Enu<<",sig/E= "<<DIStsig<<G4endl;
-    //...................
-    //DIStsig*=sik;
-    //G4cout<<"***qelas*** E="<<Enu<<",sig= "<<DIStsig<<G4endl;
-    //===== End of the choice
-		} // End of the Enery LOOP
-  // int np=0;
-  //for(int m=0; m<2; m++)
-  //{
-  //  //fileNuMuX<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
-  //  //G4cout<<"**** A_high="<<m<<G4endl;
-  //  np=0;
-  //  int nC=14;
-	 //  for(G4int en=0; en<nC; en++)
-	 //  {
-  //    //G4double sig=1.;
-  //    np++;
-  //    //if(np==7)           // Write by 7 number in brackets
-  //    //{
-  //    //  if(en==nC-1) fileNuMuX<<sig<<"};"<<G4endl;
-  //    //  else         fileNuMuX<<sig<<","<<G4endl<<"    ";
-	 //    //}
-  //    //else           fileNuMuX<<sig<<",";
-  //    //if(np==7) np=0;
-  //  } // End of the point LOOP
-  //} // End of the isotop LOOP
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuNQEQ2Integration.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuNQEQ2Integration.cc	(revision 968)
+++ 	(revision )
@@ -1,183 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS(Q2) fixed step integration
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-2005
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-//#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-// All calculations have been done for C12 nucleus
-
-double nuE(double E) { return .9673/(1.+.323/E/E)/std::pow(E,.78);} // (E is in GeV)
-
-double nuX(double E, double r, double p) // (E is in GeV, r=Q2/Q2max, p=1.+nuE(E))
-{
-  double y=p-r;
-  double p3=(.3088+.0012352*E)/(1.+1.836/E/E);
-  return std::pow(y,6)*(r+p3)/(p3*r+1.);
-}
-
-double anuE(double E) { return .875/(1.+.2977/E/E)/std::pow(E,.78);} // (E is in GeV)
-
-double anuX(double E, double r, double p) // (E is in GeV, r=Q2/Q2max, p=1+anuE(E))
-{
-		double	E2=E*E;
-  double y=p-r;
-  double p3=(13.88+.9373*(1.+.000033*E2)*std::sqrt(E))/(1.+(10.12+1.532/E2)/E);
-  return std::pow(y,6)*(r+p3)/(p3*r+1.);
-}
-
-int main()
-{
-  const double eps=.00000001;    // relative accuracy of the integral calculation
-  const double mmu=.105658369;   // mu meson mass in GeV
-  const double mmu2=mmu*mmu;     // m_mu^2 in GeV^2
-  const double MN=.931494043;    // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  const double Emin=mmu+mmu2/(MN+MN); // the threshold energy in GeV
-  const double Emax=390.;        // the maximum energy in GeV
-  const double lEmi=std::log(Emin);   // logarithm of the threshold energy in GeV
-  const double lEma=std::log(Emax);   // logarithm of the maximum energy in GeV
-  const int    power=7;          // power for the magic variable (E-dependent)
-  const double pconv=1./power;   // conversion power for the magic variable
-  //           =========
-  const int niX=20;              // number of points for X table
-  const int liX=niX-1;           // the last index for X table
-  const int niE=20;              // number of points for E table
-  double Xl[niE][niX];           // reversed table
-  double inl[niE][niX];          // direct table
-  // ----------------- nu/anu switch ------------------------------
-  //bool nu=true;
-  bool nu=false;
-  // --------------------------------------------------------------
-  double dE=(lEma-lEmi)/niE;
-  double hE=dE/2;
-  int ne=0;
-  for(double len=lEmi+hE; len<lEma; len+=dE)
-  {
-    //G4cout<<"log(E)="<<len<<G4endl;
-    double en=std::exp(len);
-    double shift=0.;
-    if(nu) shift=1.+nuE(en);
-    else   shift=1.+anuE(en);
-    double Xma=std::pow(shift,power);
-    double Xmi=std::pow((shift-1.),power);
-    int    nX=8;
-    double DISmsig=0.;
-    double DIStsig=1.;
-    while(std::fabs(DIStsig-DISmsig)/DIStsig>eps)
-    {
-      DISmsig=DIStsig;
-      DIStsig=0.;
-      nX*=2;
-      double dX=(Xma-Xmi)/nX;
-      double hX=dX/2;
-      for(double X=Xmi+hX; X<Xma; X+=dX)
-				  {
-        double r=shift-std::pow(X,pconv); // the same for nu and anu
-        if(nu) DIStsig+=nuX(en,r,shift);  // neutrino
-        else   DIStsig+=anuX(en,r,shift); // anti-neutrino
-      }
-      DIStsig*=dX;
-      //G4cout<<"E="<<en<<",nX="<<nX<<",i="<<DIStsig<<",m="<<DISmsig<<",r="
-      //      <<(DIStsig-DISmsig)/DIStsig<<G4endl;
-		  }
-    //G4cout<<"E="<<en<<", Total: nX="<<nX<<", Int="<<DIStsig<<G4endl;
-    // ***************** Calculate the reversed table *****************
-    double dInt=DIStsig/niX;             // Step for the integral
-    DIStsig=0.;
-    DISmsig=0.;
-    nX*=2;
-    double dX=(Xma-Xmi)/nX;      // New step for the multiplied number of nodes
-    dInt/=dX;                    // To avoid multiplication (nX is fixed)
-    double sum=dInt;
-    int nn=0;
-    double hX=dX/2;
-    G4cout<<"E="<<en<<", Xl_min="<<Xmi<<G4endl;
-    for(double X=Xmi+hX; X<Xma; X+=dX)
-		  {
-      double r=shift-std::pow(X,pconv); // the same for nu and anu
-      if(nu) DIStsig+=nuX(en,r,shift);  // neutrino
-      else   DIStsig+=anuX(en,r,shift); // anti-neutrino
-      if(DIStsig>sum+eps)
-		  		{
-        inl[ne][nn]=sum*dX;
-        Xl[ne][nn]=X-(DIStsig-sum)*dX/(DIStsig-DISmsig);
-        G4cout<<"E="<<en<<", sum="<<inl[ne][nn]<<", Xl["<<nn<<"]="<<Xl[ne][nn]<<G4endl;
-        nn++;
-        sum+=dInt;
-      }
-      DISmsig=DIStsig;
-    }
-    inl[ne][nn]=sum*dX;
-    Xl[ne][nn]=Xma;
-    G4cout<<"E="<<en<<", sum="<<inl[ne][nn]<<", Xl["<<nn<<"]="<<Xma<<G4endl;
-    // ************* The following is just a test (better than .5%) ************
-    //for(int i=1; i<niX; i++)
-    //{
-    //  DIStsig=0.;
-    //  double Xm=Xl[ne][i];
-    //  double dX=(Xm-Xmi)/nX;
-    //  double hX=dX/2;
-    //  for(double X=Xmi+hX; X<Xm; X+=dX)
-		  //		{
-    //    double r=shift-pow(X,pconv); // the same for nu and anu
-    //    if(nu) DIStsig+=nuX(en,r,shift);  // neutrino
-    //    else   DIStsig+=anuX(en,r,shift); // anti-neutrino
-    //  }
-    //  DIStsig*=dX;
-    //  G4cout<<"i="<<i<<", v="<<DIStsig<<", d="<<fabs(DIStsig-inl[ne][i])/DIStsig<<G4endl;
-    //}
-    // ***************** Calculate the direct table *****************
-    dX=(Xma-Xmi)/niX;
-    double nor=inl[ne][liX]/niX;
-    //G4cout<<"E="<<en<<", i=0, Xm="<<Xmi<<", I=0."<<G4endl;
-    for(int i=1; i<=niX; i++)
-    {
-      DIStsig=0.;
-      double Xm=Xmi+dX*i;
-      double rX=(Xm-Xmi)/nX;
-      double hX=rX/2;
-      for(double X=Xmi+hX; X<Xm; X+=rX)
-		  		{
-        double r=shift-std::pow(X,pconv); // the same for nu and anu
-        if(nu) DIStsig+=nuX(en,r,shift);  // neutrino
-        else   DIStsig+=anuX(en,r,shift); // anti-neutrino
-      }
-      DIStsig*=rX/nor;
-      //////////G4cout<<"E="<<en<<", i="<<i<<", Xm="<<Xm<<", I="<<DIStsig<<G4endl;
-    }
-    ne++;
-  } // End of the big loop over log(E)
-  G4cout<<"End"<<G4endl;
-  return 0;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuQEQ2Integration.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuQEQ2Integration.cc	(revision 968)
+++ 	(revision )
@@ -1,305 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS(Q2) fixed step integration
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-2005
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-//#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-// All calculations have been done for C12 nucleus
-
-double nuQ2(double Q2) // (Q2 is in GeV^2)
-{
-  //  x      *     *     *     *     *     *     *     *      *     *
-  // -1-    -2-   -3-   -4-   -5-   -6-   -7-   -8-   -9-   -10-  -11-
-  // 2.068 2.634 4.333 3.000 10.59 .0011 18.74 1.484 134950 .0755 4.5
-  G4double y=Q2/3;
-  G4double z=std::pow((Q2/2.634),4.333);
-		// The first (1.+Q2)**4.5 fuctor is just normalization, which payed back by 1/(1+Q2)^3.5
-  G4double f=std::pow((1.+Q2),4.5)*
-				(1.+z*(1.+18.74*std::exp(-Q2/1.484)-134950*std::exp(-Q2/.0755)))/std::pow((1.+y+.0011*y*y),10.59);
-  //G4cout<<"Q2="<<Q2<<", y="<<y<<", z="<<z<<", f="<<f<<G4endl;
-  return f;
-}
-
-double anuQ2(double Q2)
-{
-  //    x     *     *     *     *     *     *     *     -     *
-  //   -1-   -2-   -3-   -4-   -5-   -6-   -7-   -8-   -9-   -10-
-  // .2000 .0100 .3000 .4000 3.900 .0150 .2000 6.500 0.000 .00000001
-  G4double  y=Q2/.4;
-  G4double  y2=y*y;
-  G4double  z=std::pow((Q2/.01),-.3);
-		// The first (1.+Q2)**6.5 fuctor is just normalization, which payed back by 1/(1+Q2)^5.5
-  G4double f=(1.+z)*std::pow((1+Q2),6.5)*std::pow(Q2,-.2)/std::pow((1.+y+(.015+.00000001*y2)*y2),3.9);
-  //G4cout<<"Q2="<<Q2<<", y="<<y<<", z="<<z<<", f="<<f<<G4endl;
-  return f;
-}
-
-int main()
-{
-  const double eps=.000001;      // relative accuracy of the integral calculation
-  //           =========
-  const double mmu=.105658369;  // mu meson mass in GeV
-  const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
-  const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
-  //const double mtau=1.777;      // tau meson mass in GeV
-  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
-  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
-  //const double mp=.93827203;    // proton mass in GeV
-  //const double mn=.93956536;    // neutron mass in GeV
-  //const double md=1.87561282;   // deuteron mass in GeV
-  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
-  //const double MD=1.232;        // proton mass in GeV
-  //const double mp2=mp*mp;       // m_p^2 in GeV^2
-  const double MN2=MN*MN;       // M_N^2 in GeV^2
-  const double dMN=MN+MN;       // 2*M_N in GeV
-  //const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
-  //const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
-  //const double EminE=me+me2/dMN;// Threshold for muon production
-  const double Emin=mmu+mmu2/dMN;  // Threshold for muon production
-  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
-  //
-  //const double mc=.3;           // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  //const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  //const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
-  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
-  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
-  // _____ Begin of Test Area
-  //Genfun::LogGamma logGamma;
-  //double n=4.9;
-		//double g=exp(logGamma(n));
-  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
-  // ^^^^^ End of Test Area
-  //
-  const int niQ2=100;           // Number of points in the Q2 integration
-  const int diQ2=niQ2+1;        // Dimention for arrays
-  double Xl[diQ2];
-  double inl[diQ2];
-  double Enu=Emin;              // Initial E=minE  
-  double Emo=Emin+Emin;         // First new
-  double Emv=0.;
-  double Q2i=1.;
-  double Q2o=1.;
-  double Q2v=0.;
-		int nit=0;
-  while (std::fabs(Q2v-Q2o)>eps && nit<10)
-  {
-    nit++;
-    double Emm=Emo;               // Candidate for the next point is predefined
-    double dEnu=Enu+Enu;          // doubled energy of nu/anu
-    double Enu2=Enu*Enu;          // squared energy of nu/anu
-    double ME=Enu*MN;             // M*E
-    double dEMN=(dEnu+MN)*ME;
-    double MEm=ME-hmmu2;
-    double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
-    double E2M=MN*Enu2-(Enu+MN)*hmmu2;
-    //double ymax=(E2M+sqE)/dEMN;
-    double ymin=(E2M-sqE)/dEMN;
-    double rmin=1.-ymin;
-    double rhm2E=hmmu2/Enu2;
-    Q2i=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // minQ2 for Enu
-    //G4cout<<"MinSearch: E="<<Enu<<": Q2="<<Q2i<<G4endl;
-    if(Emv>0.0001) // Not initialization (not first two) Use three points
-				{
-      //double d=Enu*Emo*(Enu-Emo)+Enu*Emv*(Emv-Enu)+Emo*Emv*(Emo-Emv);
-      double a=Q2i*Emo-Enu*Q2o+Q2v*Enu-Emv*Q2i+Q2o*Emv-Emo*Q2v;
-						double b=Enu*Enu*(Q2o-Q2v)+Emo*Emo*(Q2v-Q2i)+Emv*Emv*(Q2i-Q2o);
-      //double a
-				  Emm=-b/(a+a);
-      if(Q2v<Q2i)                 // swap to make Q2v the biggest
-				  {
-        double q2=Q2v;
-        double en=Emv;
-        Q2v=Q2i;
-        Emv=Enu;
-        Q2i=q2;
-        Enu=en;
-      }
-      if(Q2v<Q2o)                 // swap to make Q2v the biggest
-				  {
-        double q2=Q2v;
-        double en=Emv;
-        Q2v=Q2o;
-        Emv=Emo;
-        Q2o=q2;
-        Emo=en;
-      }
-    }
-    else if(Q2o>0.0000000001) Emm=Enu+Emin; // the second step
-    //G4cout<<"***Enu="<<Enu<<", Emm="<<Emm<<", Emin="<<Emin<<G4endl;
-    Emv=Emo;
-    Q2v=Q2o;
-    Emo=Enu;
-    Q2o=Q2i;
-    Enu=Emm;
-    //G4cout<<"___Q2i="<<Q2i<<": Q2o="<<Q2o<<", Q2v="<<Q2v<<G4endl;
-  }
-  double Q2min=Q2i-eps;
-  double Q2max=600.;            // covers the calculated region
-  // ----------------- nu/anu switch ------------------------------
-  //bool nu=true;
-  bool nu=false;
-  // -----------------------------**************************-------
-  // *************** Convert to the magic variable ****************
-  double Xmin=0.;
-  double Xmax=0.;
-		if(nu)
-  {
-    Xmin=std::pow(1+Q2max,-3.5);
-    Xmax=std::pow(1+Q2min,-3.5);
-		}  
-  else
-  {
-    Xmin=std::pow(1+Q2max,-5.5);
-    Xmax=std::pow(1+Q2min,-5.5);
-		}
-  G4cout<<"Ei="<<Enu<<":Q2i="<<Q2min<<",Q2a="<<Q2max<<",Xmi="<<Xmin<<",Xma="<<Xmax<<G4endl;
-  int    nQ2=8;                 // nitial #Of points for the overall integration
-  double DISmsig=0.;
-  double DIStsig=1.;
-		double Q2=0.;                 // Prototype for the integration 
-  while(std::fabs(DIStsig-DISmsig)/DIStsig>eps)
-  {
-    DISmsig=DIStsig;
-    DIStsig=0.;
-    nQ2*=2;
-    double dX=(Xmax-Xmin)/nQ2;
-    double hX=dX/2;
-    for(double X=Xmin+hX; X<Xmax; X+=dX)
-				{
-      if(nu) // neutrino
-						{
-        Q2=std::pow(X,-.2857142); // -1./3.5
-        DIStsig+=nuQ2(Q2);
-      }
-      else  // anti-neutrino
-						{
-        Q2=std::pow(X,-.1818182); // -1./5.5
-        DIStsig+=anuQ2(Q2);
-      }
-    }
-    DIStsig*=dX;
-    //G4cout<<"n="<<nQ2<<",i="<<DIStsig<<",m="<<DISmsig<<",r="<<(DIStsig-DISmsig)/DIStsig
-				//		    <<G4endl;
-		}
-  G4cout<<"Total: nQ2="<<nQ2<<", Int="<<DIStsig<<G4endl;
-  // ***************** Calculate the reversed table *****************
-  double dInt=DIStsig/niQ2;             // Step for the integral
-  DIStsig=0.;
-  DISmsig=0.;
-  nQ2*=2;
-  double dX=(Xmax-Xmin)/nQ2;
-  dInt/=dX;                     // To avoid multiplication (nQ2 is fixed)
-  double sum=dInt;
-  int nn=1;
-  Xl[0]=0.;
-  double hX=dX/2;
-  for(double X=Xmin+hX; X<Xmax; X+=dX)
-		{
-    if(nu) // neutrino
-				{
-      Q2=std::pow(X,-.2857142); // -1./3.5
-      DIStsig+=nuQ2(Q2);
-    }
-    else  // anti-neutrino
-				{
-      Q2=std::pow(X,-.1818182); // -1./5.5
-      DIStsig+=anuQ2(Q2);
-    }
-    if(DIStsig>sum+eps)
-				{
-      inl[nn]=sum*dX;
-      Xl[nn]=X-(DIStsig-sum)*dX/(DIStsig-DISmsig);
-      G4cout<<"sum="<<inl[nn]<<", Xl["<<nn<<"]="<<Xl[nn]<<G4endl;
-      nn++;
-      sum+=dInt;
-    }
-    DISmsig=DIStsig;
-  }
-  inl[nn]=sum*dX;
-  Xl[nn]=Xmax;
-  G4cout<<"sum="<<inl[nn]<<", Xl["<<nn<<"]="<<Xmax<<G4endl;
-  // ************* The following is just a test (better than .5%) ************
-  for(int i=1; i<=niQ2; i++)
-  {
-    DIStsig=0.;
-    double Xm=Xl[i];
-    double dX=(Xm-Xmin)/nQ2;
-    double hX=dX/2;
-    for(double X=Xmin+hX; X<Xm; X+=dX)
-				{
-      if(nu) // neutrino
-						{
-        Q2=std::pow(X,-.2857142); // -1./3.5
-        DIStsig+=nuQ2(Q2);
-      }
-      else  // anti-neutrino
-						{
-        Q2=std::pow(X,-.1818182); // -1./5.5
-        DIStsig+=anuQ2(Q2);
-      }
-    }
-    DIStsig*=dX;
-    G4cout<<"i="<<i<<", v="<<DIStsig<<", d="<<std::fabs(DIStsig-inl[i])/DIStsig<<G4endl;
-  }
-  // ***************** Calculate the direct table *****************
-  dX=(Xmax-Xmin)/niQ2;
-  double nor=inl[niQ2]/niQ2;
-  for(int i=1; i<=niQ2; i++)
-  {
-    DIStsig=0.;
-    double Xm=Xmin+dX*i;
-    double rX=(Xm-Xmin)/nQ2;
-    double hX=rX/2;
-    for(double X=Xmin+hX; X<Xm; X+=rX)
-				{
-      if(nu) // neutrino
-						{
-        Q2=std::pow(X,-.2857142); // -1./3.5
-        DIStsig+=nuQ2(Q2);
-      }
-      else  // anti-neutrino
-						{
-        Q2=std::pow(X,-.1818182); // -1./5.5
-        DIStsig+=anuQ2(Q2);
-      }
-    }
-    DIStsig*=rX/nor;
-    G4cout<<"i="<<i<<", Xm="<<Xm<<", I="<<DIStsig<<G4endl;
-  }
-  G4cout<<"End"<<G4endl;
-  DIStsig=0.;
-  return 0;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXIntegration_Q2_QE.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXIntegration_Q2_QE.cc	(revision 968)
+++ 	(revision )
@@ -1,312 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS (x,Q2) approximation is integrated over x
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 30-Sept-05
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-void strucf(int A, double x, double Q2, double& f2, double& xf3, double& fL)
-{
-  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN2=MN*MN;       // M_N^2 in GeV^2
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double Wt=MN+mpi;       // Delta threshold
-  //const double W2t=Wt*Wt;       // Squared Delta threshold
-  const Genfun::LogGamma lGam;
-  static int mA=0;
-  static double mQ2=0., mN, mD, mDel, mU2, mU3, mV;
-  //static double mUU;
-  double N=3., D=0., Del=0., U2=0., U3=0., V=0.;
-  //double UU=0.;
-  if(A==mA && Q2==mQ2)         // Associative memory for acceleration
-		{
-    N  =mN;
-    D  =mD;
-    Del=mDel;
-    U2 =mU2;
-    U3 =mU3;
-    //UU =mUU;
-    V  =mV;
-  }
-  else
-  {
-		  double r=0.;
-    double max=1.;
-    double H=1.22;
-    if(A==1)                   // Proton
-		  {
-      r=std::sqrt(Q2/1.66);
-      max=.5;
-    }
-    else if(A<13)              // Light nuclei
-		  {
-      double f=Q2/4.62;
-      r=f*f;
-      max=.3;
-      if(A>2) H=1.;
-    }
-    else if(A>0)               // Heavy nuclei
-		  {
-      double f=Q2/3.4;
-      double ff=f*f;
-      r=ff*ff;
-      max=.5;
-      H=1.;
-    }
-    else G4cout<<"strucf: A="<<A<<" <= 0"<< G4endl;
-    //
-    N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
-    Del=(1.+r)/(12.5+r/max);
-    double S=std::pow(1.+.6/Q2,-1.-Del);
-    D=H*S*(1.-.5*S);
-    V=3*(1.-D)*(N-1.);
-    double cc=Q2/.08;
-    double cc2=cc*cc;
-    //double C=(1.+cc2)/(1.+cc2/.24);              // Weak?
-    double C=(1.+cc2)/(1.+cc2/.24)/(1.+Q2/21.6); // EM
-    double c3=C+C+C;
-    double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
-    U2=(c3+N-3.)*uu;
-    U3=c3*uu;
-    //UU=uu+uu+uu; // @@
-    mA  = A;
-    mQ2 = Q2;
-    mN  = N;
-    mD  = D;
-    mDel=Del;
-    mU2 =U2;
-    mU3 =U3;
-    //mUU =UU; // @@
-    mV  =V;
-  }  
-  // From here the Q2 coefficients are used
-  double x1=std::pow(1.-x,N-2.);
-  double pp=D*std::pow(x,-Del)*x1;
-  double dir=V*x*x1;
-  double per=U2*pp;
-		f2 = per + dir;
-  //double W2=MN2-MN2*x+Q2/x-Q2;
-  //if(W2<W2t)
-		//{
-  //  per=UU*pp;
-  //  xf3= per+dir;
-  //}
-  //else
-  xf3= U3*pp+dir;
-  fL = per/4.;
-  return;
-}
-
-void getFun(int A, double lx, double Q2, double* f)
-{
-  double f2=0., xf3=0., fL=0.;
-  if (lx>0.5) G4cerr<<"***getFun: ln(x)="<<lx<<">.5"<<G4endl;
-  double x=std::exp(lx);
-  double x2=x*x;
-  strucf(A, x, Q2, f2, xf3, fL);
-  f[0]=f2;             // direct part
-  f[1]=(-f2+xf3)/x;    // *y (neutrino) part
-  f[2]=(-f2-xf3)/x;    // *y (anti-neutrino) part
-  f[3]=(f2-fL-xf3)/x2; // *y2 (neutrino) part
-  f[4]=(f2-fL+xf3)/x2; // *y2 (anti-neutrino) part
-}
-
-int main()
-{
-  const double reps=.001;        // relative accuracy of the total Q2 integral calculation
-  const double xeps=.0001;      // relative accuracy of the total X integral calculation
-  //           =========
-  const double GF=1.16637e-5;   // Fermi constant in GeV^-2
-  const double GF2=GF*GF;       // Squared Fermi constant in GeV^-4
-  const double MW=80.425;       // Mass of W-boson in GeV
-  const double MW2=MW*MW;       // Squared mass of W-boson in GeV^2
-  const double MW4=MW2*MW2;     // Quadro mass of W-boson in GeV^4
-  const double hc2=38937932300.;// (hc)^2 in GeV^2*10^-38cm2 to convert GeV^-2 to 10^-38cm2
-  const double pif=3.14159265*4;// 4pi
-  const double sik=GF2*hc2/pif; // precalculated coefficient
-  //const double mpi=.1349766;    // pi0 meson mass in GeV
-  const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
-  //const double me=.00051099892; // electron mass in GeV
-  //const double me2=me*me;       // m_e^2 in GeV^2
-  //const double hme2=me2/2;      // .5*m_e^2 in GeV^2
-  const double mmu=.105658369;  // mu meson mass in GeV
-  const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
-  const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
-  //const double mtau=1.777;      // tau meson mass in GeV
-  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
-  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
-  //const double mp=.93827203;    // proton mass in GeV
-  //const double mn=.93956536;    // neutron mass in GeV
-  //const double md=1.87561282;   // deuteron mass in GeV
-  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
-  //const double MD=1.232;        // proton mass in GeV
-  //const double mp2=mp*mp;       // m_p^2 in GeV^2
-  const double MN2=MN*MN;       // M_N^2 in GeV^2
-  const double dMN=MN+MN;       // 2*M_N in GeV
-  const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
-  const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
-  //const double EminE=me+me2/dMN;// Threshold for muon production
-  const double EminMu=mmu+mmu2/dMN;  // Threshold for muon production
-  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
-  //
-  //const double mc=.3;           // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
-  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
-  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
-  // _____ Begin of Test Area
-  //Genfun::LogGamma logGamma;
-  //double n=4.9;
-		//double g=exp(logGamma(n));
-  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
-  // ^^^^^ End of Test Area
-  //
-  double f[5];                    // A working array
-  int    A=12;                    // Neucleus for which calculations should be done
-  double lEnuMin=std::log(EminMu);     // Log of Minimum energy of neutrino
-  double lEnuMax=std::log(300.);       // Log of Maximum energy of neutrino
-  int    nE=20;
-  double dlE=(lEnuMax-lEnuMin)/nE;
-  lEnuMin+=dlE/10;
-  lEnuMax+=dlE/5;
-  G4cout<<"Emin="<<EminMu<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
-  for(double lEnu=lEnuMin; lEnu<lEnuMax; lEnu+=dlE)
-		{
-    double Enu=std::exp(lEnu);         // Energy of neutrino/anti-neutrino
-    double dEnu=Enu+Enu;          // doubled energy of nu/anu
-    double Enu2=Enu*Enu;          // squared energy of nu/anu
-    double Emu=Enu-mmu;           // Free Energy of neutrino/anti-neutrino
-    double Emu2=Emu*Emu;          // squared energy of nu/anu
-    double ME=Enu*MN;             // M*E
-    double dME=ME+ME;             // 2*M*E
-    double DIStsig=1.;            // Total curent DIS cross-section to be integrated
-    double DISmsig=1.e20;         // Total remembered DIS cross-section
-    double dEMN=(dEnu+MN)*ME;
-    double MEm=ME-hmmu2;
-    double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
-    double E2M=MN*Enu2-(Enu+MN)*hmmu2;
-    double ymax=(E2M+sqE)/dEMN;
-    double ymin=(E2M-sqE)/dEMN;
-    double rmin=1.-ymin;
-    double rhm2E=hmmu2/Enu2;
-    double Q2min=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
-    double Q2max=dME*ymax;
-    // *** Additional notQE limit ****
-    //double Q2nqe=Emu*dMN-mcV;
-    //if(Q2max>Q2nqe) Q2max=Q2nqe; // !! only for not QE
-    // *** End of notQE ^^^^^
-    int    nQ2=8;
-    //G4cout<<"*** E="<<Enu<<", Q2i="<<Q2min<<" < Q2a="<<Q2max<<", yi="<<ymin<<" < ya="
-    //      <<ymax<<G4endl;
-    //while(fabs(DIStsig-DISmsig)/DIStsig>reps)
-    //{
-      DISmsig=DIStsig;
-      DIStsig=0.;
-      nQ2*=2;
-      double dQ2=(Q2max-Q2min)/nQ2;
-      if(dQ2>0) for(double Q2=Q2min+dQ2/2; Q2<Q2max; Q2+=dQ2)
-				  {
-        double DISxint=1.;        // Curent DIS x-integral
-        double DISmint=1.e20;     // Remembered DIS x-integral
-        double Q2M=Q2+MW2;
-        double dik=MW4/Q2M/Q2M;
-        double qmc=Q2+mcV;
-        double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2); // Quasielastic boundary
-        //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c)
-        //double xN=Q2/dME;
-        double xN=Q2/MN/(Emu+std::sqrt(Emu2+Q2));
-        //double lXmin=log(xN/ymax);
-        double lXmin=std::log(xN);
-        double lXmax=0.;               // QES is in DIS
-        // ****** QE ********
-        if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES
-        //if(lXQES<lXmax) lXmax=lXQES;   // A cut which excludes QES (see above "notQE")
-        // *** End of QE^^^^^
-        //double lXmax=lXQES;          // Cut off quasielastic
-        int    nX=8;
-        while(std::fabs(DISxint-DISmint)/DISxint>xeps)
-        {
-          DISmint=DISxint;
-          DISxint=0.;
-          nX*=2;
-          double dlX=(lXmax-lXmin)/nX;
-          for(double lX=lXmin+dlX/2; lX<lXmax; lX+=dlX)
-				      {
-            getFun(A, lX, Q2, f);
-            DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
-            //DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
-            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
-          }
-          DISxint*=dlX;
-          //G4cout<<"--- E="<<Enu<<" --- Q2="<<Q2<<" --- nX="<<nX<<", iX="<<DISxint
-          //      <<", mX="<<DISmint<<", rX="<<(DISxint-DISmint)/DISxint<<G4endl;
-        }
-        //G4cout<<"(E="<<Enu<<"), r="<<Q2/Q2mx<<" "<<DISxint*dik*sik<<G4endl;//NoNormedByEn
-        // QE
-        G4cout<<"(E="<<Enu<<"), Q2="<<Q2<<" "<<DISxint*dik*sik<<G4endl;//Not Normed By En
-        // notQE
-        //G4cout<<"(E="<<Enu<<"), r="<<Q2/Q2max<<" "<<DISxint*dik*sik/Enu<<G4endl;//NormByE
-        DIStsig+=DISxint*dik;
-      }
-      DIStsig*=dQ2;
-      //G4cout<<"=== E="<<Enu<<" ===> nQ="<<nQ2<<", iQ="<<DIStsig<<", mQ="<<DISmsig
-						//      <<", rQ="<<(DIStsig-DISmsig)/DIStsig<<G4endl;
-				//}
-    DIStsig*=sik/Enu;
-    G4cout<<"*** E="<<Enu<<",sig/E="<<DIStsig<<G4endl;
-		} // End of the Enery LOOP
-  // int np=0;
-  //for(int m=0; m<2; m++)
-  //{
-  //  //fileNuMuX<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
-  //  //G4cout<<"**** A_high="<<m<<G4endl;
-  //  np=0;
-  //  int nC=14;
-	 //  for(G4int en=0; en<nC; en++)
-	 //  {
-  //    //G4double sig=1.;
-  //    np++;
-  //    //if(np==7)           // Write by 7 number in brackets
-  //    //{
-  //    //  if(en==nC-1) fileNuMuX<<sig<<"};"<<G4endl;
-  //    //  else         fileNuMuX<<sig<<","<<G4endl<<"    ";
-	 //    //}
-  //    //else           fileNuMuX<<sig<<",";
-  //    //if(np==7) np=0;
-  //  } // End of the point LOOP
-  //} // End of the isotop LOOP
-  return 0;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXQ2Integration.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXQ2Integration.cc	(revision 968)
+++ 	(revision )
@@ -1,310 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS (x,Q2) approximation is integrated over x
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 30-Sept-05
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-void strucf(int A, double x, double Q2, double& f2, double& xf3, double& fL)
-{
-  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN2=MN*MN;       // M_N^2 in GeV^2
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double Wt=MN+mpi;       // Delta threshold
-  //const double W2t=Wt*Wt;       // Squared Delta threshold
-  const Genfun::LogGamma lGam;
-  static int mA=0;
-  static double mQ2=0., mN, mD, mDel, mU2, mU3, mV;
-  //static double mUU;
-  double N=3., D=0., Del=0., U2=0., U3=0., V=0.;
-  //double UU=0.;
-  if(A==mA && Q2==mQ2)         // Associative memory for acceleration
-		{
-    N  =mN;
-    D  =mD;
-    Del=mDel;
-    U2 =mU2;
-    U3 =mU3;
-    //UU =mUU;
-    V  =mV;
-  }
-  else
-  {
-		  //double r=0.;
-    double min=.077;             // Delta(0)
-				double max=.47;              // Hard Pomeron
-    //double H=1.22;
-    if(Q2<=0.) G4cout<<"strucf: Q2="<<Q2<<" <= 0"<< G4endl;
-    double Q=std::sqrt(Q2);
-    //if(A==1)                   // Proton
-		  //{
-    //  r=std::sqrt(Q2/1.66);
-    //  max=.5;
-    //}
-    //else if(A<13)              // Light nuclei
-		  //{
-    //  double f=Q2/4.62;
-    //  r=f*f;
-    //  max=.3;
-    //  if(A>2) H=1.;
-    //}
-    //else if(A>0)               // Heavy nuclei
-		  //{
-    //  double f=Q2/3.4;
-    //  double ff=f*f;
-    //  r=ff*ff;
-    //  max=.5;
-    //  H=1.;
-    //}
-    //else G4cout<<"strucf: A="<<A<<" <= 0"<< G4endl;
-    //
-    N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
-    Del=min+(max-min)/(1.+5./Q);
-    D=0.68*std::pow(1.+.145/Q2,-1.-Del); // 0.68=1-0.34, m2=.145 GeV2
-    V=3*(1.-D)*(N-1.);
-    double c3=.75;                       // 3*0.25
-    double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
-    U2=(c3+N-3.)*uu;
-    U3=c3*uu;
-    mA  = A;
-    mQ2 = Q2;
-    mN  = N;
-    mD  = D;
-    mDel=Del;
-    mU2 =U2;
-    mU3 =U3;
-    //mUU =UU; // @@
-    mV  =V;
-  }  
-  // From here the Q2 coefficients are used
-  double x1=std::pow(1.-x,N-2.);
-  double pp=D*std::pow(x,-Del)*x1;
-  double dir=V*x*x1;
-  double per=U2*pp;
-		f2 = per + dir;
-  //double W2=MN2-MN2*x+Q2/x-Q2;
-  //if(W2<W2t)
-		//{
-  //  per=UU*pp;
-  //  xf3= per+dir;
-  //}
-  //else
-  xf3= U3*pp+dir;
-  fL = per/5.;                             // 20%
-  return;
-}
-
-void getFun(int A, double lx, double Q2, double* f)
-{
-  double f2=0., xf3=0., fL=0.;
-  if (lx>0.5) G4cerr<<"***getFun: ln(x)="<<lx<<">.5"<<G4endl;
-  double x=std::exp(lx);
-  double x2=x*x;
-  strucf(A, x, Q2, f2, xf3, fL);
-  f[0]=f2;             // direct part
-  f[1]=(-f2+xf3)/x;    // *y (neutrino) part
-  f[2]=(-f2-xf3)/x;    // *y (anti-neutrino) part
-  f[3]=(f2-fL-xf3)/x2; // *y2 (neutrino) part
-  f[4]=(f2-fL+xf3)/x2; // *y2 (anti-neutrino) part
-}
-
-int main()
-{
-  const double reps=.001;        // relative accuracy of the total Q2 integral calculation
-  const double xeps=.0001;      // relative accuracy of the total X integral calculation
-  //           =========
-  const double GF=1.16637e-5;   // Fermi constant in GeV^-2
-  const double GF2=GF*GF;       // Squared Fermi constant in GeV^-4
-  const double MW=80.425;       // Mass of W-boson in GeV
-  const double MW2=MW*MW;       // Squared mass of W-boson in GeV^2
-  const double MW4=MW2*MW2;     // Quadro mass of W-boson in GeV^4
-  const double hc2=38937932300.;// (hc)^2 in GeV^2*10^-38cm2 to convert GeV^-2 to 10^-38cm2
-  const double pif=3.14159265*4;// 4pi
-  const double sik=GF2*hc2/pif; // precalculated coefficient
-  //const double mpi=.1349766;    // pi0 meson mass in GeV
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
-  //const double me=.00051099892; // electron mass in GeV
-  //const double me2=me*me;       // m_e^2 in GeV^2
-  //const double hme2=me2/2;      // .5*m_e^2 in GeV^2
-  const double mmu=.105658369;  // mu meson mass in GeV
-  const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
-  const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
-  //const double mtau=1.777;      // tau meson mass in GeV
-  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
-  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
-  const double mp=.93827203;    // proton mass in GeV
-  const double mn=.93956536;    // neutron mass in GeV
-  //const double md=1.87561282;   // deuteron mass in GeV
-  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
-  const double MD=1.232;        // proton mass in GeV
-  //const double mp2=mp*mp;       // m_p^2 in GeV^2
-  const double MN2=MN*MN;       // M_N^2 in GeV^2
-  const double dMN=MN+MN;       // 2*M_N in GeV
-  const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
-  const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
-  //const double EminE=me+me2/dMN;// Threshold for muon production
-  const double EminMu=mmu+mmu2/dMN;  // Threshold for muon production
-  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
-  //
-  //const double mc=.3;           // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  const double mc=std::sqrt((MN*MN+MD*MD)/2)-MN; // Squared mean between N and \Delta
-  //const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
-  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
-  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
-  // _____ Begin of Test Area
-  //Genfun::LogGamma logGamma;
-  //double n=4.9;
-		//double g=exp(logGamma(n));
-  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
-  // ^^^^^ End of Test Area
-  //
-  double f[5];                    // A working array
-  int    A=12;                    // Neucleus for which calculations should be done
-  double lEnuMin=0;               // LogLog of Minimum energy of neutrino
-  double lEnuMax=std::log(1.+std::log(300./EminMu)); // LogLog of MaximumEnergy of neutrino
-  int    nE=49;                   // Number of points
-  double dlE=(lEnuMax-lEnuMin)/nE;
-  lEnuMin+=dlE/10;
-  lEnuMax+=dlE/5;
-  G4cout<<"Emin="<<EminMu<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
-  for(double lEnu=lEnuMin; lEnu<lEnuMax; lEnu+=dlE)
-		{
-    double Enu=std::exp(std::exp(lEnu)-1.)*EminMu; // Energy of neutrino/anti-neutrino
-    double dEnu=Enu+Enu;          // doubled energy of nu/anu
-    double Enu2=Enu*Enu;          // squared energy of nu/anu
-    double Emu=Enu-mmu;           // Free Energy of neutrino/anti-neutrino
-    double Emu2=Emu*Emu;          // squared energy of nu/anu
-    double ME=Enu*MN;             // M*E
-    double dME=ME+ME;             // 2*M*E
-    double DIStsig=1.;            // Total curent DIS cross-section to be integrated
-    double DISmsig=1.e20;         // Total remembered DIS cross-section
-    double dEMN=(dEnu+MN)*ME;
-    double MEm=ME-hmmu2;
-    double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
-    double E2M=MN*Enu2-(Enu+MN)*hmmu2;
-    double ymax=(E2M+sqE)/dEMN;
-    double ymin=(E2M-sqE)/dEMN;
-    double rmin=1.-ymin;
-    double rhm2E=hmmu2/Enu2;
-    double Q2min=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
-    double Q2max=dME*ymax;
-    int    nQ2=8;
-    //G4cout<<"*** E="<<Enu<<", Q2i="<<Q2min<<" < Q2a="<<Q2max<<", yi="<<ymin<<" < ya="
-    //      <<ymax<<G4endl;
-    while(std::fabs(DIStsig-DISmsig)/DIStsig>reps)
-    {
-      DISmsig=DIStsig;
-      DIStsig=0.;
-      nQ2*=2;
-      double dQ2=(Q2max-Q2min)/nQ2;
-      for(double Q2=Q2min+dQ2/2; Q2<Q2max; Q2+=dQ2)
-				  {
-        double DISxint=1.;        // Curent DIS x-integral
-        double DISmint=1.e20;     // Remembered DIS x-integral
-        double Q2M=Q2+MW2;
-        double dik=MW4/Q2M/Q2M;
-        double qmc=Q2+mcV;
-        double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2);//QuasielasticBoundary
-        //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c)
-        //double xN=Q2/dME;
-        double xN=Q2/MN/(Emu+std::sqrt(Emu2+Q2));
-        //double lXmin=log(xN/ymax);
-        double lXmin=std::log(xN);
-        // ****** QE ********
-        if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES >>>>>>>>>>>>>>>>
-        // *** End of QE^^^^^
-        double lXmax=0.;               // QES is in DIS
-        //double lXmax=lXQES;          // Cut off quasielastic
-        int    nX=8;
-        while(std::fabs(DISxint-DISmint)/DISxint>xeps)
-        {
-          DISmint=DISxint;
-          DISxint=0.;
-          nX*=2;
-          double dlX=(lXmax-lXmin)/nX;
-          for(double lX=lXmin+dlX/2; lX<lXmax; lX+=dlX)
-				      {
-            getFun(A, lX, Q2, f);
-            // ***** Neutrino/Antineutrino switch ******>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-            //DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
-            DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
-            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
-          }
-          DISxint*=dlX;
-          //G4cout<<"--- E="<<Enu<<" --- Q2="<<Q2<<" --- nX="<<nX<<", iX="<<DISxint
-          //      <<", mX="<<DISmint<<", rX="<<(DISxint-DISmint)/DISxint<<G4endl;
-        }
-        //G4cout<<"(E="<<Enu<<"), Q2="<<Q2<<", I="<<DISxint/dik/dik<<G4endl;
-        DIStsig+=DISxint*dik;
-      }
-      DIStsig*=dQ2;
-      //G4cout<<"=== E="<<Enu<<" ===> nQ="<<nQ2<<", iQ="<<DIStsig<<", mQ="<<DISmsig
-						//      <<", rQ="<<(DIStsig-DISmsig)/DIStsig<<G4endl;
-    }
-    //===== tot/qe choice ====
-    DIStsig*=sik/Enu;
-    //G4cout<<"***total-neutrino*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
-    G4cout<<"***total-antineutrino*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
-    //G4cout<<"***quasiel-nu*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
-    //G4cout<<"***quasiel-antinu*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
-    //...................
-    //DIStsig*=sik;
-    //G4cout<<"***qelas*** E="<<Enu<<",sig="<<DIStsig<<G4endl;
-    //===== End of the choice
-		} // End of the Enery LOOP
-  // int np=0;
-  //for(int m=0; m<2; m++)
-  //{
-  //  //fileNuMuX<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
-  //  //G4cout<<"**** A_high="<<m<<G4endl;
-  //  np=0;
-  //  int nC=14;
-	 //  for(G4int en=0; en<nC; en++)
-	 //  {
-  //    //G4double sig=1.;
-  //    np++;
-  //    //if(np==7)           // Write by 7 number in brackets
-  //    //{
-  //    //  if(en==nC-1) fileNuMuX<<sig<<"};"<<G4endl;
-  //    //  else         fileNuMuX<<sig<<","<<G4endl<<"    ";
-	 //    //}
-  //    //else           fileNuMuX<<sig<<",";
-  //    //if(np==7) np=0;
-  //  } // End of the point LOOP
-  //} // End of the isotop LOOP
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXQ2Integration_N.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXQ2Integration_N.cc	(revision 968)
+++ 	(revision )
@@ -1,284 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS (x,Q2) approximation is integrated over x & Q2
-// .....................................................[ for a <nucleon> ]
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Oct-07
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-void strucf(int A, double x, double Q2, double& f2, double& xf3, double& fL)
-{
-  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN2=MN*MN;       // M_N^2 in GeV^2
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double Wt=MN+mpi;       // Delta threshold
-  //const double W2t=Wt*Wt;       // Squared Delta threshold
-  const Genfun::LogGamma lGam;
-  static int mA=0;
-  static double mQ2=0., mN, mD, mDel, mU2, mU3, mV;
-  //static double mUU;
-  double N=3., D=0., Del=0., U2=0., U3=0., V=0.;
-  //double UU=0.;
-  if(A==mA && Q2==mQ2)         // Associative memory for acceleration
-		{
-    N  =mN;
-    D  =mD;
-    Del=mDel;
-    U2 =mU2;
-    U3 =mU3;
-    //UU =mUU;
-    V  =mV;
-  }
-  else
-  {
-    double Q=std::sqrt(Q2);
-    N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
-    Del=.077+.393/(1.+5./Q);
-    D=.68*std::pow(1.+.145/Q2,-1.-Del);
-    V=3*(N-1.);
-    double c3=.75;
-    double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
-    U2=(c3+N-3.)*uu;
-    U3=c3*uu;
-    //UU=uu+uu+uu; // @@
-    mA  = A;
-    mQ2 = Q2;
-    mN  = N;
-    mD  = D;
-    mDel=Del;
-    mU2 =U2;
-    mU3 =U3;
-    //mUU =UU; // @@
-    mV  =V;
-  }  
-  // From here the Q2 coefficients are used
-  double x1=std::pow(1.-x,N-2.);
-  double pp=D*std::pow(x,-Del)*x1;
-  double dir=(1-D)*V*x*x1;
-  double per=U2*pp;
-		f2 = per + dir;
-  //double W2=MN2-MN2*x+Q2/x-Q2;
-  //if(W2<W2t)
-		//{
-  //  per=UU*pp;
-  //  xf3= per+dir;
-  //}
-  //else
-  xf3= U3*pp+dir;
-  fL = per/5.;
-  return;
-}
-
-void getFun(int A, double lx, double Q2, double* f)
-{
-  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  const double MN2=MN*MN;   // Squared Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  double f2=0., xf3=0., fL=0.;
-  if (lx>0.5) G4cerr<<"***getFun: ln(x)="<<lx<<">.5"<<G4endl;
-  double x=std::exp(lx);
-  double x2=x*x;
-  double c=0.;
-  if(Q2>.0001) c=x2*MN2/Q2;
-  double c2=1+c+c;     // Q^2/nu^2 correction
-  strucf(A, x, Q2, f2, xf3, fL);
-  f[0]=f2;             // direct part
-  f[1]=(-f2+xf3)/x;    // *y (neutrino) part
-  f[2]=(-f2-xf3)/x;    // *y (anti-neutrino) part
-  f[3]=(f2*c2-fL-xf3)/x2; // *y2 (neutrino) part
-  f[4]=(f2*c2-fL+xf3)/x2; // *y2 (anti-neutrino) part
-}
-
-int main()
-{
-  const double reps=.001;        // relative accuracy of the total Q2 integral calculation
-  const double xeps=.0001;      // relative accuracy of the total X integral calculation
-  //           =========
-  const double GF=1.16637e-5;   // Fermi constant in GeV^-2
-  const double GF2=GF*GF;       // Squared Fermi constant in GeV^-4
-  const double MW=80.425;       // Mass of W-boson in GeV
-  const double MW2=MW*MW;       // Squared mass of W-boson in GeV^2
-  const double MW4=MW2*MW2;     // Quadro mass of W-boson in GeV^4
-  const double hc2=38937932300.;// (hc)^2 in GeV^2*10^-38cm2 to convert GeV^-2 to 10^-38cm2
-  const double pif=3.14159265*4;// 4pi
-  const double sik=GF2*hc2/pif; // precalculated coefficient
-  //const double mpi=.1349766;    // pi0 meson mass in GeV
-  //const double mpi=.13957018;   // charged pi meson mass in GeV
-  //const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
-  //const double me=.00051099892; // electron mass in GeV
-  //const double me2=me*me;       // m_e^2 in GeV^2
-  //const double hme2=me2/2;      // .5*m_e^2 in GeV^2
-  const double mmu=.105658369;  // mu meson mass in GeV
-  const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
-  const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
-  //const double mtau=1.777;      // tau meson mass in GeV
-  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
-  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
-  //const double mp=.93827203;    // proton mass in GeV
-  //const double mn=.93956536;    // neutron mass in GeV
-  //const double md=1.87561282;   // deuteron mass in GeV
-  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
-  //const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
-  //const double MD=1.232;        // proton mass in GeV
-  //const double mp2=mp*mp;       // m_p^2 in GeV^2
-  const double MN2=MN*MN;       // M_N^2 in GeV^2
-  const double dMN=MN+MN;       // 2*M_N in GeV
-  const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
-  const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
-  //const double EminE=me+me2/dMN;// Threshold for muon production
-  const double EminMu=mmu+mmu2/dMN;  // Threshold for muon production
-  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
-  //
-  const double mc=.261;           // parameter of W2>(M_N^2+M_D^2)/2 cut for QuasiEl/Delta
-  //const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
-  const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
-  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
-  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
-  // _____ Begin of Test Area
-  //Genfun::LogGamma logGamma;
-  //double n=4.9;
-		//double g=exp(logGamma(n));
-  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
-  // ^^^^^ End of Test Area
-  //
-  double f[5];                    // A working array
-  int    A=12;                    // Neucleus for which calculations should be done
-  double lEnuMin=0;               // LogLog of Minimum energy of neutrino
-  double lEnuMax=std::log(1.+std::log(300./EminMu)); // LogLog of MaximumEnergy of neutrino
-  int    nE=63;
-  double dlE=(lEnuMax-lEnuMin)/nE;
-  lEnuMin+=dlE/10;
-  lEnuMax+=dlE/5;
-  G4cout<<"Emin="<<EminMu<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
-  for(double lEnu=lEnuMin; lEnu<lEnuMax; lEnu+=dlE)
-		{
-    double Enu=std::exp(std::exp(lEnu)-1.)*EminMu; // Energy of neutrino/anti-neutrino
-    double dEnu=Enu+Enu;          // doubled energy of nu/anu
-    double Enu2=Enu*Enu;          // squared energy of nu/anu
-    double Emu=Enu-mmu;           // Free Energy of neutrino/anti-neutrino
-    double Emu2=Emu*Emu;          // squared energy of nu/anu
-    double ME=Enu*MN;             // M*E
-    double dME=ME+ME;             // 2*M*E
-    double DIStsig=1.;            // Total curent DIS cross-section to be integrated
-    double DISmsig=1.e20;         // Total remembered DIS cross-section
-    double dEMN=(dEnu+MN)*ME;
-    double MEm=ME-hmmu2;
-    double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
-    double E2M=MN*Enu2-(Enu+MN)*hmmu2;
-    double ymax=(E2M+sqE)/dEMN;
-    double ymin=(E2M-sqE)/dEMN;
-    double rmin=1.-ymin;
-    double rhm2E=hmmu2/Enu2;
-    double Q2min=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
-    double Q2max=dME*ymax;
-    int    nQ2=8;
-    //G4cout<<"*** E="<<Enu<<", Q2i="<<Q2min<<" < Q2a="<<Q2max<<", yi="<<ymin<<" < ya="
-    //      <<ymax<<G4endl;
-    while(std::fabs(DIStsig-DISmsig)/DIStsig>reps)
-    {
-      DISmsig=DIStsig;
-      DIStsig=0.;
-      nQ2*=2;
-      double dQ2=(Q2max-Q2min)/nQ2;
-      for(double Q2=Q2min+dQ2/2; Q2<Q2max; Q2+=dQ2)
-				  {
-        double DISxint=1.;        // Curent DIS x-integral
-        double DISmint=1.e20;     // Remembered DIS x-integral
-        double Q2M=Q2+MW2;
-        double dik=MW4/Q2M/Q2M;
-        double qmc=Q2+mcV;
-        double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2); // QuasielastBoundary
-        //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c)
-        //double xN=Q2/dME;
-        double xN=Q2/MN/(Emu+std::sqrt(Emu2+Q2));
-        //double lXmin=log(xN/ymax);
-        double lXmin=std::log(xN);
-        // ****** QE ********
-        if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES
-        // *** End of QE^^^^^
-        double lXmax=0.;               // QES is in DIS
-        //double lXmax=lXQES;          // Cut off quasielastic
-        int    nX=8;
-        while(std::fabs(DISxint-DISmint)/DISxint>xeps)
-        {
-          DISmint=DISxint;
-          DISxint=0.;
-          nX*=2;
-          double dlX=(lXmax-lXmin)/nX;
-          for(double lX=lXmin+dlX/2; lX<lXmax; lX+=dlX)
-				      {
-            getFun(A, lX, Q2, f);
-            DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
-            //DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
-            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
-          }
-          DISxint*=dlX;
-          //G4cout<<"--- E="<<Enu<<" --- Q2="<<Q2<<" --- nX="<<nX<<", iX="<<DISxint
-          //      <<", mX="<<DISmint<<", rX="<<(DISxint-DISmint)/DISxint<<G4endl;
-        }
-        //G4cout<<"(E="<<Enu<<"), Q2="<<Q2<<", I="<<DISxint/dik/dik<<G4endl;
-        DIStsig+=DISxint*dik;
-      }
-      DIStsig*=dQ2;
-      //G4cout<<"=== E="<<Enu<<" ===> nQ="<<nQ2<<", iQ="<<DIStsig<<", mQ="<<DISmsig
-						//      <<", rQ="<<(DIStsig-DISmsig)/DIStsig<<G4endl;
-    }
-    //===== tot/qe choice ====
-    //DIStsig*=sik/Enu;
-    //G4cout<<"***total*** E="<<Enu<<",sig/E= "<<DIStsig<<G4endl;
-    //...................
-    DIStsig*=sik;
-    G4cout<<"***qelas*** E="<<Enu<<",sig= "<<DIStsig<<G4endl;
-    //===== End of the choice
-		} // End of the Enery LOOP
-  // int np=0;
-  //for(int m=0; m<2; m++)
-  //{
-  //  //fileNuMuX<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
-  //  //G4cout<<"**** A_high="<<m<<G4endl;
-  //  np=0;
-  //  int nC=14;
-	 //  for(G4int en=0; en<nC; en++)
-	 //  {
-  //    //G4double sig=1.;
-  //    np++;
-  //    //if(np==7)           // Write by 7 number in brackets
-  //    //{
-  //    //  if(en==nC-1) fileNuMuX<<sig<<"};"<<G4endl;
-  //    //  else         fileNuMuX<<sig<<","<<G4endl<<"    ";
-	 //    //}
-  //    //else           fileNuMuX<<sig<<",";
-  //    //if(np==7) np=0;
-  //  } // End of the point LOOP
-  //} // End of the isotop LOOP
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXYIntegration.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXYIntegration.cc	(revision 968)
+++ 	(revision )
@@ -1,241 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS (x,Q2) approximation is integrated over x
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 30-Sept-05
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-void strucf(int A, double x, double ydME, double& f2, double& xf3, double& fL)
-{
-  const Genfun::LogGamma lGam;
-  double N=3., D=0., Del=0., U2=0., U3=0., V=0.;
-  double Q2=x*ydME;
-  double r=0.;
-  double max=1.;
-  double H=1.22;
-  if(A==1)                   // Proton
-  {
-    r=std::sqrt(Q2/1.66);
-    max=.5;
-  }
-  else if(A<13)              // Light nuclei
-  {
-    double f=Q2/4.62;
-    r=f*f;
-    max=.3;
-    if(A>2) H=1.;
-  }
-  else if(A>0)               // Heavy nuclei
-  {
-    double f=Q2/3.4;
-    double ff=f*f;
-    r=ff*ff;
-    max=.5;
-    H=1.;
-  }
-  else G4cout<<"strucf: A="<<A<<" <= 0"<< G4endl;
-  //
-  N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
-  Del=(1.+r)/(12.5+r/max);
-  double S=std::pow(1.+.6/Q2,-1.-Del);
-  D=H*S*(1.-.5*S);
-  V=3*(1.-D)*(N-1.);
-  double cc=Q2/.08;
-  double cc2=cc*cc;
-  double C=(1.+cc2)/(1.+cc2/.24)/(1.+Q2/21.6);
-  double c3=C+C+C;
-  double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
-  U2=(c3+N-3.)*uu;
-  U3=c3*uu;
-  // From here the Q2 coefficients are used
-  double x1=std::pow(1.-x,N-2.);
-  double pp=D*std::pow(x,-Del)*x1;
-  double dir=V*x*x1;
-  double per=U2*pp;
-		f2 = per + dir;
-  xf3= U3*pp+dir;
-  fL = per/4.;
-  return;
-}
-
-void getFun(int A, double lx, double ydME, double* f)
-{
-  double f2=0., xf3=0., fL=0.;
-  if (lx>0.5) G4cerr<<"***getFun: ln(x)="<<lx<<">.5"<<G4endl;
-  double x=std::exp(lx);
-  strucf(A, x, ydME, f2, xf3, fL);
-  f[0]=x*f2;          // direct part
-  f[1]=x*(-f2+xf3);   // *y (neutrino) part
-  f[2]=x*(-f2-xf3);   // *y (anti-neutrino) part
-  f[3]=x*(f2-fL-xf3); // *y2 (neutrino) part
-  f[4]=x*(f2-fL+xf3); // *y2 (anti-neutrino) part
-}
-
-int main()
-{
-  const double reps=.01;        // relative accuracy of the total Q2 integral calculation
-  const double xeps=.0001;      // relative accuracy of the total X integral calculation
-  //           =========
-  const double GF=1.16637e-5;   // Fermi constant in GeV^-2
-  const double GF2=GF*GF;       // Squared Fermi constant in GeV^-4
-  const double MW=80.425;       // Mass of W-boson in GeV
-  const double MW2=MW*MW;       // Squared mass of W-boson in GeV^2
-  const double MW4=MW2*MW2;     // Quadro mass of W-boson in GeV^4
-  const double hc2=38937932300.;// (hc)^2 in GeV^2*10^-38cm2 to convert GeV^-2 to 10^-38cm2
-  const double pif=3.14159265*4;// 4pi
-  const double sik=GF2*MW4*hc2/pif; // precalculated coefficient
-  const double mpi=.1349766;    // pi0 meson mass in GeV
-  const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
-  const double me=.00051099892; // electron mass in GeV
-  const double me2=me*me;       // m_e^2 in GeV^2
-  const double hme2=me2/2;      // .5*m_e^2 in GeV^2
-  const double mmu=.105658369;  // mu meson mass in GeV
-  const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
-  const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
-  const double mtau=1.777;      // tau meson mass in GeV
-  const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
-  const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
-  const double mp=.93827203;    // proton mass in GeV
-  const double mn=.93956536;    // neutron mass in GeV
-  const double md=1.87561282;   // deuteron mass in GeV
-  const double MN=.931494043;   // Nucleon mass (inside nucleus) in GeV (atomic mass unit)
-  const double mp2=mp*mp;       // m_p^2 in GeV^2
-  const double MN2=MN*MN;       // M_N^2 in GeV^2
-  const double dMN=MN+MN;       // 2*M_N in GeV
-  const double EminE=me+me2/dMN;// Threshold for muon production
-  const double EminMu=mmu+mmu2/dMN;  // Threshold for muon production
-  const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
-  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
-  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
-  // _____ Begin of Test Area
-  //Genfun::LogGamma logGamma;
-  //double n=4.9;
-		//double g=exp(logGamma(n));
-  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
-  // ^^^^^ End of Test Area
-  //
-  double f[5];                    // A working array
-  int    A=12;                    // Neucleus for which calculations should be done
-  double lEnuMin=std::log(EminMu);     // Log of Minimum energy of neutrino
-  double lEnuMax=std::log(300.);       // Log of Maximum energy of neutrino
-  int    nE=20;
-  double dlE=(lEnuMax-lEnuMin)/nE;
-  G4cout<<"Emin="<<EminMu<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
-  for(double lEnu=lEnuMin+dlE/2; lEnu<lEnuMax; lEnu+=dlE)
-		{
-    double Enu=std::exp(lEnu);         // Energy of neutrino/anti-neutrino
-    double dEnu=Enu+Enu;          // doubled energy of nu/anu
-    double Enu2=Enu*Enu;          // squared energy of nu/anu
-    double ME=Enu*MN;             // M*E
-    double dME=ME+ME;             // 2*M*E
-    double DIStsig=1.;            // Total curent DIS cross-section to be integrated
-    double DISmsig=1.e20;         // Total remembered DIS cross-section
-    double dEMN=(dEnu+MN)*ME;
-    double MEm=ME-hmmu2;
-    double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
-    double E2M=MN*Enu2-(Enu+MN)*hmmu2;
-    double ymax=(E2M+sqE)/dEMN;
-    double ymin=(E2M-sqE)/dEMN;
-    // The Q2min/max calculation is not necessary in this method
-    //double rmin=1.-ymin;
-    //double rhm2E=hmmu2/Enu2;
-    //double Q2min=(Enu2+Enu2)*(rmin-rhm2E-sqrt(rmin*rmin-rhm2E-rhm2E));
-    //double Q2max=dME*ymax;
-    int    nY=8;
-    //G4cout<<"*** E="<<Enu<<", N="<<nY<<", yi="<<ymin<<" < ya="<<ymax<<G4endl;
-    while(std::fabs(DIStsig-DISmsig)/DIStsig>reps)
-    {
-      DISmsig=DIStsig;
-      DIStsig=0.;
-      nY*=2;
-      double dY=(ymax-ymin)/nY;
-      for(double y=ymin+dY/2; y<ymax; y+=dY)
-				  {
-        double y2=y*y;
-        double DISxint=1.;        // Curent DIS x-integral
-        double DISmint=1.e20;     // Remembered DIS x-integral
-        double lXmin=std::log((std::sqrt(Enu2+mmu2)-Enu)/MN);
-        double lXmax=0.;          // Temporary, while QES is in DIS
-        double ydME=y*dME;
-        int    nX=8;
-        while(std::fabs(DISxint-DISmint)/DISxint>xeps)
-        {
-          DISmint=DISxint;
-          DISxint=0.;
-          nX*=2;
-          double dlX=(lXmax-lXmin)/nX;
-          for(double lX=lXmin+dlX/2; lX<lXmax; lX+=dlX)
-				      {
-            getFun(A, lX, ydME, f);
-            double dik=std::exp(lX)*ydME+MW2;
-            //DISxint+=(f[0]+f[0]+(y+y)*f[1]+y2*f[3])/dik/dik; // neutrino
-            DISxint+=(f[0]+f[0]+(y+y)*f[2]+y2*f[4])/dik/dik; // anti-neutrino
-            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
-          }
-          DISxint*=dlX;
-          //G4cout<<"--- E="<<Enu<<" --- Q2="<<Q2<<" --- nX="<<nX<<", iX="<<DISxint
-          //      <<", mX="<<DISmint<<", rX="<<(DISxint-DISmint)/DISxint<<G4endl;
-        }
-        //G4cout<<"(E="<<Enu<<"), Q2="<<Q2<<", I="<<DISxint/dik/dik<<G4endl;
-        DIStsig+=DISxint;
-      }
-      DIStsig*=dY;
-      //G4cout<<"=== E="<<Enu<<" ===> nQ="<<nQ2<<", iQ="<<DIStsig<<", mQ="<<DISmsig
-						//      <<", rQ="<<(DIStsig-DISmsig)/DIStsig<<G4endl;
-    }
-    DIStsig*=sik*dMN; // Jakobian is 2ME => instead of 1/E must be *2M
-    G4cout<<"*** E="<<Enu<<",sig/E="<<DIStsig<<G4endl;
-		} // End of the Enery LOOP
-  // int np=0;
-  //for(int m=0; m<2; m++)
-  //{
-  //  //fileNuMuX<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
-  //  //G4cout<<"**** A_high="<<m<<G4endl;
-  //  np=0;
-  //  int nC=14;
-	 //  for(G4int en=0; en<nC; en++)
-	 //  {
-  //    //G4double sig=1.;
-  //    np++;
-  //    //if(np==7)           // Write by 7 number in brackets
-  //    //{
-  //    //  if(en==nC-1) fileNuMuX<<sig<<"};"<<G4endl;
-  //    //  else         fileNuMuX<<sig<<","<<G4endl<<"    ";
-	 //    //}
-  //    //else           fileNuMuX<<sig<<",";
-  //    //if(np==7) np=0;
-  //  } // End of the point LOOP
-  //} // End of the isotop LOOP
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4PhotoNuclearCalculation.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4PhotoNuclearCalculation.cc	(revision 968)
+++ 	(revision )
@@ -1,261 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: PhotoNuclearCalculation of gamma+A cross sections
-// according to the paper M.V.Kossov "Approximation of photonuclear
-// interaction cross sections" (EPJA,2002, be published). The cross
-// sections are calculated according to the exact complicated formulas
-// and tabulated for the fast response. The total energy range is
-// subdivided in three regions:
-// a) The GDR region (from the hadron production threshold to 106 MeV)
-//    covering 46 nuclei: H2,He4,Li6,Li7,Be9,C12,N14,N15,O16,F19,Na23,
-//    Mg24,Al27,Si28,S32,S34,(Ar40,Ca40),Fe54,Mn55,Fe56,Ni58,Co59,Cu,Zn
-//    Se76,Se82,Ag,Cd,Sn,I,Sm154,Gd156,Tb159,Ho165,Er168,Yb174,Hf178,Hf180,
-//    Ta181,W184,W186,Au197,Tl,Pb,Bi209,Th232,U235,U238,Pu239. For these
-//    isotops the approximation fits the existing measurements. For
-//    them the basic functions are calculated. For other isotops the
-//    approximation is interpolated linearly with A.
-// b) The resonance region (from 106 MeV to 50 GeV)
-//    covering 14 nuclei: H1,H2,He3,He4,Li6,Li7,Be,C,O,Al,Cu,Sn,Pb,U.
-//    For these isotops the melting of resonances in nuclear matter is
-//    measured and approximated. For them the basic functions are
-//    calculated. For other isotops the approximation is interpolated
-//    linearly with A.
-// c) The high energy region is an A-dependent function for E>100 GeV,
-//    taking into account nuclear shaddowing.
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-May-02
-// The lust update: M.V. Kossov, CERN/ITEP(Moscow) 17-June-02
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-
-int main()
-{
-  const G4int mN=49;                // A#of GDR basic nuclei
-  const G4int mC=105;               // A#of Resonance points in E (each MeV from 2 to 106)
-  const G4double mA[mN]={
-	2.,4.,6.,7.,9.,12.,14.,15.,16.,19.,23.,24.,27.,28.,32.,34.,40.,54.,
-    55.,56.,58.7,58.9,63.5,65.4,76.,82.,107.9,112.4,118.7,126.9,154.,156.,159.,165.,
-    168.,174.,178.,180.,181.,184.,186.,197.,204.4,207.2,209.,232.,235.,
-    238.,239.};
-  const G4int nN=14;                // A#of Resonance basic nuclei
-  const G4int nC=224;               // A#of Resonance points in lnE
-  const G4double nA[nN]={1.,2.,3.,4.,6.,7.,9.,12.,16.,27.,63.5,118.7,207.2,238.};
-  const G4double r4[6]={0.,3.38,3.51,3.23,3.57,3.60};
-  const G4double t4[6]={0.,3.14,3.09,2.80,3.25,3.17};
-  const G4double r8[mN]={
-    0.00,0.00,0.00,2.57,3.26,3.52,3.51,3.51,3.49,3.05,3.20,3.25,3.40,3.42,3.38,3.50,3.51,
-    3.59,3.45,3.46,3.52,3.56,3.40,3.40,3.29,3.38,3.50,3.49,3.41,3.46,3.40,3.37,3.40,3.29,
-    3.37,3.21,3.45,3.45,3.39,3.36,3.36,3.39,3.47,3.40,3.36,3.31,3.30,3.37,3.36};
-  const G4double t8[mN]={
-    0.00,0.00,0.00,2.48,3.05,3.11,3.09,3.08,3.10,2.70,2.90,2.95,3.00,2.97,2.98,3.00,3.05,
-    2.94,2.93,2.86,2.99,3.01,2.85,2.84,2.72,2.76,2.81,2.77,2.72,2.75,2.78,2.74,2.72,2.72,
-    2.70,2.74,2.70,2.70,2.64,2.63,2.61,2.59,2.63,2.57,2.63,2.62,2.62,2.67,2.65};
-  const G4double iE=std::log(106.);      // Start logarithm energy
-  const G4double fE=std::log(50000.);    // Finish logarithm energy (each 2.75 percent)
-  const G4double dE=(fE-iE)/(nC-1); // Step in logarithm energy
-  std::ofstream fileGDR("GDR.out", std::ios::out);
-  fileGDR.setf( std::ios::scientific, std::ios::floatfield );
-  std::ofstream fileRes("Res.out", std::ios::out);
-  fileRes.setf( std::ios::scientific, std::ios::floatfield );
-  G4int np=0;
-  for(G4int m=0; m<mN; m++)
-  {
-    G4double A=mA[m];
-    G4double lnA=std::log(A);
-    G4double A2=A*A;
-    G4double red=1.+16/A2/A2;
-    G4double rho1=(3.2+.75*lnA)/red;
-    G4double tau1=(6.6-.5*lnA)/red;
-	if(A==2.)
-	{
-	  rho1=1.86;
-	  tau1=1.2;
-	}
-	G4double rho2=(4.+.125*lnA)/red;
-	G4double tau2=3.4/red;
-	if(A==6.)
-	{
-	  rho2=2.9;
-	  tau2=2.32;
-	}
-	else if(A==2.)
-	{
-	  rho2=2.11;
-	  tau2=1.5;
-	}
-	G4double rho4=3.8+.05*lnA;
-	G4double tau4=3.8-.25*lnA;
-	if(A<13.)
-	{
-	  rho4=r4[m];
-	  tau4=t4[m];
-	}
-	G4double rho8=r8[m];
-	G4double tau8=t8[m];
-    G4double tr=5.13-.00075*A; // Resonance threshold Position
-    //if(A==1.)tr=5.24;
-	G4double rr=11.;           // Resonance threshold Slope
-	if(A<2.5)rr=25.;
-	G4double dw=.056+lnA*(.03-.001*lnA); // Delta width
-	G4double du=5.82-.07/(1.+.003*A2);   // Delta position
-    G4double da=.39*A;                   // Delta amplitude
-    if(A==2.)da=.88;
-    //else if(A==1.) da=.55;
-    G4double ekin=1.;
-    fileGDR<<"  static const G4double SL"<<m<<"[nL]={"<<G4endl<<"    ";
-    G4cout<<"**** A_low="<<A<<G4endl;
-    np=0;
-	for(G4int em=0; em<mC; em++)
-	{
-      ekin+=1.;
-      G4double z=std::log(ekin);
-      G4double ds=z-du;
-	  G4double sigd=da/(1.+ds*ds/dw)/(1.+std::exp(rr*(tr-z))); // Delta contribution
-	  G4double g1=std::exp(rho1-z)/(1.+std::exp(3*(tau1-z)));
-	  G4double g2=std::exp(2*(rho2-z))/(1.+std::exp(6*(tau2-z)));
-      G4double g4=0.;
-      if(A>3.5)g4=std::exp(4*(rho4-z))/(1.+std::exp(12*(tau4-z)));
-      G4double g8=0.;
-      if(A>6.5)g8=std::exp(8*(rho8-z))/(1.+std::exp(24*(tau8-z)));
-      G4double sig=sigd+g1+g2+g4+g8;
-      np++;
-      if(np==7)
-      {
-        if(em==mC-1) fileGDR<<sig<<"};"<<G4endl;
-        else         fileGDR<<sig<<","<<G4endl<<"    ";
-	  }
-      else      fileGDR<<sig<<",";
-      if(np==7) np=0;
-    } // End of the point LOOP
-  } // End of the isotop LOOP
-  // ===== High Energy Calculations =======================
-  //G4double shd=1.0663-.0023*log(2.);   // HE PomShadowing(D)
-  np=0;
-  for(G4int n=0; n<nN; n++)
-  {
-    G4double A=nA[n];
-    G4double A2=A*A;
-    fileRes<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
-    G4cout<<"**** A_high="<<A<<G4endl;
-    G4double lnA=std::log(A);
-    G4double slA=std::sqrt(lnA);
-    G4double red=1.+16/A2/A2;
-    G4double rho1=(3.2+.75*lnA)/red;
-    G4double tau1=(6.6-.5*lnA)/red;
-	if(A==2.)
-	{
-	  rho1=1.86;
-	  tau1=1.2;
-	}
-	G4double rho2=(4.+.125*lnA)/red;
-	G4double tau2=3.4/red;
-	if(A==6.)
-	{
-	  rho2=2.9;
-	  tau2=2.32;
-	}
-	else if(A==2.)
-	{
-	  rho2=2.11;
-	  tau2=1.5;
-	}
-    G4double rho4=6.27;
-    G4double tau4=7.25;
-    G4double rho8=6.66;
-    G4double tau8=6.90;
-    if(A==2.)
-	{
-     rho4=6.2;
-     tau4=7.1;
-     rho8=6.62;
-     tau8=6.91;
-    }
-    G4double tr=5.13-.00075*A; // Resonance threshold Position
-    if(A==1.)tr=5.24;
-	G4double rr=11.;           // Resonance threshold Slope
-	if(A<2.5)rr=25.;
-	// Delta
-	G4double dw=.056+lnA*(.03-.001*lnA); // Delta width
-	G4double du=5.82-.07/(1.+.003*A2);   // Delta position
-    G4double da=.39*A;                   // Delta amplitude
-    if(A==2.)da=.88;
-    else if(A==1.) da=.55;
-	// High Resonance
-	G4double hw=.045+.04*slA*slA*slA;    // HighR width
-	G4double hu=6.496+.042*lnA;          // HighR position
-    if(A==1.)hu=6.57;
-    else if(A==2.)hu=6.575;
-    G4double ha=.223;                    // HighR amplitude
-    if(A>2.5)ha=.16*A/slA;
-    else if(A==2.) ha=.348;
-    G4double sp=A*(1.-.072*lnA);         // HE TotShadowing
-    G4double sh=1.0663-.0023*lnA;        // HE PomShadowing
-    if(A==1.)sh=1.07;
-	// -- Loop over energies ---
-    G4double z=iE-dE;
-    np=0;
-	for(G4int en=0; en<nC; en++)
-	{
-      z+=dE;
-      G4double ds=z-du;
-      G4double fr=1.+std::exp(rr*(tr-z));
-	  G4double sigd=da/(1.+ds*ds/dw); // Delta contribution
-      G4double hs=z-hu;
-	  G4double sigh=ha/(1.+hs*hs/hw); // HighR contribution
-	  G4double g1=0.;
-      if(A>1.5)g1=std::exp(rho1-z)/(1.+std::exp(3*(tau1-z)));
-	  G4double g2=0.;
-      if(A>1.5)g2=std::exp(2*(rho2-z))/(1.+std::exp(6*(tau2-z)));
-      G4double g4=0.;
-      if(A<2.5)g4=std::exp(4*(rho4-z))/(1.+std::exp(12*(tau4-z)));
-      G4double g8=0.;
-      if(A<2.5)g8=std::exp(8*(rho8-z))/(1.+std::exp(24*(tau8-z)));
-      G4double hp=.0375*(z-16.5)+sh*std::exp(-.11*z);
-      G4double fp=hp/(1.+std::exp(4*(7.-z)));
-      G4double sig=(sigd+sigh)/fr+g1+g2+g4+g8+sp*fp;
-      np++;
-      if(np==7)
-      {
-        if(en==nC-1) fileRes<<sig<<"};"<<G4endl;
-        else         fileRes<<sig<<","<<G4endl<<"    ";
-	  }
-      else      fileRes<<sig<<",";
-      if(np==7) np=0;
-      //if(en==nC)
-	  //{
-      //  G4double fsig=sp*(.0375*(z-16.5)+shd*exp(-.11*z));
-      //  fileRes<<">>> fun="<<fsig<<G4endl;
-	  //}
-    } // End of the point LOOP
-  } // End of the isotop LOOP
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4PhotoNuclearCrossSectionTest.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4PhotoNuclearCrossSectionTest.cc	(revision 968)
+++ 	(revision )
@@ -1,84 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-#include "G4ios.hh"
-#include "G4Gamma.hh"
-#include "G4Element.hh"
-#include "G4PhotoNuclearCrossSection.hh"
-
-int main()
-{
-  G4PhotoNuclearCrossSection gammaACrossSection;
-  const G4int nN=6;                          // A#of plots (nuclei in the test)
-  const G4int nC=227;                        // A#of points in the plots (energy values)
-  ofstream fileH1("h1.out", ios::out);
-  fileH1.setf( ios::scientific, ios::floatfield );
-  ofstream fileH2("h2.out", ios::out);
-  fileH2.setf( ios::scientific, ios::floatfield );
-  ofstream fileC12("c12.out", ios::out);
-  fileC12.setf( ios::scientific, ios::floatfield );
-  ofstream fileAl27("al27.out", ios::out);
-  fileAl27.setf( ios::scientific, ios::floatfield );
-  ofstream fileCu("cu.out", ios::out);
-  fileCu.setf( ios::scientific, ios::floatfield );
-  ofstream filePb("pb.out", ios::out);
-  filePb.setf( ios::scientific, ios::floatfield );
-  G4double low[nN]={137. , 2.22, 12.5 , 8.3 , 7.  , 7.35 };
-  G4double high[nN]={4.e7, 4.e7, 4.e7, 4.e7, 4.e7, 4.e7};
-  G4Element* theElement[nN]={new G4Element("Hydrogen", "H", 1, 1.*g/mole),
-							 new G4Element("Deuterium", "D", 1, 2.*g/mole),
-							 new G4Element("Carbon", "C", 6, 12.*g/mole),
-							 new G4Element("Aluminum", "Al", 13, 27.*g/mole),
-							 new G4Element("Copper", "Cu", 29, 64.*g/mole),
-							 new G4Element("Lead", "Pb", 82, 208.*g/mole)};
-  G4ParticleDefinition* theParticleDefinition = G4Gamma::GammaDefinition();
-  //G4VProcess* theProcess = new G4PhotonInelasticProcess("Inelastic"); // @@ Make it?! 
-  G4DynamicParticle* theDynamicParticle;
-  G4double sig=0.;
-  for(G4int n=0; n<6; n++)
-  {
-	G4double lekin = std::log(low[n]);
-	G4double dlekin= std::exp((std::log(std::log(high[n]))-std::log(lekin))/(nC-1));
-	lekin /= dlekin;
-	//G4cout<<"n="<<n<<", low="<<low[n]<<", high="<<high[n]<<", np="<<nC<<", d="<<dlekin<<G4endl;
-	for(G4int ll=0; ll<nC; ll++)
-	{
-	  lekin*=dlekin;
-      G4double ekin=std::exp(lekin);
-	  theDynamicParticle = new G4DynamicParticle(theParticleDefinition,
-												 G4ParticleMomentum(1.,0.,0.), ekin*MeV);
-	  sig = gammaACrossSection.GetCrossSection(theDynamicParticle,theElement[n])/millibarn;
-	  delete theDynamicParticle;
-      //G4cout<<"n="<<n<<", e="<<ekin<<", sigma="<<sig<<G4endl;
-      if(!n) fileH1<<ekin<<" "<<sig<<G4endl;
-      else if(n==1) fileH2<<ekin<<" "<<sig<<G4endl;
-      else if(n==2) fileC12<<ekin<<" "<<sig<<G4endl;
-      else if(n==3) fileAl27<<ekin<<" "<<sig<<G4endl;
-      else if(n==4) fileCu<<ekin<<" "<<sig<<G4endl;
-      else if(n==5) filePb<<ekin<<" "<<sig<<G4endl;
-    } // End of the point LOOP
-  } // End of the plot LOOP
-  return EXIT_SUCCESS;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4PolinomTabulation.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4PolinomTabulation.cc	(revision 968)
+++ 	(revision )
@@ -1,204 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS(Q2) fixed step integration
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-2005
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-//#include <CLHEP/GenericFunctions/LogGamma.hh>
-
-// Solution for randomization of the high order polinoms is tabulated
-
-double poli(double x, int n) // f=n*x^(n-1)-(n-1)*x^n=R(0-1), n>2, as n=2 has analSolution
-{
-  double p=x;
-  int n1=n-1;
-  for(int i=1; i<n1; i++) p*=x;
-  return p*(n-n1*x);
-}
-
-double dpoli(double x, int n) // df=n*(n-1)*x^(n-2)-(n-1)*n*x^(n-1), n>2 (n=2 has analSol)
-{
-  double p=x;
-  int n1=n-1;
-  if(n>3) for(int i=2; i<n1; i++) p*=x;
-  return p*n*n1*(1.-x);
-}
-// One can define for tabulation more similar functions
-
-int main()
-{
-  const double eps=.0000001;      // relative accuracy of the integral calculation
-  //           =============
-  int nSub=10;
-  double dSub=1./nSub;
-  for(int n=3; n<12; n++)
-		{
-    G4cout<<"n="<<n<<G4endl;
-    G4double s=0.;
-    for(int i=1; i<nSub; i++)
-    {
-      s+=dSub;
-      G4double r=std::pow(s,(n-1.5)/1.5);
-      G4double x=0.5;
-		    if(n==3)
-		    {
-        if    (r==0.5) x=0.5;
-        else if(r<0.5) x=std::sqrt(r+r)*(.5+.1579*(r-.5));
-        else           x=1.-std::sqrt(2.-r-r)*(.5+.1579*(.5-r));
-      }
-      else
-		    {
-        G4int n1=n-1;
-        G4double r1=n1;
-        G4double r2=r1-1.;
-        G4double rr=r2/r1;
-        G4double rp=std::pow(rr,n1);
-        G4double p2=rp+rp;
-        if  (r==rr)  x=p2;
-        else
-				    {
-          if(r<rr)
-          {
-								    G4double pr=0.;
-								    G4double pra=0.;
-            if(n>7)
-								    {
-              if(n>9)
-								      {
-                if(n>10)                         // >10(11)
-                {
-                  pr=.614/std::pow((n+1+1.25),.75);
-                  pra=.915/std::pow((n+1+6.7),1.75);
-                }
-												    else                             // 10
-                {
-                  pr=.09945;
-                  pra=.00667;
-                }
-              }
-              else
-								      {
-                if(n>8)                          // 9
-                {
-                  pr=.1064;
-                  pra=.00741;
-                }
-												    else                             // 8
-                {
-                  pr=.11425;
-                  pra=.00828;
-                }
-              }
-            }
-            else
-								    {
-              if(n>5)
-								      {
-                if(n>6)                          // 7
-                {
-                  pr=.12347;
-                  pra=.00926;
-                }
-												    else                             // 6
-                {
-                  pr=.13405;
-                  pra=.01027;
-                }
-              }
-              else
-								      {
-                if(n>4)                          // 5
-                {
-                  pr=.1454;
-                  pra=.01112;
-                }
-												    else                             // 4
-                {
-                  pr=.15765;
-                  pra=.00965;
-                }
-              }
-            }
-            x=std::pow((r/p2),(1.-rr+pra))*(rr+pr*(r-p2));
-          }
-          else
-          {
-								    G4double sr=0.;
-            if(n>7)
-								    {
-              if(n>9)
-								      {
-                if(n>10) sr=.86/(n+1+1.05);      // >10(11)
-												    else     sr=.0774;               // 10
-              }
-              else
-								      {
-                if(n>8) sr=.0849;                // 9
-												    else    sr=.0938;                // 8
-              }
-            }
-            else
-								    {
-              if(n>5)
-								      {
-                if(n>6) sr=.1047;                // 7
-												    else    sr=.1179;                // 6
-              }
-              else
-								      {
-                if(n>7) sr=.1339;                // 5
-												    else    sr=.15135;               // 4
-              }
-            }
-            x=1.-std::sqrt((1.-r)/(1.-p2))*(1.-rr+sr*(p2-r));
-          }
-        }
-      }
-      G4double dx=x;
-      G4double f=poli(x,n)-r;
-      G4int it=0;
-      while(std::fabs(f/r)>eps)
-						{
-        it++;
-        G4double df=dpoli(x,n);
-        //G4cout<<"n="<<n<<", r="<<r<<", f="<<f<<", d="<<df<<", x="<<x<<", i="<<it<<G4endl;
-        x-=f/df;
-        f=poli(x,n)-r;
-      }
-      dx=std::fabs(dx-x);
-      G4double d=std::fabs(f/r);
-      G4cout<<"n="<<n<<",r="<<r<<", Final: f="<<d<<",x="<<x<<",d="<<dx<<",it="<<it<<G4endl;
-    } // End of loop over r
-  } // End of loop over n
-  return 0;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4SurfaceCalculation.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4SurfaceCalculation.cc	(revision 968)
+++ 	(revision )
@@ -1,125 +1,0 @@
-//
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-// G4 Tools program: NuMu DIS(Q2) fixed step integration
-// .....................................................
-// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-2005
-// 
-//=====================================================================
-#include "globals.hh"
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include "G4ios.hh"
-#include "Randomize.hh"
-#include "G4Timer.hh"
-#include "time.h"
-
-// All calculations have been done for C12 nucleus
-
-double anuX(double E, double r, double p) // (E is in GeV, r=Q2/Q2max, p=1+anuE(E))
-{
-		double	E2=E*E;
-  double y=p-r;
-  double p3=(13.88+.9373*(1.+.000033*E2)*std::sqrt(E))/(1.+(10.12+1.532/E2)/E);
-  return std::pow(y,6)*(r+p3)/(p3*r+1.);
-}
-
-int main()
-{
-  G4double PI=3.14159265;        // \pi
-  G4double A=.5;                 // half side of the tested cube
-  G4double R=.5;                 // radius of the tested sphere
-  G4double B=1.;                 // half width of the mother cube
-  G4double S=6*(B+B)*(B+B);      // 6 sides 2x2
-  G4double SC=6*(A+A)*(A+A);     // The target surface of the test cube
-  G4double SS=4*PI*R*R;          // The target surface of the test sphere
-  G4double c=0.;                 // summ for the cube							
-  G4double s=0.;                 // summ for the sphere					
-  G4double m=0.;                 // summ for the mother cube
-  G4double cf=0.;                // reversed length summ for the cube							
-  G4double sf=0.;                // reversed length summ for the sphere					
-  G4double mf=0.;                // reversed length summ for the mother cube
-  G4double cr=0.;                // length summ for the cube							
-  G4double sr=0.;                // length summ for the sphere					
-  G4double mr=0.;                // length summ for the mother cube
-  G4double cm=B;                 // min length for the cube							
-  G4double sm=B;                 // min length for the sphere					
-  G4double mm=B;                 // min length for the mother cube
-  G4int nEv=10000000;            // calculation statistics
-  G4double f=.99;                // scale factor
-  G4double Af=A*f;               // Scaled A
-  G4double Rf=R*f;               // Scaled R
-  G4double Bf=B*f;               // Scaled B
-  G4Timer* timer = new G4Timer();
-  timer->Start();
-  for(G4int i=0; i<nEv; i++)
-		{
-    // Randomize coordinates within (+/-1,+/-1,+/-1) cube (S=6+4=24)
-    G4double x=G4UniformRand();
-    x=B*(1.-x-x);                    // Can be x=B*x, but B=1 is a hilf width of mother cube
-    G4double y=G4UniformRand();
-    y=B*(1.-y-y);
-    G4double z=G4UniformRand();
-    z=B*(1.-z-z);
-    G4double r2=x*x+y*y+z*z;
-    G4double r=std::sqrt(r2);
-    G4double w=1./r;
-    G4double ax=std::fabs(x);
-    G4double ay=std::fabs(y);
-    G4double az=std::fabs(z);
-    if(r<R && r>Rf)
-    {
-      s+=1.;
-      sf+=w;
-      sr+=r;
-      if(r<sm) sm=r;
-    }
-    if(ax<A && ay<A && az<A && (ax>Af || ay>Af || az>Af))
-    {
-      c+=1.;
-      cf+=w;
-      cr+=r;
-      if(r<cm) cm=r;
-    }
-    if(ax>Bf || ay>Bf || az>Bf)
-    {
-      m+=1.;
-      mf+=w;
-      mr+=r;
-      if(r<mm) mm=r;
-    }
-  }
-  G4double N=S*mm/m;
-  //G4double N=S*m*m/mf;
-  //G4double N=S/mf;
-  timer->Stop();
-  G4cout<<"SurfaceCalc: nE="<<nEv<<", time="<<timer->GetUserElapsed()<<" s, f="<<f<<G4endl;
-  delete timer;
-  G4cout<<"SurfaceCalculation: Cube="<<N*c/cm/SC<<", Sphere="<<N*s/sm/SS<<G4endl;
-  //G4cout<<"SurfaceCalculation: Cube="<<N*cf/SC<<", Sphere="<<N*sf/SS<<G4endl;
-  //G4cout<<"SurfaceCalculation: Cube="<<N*cr/c/c/SC<<", Sphere="<<N*sr/s/s/SS<<G4endl;
-}
Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/GNUmakefile
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/GNUmakefile	(revision 968)
+++ 	(revision )
@@ -1,57 +1,0 @@
-# $Id: GNUmakefile,v 1.2 2006/06/02 12:37:30 mkossov Exp $
-# ----------------------------------------------------------------
-# Makes the calculation program in environment variable G4TARGET.
-# ----------------------------------------------------------------
-
-ifndef G4TARGET
-  G4TARGET := $(TESTTARGET)
-endif
-
-ifndef G4INSTALL
-  G4INSTALL = ../../../../..
-endif
-
-G4EXEC_BUILD = true
-
-include $(G4INSTALL)/config/architecture.gmk
-
-# Override some variables for binmake.gmk.
-#
-INCFLAGS := \
-            -I$(G4BASE)/global/management/include \
-            -I$(G4BASE)/global/HEPRandom/include \
-            -I$(G4BASE)/global/HEPGeometry/include \
-            -I$(G4BASE)/geometry/management/include \
-            -I$(G4BASE)/geometry/volumes/include \
-            -I$(G4BASE)/geometry/magneticfield/include \
-            -I$(G4BASE)/geometry/solids/CSG/include \
-            -I$(G4BASE)/event/include \
-            -I$(G4BASE)/graphics_reps/include \
-            -I$(G4BASE)/track/include \
-            -I$(G4BASE)/tracking/include \
-            -I$(G4BASE)/processes/management/include \
-            -I$(G4BASE)/materials/include \
-            -I$(G4BASE)/particles/management/include \
-            -I$(G4BASE)/particles/bosons/include \
-            -I$(G4BASE)/particles/leptons/include \
-            -I$(G4BASE)/particles/hadrons/mesons/include \
-            -I$(G4BASE)/particles/hadrons/barions/include \
-            -I$(G4BASE)/particles/hadrons/ions/include \
-            -I$(G4BASE)/processes/general/include \
-            -I$(G4BASE)/processes/management/include \
-            -I$(G4BASE)/processes/hadronic/include \
-            -I$(G4BASE)/processes/hadronic/processes/include \
-            -I$(G4BASE)/processes/hadronic/management/include \
-            -I$(G4BASE)/processes/hadronic/util/include \
-            -I$(G4BASE)/processes/hadronic/cross_sections/include \
-            -I$(G4BASE)/processes/hadronic/models/neutron_hp/include \
-            -I$(G4BASE)/processes/hadronic/models/chiral_inv_phase_space/body/include \
-            -I$(G4BASE)/processes/electromagnetic/standard/include \
-            -I$(G4BASE)/processes/electromagnetic/muons/include \
-            -I$(G4BASE)/processes/electromagnetic/utils/include \
-            -I$(G4BASE)/processes/electromagnetic/xrays/include \
-            -I$(G4BASE)/run/include \
-            -I$(G4BASE)/digits+hits/include \
-            -I$(G4BASE)/intercoms/include
-
-include $(G4INSTALL)/config/binmake.gmk
