G4INCLCrossSectionsMultiPionsAndResonances.cc

이 파일의 문서화 페이지로 가기

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
// INCL++ intra-nuclear cascade model
// Alain Boudard, CEA-Saclay, France
// Joseph Cugnon, University of Liege, Belgium
// Jean-Christophe David, CEA-Saclay, France
// Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
// Sylvie Leray, CEA-Saclay, France
// Davide Mancusi, CEA-Saclay, France
//
#define INCLXX_IN_GEANT4_MODE 1

#include "globals.hh"

#include "G4INCLCrossSectionsMultiPionsAndResonances.hh"
#include "G4INCLKinematicsUtils.hh"
#include "G4INCLParticleTable.hh"
// #include <cassert>

namespace G4INCL {
  
  template<G4int N>
    struct BystrickyEvaluator {
    static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients<N> const &coeffs) {
      const G4double pMeV = pLab*1E3;
      const G4double ekin=std::sqrt(ParticleTable::effectiveNucleonMass2+pMeV*pMeV)-ParticleTable::effectiveNucleonMass;
      const G4double xrat=ekin*oneOverThreshold;
      const G4double x=std::log(xrat);
      return HornerEvaluator<N>::eval(x, coeffs) * x * std::exp(-0.5*x);
    }
    };
  
  const G4int CrossSectionsMultiPionsAndResonances::nMaxPiNN = 4;
  const G4int CrossSectionsMultiPionsAndResonances::nMaxPiPiN = 4;
  
  const G4double CrossSectionsMultiPionsAndResonances::s11pzOOT = 0.0035761542037692665889;
  const G4double CrossSectionsMultiPionsAndResonances::s01ppOOT = 0.003421025623481919853;
  const G4double CrossSectionsMultiPionsAndResonances::s01pzOOT = 0.0035739814152966403123;
  const G4double CrossSectionsMultiPionsAndResonances::s11pmOOT = 0.0034855350296270480281;
  const G4double CrossSectionsMultiPionsAndResonances::s12pmOOT = 0.0016672224074691565119;
  const G4double CrossSectionsMultiPionsAndResonances::s12ppOOT = 0.0016507643038726931312;
  const G4double CrossSectionsMultiPionsAndResonances::s12zzOOT = 0.0011111111111111111111;
  const G4double CrossSectionsMultiPionsAndResonances::s02pzOOT = 0.00125;
  const G4double CrossSectionsMultiPionsAndResonances::s02pmOOT = 0.0016661112962345883443;
  const G4double CrossSectionsMultiPionsAndResonances::s12mzOOT = 0.0017047391749062392793;
  
  CrossSectionsMultiPionsAndResonances::CrossSectionsMultiPionsAndResonances() :
    s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
    s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
    s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
    s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
    s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
    s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
    s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
    s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
    s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
    s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
  {
  }
  
  G4double CrossSectionsMultiPionsAndResonances::total(Particle const * const p1, Particle const * const p2) {
    G4double inelastic;
    if(p1->isNucleon() && p2->isNucleon()) {
      return CrossSectionsMultiPions::NNTot(p1, p2);
    } else if((p1->isNucleon() && p2->isDelta()) ||
          (p1->isDelta() && p2->isNucleon())) {
      inelastic = CrossSectionsMultiPions::NDeltaToNN(p1, p2);
    } else if((p1->isNucleon() && p2->isPion()) ||
          (p1->isPion() && p2->isNucleon())) {
      return CrossSectionsMultiPions::piNTot(p1,p2);
    } else if((p1->isNucleon() && p2->isEta()) ||
          (p1->isEta() && p2->isNucleon())) {
      inelastic = etaNToPiN(p1,p2) + etaNToPiPiN(p1,p2);
    } else if((p1->isNucleon() && p2->isOmega()) ||
          (p1->isOmega() && p2->isNucleon())) {
      inelastic = omegaNInelastic(p1,p2);
    } else if((p1->isNucleon() && p2->isEtaPrime()) ||
          (p1->isEtaPrime() && p2->isNucleon())) {
      inelastic = etaPrimeNToPiN(p1,p2);
    } else {
      inelastic = 0.;
    }
    
    return inelastic + elastic(p1, p2);
  } 

  
  G4double CrossSectionsMultiPionsAndResonances::elastic(Particle const * const p1, Particle const * const p2) {
    if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta())){
      return CrossSectionsMultiPions::elastic(p1, p2);
    }
    else if ((p1->isNucleon() && p2->isPion()) || (p2->isNucleon() && p1->isPion())){
      return CrossSectionsMultiPions::elastic(p1, p2);
    }
    else if ((p1->isNucleon() && p2->isEta()) || (p2->isNucleon() && p1->isEta())){
      return etaNElastic(p1, p2);
    }
    else if ((p1->isNucleon() && p2->isOmega()) || (p2->isNucleon() && p1->isOmega())){
      return omegaNElastic(p1, p2);
    }
    else {
      return 0.0;
    }
  }
  
  
  G4double CrossSectionsMultiPionsAndResonances::piNToxPiN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
      //
      //     pion-Nucleon producing xpi pions cross sections (corrected due to eta and omega)
      //
