Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/utils/src/G4INCLParticleTable.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/utils/src/G4INCLParticleTable.cc (Version 11.3.0) and /processes/hadronic/models/inclxx/utils/src/G4INCLParticleTable.cc (Version 10.3.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 #include "G4INCLParticleTable.hh"                  38 #include "G4INCLParticleTable.hh"
 39 #include "G4INCLNuclearMassTable.hh"               39 #include "G4INCLNuclearMassTable.hh"
 40 #include <algorithm>                               40 #include <algorithm>
 41 // #include <cassert>                              41 // #include <cassert>
 42 #include <cmath>                                   42 #include <cmath>
 43 #include <cctype>                                  43 #include <cctype>
 44 #include <sstream>                                 44 #include <sstream>
 45 #ifdef INCLXX_IN_GEANT4_MODE                       45 #ifdef INCLXX_IN_GEANT4_MODE
 46 #include "G4SystemOfUnits.hh"                      46 #include "G4SystemOfUnits.hh"
 47 #endif                                             47 #endif
 48                                                    48 
 49 #ifdef INCLXX_IN_GEANT4_MODE                       49 #ifdef INCLXX_IN_GEANT4_MODE
 50 #include "G4PhysicalConstants.hh"                  50 #include "G4PhysicalConstants.hh"
 51 #include "G4SystemOfUnits.hh"                      51 #include "G4SystemOfUnits.hh"
 52 #endif                                             52 #endif
 53                                                    53 
 54 namespace G4INCL {                                 54 namespace G4INCL {
 55                                                    55 
 56   namespace ParticleTable {                        56   namespace ParticleTable {
 57                                                    57 
 58     namespace {                                    58     namespace {
 59                                                    59 
 60       /// \brief Static instance of the Natura     60       /// \brief Static instance of the NaturalIsotopicAbundances class
 61       const NaturalIsotopicDistributions *theN     61       const NaturalIsotopicDistributions *theNaturalIsotopicDistributions = NULL;
 62                                                    62 
 63       const G4double theINCLNucleonMass = 938.     63       const G4double theINCLNucleonMass = 938.2796;
 64       const G4double theINCLPionMass = 138.0;      64       const G4double theINCLPionMass = 138.0;
 65       const G4double theINCLLambdaMass = 1115. << 
 66       //const G4double theINCLSigmaMass = 1197 << 
 67       //const G4double theINCLKaonMass = 497.6 << 
 68       const G4double theINCLEtaMass = 547.862;     65       const G4double theINCLEtaMass = 547.862;
 69       const G4double theINCLOmegaMass = 782.65     66       const G4double theINCLOmegaMass = 782.65;
 70       const G4double theINCLEtaPrimeMass = 957     67       const G4double theINCLEtaPrimeMass = 957.78;
 71       const G4double theINCLPhotonMass = 0.0;      68       const G4double theINCLPhotonMass = 0.0;
 72       G4ThreadLocal G4double protonMass = 0.0;     69       G4ThreadLocal G4double protonMass = 0.0;
 73       G4ThreadLocal G4double neutronMass = 0.0     70       G4ThreadLocal G4double neutronMass = 0.0;
 74       G4ThreadLocal G4double piPlusMass = 0.0;     71       G4ThreadLocal G4double piPlusMass = 0.0;
 75       G4ThreadLocal G4double piMinusMass = 0.0     72       G4ThreadLocal G4double piMinusMass = 0.0;
 76       G4ThreadLocal G4double piZeroMass = 0.0;     73       G4ThreadLocal G4double piZeroMass = 0.0;
 77       G4ThreadLocal G4double SigmaPlusMass = 0 << 
 78       G4ThreadLocal G4double SigmaZeroMass = 0 << 
 79       G4ThreadLocal G4double SigmaMinusMass =  << 
 80       G4ThreadLocal G4double LambdaMass = 0.0; << 
 81       G4ThreadLocal G4double XiMinusMass = 0.0 << 
 82       G4ThreadLocal G4double XiZeroMass = 0.0; << 
 83       G4ThreadLocal G4double antiProtonMass =  << 
 84       G4ThreadLocal G4double antiNeutronMass = << 
 85       G4ThreadLocal G4double antiSigmaPlusMass << 
 86       G4ThreadLocal G4double antiSigmaZeroMass << 
 87       G4ThreadLocal G4double antiSigmaMinusMas << 
 88       G4ThreadLocal G4double antiLambdaMass =  << 
 89       G4ThreadLocal G4double antiXiMinusMass = << 
 90       G4ThreadLocal G4double antiXiZeroMass =  << 
 91       G4ThreadLocal G4double KPlusMass = 0.0;  << 
 92       G4ThreadLocal G4double KZeroMass = 0.0;  << 
 93       G4ThreadLocal G4double KZeroBarMass = 0. << 
 94       G4ThreadLocal G4double KShortMass = 0.0; << 
 95       G4ThreadLocal G4double KLongMass = 0.0;  << 
 96       G4ThreadLocal G4double KMinusMass = 0.0; << 
 97       G4ThreadLocal G4double etaMass = 0.0;        74       G4ThreadLocal G4double etaMass = 0.0;
 98       G4ThreadLocal G4double omegaMass = 0.0;      75       G4ThreadLocal G4double omegaMass = 0.0;
 99       G4ThreadLocal G4double etaPrimeMass = 0.     76       G4ThreadLocal G4double etaPrimeMass = 0.0;
100       G4ThreadLocal G4double photonMass = 0.0;     77       G4ThreadLocal G4double photonMass = 0.0;
101                                                    78 
102       // Hard-coded values of the real particl     79       // Hard-coded values of the real particle masses (MeV/c^2)
103       G4ThreadLocal G4double theRealProtonMass     80       G4ThreadLocal G4double theRealProtonMass = 938.27203;
104       G4ThreadLocal G4double theRealNeutronMas     81       G4ThreadLocal G4double theRealNeutronMass = 939.56536;
105       G4ThreadLocal G4double theRealChargedPiM     82       G4ThreadLocal G4double theRealChargedPiMass = 139.57018;
106       G4ThreadLocal G4double theRealPiZeroMass     83       G4ThreadLocal G4double theRealPiZeroMass = 134.9766;
107       G4ThreadLocal G4double theRealLambdaMass << 
108       G4ThreadLocal G4double theRealSigmaPlusM << 
109       G4ThreadLocal G4double theRealSigmaZeroM << 
110       G4ThreadLocal G4double theRealSigmaMinus << 
111       G4ThreadLocal G4double theRealAntiProton << 
112       G4ThreadLocal G4double theRealXiMinusMas << 
113       G4ThreadLocal G4double theRealXiZeroMass << 
114       G4ThreadLocal G4double theRealAntiNeutro << 
115       G4ThreadLocal G4double theRealAntiLambda << 
116       G4ThreadLocal G4double theRealAntiSigmaP << 
117       G4ThreadLocal G4double theRealAntiSigmaZ << 
118       G4ThreadLocal G4double theRealAntiSigmaM << 
119       G4ThreadLocal G4double theRealAntiXiMinu << 
120       G4ThreadLocal G4double theRealAntiXiZero << 
121       G4ThreadLocal G4double theRealChargedKao << 
122       G4ThreadLocal G4double theRealNeutralKao << 
123       G4ThreadLocal G4double theRealEtaMass =      84       G4ThreadLocal G4double theRealEtaMass = 547.862;
124       G4ThreadLocal G4double theRealOmegaMass      85       G4ThreadLocal G4double theRealOmegaMass = 782.65;
125       G4ThreadLocal G4double theRealEtaPrimeMa     86       G4ThreadLocal G4double theRealEtaPrimeMass = 957.78;
126       G4ThreadLocal G4double theRealPhotonMass     87       G4ThreadLocal G4double theRealPhotonMass = 0.0;
127                                                << 
128       // Width (second)                        << 
129       const G4double theChargedPiWidth = 2.603 << 
130       const G4double thePiZeroWidth = 8.52e-17 << 
131       const G4double theEtaWidth = 5.025e-19;  << 
132       const G4double theOmegaWidth = 7.7528e-2 << 
133       const G4double theEtaPrimeWidth = 3.3243 << 
134       const G4double theChargedKaonWidth = 1.2 << 
135       const G4double theKShortWidth = 8.954e-1 << 
136       const G4double theKLongWidth = 5.116e-08 << 
137       const G4double theLambdaWidth = 2.632e-1 << 
138       const G4double theSigmaPlusWidth = 8.018 << 
139       const G4double theSigmaZeroWidth = 7.4e- << 
140       const G4double theSigmaMinusWidth = 1.47 << 
141       //const G4double theXiMinusWidth = 1.639 << 
142       //const G4double theXiZeroWidth = 2.90e- << 
143       //const G4double theAntiLambdaWidth = 2. << 
144       //const G4double theAntiSigmaPlusWidth = << 
145       //const G4double theAntiSigmaZeroWidth = << 
146       //const G4double theAntiSigmaMinusWidth  << 
147       //const G4double theAntiXiMinusWidth = 1 << 
148       //const G4double theAntiXiZeroWidth = 2. << 
149       G4ThreadLocal G4double piPlusWidth = 0.0 << 
150       G4ThreadLocal G4double piMinusWidth = 0. << 
151       G4ThreadLocal G4double piZeroWidth = 0.0 << 
152       G4ThreadLocal G4double etaWidth = 0.0;   << 
153       G4ThreadLocal G4double omegaWidth = 0.0; << 
154       G4ThreadLocal G4double etaPrimeWidth = 0 << 
155       G4ThreadLocal G4double LambdaWidth = 0.0 << 
156       G4ThreadLocal G4double SigmaPlusWidth =  << 
157       G4ThreadLocal G4double SigmaZeroWidth =  << 
158       G4ThreadLocal G4double SigmaMinusWidth = << 
159       G4ThreadLocal G4double KPlusWidth = 0.0; << 
160       G4ThreadLocal G4double KMinusWidth = 0.0 << 
161       G4ThreadLocal G4double KShortWidth = 0.0 << 
162       G4ThreadLocal G4double KLongWidth = 0.0; << 
163       G4ThreadLocal G4double XiMinusWidth = 0. << 
164       G4ThreadLocal G4double XiZeroWidth = 0.0 << 
165       G4ThreadLocal G4double antiLambdaWidth = << 
166       G4ThreadLocal G4double antiSigmaZeroWidt << 
167       G4ThreadLocal G4double antiSigmaMinusWid << 
168       G4ThreadLocal G4double antiSigmaPlusWidt << 
169       G4ThreadLocal G4double antiXiZeroWidth = << 
170       G4ThreadLocal G4double antiXiMinusWidth  << 
171                                                    88 
                                                   >>  89     // Width (second)
                                                   >>  90     const G4double theChargedPiWidth = 2.6033e-08;
                                                   >>  91     const G4double thePiZeroWidth = 8.52e-17;
                                                   >>  92     const G4double theEtaWidth = 5.025e-19; // 1.31 keV
                                                   >>  93     const G4double theOmegaWidth = 7.7528e-23; // 8.49 MeV
                                                   >>  94     const G4double theEtaPrimeWidth = 3.3243e-21; // 0.198 MeV
                                                   >>  95     G4ThreadLocal G4double piPlusWidth = 0.0;
                                                   >>  96     G4ThreadLocal G4double piMinusWidth = 0.0;
                                                   >>  97     G4ThreadLocal G4double piZeroWidth = 0.0;
                                                   >>  98     G4ThreadLocal G4double etaWidth = 0.0;
                                                   >>  99     G4ThreadLocal G4double omegaWidth = 0.0;
                                                   >> 100     G4ThreadLocal G4double etaPrimeWidth = 0.0;
                                                   >> 101     
                                                   >> 102     
172       const G4int mediumNucleiTableSize = 30;     103       const G4int mediumNucleiTableSize = 30;
173                                                   104 
174       const G4double mediumDiffuseness[mediumN    105       const G4double mediumDiffuseness[mediumNucleiTableSize] =
175       {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.69 << 106       {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.77,1.71,
176         1.69,1.72,1.635,1.730,1.81,1.833,1.798 << 107         1.69,1.69,1.635,1.730,1.81,1.833,1.798,
177         1.93,0.567,0.571, 0.560,0.549,0.550,0. << 108         1.841,0.567,0.571, 0.560,0.549,0.550,0.551,
178         0.580,0.575,0.569,0.537,0.0,0.0};         109         0.580,0.575,0.569,0.537,0.0,0.0};
179       const G4double mediumRadius[mediumNuclei    110       const G4double mediumRadius[mediumNucleiTableSize] =
180       {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0    111       {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0.631,0.838,
181         0.811,0.84,1.403,1.335,1.25,1.544,1.49 << 112         0.811,1.07,1.403,1.335,1.25,1.544,1.498,1.513,
182         2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3    113         2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3.03,2.84,
183         3.14,0.0,0.0};                            114         3.14,0.0,0.0};
184                                                   115 
185       const G4double positionRMS[clusterTableZ    116       const G4double positionRMS[clusterTableZSize][clusterTableASize] = {
186         /*     A=   0     1     2     3     4     117         /*     A=   0     1     2     3     4     5     6     7     8     9    10    11    12 */
187         /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1.    118         /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0},
188         /* Z=1 */ {-1.0, -1.0, 2.10, 1.80, 1.7    119         /* Z=1 */ {-1.0, -1.0, 2.10, 1.80, 1.70, 1.83, 2.60, 2.50, -1.0, -1.0, -1.0, -1.0, -1.0},
189         /* Z=2 */ {-1.0, -1.0, -1.0, 1.80, 1.6    120         /* Z=2 */ {-1.0, -1.0, -1.0, 1.80, 1.68, 1.70, 2.60, 2.50, 2.50, 2.50, 2.50, -1.0, -1.0},
190         /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 1.7    121         /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 1.70, 1.83, 2.56, 2.40, 2.50, 2.50, 2.50, 2.50, 2.50},
191         /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.    122         /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.60, 2.50, 2.50, 2.51, 2.50, 2.50, 2.50},
192         /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.    123         /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50, 2.50, 2.45, 2.40, 2.50},
193         /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.    124         /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50, 2.50, 2.47},
194         /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1.    125         /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50},
195         /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1.    126         /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50}
196       };                                          127       };
197                                                   128 
198       const G4double momentumRMS[clusterTableZ    129       const G4double momentumRMS[clusterTableZSize][clusterTableASize] = {
199         /*     A=   0     1     2     3     4     130         /*     A=   0     1     2     3     4     5     6     7     8     9    10    11    12 */
200         /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1.    131         /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0},
201         /* Z=1 */ {-1.0, -1.0, 77.0, 110., 153    132         /* Z=1 */ {-1.0, -1.0, 77.0, 110., 153., 100., 100., 100., -1.0, -1.0, -1.0, -1.0, -1.0},
202         /* Z=2 */ {-1.0, -1.0, -1.0, 110., 153    133         /* Z=2 */ {-1.0, -1.0, -1.0, 110., 153., 100., 100., 100., 100., 100., 100., -1.0, -1.0},
203         /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 153    134         /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 153., 100., 100., 100., 100., 100., 100., 100., 100.},
204         /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.    135         /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100., 100., 100.},
205         /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.    136         /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100., 100., 100.},
206         /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.    137         /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100.},
207         /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1.    138         /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100.},
208         /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1.    139         /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100.}
209       };                                          140       };
210                                                   141 
211       const G4int elementTableSize = 113; // u    142       const G4int elementTableSize = 113; // up to Cn
212                                                   143 
213       /// \brief Table of chemical element nam    144       /// \brief Table of chemical element names
214       const std::string elementTable[elementTa    145       const std::string elementTable[elementTableSize] = {
215         "",                                       146         "",
216         "H",                                      147         "H",
217         "He",                                     148         "He",
218         "Li",                                     149         "Li",
219         "Be",                                     150         "Be",
220         "B",                                      151         "B",
221         "C",                                      152         "C",
222         "N",                                      153         "N",
223         "O",                                      154         "O",
224         "F",                                      155         "F",
225         "Ne",                                     156         "Ne",
226         "Na",                                     157         "Na",
227         "Mg",                                     158         "Mg",
228         "Al",                                     159         "Al",
229         "Si",                                     160         "Si",
230         "P",                                      161         "P",
231         "S",                                      162         "S",
232         "Cl",                                     163         "Cl",
233         "Ar",                                     164         "Ar",
234         "K",                                      165         "K",
235         "Ca",                                     166         "Ca",
236         "Sc",                                     167         "Sc",
237         "Ti",                                     168         "Ti",
238         "V",                                      169         "V",
239         "Cr",                                     170         "Cr",
240         "Mn",                                     171         "Mn",
241         "Fe",                                     172         "Fe",
242         "Co",                                     173         "Co",
243         "Ni",                                     174         "Ni",
244         "Cu",                                     175         "Cu",
245         "Zn",                                     176         "Zn",
246         "Ga",                                     177         "Ga",
247         "Ge",                                     178         "Ge",
248         "As",                                     179         "As",
249         "Se",                                     180         "Se",
250         "Br",                                     181         "Br",
251         "Kr",                                     182         "Kr",
252         "Rb",                                     183         "Rb",
253         "Sr",                                     184         "Sr",
254         "Y",                                      185         "Y",
255         "Zr",                                     186         "Zr",
256         "Nb",                                     187         "Nb",
257         "Mo",                                     188         "Mo",
258         "Tc",                                     189         "Tc",
259         "Ru",                                     190         "Ru",
260         "Rh",                                     191         "Rh",
261         "Pd",                                     192         "Pd",
262         "Ag",                                     193         "Ag",
263         "Cd",                                     194         "Cd",
264         "In",                                     195         "In",
265         "Sn",                                     196         "Sn",
266         "Sb",                                     197         "Sb",
267         "Te",                                     198         "Te",
268         "I",                                      199         "I",
269         "Xe",                                     200         "Xe",
270         "Cs",                                     201         "Cs",
271         "Ba",                                     202         "Ba",
272         "La",                                     203         "La",
273         "Ce",                                     204         "Ce",
274         "Pr",                                     205         "Pr",
275         "Nd",                                     206         "Nd",
276         "Pm",                                     207         "Pm",
277         "Sm",                                     208         "Sm",
278         "Eu",                                     209         "Eu",
279         "Gd",                                     210         "Gd",
280         "Tb",                                     211         "Tb",
281         "Dy",                                     212         "Dy",
282         "Ho",                                     213         "Ho",
283         "Er",                                     214         "Er",
284         "Tm",                                     215         "Tm",
285         "Yb",                                     216         "Yb",
286         "Lu",                                     217         "Lu",
287         "Hf",                                     218         "Hf",
288         "Ta",                                     219         "Ta",
289         "W",                                      220         "W",
290         "Re",                                     221         "Re",
291         "Os",                                     222         "Os",
292         "Ir",                                     223         "Ir",
293         "Pt",                                     224         "Pt",
294         "Au",                                     225         "Au",
295         "Hg",                                     226         "Hg",
296         "Tl",                                     227         "Tl",
297         "Pb",                                     228         "Pb",
298         "Bi",                                     229         "Bi",
299         "Po",                                     230         "Po",
300         "At",                                     231         "At",
301         "Rn",                                     232         "Rn",
302         "Fr",                                     233         "Fr",
303         "Ra",                                     234         "Ra",
304         "Ac",                                     235         "Ac",
305         "Th",                                     236         "Th",
306         "Pa",                                     237         "Pa",
307         "U",                                      238         "U",
308         "Np",                                     239         "Np",
309         "Pu",                                     240         "Pu",
310         "Am",                                     241         "Am",
311         "Cm",                                     242         "Cm",
312         "Bk",                                     243         "Bk",
313         "Cf",                                     244         "Cf",
314         "Es",                                     245         "Es",
315         "Fm",                                     246         "Fm",
316         "Md",                                     247         "Md",
317         "No",                                     248         "No",
318         "Lr",                                     249         "Lr",
319         "Rf",                                     250         "Rf",
320         "Db",                                     251         "Db",
321         "Sg",                                     252         "Sg",
322         "Bh",                                     253         "Bh",
323         "Hs",                                     254         "Hs",
324         "Mt",                                     255         "Mt",
325         "Ds",                                     256         "Ds",
326         "Rg",                                     257         "Rg",
327         "Cn"                                      258         "Cn"
328       };                                          259       };
329                                                   260 
330       /// \brief Digit names to compose IUPAC     261       /// \brief Digit names to compose IUPAC element names
331       const std::string elementIUPACDigits = "    262       const std::string elementIUPACDigits = "nubtqphsoe";
332                                                   263 
333 #define INCL_DEFAULT_SEPARATION_ENERGY 6.83       264 #define INCL_DEFAULT_SEPARATION_ENERGY 6.83
334       const G4double theINCLProtonSeparationEn    265       const G4double theINCLProtonSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
335       const G4double theINCLNeutronSeparationE    266       const G4double theINCLNeutronSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
336       const G4double theINCLLambdaSeparationEn << 
337       //const G4double theINCLantiProtonSepara << 
338       const G4double theINCLantiProtonSeparati << 
339       G4ThreadLocal G4double protonSeparationE    267       G4ThreadLocal G4double protonSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
340       G4ThreadLocal G4double neutronSeparation    268       G4ThreadLocal G4double neutronSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
341       G4ThreadLocal G4double lambdaSeparationE << 
342       //G4ThreadLocal G4double antiprotonSepar << 
343       //G4ThreadLocal G4double antiprotonSepar << 
344 #undef INCL_DEFAULT_SEPARATION_ENERGY             269 #undef INCL_DEFAULT_SEPARATION_ENERGY
345                                                   270 
346       G4ThreadLocal G4double rpCorrelationCoef    271       G4ThreadLocal G4double rpCorrelationCoefficient[UnknownParticle];
347                                                   272 
348       G4ThreadLocal G4double neutronSkin = 0.0    273       G4ThreadLocal G4double neutronSkin = 0.0;
349       G4ThreadLocal G4double neutronHalo = 0.0    274       G4ThreadLocal G4double neutronHalo = 0.0;
350                                                   275 
351 #ifdef INCLXX_IN_GEANT4_MODE                      276 #ifdef INCLXX_IN_GEANT4_MODE
352       G4ThreadLocal G4IonTable *theG4IonTable;    277       G4ThreadLocal G4IonTable *theG4IonTable;
353 #endif                                            278 #endif
354                                                   279 
355       /// \brief Default value for constant Fe    280       /// \brief Default value for constant Fermi momentum
356       G4ThreadLocal G4double constantFermiMome    281       G4ThreadLocal G4double constantFermiMomentum = 0.0;
357                                                   282 
358       /// \brief Transform a IUPAC char to an     283       /// \brief Transform a IUPAC char to an char representing an integer digit
359       char iupacToInt(char c) {                   284       char iupacToInt(char c) {
360         return (char)(((G4int)'0')+elementIUPA    285         return (char)(((G4int)'0')+elementIUPACDigits.find(c));
361       }                                           286       }
362                                                   287 
363       /// \brief Transform an integer digit (r    288       /// \brief Transform an integer digit (represented by a char) to a IUPAC char
364       char intToIUPAC(char n) { return element    289       char intToIUPAC(char n) { return elementIUPACDigits.at(n); }
365                                                   290 
366       /// \brief Get the singleton instance of    291       /// \brief Get the singleton instance of the natural isotopic distributions
367       const NaturalIsotopicDistributions *getN    292       const NaturalIsotopicDistributions *getNaturalIsotopicDistributions() {
368         if(!theNaturalIsotopicDistributions)      293         if(!theNaturalIsotopicDistributions)
369           theNaturalIsotopicDistributions = ne    294           theNaturalIsotopicDistributions = new NaturalIsotopicDistributions;
370         return theNaturalIsotopicDistributions    295         return theNaturalIsotopicDistributions;
371       }                                           296       }
372                                                   297 
373     } // namespace                                298     } // namespace
374                                                   299 
375     void initialize(Config const * const theCo    300     void initialize(Config const * const theConfig /*=0*/) {
376       protonMass = theINCLNucleonMass;            301       protonMass = theINCLNucleonMass;
377       neutronMass = theINCLNucleonMass;        << 302       neutronMass = theINCLNucleonMass;
378       piPlusMass = theINCLPionMass;               303       piPlusMass = theINCLPionMass;
379       piMinusMass = theINCLPionMass;              304       piMinusMass = theINCLPionMass;
380       piZeroMass = theINCLPionMass;               305       piZeroMass = theINCLPionMass;
381                                                << 
382       etaMass = theINCLEtaMass;                   306       etaMass = theINCLEtaMass;
383       omegaMass = theINCLOmegaMass;               307       omegaMass = theINCLOmegaMass;
384       etaPrimeMass = theINCLEtaPrimeMass;         308       etaPrimeMass = theINCLEtaPrimeMass;
385       photonMass = theINCLPhotonMass;             309       photonMass = theINCLPhotonMass;
386                                                   310 
387       SigmaPlusMass = theRealSigmaPlusMass;    << 
388       SigmaMinusMass = theRealSigmaMinusMass;  << 
389       SigmaZeroMass = theRealSigmaZeroMass;    << 
390       LambdaMass = theINCLLambdaMass;          << 
391       KPlusMass = theRealChargedKaonMass;      << 
392       KZeroMass = theRealNeutralKaonMass;      << 
393       KZeroBarMass = theRealNeutralKaonMass;   << 
394       KShortMass = theRealNeutralKaonMass;     << 
395       KLongMass = theRealNeutralKaonMass;      << 
396       KMinusMass = theRealChargedKaonMass;     << 
397                                                << 
398       antiProtonMass = theRealAntiProtonMass;  << 
399       XiZeroMass = theRealXiZeroMass;          << 
400       XiMinusMass = theRealXiMinusMass;        << 
401       antiNeutronMass = theRealAntiNeutronMass << 
402       antiSigmaPlusMass = theRealAntiSigmaPlus << 
403       antiSigmaMinusMass = theRealAntiSigmaMin << 
404       antiSigmaZeroMass = theRealAntiSigmaZero << 
405       antiLambdaMass = theRealAntiLambdaMass;  << 
406       antiXiZeroMass = theRealAntiXiZeroMass;  << 
407       antiXiMinusMass = theRealAntiXiMinusMass << 
408                                                << 
409       if(theConfig && theConfig->getUseRealMas    311       if(theConfig && theConfig->getUseRealMasses()) {
410         getTableMass = getRealMass;               312         getTableMass = getRealMass;
411         getTableParticleMass = getRealMass;       313         getTableParticleMass = getRealMass;
412       } else {                                    314       } else {
413         getTableMass = getINCLMass;               315         getTableMass = getINCLMass;
414         getTableParticleMass = getINCLMass;       316         getTableParticleMass = getINCLMass;
415       }                                           317       }
416                                                   318 
417 #ifndef INCLXX_IN_GEANT4_MODE                     319 #ifndef INCLXX_IN_GEANT4_MODE
418       std::string dataFilePath;                   320       std::string dataFilePath;
419       if(theConfig)                               321       if(theConfig)
420         dataFilePath = theConfig->getINCLXXDat    322         dataFilePath = theConfig->getINCLXXDataFilePath();
421       NuclearMassTable::initialize(dataFilePat    323       NuclearMassTable::initialize(dataFilePath, getRealMass(Proton), getRealMass(Neutron));
422 #endif                                            324 #endif
423                                                   325 
424 #ifdef INCLXX_IN_GEANT4_MODE                      326 #ifdef INCLXX_IN_GEANT4_MODE
425       G4ParticleTable *theG4ParticleTable = G4    327       G4ParticleTable *theG4ParticleTable = G4ParticleTable::GetParticleTable();
426       theG4IonTable = theG4ParticleTable->GetI    328       theG4IonTable = theG4ParticleTable->GetIonTable();
427       theRealProtonMass = theG4ParticleTable->    329       theRealProtonMass = theG4ParticleTable->FindParticle("proton")->GetPDGMass() / MeV;
428       theRealNeutronMass = theG4ParticleTable- << 330       theRealNeutronMass = theG4ParticleTable->FindParticle("neutron")->GetPDGMass() / MeV;
429       theRealChargedPiMass = theG4ParticleTabl    331       theRealChargedPiMass = theG4ParticleTable->FindParticle("pi+")->GetPDGMass() / MeV;
430       theRealPiZeroMass = theG4ParticleTable->    332       theRealPiZeroMass = theG4ParticleTable->FindParticle("pi0")->GetPDGMass() / MeV;
431                                                << 
432       theRealEtaMass = theG4ParticleTable->Fin    333       theRealEtaMass = theG4ParticleTable->FindParticle("eta")->GetPDGMass() / MeV;
433       theRealOmegaMass = theG4ParticleTable->F    334       theRealOmegaMass = theG4ParticleTable->FindParticle("omega")->GetPDGMass() / MeV;
434       theRealEtaPrimeMass = theG4ParticleTable    335       theRealEtaPrimeMass = theG4ParticleTable->FindParticle("eta_prime")->GetPDGMass() / MeV;
435       theRealPhotonMass = theG4ParticleTable->    336       theRealPhotonMass = theG4ParticleTable->FindParticle("gamma")->GetPDGMass() / MeV;
436                                                << 
437       theRealSigmaPlusMass = theG4ParticleTabl << 
438       theRealSigmaZeroMass = theG4ParticleTabl << 
439       theRealSigmaMinusMass = theG4ParticleTab << 
440       theRealLambdaMass = theG4ParticleTable-> << 
441       theRealChargedKaonMass = theG4ParticleTa << 
442       theRealNeutralKaonMass = theG4ParticleTa << 
443                                                << 
444       theRealAntiProtonMass = theG4ParticleTab << 
445       theRealAntiNeutronMass = theG4ParticleTa << 
446       theRealXiZeroMass = theG4ParticleTable-> << 
447       theRealXiMinusMass = theG4ParticleTable- << 
448       theRealAntiSigmaPlusMass = theG4Particle << 
449       theRealAntiSigmaZeroMass = theG4Particle << 
450       theRealAntiSigmaMinusMass = theG4Particl << 
451       theRealAntiLambdaMass = theG4ParticleTab << 
452       theRealAntiXiZeroMass = theG4ParticleTab << 
453       theRealAntiXiMinusMass = theG4ParticleTa << 
454 #endif                                            337 #endif
455                                                   338 
456       minDeltaMass = theRealNeutronMass + theR    339       minDeltaMass = theRealNeutronMass + theRealChargedPiMass + 0.5;
457       minDeltaMass2 = minDeltaMass*minDeltaMas    340       minDeltaMass2 = minDeltaMass*minDeltaMass;
458       minDeltaMassRndm = std::atan((minDeltaMa    341       minDeltaMassRndm = std::atan((minDeltaMass-effectiveDeltaMass)*2./effectiveDeltaWidth);
459                                                   342 
460       piPlusWidth   = theChargedPiWidth;       << 343     piPlusWidth   = theChargedPiWidth;
461       piMinusWidth  = theChargedPiWidth;       << 344     piMinusWidth  = theChargedPiWidth;
462       piZeroWidth   = thePiZeroWidth;          << 345     piZeroWidth   = thePiZeroWidth;
463       etaWidth      = theEtaWidth;             << 346     etaWidth      = theEtaWidth;
464       omegaWidth    = theOmegaWidth;           << 347     omegaWidth    = theOmegaWidth;
465       etaPrimeWidth = theEtaPrimeWidth;        << 348     etaPrimeWidth = theEtaPrimeWidth;
466                                                << 349     
467       SigmaMinusWidth = theSigmaMinusWidth;    << 350     
468       SigmaPlusWidth = theSigmaPlusWidth;      << 
469       SigmaZeroWidth = theSigmaZeroWidth;      << 
470       LambdaWidth = theLambdaWidth;            << 
471       KPlusWidth = theChargedKaonWidth;        << 
472       KMinusWidth = theChargedKaonWidth;       << 
473       KShortWidth = theKShortWidth;            << 
474       KLongWidth = theKLongWidth;              << 
475                                                << 
476       // Initialise HFB tables                 << 
477 #ifdef INCLXX_IN_GEANT4_MODE                   << 
478         HFB::initialize();                     << 
479 #else                                          << 
480         HFB::initialize(dataFilePath);         << 
481 #endif                                         << 
482                                                << 
483       // Initialise the separation-energy func    351       // Initialise the separation-energy function
484       if(!theConfig || theConfig->getSeparatio    352       if(!theConfig || theConfig->getSeparationEnergyType()==INCLSeparationEnergy)
485         getSeparationEnergy = getSeparationEne    353         getSeparationEnergy = getSeparationEnergyINCL;
486       else if(theConfig->getSeparationEnergyTy    354       else if(theConfig->getSeparationEnergyType()==RealSeparationEnergy)
487         getSeparationEnergy = getSeparationEne    355         getSeparationEnergy = getSeparationEnergyReal;
488       else if(theConfig->getSeparationEnergyTy    356       else if(theConfig->getSeparationEnergyType()==RealForLightSeparationEnergy)
489         getSeparationEnergy = getSeparationEne    357         getSeparationEnergy = getSeparationEnergyRealForLight;
490       else {                                      358       else {
491         INCL_FATAL("Unrecognized separation-en    359         INCL_FATAL("Unrecognized separation-energy type in ParticleTable initialization: " << theConfig->getSeparationEnergyType() << '\n');
492         return;                                   360         return;
493       }                                           361       }
494                                                   362 
495       // Initialise the Fermi-momentum functio    363       // Initialise the Fermi-momentum function
496       if(!theConfig || theConfig->getFermiMome    364       if(!theConfig || theConfig->getFermiMomentumType()==ConstantFermiMomentum) {
497         getFermiMomentum = ParticleTable::getF    365         getFermiMomentum = ParticleTable::getFermiMomentumConstant;
498         if(theConfig) {                           366         if(theConfig) {
499           const G4double aFermiMomentum = theC    367           const G4double aFermiMomentum = theConfig->getFermiMomentum();
500           if(aFermiMomentum>0.)                   368           if(aFermiMomentum>0.)
501             constantFermiMomentum = aFermiMome    369             constantFermiMomentum = aFermiMomentum;
502           else                                    370           else
503             constantFermiMomentum = PhysicalCo    371             constantFermiMomentum = PhysicalConstants::Pf;
504         } else {                                  372         } else {
505           constantFermiMomentum = PhysicalCons    373           constantFermiMomentum = PhysicalConstants::Pf;
506         }                                         374         }
507       } else if(theConfig->getFermiMomentumTyp    375       } else if(theConfig->getFermiMomentumType()==ConstantLightFermiMomentum)
508         getFermiMomentum = ParticleTable::getF    376         getFermiMomentum = ParticleTable::getFermiMomentumConstantLight;
509       else if(theConfig->getFermiMomentumType(    377       else if(theConfig->getFermiMomentumType()==MassDependentFermiMomentum)
510         getFermiMomentum = ParticleTable::getF    378         getFermiMomentum = ParticleTable::getFermiMomentumMassDependent;
511       else {                                      379       else {
512         INCL_FATAL("Unrecognized Fermi-momentu    380         INCL_FATAL("Unrecognized Fermi-momentum type in ParticleTable initialization: " << theConfig->getFermiMomentumType() << '\n');
513         return;                                   381         return;
514       }                                           382       }
515                                                   383 
516       // Initialise the r-p correlation coeffi    384       // Initialise the r-p correlation coefficients
517       std::fill(rpCorrelationCoefficient, rpCo    385       std::fill(rpCorrelationCoefficient, rpCorrelationCoefficient + UnknownParticle, 1.);
518       if(theConfig) {                             386       if(theConfig) {
519         rpCorrelationCoefficient[Proton] = the    387         rpCorrelationCoefficient[Proton] = theConfig->getRPCorrelationCoefficient(Proton);
520         rpCorrelationCoefficient[Neutron] = th    388         rpCorrelationCoefficient[Neutron] = theConfig->getRPCorrelationCoefficient(Neutron);
521       }                                           389       }
522                                                   390 
523       // Initialise the neutron-skin parameter    391       // Initialise the neutron-skin parameters
524       if(theConfig) {                             392       if(theConfig) {
525         neutronSkin = theConfig->getNeutronSki    393         neutronSkin = theConfig->getNeutronSkin();
526         neutronHalo = theConfig->getNeutronHal    394         neutronHalo = theConfig->getNeutronHalo();
527       }                                           395       }
528                                                   396 
529     }                                             397     }
530                                                   398 
531     G4int getIsospin(const ParticleType t) {      399     G4int getIsospin(const ParticleType t) {
532       // Actually this is the 3rd component of    400       // Actually this is the 3rd component of isospin (I_z) multiplied by 2!
533       if(t == Proton) {                           401       if(t == Proton) {
534         return 1;                                 402         return 1;
535       } else if(t == Neutron) {                   403       } else if(t == Neutron) {
536         return -1;                                404         return -1;
537       } else if(t == PiPlus) {                    405       } else if(t == PiPlus) {
538         return 2;                                 406         return 2;
539       } else if(t == PiMinus) {                   407       } else if(t == PiMinus) {
540         return -2;                                408         return -2;
541       } else if(t == PiZero) {                    409       } else if(t == PiZero) {
542         return 0;                                 410         return 0;
543       } else if(t == DeltaPlusPlus) {             411       } else if(t == DeltaPlusPlus) {
544         return 3;                                 412         return 3;
545       } else if(t == DeltaPlus) {                 413       } else if(t == DeltaPlus) {
546         return 1;                                 414         return 1;
547       } else if(t == DeltaZero) {                 415       } else if(t == DeltaZero) {
548         return -1;                                416         return -1;
549       } else if(t == DeltaMinus) {                417       } else if(t == DeltaMinus) {
550         return -3;                             << 418     return -3;
551       } else if(t == Lambda) {                 << 
552         return 0;                              << 
553       } else if(t == SigmaPlus) {              << 
554         return 2;                              << 
555       } else if(t == SigmaZero) {              << 
556         return 0;                              << 
557       } else if(t == SigmaMinus) {             << 
558         return -2;                             << 
559       } else if(t == KPlus) {                  << 
560         return 1;                              << 
561       } else if(t == KZero) {                  << 
562         return -1;                             << 
563       } else if(t == KZeroBar) {               << 
564         return 1;                              << 
565       } else if(t == KShort) {                 << 
566         return 0;                              << 
567       } else if(t == KLong) {                  << 
568         return 0;                              << 
569       } else if(t == KMinus) {                 << 
570         return -1;                             << 
571       } else if(t == Eta) {                       419       } else if(t == Eta) {
572         return 0;                              << 420     return 0;
573       } else if(t == Omega) {                     421       } else if(t == Omega) {
574         return 0;                              << 422     return 0;
575       } else if(t == EtaPrime) {                  423       } else if(t == EtaPrime) {
576         return 0;                              << 424     return 0;
577       } else if(t == Photon) {                    425       } else if(t == Photon) {
578         return 0;                              << 426     return 0;
579       } else if(t == antiProton) {             << 
580         return -1;                             << 
581       } else if(t == XiMinus) {                << 
582         return -1;                             << 
583       } else if(t == XiZero) {                 << 
584         return 1;                              << 
585       } else if(t == antiNeutron) {            << 
586         return 1;                              << 
587       } else if(t == antiLambda) {             << 
588         return 0;                              << 
589       } else if(t == antiSigmaPlus) {          << 
590         return -2;                             << 
591       } else if(t == antiSigmaZero) {          << 
592         return 0;                              << 
593       } else if(t == antiSigmaMinus) {         << 
594         return 2;                              << 
595       } else if(t == antiXiMinus) {            << 
596         return 1;                              << 
597       } else if(t == antiXiZero) {             << 
598         return -1;                             << 
599       }                                           427       }
                                                   >> 428 
600       INCL_ERROR("Requested isospin of an unkn    429       INCL_ERROR("Requested isospin of an unknown particle!");
601       return -10; // Unknown                      430       return -10; // Unknown
602     }                                             431     }
603                                                   432 
604     std::string getShortName(const ParticleSpe << 433     std::string getShortName(const ParticleSpecies &s) {
605       if(sp.theType==Composite && sp.theS == 0 << 434       if(s.theType==Composite)
606         return getShortName(sp.theA,sp.theZ);  << 435         return getShortName(s.theA,s.theZ);
607       else if(sp.theType==Composite)           << 
608         return getName(sp.theA,sp.theZ,sp.theS << 
609       else                                        436       else
610         return getShortName(sp.theType);       << 437         return getShortName(s.theType);
611     }                                             438     }
612                                                << 439 
613     std::string getName(const ParticleSpecies  << 440     std::string getName(const ParticleSpecies &s) {
614       if(sp.theType==Composite && sp.theS == 0 << 441       if(s.theType==Composite)
615         return getName(sp.theA,sp.theZ);       << 442         return getName(s.theA,s.theZ);
616       else if(sp.theType==Composite)           << 
617         return getName(sp.theA,sp.theZ,sp.theS << 
618       else                                        443       else
619         return getName(sp.theType);            << 444         return getName(s.theType);
620     }                                             445     }
621                                                   446 
622     std::string getName(const G4int A, const G    447     std::string getName(const G4int A, const G4int Z) {
623       std::stringstream stream;                   448       std::stringstream stream;
624       stream << getElementName(Z) << "-" << A;    449       stream << getElementName(Z) << "-" << A;
625       return stream.str();                        450       return stream.str();
626     }                                             451     }
627                                                   452 
628     std::string getName(const G4int A, const G << 
629       std::stringstream stream;                << 
630       if(S >= 0) // S < 0 for hypernuclei      << 
631         return getName(A, Z);                  << 
632       else if(S == -1)                         << 
633         stream << getElementName(Z) << "-" <<  << 
634       else                                     << 
635         stream << getElementName(Z) << "-" <<  << 
636       return stream.str();                     << 
637     }                                          << 
638                                                << 
639     std::string getShortName(const G4int A, co    453     std::string getShortName(const G4int A, const G4int Z) {
640       std::stringstream stream;                   454       std::stringstream stream;
641       stream << getElementName(Z);                455       stream << getElementName(Z);
642       if(A>0)                                     456       if(A>0)
643         stream << A;                              457         stream << A;
644       return stream.str();                        458       return stream.str();
645     }                                             459     }
646                                                   460 
647     std::string getName(const ParticleType p)     461     std::string getName(const ParticleType p) {
648       if(p == G4INCL::Proton) {                   462       if(p == G4INCL::Proton) {
649         return std::string("proton");             463         return std::string("proton");
650       } else if(p == G4INCL::Neutron) {           464       } else if(p == G4INCL::Neutron) {
651         return std::string("neutron");            465         return std::string("neutron");
652       } else if(p == G4INCL::DeltaPlusPlus) {     466       } else if(p == G4INCL::DeltaPlusPlus) {
653         return std::string("delta++");            467         return std::string("delta++");
654       } else if(p == G4INCL::DeltaPlus) {         468       } else if(p == G4INCL::DeltaPlus) {
655         return std::string("delta+");             469         return std::string("delta+");
656       } else if(p == G4INCL::DeltaZero) {         470       } else if(p == G4INCL::DeltaZero) {
657         return std::string("delta0");             471         return std::string("delta0");
658       } else if(p == G4INCL::DeltaMinus) {        472       } else if(p == G4INCL::DeltaMinus) {
659         return std::string("delta-");             473         return std::string("delta-");
660       } else if(p == G4INCL::PiPlus) {            474       } else if(p == G4INCL::PiPlus) {
661         return std::string("pi+");                475         return std::string("pi+");
662       } else if(p == G4INCL::PiZero) {            476       } else if(p == G4INCL::PiZero) {
663         return std::string("pi0");                477         return std::string("pi0");
664       } else if(p == G4INCL::PiMinus) {           478       } else if(p == G4INCL::PiMinus) {
665         return std::string("pi-");                479         return std::string("pi-");
666       } else if(p == G4INCL::Lambda) {         << 
667         return std::string("lambda");          << 
668       } else if(p == G4INCL::SigmaPlus) {      << 
669         return std::string("sigma+");          << 
670       } else if(p == G4INCL::SigmaZero) {      << 
671         return std::string("sigma0");          << 
672       } else if(p == G4INCL::SigmaMinus) {     << 
673         return std::string("sigma-");          << 
674       } else if(p == G4INCL::antiProton) {     << 
675         return std::string("antiproton");      << 
676       } else if(p == G4INCL::XiMinus) {        << 
677         return std::string("xi-");             << 
678       } else if(p == G4INCL::XiZero) {         << 
679         return std::string("xi0");             << 
680       } else if(p == G4INCL::antiNeutron) {    << 
681         return std::string("antineutron");     << 
682       } else if(p == G4INCL::antiSigmaPlus) {  << 
683         return std::string("antisigma+");      << 
684       } else if(p == G4INCL::antiSigmaZero) {  << 
685         return std::string("antisigma0");      << 
686       } else if(p == G4INCL::antiSigmaMinus) { << 
687         return std::string("antisigma-");      << 
688       } else if(p == G4INCL::antiLambda) {     << 
689         return std::string("antilambda");      << 
690       } else if(p == G4INCL::antiXiMinus) {    << 
691         return std::string("antixi-");         << 
692       } else if(p == G4INCL::antiXiZero) {     << 
693         return std::string("antixi0");         << 
694       } else if(p == G4INCL::KPlus) {          << 
695         return std::string("kaon+");           << 
696       } else if(p == G4INCL::KZero) {          << 
697         return std::string("kaon0");           << 
698       } else if(p == G4INCL::KZeroBar) {       << 
699         return std::string("kaon0bar");        << 
700       } else if(p == G4INCL::KMinus) {         << 
701         return std::string("kaon-");           << 
702       } else if(p == G4INCL::KShort) {         << 
703         return std::string("kaonshort");       << 
704       } else if(p == G4INCL::KLong) {          << 
705         return std::string("kaonlong");        << 
706       } else if(p == G4INCL::Composite) {         480       } else if(p == G4INCL::Composite) {
707         return std::string("composite");          481         return std::string("composite");
708       } else if(p == G4INCL::Eta) {               482       } else if(p == G4INCL::Eta) {
709         return std::string("eta");                483         return std::string("eta");
710       } else if(p == G4INCL::Omega) {             484       } else if(p == G4INCL::Omega) {
711         return std::string("omega");              485         return std::string("omega");
712       } else if(p == G4INCL::EtaPrime) {          486       } else if(p == G4INCL::EtaPrime) {
713         return std::string("etaprime");           487         return std::string("etaprime");
714       } else if(p == G4INCL::Photon) {            488       } else if(p == G4INCL::Photon) {
715         return std::string("photon");          << 489         return std::string("photon");
716       }                                           490       }
717       return std::string("unknown");              491       return std::string("unknown");
718     }                                             492     }
719                                                   493 
720     std::string getShortName(const ParticleTyp    494     std::string getShortName(const ParticleType p) {
721       if(p == G4INCL::Proton) {                   495       if(p == G4INCL::Proton) {
722         return std::string("p");                  496         return std::string("p");
723       } else if(p == G4INCL::Neutron) {           497       } else if(p == G4INCL::Neutron) {
724         return std::string("n");                  498         return std::string("n");
725       } else if(p == G4INCL::DeltaPlusPlus) {     499       } else if(p == G4INCL::DeltaPlusPlus) {
726         return std::string("d++");                500         return std::string("d++");
727       } else if(p == G4INCL::DeltaPlus) {         501       } else if(p == G4INCL::DeltaPlus) {
728         return std::string("d+");                 502         return std::string("d+");
729       } else if(p == G4INCL::DeltaZero) {         503       } else if(p == G4INCL::DeltaZero) {
730         return std::string("d0");                 504         return std::string("d0");
731       } else if(p == G4INCL::DeltaMinus) {        505       } else if(p == G4INCL::DeltaMinus) {
732         return std::string("d-");                 506         return std::string("d-");
733       } else if(p == G4INCL::PiPlus) {            507       } else if(p == G4INCL::PiPlus) {
734         return std::string("pi+");                508         return std::string("pi+");
735       } else if(p == G4INCL::PiZero) {            509       } else if(p == G4INCL::PiZero) {
736         return std::string("pi0");                510         return std::string("pi0");
737       } else if(p == G4INCL::PiMinus) {           511       } else if(p == G4INCL::PiMinus) {
738         return std::string("pi-");                512         return std::string("pi-");
739       } else if(p == G4INCL::Lambda) {         << 
740         return std::string("l");               << 
741       } else if(p == G4INCL::SigmaPlus) {      << 
742         return std::string("s+");              << 
743       } else if(p == G4INCL::SigmaZero) {      << 
744         return std::string("s0");              << 
745       } else if(p == G4INCL::SigmaMinus) {     << 
746         return std::string("s-");              << 
747       } else if(p == G4INCL::antiProton) {     << 
748         return std::string("pb");              << 
749       } else if(p == G4INCL::XiMinus) {        << 
750         return std::string("x-");              << 
751       } else if(p == G4INCL::XiZero) {         << 
752         return std::string("x0");              << 
753       } else if(p == G4INCL::antiNeutron) {    << 
754         return std::string("nb");              << 
755       } else if(p == G4INCL::antiSigmaPlus) {  << 
756         return std::string("s+b");             << 
757       } else if(p == G4INCL::antiSigmaZero) {  << 
758         return std::string("s0b");             << 
759       } else if(p == G4INCL::antiSigmaMinus) { << 
760         return std::string("s-b");             << 
761       } else if(p == G4INCL::antiLambda) {     << 
762         return std::string("lb");              << 
763       } else if(p == G4INCL::antiXiMinus) {    << 
764         return std::string("x-b");             << 
765       } else if(p == G4INCL::antiXiZero) {     << 
766         return std::string("x0b");             << 
767       } else if(p == G4INCL::KPlus) {          << 
768         return std::string("k+");              << 
769       } else if(p == G4INCL::KZero) {          << 
770         return std::string("k0");              << 
771       } else if(p == G4INCL::KZeroBar) {       << 
772         return std::string("k0b");             << 
773       } else if(p == G4INCL::KMinus) {         << 
774         return std::string("k-");              << 
775       } else if(p == G4INCL::KShort) {         << 
776         return std::string("ks");              << 
777       } else if(p == G4INCL::KLong) {          << 
778         return std::string("kl");              << 
779       } else if(p == G4INCL::Composite) {         513       } else if(p == G4INCL::Composite) {
780         return std::string("comp");               514         return std::string("comp");
781       } else if(p == G4INCL::Eta) {               515       } else if(p == G4INCL::Eta) {
782         return std::string("eta");                516         return std::string("eta");
783       } else if(p == G4INCL::Omega) {             517       } else if(p == G4INCL::Omega) {
784         return std::string("omega");              518         return std::string("omega");
785       } else if(p == G4INCL::EtaPrime) {          519       } else if(p == G4INCL::EtaPrime) {
786         return std::string("etap");               520         return std::string("etap");
787       } else if(p == G4INCL::Photon) {            521       } else if(p == G4INCL::Photon) {
788         return std::string("photon");             522         return std::string("photon");
789       }                                           523       }
790       return std::string("unknown");              524       return std::string("unknown");
791     }                                             525     }
792                                                   526 
793     G4double getINCLMass(const ParticleType pt    527     G4double getINCLMass(const ParticleType pt) {
794       if(pt == Proton) {                          528       if(pt == Proton) {
795         return protonMass;                        529         return protonMass;
796       } else if(pt == Neutron) {                  530       } else if(pt == Neutron) {
797         return neutronMass;                       531         return neutronMass;
798       } else if(pt == PiPlus) {                   532       } else if(pt == PiPlus) {
799         return piPlusMass;                        533         return piPlusMass;
800       } else if(pt == PiMinus) {                  534       } else if(pt == PiMinus) {
801         return piMinusMass;                       535         return piMinusMass;
802       } else if(pt == PiZero) {                   536       } else if(pt == PiZero) {
803         return piZeroMass;                        537         return piZeroMass;
804       } else if(pt == SigmaPlus) {             << 
805         return SigmaPlusMass;                  << 
806       } else if(pt == SigmaMinus) {            << 
807         return SigmaMinusMass;                 << 
808       } else if(pt == SigmaZero) {             << 
809         return SigmaZeroMass;                  << 
810       } else if(pt == Lambda) {                << 
811         return LambdaMass;                     << 
812       } else if(pt == antiProton) {            << 
813         return antiProtonMass;                 << 
814       } else if(pt == XiMinus) {               << 
815         return XiMinusMass;                    << 
816       } else if(pt == XiZero) {                << 
817         return XiZeroMass;                     << 
818       } else if(pt == antiNeutron) {           << 
819         return antiNeutronMass;                << 
820       } else if(pt == antiSigmaPlus) {         << 
821         return antiSigmaPlusMass;              << 
822       } else if(pt == antiSigmaMinus) {        << 
823         return antiSigmaMinusMass;             << 
824       } else if(pt == antiSigmaZero) {         << 
825         return antiSigmaZeroMass;              << 
826       } else if(pt == antiLambda) {            << 
827         return antiLambdaMass;                 << 
828       } else if(pt == antiXiMinus) {           << 
829         return antiXiMinusMass;                << 
830       } else if(pt == antiXiZero) {            << 
831         return antiXiZeroMass;                 << 
832       } else if(pt == KPlus) {                 << 
833         return KPlusMass;                      << 
834       } else if(pt == KZero) {                 << 
835         return KZeroMass;                      << 
836       } else if(pt == KZeroBar) {              << 
837         return KZeroBarMass;                   << 
838       } else if(pt == KMinus) {                << 
839         return KMinusMass;                     << 
840       } else if(pt == KShort) {                << 
841         return KShortMass;                     << 
842       } else if(pt == KLong) {                 << 
843         return KLongMass;                      << 
844       } else if(pt == Eta) {                      538       } else if(pt == Eta) {
845         return etaMass;                           539         return etaMass;
846       } else if(pt == Omega) {                    540       } else if(pt == Omega) {
847         return omegaMass;                         541         return omegaMass;
848       } else if(pt == EtaPrime) {                 542       } else if(pt == EtaPrime) {
849         return etaPrimeMass;                      543         return etaPrimeMass;
850       } else if(pt == Photon) {                   544       } else if(pt == Photon) {
851         return photonMass;                     << 545         return photonMass;
852       } else {                                    546       } else {
853         INCL_ERROR("getMass : Unknown particle    547         INCL_ERROR("getMass : Unknown particle type." << '\n');
854         return 0.0;                               548         return 0.0;
855       }                                           549       }
856     }                                             550     }
857                                                << 551 
858     G4double getRealMass(const ParticleType t)    552     G4double getRealMass(const ParticleType t) {
859       switch(t) {                                 553       switch(t) {
860         case Proton:                              554         case Proton:
861           return theRealProtonMass;               555           return theRealProtonMass;
862           break;                                  556           break;
863         case Neutron:                             557         case Neutron:
864           return theRealNeutronMass;              558           return theRealNeutronMass;
865           break;                                  559           break;
866         case PiPlus:                              560         case PiPlus:
867         case PiMinus:                             561         case PiMinus:
868           return theRealChargedPiMass;            562           return theRealChargedPiMass;
869           break;                                  563           break;
870         case PiZero:                              564         case PiZero:
871           return theRealPiZeroMass;               565           return theRealPiZeroMass;
872           break;                                  566           break;
873         case Eta:                                 567         case Eta:
874           return theRealEtaMass;                  568           return theRealEtaMass;
875           break;                                  569           break;
876         case Omega:                               570         case Omega:
877           return theRealOmegaMass;                571           return theRealOmegaMass;
878           break;                                  572           break;
879         case EtaPrime:                            573         case EtaPrime:
880           return theRealEtaPrimeMass;             574           return theRealEtaPrimeMass;
881           break;                                  575           break;
882         case Photon:                              576         case Photon:
883           return theRealPhotonMass;               577           return theRealPhotonMass;
884           break;                               << 
885         case Lambda:                           << 
886           return theRealLambdaMass;            << 
887           break;                               << 
888         case KPlus:                            << 
889         case KMinus:                           << 
890           return theRealChargedKaonMass;       << 
891           break;                               << 
892         case KZero:                            << 
893         case KZeroBar:                         << 
894         case KShort:                           << 
895         case KLong:                            << 
896           return theRealNeutralKaonMass;       << 
897           break;                               << 
898         case SigmaPlus:                        << 
899           return theRealSigmaPlusMass;         << 
900           break;                               << 
901         case SigmaZero:                        << 
902           return theRealSigmaZeroMass;         << 
903           break;                               << 
904         case SigmaMinus:                       << 
905           return theRealSigmaMinusMass;        << 
906           break;                               << 
907         case antiProton:                       << 
908           return theRealAntiProtonMass;        << 
909           break;                               << 
910         case XiMinus:                          << 
911           return theRealXiMinusMass;           << 
912           break;                               << 
913         case XiZero:                           << 
914           return theRealXiZeroMass;            << 
915           break;                               << 
916         case antiNeutron:                      << 
917           return theRealAntiNeutronMass;       << 
918           break;                               << 
919         case antiSigmaPlus:                    << 
920           return theRealAntiSigmaPlusMass;     << 
921           break;                                  578           break;
922         case antiSigmaZero:                    << 
923           return theRealAntiSigmaZeroMass;     << 
924           break;                               << 
925         case antiSigmaMinus:                   << 
926           return theRealAntiSigmaMinusMass;    << 
927           break;                               << 
928         case antiXiMinus:                      << 
929           return theRealAntiXiMinusMass;       << 
930           break;                               << 
931         case antiXiZero:                       << 
932           return theRealAntiXiZeroMass;        << 
933           break;                               << 
934         case antiLambda:                       << 
935           return theRealAntiLambdaMass;        << 
936           break;                               << 
937         default:                                  579         default:
938           INCL_ERROR("Particle::getRealMass :     580           INCL_ERROR("Particle::getRealMass : Unknown particle type." << '\n');
939           return 0.0;                             581           return 0.0;
940           break;                                  582           break;
941       }                                           583       }
942     }                                             584     }
943                                                << 585 
944     G4double getRealMass(const G4int A, const  << 586     G4double getRealMass(const G4int A, const G4int Z) {
945 // assert(A>=0);                                  587 // assert(A>=0);
946       // For nuclei with Z<0 or Z>A, assume th    588       // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
947       if(Z<0 && S<0)                           << 
948         return (A+S)*theRealNeutronMass - S*La << 
949       else if(Z>A && S<0)                      << 
950         return (A+S)*theRealProtonMass - S*Lam << 
951       if(Z<0)                                     589       if(Z<0)
952         return (A)*theRealNeutronMass - Z*getR << 590         return A*neutronMass - Z*getRealMass(PiMinus);
953       else if(Z>A)                                591       else if(Z>A)
954         return (A)*theRealProtonMass + (A-Z)*g << 592         return A*protonMass + (A-Z)*getRealMass(PiPlus);
955       else if(Z==0 && S==0)                    << 593       else if(Z==0)
956         return A*theRealNeutronMass;           << 594         return A*getRealMass(Neutron);
957       else if(A==Z)                               595       else if(A==Z)
958         return A*theRealProtonMass;            << 596         return A*getRealMass(Proton);
959       else if(Z==0 && S<0)                     << 
960         return (A+S)*theRealNeutronMass-S*Lamb << 
961       else if(A>1) {                              597       else if(A>1) {
962 #ifndef INCLXX_IN_GEANT4_MODE                     598 #ifndef INCLXX_IN_GEANT4_MODE
963         return ::G4INCL::NuclearMassTable::get << 599         return ::G4INCL::NuclearMassTable::getMass(A,Z);
964 #else                                             600 #else
965         if(S<0) return theG4IonTable->GetNucle << 601         return theG4IonTable->GetNucleusMass(Z,A) / MeV;
966         else    return theG4IonTable->GetNucle << 
967 #endif                                            602 #endif
968       } else                                      603       } else
969         return 0.;                                604         return 0.;
970     }                                             605     }
971                                                   606 
972     G4double getINCLMass(const G4int A, const  << 607     G4double getINCLMass(const G4int A, const G4int Z) {
973 // assert(A>=0);                                  608 // assert(A>=0);
974       // For nuclei with Z<0 or Z>A, assume th    609       // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
975       // Note that S<0 for lambda              << 610       if(Z<0)
976       if(Z<0 && S<0)                           << 611         return A*neutronMass - Z*getINCLMass(PiMinus);
977         return (A+S)*neutronMass - S*LambdaMas << 
978       else if(Z>A && S<0)                      << 
979         return (A+S)*protonMass - S*LambdaMass << 
980       else if(Z<0)                             << 
981         return (A)*neutronMass - Z*getINCLMass << 
982       else if(Z>A)                                612       else if(Z>A)
983         return (A)*protonMass + (A-Z)*getINCLM << 613         return A*protonMass + (A-Z)*getINCLMass(PiPlus);
984       else if(A>1 && S<0)                      << 
985         return Z*(protonMass - protonSeparatio << 
986       else if(A>1)                                614       else if(A>1)
987         return Z*(protonMass - protonSeparatio    615         return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy);
988       else if(A==1 && Z==0 && S==0)            << 616       else if(A==1 && Z==0)
989         return getINCLMass(Neutron);              617         return getINCLMass(Neutron);
990       else if(A==1 && Z==1 && S==0)            << 618       else if(A==1 && Z==1)
991         return getINCLMass(Proton);               619         return getINCLMass(Proton);
992       else if(A==1 && Z==0 && S==-1)           << 
993         return getINCLMass(Lambda);            << 
994       else                                        620       else
995         return 0.;                                621         return 0.;
996     }                                             622     }
997                                                   623 
998     G4double getTableQValue(const G4int A1, co << 624     G4double getTableQValue(const G4int A1, const G4int Z1, const G4int A2, const G4int Z2) {
999       return getTableMass(A1,Z1,S1) + getTable << 625       return getTableMass(A1,Z1) + getTableMass(A2,Z2) - getTableMass(A1+A2,Z1+Z2);
1000     }                                            626     }
1001                                                  627 
1002     G4double getTableQValue(const G4int A1, c << 628     G4double getTableQValue(const G4int A1, const G4int Z1, const G4int A2, const G4int Z2, const G4int A3, const G4int Z3) {
1003       return getTableMass(A1,Z1,S1) + getTabl << 629       return getTableMass(A1,Z1) + getTableMass(A2,Z2) - getTableMass(A3,Z3) - getTableMass(A1+A2-A3,Z1+Z2-Z3);
1004     }                                            630     }
1005                                                  631 
1006     G4double getTableSpeciesMass(const Partic    632     G4double getTableSpeciesMass(const ParticleSpecies &p) {
1007       if(p.theType == Composite)                 633       if(p.theType == Composite)
1008         return (*getTableMass)(p.theA, p.theZ << 634         return (*getTableMass)(p.theA, p.theZ);
1009       else                                       635       else
1010         return (*getTableParticleMass)(p.theT    636         return (*getTableParticleMass)(p.theType);
1011     }                                            637     }
1012                                                  638 
1013     G4int getMassNumber(const ParticleType t) << 639     G4int getMassNumber(const ParticleType t) {
1014                                               << 
1015       switch(t) {                                640       switch(t) {
1016         case Proton:                             641         case Proton:
1017         case Neutron:                            642         case Neutron:
1018         case DeltaPlusPlus:                      643         case DeltaPlusPlus:
1019         case DeltaPlus:                          644         case DeltaPlus:
1020         case DeltaZero:                          645         case DeltaZero:
1021         case DeltaMinus:                         646         case DeltaMinus:
1022         case SigmaPlus:                       << 
1023         case SigmaZero:                       << 
1024         case SigmaMinus:                      << 
1025         case Lambda:                          << 
1026         case XiZero:                          << 
1027         case XiMinus:                         << 
1028           return 1;                              647           return 1;
1029           break;                              << 
1030         case antiProton:                      << 
1031         case antiNeutron:                     << 
1032         case antiSigmaPlus:                   << 
1033         case antiSigmaZero:                   << 
1034         case antiSigmaMinus:                  << 
1035         case antiLambda:                      << 
1036         case antiXiZero:                      << 
1037         case antiXiMinus:                     << 
1038           return -1;                          << 
1039           break;                                 648           break;
1040         case PiPlus:                             649         case PiPlus:
1041         case PiMinus:                            650         case PiMinus:
1042         case PiZero:                             651         case PiZero:
1043         case KPlus:                           << 
1044         case KZero:                           << 
1045         case KZeroBar:                        << 
1046         case KShort:                          << 
1047         case KLong:                           << 
1048         case KMinus:                          << 
1049         case Eta:                                652         case Eta:
1050         case Omega:                              653         case Omega:
1051         case EtaPrime:                           654         case EtaPrime:
1052         case Photon:                             655         case Photon:
1053           return 0;                              656           return 0;
1054           break;                                 657           break;
1055         default:                                 658         default:
1056           return 0;                              659           return 0;
1057           break;                                 660           break;
1058       }                                          661       }
1059     }                                            662     }
1060                                                  663 
1061     G4int getChargeNumber(const ParticleType     664     G4int getChargeNumber(const ParticleType t) {
1062       switch(t) {                                665       switch(t) {
1063         case DeltaPlusPlus:                      666         case DeltaPlusPlus:
1064           return 2;                              667           return 2;
1065           break;                                 668           break;
1066         case Proton:                             669         case Proton:
1067         case DeltaPlus:                          670         case DeltaPlus:
1068         case PiPlus:                             671         case PiPlus:
1069         case SigmaPlus:                       << 
1070         case KPlus:                           << 
1071         case antiSigmaMinus:                  << 
1072         case antiXiMinus:                     << 
1073           return 1;                              672           return 1;
1074           break;                              << 673           break;
1075         case Neutron:                            674         case Neutron:
1076         case DeltaZero:                          675         case DeltaZero:
1077         case PiZero:                             676         case PiZero:
1078         case SigmaZero:                       << 
1079         case Lambda:                          << 
1080         case KZero:                           << 
1081         case KZeroBar:                        << 
1082         case KShort:                          << 
1083         case KLong:                           << 
1084         case Eta:                                677         case Eta:
1085         case Omega:                              678         case Omega:
1086         case EtaPrime:                           679         case EtaPrime:
1087         case Photon:                             680         case Photon:
1088         case XiZero:                          << 
1089         case antiNeutron:                     << 
1090         case antiLambda:                      << 
1091         case antiSigmaZero:                   << 
1092         case antiXiZero:                      << 
1093           return 0;                              681           return 0;
1094           break;                                 682           break;
1095         case DeltaMinus:                         683         case DeltaMinus:
1096         case PiMinus:                            684         case PiMinus:
1097         case SigmaMinus:                      << 
1098         case KMinus:                          << 
1099         case antiProton:                      << 
1100         case XiMinus:                         << 
1101         case antiSigmaPlus:                   << 
1102           return -1;                             685           return -1;
1103           break;                              << 
1104         default:                              << 
1105           return 0;                           << 
1106           break;                                 686           break;
1107       }                                       << 
1108     }                                         << 
1109                                               << 
1110     G4int getStrangenessNumber(const Particle << 
1111       switch(t) {                             << 
1112         case DeltaPlusPlus:                   << 
1113         case DeltaPlus:                       << 
1114         case DeltaZero:                       << 
1115         case DeltaMinus:                      << 
1116         case Proton:                          << 
1117         case Neutron:                         << 
1118         case PiPlus:                          << 
1119         case PiZero:                          << 
1120         case PiMinus:                         << 
1121         case Eta:                             << 
1122         case Omega:                           << 
1123         case EtaPrime:                        << 
1124         case Photon:                          << 
1125         case antiProton:                      << 
1126         case antiNeutron:                     << 
1127           return 0;                           << 
1128           break;                              << 
1129         case XiMinus:                         << 
1130         case XiZero:                          << 
1131         case antiXiMinus:                     << 
1132         case antiXiZero:                      << 
1133           return 2;                           << 
1134           break;                              << 
1135         case antiLambda:                      << 
1136         case antiSigmaPlus:                   << 
1137         case antiSigmaZero:                   << 
1138         case antiSigmaMinus:                  << 
1139           return 1;                           << 
1140           break;                              << 
1141         case Lambda:                          << 
1142         case SigmaPlus:                       << 
1143         case SigmaZero:                       << 
1144         case SigmaMinus:                      << 
1145         case KZeroBar:                        << 
1146         case KMinus:                          << 
1147           return -1;                          << 
1148           break;                              << 
1149         case KPlus:                           << 
1150         case KZero:                           << 
1151           return 1;                           << 
1152           break;                              << 
1153         case KShort:                          << 
1154           return 0;                           << 
1155           break;                              << 
1156         case KLong:                           << 
1157           return 0;                           << 
1158           break;                              << 
1159         default:                                 687         default:
1160           return 0;                              688           return 0;
1161           break;                                 689           break;
1162       }                                          690       }
1163     }                                            691     }
1164                                                  692 
1165     G4double getNuclearRadius(const ParticleT    693     G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z) {
1166 // assert(A>=0);                                 694 // assert(A>=0);
1167       if(A > 19 || (A < 6 && A >= 2)) {       << 695       if(A >= 19 || (A < 6 && A >= 2)) {
1168         // For large (Woods-Saxon or Modified    696         // For large (Woods-Saxon or Modified Harmonic Oscillator) or small
1169         // (Gaussian) nuclei, the radius para    697         // (Gaussian) nuclei, the radius parameter is just the nuclear radius
1170         return getRadiusParameter(t,A,Z);        698         return getRadiusParameter(t,A,Z);
1171       } else if(A < clusterTableASize && Z>=0    699       } else if(A < clusterTableASize && Z>=0 && Z < clusterTableZSize && A >= 6) {
1172         const G4double thisRMS = positionRMS[    700         const G4double thisRMS = positionRMS[Z][A];
1173         if(thisRMS>0.0)                          701         if(thisRMS>0.0)
1174           return thisRMS;                        702           return thisRMS;
1175         else {                                   703         else {
1176           INCL_DEBUG("getNuclearRadius: Radiu    704           INCL_DEBUG("getNuclearRadius: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
1177                      << "returning radius for    705                      << "returning radius for C12");
1178           return positionRMS[6][12];             706           return positionRMS[6][12];
1179         }                                        707         }
1180       } else if(A <= 19) {                    << 708       } else if(A < 19) {
1181         const G4double theRadiusParameter = g    709         const G4double theRadiusParameter = getRadiusParameter(t, A, Z);
1182         const G4double theDiffusenessParamete    710         const G4double theDiffusenessParameter = getSurfaceDiffuseness(t, A, Z);
1183         // The formula yields the nuclear RMS    711         // The formula yields the nuclear RMS radius based on the parameters of
1184         // the nuclear-density function          712         // the nuclear-density function
1185         return 1.225*theDiffusenessParameter*    713         return 1.225*theDiffusenessParameter*
1186           std::sqrt((2.+5.*theRadiusParameter    714           std::sqrt((2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter));
1187       } else {                                   715       } else {
1188         INCL_ERROR("getNuclearRadius: No radi    716         INCL_ERROR("getNuclearRadius: No radius for nucleus A = " << A << " Z = " << Z << '\n');
1189         return 0.0;                              717         return 0.0;
1190       }                                          718       }
1191     }                                            719     }
1192                                                  720 
1193     G4double getLargestNuclearRadius(const G4    721     G4double getLargestNuclearRadius(const G4int A, const G4int Z) {
1194       return Math::max(getNuclearRadius(Proto    722       return Math::max(getNuclearRadius(Proton, A, Z), getNuclearRadius(Neutron, A, Z));
1195     }                                            723     }
1196                                                  724 
1197     G4double getRadiusParameter(const Particl    725     G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z) {
1198 // assert(A>0);                                  726 // assert(A>0);
1199       if(A > 19) {                            << 727       if(A >= 28) {
1200         // radius fit for lambdas             << 
1201         if(t==Lambda){                        << 
1202          G4double r0 = (1.128+0.439*std::pow( << 
1203          return r0;                           << 
1204         }                                     << 
1205         // phenomenological radius fit           728         // phenomenological radius fit
1206         G4double r0 = (2.745e-4 * A + 1.063)     729         G4double r0 = (2.745e-4 * A + 1.063) * std::pow(A, 1.0/3.0);
1207         // HFB calculations                   << 
1208         if(getRPCorrelationCoefficient(t)<1.) << 
1209          G4double r0hfb = HFB::getRadiusParam << 
1210          if(r0hfb>0.)r0 = r0hfb;              << 
1211         }                                     << 
1212         //                                    << 
1213         if(t==Neutron)                           730         if(t==Neutron)
1214           r0 += neutronSkin;                     731           r0 += neutronSkin;
1215         return r0;                               732         return r0;
1216       } else if(A < 6 && A >= 2) {               733       } else if(A < 6 && A >= 2) {
1217         if(Z<clusterTableZSize && Z>=0) {        734         if(Z<clusterTableZSize && Z>=0) {
1218           const G4double thisRMS = positionRM    735           const G4double thisRMS = positionRMS[Z][A];
1219           if(thisRMS>0.0)                        736           if(thisRMS>0.0)
1220             return thisRMS;                      737             return thisRMS;
1221           else {                                 738           else {
1222             INCL_DEBUG("getRadiusParameter: R    739             INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
1223                        << "returning radius f    740                        << "returning radius for C12");
1224             return positionRMS[6][12];           741             return positionRMS[6][12];
1225           }                                      742           }
1226         } else {                                 743         } else {
1227           INCL_DEBUG("getRadiusParameter: Rad    744           INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
1228                      << "returning radius for    745                      << "returning radius for C12");
1229           return positionRMS[6][12];             746           return positionRMS[6][12];
1230         }                                        747         }
1231       } else if(A <= 19 && A >= 6) {          << 748       } else if(A < 28 && A >= 6) {
1232         if(t==Lambda){                        << 
1233          G4double r0 = (1.128+0.439*std::pow( << 
1234          return r0;                           << 
1235         }                                     << 
1236         // HFB calculations                   << 
1237         if(getRPCorrelationCoefficient(t)<1.) << 
1238          G4double r0hfb = HFB::getSurfaceDiff << 
1239          if(r0hfb>0.)return r0hfb;            << 
1240         }                                     << 
1241         return mediumRadius[A-1];                749         return mediumRadius[A-1];
1242         //      return 1.581*mediumDiffusenes    750         //      return 1.581*mediumDiffuseness[A-1]*(2.+5.*mediumRadius[A-1])/(2.+3.*mediumRadius[A-1]);
1243       } else {                                   751       } else {
1244         INCL_ERROR("getRadiusParameter: No ra    752         INCL_ERROR("getRadiusParameter: No radius for nucleus A = " << A << " Z = " << Z << '\n');
1245         return 0.0;                              753         return 0.0;
1246       }                                          754       }
1247     }                                            755     }
1248                                                  756 
1249     G4double getMaximumNuclearRadius(const Pa    757     G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z) {
1250       const G4double XFOISA = 8.0;               758       const G4double XFOISA = 8.0;
1251       if(A > 19) {                            << 759       if(A >= 19) {
1252         return getNuclearRadius(t,A,Z) + XFOI    760         return getNuclearRadius(t,A,Z) + XFOISA * getSurfaceDiffuseness(t,A,Z);
1253       } else if(A <= 19 && A >= 6) {          << 761       } else if(A < 19 && A >= 6) {
1254         return 5.5 + 0.3 * (G4double(A) - 6.0    762         return 5.5 + 0.3 * (G4double(A) - 6.0)/12.0;
1255       } else if(A >= 2) {                        763       } else if(A >= 2) {
1256         return getNuclearRadius(t, A, Z) + 4.    764         return getNuclearRadius(t, A, Z) + 4.5;
1257       } else {                                   765       } else {
1258         INCL_ERROR("getMaximumNuclearRadius :    766         INCL_ERROR("getMaximumNuclearRadius : No maximum radius for nucleus A = " << A << " Z = " << Z << '\n');
1259         return 0.0;                              767         return 0.0;
1260       }                                          768       }
1261     }                                            769     }
1262                                                  770 
1263     G4double getSurfaceDiffuseness(const Part    771     G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z) {
1264       if(A > 19) {                            << 772       if(A >= 28) {
1265         // phenomenological fit               << 
1266         G4double a = 1.63e-4 * A + 0.510;        773         G4double a = 1.63e-4 * A + 0.510;
1267         // HFB calculations                   << 
1268         if(getRPCorrelationCoefficient(t)<1.) << 
1269           G4double ahfb = HFB::getSurfaceDiff << 
1270           if(ahfb>0.)a=ahfb;                  << 
1271         }                                     << 
1272         //                                    << 
1273         if(t==Lambda){                        << 
1274         // Like for neutrons                  << 
1275           G4double ahfb = HFB::getSurfaceDiff << 
1276           if(ahfb>0.)a=ahfb;                  << 
1277         }                                     << 
1278         if(t==Neutron)                           774         if(t==Neutron)
1279           a += neutronHalo;                      775           a += neutronHalo;
1280         return a;                                776         return a;
1281       } else if(A <= 19 && A >= 6) {          << 777       } else if(A < 28 && A >= 19) {
1282         // HFB calculations                   << 778         return mediumDiffuseness[A-1];
1283         if(getRPCorrelationCoefficient(t)<1.) << 779       } else if(A < 19 && A >= 6) {
1284           G4double ahfb = HFB::getRadiusParam << 
1285           if(ahfb>0.)return ahfb;             << 
1286         }                                     << 
1287         return mediumDiffuseness[A-1];           780         return mediumDiffuseness[A-1];
1288       } else if(A < 6 && A >= 2) {               781       } else if(A < 6 && A >= 2) {
1289         INCL_ERROR("getSurfaceDiffuseness: wa    782         INCL_ERROR("getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << '\n');
1290         return 0.0;                              783         return 0.0;
1291       } else {                                   784       } else {
1292         INCL_ERROR("getSurfaceDiffuseness: No    785         INCL_ERROR("getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << '\n');
1293         return 0.0;                              786         return 0.0;
1294       }                                          787       }
1295     }                                            788     }
1296                                                  789 
1297     G4double getMomentumRMS(const G4int A, co    790     G4double getMomentumRMS(const G4int A, const G4int Z) {
1298 // assert(Z>=0 && A>=0 && Z<=A);                 791 // assert(Z>=0 && A>=0 && Z<=A);
1299       return getFermiMomentum(A,Z) * Math::sq    792       return getFermiMomentum(A,Z) * Math::sqrtThreeFifths;
1300     }                                            793     }
1301                                                  794 
1302     G4double getSeparationEnergyINCL(const Pa    795     G4double getSeparationEnergyINCL(const ParticleType t, const G4int /*A*/, const G4int /*Z*/) {
1303       if(t==Proton)                              796       if(t==Proton)
1304         return theINCLProtonSeparationEnergy;    797         return theINCLProtonSeparationEnergy;
1305       else if(t==Neutron)                        798       else if(t==Neutron)
1306         return theINCLNeutronSeparationEnergy    799         return theINCLNeutronSeparationEnergy;
1307       else if(t==Lambda)                      << 
1308         return theINCLLambdaSeparationEnergy; << 
1309       else if(t==antiProton)                  << 
1310         return theINCLantiProtonSeparationEne << 
1311       else {                                     800       else {
1312         INCL_ERROR("ParticleTable::getSeparat    801         INCL_ERROR("ParticleTable::getSeparationEnergyINCL : Unknown particle type." << '\n');
1313         return 0.0;                              802         return 0.0;
1314       }                                          803       }
1315     }                                            804     }
1316                                                  805 
1317     G4double getSeparationEnergyReal(const Pa    806     G4double getSeparationEnergyReal(const ParticleType t, const G4int A, const G4int Z) {
1318       // Real separation energies for all nuc    807       // Real separation energies for all nuclei
1319       if(t==Proton)                              808       if(t==Proton)
1320         return (*getTableParticleMass)(Proton << 809         return (*getTableParticleMass)(Proton) + (*getTableMass)(A-1,Z-1) - (*getTableMass)(A,Z);
1321       else if(t==Neutron)                        810       else if(t==Neutron)
1322         return (*getTableParticleMass)(Neutro << 811         return (*getTableParticleMass)(Neutron) + (*getTableMass)(A-1,Z) - (*getTableMass)(A,Z);
1323       else if(t==Lambda)                      << 
1324         return (*getTableParticleMass)(Lambda << 
1325       else {                                     812       else {
1326         INCL_ERROR("ParticleTable::getSeparat    813         INCL_ERROR("ParticleTable::getSeparationEnergyReal : Unknown particle type." << '\n');
1327         return 0.0;                              814         return 0.0;
1328       }                                          815       }
1329     }                                            816     }
1330                                                  817 
1331     G4double getSeparationEnergyRealForLight(    818     G4double getSeparationEnergyRealForLight(const ParticleType t, const G4int A, const G4int Z) {
1332       // Real separation energies for light n    819       // Real separation energies for light nuclei, fixed values for heavy nuclei
1333       if(Z<clusterTableZSize && A<clusterTabl    820       if(Z<clusterTableZSize && A<clusterTableASize)
1334         return getSeparationEnergyReal(t, A,     821         return getSeparationEnergyReal(t, A, Z);
1335       else                                       822       else
1336         return getSeparationEnergyINCL(t, A,     823         return getSeparationEnergyINCL(t, A, Z);
1337     }                                            824     }
1338                                                  825 
1339     G4double getProtonSeparationEnergy() { re    826     G4double getProtonSeparationEnergy() { return protonSeparationEnergy; }
1340                                                  827 
1341     G4double getNeutronSeparationEnergy() { r    828     G4double getNeutronSeparationEnergy() { return neutronSeparationEnergy; }
1342                                                  829 
1343     G4double getLambdaSeparationEnergy() { re << 830     void setProtonSeparationEnergy(const G4double s) { protonSeparationEnergy = s; }
1344                                               << 
1345     void setProtonSeparationEnergy(const G4do << 
1346                                                  831 
1347     void setNeutronSeparationEnergy(const G4d << 832     void setNeutronSeparationEnergy(const G4double s) { neutronSeparationEnergy  = s; }
1348                                               << 
1349     void setLambdaSeparationEnergy(const G4do << 
1350                                                  833 
1351     std::string getElementName(const G4int Z)    834     std::string getElementName(const G4int Z) {
1352       if(Z<1) {                                  835       if(Z<1) {
1353         INCL_WARN("getElementName called with    836         INCL_WARN("getElementName called with Z<1" << '\n');
1354         return elementTable[0];                  837         return elementTable[0];
1355       } else if(Z<elementTableSize)              838       } else if(Z<elementTableSize)
1356         return elementTable[Z];                  839         return elementTable[Z];
1357       else                                       840       else
1358         return getIUPACElementName(Z);           841         return getIUPACElementName(Z);
1359     }                                            842     }
1360                                                  843 
1361     std::string getIUPACElementName(const G4i    844     std::string getIUPACElementName(const G4int Z) {
1362       std::stringstream elementStream;           845       std::stringstream elementStream;
1363       elementStream << Z;                        846       elementStream << Z;
1364       std::string elementName = elementStream    847       std::string elementName = elementStream.str();
1365       std::transform(elementName.begin(), ele    848       std::transform(elementName.begin(), elementName.end(), elementName.begin(), intToIUPAC);
1366       elementName[0] = (char)std::toupper(ele << 849       elementName[0] = std::toupper(elementName.at(0));
1367       return elementName;                        850       return elementName;
1368     }                                            851     }
1369                                                  852 
1370     G4int parseElement(std::string pS) {         853     G4int parseElement(std::string pS) {
1371       // Normalize the element name              854       // Normalize the element name
1372       std::transform(pS.begin(), pS.end(), pS    855       std::transform(pS.begin(), pS.end(), pS.begin(), ::tolower);
1373       pS[0] = (char)std::toupper(pS[0]);      << 856       pS[0] = ::toupper(pS[0]);
1374                                                  857 
1375       const std::string *iter = std::find(ele    858       const std::string *iter = std::find(elementTable, elementTable+elementTableSize, pS);
1376       if(iter != elementTable+elementTableSiz    859       if(iter != elementTable+elementTableSize)
1377         return G4int(iter - elementTable);    << 860         return iter - elementTable;
1378       else                                       861       else
1379         return ParticleTable::parseIUPACEleme    862         return ParticleTable::parseIUPACElement(pS);
1380     }                                            863     }
1381                                                  864 
1382     G4int parseIUPACElement(std::string const << 865     G4int parseIUPACElement(std::string const &s) {
1383       // Normalise to lower case                 866       // Normalise to lower case
1384       std::string elementName(sel);           << 867       std::string elementName(s);
1385       std::transform(elementName.begin(), ele    868       std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower);
1386       // Return 0 if the element name contain    869       // Return 0 if the element name contains anything but IUPAC digits
1387       if(elementName.find_first_not_of(elemen    870       if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos)
1388         return 0;                                871         return 0;
1389       std::transform(elementName.begin(), ele    872       std::transform(elementName.begin(), elementName.end(), elementName.begin(), iupacToInt);
1390       std::stringstream elementStream(element    873       std::stringstream elementStream(elementName);
1391       G4int Z;                                   874       G4int Z;
1392       elementStream >> Z;                        875       elementStream >> Z;
1393       return Z;                                  876       return Z;
1394     }                                            877     }
1395                                                  878 
1396     IsotopicDistribution const &getNaturalIso    879     IsotopicDistribution const &getNaturalIsotopicDistribution(const G4int Z) {
1397       return getNaturalIsotopicDistributions(    880       return getNaturalIsotopicDistributions()->getIsotopicDistribution(Z);
1398     }                                            881     }
1399                                                  882 
1400     G4int drawRandomNaturalIsotope(const G4in    883     G4int drawRandomNaturalIsotope(const G4int Z) {
1401       return getNaturalIsotopicDistributions(    884       return getNaturalIsotopicDistributions()->drawRandomIsotope(Z);
1402     }                                            885     }
1403                                                  886 
1404     G4double getFermiMomentumConstant(const G    887     G4double getFermiMomentumConstant(const G4int /*A*/, const G4int /*Z*/) {
1405       return constantFermiMomentum;              888       return constantFermiMomentum;
1406     }                                            889     }
1407                                                  890 
1408     G4double getFermiMomentumConstantLight(co    891     G4double getFermiMomentumConstantLight(const G4int A, const G4int Z) {
1409 // assert(Z>0 && A>0 && Z<=A);                   892 // assert(Z>0 && A>0 && Z<=A);
1410       if(Z<clusterTableZSize && A<clusterTabl    893       if(Z<clusterTableZSize && A<clusterTableASize) {
1411         const G4double rms = momentumRMS[Z][A    894         const G4double rms = momentumRMS[Z][A];
1412         return ((rms>0.) ? rms : momentumRMS[    895         return ((rms>0.) ? rms : momentumRMS[6][12]) * Math::sqrtFiveThirds;
1413       } else                                     896       } else
1414         return getFermiMomentumConstant(A,Z);    897         return getFermiMomentumConstant(A,Z);
1415     }                                            898     }
1416                                                  899 
1417     G4double getFermiMomentumMassDependent(co    900     G4double getFermiMomentumMassDependent(const G4int A, const G4int /*Z*/) {
1418 // assert(A>0);                                  901 // assert(A>0);
1419       static const G4double alphaParam = 259.    902       static const G4double alphaParam = 259.416; // MeV/c
1420       static const G4double betaParam  = 152.    903       static const G4double betaParam  = 152.824; // MeV/c
1421       static const G4double gammaParam = 9.51    904       static const G4double gammaParam = 9.5157E-2;
1422       return alphaParam - betaParam*std::exp(    905       return alphaParam - betaParam*std::exp(-gammaParam*((G4double)A));
1423     }                                            906     }
1424                                                  907 
1425     G4double getRPCorrelationCoefficient(cons    908     G4double getRPCorrelationCoefficient(const ParticleType t) {
1426 // assert(t==Proton || t==Neutron || t==Lambd << 909 // assert(t==Proton || t==Neutron);
1427       return rpCorrelationCoefficient[t];        910       return rpCorrelationCoefficient[t];
1428     }                                            911     }
1429                                                  912 
1430     G4double getNeutronSkin() { return neutro    913     G4double getNeutronSkin() { return neutronSkin; }
1431                                                  914 
1432     G4double getNeutronHalo() { return neutro    915     G4double getNeutronHalo() { return neutronHalo; }
1433                                                  916 
1434     G4ThreadLocal G4double minDeltaMass = 0.;    917     G4ThreadLocal G4double minDeltaMass = 0.;
1435     G4ThreadLocal G4double minDeltaMass2 = 0.    918     G4ThreadLocal G4double minDeltaMass2 = 0.;
1436     G4ThreadLocal G4double minDeltaMassRndm =    919     G4ThreadLocal G4double minDeltaMassRndm = 0.;
1437     G4ThreadLocal NuclearMassFn getTableMass     920     G4ThreadLocal NuclearMassFn getTableMass = NULL;
1438     G4ThreadLocal ParticleMassFn getTablePart    921     G4ThreadLocal ParticleMassFn getTableParticleMass = NULL;
1439     G4ThreadLocal SeparationEnergyFn getSepar    922     G4ThreadLocal SeparationEnergyFn getSeparationEnergy = NULL;
1440     G4ThreadLocal FermiMomentumFn getFermiMom    923     G4ThreadLocal FermiMomentumFn getFermiMomentum = NULL;
1441                                                  924 
1442     ParticleType getPionType(const G4int isos    925     ParticleType getPionType(const G4int isosp) {
1443 // assert(isosp == -2 || isosp == 0 || isosp     926 // assert(isosp == -2 || isosp == 0 || isosp == 2);
1444         if (isosp == -2) {                       927         if (isosp == -2) {
1445             return PiMinus;                      928             return PiMinus;
1446         }                                        929         }
1447         else if (isosp == 0) {                   930         else if (isosp == 0) {
1448             return PiZero;                       931             return PiZero;
1449         }                                        932         }
1450         else {                                   933         else {
1451             return PiPlus;                       934             return PiPlus;
1452         }                                        935         }
1453     }                                            936     }
1454                                                  937 
1455     ParticleType getNucleonType(const G4int i    938     ParticleType getNucleonType(const G4int isosp) {
1456 // assert(isosp == -1 || isosp == 1);            939 // assert(isosp == -1 || isosp == 1);
1457         if (isosp == -1) {                       940         if (isosp == -1) {
1458             return Neutron;                      941             return Neutron;
1459         }                                        942         }
1460         else {                                   943         else {
1461             return Proton;                       944             return Proton;
1462         }                                        945         }
1463     }                                            946     }
1464                                                  947 
1465     ParticleType getDeltaType(const G4int iso    948     ParticleType getDeltaType(const G4int isosp) {
1466 // assert(isosp == -3 || isosp == -1 || isosp    949 // assert(isosp == -3 || isosp == -1 || isosp == 1 || isosp == 3);
1467         if (isosp == -3) {                       950         if (isosp == -3) {
1468             return DeltaMinus;                   951             return DeltaMinus;
1469         }                                        952         }
1470         else if (isosp == -1) {                  953         else if (isosp == -1) {
1471             return DeltaZero;                    954             return DeltaZero;
1472         }                                        955         }
1473         else if (isosp == 1) {                   956         else if (isosp == 1) {
1474             return DeltaPlus;                    957             return DeltaPlus;
1475         }                                        958         }
1476         else {                                   959         else {
1477             return DeltaPlusPlus;                960             return DeltaPlusPlus;
1478         }                                        961         }
1479     }                                            962     }
1480                                                  963 
1481     ParticleType getSigmaType(const G4int iso << 964   G4double getWidth(const ParticleType pt) {
1482 // assert(isosp == -2 || isosp == 0 || isosp  << 965 // assert(pt == PiPlus || pt == PiMinus || pt == PiZero || pt == Eta || pt == Omega || pt == EtaPrime);
1483         if (isosp == -2) {                    << 966       if(pt == PiPlus) {
1484             return SigmaMinus;                << 967         return piPlusWidth;
1485         }                                     << 968       } else if(pt == PiMinus) {
1486         else if (isosp == 0) {                << 969         return piMinusWidth;
1487             return SigmaZero;                 << 970       } else if(pt == PiZero) {
1488         }                                     << 971         return piZeroWidth;
1489         else {                                << 972       } else if(pt == Eta) {
1490             return SigmaPlus;                 << 973         return etaWidth;
1491         }                                     << 974       } else if(pt == Omega) {
1492     }                                         << 975         return omegaWidth;
1493                                               << 976       } else if(pt == EtaPrime) {
1494     ParticleType getXiType(const G4int isosp) << 977         return etaPrimeWidth;
1495 // assert(isosp == -1 || isosp == 1);         << 978       } else {
1496         if (isosp == -1) {                    << 979         INCL_ERROR("getWidth : Unknown particle type." << '\n');
1497             return XiMinus;                   << 980         return 0.0;
1498         }                                     << 981       }
1499         else {                                << 982     }
1500             return XiZero;                    << 983     
1501         }                                     << 984     
1502     }                                         << 
1503                                               << 
1504   /*ParticleType getAntiNucleonType(const G4i << 
1505 // assert(isosp == -1); //|| isosp == 1       << 
1506       if (isosp == -1) {                      << 
1507           return antiProton;                  << 
1508       }                                       << 
1509       else {                                  << 
1510           return antiNeutron;                 << 
1511       }                                       << 
1512     }*/                                       << 
1513                                               << 
1514     ParticleType getAntiSigmaType(const G4int << 
1515 // assert(isosp == -2 || isosp == 0 || isosp  << 
1516         if (isosp == -2) {                    << 
1517             return antiSigmaPlus;             << 
1518         }                                     << 
1519         else if (isosp == 0) {                << 
1520             return antiSigmaZero;             << 
1521         }                                     << 
1522         else {                                << 
1523             return antiSigmaMinus;            << 
1524         }                                     << 
1525     }                                         << 
1526                                               << 
1527     ParticleType getAntiXiType(const G4int is << 
1528 // assert(isosp == -1 || isosp == 1);         << 
1529         if (isosp == -1) {                    << 
1530             return antiXiZero;                << 
1531         }                                     << 
1532         else {                                << 
1533             return antiXiMinus;               << 
1534         }                                     << 
1535     }                                         << 
1536                                               << 
1537     ParticleType getKaonType(const G4int isos << 
1538 // assert(isosp == -1 || isosp == 1);         << 
1539         if (isosp == -1) {                    << 
1540             return KZero;                     << 
1541         }                                     << 
1542         else {                                << 
1543             return KPlus;                     << 
1544         }                                     << 
1545     }                                         << 
1546                                               << 
1547     ParticleType getAntiKaonType(const G4int  << 
1548 // assert(isosp == -1 || isosp == 1);         << 
1549         if (isosp == -1) {                    << 
1550             return KMinus;                    << 
1551         }                                     << 
1552         else {                                << 
1553             return KZeroBar;                  << 
1554         }                                     << 
1555     }                                         << 
1556                                               << 
1557     G4double getWidth(const ParticleType pt)  << 
1558 // assert(pt == PiPlus || pt == PiMinus || pt << 
1559           if(pt == PiPlus) {                  << 
1560               return piPlusWidth;             << 
1561           } else if(pt == PiMinus) {          << 
1562               return piMinusWidth;            << 
1563           } else if(pt == PiZero) {           << 
1564               return piZeroWidth;             << 
1565           } else if(pt == Eta) {              << 
1566               return etaWidth;                << 
1567           } else if(pt == Omega) {            << 
1568               return omegaWidth;              << 
1569           } else if(pt == EtaPrime) {         << 
1570               return etaPrimeWidth;           << 
1571           } else if(pt == SigmaPlus) {        << 
1572               return SigmaPlusWidth;          << 
1573           } else if(pt == SigmaZero) {        << 
1574               return SigmaZeroWidth;          << 
1575           } else if(pt == SigmaMinus) {       << 
1576               return SigmaMinusWidth;         << 
1577           } else if(pt == KPlus) {            << 
1578               return KPlusWidth;              << 
1579           } else if(pt == KMinus) {           << 
1580               return KMinusWidth;             << 
1581           } else if(pt == KShort) {           << 
1582               return KShortWidth;             << 
1583           } else if(pt == KLong) {            << 
1584               return KLongWidth;              << 
1585           } else if(pt == Lambda) {           << 
1586               return LambdaWidth;             << 
1587           } else if(pt == XiMinus) {          << 
1588               return XiMinusWidth;            << 
1589           } else if(pt == XiZero) {           << 
1590               return XiZeroWidth;             << 
1591           } else if(pt == antiSigmaPlus) {    << 
1592               return antiSigmaPlusWidth;      << 
1593           } else if(pt == antiSigmaZero) {    << 
1594               return antiSigmaZeroWidth;      << 
1595           } else if(pt == antiSigmaMinus) {   << 
1596               return antiSigmaMinusWidth;     << 
1597           } else if(pt == antiLambda) {       << 
1598               return antiLambdaWidth;         << 
1599           } else if(pt == antiXiMinus) {      << 
1600               return antiXiMinusWidth;        << 
1601           } else if(pt == antiXiZero) {       << 
1602               return antiXiZeroWidth;         << 
1603           } else {                            << 
1604               INCL_ERROR("getWidth : Unknown  << 
1605               return 0.0;                     << 
1606           }                                   << 
1607       }                                       << 
1608                                               << 
1609   } // namespace ParticleTable                   985   } // namespace ParticleTable
1610 } // namespace G4INCL                            986 } // namespace G4INCL
1611                                                  987 
1612                                                  988