Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLCrossSectionsMultiPionsAndResonances.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/inclxx/incl_physics/src/G4INCLCrossSectionsMultiPionsAndResonances.cc (Version 11.3.0) and /processes/hadronic/models/inclxx/incl_physics/src/G4INCLCrossSectionsMultiPionsAndResonances.cc (Version 10.4.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // INCL++ intra-nuclear cascade model              26 // INCL++ intra-nuclear cascade model
 27 // Alain Boudard, CEA-Saclay, France               27 // Alain Boudard, CEA-Saclay, France
 28 // Joseph Cugnon, University of Liege, Belgium     28 // Joseph Cugnon, University of Liege, Belgium
 29 // Jean-Christophe David, CEA-Saclay, France       29 // Jean-Christophe David, CEA-Saclay, France
 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H     30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
 31 // Sylvie Leray, CEA-Saclay, France                31 // Sylvie Leray, CEA-Saclay, France
 32 // Davide Mancusi, CEA-Saclay, France              32 // Davide Mancusi, CEA-Saclay, France
 33 //                                                 33 //
 34 #define INCLXX_IN_GEANT4_MODE 1                    34 #define INCLXX_IN_GEANT4_MODE 1
 35                                                    35 
 36 #include "globals.hh"                              36 #include "globals.hh"
 37                                                    37 
 38 /** \file G4INCLCrossSectionsMultiPionsAndReso     38 /** \file G4INCLCrossSectionsMultiPionsAndResonances.cc
 39  * \brief Multipion and mesonic Resonances cro     39  * \brief Multipion and mesonic Resonances cross sections
 40  *                                                 40  *
 41  * \date 4th February 2014                         41  * \date 4th February 2014
 42  * \author Jean-Christophe David                   42  * \author Jean-Christophe David
 43  */                                                43  */
 44                                                    44 
 45 #include "G4INCLCrossSectionsMultiPionsAndReso     45 #include "G4INCLCrossSectionsMultiPionsAndResonances.hh"
 46 #include "G4INCLKinematicsUtils.hh"                46 #include "G4INCLKinematicsUtils.hh"
 47 #include "G4INCLParticleTable.hh"                  47 #include "G4INCLParticleTable.hh"
 48 // #include <cassert>                              48 // #include <cassert>
 49                                                    49 
 50 namespace G4INCL {                                 50 namespace G4INCL {
 51                                                <<  51   
 52     template<G4int N>                          <<  52   template<G4int N>
 53     struct BystrickyEvaluator {                    53     struct BystrickyEvaluator {
 54         static G4double eval(const G4double pL <<  54     static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients<N> const &coeffs) {
 55             const G4double pMeV = pLab*1E3;    <<  55       const G4double pMeV = pLab*1E3;
 56             const G4double ekin=std::sqrt(Part <<  56       const G4double ekin=std::sqrt(ParticleTable::effectiveNucleonMass2+pMeV*pMeV)-ParticleTable::effectiveNucleonMass;
 57             const G4double xrat=ekin*oneOverTh <<  57       const G4double xrat=ekin*oneOverThreshold;
 58             const G4double x=std::log(xrat);   <<  58       const G4double x=std::log(xrat);
 59             return HornerEvaluator<N>::eval(x, <<  59       return HornerEvaluator<N>::eval(x, coeffs) * x * std::exp(-0.5*x);
 60         }                                      <<  60     }
 61     };                                             61     };
 62                                                <<  62   
 63     const G4int CrossSectionsMultiPionsAndReso <<  63   const G4int CrossSectionsMultiPionsAndResonances::nMaxPiNN = 4;
 64     const G4int CrossSectionsMultiPionsAndReso <<  64   const G4int CrossSectionsMultiPionsAndResonances::nMaxPiPiN = 4;
 65                                                <<  65   
 66     const G4double CrossSectionsMultiPionsAndR <<  66   const G4double CrossSectionsMultiPionsAndResonances::s11pzOOT = 0.0035761542037692665889;
 67     const G4double CrossSectionsMultiPionsAndR <<  67   const G4double CrossSectionsMultiPionsAndResonances::s01ppOOT = 0.003421025623481919853;
 68     const G4double CrossSectionsMultiPionsAndR <<  68   const G4double CrossSectionsMultiPionsAndResonances::s01pzOOT = 0.0035739814152966403123;
 69     const G4double CrossSectionsMultiPionsAndR <<  69   const G4double CrossSectionsMultiPionsAndResonances::s11pmOOT = 0.0034855350296270480281;
 70     const G4double CrossSectionsMultiPionsAndR <<  70   const G4double CrossSectionsMultiPionsAndResonances::s12pmOOT = 0.0016672224074691565119;
 71     const G4double CrossSectionsMultiPionsAndR <<  71   const G4double CrossSectionsMultiPionsAndResonances::s12ppOOT = 0.0016507643038726931312;
 72     const G4double CrossSectionsMultiPionsAndR <<  72   const G4double CrossSectionsMultiPionsAndResonances::s12zzOOT = 0.0011111111111111111111;
 73     const G4double CrossSectionsMultiPionsAndR <<  73   const G4double CrossSectionsMultiPionsAndResonances::s02pzOOT = 0.00125;
 74     const G4double CrossSectionsMultiPionsAndR <<  74   const G4double CrossSectionsMultiPionsAndResonances::s02pmOOT = 0.0016661112962345883443;
 75     const G4double CrossSectionsMultiPionsAndR <<  75   const G4double CrossSectionsMultiPionsAndResonances::s12mzOOT = 0.0017047391749062392793;
 76                                                <<  76   
 77     CrossSectionsMultiPionsAndResonances::Cros <<  77   CrossSectionsMultiPionsAndResonances::CrossSectionsMultiPionsAndResonances() :
 78     s11pzHC(-2.228000000000294018,8.7560000000     78     s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
 79     s01ppHC(2.0570000000126518344,-6.029000000     79     s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
 80     s01pzHC(0.18030000000000441851,7.870099999     80     s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
 81     s11pmHC(0.20590000000000031866,3.345099999     81     s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
 82     s12pmHC(-0.77235999999999901328,4.26265999     82     s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
 83     s12ppHC(-0.75724999999999975664,2.09343999     83     s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
 84     s12zzHC(-0.89599999999996965072,7.88299999     84     s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
 85     s02pzHC(-1.0579999999999967036,11.11399999     85     s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
 86     s02pmHC(2.4009000000012553286,-7.768000000     86     s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
 87     s12mzHC(-0.21858699999999976269,1.91489999     87     s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
 88     {                                          <<  88   {
 89     }                                          <<  89   }
 90                                                <<  90   
 91     G4double CrossSectionsMultiPionsAndResonan <<  91   G4double CrossSectionsMultiPionsAndResonances::total(Particle const * const p1, Particle const * const p2) {
 92         G4double inelastic;                    <<  92     G4double inelastic;
 93         if(p1->isNucleon() && p2->isNucleon()) <<  93     if(p1->isNucleon() && p2->isNucleon()) {
 94             return CrossSectionsMultiPions::NN <<  94       return CrossSectionsMultiPions::NNTot(p1, p2);
 95         } else if((p1->isNucleon() && p2->isDe <<  95     } else if((p1->isNucleon() && p2->isDelta()) ||
 96                   (p1->isDelta() && p2->isNucl <<  96           (p1->isDelta() && p2->isNucleon())) {
 97             inelastic = CrossSectionsMultiPion <<  97       inelastic = CrossSectionsMultiPions::NDeltaToNN(p1, p2);
 98         } else if((p1->isNucleon() && p2->isPi <<  98     } else if((p1->isNucleon() && p2->isPion()) ||
 99                   (p1->isPion() && p2->isNucle <<  99           (p1->isPion() && p2->isNucleon())) {
100             return CrossSectionsMultiPions::pi << 100       return CrossSectionsMultiPions::piNTot(p1,p2);
101         } else if((p1->isNucleon() && p2->isEt << 101     } else if((p1->isNucleon() && p2->isEta()) ||
102                   (p1->isEta() && p2->isNucleo << 102           (p1->isEta() && p2->isNucleon())) {
103             inelastic = etaNToPiN(p1,p2) + eta << 103       inelastic = etaNToPiN(p1,p2) + etaNToPiPiN(p1,p2);
104         } else if((p1->isNucleon() && p2->isOm << 104     } else if((p1->isNucleon() && p2->isOmega()) ||
105                   (p1->isOmega() && p2->isNucl << 105           (p1->isOmega() && p2->isNucleon())) {
106             inelastic = omegaNInelastic(p1,p2) << 106       inelastic = omegaNInelastic(p1,p2);
107         } else if((p1->isNucleon() && p2->isEt << 107     } else if((p1->isNucleon() && p2->isEtaPrime()) ||
108                   (p1->isEtaPrime() && p2->isN << 108           (p1->isEtaPrime() && p2->isNucleon())) {
109             inelastic = etaPrimeNToPiN(p1,p2); << 109       inelastic = etaPrimeNToPiN(p1,p2);
110         } else {                               << 110     } else {
111             inelastic = 0.;                    << 111       inelastic = 0.;
112         }                                      << 112     }
113                                                << 113     
114         return inelastic + elastic(p1, p2);    << 114     return inelastic + elastic(p1, p2);
115     }                                          << 115   } 
116                                                << 116 
117                                                << 117   
118     G4double CrossSectionsMultiPionsAndResonan << 118   G4double CrossSectionsMultiPionsAndResonances::elastic(Particle const * const p1, Particle const * const p2) {
119         if((p1->isNucleon()||p1->isDelta()) && << 119     if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta())){
120             return CrossSectionsMultiPions::el << 120       return CrossSectionsMultiPions::elastic(p1, p2);
121         }                                      << 121     }
122         else if ((p1->isNucleon() && p2->isPio << 122     else if ((p1->isNucleon() && p2->isPion()) || (p2->isNucleon() && p1->isPion())){
123             return CrossSectionsMultiPions::el << 123       return CrossSectionsMultiPions::elastic(p1, p2);
124         }                                      << 124     }
125         else if ((p1->isNucleon() && p2->isEta << 125     else if ((p1->isNucleon() && p2->isEta()) || (p2->isNucleon() && p1->isEta())){
126             return etaNElastic(p1, p2);        << 126       return etaNElastic(p1, p2);
127         }                                      << 127     }
128         else if ((p1->isNucleon() && p2->isOme << 128     else if ((p1->isNucleon() && p2->isOmega()) || (p2->isNucleon() && p1->isOmega())){
129             return omegaNElastic(p1, p2);      << 129       return omegaNElastic(p1, p2);
130         }                                      << 130     }
131         else {                                 << 131     else {
132             return 0.0;                        << 132       return 0.0;
133         }                                      << 133     }
134     }                                          << 134   }
135                                                << 135   
136                                                << 136   
137     G4double CrossSectionsMultiPionsAndResonan << 137   G4double CrossSectionsMultiPionsAndResonances::piNToxPiN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
138         //                                     << 138       //
139         //     pion-Nucleon producing xpi pion << 139       //     pion-Nucleon producing xpi pions cross sections (corrected due to eta and omega)
140         //                                     << 140       //
141 // assert(xpi>1 && xpi<=nMaxPiPiN);               141 // assert(xpi>1 && xpi<=nMaxPiPiN);
142 // assert((particle1->isNucleon() && particle2    142 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
143                                                   143 
144         const G4double oldXS2Pi=CrossSectionsM << 144       const G4double oldXS2Pi=CrossSectionsMultiPions::piNToxPiN(2,particle1, particle2);
145         const G4double oldXS3Pi=CrossSectionsM << 145       const G4double oldXS3Pi=CrossSectionsMultiPions::piNToxPiN(3,particle1, particle2);
146         const G4double oldXS4Pi=CrossSectionsM << 146       const G4double oldXS4Pi=CrossSectionsMultiPions::piNToxPiN(4,particle1, particle2);
147         const G4double xsEta=piNToEtaN(particl << 147       const G4double xsEta=piNToEtaN(particle1, particle2);
148         const G4double xsOmega=piNToOmegaN(par << 148       const G4double xsOmega=piNToOmegaN(particle1, particle2);
149         G4double newXS2Pi=0.;                  << 149       G4double newXS2Pi=0.; 
150         G4double newXS3Pi=0.;                  << 150       G4double newXS3Pi=0.; 
151         G4double newXS4Pi=0.;                  << 151       G4double newXS4Pi=0.;
152                                                << 152       
153         if (xpi == 2) {                        << 153       if (xpi == 2) {
154             if (oldXS4Pi != 0.)                << 154         if (oldXS4Pi != 0.)
155                 newXS2Pi=oldXS2Pi;             << 155           newXS2Pi=oldXS2Pi;
156             else if (oldXS3Pi != 0.) {         << 156         else if (oldXS3Pi != 0.) {
157                 newXS3Pi=oldXS3Pi-xsEta-xsOmeg << 157           newXS3Pi=oldXS3Pi-xsEta-xsOmega;
158                 if (newXS3Pi < 1.e-09)         << 158           if (newXS3Pi < 1.e-09)
159                     newXS2Pi=oldXS2Pi-(xsEta+x << 159             newXS2Pi=oldXS2Pi-(xsEta+xsOmega-oldXS3Pi);
160                 else                           << 160           else
161                     newXS2Pi=oldXS2Pi;         << 161             newXS2Pi=oldXS2Pi;
162             }                                  << 162         }
163             else {                             << 163         else { 
164                 newXS2Pi=oldXS2Pi-xsEta-xsOmeg << 164           newXS2Pi=oldXS2Pi-xsEta-xsOmega;
165             if (newXS2Pi < 1.e-09)             << 165           if (newXS2Pi < 1.e-09)
166                     newXS2Pi=0.;               << 166             newXS2Pi=0.;
167             }                                  << 167         }
168             return newXS2Pi;                   << 168         return newXS2Pi;
169         }                                      << 169       }                     
170         else if (xpi == 3) {                   << 170       else if (xpi == 3) {
171             if (oldXS4Pi != 0.) {              << 171         if (oldXS4Pi != 0.) {
172                 newXS4Pi=oldXS4Pi-xsEta-xsOmeg << 172           newXS4Pi=oldXS4Pi-xsEta-xsOmega;
173                 if (newXS4Pi < 1.e-09)         << 173           if (newXS4Pi < 1.e-09)
174                     newXS3Pi=oldXS3Pi-(xsEta+x << 174             newXS3Pi=oldXS3Pi-(xsEta+xsOmega-oldXS4Pi);
175                 else                           << 175           else
176                     newXS3Pi=oldXS3Pi;         << 176             newXS3Pi=oldXS3Pi;
177             }                                  << 177         }
178             else {                             << 178         else { 
179                 newXS3Pi=oldXS3Pi-xsEta-xsOmeg << 179           newXS3Pi=oldXS3Pi-xsEta-xsOmega;
180                 if (newXS3Pi < 1.e-09)         << 180           if (newXS3Pi < 1.e-09)
181                     newXS3Pi=0.;               << 181             newXS3Pi=0.;
182             }                                  << 182         }
183             return newXS3Pi;                   << 183         return newXS3Pi;
184         }                                      << 184       }                 
185         else if (xpi == 4) {                   << 185       else if (xpi == 4) {
186             newXS4Pi=oldXS4Pi-xsEta-xsOmega;   << 186         newXS4Pi=oldXS4Pi-xsEta-xsOmega;
187             if (newXS4Pi < 1.e-09)             << 187         if (newXS4Pi < 1.e-09)
188                 newXS4Pi=0.;                   << 188           newXS4Pi=0.;
189             return newXS4Pi;                   << 189         return newXS4Pi;
190         }                                      << 190       }     
191         else // should never reach this point  << 191       else // should never reach this point
192             return 0.;                         << 192         return 0.;
193     }                                          << 193     }
194                                                << 194   
195     G4double CrossSectionsMultiPionsAndResonan << 195   G4double CrossSectionsMultiPionsAndResonances::piNToEtaN(Particle const * const particle1, Particle const * const particle2) {
196         //                                     << 196     //
197         //     Pion-Nucleon producing Eta cros << 197     //     Pion-Nucleon producing Eta cross sections
198         //                                     << 198     //
199 // assert((particle1->isNucleon() && particle2    199 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
200                                                << 200     
201         G4double sigma;                        << 201     G4double sigma;
202         sigma=piMinuspToEtaN(particle1,particl << 202     sigma=piMinuspToEtaN(particle1,particle2);
203                                                << 203     
204         const G4int isoin = ParticleTable::get << 204     const G4int isoin = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
205                                                << 205     
206         if (isoin == -1) {                     << 206     if (isoin == -1) {
207             if ((particle1->getType()) == Prot << 207       if ((particle1->getType()) == Proton || (particle2->getType()) == Proton) return sigma;
208             else return 0.5 * sigma;           << 208       else return 0.5 * sigma;
209         }                                      << 209     }
210         else if (isoin == 1) {                 << 210     else if (isoin == 1) {
211             if ((particle1->getType()) == Neut << 211       if ((particle1->getType()) == Neutron || (particle2->getType()) == Neutron) return sigma;
212             else return 0.5 * sigma;           << 212       else return 0.5 * sigma;
213         }                                      << 213     }
214         else return 0. ; // should never retur << 214     else return 0. ; // should never return 0. (?) // pi+ p and pi- n return 0.
215                                                << 215     
216 //        return sigma;                        << 216 //    return sigma;
217     }                                          << 217   }
218                                                << 218   
219     G4double CrossSectionsMultiPionsAndResonan << 219   G4double CrossSectionsMultiPionsAndResonances::piNToOmegaN(Particle const * const particle1, Particle const * const particle2) {
220         //                                     << 220     //
221         //     Pion-Nucleon producing Omega cr << 221     //     Pion-Nucleon producing Omega cross sections
222         //                                     << 222     //
223 // assert((particle1->isNucleon() && particle2    223 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
224                                                << 224     
225         G4double sigma;                        << 225     G4double sigma;
226         sigma=piMinuspToOmegaN(particle1,parti << 226     sigma=piMinuspToOmegaN(particle1,particle2);
227                                                << 227     
228         const G4int isoin = ParticleTable::get << 228     const G4int isoin = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
229                                                << 229     
230         if (isoin == -1) {                     << 230   if (isoin == -1) {
231             if ((particle1->getType()) == Prot << 231       if ((particle1->getType()) == Proton || (particle2->getType()) == Proton) return sigma;
232             else return 0.5 * sigma;           << 232       else return 0.5 * sigma;
233         }                                      << 233     }
234         else if (isoin == 1) {                 << 234     else if (isoin == 1) {
235             if ((particle1->getType()) == Neut << 235       if ((particle1->getType()) == Neutron || (particle2->getType()) == Neutron) return sigma;
236             else return 0.5 * sigma;           << 236       else return 0.5 * sigma;
237         }                                      << 237     }
238         else return 0. ; // should never retur << 238     else return 0. ; // should never return 0. (?) // pi+ p and pi- n return 0.
239                                                << 239     
240 //        return sigma;                        << 240 //    return sigma;
241     }                                          << 241   }
242                                                << 242   
243 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT    243 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
244     G4double CrossSectionsMultiPionsAndResonan    244     G4double CrossSectionsMultiPionsAndResonances::piNToEtaPrimeN(Particle const * const /*particle1*/, Particle const * const /*particle2*/) {
245 #else                                             245 #else
246     G4double CrossSectionsMultiPionsAndResonan << 246   G4double CrossSectionsMultiPionsAndResonances::piNToEtaPrimeN(Particle const * const particle1, Particle const * const particle2) {
247 #endif                                            247 #endif
248         //                                     << 248     //
249         //     Pion-Nucleon producing EtaPrime << 249     //     Pion-Nucleon producing EtaPrime cross sections
250         //                                     << 250     //
251 // assert((particle1->isNucleon() && particle2    251 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
252                                                << 252     
253         return 0.;                             << 253     return 0.;
254     }                                          << 254   }
255                                                << 255   
256     G4double CrossSectionsMultiPionsAndResonan << 256   G4double CrossSectionsMultiPionsAndResonances::etaNToPiN(Particle const * const particle1, Particle const * const particle2) {
257         //                                     << 257  //
258         //     Eta-Nucleon producing Pion cros << 258  //     Eta-Nucleon producing Pion cross sections
259         //                                     << 259  //
260 // assert((particle1->isNucleon() && particle2    260 // assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
261                                                << 261  
262         const Particle *eta;                   << 262  const Particle *eta;
263         const Particle *nucleon;               << 263  const Particle *nucleon;
264                                                << 264  
265         if (particle1->isEta()) {              << 265  if (particle1->isEta()) {
266          eta = particle1;                      << 266   eta = particle1;
267          nucleon = particle2;                  << 267   nucleon = particle2;
268         }                                      << 268  }
269         else {                                 << 269  else {
270          eta = particle2;                      << 270   eta = particle2;
271          nucleon = particle1;                  << 271   nucleon = particle1;
272         }                                      << 272  }
273                                                << 273  
274         const G4double pLab = KinematicsUtils: << 274  const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
275         G4double sigma=0.;                     << 275  G4double sigma=0.;   
276                                                << 276  
277         if (pLab <= 574.)                      << 277  if (pLab <= 574.) 
278          sigma= 1.511147E-13*std::pow(pLab,6)- << 278   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;
279         else if (pLab <= 850.)                 << 279  else if (pLab <= 850.) 
280          sigma= -8.00018E-14*std::pow(pLab,6)+ << 280   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;
281         else if (pLab <= 1300.)                << 281  else if (pLab <= 1300.)
282          sigma= 6.56364E-09*std::pow(pLab,3)-  << 282   sigma= 6.56364E-09*std::pow(pLab,3)- 2.07653E-05*std::pow(pLab,2)+ 1.84148E-02*pLab- 1.70427E+00;  
283         else {                                 << 283  else {
284             G4double ECM=KinematicsUtils::tota << 284   G4double ECM=KinematicsUtils::totalEnergyInCM(eta, nucleon);
285             G4double massPiZero=ParticleTable: << 285     G4double massPiZero=ParticleTable::getINCLMass(PiZero);
286             G4double massPiMinus=ParticleTable << 286     G4double massPiMinus=ParticleTable::getINCLMass(PiMinus);
287             G4double massProton=ParticleTable: << 287     G4double massProton=ParticleTable::getINCLMass(Proton);
288             G4double masseta;                  << 288     G4double masseta;
289             G4double massnucleon;              << 289     G4double massnucleon;
290             G4double pCM_eta;                  << 290     G4double pCM_eta;
291             masseta=eta->getMass();            << 291   masseta=eta->getMass();
292             massnucleon=nucleon->getMass();    << 292   massnucleon=nucleon->getMass();
293             pCM_eta=KinematicsUtils::momentumI << 293   pCM_eta=KinematicsUtils::momentumInCM(ECM, masseta, massnucleon);    
294             G4double pCM_PiZero=KinematicsUtil << 294   G4double pCM_PiZero=KinematicsUtils::momentumInCM(ECM, massPiZero, massProton);
295             G4double pCM_PiMinus=KinematicsUti << 295   G4double pCM_PiMinus=KinematicsUtils::momentumInCM(ECM, massPiMinus, massProton); // = pCM_PiPlus (because massPiMinus = massPiPlus) 
296             sigma = (piMinuspToEtaN(ECM)/2.) * << 296   sigma = (piMinuspToEtaN(ECM)/2.) * std::pow((pCM_PiZero/pCM_eta), 2) + piMinuspToEtaN(ECM) * std::pow((pCM_PiMinus/pCM_eta), 2);      
297         }                                      << 297  } 
298         if (sigma < 0.) sigma=0.;              << 298  if (sigma < 0.) sigma=0.;
299                                                << 299  return sigma;
300         return sigma;                          << 300 }
301     }                                          << 
302                                                   301 
303     G4double CrossSectionsMultiPionsAndResonan << 302   G4double CrossSectionsMultiPionsAndResonances::etaNToPiPiN(Particle const * const particle1, Particle const * const particle2) {
304         //                                     << 303       //
305         //     Eta-Nucleon producing Two Pions << 304       //     Eta-Nucleon producing Two Pions cross sections
306         //                                     << 305       //
307 // assert((particle1->isNucleon() && particle2    306 // assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
308                                                << 307             
309         G4double sigma=0.;                     << 308       G4double sigma=0.;    
310                                                << 309       
311         const Particle *eta;                   << 310       const Particle *eta;
312         const Particle *nucleon;               << 311       const Particle *nucleon;
313                                                << 312 
314         if (particle1->isEta()) {              << 313       if (particle1->isEta()) {
315             eta = particle1;                   << 314         eta = particle1;
316             nucleon = particle2;               << 315         nucleon = particle2;
317         }                                      << 316       }
318         else {                                 << 317      else {
319             eta = particle2;                   << 318         eta = particle2;
320             nucleon = particle1;               << 319         nucleon = particle1;
321         }                                      << 320       }
322                                                << 321         
323         const G4double pLab = KinematicsUtils: << 322      const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
324                                                << 323       
325         if (pLab < 450.)                       << 324    if (pLab < 450.) 
326             sigma = 2.01854221E-13*std::pow(pL << 325         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;       
327         else if (pLab < 600.)                  << 326       else if (pLab < 600.) 
328             sigma = 2.01854221E-13*std::pow(45 << 327         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;       
329         else if (pLab <= 1300.)                << 328       else if (pLab <= 1300.) 
330             sigma = -6.32793049e-16*std::pow(p << 329         sigma = -6.32793049e-16*std::pow(pLab,6) + 3.95985900e-12*std::pow(pLab,5) - 1.01727714e-8*std::pow(pLab,4) + 
331             1.37055547e-05*std::pow(pLab,3) -  << 330                  1.37055547e-05*std::pow(pLab,3) - 1.01830486e-02*std::pow(pLab,2) + 3.93492126*pLab - 609.447145;
332         else                                   << 331    else 
333             sigma = etaNToPiN(particle1,partic << 332         sigma = etaNToPiN(particle1,particle2);
334                                                << 333       
335         if (sigma < 0.) sigma = 0.;            << 334       if (sigma < 0.) sigma = 0.;
336         return sigma; // Parameterization from << 335       return sigma; // Parameterization from the ANL-Osaka DCC model [PRC88(2013)035209] - eta p --> "pi+pi0 n" + "pi0 pi0 p" total XS
337     }                                          << 336     }
338                                                << 337     
339                                                << 338     
340     G4double CrossSectionsMultiPionsAndResonan << 339   G4double CrossSectionsMultiPionsAndResonances::etaNElastic(Particle const * const particle1, Particle const * const particle2) {
341         //                                     << 340       //
342         //     Eta-Nucleon elastic cross secti << 341       //     Eta-Nucleon elastic cross sections
343         //                                     << 342       //
344 // assert((particle1->isNucleon() && particle2    343 // assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
345                                                << 344       
346         G4double sigma=0.;                     << 345       G4double sigma=0.;    
347                                                << 346       
348         const Particle *eta;                   << 347       const Particle *eta;
349         const Particle *nucleon;               << 348       const Particle *nucleon;
350                                                << 349       
351         if (particle1->isEta()) {              << 350       if (particle1->isEta()) {
352             eta = particle1;                   << 351         eta = particle1;
353             nucleon = particle2;               << 352         nucleon = particle2;
354         }                                      << 353       }
355         else {                                 << 354      else {
356             eta = particle2;                   << 355         eta = particle2;
357             nucleon = particle1;               << 356         nucleon = particle1;
358         }                                      << 357       }
359                                                << 358       
360         const G4double pLab = KinematicsUtils: << 359      const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
361                                                << 360       
362         if (pLab < 700.)                       << 361    if (pLab < 700.) 
363             sigma = 3.6838e-15*std::pow(pLab,6 << 362         sigma = 3.6838e-15*std::pow(pLab,6) - 9.7815e-12*std::pow(pLab,5) + 9.7914e-9*std::pow(pLab,4) - 
364             4.3222e-06*std::pow(pLab,3) + 7.91 << 363                 4.3222e-06*std::pow(pLab,3) + 7.9188e-04*std::pow(pLab,2) - 1.8379e-01*pLab + 84.965;     
365         else if (pLab < 1400.)                 << 364       else if (pLab < 1400.) 
366             sigma = 3.562630e-16*std::pow(pLab << 365         sigma = 3.562630e-16*std::pow(pLab,6) - 2.384766e-12*std::pow(pLab,5) + 6.601312e-9*std::pow(pLab,4) - 
367             9.667078e-06*std::pow(pLab,3) + 7. << 366                 9.667078e-06*std::pow(pLab,3) + 7.894845e-03*std::pow(pLab,2) - 3.409200*pLab + 609.8501;
368         else if (pLab < 2025.)                 << 367       else if (pLab < 2025.)
369             sigma = -1.041950E-03*pLab + 2.110 << 368         sigma = -1.041950E-03*pLab + 2.110529E+00;
370         else                                   << 369       else  
371             sigma=0.;                          << 370        sigma=0.;
372                                                << 371         
373         if (sigma < 0.) sigma = 0.;            << 372       if (sigma < 0.) sigma = 0.;
374             return sigma; // Parameterization  << 373       return sigma; // Parameterization from the ANL-Osaka DCC model [PRC88(2013)035209]
375         }                                      << 374     }
376                                                << 375 
377     G4double CrossSectionsMultiPionsAndResonan << 376   G4double CrossSectionsMultiPionsAndResonances::omegaNInelastic(Particle const * const particle1, Particle const * const particle2) {
378         //                                     << 377     //
379         //     Omega-Nucleon inelastic cross s << 378     //     Omega-Nucleon inelastic cross sections
380         //                                     << 379     //
381 // assert((particle1->isNucleon() && particle2    380 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
382                                                << 381     
383         G4double sigma=0.;                     << 382     G4double sigma=0.;    
384                                                << 383     
385         const Particle *omega;                 << 384     const Particle *omega;
386         const Particle *nucleon;               << 385     const Particle *nucleon;
387                                                << 386     
388         if (particle1->isOmega()) {            << 387     if (particle1->isOmega()) {
389             omega = particle1;                 << 388       omega = particle1;
390             nucleon = particle2;               << 389       nucleon = particle2;
391             }                                  << 390       }
392         else {                                 << 391     else {
393             omega = particle2;                 << 392       omega = particle2;
394             nucleon = particle1;               << 393       nucleon = particle1;
395         }                                      << 394     }
396                                                << 395     
397         const G4double pLab = KinematicsUtils: << 396     const G4double pLab = KinematicsUtils::momentumInLab(omega, nucleon)/1000.; // GeV/c
398                                                << 397     
399         sigma = 20. + 4.0/pLab;    // Eq.(24)  << 398     sigma = 20. + 4.0/pLab; // Eq.(24) in G.I. Lykasov et al., EPJA 6, 71-81 (1999)
400                                                << 399     
401         return sigma;                          << 400     return sigma;
402     }                                          << 401   }
403                                                   402 
404                                                   403   
405     G4double CrossSectionsMultiPionsAndResonan << 404   G4double CrossSectionsMultiPionsAndResonances::omegaNElastic(Particle const * const particle1, Particle const * const particle2) {
406         //                                     << 405     //
407         //     Omega-Nucleon elastic cross sec << 406     //     Omega-Nucleon elastic cross sections
408         //                                     << 407     //
409 // assert((particle1->isNucleon() && particle2    408 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
410                                                << 409     
411         G4double sigma=0.;                     << 410     G4double sigma=0.;    
412                                                << 411     
413         const Particle *omega;                 << 412     const Particle *omega;
414         const Particle *nucleon;               << 413     const Particle *nucleon;
415                                                << 414     
416         if (particle1->isOmega()) {            << 415     if (particle1->isOmega()) {
417             omega = particle1;                 << 416       omega = particle1;
418             nucleon = particle2;               << 417       nucleon = particle2;
419         }                                      << 418     }
420         else {                                 << 419     else {
421             omega = particle2;                 << 420       omega = particle2;
422             nucleon = particle1;               << 421       nucleon = particle1;
423         }                                      << 422     }
424                                                << 423     
425         const G4double pLab = KinematicsUtils: << 424     const G4double pLab = KinematicsUtils::momentumInLab(omega, nucleon)/1000.; // GeV/c
426                                                << 425     
427         sigma = 5.4 + 10.*std::exp(-0.6*pLab); << 426     sigma = 5.4 + 10.*std::exp(-0.6*pLab);  // Eq.(21) in G.I. Lykasov et al., EPJA 6, 71-81 (1999)
428                                                << 427     
429         return sigma;                          << 428     return sigma;
430     }                                          << 429   }
431                                                   430   
432                                                << 431     
433     G4double CrossSectionsMultiPionsAndResonan << 432   G4double CrossSectionsMultiPionsAndResonances::omegaNToPiN(Particle const * const particle1, Particle const * const particle2) {
434         //                                     << 433     //
435         //     Omega-Nucleon producing Pion cr << 434     //     Omega-Nucleon producing Pion cross sections
436         //                                     << 435     //
437 // assert((particle1->isNucleon() && particle2    436 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
438                                                << 437     
439         G4double ECM=KinematicsUtils::totalEne << 438     G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
440                                                << 439     
441         G4double massPiZero=ParticleTable::get << 440     G4double massPiZero=ParticleTable::getINCLMass(PiZero);
442         G4double massPiMinus=ParticleTable::ge << 441     G4double massPiMinus=ParticleTable::getINCLMass(PiMinus);
443         G4double massProton=ParticleTable::get << 442     G4double massProton=ParticleTable::getINCLMass(Proton);
444                                                << 443     
445         G4double massomega;                    << 444     G4double massomega;
446         G4double massnucleon;                  << 445     G4double massnucleon;
447         G4double pCM_omega;                    << 446     G4double pCM_omega;
448         G4double pLab_omega;                   << 447     G4double pLab_omega;
449                                                << 448     
450         G4double sigma=0.;                     << 449     G4double sigma=0.;    
451                                                << 450     
452         if (particle1->isOmega()) {            << 451     if (particle1->isOmega()) {
453             massomega=particle1->getMass();    << 452       massomega=particle1->getMass();
454             massnucleon=particle2->getMass();  << 453       massnucleon=particle2->getMass();
455         }                                      << 454     }
456         else {                                 << 455     else {
457             massomega=particle2->getMass();    << 456       massomega=particle2->getMass();
458             massnucleon=particle1->getMass();  << 457       massnucleon=particle1->getMass();
459         }                                      << 458     }
460         pCM_omega=KinematicsUtils::momentumInC << 459     pCM_omega=KinematicsUtils::momentumInCM(ECM, massomega, massnucleon);   
461         pLab_omega=KinematicsUtils::momentumIn << 460     pLab_omega=KinematicsUtils::momentumInLab(ECM*ECM, massomega, massnucleon);   
462                                                << 461     
463         G4double pCM_PiZero=KinematicsUtils::m << 462     G4double pCM_PiZero=KinematicsUtils::momentumInCM(ECM, massPiZero, massProton);
464         G4double pCM_PiMinus=KinematicsUtils:: << 463     G4double pCM_PiMinus=KinematicsUtils::momentumInCM(ECM, massPiMinus, massProton); // = pCM_PiPlus (because massPiMinus = massPiPlus)
465                                                << 464     
466         sigma = (piMinuspToOmegaN(ECM)/2.) * s << 465     sigma = (piMinuspToOmegaN(ECM)/2.) * std::pow((pCM_PiZero/pCM_omega), 2) 
467         + piMinuspToOmegaN(ECM) * std::pow((pC << 466     + piMinuspToOmegaN(ECM) * std::pow((pCM_PiMinus/pCM_omega), 2);     
468                                                   467   
469         if (sigma > omegaNInelastic(particle1, << 468   if (sigma > omegaNInelastic(particle1, particle2) || (pLab_omega < 200.)) {
470 //        if (sigma > omegaNInelastic(particle << 469 //  if (sigma > omegaNInelastic(particle1, particle2)) {
471             sigma = omegaNInelastic(particle1, << 470    sigma = omegaNInelastic(particle1, particle2);
472         }                                      << 471   }  
473                                                   472   
474         return sigma;                          << 473     return sigma;
475     }                                          << 474   }
476                                                << 475   
477                                                   476   
478     G4double CrossSectionsMultiPionsAndResonan << 477   G4double CrossSectionsMultiPionsAndResonances::omegaNToPiPiN(Particle const * const particle1, Particle const * const particle2) {
479         //                                     << 478     //
480         //     Omega-Nucleon producing 2 PionS << 479     //     Omega-Nucleon producing 2 PionS cross sections
481         //                                     << 480     //
482 // assert((particle1->isNucleon() && particle2    481 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
483                                                << 482     
484         G4double sigma=0.;                     << 483     G4double sigma=0.;    
485                                                << 484     
486         sigma = omegaNInelastic(particle1,part << 485     sigma = omegaNInelastic(particle1,particle2) - omegaNToPiN(particle1,particle2) ;
487                                                << 486     
488         return sigma;                          << 487     return sigma;
489     }                                          << 488   }
490                                                   489   
491                                                << 490   
492 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT    491 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
493   G4double CrossSectionsMultiPionsAndResonance    492   G4double CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(Particle const * const /*particle1*/, Particle const * const /*particle2*/) {
494 #else                                             493 #else
495   G4double CrossSectionsMultiPionsAndResonance    494   G4double CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(Particle const * const particle1, Particle const * const particle2) {
496 #endif                                            495 #endif
497         //                                     << 496     //
498         //     EtaPrime-Nucleon producing Pion << 497     //     EtaPrime-Nucleon producing Pion cross sections
499         //                                     << 498     //
500 // assert((particle1->isNucleon() && particle2    499 // assert((particle1->isNucleon() && particle2->isEtaPrime()) || (particle1->isEtaPrime() && particle2->isNucleon()));
501                                                << 500     
502         return 0.;                             << 501     return 0.;
503     }                                          << 502   } 
504                                                << 503   
505     G4double CrossSectionsMultiPionsAndResonan << 504   G4double CrossSectionsMultiPionsAndResonances::piMinuspToEtaN(Particle const * const particle1, Particle const * const particle2) {
506         //                                     << 505     //
507         //     Pion-Nucleon producing Eta cros << 506     //     Pion-Nucleon producing Eta cross sections
508         //                                     << 507     //
509 // assert((particle1->isNucleon() && particle2    508 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
510                                                << 509     
511         G4double masspion;                     << 510     G4double masspion;
512         G4double massnucleon;                  << 511     G4double massnucleon;
513         if (particle1->isPion()) {             << 512     if (particle1->isPion()) {
514             masspion=particle1->getMass();     << 513         masspion=particle1->getMass();
515             massnucleon=particle2->getMass();  << 514         massnucleon=particle2->getMass();
516         }                                      << 515     }
517         else {                                 << 516     else {
518             masspion=particle2->getMass();     << 517         masspion=particle2->getMass();
519             massnucleon=particle1->getMass();  << 518         massnucleon=particle1->getMass();
520         }                                      << 519     }
521                                                << 520     
522         G4double ECM=KinematicsUtils::totalEne << 521     G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
523         G4double plab=KinematicsUtils::momentu << 522     G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.; // GeV/c
524                                                << 523     
525         G4double sigma;                        << 524     G4double sigma;
526                                                << 525     
527 // new parameterization (JCD) - end of februar    526 // new parameterization (JCD) - end of february 2016
528         if (ECM < 1486.5) sigma=0.;            << 527     if (ECM < 1486.5) sigma=0.;
529         else                                   << 528     else
530         {                                      << 529     {
531             if (ECM < 1535.)                   << 530       if (ECM < 1535.)
532             {                                  << 531       {
533                 sigma = -0.0000003689197974814 << 532         sigma = -0.0000003689197974814*std::pow(ECM,4) + 0.002260193900097*std::pow(ECM,3) - 5.193105877187*std::pow(ECM,2) + 5303.505273919*ECM - 2031265.900648;
534             }                                  << 533       }
535             else if (ECM < 1670.)              << 534       else if (ECM < 1670.)
536             {                                  << 535       {
537                 sigma = -0.0000000337986446*st << 536         sigma = -0.0000000337986446*std::pow(ECM,4) + 0.000218279989*std::pow(ECM,3) - 0.528276144*std::pow(ECM,2) + 567.828367*ECM - 228709.42;
538             }                                  << 537       }
539             else if (ECM < 1714.)              << 538       else if (ECM < 1714.)
540             {                                  << 539       {
541                 sigma =  0.000003737765*std::p << 540         sigma =  0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
542             }                                  << 541       }
543             else sigma=1.47*std::pow(plab, -1. << 542       else sigma=1.47*std::pow(plab, -1.68);
544         }                                      << 543     }   
545 //                                             << 544 //    
546                                                << 545 
547         return sigma;                          << 546     return sigma;
548     }                                          << 547   }
549                                                << 548   
550     G4double CrossSectionsMultiPionsAndResonan << 549   G4double CrossSectionsMultiPionsAndResonances::piMinuspToEtaN(const G4double ECM) {
551         //                                     << 550     //
552         //     Pion-Nucleon producing Eta cros << 551     //     Pion-Nucleon producing Eta cross sections
553         //                                     << 552     //
554                                                << 553     
555         const G4double masspion = ParticleTabl << 554   const G4double masspion = ParticleTable::getRealMass(PiMinus);    
556         const G4double massnucleon = ParticleT << 555   const G4double massnucleon = ParticleTable::getRealMass(Proton);    
557                                                   556   
558         G4double plab=KinematicsUtils::momentu << 557   G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.;  // GeV/c
559                                                << 558     
560         G4double sigma;                        << 559     G4double sigma;
561                                                << 560     
562 // new parameterization (JCD) - end of februar    561 // new parameterization (JCD) - end of february 2016
563         if (ECM < 1486.5) sigma=0.;            << 562     if (ECM < 1486.5) sigma=0.;
564         else                                   << 563     else
565         {                                      << 564     {
566             if (ECM < 1535.)                   << 565       if (ECM < 1535.)
567             {                                  << 566       {
568                 sigma = -0.0000003689197974814 << 567         sigma = -0.0000003689197974814*std::pow(ECM,4) + 0.002260193900097*std::pow(ECM,3) - 5.193105877187*std::pow(ECM,2) + 5303.505273919*ECM - 2031265.900648;
569             }                                  << 568       }
570             else if (ECM < 1670.)              << 569       else if (ECM < 1670.)
571             {                                  << 570       {
572                 sigma = -0.0000000337986446*st << 571         sigma = -0.0000000337986446*std::pow(ECM,4) + 0.000218279989*std::pow(ECM,3) - 0.528276144*std::pow(ECM,2) + 567.828367*ECM - 228709.42;
573             }                                  << 572       }
574             else if (ECM < 1714.)              << 573       else if (ECM < 1714.)
575             {                                  << 574       {
576                 sigma =  0.000003737765*std::p << 575         sigma =  0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
577             }                                  << 576       }
578             else sigma=1.47*std::pow(plab, -1. << 577       else sigma=1.47*std::pow(plab, -1.68);
579         }                                      << 578     }   
580                                                << 579 //    
581         return sigma;                          << 580 
582     }                                          << 581     return sigma;
583                                                << 582   }
584     G4double CrossSectionsMultiPionsAndResonan << 583 
585         //                                     << 584   G4double CrossSectionsMultiPionsAndResonances::piMinuspToOmegaN(Particle const * const particle1, Particle const * const particle2) {
586         //     Pion-Nucleon producing Omega cr << 585     //
587         //                                     << 586     //     Pion-Nucleon producing Omega cross sections
                                                   >> 587     //
588 // assert((particle1->isNucleon() && particle2    588 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
589 //jcd to be removed                               589 //jcd to be removed
590 //        return 0.;                           << 590 //  return 0.;
591 //jcd                                             591 //jcd    
592                                                << 592     
593 //        G4double param=1.095 ; // Deneye (Th << 593 //    G4double param=1.095 ; // Deneye (Thesis)
594         G4double param=1.0903 ; // JCD (thresh << 594     G4double param=1.0903 ; // JCD (threshold taken into account)
595                                                   595   
596         G4double masspion;                     << 596   G4double masspion;
597         G4double massnucleon;                  << 597     G4double massnucleon;
598         if (particle1->isPion()) {             << 598     if (particle1->isPion()) {
599             masspion=particle1->getMass();     << 599         masspion=particle1->getMass();
600             massnucleon=particle2->getMass();  << 600         massnucleon=particle2->getMass();
601         }                                      << 601     }
602         else {                                 << 602     else {
603             masspion=particle2->getMass();     << 603         masspion=particle2->getMass();
604             massnucleon=particle1->getMass();  << 604         massnucleon=particle1->getMass();
605         }                                      << 605     }
606         G4double ECM=KinematicsUtils::totalEne << 606     G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
607         G4double plab=KinematicsUtils::momentu << 607     G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.;  // GeV/c
608                                                << 608     
609         G4double sigma;                        << 609     G4double sigma;
610         if (plab < param) sigma=0.;            << 610     if (plab < param) sigma=0.;
611         else sigma=13.76*(plab-param)/(std::po << 611     else sigma=13.76*(plab-param)/(std::pow(plab, 3.33) - 1.07); // Phys. Rev. C 41, 1701–1718 (1990)
612                                                   612   
613         return sigma;                          << 613     return sigma;
614 }                                                 614 }
615     G4double CrossSectionsMultiPionsAndResonan << 615       G4double CrossSectionsMultiPionsAndResonances::piMinuspToOmegaN(const G4double ECM) {
616         //                                     << 616         //
617         //     Pion-Nucleon producing Omega cr << 617         //     Pion-Nucleon producing Omega cross sections
618         //                                     << 618         //
619 //jcd to be removed                               619 //jcd to be removed
620 //    return 0.;                                  620 //    return 0.;
621 //jcd                                             621 //jcd    
622                                                << 622         
623 //        G4double param=1.095 ; // Deneye (Th << 623 //    G4double param=1.095 ; // Deneye (Thesis)
624         G4double param=1.0903 ; // JCD (thresh << 624     G4double param=1.0903 ; // JCD (threshold taken into account)
625                                                << 625 
626         const G4double masspion = ParticleTabl << 626         const G4double masspion = ParticleTable::getRealMass(PiMinus);    
627         const G4double massnucleon = ParticleT << 627         const G4double massnucleon = ParticleTable::getRealMass(Proton);    
628                                                << 628 
629         G4double plab=KinematicsUtils::momentu << 629         G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.;  // GeV/c
630                                                << 630         
631         G4double sigma;                        << 631         G4double sigma;
632         if (plab < param) sigma=0.;            << 632         if (plab < param) sigma=0.;
633         else sigma=13.76*(plab-param)/(std::po << 633         else sigma=13.76*(plab-param)/(std::pow(plab, 3.33)-1.07);
634                                                << 634         
635         return sigma;                          << 635         return sigma;
636     }                                          << 636 }
637                                                   637 
638     G4double CrossSectionsMultiPionsAndResonan    638     G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaIso(const G4double ener, const G4int iso) {
639                                                << 639     
640         const G4double Ecm=0.001*ener;         << 640     const G4double Ecm=0.001*ener;
641         G4double sNNEta; // pp->pp+eta(+X)     << 641     G4double sNNEta; // pp->pp+eta(+X)
642         G4double sNNEta1; // np->np+eta(+X)    << 642     G4double sNNEta1; // np->np+eta(+X)
643         G4double sNNEta2; // np->d+eta (d will << 643     G4double sNNEta2; // np->d+eta (d will be considered as np - How far is this right?)
644         G4double x=Ecm*Ecm/5.88;               << 644     G4double x=Ecm*Ecm/5.88;
645                                                << 645                 
646 //jcd                                             646 //jcd
647         if (Ecm >= 3.05) {                     << 647   if (Ecm >= 3.05) {
648             sNNEta = 2.5*std::pow((x-1.),1.47) << 648      sNNEta = 2.5*std::pow((x-1.),1.47)*std::pow(x,-1.25)*1000.;
649         }                                      << 649   }      
650         else if(Ecm >= 2.6) {                  << 650   else if(Ecm >= 2.6) {
651             sNNEta = -327.29*Ecm*Ecm*Ecm + 287 << 651      sNNEta = -327.29*Ecm*Ecm*Ecm + 2870.*Ecm*Ecm - 7229.3*Ecm + 5273.3;
652             if (sNNEta <= NNToNNEtaExcluIso(en << 652      if (sNNEta <= NNToNNEtaExcluIso(ener, 2)*1000.) sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
653         }                                      << 653   }
654         else {                                 << 654   else {
655             sNNEta = NNToNNEtaExcluIso(ener, 2 << 655      sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
656         }                                      << 656   }
657 //jcd                                             657 //jcd
658         if (sNNEta < 1.e-9) sNNEta = 0.;       << 658   if (sNNEta < 1.e-9) sNNEta = 0.;
659                                                << 
660         if (iso != 0) {                        << 
661             return sNNEta/1000.; // parameteri << 
662         }                                      << 
663                                                << 
664         if(Ecm >= 6.25) {                      << 
665             sNNEta1 = sNNEta;                  << 
666         }                                      << 
667         else if (Ecm >= 2.6) {                 << 
668             sNNEta1 = sNNEta*std::exp(-(-5.531 << 
669         }                                      << 
670         else if (Ecm >= 2.525) { // = exclusiv << 
671             sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ec << 
672         }                                      << 
673         else { // = exclusive pn               << 
674             sNNEta1 = 17570.217219*Ecm*Ecm - 8 << 
675         }                                      << 
676                                                << 
677         sNNEta2 = -10220.89518466*Ecm*Ecm+5122 << 
678         if (sNNEta2 < 0.) sNNEta2=0.;          << 
679                                                << 
680         sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;   << 
681                                                << 
682         G4double Mn=ParticleTable::getRealMass << 
683         G4double Mp=ParticleTable::getRealMass << 
684         G4double Meta=ParticleTable::getRealMa << 
685         if (sNNEta < 1.e-9 || Ecm < Mn+Mp+Meta << 
686                                                << 
687         return sNNEta/1000.; // parameterizati << 
688     }                                          << 
689                                                   659 
                                                   >> 660   if (iso != 0) {
                                                   >> 661      return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
                                                   >> 662   }
                                                   >> 663 
                                                   >> 664   if(Ecm >= 6.25) {
                                                   >> 665      sNNEta1 = sNNEta;
                                                   >> 666   }
                                                   >> 667   else if (Ecm >= 2.6) {
                                                   >> 668      sNNEta1 = sNNEta*std::exp(-(-5.53151576/Ecm+0.8850425));
                                                   >> 669   }
                                                   >> 670   else if (Ecm >= 2.525) { // = exclusive pn
                                                   >> 671       sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ecm + 56581.54*Ecm*Ecm*Ecm - 270212.6*Ecm*Ecm + 571650.6*Ecm - 451091.6;
                                                   >> 672   }
                                                   >> 673   else { // = exclusive pn
                                                   >> 674       sNNEta1 = 17570.217219*Ecm*Ecm - 84910.985402*Ecm + 102585.55847;
                                                   >> 675   }
                                                   >> 676 
                                                   >> 677   sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
                                                   >> 678   if (sNNEta2 < 0.) sNNEta2=0.;
                                                   >> 679 
                                                   >> 680   sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
                                                   >> 681 
                                                   >> 682   G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
                                                   >> 683   G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
                                                   >> 684   G4double Meta=ParticleTable::getRealMass(Eta)/1000.;
                                                   >> 685     if (sNNEta < 1.e-9 || Ecm < Mn+Mp+Meta) sNNEta = 0.;
                                                   >> 686     
                                                   >> 687     return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
                                                   >> 688   }
                                                   >> 689 
                                                   >> 690   
                                                   >> 691   G4double CrossSectionsMultiPionsAndResonances::NNToNNEta(Particle const * const particle1, Particle const * const particle2) {
                                                   >> 692     
                                                   >> 693   const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
                                                   >> 694   const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
                                                   >> 695     
                                                   >> 696     if (iso != 0) {
                                                   >> 697    return NNToNNEtaIso(ener, iso);
                                                   >> 698     }
                                                   >> 699     else {
                                                   >> 700       return 0.5*(NNToNNEtaIso(ener, 0)+NNToNNEtaIso(ener, 2));
                                                   >> 701     }
                                                   >> 702     }
                                                   >> 703   
                                                   >> 704   G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaExcluIso(const G4double ener, const G4int iso) {
                                                   >> 705         
                                                   >> 706   const G4double Ecm=0.001*ener;
                                                   >> 707   G4double sNNEta; // pp->pp+eta
                                                   >> 708   G4double sNNEta1; // np->np+eta
                                                   >> 709   G4double sNNEta2; // np->d+eta (d wil be considered as np - How far is this right?)
                                                   >> 710         
                                                   >> 711   if(Ecm>=3.875) { // By hand (JCD)
                                                   >> 712      sNNEta = -13.008*Ecm*Ecm + 84.531*Ecm + 36.234;
                                                   >> 713   }
                                                   >> 714   else if(Ecm>=2.725) { // By hand (JCD)
                                                   >> 715      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;
                                                   >> 716   }
                                                   >> 717   else if(Ecm>=2.575) { // By hand (JCD)
                                                   >> 718      sNNEta = -2640.3*Ecm*Ecm + 14692*Ecm - 20225;
                                                   >> 719   }
                                                   >> 720   else {
                                                   >> 721      sNNEta = -147043.497285*std::pow(Ecm,4) + 1487222.5438123*std::pow(Ecm,3) - 5634399.900744*std::pow(Ecm,2) + 9477290.199378*Ecm - 5972174.353438;
                                                   >> 722   }
                                                   >> 723         
                                                   >> 724    G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
                                                   >> 725    G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
                                                   >> 726    G4double Meta=ParticleTable::getRealMass(Eta)/1000.;
                                                   >> 727    G4double Thr0=0.;
                                                   >> 728    if (iso > 0) {
                                                   >> 729     Thr0=2.*Mp+Meta;
                                                   >> 730    }
                                                   >> 731    else if (iso < 0) {
                                                   >> 732     Thr0=2.*Mn+Meta;
                                                   >> 733    }
                                                   >> 734    else {
                                                   >> 735     Thr0=Mn+Mp+Meta;
                                                   >> 736    }
                                                   >> 737 
                                                   >> 738    if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.;  // Thr0: Ecm threshold
                                                   >> 739         
                                                   >> 740    if (iso != 0) {
                                                   >> 741       return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
                                                   >> 742    }
                                                   >> 743         
                                                   >> 744    if(Ecm>=3.9) {
                                                   >> 745       sNNEta1 = sNNEta;
                                                   >> 746    }
                                                   >> 747    else if (Ecm >= 3.5) {
                                                   >> 748       sNNEta1 = -1916.2*Ecm*Ecm*Ecm + 21556.0*Ecm*Ecm - 80828.0*Ecm + 101200.0;
                                                   >> 749    }
                                                   >> 750    else if (Ecm >= 2.525) {
                                                   >> 751       sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ecm + 56581.54*Ecm*Ecm*Ecm - 270212.6*Ecm*Ecm + 571650.6*Ecm - 451091.6;
                                                   >> 752    }
                                                   >> 753    else {
                                                   >> 754       sNNEta1 = 17570.217219*Ecm*Ecm - 84910.985402*Ecm + 102585.55847;
                                                   >> 755   }
690                                                   756     
691     G4double CrossSectionsMultiPionsAndResonan << 757    sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
                                                   >> 758    if (sNNEta2 < 0.) sNNEta2=0.;
                                                   >> 759     
                                                   >> 760    sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
                                                   >> 761    if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.;  // Thr0: Ecm threshold
                                                   >> 762         
                                                   >> 763   return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
692                                                   764 
693         const G4double ener=KinematicsUtils::t << 765 }
694         const G4int iso=ParticleTable::getIsos << 766       
                                                   >> 767   G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaExclu(Particle const * const particle1, Particle const * const particle2) {
                                                   >> 768         
                                                   >> 769    const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
                                                   >> 770    const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
                                                   >> 771         
                                                   >> 772    if (iso != 0) {
                                                   >> 773       return NNToNNEtaExcluIso(ener, iso);
                                                   >> 774    }
                                                   >> 775    else {
                                                   >> 776       return 0.5*(NNToNNEtaExcluIso(ener, 0)+NNToNNEtaExcluIso(ener, 2));
                                                   >> 777    }
                                                   >> 778 }
695                                                   779 
696         if (iso != 0) {                        << 780    
697             return NNToNNEtaIso(ener, iso);    << 781    G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaIso(const G4double ener, const G4int iso) {
698         }                                      << 782     
699         else {                                 << 783     const G4double Ecm=0.001*ener;
700             return 0.5*(NNToNNEtaIso(ener, 0)+ << 784     G4double sNNOmega; // pp->pp+eta(+X)
701         }                                      << 785     G4double sNNOmega1; // np->np+eta(+X)
                                                   >> 786     G4double x=Ecm*Ecm/7.06;
                                                   >> 787     
                                                   >> 788     if(Ecm>4.0) {
                                                   >> 789      sNNOmega = 2.5*std::pow(x-1, 1.47)*std::pow(x, -1.11) ;
702     }                                             790     }
703                                                << 791     else if(Ecm>2.802) { // 2802 MeV -> threshold to open inclusive (based on multipion threshold and omega mass) 
704     G4double CrossSectionsMultiPionsAndResonan << 792      sNNOmega = (568.5254*Ecm*Ecm - 2694.045*Ecm + 3106.247)/1000.;
705                                                << 793      if (sNNOmega <= NNToNNOmegaExcluIso(ener, 2)) sNNOmega = NNToNNOmegaExcluIso(ener, 2);
706         const G4double Ecm=0.001*ener;         << 
707         G4double sNNEta; // pp->pp+eta         << 
708         G4double sNNEta1; // np->np+eta        << 
709         G4double sNNEta2; // np->d+eta (d wil  << 
710                                                << 
711         if(Ecm>=3.875) { // By hand (JCD)      << 
712             sNNEta = -13.008*Ecm*Ecm + 84.531* << 
713         }                                      << 
714         else if(Ecm>=2.725) { // By hand (JCD) << 
715             sNNEta = -913.2809*std::pow(Ecm,5) << 
716         }                                      << 
717         else if(Ecm>=2.575) { // By hand (JCD) << 
718             sNNEta = -2640.3*Ecm*Ecm + 14692*E << 
719         }                                      << 
720         else {                                 << 
721             sNNEta = -147043.497285*std::pow(E << 
722         }                                      << 
723                                                << 
724         G4double Mn=ParticleTable::getRealMass << 
725         G4double Mp=ParticleTable::getRealMass << 
726         G4double Meta=ParticleTable::getRealMa << 
727         G4double Thr0=0.;                      << 
728         if (iso > 0) {                         << 
729             Thr0=2.*Mp+Meta;                   << 
730         }                                      << 
731         else if (iso < 0) {                    << 
732             Thr0=2.*Mn+Meta;                   << 
733         }                                      << 
734         else {                                 << 
735             Thr0=Mn+Mp+Meta;                   << 
736         }                                      << 
737                                                << 
738         if (sNNEta < 1.e-9 || Ecm < Thr0) sNNE << 
739                                                << 
740         if (iso != 0) {                        << 
741             return sNNEta/1000.; // parameteri << 
742         }                                      << 
743                                                << 
744         if(Ecm>=3.9) {                         << 
745             sNNEta1 = sNNEta;                  << 
746         }                                      << 
747         else if (Ecm >= 3.5) {                 << 
748             sNNEta1 = -1916.2*Ecm*Ecm*Ecm + 21 << 
749         }                                      << 
750         else if (Ecm >= 2.525) {               << 
751             sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ec << 
752         }                                      << 
753         else {                                 << 
754             sNNEta1 = 17570.217219*Ecm*Ecm - 8 << 
755         }                                      << 
756                                                << 
757         sNNEta2 = -10220.89518466*Ecm*Ecm+5122 << 
758         if (sNNEta2 < 0.) sNNEta2=0.;          << 
759                                                << 
760         sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;   << 
761         if (sNNEta < 1.e-9 || Ecm < Thr0) sNNE << 
762                                                << 
763         return sNNEta/1000.; // parameterizati << 
764                                                << 
765     }                                             794     }
766                                                << 795     else {
767     G4double CrossSectionsMultiPionsAndResonan << 796      sNNOmega = NNToNNOmegaExcluIso(ener, 2);
768                                                << 
769         const G4double ener=KinematicsUtils::t << 
770         const G4int iso=ParticleTable::getIsos << 
771                                                << 
772         if (iso != 0) {                        << 
773             return NNToNNEtaExcluIso(ener, iso << 
774         }                                      << 
775         else {                                 << 
776             return 0.5*(NNToNNEtaExcluIso(ener << 
777         }                                      << 
778     }                                             797     }
779                                                << 798     
780                                                << 799     if (sNNOmega < 1.e-9) sNNOmega = 0.;
781     G4double CrossSectionsMultiPionsAndResonan << 800     
782                                                << 801     if (iso != 0) {
783         const G4double Ecm=0.001*ener;         << 802      return sNNOmega; 
784         G4double sNNOmega; // pp->pp+eta(+X)   << 
785         G4double sNNOmega1; // np->np+eta(+X)  << 
786         G4double x=Ecm*Ecm/7.06;               << 
787                                                << 
788         if(Ecm>4.0) {                          << 
789             sNNOmega = 2.5*std::pow(x-1, 1.47) << 
790         }                                      << 
791         else if(Ecm>2.802) { // 2802 MeV -> th << 
792             sNNOmega = (568.5254*Ecm*Ecm - 269 << 
793         if (sNNOmega <= NNToNNOmegaExcluIso(en << 
794         }                                      << 
795         else {                                 << 
796             sNNOmega = NNToNNOmegaExcluIso(ene << 
797         }                                      << 
798                                                << 
799         if (sNNOmega < 1.e-9) sNNOmega = 0.;   << 
800                                                << 
801         if (iso != 0) {                        << 
802         return sNNOmega;                       << 
803         }                                      << 
804                                                << 
805         sNNOmega1 = 3.*sNNOmega; // 3.0: ratio << 
806                                                << 
807         sNNOmega = 2.*sNNOmega1-sNNOmega;      << 
808                                                << 
809         if (sNNOmega < 1.e-9) sNNOmega = 0.;   << 
810                                                << 
811         return sNNOmega;                       << 
812     }                                             803     }
813                                                << 804     
814                                                << 805     sNNOmega1 = 3.*sNNOmega; // 3.0: ratio pn/pp (5 from theory; 2 from experiments)
815     G4double CrossSectionsMultiPionsAndResonan << 806         
816                                                << 807     sNNOmega = 2.*sNNOmega1-sNNOmega;
817         const G4double ener=KinematicsUtils::t << 808     
818         const G4int iso=ParticleTable::getIsos << 809     if (sNNOmega < 1.e-9) sNNOmega = 0.;
                                                   >> 810     
                                                   >> 811     return sNNOmega; 
                                                   >> 812    }
                                                   >> 813    
                                                   >> 814    
                                                   >> 815    G4double CrossSectionsMultiPionsAndResonances::NNToNNOmega(Particle const * const particle1, Particle const * const particle2) {
                                                   >> 816     
                                                   >> 817     const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
                                                   >> 818     const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
819 //jcd to be removed                               819 //jcd to be removed
820 //        return 0.;                           << 820 //    return 0.;
821 //jcd                                             821 //jcd    
822         if (iso != 0) {                        << 822     if (iso != 0) {
823             return NNToNNOmegaIso(ener, iso);  << 823      return NNToNNOmegaIso(ener, iso);
824         }                                      << 
825         else {                                 << 
826             return 0.5*(NNToNNOmegaIso(ener, 0 << 
827         }                                      << 
828     }                                             824     }
829                                                << 825     else {
830     G4double CrossSectionsMultiPionsAndResonan << 826      return 0.5*(NNToNNOmegaIso(ener, 0)+NNToNNOmegaIso(ener, 2));
831                                                << 
832         const G4double Ecm=0.001*ener;         << 
833         G4double sNNOmega; // pp->pp+eta       << 
834         G4double sNNOmega1; // np->np+eta      << 
835         G4double sthroot=std::sqrt(7.06);      << 
836                                                << 
837         if(Ecm>=3.0744) { // By hand (JCD)     << 
838             sNNOmega = 330.*(Ecm-sthroot)/(1.0 << 
839         }                                      << 
840         else if(Ecm>=2.65854){                 << 
841             sNNOmega = -1208.09757*std::pow(Ec << 
842         }                                      << 
843         else {                                 << 
844             sNNOmega = 0. ;                    << 
845         }                                      << 
846                                                << 
847         G4double Mn=ParticleTable::getRealMass << 
848         G4double Mp=ParticleTable::getRealMass << 
849         G4double Momega=ParticleTable::getReal << 
850         G4double Thr0=0.;                      << 
851         if (iso > 0) {                         << 
852             Thr0=2.*Mp+Momega;                 << 
853         }                                      << 
854         else if (iso < 0) {                    << 
855             Thr0=2.*Mn+Momega;                 << 
856         }                                      << 
857         else {                                 << 
858             Thr0=Mn+Mp+Momega;                 << 
859         }                                      << 
860                                                << 
861         if (sNNOmega < 1.e-9 || Ecm < Thr0) sN << 
862                                                << 
863         if (iso != 0) {                        << 
864             return sNNOmega/1000.; // paramete << 
865         }                                      << 
866                                                << 
867         sNNOmega1 = 3*sNNOmega; // 3.0: ratio  << 
868                                                << 
869         sNNOmega = 2*sNNOmega1-sNNOmega;       << 
870         if (sNNOmega < 1.e-9 || Ecm < Thr0) sN << 
871                                                << 
872         return sNNOmega/1000.; // parameteriza << 
873     }                                             827     }
874                                                << 828    }
875     G4double CrossSectionsMultiPionsAndResonan << 829    
                                                   >> 830    G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaExcluIso(const G4double ener, const G4int iso) {
                                                   >> 831         
                                                   >> 832         const G4double Ecm=0.001*ener;
                                                   >> 833         G4double sNNOmega; // pp->pp+eta
                                                   >> 834         G4double sNNOmega1; // np->np+eta
                                                   >> 835     G4double sthroot=std::sqrt(7.06);
                                                   >> 836         
                                                   >> 837         if(Ecm>=3.0744) { // By hand (JCD)
                                                   >> 838           sNNOmega = 330.*(Ecm-sthroot)/(1.05+std::pow((Ecm-sthroot),2));
                                                   >> 839         }
                                                   >> 840         else if(Ecm>=2.65854){
                                                   >> 841           sNNOmega = -1208.09757*std::pow(Ecm,3) + 10773.3322*std::pow(Ecm,2) - 31661.0223*Ecm + 30728.7241 ;
                                                   >> 842         }
                                                   >> 843         else {
                                                   >> 844           sNNOmega = 0. ;
                                                   >> 845         }
                                                   >> 846         
                                                   >> 847     G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
                                                   >> 848     G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
                                                   >> 849     G4double Momega=ParticleTable::getRealMass(Omega)/1000.;
                                                   >> 850     G4double Thr0=0.;
                                                   >> 851     if (iso > 0) {
                                                   >> 852      Thr0=2.*Mp+Momega;
                                                   >> 853     }
                                                   >> 854     else if (iso < 0) {
                                                   >> 855      Thr0=2.*Mn+Momega;
                                                   >> 856     }
                                                   >> 857     else {
                                                   >> 858      Thr0=Mn+Mp+Momega;
                                                   >> 859     }
                                                   >> 860     
                                                   >> 861     if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.;  // Thr0: Ecm threshold
                                                   >> 862         
                                                   >> 863         if (iso != 0) {
                                                   >> 864           return sNNOmega/1000.; // parameterization in microbarn (not millibarn)!
                                                   >> 865         }
                                                   >> 866         
                                                   >> 867     sNNOmega1 = 3*sNNOmega; // 3.0: ratio pn/pp
                                                   >> 868         
                                                   >> 869         sNNOmega = 2*sNNOmega1-sNNOmega;
                                                   >> 870         if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.;
                                                   >> 871         
                                                   >> 872         return sNNOmega/1000.; // parameterization in microbarn (not millibarn)!
                                                   >> 873       }
                                                   >> 874       
                                                   >> 875    G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaExclu(Particle const * const particle1, Particle const * const particle2) {
876 //jcd to be removed                               876 //jcd to be removed
877 //        return 0.;                           << 877 //    return 0.;
878 //jcd                                             878 //jcd    
879                                                << 879         
880         const G4double ener=KinematicsUtils::t << 880         const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
881         const G4int iso=ParticleTable::getIsos << 881         const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
882                                                << 882         
883         if (iso != 0) {                        << 883         if (iso != 0) {
884             return NNToNNOmegaExcluIso(ener, i << 884           return NNToNNOmegaExcluIso(ener, iso);
885         }                                      << 885         }
886         else {                                 << 886         else {
887             return 0.5*(NNToNNOmegaExcluIso(en << 887           return 0.5*(NNToNNOmegaExcluIso(ener, 0)+NNToNNOmegaExcluIso(ener, 2));
888         }                                      << 888         }
889     }                                          << 889       }
890                                                << 890    
891                                                << 891    
892     G4double CrossSectionsMultiPionsAndResonan << 892     G4double CrossSectionsMultiPionsAndResonances::NNToxPiNN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
893         //                                     << 893     //
894         //     Nucleon-Nucleon producing xpi p << 894     //     Nucleon-Nucleon producing xpi pions cross sections
895         //                                     << 895     //
896 // assert(xpi>0 && xpi<=nMaxPiNN);                896 // assert(xpi>0 && xpi<=nMaxPiNN);
897 // assert(particle1->isNucleon() && particle2-    897 // assert(particle1->isNucleon() && particle2->isNucleon());
898                                                << 898     
899         G4double oldXS1Pi=CrossSectionsMultiPi << 899     G4double oldXS1Pi=CrossSectionsMultiPions::NNToxPiNN(1,particle1, particle2);
900         G4double oldXS2Pi=CrossSectionsMultiPi << 900     G4double oldXS2Pi=CrossSectionsMultiPions::NNToxPiNN(2,particle1, particle2);
901         G4double oldXS3Pi=CrossSectionsMultiPi << 901     G4double oldXS3Pi=CrossSectionsMultiPions::NNToxPiNN(3,particle1, particle2);
902         G4double oldXS4Pi=CrossSectionsMultiPi << 902     G4double oldXS4Pi=CrossSectionsMultiPions::NNToxPiNN(4,particle1, particle2);
903         G4double xsEtaOmega=NNToNNEta(particle << 903     G4double xsEtaOmega=NNToNNEta(particle1, particle2)+NNToNNOmega(particle1, particle2);
904         G4double newXS1Pi=0.;                  << 904     G4double newXS1Pi=0.; 
905         G4double newXS2Pi=0.;                  << 905     G4double newXS2Pi=0.; 
906         G4double newXS3Pi=0.;                  << 906     G4double newXS3Pi=0.; 
907         G4double newXS4Pi=0.;                  << 907     G4double newXS4Pi=0.; 
908                                                << 908       
909         if (xpi == 1) {                        << 909     if (xpi == 1) {
910             if (oldXS4Pi != 0. || oldXS3Pi !=  << 910       if (oldXS4Pi != 0. || oldXS3Pi != 0.)
911                 newXS1Pi=oldXS1Pi;             << 911         newXS1Pi=oldXS1Pi;
912             else if (oldXS2Pi != 0.) {         << 912       else if (oldXS2Pi != 0.) {
913                 newXS2Pi=oldXS2Pi-xsEtaOmega;  << 913         newXS2Pi=oldXS2Pi-xsEtaOmega;
914                 if (newXS2Pi < 0.)             << 914        if (newXS2Pi < 0.)
915                     newXS1Pi=oldXS1Pi-(xsEtaOm << 915          newXS1Pi=oldXS1Pi-(xsEtaOmega-oldXS2Pi);
916                 else                           << 916         else
917                     newXS1Pi=oldXS1Pi;         << 917           newXS1Pi=oldXS1Pi;
918             }                                  << 918       }
919             else                               << 919       else 
920                 newXS1Pi=oldXS1Pi-xsEtaOmega;  << 920         newXS1Pi=oldXS1Pi-xsEtaOmega;
921             return newXS1Pi;                   << 921      return newXS1Pi;
922         }                                      << 922     }
923         else if (xpi == 2) {                   << 923     else if (xpi == 2) {
924             if (oldXS4Pi != 0.)                << 924       if (oldXS4Pi != 0.)
925             newXS2Pi=oldXS2Pi;                 << 925         newXS2Pi=oldXS2Pi;
926             else if (oldXS3Pi != 0.) {         << 926       else if (oldXS3Pi != 0.) {
927                 newXS3Pi=oldXS3Pi-xsEtaOmega;  << 927         newXS3Pi=oldXS3Pi-xsEtaOmega;
928                 if (newXS3Pi < 0.)             << 928         if (newXS3Pi < 0.)
929                     newXS2Pi=oldXS2Pi-(xsEtaOm << 929           newXS2Pi=oldXS2Pi-(xsEtaOmega-oldXS3Pi);
930                 else                           << 930         else
931                     newXS2Pi=oldXS2Pi;         << 931           newXS2Pi=oldXS2Pi;
932             }                                  << 932       }
933             else {                             << 933       else { 
934                 newXS2Pi=oldXS2Pi-xsEtaOmega;  << 934         newXS2Pi=oldXS2Pi-xsEtaOmega;
935                 if (newXS2Pi < 0.)             << 935         if (newXS2Pi < 0.)
936                     newXS2Pi=0.;               << 936           newXS2Pi=0.;
937             }                                  << 937       }
938             return newXS2Pi;                   << 938       return newXS2Pi;
939         }                                      << 939     }
940         else if (xpi == 3) {                   << 940     else if (xpi == 3) {
941             if (oldXS4Pi != 0.) {              << 941       if (oldXS4Pi != 0.) {
942                 newXS4Pi=oldXS4Pi-xsEtaOmega;  << 942         newXS4Pi=oldXS4Pi-xsEtaOmega;
943                 if (newXS4Pi < 0.)             << 943        if (newXS4Pi < 0.)
944                     newXS3Pi=oldXS3Pi-(xsEtaOm << 944         newXS3Pi=oldXS3Pi-(xsEtaOmega-oldXS4Pi);
945                 else                           << 945        else
946                     newXS3Pi=oldXS3Pi;         << 946         newXS3Pi=oldXS3Pi;
947             }                                  << 947       }
948             else {                             << 948       else { 
949                 newXS3Pi=oldXS3Pi-xsEtaOmega;  << 949         newXS3Pi=oldXS3Pi-xsEtaOmega;       
950                 if (newXS3Pi < 0.)             << 950         if (newXS3Pi < 0.)
951                     newXS3Pi=0.;               << 951           newXS3Pi=0.;
952             }                                  << 952       }
953             return newXS3Pi;                   << 953       return newXS3Pi;
954         }                                      << 954     }
955         else if (xpi == 4) {                   << 955     else if (xpi == 4) {
956             newXS4Pi=oldXS4Pi-xsEtaOmega;      << 956       newXS4Pi=oldXS4Pi-xsEtaOmega;
957             if (newXS4Pi < 0.)                 << 957       if (newXS4Pi < 0.)
958                 newXS4Pi=0.;                   << 958         newXS4Pi=0.;
959             return newXS4Pi;                   << 959     return newXS4Pi;
960         }                                      << 960     }
961                                                << 961 
962         else // should never reach this point  << 962     else // should never reach this point
963         return 0.;                             << 963       return 0.;
964     }                                          << 964   }
965                                                << 965 
966                                                << 966 
967     G4double CrossSectionsMultiPionsAndResonan << 967       G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaOnePi(Particle const * const particle1, Particle const * const particle2) {
968         // Cross section for nucleon-nucleon p << 968         // Cross section for nucleon-nucleon producing one eta and one pion
969                                                << 969         
970         const G4int iso=ParticleTable::getIsos << 970         const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
971         if (iso!=0)                            << 971         if (iso!=0)
972             return 0.;                         << 972           return 0.;
973                                                << 973         
974         const G4double ener=KinematicsUtils::t << 974         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)
975         if (ener < 2018.563) return 0.;        << 975         if (ener < 2018.563) return 0.;
976                                                << 976         
977         const G4double xsiso2=CrossSectionsMul << 977         const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
978         const G4double xsiso0=CrossSectionsMul << 978         const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
979                                                << 979 
980         return 0.25*(CrossSectionsMultiPions:: << 980         return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
981     }                                          << 981       }
982                                                << 982       
983     G4double CrossSectionsMultiPionsAndResonan << 983       G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
984         const G4double ener=KinematicsUtils::t << 984         const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
985         if (ener < 2018.563) return 0.;        << 985         if (ener < 2018.563) return 0.;
986         const G4int iso=ParticleTable::getIsos << 986         const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
987                                                << 987         
988         const G4double xsiso2=CrossSectionsMul << 988         const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
989         if (iso != 0)                          << 989         if (iso != 0)
990             return CrossSectionsMultiPions::NN << 990           return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
991         else {                                 << 991         else {
992             const G4double xsiso0=CrossSection << 992           const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
993             return 0.5*(CrossSectionsMultiPion << 993           return 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
994         }                                      << 994         }
995     }                                          << 995       }
996                                                << 996       
997     G4double CrossSectionsMultiPionsAndResonan << 997       G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaTwoPi(Particle const * const particle1, Particle const * const particle2) {
998         //                                     << 998         //
999         //     Nucleon-Nucleon producing one e << 999         //     Nucleon-Nucleon producing one eta and two pions 
1000         //                                    << 1000         //
1001         const G4double ener=KinematicsUtils:: << 1001         const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1002         if (ener < 2018.563) return 0.;       << 1002         if (ener < 2018.563) return 0.;
1003         const G4int iso=ParticleTable::getIso << 1003         const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1004                                               << 1004         
1005                                               << 1005         
1006         const G4double xsiso2=CrossSectionsMu << 1006         const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1007         if (iso != 0) {                       << 1007         if (iso != 0) {
1008             return CrossSectionsMultiPions::N << 1008           return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1009         }                                     << 1009         }
1010         else {                                << 1010         else {
1011             const G4double xsiso0=CrossSectio << 1011           const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1012             return 0.5*(CrossSectionsMultiPio << 1012           return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
1013         }                                     << 1013         }
1014     }                                         << 1014       }
1015                                               << 1015       
1016     G4double CrossSectionsMultiPionsAndResona << 1016       G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaThreePi(Particle const * const particle1, Particle const * const particle2) {
1017         //                                    << 1017         //
1018         //     Nucleon-Nucleon producing one  << 1018         //     Nucleon-Nucleon producing one eta and three pions
1019         //                                    << 1019         //
1020                                               << 1020         
1021         const G4double ener=KinematicsUtils:: << 1021         const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1022         if (ener < 2018.563) return 0.;       << 1022         if (ener < 2018.563) return 0.;
1023         const G4int iso=ParticleTable::getIso << 1023         const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1024                                               << 1024         
1025                                               << 1025         
1026         const G4double xsiso2=CrossSectionsMu << 1026         const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1027         const G4double xs1pi2=CrossSectionsMu << 1027         const G4double xs1pi2=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2);
1028         const G4double xs2pi2=CrossSectionsMu << 1028         const G4double xs2pi2=CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1029         if (iso != 0)                         << 1029         if (iso != 0)
1030             return CrossSectionsMultiPions::N << 1030           return CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
1031         else {                                << 1031         else {
1032             const G4double xsiso0=CrossSectio << 1032           const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1033             const G4double xs1pi0=CrossSectio << 1033           const G4double xs1pi0=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0);
1034             const G4double xs2pi0=CrossSectio << 1034           const G4double xs2pi0=CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0);
1035             return 0.5*(CrossSectionsMultiPio << 1035           return 0.5*(CrossSectionsMultiPions::NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
1036         }                                     << 1036         }
1037     }                                         << 1037       }
1038                                               << 1038       
1039     G4double CrossSectionsMultiPionsAndResona << 1039       G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaFourPi(Particle const * const particle1, Particle const * const particle2) {
1040         //                                    << 1040         //
1041         //     Nucleon-Nucleon producing one  << 1041         //     Nucleon-Nucleon producing one eta and four pions
1042         //                                    << 1042         //
1043                                               << 1043         
1044         const G4double ener=KinematicsUtils:: << 1044         const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1045         if (ener < 2018.563) return 0.;       << 1045         if (ener < 2018.563) return 0.;
1046         const G4double s = ener*ener;         << 1046         const G4double s = ener*ener;
1047         const G4int i = ParticleTable::getIso << 1047         const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1048         G4double xsinelas;                    << 1048     G4double xsinelas;
1049         if (i!=0)                             << 1049         if (i!=0) 
1050             xsinelas=CrossSectionsMultiPions: << 1050      xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1051         else                                  << 1051     else 
1052             xsinelas=0.5*(CrossSectionsMultiP << 1052           xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1053         if (xsinelas <= 1.e-9) return 0.;     << 1053         if (xsinelas <= 1.e-9) return 0.;
1054         G4double ratio=(NNToNNEta(particle1,  << 1054        G4double ratio=(NNToNNEta(particle1, particle2)-NNToNNEtaExclu(particle1, particle2))/xsinelas;
1055         if(s<6.25E6)                          << 1055         if(s<6.25E6)
1056             return 0.;                        << 1056           return 0.;
1057         const G4double sigma = NNToNNEta(part << 1057         const G4double sigma = NNToNNEta(particle1, particle2) - NNToNNEtaExclu(particle1, particle2) - ratio*(NNToNNEtaOnePiOrDelta(particle1, particle2) + NNToNNEtaTwoPi(particle1, particle2) + NNToNNEtaThreePi(particle1, particle2));
1058         return ((sigma>1.e-9) ? sigma : 0.);  << 1058         return ((sigma>1.e-9) ? sigma : 0.);
1059     }                                         << 1059       }
1060                                               << 1060       
1061     G4double CrossSectionsMultiPionsAndResona << 1061       G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaxPi(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
1062         //                                    << 1062         //
1063         //     Nucleon-Nucleon producing one  << 1063         //     Nucleon-Nucleon producing one eta and xpi pions
1064         //                                    << 1064         //
1065 // assert(xpi>0 && xpi<=nMaxPiNN);               1065 // assert(xpi>0 && xpi<=nMaxPiNN);
1066 // assert(particle1->isNucleon() && particle2    1066 // assert(particle1->isNucleon() && particle2->isNucleon());
                                                   >> 1067         
                                                   >> 1068         const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
                                                   >> 1069         if (ener < 2018.563) return 0.;
                                                   >> 1070         const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
                                                   >> 1071     G4double xsinelas;
                                                   >> 1072         if (i!=0) 
                                                   >> 1073      xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
                                                   >> 1074     else 
                                                   >> 1075           xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
                                                   >> 1076         if (xsinelas <= 1.e-9) return 0.;
                                                   >> 1077        G4double ratio=(NNToNNEta(particle1, particle2)-NNToNNEtaExclu(particle1, particle2))/xsinelas;
                                                   >> 1078 
                                                   >> 1079         if (xpi == 1)
                                                   >> 1080           return NNToNNEtaOnePi(particle1, particle2)*ratio;
                                                   >> 1081         else if (xpi == 2)
                                                   >> 1082           return NNToNNEtaTwoPi(particle1, particle2)*ratio;
                                                   >> 1083         else if (xpi == 3)
                                                   >> 1084           return NNToNNEtaThreePi(particle1, particle2)*ratio;
                                                   >> 1085         else if (xpi == 4)
                                                   >> 1086           return NNToNNEtaFourPi(particle1, particle2);
                                                   >> 1087         else // should never reach this point
                                                   >> 1088           return 0.;
                                                   >> 1089       }
1067                                                  1090 
1068         const G4double ener=KinematicsUtils:: << 1091       
1069         if (ener < 2018.563) return 0.;       << 1092       G4double CrossSectionsMultiPionsAndResonances::NNToNDeltaEta(Particle const * const p1, Particle const * const p2) {
1070         const G4int i = ParticleTable::getIso << 
1071         G4double xsinelas;                    << 
1072         if (i!=0)                             << 
1073             xsinelas=CrossSectionsMultiPions: << 
1074         else                                  << 
1075             xsinelas=0.5*(CrossSectionsMultiP << 
1076         if (xsinelas <= 1.e-9) return 0.;     << 
1077         G4double ratio=(NNToNNEta(particle1,  << 
1078                                               << 
1079         if (xpi == 1)                         << 
1080             return NNToNNEtaOnePi(particle1,  << 
1081         else if (xpi == 2)                    << 
1082             return NNToNNEtaTwoPi(particle1,  << 
1083         else if (xpi == 3)                    << 
1084             return NNToNNEtaThreePi(particle1 << 
1085         else if (xpi == 4)                    << 
1086             return NNToNNEtaFourPi(particle1, << 
1087         else // should never reach this point << 
1088             return 0.;                        << 
1089     }                                         << 
1090                                               << 
1091                                               << 
1092     G4double CrossSectionsMultiPionsAndResona << 
1093 // assert(p1->isNucleon() && p2->isNucleon())    1093 // assert(p1->isNucleon() && p2->isNucleon());
1094         const G4int i = ParticleTable::getIso << 1094     const G4int i = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1095         const G4double ener=KinematicsUtils:: << 1095         const G4double ener=KinematicsUtils::totalEnergyInCM(p1, p2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1096         if (ener < 2018.563) return 0.;       << 1096         if (ener < 2018.563) return 0.;
1097         G4double xsinelas;                    << 1097     G4double xsinelas;
1098         if (i!=0)                             << 1098         if (i!=0) 
1099             xsinelas=CrossSectionsMultiPions: << 1099      xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1100         else                                  << 1100     else 
1101             xsinelas=0.5*(CrossSectionsMultiP << 1101           xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1102         if (xsinelas <= 1.e-9) return 0.;     << 1102         if (xsinelas <= 1.e-9) return 0.;
1103         G4double ratio=(NNToNNEta(p1, p2)-NNT << 1103        G4double ratio=(NNToNNEta(p1, p2)-NNToNNEtaExclu(p1, p2))/xsinelas;
1104         G4double sigma = NNToNNEtaOnePiOrDelt << 1104     G4double sigma = NNToNNEtaOnePiOrDelta(p1, p2)*ratio;
1105         if(i==0)                              << 1105     if(i==0)
1106             sigma *= 0.5;                     << 1106           sigma *= 0.5;
1107         return sigma;                         << 1107     return sigma;
1108     }                                         << 1108       }
1109                                               << 
1110                                               << 
1111     G4double CrossSectionsMultiPionsAndResona << 
1112         // Cross section for nucleon-nucleon  << 
1113                                               << 
1114         const G4int iso=ParticleTable::getIso << 
1115         if (iso!=0)                           << 
1116             return 0.;                        << 
1117                                               << 
1118         const G4double ener=KinematicsUtils:: << 
1119         if (ener < 2018.563) return 0.;       << 
1120                                               << 
1121         const G4double xsiso2=CrossSectionsMu << 
1122         const G4double xsiso0=CrossSectionsMu << 
1123                                                  1109 
1124         return 0.25*(CrossSectionsMultiPions: << 1110   
1125     }                                         << 1111   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaOnePi(Particle const * const particle1, Particle const * const particle2) {
1126                                               << 1112    // Cross section for nucleon-nucleon producing one omega and one pion
1127     G4double CrossSectionsMultiPionsAndResona << 1113    
1128         const G4double ener=KinematicsUtils:: << 1114    const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1129         if (ener < 2018.563) return 0.;       << 1115    if (iso!=0)
1130         const G4int iso=ParticleTable::getIso << 1116     return 0.;
1131                                               << 1117    
1132         const G4double xsiso2=CrossSectionsMu << 1118    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)
1133         if (iso != 0)                         << 1119    if (ener < 2018.563) return 0.;
1134             return CrossSectionsMultiPions::N << 1120    
1135         else {                                << 1121    const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1136             const G4double xsiso0=CrossSectio << 1122    const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1137             return 0.5*(CrossSectionsMultiPio << 1123    
1138         }                                     << 1124    return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
1139     }                                         << 1125   }
1140                                               << 1126   
1141     G4double CrossSectionsMultiPionsAndResona << 1127   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
1142         //                                    << 1128    const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1143         //     Nucleon-Nucleon producing one  << 1129    if (ener < 2018.563) return 0.;
1144         //                                    << 1130    const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1145         const G4double ener=KinematicsUtils:: << 1131    
1146         if (ener < 2018.563) return 0.;       << 1132    const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1147         const G4int iso=ParticleTable::getIso << 1133    if (iso != 0)
1148                                               << 1134     return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
1149         const G4double xsiso2=CrossSectionsMu << 1135    else {
1150         if (iso != 0) {                       << 1136     const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1151             return CrossSectionsMultiPions::N << 1137     return 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
1152         }                                     << 1138    }
1153         else {                                << 1139   }
1154             const G4double xsiso0=CrossSectio << 1140   
1155             return 0.5*(CrossSectionsMultiPio << 1141   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaTwoPi(Particle const * const particle1, Particle const * const particle2) {
1156         }                                     << 1142    //
1157     }                                         << 1143    //     Nucleon-Nucleon producing one omega and two pions 
1158                                               << 1144    //
1159     G4double CrossSectionsMultiPionsAndResona << 1145    const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1160         //                                    << 1146    if (ener < 2018.563) return 0.;
1161         //     Nucleon-Nucleon producing one  << 1147    const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1162         //                                    << 1148    
1163                                               << 1149    
1164         const G4double ener=KinematicsUtils:: << 1150    const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1165         if (ener < 2018.563) return 0.;       << 1151    if (iso != 0) {
1166         const G4int iso=ParticleTable::getIso << 1152     return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1167                                               << 1153    }
1168                                               << 1154    else {
1169         const G4double xsiso2=CrossSectionsMu << 1155     const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1170         const G4double xs1pi2=CrossSectionsMu << 1156     return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
1171         const G4double xs2pi2=CrossSectionsMu << 1157    }
1172         if (iso != 0)                         << 1158   }
1173             return CrossSectionsMultiPions::N << 1159   
1174         else {                                << 1160   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaThreePi(Particle const * const particle1, Particle const * const particle2) {
1175             const G4double xsiso0=CrossSectio << 1161    //
1176             const G4double xs1pi0=CrossSectio << 1162    //     Nucleon-Nucleon producing one omega and three pions
1177             const G4double xs2pi0=CrossSectio << 1163    //
1178             return 0.5*(CrossSectionsMultiPio << 1164    
1179         }                                     << 1165    const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1180     }                                         << 1166    if (ener < 2018.563) return 0.;
1181                                               << 1167    const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1182     G4double CrossSectionsMultiPionsAndResona << 1168    
1183         //                                    << 1169    
1184         //     Nucleon-Nucleon producing one  << 1170    const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1185         //                                    << 1171    const G4double xs1pi2=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2);
                                                   >> 1172    const G4double xs2pi2=CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
                                                   >> 1173    if (iso != 0)
                                                   >> 1174     return CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
                                                   >> 1175    else {
                                                   >> 1176     const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
                                                   >> 1177     const G4double xs1pi0=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0);
                                                   >> 1178     const G4double xs2pi0=CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0);
                                                   >> 1179     return 0.5*(CrossSectionsMultiPions::NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
                                                   >> 1180    }
                                                   >> 1181   }
                                                   >> 1182   
                                                   >> 1183   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaFourPi(Particle const * const particle1, Particle const * const particle2) {
                                                   >> 1184    //
                                                   >> 1185    //     Nucleon-Nucleon producing one omega and four pions
                                                   >> 1186    //
1186 //jcd to be removed                              1187 //jcd to be removed
1187 //        return 0.;                          << 1188 //   return 0.;
1188 //jcd                                            1189 //jcd    
1189                                               << 1190    
1190         const G4double ener=KinematicsUtils:: << 1191    const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1191         if (ener < 2018.563) return 0.;       << 1192    if (ener < 2018.563) return 0.;
1192         const G4double s = ener*ener;         << 1193    const G4double s = ener*ener;
1193         const G4int i = ParticleTable::getIso << 1194    const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1194         G4double xsinelas;                    << 1195    G4double xsinelas;
1195         if (i!=0)                             << 1196    if (i!=0) 
1196             xsinelas=CrossSectionsMultiPions: << 1197     xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1197         else                                  << 1198    else 
1198             xsinelas=0.5*(CrossSectionsMultiP << 1199     xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1199         if (xsinelas <= 1.e-9) return 0.;     << 1200    if (xsinelas <= 1.e-9) return 0.;
1200         G4double ratio=(NNToNNOmega(particle1 << 1201    G4double ratio=(NNToNNOmega(particle1, particle2)-NNToNNOmegaExclu(particle1, particle2))/xsinelas;
1201         if(s<6.25E6)                          << 1202    if(s<6.25E6)
1202             return 0.;                        << 1203     return 0.;
1203         const G4double sigma = NNToNNOmega(pa << 1204    const G4double sigma = NNToNNOmega(particle1, particle2) - NNToNNOmegaExclu(particle1, particle2) - ratio*(NNToNNOmegaOnePiOrDelta(particle1, particle2) + NNToNNOmegaTwoPi(particle1, particle2) + NNToNNOmegaThreePi(particle1, particle2));
1204         return ((sigma>1.e-9) ? sigma : 0.);  << 1205    return ((sigma>1.e-9) ? sigma : 0.);
1205     }                                         << 1206   }
1206                                               << 1207   
1207     G4double CrossSectionsMultiPionsAndResona << 1208   G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaxPi(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
1208         //                                    << 1209    //
1209         //     Nucleon-Nucleon producing one  << 1210    //     Nucleon-Nucleon producing one omega and xpi pions
1210         //                                    << 1211    //
1211 // assert(xpi>0 && xpi<=nMaxPiNN);               1212 // assert(xpi>0 && xpi<=nMaxPiNN);
1212 // assert(particle1->isNucleon() && particle2    1213 // assert(particle1->isNucleon() && particle2->isNucleon());
1213 //jcd to be removed                              1214 //jcd to be removed
1214 //        return 0.;                          << 1215 //   return 0.;
1215 //jcd                                            1216 //jcd    
1216                                               << 1217    
1217         const G4double ener=KinematicsUtils:: << 1218    const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1218         if (ener < 2018.563) return 0.;       << 1219    if (ener < 2018.563) return 0.;
1219         const G4int i = ParticleTable::getIso << 1220    const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1220         G4double xsinelas;                    << 1221    G4double xsinelas;
1221         if (i!=0)                             << 1222    if (i!=0) 
1222             xsinelas=CrossSectionsMultiPions: << 1223     xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1223         else                                  << 1224    else 
1224             xsinelas=0.5*(CrossSectionsMultiP << 1225     xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1225         if (xsinelas <= 1.e-9) return 0.;     << 1226    if (xsinelas <= 1.e-9) return 0.;
1226         G4double ratio=(NNToNNOmega(particle1 << 1227    G4double ratio=(NNToNNOmega(particle1, particle2)-NNToNNOmegaExclu(particle1, particle2))/xsinelas;
1227                                               << 1228    
1228         if (xpi == 1)                         << 1229    if (xpi == 1)
1229             return NNToNNOmegaOnePi(particle1 << 1230     return NNToNNOmegaOnePi(particle1, particle2)*ratio;
1230         else if (xpi == 2)                    << 1231    else if (xpi == 2)
1231             return NNToNNOmegaTwoPi(particle1 << 1232     return NNToNNOmegaTwoPi(particle1, particle2)*ratio;
1232         else if (xpi == 3)                    << 1233    else if (xpi == 3)
1233             return NNToNNOmegaThreePi(particl << 1234     return NNToNNOmegaThreePi(particle1, particle2)*ratio;
1234         else if (xpi == 4)                    << 1235    else if (xpi == 4)
1235             return NNToNNOmegaFourPi(particle << 1236     return NNToNNOmegaFourPi(particle1, particle2);
1236         else // should never reach this point << 1237    else // should never reach this point
1237             return 0.;                        << 1238     return 0.;
1238     }                                         << 1239   }
1239                                               << 1240   
1240                                               << 1241   
1241     G4double CrossSectionsMultiPionsAndResona << 1242   G4double CrossSectionsMultiPionsAndResonances::NNToNDeltaOmega(Particle const * const p1, Particle const * const p2) {
1242 // assert(p1->isNucleon() && p2->isNucleon())    1243 // assert(p1->isNucleon() && p2->isNucleon());
1243 //jcd to be removed                              1244 //jcd to be removed
1244 //        return 0.;                          << 1245 //   return 0.;
1245 //jcd                                            1246 //jcd    
1246         const G4int i = ParticleTable::getIso << 1247    const G4int i = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1247         const G4double ener=KinematicsUtils:: << 1248    const G4double ener=KinematicsUtils::totalEnergyInCM(p1, p2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1248         if (ener < 2018.563) return 0.;       << 1249    if (ener < 2018.563) return 0.;
1249         G4double xsinelas;                    << 1250    G4double xsinelas;
1250         if (i!=0)                             << 1251    if (i!=0) 
1251             xsinelas=CrossSectionsMultiPions: << 1252     xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1252         else                                  << 1253    else 
1253             xsinelas=0.5*(CrossSectionsMultiP << 1254     xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1254         if (xsinelas <= 1.e-9) return 0.;     << 1255    if (xsinelas <= 1.e-9) return 0.;
1255         G4double ratio=(NNToNNOmega(p1, p2)-N << 1256    G4double ratio=(NNToNNOmega(p1, p2)-NNToNNOmegaExclu(p1, p2))/xsinelas;
1256         G4double sigma = NNToNNOmegaOnePiOrDe << 1257    G4double sigma = NNToNNOmegaOnePiOrDelta(p1, p2)*ratio;
1257         if(i==0)                              << 1258    if(i==0)
1258             sigma *= 0.5;                     << 1259     sigma *= 0.5;
1259         return sigma;                         << 1260    return sigma;
1260     }                                         << 1261   }
1261                                               << 1262   
1262                                               << 1263   
1263 } // namespace G4INCL                            1264 } // namespace G4INCL
1264                                                  1265 
1265                                                  1266