// assert(xpi>1 && xpi<=nMaxPiPiN);
// assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));

      const G4double oldXS2Pi=CrossSectionsMultiPions::piNToxPiN(2,particle1, particle2);
      const G4double oldXS3Pi=CrossSectionsMultiPions::piNToxPiN(3,particle1, particle2);
      const G4double oldXS4Pi=CrossSectionsMultiPions::piNToxPiN(4,particle1, particle2);
      const G4double xsEta=piNToEtaN(particle1, particle2);
      const G4double xsOmega=piNToOmegaN(particle1, particle2);
      G4double newXS2Pi=0.; 
      G4double newXS3Pi=0.; 
      G4double newXS4Pi=0.;
      
      if (xpi == 2) {
        if (oldXS4Pi != 0.)
          newXS2Pi=oldXS2Pi;
        else if (oldXS3Pi != 0.) {
          newXS3Pi=oldXS3Pi-xsEta-xsOmega;
          if (newXS3Pi < 1.e-09)
            newXS2Pi=oldXS2Pi-(xsEta+xsOmega-oldXS3Pi);
          else
            newXS2Pi=oldXS2Pi;
        }
        else { 
          newXS2Pi=oldXS2Pi-xsEta-xsOmega;
          if (newXS2Pi < 1.e-09)
            newXS2Pi=0.;
        }
        return newXS2Pi;
      }                     
      else if (xpi == 3) {
        if (oldXS4Pi != 0.) {
          newXS4Pi=oldXS4Pi-xsEta-xsOmega;
          if (newXS4Pi < 1.e-09)
            newXS3Pi=oldXS3Pi-(xsEta+xsOmega-oldXS4Pi);
          else
            newXS3Pi=oldXS3Pi;
        }
        else { 
          newXS3Pi=oldXS3Pi-xsEta-xsOmega;
          if (newXS3Pi < 1.e-09)
            newXS3Pi=0.;
        }
        return newXS3Pi;
      }                 
      else if (xpi == 4) {
        newXS4Pi=oldXS4Pi-xsEta-xsOmega;
        if (newXS4Pi < 1.e-09)
          newXS4Pi=0.;
        return newXS4Pi;
      }     
      else // should never reach this point
        return 0.;
    }
  
  G4double CrossSectionsMultiPionsAndResonances::piNToEtaN(Particle const * const particle1, Particle const * const particle2) {
    //
    //     Pion-Nucleon producing Eta cross sections
    //
// assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
    
    G4double sigma;
    sigma=piMinuspToEtaN(particle1,particle2);
    
    const G4int isoin = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
    
    if (isoin == -1) {
      if ((particle1->getType()) == Proton || (particle2->getType()) == Proton) return sigma;
      else return 0.5 * sigma;
    }
    else if (isoin == 1) {
      if ((particle1->getType()) == Neutron || (particle2->getType()) == Neutron) return sigma;
      else return 0.5 * sigma;
    }
    else return 0. ; // should never return 0. (?) // pi+ p and pi- n return 0.
    
//    return sigma;
  }
  
  G4double CrossSectionsMultiPionsAndResonances::piNToOmegaN(Particle const * const particle1, Particle const * const particle2) {
    //
    //     Pion-Nucleon producing Omega cross sections
    //
// assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
    
    G4double sigma;
    sigma=piMinuspToOmegaN(particle1,particle2);
    
    const G4int isoin = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
    
  if (isoin == -1) {
      if ((particle1->getType()) == Proton || (particle2->getType()) == Proton) return sigma;
      else return 0.5 * sigma;
    }
    else if (isoin == 1) {
      if ((particle1->getType()) == Neutron || (particle2->getType()) == Neutron) return sigma;
      else return 0.5 * sigma;
    }
    else return 0. ; // should never return 0. (?) // pi+ p and pi- n return 0.
    
//    return sigma;
  }
  
#if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
    G4double CrossSectionsMultiPionsAndResonances::piNToEtaPrimeN(Particle const * const /*particle1*/, Particle const * const /*particle2*/) {
#else
  G4double CrossSectionsMultiPionsAndResonances::piNToEtaPrimeN(Particle const * const particle1, Particle const * const particle2) {
#endif
    //
    //     Pion-Nucleon producing EtaPrime cross sections
    //
// assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
    
    return 0.;
  }
  
  G4double CrossSectionsMultiPionsAndResonances::etaNToPiN(Particle const * const particle1, Particle const * const particle2) {
 //
 //     Eta-Nucleon producing Pion cross sections
 //
// assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
 
 const Particle *eta;
 const Particle *nucleon;
 
 if (particle1->isEta()) {
  eta = particle1;
  nucleon = particle2;
 }
 else {
  eta = particle2;
  nucleon = particle1;
 }
 
 const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
 G4double sigma=0.;   
 
 if (pLab <= 574.) 
  sigma= 1.511147E-13*std::pow(pLab,6)- 3.603636E-10*std::pow(pLab,5)+ 3.443487E-07*std::pow(pLab,4)- 1.681980E-04*std::pow(pLab,3)+ 4.437913E-02*std::pow(pLab,2)- 6.172108E+00*pLab+ 4.031449E+02;
 else if (pLab <= 850.) 
  sigma= -8.00018E-14*std::pow(pLab,6)+ 3.50041E-10*std::pow(pLab,5)- 6.33891E-07*std::pow(pLab,4)+ 6.07658E-04*std::pow(pLab,3)- 3.24936E-01*std::pow(pLab,2)+ 9.18098E+01*pLab- 1.06943E+04;
 else if (pLab <= 1300.)
  sigma= 6.56364E-09*std::pow(pLab,3)- 2.07653E-05*std::pow(pLab,2)+ 1.84148E-02*pLab- 1.70427E+00;  
 else {
  G4double ECM=KinematicsUtils::totalEnergyInCM(eta, nucleon);
    G4double massPiZero=ParticleTable::getINCLMass(PiZero);
    G4double massPiMinus=ParticleTable::getINCLMass(PiMinus);
    G4double massProton=ParticleTable::getINCLMass(Proton);
    G4double masseta;
    G4double massnucleon;
    G4double pCM_eta;
  masseta=eta->getMass();
  massnucleon=nucleon->getMass();
  pCM_eta=KinematicsUtils::momentumInCM(ECM, masseta, massnucleon);    
  G4double pCM_PiZero=KinematicsUtils::momentumInCM(ECM, massPiZero, massProton);
  G4double pCM_PiMinus=KinematicsUtils::momentumInCM(ECM, massPiMinus, massProton); // = pCM_PiPlus (because massPiMinus = massPiPlus) 
  sigma = (piMinuspToEtaN(ECM)/2.) * std::pow((pCM_PiZero/pCM_eta), 2) + piMinuspToEtaN(ECM) * std::pow((pCM_PiMinus/pCM_eta), 2);      
 } 
 if (sigma < 0.) sigma=0.;
 return sigma;
}

  G4double CrossSectionsMultiPionsAndResonances::etaNToPiPiN(Particle const * const particle1, Particle const * const particle2) {
      //
      //     Eta-Nucleon producing Two Pions cross sections
      //
// assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
            
      G4double sigma=0.;    
      
      const Particle *eta;
      const Particle *nucleon;

      if (particle1->isEta()) {
        eta = particle1;
        nucleon = particle2;
      }
     else {
        eta = particle2;
        nucleon = particle1;
      }
        
     const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
      
   if (pLab < 450.) 
        sigma = 2.01854221E-13*std::pow(pLab,6) - 3.49750459E-10*std::pow(pLab,5) + 2.46011585E-07*std::pow(pLab,4) - 9.01422901E-05*std::pow(pLab,3) + 1.83382964E-02*std::pow(pLab,2) - 2.03113098E+00*pLab + 1.10358550E+02;       
      else if (pLab < 600.) 
        sigma = 2.01854221E-13*std::pow(450.,6) - 3.49750459E-10*std::pow(450.,5) + 2.46011585E-07*std::pow(450.,4) - 9.01422901E-05*std::pow(450.,3) + 1.83382964E-02*std::pow(450.,2) - 2.03113098E+00*450. + 1.10358550E+02;       
      else if (pLab <= 1300.) 
        sigma = -6.32793049e-16*std::pow(pLab,6) + 3.95985900e-12*std::pow(pLab,5) - 1.01727714e-8*std::pow(pLab,4) + 
                 1.37055547e-05*std::pow(pLab,3) - 1.01830486e-02*std::pow(pLab,2) + 3.93492126*pLab - 609.447145;
   else 
        sigma = etaNToPiN(particle1,particle2);
      
      if (sigma < 0.) sigma = 0.;
      return sigma; // Parameterization from the ANL-Osaka DCC model [PRC88(2013)035209] - eta p --> "pi+pi0 n" + "pi0 pi0 p" total XS
    }
    
    
  G4double CrossSectionsMultiPionsAndResonances::etaNElastic(Particle const * const particle1, Particle const * const particle2) {
      //
      //     Eta-Nucleon elastic cross sections
      //
// assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
      
      G4double sigma=0.;    
      
      const Particle *eta;
      const Particle *nucleon;
      
      if (particle1->isEta()) {
        eta = particle1;
        nucleon = particle2;
      }
     else {
        eta = particle2;
        nucleon = particle1;
      }
      
     const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
      
   if (pLab < 700.) 
        sigma = 3.6838e-15*std::pow(pLab,6) - 9.7815e-12*std::pow(pLab,5) + 9.7914e-9*std::pow(pLab,4) - 
                4.3222e-06*std::pow(pLab,3) + 7.9188e-04*std::pow(pLab,2) - 1.8379e-01*pLab + 84.965;     
      else if (pLab < 1400.) 
        sigma = 3.562630e-16*std::pow(pLab,6) - 2.384766e-12*std::pow(pLab,5) + 6.601312e-9*std::pow(pLab,4) - 
                9.667078e-06*std::pow(pLab,3) + 7.894845e-03*std::pow(pLab,2) - 3.409200*pLab + 609.8501;
      else if (pLab < 2025.)
        sigma = -1.041950E-03*pLab + 2.110529E+00;
      else  
       sigma=0.;
        
      if (sigma < 0.) sigma = 0.;
      return sigma; // Parameterization from the ANL-Osaka DCC model [PRC88(2013)035209]
    }

  G4double CrossSectionsMultiPionsAndResonances::omegaNInelastic(Particle const * const particle1, Particle const * const particle2) {
    //
    //     Omega-Nucleon inelastic cross sections
    //
// assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
    
    G4double sigma=0.;    
    
    const Particle *omega;
    const Particle *nucleon;
    
    if (particle1->isOmega()) {
      omega = particle1;
      nucleon = particle2;
      }
    else {
      omega = particle2;
      nucleon = particle1;
    }
    
    const G4double pLab = KinematicsUtils::momentumInLab(omega, nucleon)/1000.; // GeV/c
    
    sigma = 20. + 4.0/pLab; // Eq.(24) in G.I. Lykasov et al., EPJA 6, 71-81 (1999)
    
    return sigma;
  }

  
  G4double CrossSectionsMultiPionsAndResonances::omegaNElastic(Particle const * const particle1, Particle const * const particle2) {
    //
    //     Omega-Nucleon elastic cross sections
    //
// assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
    
    G4double sigma=0.;    
    
    const Particle *omega;
    const Particle *nucleon;
    
    if (particle1->isOmega()) {
      omega = particle1;
      nucleon = particle2;
    }
    else {
      omega = particle2;
      nucleon = particle1;
    }
    
    const G4double pLab = KinematicsUtils::momentumInLab(omega, nucleon)/1000.; // GeV/c
    
    sigma = 5.4 + 10.*std::exp(-0.6*pLab);  // Eq.(21) in G.I. Lykasov et al., EPJA 6, 71-81 (1999)
    
    return sigma;
  }
  
    
  G4double CrossSectionsMultiPionsAndResonances::omegaNToPiN(Particle const * const particle1, Particle const * const particle2) {
    //
    //     Omega-Nucleon producing Pion cross sections
    //
// assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
    
    G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
    
    G4double massPiZero=ParticleTable::getINCLMass(PiZero);
    G4double massPiMinus=ParticleTable::getINCLMass(PiMinus);
    G4double massProton=ParticleTable::getINCLMass(Proton);
    
    G4double massomega;
    G4double massnucleon;
    G4double pCM_omega;
    G4double pLab_omega;
    
    G4double sigma=0.;    
    
    if (particle1->isOmega()) {
      massomega=particle1->getMass();
      massnucleon=particle2->getMass();
    }
    else {
      massomega=particle2->getMass();
      massnucleon=particle1->getMass();
    }
    pCM_omega=KinematicsUtils::momentumInCM(ECM, massomega, massnucleon);   
    pLab_omega=KinematicsUtils::momentumInLab(ECM*ECM, massomega, massnucleon);   
    
    G4double pCM_PiZero=KinematicsUtils::momentumInCM(ECM, massPiZero, massProton);
    G4double pCM_PiMinus=KinematicsUtils::momentumInCM(ECM, massPiMinus, massProton); // = pCM_PiPlus (because massPiMinus = massPiPlus)
    
    sigma = (piMinuspToOmegaN(ECM)/2.) * std::pow((pCM_PiZero/pCM_omega), 2) 
    + piMinuspToOmegaN(ECM) * std::pow((pCM_PiMinus/pCM_omega), 2);     
  
  if (sigma > omegaNInelastic(particle1, particle2) || (pLab_omega < 200.)) {
//  if (sigma > omegaNInelastic(particle1, particle2)) {
   sigma = omegaNInelastic(particle1, particle2);
  }  
  
    return sigma;
  }
  
  
  G4double CrossSectionsMultiPionsAndResonances::omegaNToPiPiN(Particle const * const particle1, Particle const * const particle2) {
    //
    //     Omega-Nucleon producing 2 PionS cross sections
    //
// assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
    
    G4double sigma=0.;    
    
    sigma = omegaNInelastic(particle1,particle2) - omegaNToPiN(particle1,particle2) ;
    
    return sigma;
  }
  
  
#if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
  G4double CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(Particle const * const /*particle1*/, Particle const * const /*particle2*/) {
#else
  G4double CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(Particle const * const particle1, Particle const * const particle2) {
#endif
    //
    //     EtaPrime-Nucleon producing Pion cross sections
    //
// assert((particle1->isNucleon() && particle2->isEtaPrime()) || (particle1->isEtaPrime() && particle2->isNucleon()));
    
    return 0.;
  } 
  
  G4double CrossSectionsMultiPionsAndResonances::piMinuspToEtaN(Particle const * const particle1, Particle const * const particle2) {
    //
    //     Pion-Nucleon producing Eta cross sections
    //
// assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
    
    G4double masspion;
    G4double massnucleon;
    if (particle1->isPion()) {
        masspion=particle1->getMass();
        massnucleon=particle2->getMass();
    }
    else {
        masspion=particle2->getMass();
        massnucleon=particle1->getMass();
    }
    
    G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
    G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.; // GeV/c
    
    G4double sigma;
    
// new parameterization (JCD) - end of february 2016
    if (ECM < 1486.5) sigma=0.;
    else
    {
      if (ECM < 1535.)
      {
        sigma = -0.0000003689197974814*std::pow(ECM,4) + 0.002260193900097*std::pow(ECM,3) - 5.193105877187*std::pow(ECM,2) + 5303.505273919*ECM - 2031265.900648;
      }
      else if (ECM < 1670.)
      {
        sigma = -0.0000000337986446*std::pow(ECM,4) + 0.000218279989*std::pow(ECM,3) - 0.528276144*std::pow(ECM,2) + 567.828367*ECM - 228709.42;
      }
      else if (ECM < 1714.)
      {
        sigma =  0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
      }
      else sigma=1.47*std::pow(plab, -1.68);
    }   
//    

    return sigma;
  }
  
  G4double CrossSectionsMultiPionsAndResonances::piMinuspToEtaN(const G4double ECM) {
    //
    //     Pion-Nucleon producing Eta cross sections
    //
    
  const G4double masspion = ParticleTable::getRealMass(PiMinus);    
  const G4double massnucleon = ParticleTable::getRealMass(Proton);    
  
  G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.;  // GeV/c
    
    G4double sigma;
    
// new parameterization (JCD) - end of february 2016
    if (ECM < 1486.5) sigma=0.;
    else
    {
      if (ECM < 1535.)
      {
        sigma = -0.0000003689197974814*std::pow(ECM,4) + 0.002260193900097*std::pow(ECM,3) - 5.193105877187*std::pow(ECM,2) + 5303.505273919*ECM - 2031265.900648;
      }
      else if (ECM < 1670.)
      {
        sigma = -0.0000000337986446*std::pow(ECM,4) + 0.000218279989*std::pow(ECM,3) - 0.528276144*std::pow(ECM,2) + 567.828367*ECM - 228709.42;
      }
      else if (ECM < 1714.)
      {
        sigma =  0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
      }
      else sigma=1.47*std::pow(plab, -1.68);
    }   
//    

    return sigma;
  }

  G4double CrossSectionsMultiPionsAndResonances::piMinuspToOmegaN(Particle const * const particle1, Particle const * const particle2) {
    //
    //     Pion-Nucleon producing Omega cross sections
    //
// assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
//jcd to be removed
//  return 0.;
//jcd    
    
//    G4double param=1.095 ; // Deneye (Thesis)
    G4double param=1.0903 ; // JCD (threshold taken into account)
  
  G4double masspion;
    G4double massnucleon;
    if (particle1->isPion()) {
        masspion=particle1->getMass();
        massnucleon=particle2->getMass();
    }
    else {
        masspion=particle2->getMass();
        massnucleon=particle1->getMass();
    }
    G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
    G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.;  // GeV/c
    
    G4double sigma;
    if (plab < param) sigma=0.;
    else sigma=13.76*(plab-param)/(std::pow(plab, 3.33) - 1.07); // Phys. Rev. C 41, 1701–1718 (1990)
  
    return sigma;
}
      G4double CrossSectionsMultiPionsAndResonances::piMinuspToOmegaN(const G4double ECM) {
        //
        //     Pion-Nucleon producing Omega cross sections
        //
//jcd to be removed
//    return 0.;
//jcd    
        
//    G4double param=1.095 ; // Deneye (Thesis)
    G4double param=1.0903 ; // JCD (threshold taken into account)

        const G4double masspion = ParticleTable::getRealMass(PiMinus);    
        const G4double massnucleon = ParticleTable::getRealMass(Proton);    

        G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.;  // GeV/c
        
        G4double sigma;
        if (plab < param) sigma=0.;
        else sigma=13.76*(plab-param)/(std::pow(plab, 3.33)-1.07);
        
        return sigma;
}

    G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaIso(const G4double ener, const G4int iso) {
    
    const G4double Ecm=0.001*ener;
    G4double sNNEta; // pp->pp+eta(+X)
    G4double sNNEta1; // np->np+eta(+X)
    G4double sNNEta2; // np->d+eta (d will be considered as np - How far is this right?)
    G4double x=Ecm*Ecm/5.88;
                
//jcd
  if (Ecm >= 3.05) {
     sNNEta = 2.5*std::pow((x-1.),1.47)*std::pow(x,-1.25)*1000.;
  }      
  else if(Ecm >= 2.6) {
     sNNEta = -327.29*Ecm*Ecm*Ecm + 2870.*Ecm*Ecm - 7229.3*Ecm + 5273.3;
     if (sNNEta <= NNToNNEtaExcluIso(ener, 2)*1000.) sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
  }
  else {
     sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
  }
//jcd
  if (sNNEta < 1.e-9) sNNEta = 0.;

  if (iso != 0) {
     return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
  }

  if(Ecm >= 6.25) {
     sNNEta1 = sNNEta;
  }
  else if (Ecm >= 2.6) {
     sNNEta1 = sNNEta*std::exp(-(-5.53151576/Ecm+0.8850425));
  }
  else if (Ecm >= 2.525) { // = exclusive pn
      sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ecm + 56581.54*Ecm*Ecm*Ecm - 270212.6*Ecm*Ecm + 571650.6*Ecm - 451091.6;
  }
  else { // = exclusive pn
      sNNEta1 = 17570.217219*Ecm*Ecm - 84910.985402*Ecm + 102585.55847;
  }

  sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
  if (sNNEta2 < 0.) sNNEta2=0.;

  sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;

  G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
  G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
  G4double Meta=ParticleTable::getRealMass(Eta)/1000.;
    if (sNNEta < 1.e-9 || Ecm < Mn+Mp+Meta) sNNEta = 0.;
    
    return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
  }

  
  G4double CrossSectionsMultiPionsAndResonances::NNToNNEta(Particle const * const particle1, Particle const * const particle2) {
    
  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
    
    if (iso != 0) {
   return NNToNNEtaIso(ener, iso);
    }
    else {
      return 0.5*(NNToNNEtaIso(ener, 0)+NNToNNEtaIso(ener, 2));
    }
    }
  
  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaExcluIso(const G4double ener, const G4int iso) {
        
  const G4double Ecm=0.001*ener;
  G4double sNNEta; // pp->pp+eta
  G4double sNNEta1; // np->np+eta
  G4double sNNEta2; // np->d+eta (d wil be considered as np - How far is this right?)
        
  if(Ecm>=3.875) { // By hand (JCD)
     sNNEta = -13.008*Ecm*Ecm + 84.531*Ecm + 36.234;
  }
  else if(Ecm>=2.725) { // By hand (JCD)
     sNNEta = -913.2809*std::pow(Ecm,5) + 15564.27*std::pow(Ecm,4) - 105054.9*std::pow(Ecm,3) + 351294.2*std::pow(Ecm,2) - 582413.9*Ecm + 383474.7;
  }
  else if(Ecm>=2.575) { // By hand (JCD)
     sNNEta = -2640.3*Ecm*Ecm + 14692*Ecm - 20225;
  }
  else {
     sNNEta = -147043.497285*std::pow(Ecm,4) + 1487222.5438123*std::pow(Ecm,3) - 5634399.900744*std::pow(Ecm,2) + 9477290.199378*Ecm - 5972174.353438;
  }
        
   G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
   G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
   G4double Meta=ParticleTable::getRealMass(Eta)/1000.;
   G4double Thr0=0.;
   if (iso > 0) {
    Thr0=2.*Mp+Meta;
   }
   else if (iso < 0) {
    Thr0=2.*Mn+Meta;
   }
   else {
    Thr0=Mn+Mp+Meta;
   }

   if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.;  // Thr0: Ecm threshold
        
   if (iso != 0) {
      return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
   }
        
   if(Ecm>=3.9) {
      sNNEta1 = sNNEta;
   }
   else if (Ecm >= 3.5) {
      sNNEta1 = -1916.2*Ecm*Ecm*Ecm + 21556.0*Ecm*Ecm - 80828.0*Ecm + 101200.0;
   }
   else if (Ecm >= 2.525) {
      sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ecm + 56581.54*Ecm*Ecm*Ecm - 270212.6*Ecm*Ecm + 571650.6*Ecm - 451091.6;
   }
   else {
      sNNEta1 = 17570.217219*Ecm*Ecm - 84910.985402*Ecm + 102585.55847;
  }
    
   sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
   if (sNNEta2 < 0.) sNNEta2=0.;
    
   sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
   if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.;  // Thr0: Ecm threshold
        
  return sNNEta/1000.; // parameterization in microbarn (not millibarn)!

}
      
  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaExclu(Particle const * const particle1, Particle const * const particle2) {
        
   const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
   const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
        
   if (iso != 0) {
      return NNToNNEtaExcluIso(ener, iso);
   }
   else {
      return 0.5*(NNToNNEtaExcluIso(ener, 0)+NNToNNEtaExcluIso(ener, 2));
   }
}

   
   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaIso(const G4double ener, const G4int iso) {
    
    const G4double Ecm=0.001*ener;
    G4double sNNOmega; // pp->pp+eta(+X)
    G4double sNNOmega1; // np->np+eta(+X)
    G4double x=Ecm*Ecm/7.06;
    
    if(Ecm>4.0) {
     sNNOmega = 2.5*std::pow(x-1, 1.47)*std::pow(x, -1.11) ;
    }
    else if(Ecm>2.802) { // 2802 MeV -> threshold to open inclusive (based on multipion threshold and omega mass) 
     sNNOmega = (568.5254*Ecm*Ecm - 2694.045*Ecm + 3106.247)/1000.;
     if (sNNOmega <= NNToNNOmegaExcluIso(ener, 2)) sNNOmega = NNToNNOmegaExcluIso(ener, 2);
    }
    else {
     sNNOmega = NNToNNOmegaExcluIso(ener, 2);
    }
    
    if (sNNOmega < 1.e-9) sNNOmega = 0.;
    
    if (iso != 0) {
     return sNNOmega; 
    }
    
    sNNOmega1 = 3.*sNNOmega; // 3.0: ratio pn/pp (5 from theory; 2 from experiments)
        
    sNNOmega = 2.*sNNOmega1-sNNOmega;
    
    if (sNNOmega < 1.e-9) sNNOmega = 0.;
    
    return sNNOmega; 
   }
   
   
   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmega(Particle const * const particle1, Particle const * const particle2) {
    
    const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
    const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
//jcd to be removed
//    return 0.;
//jcd    
    if (iso != 0) {
     return NNToNNOmegaIso(ener, iso);
    }
    else {
     return 0.5*(NNToNNOmegaIso(ener, 0)+NNToNNOmegaIso(ener, 2));
    }
   }
   
   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaExcluIso(const G4double ener, const G4int iso) {
        
        const G4double Ecm=0.001*ener;
        G4double sNNOmega; // pp->pp+eta
        G4double sNNOmega1; // np->np+eta
    G4double sthroot=std::sqrt(7.06);
        
        if(Ecm>=3.0744) { // By hand (JCD)
          sNNOmega = 330.*(Ecm-sthroot)/(1.05+std::pow((Ecm-sthroot),2));
        }
        else if(Ecm>=2.65854){
          sNNOmega = -1208.09757*std::pow(Ecm,3) + 10773.3322*std::pow(Ecm,2) - 31661.0223*Ecm + 30728.7241 ;
        }
        else {
          sNNOmega = 0. ;
        }
        
    G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
    G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
    G4double Momega=ParticleTable::getRealMass(Omega)/1000.;
    G4double Thr0=0.;
    if (iso > 0) {
     Thr0=2.*Mp+Momega;
    }
    else if (iso < 0) {
     Thr0=2.*Mn+Momega;
    }
    else {
     Thr0=Mn+Mp+Momega;
    }
    
    if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.;  // Thr0: Ecm threshold
        
        if (iso != 0) {
          return sNNOmega/1000.; // parameterization in microbarn (not millibarn)!
        }
        
    sNNOmega1 = 3*sNNOmega; // 3.0: ratio pn/pp
        
        sNNOmega = 2*sNNOmega1-sNNOmega;
        if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.;
        
        return sNNOmega/1000.; // parameterization in microbarn (not millibarn)!
      }
      
   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaExclu(Particle const * const particle1, Particle const * const particle2) {
//jcd to be removed
//    return 0.;
//jcd    
        
        const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
        const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
        
        if (iso != 0) {
          return NNToNNOmegaExcluIso(ener, iso);
        }
        else {
          return 0.5*(NNToNNOmegaExcluIso(ener, 0)+NNToNNOmegaExcluIso(ener, 2));
        }
      }
   
   
    G4double CrossSectionsMultiPionsAndResonances::NNToxPiNN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
    //
    //     Nucleon-Nucleon producing xpi pions cross sections
    //
// assert(xpi>0 && xpi<=nMaxPiNN);
// assert(particle1->isNucleon() && particle2->isNucleon());
    
    G4double oldXS1Pi=CrossSectionsMultiPions::NNToxPiNN(1,particle1, particle2);
    G4double oldXS2Pi=CrossSectionsMultiPions::NNToxPiNN(2,particle1, particle2);
    G4double oldXS3Pi=CrossSectionsMultiPions::NNToxPiNN(3,particle1, particle2);
    G4double oldXS4Pi=CrossSectionsMultiPions::NNToxPiNN(4,particle1, particle2);
    G4double xsEtaOmega=NNToNNEta(particle1, particle2)+NNToNNOmega(particle1, particle2);
    G4double newXS1Pi=0.; 
    G4double newXS2Pi=0.; 
    G4double newXS3Pi=0.; 
    G4double newXS4Pi=0.; 
      
    if (xpi == 1) {
      if (oldXS4Pi != 0. || oldXS3Pi != 0.)
        newXS1Pi=oldXS1Pi;
      else if (oldXS2Pi != 0.) {
        newXS2Pi=oldXS2Pi-xsEtaOmega;
       if (newXS2Pi < 0.)
         newXS1Pi=oldXS1Pi-(xsEtaOmega-oldXS2Pi);
        else
          newXS1Pi=oldXS1Pi;
      }
      else 
        newXS1Pi=oldXS1Pi-xsEtaOmega;
     return newXS1Pi;
    }
    else if (xpi == 2) {
      if (oldXS4Pi != 0.)
        newXS2Pi=oldXS2Pi;
      else if (oldXS3Pi != 0.) {
        newXS3Pi=oldXS3Pi-xsEtaOmega;
        if (newXS3Pi < 0.)
          newXS2Pi=oldXS2Pi-(xsEtaOmega-oldXS3Pi);
        else
          newXS2Pi=oldXS2Pi;
      }
      else { 
        newXS2Pi=oldXS2Pi-xsEtaOmega;
        if (newXS2Pi < 0.)
          newXS2Pi=0.;
      }
      return newXS2Pi;
    }
    else if (xpi == 3) {
      if (oldXS4Pi != 0.) {
        newXS4Pi=oldXS4Pi-xsEtaOmega;
       if (newXS4Pi < 0.)
        newXS3Pi=oldXS3Pi-(xsEtaOmega-oldXS4Pi);
       else
        newXS3Pi=oldXS3Pi;
      }
      else { 
        newXS3Pi=oldXS3Pi-xsEtaOmega;       
        if (newXS3Pi < 0.)
          newXS3Pi=0.;
      }
      return newXS3Pi;
    }
    else if (xpi == 4) {
      newXS4Pi=oldXS4Pi-xsEtaOmega;
      if (newXS4Pi < 0.)
        newXS4Pi=0.;
    return newXS4Pi;
    }

    else // should never reach this point
      return 0.;
  }


      G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaOnePi(Particle const * const particle1, Particle const * const particle2) {
        // Cross section for nucleon-nucleon producing one eta and one pion
        
        const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
        if (iso!=0)
          return 0.;
        
        const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta (= 2705.55 - 2018.563; 4074595.287720512986=2018.563*2018.563)
        if (ener < 2018.563) return 0.;
        
        const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
        const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);

        return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
      }
      
      G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
        const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
        if (ener < 2018.563) return 0.;
        const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
        
        const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
        if (iso != 0)
          return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
        else {
          const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
          return 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
        }
      }
      
      G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaTwoPi(Particle const * const particle1, Particle const * const particle2) {
        //
        //     Nucleon-Nucleon producing one eta and two pions 
        //
        const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
        if (ener < 2018.563) return 0.;
        const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
        
        
        const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
        if (iso != 0) {
          return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
        }
        else {
          const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
          return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
        }
      }
      
      G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaThreePi(Particle const * const particle1, Particle const * const particle2) {
        //
        //     Nucleon-Nucleon producing one eta and three pions
        //
        
        const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
        if (ener < 2018.563) return 0.;
        const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
        
        
        const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
        const G4double xs1pi2=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2);
        const G4double xs2pi2=CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
        if (iso != 0)
          return CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
        else {
          const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
          const G4double xs1pi0=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0);
          const G4double xs2pi0=CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0);
          return 0.5*(CrossSectionsMultiPions::NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
        }
      }
      
      G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaFourPi(Particle const * const particle1, Particle const * const particle2) {
        //
        //     Nucleon-Nucleon producing one eta and four pions
        //
        
        const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
        if (ener < 2018.563) return 0.;
        const G4double s = ener*ener;
        const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
    G4double xsinelas;
        if (i!=0) 
     xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
    else 
          xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
        if (xsinelas <= 1.e-9) return 0.;
       G4double ratio=(NNToNNEta(particle1, particle2)-NNToNNEtaExclu(particle1, particle2))/xsinelas;
        if(s<6.25E6)
          return 0.;
        const G4double sigma = NNToNNEta(particle1, particle2) - NNToNNEtaExclu(particle1, particle2) - ratio*(NNToNNEtaOnePiOrDelta(particle1, particle2) + NNToNNEtaTwoPi(particle1, particle2) + NNToNNEtaThreePi(particle1, particle2));
        return ((sigma>1.e-9) ? sigma : 0.);
      }
      
      G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaxPi(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
        //
        //     Nucleon-Nucleon producing one eta and xpi pions
        //
// assert(xpi>0 && xpi<=nMaxPiNN);
// assert(particle1->isNucleon() && particle2->isNucleon());
        
        const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
        if (ener < 2018.563) return 0.;
        const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
    G4double xsinelas;
        if (i!=0) 
     xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
    else 
          xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
        if (xsinelas <= 1.e-9) return 0.;
       G4double ratio=(NNToNNEta(particle1, particle2)-NNToNNEtaExclu(particle1, particle2))/xsinelas;

        if (xpi == 1)
          return NNToNNEtaOnePi(particle1, particle2)*ratio;
        else if (xpi == 2)
          return NNToNNEtaTwoPi(particle1, particle2)*ratio;
        else if (xpi == 3)
          return NNToNNEtaThreePi(particle1, particle2)*ratio;
        else if (xpi == 4)
          return NNToNNEtaFourPi(particle1, particle2);
        else // should never reach this point
          return 0.;
      }

      
      G4double CrossSectionsMultiPionsAndResonances::NNToNDeltaEta(Particle const * const p1, Particle const * const p2) {
// assert(p1->isNucleon() && p2->isNucleon());
    const G4int i = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
        const G4double ener=KinematicsUtils::totalEnergyInCM(p1, p2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
        if (ener < 2018.563) return 0.;
    G4double xsinelas;
        if (i!=0) 
     xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
    else 
          xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
        if (xsinelas <= 1.e-9) return 0.;
       G4double ratio=(NNToNNEta(p1, p2)-NNToNNEtaExclu(p1, p2))/xsinelas;
    G4double sigma = NNToNNEtaOnePiOrDelta(p1, p2)*ratio;
    if(i==0)
          sigma *= 0.5;
    return sigma;
      }

  
  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaOnePi(Particle const * const particle1, Particle const * const particle2) {
   // Cross section for nucleon-nucleon producing one omega and one pion
   
   const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
   if (iso!=0)
    return 0.;
   
   const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega (= 2802. - 2018.563; 4074595.287720512986=2018.563*2018.563)
   if (ener < 2018.563) return 0.;
   
   const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
   const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
   
   return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
  }
  
  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
   const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
   if (ener < 2018.563) return 0.;
   const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
   
   const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
   if (iso != 0)
    return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
   else {
    const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
    return 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
   }
  }
  
  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaTwoPi(Particle const * const particle1, Particle const * const particle2) {
   //
   //     Nucleon-Nucleon producing one omega and two pions 
   //
   const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
   if (ener < 2018.563) return 0.;
   const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
   
   
   const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
   if (iso != 0) {
    return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
   }
   else {
    const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
    return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
   }
  }
  
  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaThreePi(Particle const * const particle1, Particle const * const particle2) {
   //
   //     Nucleon-Nucleon producing one omega and three pions
   //
   
   const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
   if (ener < 2018.563) return 0.;
   const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
   
   
   const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
   const G4double xs1pi2=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2);
   const G4double xs2pi2=CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
   if (iso != 0)
    return CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
   else {
    const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
    const G4double xs1pi0=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0);
    const G4double xs2pi0=CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0);
    return 0.5*(CrossSectionsMultiPions::NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
   }
  }
  
  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaFourPi(Particle const * const particle1, Particle const * const particle2) {
   //
   //     Nucleon-Nucleon producing one omega and four pions
   //
//jcd to be removed
//   return 0.;
//jcd    
   
   const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
   if (ener < 2018.563) return 0.;
   const G4double s = ener*ener;
   const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
   G4double xsinelas;
   if (i!=0) 
    xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
   else 
    xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
   if (xsinelas <= 1.e-9) return 0.;
   G4double ratio=(NNToNNOmega(particle1, particle2)-NNToNNOmegaExclu(particle1, particle2))/xsinelas;
   if(s<6.25E6)
    return 0.;
   const G4double sigma = NNToNNOmega(particle1, particle2) - NNToNNOmegaExclu(particle1, particle2) - ratio*(NNToNNOmegaOnePiOrDelta(particle1, particle2) + NNToNNOmegaTwoPi(particle1, particle2) + NNToNNOmegaThreePi(particle1, particle2));
   return ((sigma>1.e-9) ? sigma : 0.);
  }
  
  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaxPi(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
   //
   //     Nucleon-Nucleon producing one omega and xpi pions
   //
// assert(xpi>0 && xpi<=nMaxPiNN);
// assert(particle1->isNucleon() && particle2->isNucleon());
//jcd to be removed
//   return 0.;
//jcd    
   
   const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
   if (ener < 2018.563) return 0.;
   const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
   G4double xsinelas;
   if (i!=0) 
    xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
   else 
    xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
   if (xsinelas <= 1.e-9) return 0.;
   G4double ratio=(NNToNNOmega(particle1, particle2)-NNToNNOmegaExclu(particle1, particle2))/xsinelas;
   
   if (xpi == 1)
    return NNToNNOmegaOnePi(particle1, particle2)*ratio;
   else if (xpi == 2)
    return NNToNNOmegaTwoPi(particle1, particle2)*ratio;
   else if (xpi == 3)
    return NNToNNOmegaThreePi(particle1, particle2)*ratio;
   else if (xpi == 4)
    return NNToNNOmegaFourPi(particle1, particle2);
   else // should never reach this point
    return 0.;
  }
  
  
  G4double CrossSectionsMultiPionsAndResonances::NNToNDeltaOmega(Particle const * const p1, Particle const * const p2) {
// assert(p1->isNucleon() && p2->isNucleon());
//jcd to be removed
//   return 0.;
//jcd    
   const G4int i = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
   const G4double ener=KinematicsUtils::totalEnergyInCM(p1, p2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
   if (ener < 2018.563) return 0.;
   G4double xsinelas;
   if (i!=0) 
    xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
   else 
    xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
   if (xsinelas <= 1.e-9) return 0.;
   G4double ratio=(NNToNNOmega(p1, p2)-NNToNNOmegaExclu(p1, p2))/xsinelas;
   G4double sigma = NNToNNOmegaOnePiOrDelta(p1, p2)*ratio;
   if(i==0)
    sigma *= 0.5;
   return sigma;
  }
  
  
} // namespace G4INCL