[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
| 27 | // $Id: G4QPDGCode.cc,v 1.55.2.1 2008/04/23 14:47:23 gcosmo Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-01-patch-02 $ |
---|
| 29 | // |
---|
| 30 | // ---------------- G4QPDGCode ---------------- |
---|
| 31 | // by Mikhail Kossov, Sept 1999. |
---|
| 32 | // class for Hadron definitions in CHIPS Model |
---|
| 33 | // ------------------------------------------------------------------- |
---|
| 34 | |
---|
| 35 | //#define debug |
---|
| 36 | //#define pdebug |
---|
| 37 | //#define qdebug |
---|
| 38 | //#define idebug |
---|
| 39 | //#define sdebug |
---|
| 40 | |
---|
| 41 | #include "G4QPDGCodeVector.hh" |
---|
| 42 | #include <cmath> |
---|
| 43 | #include <cstdlib> |
---|
| 44 | using namespace std; |
---|
| 45 | |
---|
| 46 | G4QPDGCode::G4QPDGCode(G4int PDGCode): thePDGCode(PDGCode) |
---|
| 47 | { |
---|
| 48 | #ifdef sdebug |
---|
| 49 | G4cout<<"G4QPDGCode:Constructer is called with PDGCode="<<PDGCode<<G4endl; |
---|
| 50 | #endif |
---|
| 51 | if(PDGCode) theQCode=MakeQCode(PDGCode); |
---|
| 52 | else |
---|
| 53 | { |
---|
| 54 | #ifdef sdebug |
---|
| 55 | G4cout<<"***G4QPDGCode: Constructed with PDGCode=0, QCode=-2"<<G4endl; |
---|
| 56 | #endif |
---|
| 57 | theQCode=-2; |
---|
| 58 | } |
---|
| 59 | #ifdef debug |
---|
| 60 | G4cout<<"G4QPDGCode:Constructer(PDG) PDG="<<PDGCode<<", QCode="<<theQCode<<G4endl; |
---|
| 61 | #endif |
---|
| 62 | } |
---|
| 63 | |
---|
| 64 | G4QPDGCode::G4QPDGCode(G4bool f, G4int QCode): theQCode(QCode) |
---|
| 65 | { |
---|
| 66 | if(f&&QCode<0)G4cerr<<"***G4QPDGCode::Constr. QCode="<<QCode<<G4endl; |
---|
| 67 | thePDGCode = MakePDGCode(QCode); |
---|
| 68 | #ifdef debug |
---|
| 69 | G4cout<<"G4QPDGCode::Constr: PDGCode="<<thePDGCode<<G4endl; |
---|
| 70 | #endif |
---|
| 71 | } |
---|
| 72 | |
---|
| 73 | G4QPDGCode::G4QPDGCode(G4QContent QCont) |
---|
| 74 | { |
---|
| 75 | InitByQCont(QCont); |
---|
| 76 | } |
---|
| 77 | |
---|
| 78 | G4QPDGCode::G4QPDGCode(const G4QPDGCode& rhs) |
---|
| 79 | { |
---|
| 80 | thePDGCode =rhs.thePDGCode; |
---|
| 81 | theQCode =rhs.theQCode; |
---|
| 82 | } |
---|
| 83 | |
---|
| 84 | G4QPDGCode::G4QPDGCode(G4QPDGCode* rhs) |
---|
| 85 | { |
---|
| 86 | thePDGCode =rhs->thePDGCode; |
---|
| 87 | theQCode =rhs->theQCode; |
---|
| 88 | } |
---|
| 89 | |
---|
| 90 | const G4QPDGCode& G4QPDGCode::operator=(const G4QPDGCode& rhs) |
---|
| 91 | { |
---|
| 92 | if(this != &rhs) // Beware of self assignment |
---|
| 93 | { |
---|
| 94 | thePDGCode =rhs.thePDGCode; |
---|
| 95 | theQCode =rhs.theQCode; |
---|
| 96 | } |
---|
| 97 | return *this; |
---|
| 98 | } |
---|
| 99 | |
---|
| 100 | G4QPDGCode::~G4QPDGCode() {} |
---|
| 101 | |
---|
| 102 | //G4int G4QPDGCode::nQHM=90; |
---|
| 103 | |
---|
| 104 | // Standard output for QPDGCode |
---|
| 105 | ostream& operator<<(ostream& lhs, G4QPDGCode& rhs) |
---|
| 106 | // ========================================= |
---|
| 107 | { |
---|
| 108 | lhs << "[ PDG=" << rhs.GetPDGCode() << ", Q=" << rhs.GetQCode() << "]"; |
---|
| 109 | return lhs; |
---|
| 110 | } |
---|
| 111 | |
---|
| 112 | // Standard output for const QPDGCode |
---|
| 113 | ostream& operator<<(ostream& lhs, const G4QPDGCode& rhs) |
---|
| 114 | // =============================================== |
---|
| 115 | { |
---|
| 116 | lhs << "[ PDG=" << rhs.GetPDGCode() << ", Q=" << rhs.GetQCode() << "]"; |
---|
| 117 | return lhs; |
---|
| 118 | } |
---|
| 119 | |
---|
| 120 | // Overloading of QPDGCode addition |
---|
| 121 | G4int operator+(const G4QPDGCode& lhs, const G4QPDGCode& rhs) |
---|
| 122 | // ======================================================= |
---|
| 123 | { |
---|
| 124 | G4int s = lhs.GetPDGCode(); |
---|
| 125 | return s += rhs.GetPDGCode(); |
---|
| 126 | } |
---|
| 127 | G4int operator+(const G4QPDGCode& lhs, const G4int& rhs) |
---|
| 128 | // ======================================================= |
---|
| 129 | { |
---|
| 130 | G4int s = lhs.GetPDGCode(); |
---|
| 131 | return s += rhs; |
---|
| 132 | } |
---|
| 133 | G4int operator+(const G4int& lhs, const G4QPDGCode& rhs) |
---|
| 134 | // ======================================================= |
---|
| 135 | { |
---|
| 136 | G4int s = lhs; |
---|
| 137 | return s += rhs.GetPDGCode(); |
---|
| 138 | } |
---|
| 139 | |
---|
| 140 | // Overloading of QPDGCode subtraction |
---|
| 141 | G4int operator-(const G4QPDGCode& lhs, const G4QPDGCode& rhs) |
---|
| 142 | // ======================================================= |
---|
| 143 | { |
---|
| 144 | G4int s = lhs.GetPDGCode(); |
---|
| 145 | return s -= rhs.GetPDGCode(); |
---|
| 146 | } |
---|
| 147 | G4int operator-(const G4QPDGCode& lhs, const G4int& rhs) |
---|
| 148 | // ======================================================= |
---|
| 149 | { |
---|
| 150 | G4int s = lhs.GetPDGCode(); |
---|
| 151 | return s -= rhs; |
---|
| 152 | } |
---|
| 153 | G4int operator-(const G4int& lhs, const G4QPDGCode& rhs) |
---|
| 154 | // ======================================================= |
---|
| 155 | { |
---|
| 156 | G4int s = lhs; |
---|
| 157 | return s -= rhs.GetPDGCode(); |
---|
| 158 | } |
---|
| 159 | |
---|
| 160 | // Overloading of QPDGCode multiplication |
---|
| 161 | G4int operator*(const G4QPDGCode& lhs, const G4QPDGCode& rhs) |
---|
| 162 | // ======================================================= |
---|
| 163 | { |
---|
| 164 | G4int s = lhs.GetPDGCode(); |
---|
| 165 | return s *= rhs.GetPDGCode(); |
---|
| 166 | } |
---|
| 167 | G4int operator*(const G4QPDGCode& lhs, const G4int& rhs) |
---|
| 168 | // ======================================================= |
---|
| 169 | { |
---|
| 170 | G4int s = lhs.GetPDGCode(); |
---|
| 171 | return s *= rhs; |
---|
| 172 | } |
---|
| 173 | G4int operator*(const G4int& lhs, const G4QPDGCode& rhs) |
---|
| 174 | // ======================================================= |
---|
| 175 | { |
---|
| 176 | G4int s = lhs; |
---|
| 177 | return s *= rhs.GetPDGCode(); |
---|
| 178 | } |
---|
| 179 | |
---|
| 180 | // Overloading of QPDGCode division |
---|
| 181 | G4int operator/(const G4QPDGCode& lhs, const G4QPDGCode& rhs) |
---|
| 182 | {// ======================================================= |
---|
| 183 | G4int s = lhs.GetPDGCode(); |
---|
| 184 | return s /= rhs.GetPDGCode(); |
---|
| 185 | } |
---|
| 186 | G4int operator/(const G4QPDGCode& lhs, const G4int& rhs) |
---|
| 187 | {// ================================================== |
---|
| 188 | G4int s = lhs.GetPDGCode(); |
---|
| 189 | return s /= rhs; |
---|
| 190 | } |
---|
| 191 | G4int operator/(const G4int& lhs, const G4QPDGCode& rhs) |
---|
| 192 | {// ================================================== |
---|
| 193 | G4int s = lhs; |
---|
| 194 | return s /= rhs.GetPDGCode(); |
---|
| 195 | } |
---|
| 196 | |
---|
| 197 | // Overloading of QPDGCode residual |
---|
| 198 | G4int operator%(const G4QPDGCode& lhs, const G4int& rhs) |
---|
| 199 | {// ================================================== |
---|
| 200 | G4int s = lhs.GetPDGCode(); |
---|
| 201 | return s %= rhs; |
---|
| 202 | } |
---|
| 203 | |
---|
| 204 | // TRUE if it is not RealNeutral (111,221,331 etc), FALSE if it is. |
---|
| 205 | G4bool G4QPDGCode::TestRealNeutral(const G4int& PDGCode) |
---|
| 206 | {// ================================================= |
---|
| 207 | if(PDGCode>0 && PDGCode<999) // RealNeutral are always positive && mesons |
---|
| 208 | { |
---|
| 209 | if(PDGCode==22) return false; // Photon |
---|
| 210 | G4int p=PDGCode/10; |
---|
| 211 | if(p/10==p%10) return false; // This is a RealNeutral |
---|
| 212 | } |
---|
| 213 | return true; |
---|
| 214 | } |
---|
| 215 | |
---|
| 216 | // Make a Q Code out of the PDG Code |
---|
| 217 | G4int G4QPDGCode::MakePDGCode(const G4int& QCode) |
---|
| 218 | {// =========================================== |
---|
| 219 | //static const G4int modi = 81; // Q Codes for more than di-baryon nuclei |
---|
| 220 | //static const G4int modi = 89; // Q Codes for more than di-baryon nuclei "IsoNuclei" |
---|
| 221 | static const G4int modi = 122; // Q Codes for more than quarta-baryon nuclei "Lept/Hyper" |
---|
| 222 | static G4int qC[modi] = { 11, 12, 13, 14, 15, 16, 22, 23, 24, 25, // 10 |
---|
| 223 | 37, 110, 220, 330, 111, 211, 221, 311, 321, 331, // 20 |
---|
| 224 | 2112, 2212, 3122, 3112, 3212, 3222, 3312, 3322, 113, 213, // 30 |
---|
| 225 | 223, 313, 323, 333, 1114, 2114, 2214, 2224, 3124, 3114, // 40 |
---|
| 226 | 3214, 3224, 3314, 3324, 3334, 115, 215, 225, 315, 325, // 50 |
---|
| 227 | 335, 2116, 2216, 3126, 3116, 3216, 3226, 3316, 3326, 117, // 60 |
---|
| 228 | 217, 227, 317, 327, 337, 1118, 2118, 2218, 2228, 3128, // 70 |
---|
| 229 | 3118, 3218, 3228, 3318, 3328, 3338, 119, 219, 229, 319, // 80 |
---|
| 230 | 329, 339, 90002999 , 89999003 , 90003998 , |
---|
| 231 | 89998004 , 90003999 , 89999004 , 90004998 , 89998005 , // 90 |
---|
| 232 | 90000001 , 90001000 , 91000000 , 90999001 , 91000999 , |
---|
| 233 | 91999000 , 91999999 , 92999000 , 90000002 , 90001001 , //100 |
---|
| 234 | 90002000 , 91000001 , 91001000 , 92000000 , 90999002 , |
---|
| 235 | 91001999 , 90001002 , 90002001 , 91000002 , 91001001 , //110 |
---|
| 236 | 91002000 , 92000001 , 92001000 , 90999003 , 90001003 , |
---|
| 237 | 90002002 , 90003001 , 91001002 , 91002001 , 92000002 , //120 |
---|
| 238 | 92001001 , 92002000}; |
---|
| 239 | static G4int aC[15] = {1,1000,999001,1000000,1000999,1999000,1999999, // sum 1 |
---|
| 240 | 2,1001,2000,1000001,1001000,1999001,2000000,2000999}; // sum 2 |
---|
| 241 | if (QCode<0) |
---|
| 242 | { |
---|
| 243 | G4cerr<<"***G4QPDGCode::MakePDGCode: negative Q Code ="<<QCode<<G4endl; |
---|
| 244 | return 0; |
---|
| 245 | } |
---|
| 246 | else if (QCode>=modi) |
---|
| 247 | { |
---|
| 248 | //G4int q=QCode-modi; // Starting BarNum=3 |
---|
| 249 | //G4int a=q/15+1; // BarNum/2 |
---|
| 250 | //G4int b=q%15; |
---|
| 251 | G4int q=QCode-modi; // Starting BarNum=5 |
---|
| 252 | G4int a=q/15+2; // BarNum/2 |
---|
| 253 | G4int b=q%15; |
---|
| 254 | return 90000000+a*1001+aC[b]; |
---|
| 255 | } |
---|
| 256 | return qC[theQCode]; |
---|
| 257 | } |
---|
| 258 | |
---|
| 259 | // Hadronic masses synhronized with the Geant4 hadronic masses |
---|
| 260 | G4double G4QPDGCode:: QHaM(G4int nQ) |
---|
| 261 | {// =========================== |
---|
| 262 | static G4bool iniFlag=true; |
---|
| 263 | static G4double m[nQHM]={.511, 0., 105.65837, 0., 1777., 0., 0., 91.188, 80.425, 140.00 |
---|
| 264 | ,120.000, 800., 980., 1370., 134.98, 139.57, 547.75, 497.65, 493.68, 957.78 |
---|
| 265 | ,939.5654,938.272, 1115.683, 1197.45, 1192.64, 1189.37, 1321.3, 1314.8, 775.8, 775.8 |
---|
| 266 | , 782.6, 896.1, 891.66, 1019.456, 1232., 1232., 1232., 1232., 1519.5, 1387.2 |
---|
| 267 | , 1383.7, 1382.8, 1535., 1531.8, 1672.45, 1318.3, 1318.3, 1275.4, 1432.4, 1425.6 |
---|
| 268 | , 1525., 1680., 1680., 1820., 1915., 1915., 1915., 2025., 2025., 1691. |
---|
| 269 | , 1691., 1667., 1776., 1776., 1854., 1950., 1950., 1950., 1950., 2100. |
---|
| 270 | , 2030., 2030., 2030., 2127., 2127., 2252., 2020., 2020., 2044., 2045. |
---|
| 271 | , 2045., 2297., 2170.272, 2171.565, 2464., 2464., 3108.544, 3111.13,3402.272,3403.565}; |
---|
| 272 | if(iniFlag) // Initialization of the Geant4 hadronic masses |
---|
| 273 | { |
---|
| 274 | m[ 0]= G4Electron::Electron()->GetPDGMass(); |
---|
| 275 | m[ 1]= G4NeutrinoE::NeutrinoE()->GetPDGMass(); |
---|
| 276 | m[ 2]= G4MuonMinus::MuonMinus()->GetPDGMass(); |
---|
| 277 | m[ 3]= G4NeutrinoMu::NeutrinoMu()->GetPDGMass(); |
---|
| 278 | m[ 4]= G4TauMinus::TauMinus()->GetPDGMass(); |
---|
| 279 | m[ 5]=G4NeutrinoTau::NeutrinoTau()->GetPDGMass(); |
---|
| 280 | m[14]= G4PionZero::PionZero()->GetPDGMass(); |
---|
| 281 | m[15]= G4PionMinus::PionMinus()->GetPDGMass(); |
---|
| 282 | m[16]= G4Eta::Eta()->GetPDGMass(); |
---|
| 283 | m[17]= G4KaonZero::KaonZero()->GetPDGMass(); |
---|
| 284 | m[18]= G4KaonMinus::KaonMinus()->GetPDGMass(); |
---|
| 285 | m[19]= G4EtaPrime::EtaPrime()->GetPDGMass(); |
---|
| 286 | m[20]= G4Neutron::Neutron()->GetPDGMass(); |
---|
| 287 | m[21]= G4Proton::Proton()->GetPDGMass(); |
---|
| 288 | m[22]= G4Lambda::Lambda()->GetPDGMass(); |
---|
| 289 | m[23]= G4SigmaMinus::SigmaMinus()->GetPDGMass(); |
---|
| 290 | m[24]= G4SigmaZero::SigmaZero()->GetPDGMass(); |
---|
| 291 | m[25]= G4SigmaPlus::SigmaPlus()->GetPDGMass(); |
---|
| 292 | m[26]= G4XiMinus::XiMinus()->GetPDGMass(); |
---|
| 293 | m[27]= G4XiZero::XiZero()->GetPDGMass(); |
---|
| 294 | m[44]= G4OmegaMinus::OmegaMinus()->GetPDGMass(); |
---|
| 295 | iniFlag=false; |
---|
| 296 | } |
---|
| 297 | if(nQ<0 || nQ>=nQHM) |
---|
| 298 | { |
---|
| 299 | G4cout<<"***G4QPDGCode::QHaM: negative Q-code or Q="<<nQ<<" >= nQmax = "<<nQHM<<G4endl; |
---|
| 300 | return 0.; |
---|
| 301 | } |
---|
| 302 | return m[nQ]; |
---|
| 303 | } |
---|
| 304 | |
---|
| 305 | // Make a Q Code out of the PDG Code |
---|
| 306 | G4int G4QPDGCode::MakeQCode(const G4int& PDGCode) |
---|
| 307 | {// =========================================== |
---|
| 308 | static const G4int qr[10]={0,13,19,27,33,44,50,58,64,75}; |
---|
| 309 | G4int PDGC=abs(PDGCode); // Qcode is always not negative |
---|
| 310 | G4int s=0; |
---|
| 311 | G4int z=0; |
---|
| 312 | G4int n=0; |
---|
| 313 | if (PDGC>100000000) // Not supported |
---|
| 314 | { |
---|
| 315 | #ifdef debug |
---|
| 316 | G4cout<<"***G4QPDGCode::MakeQCode: Unknown in Q-System code: "<<PDGCode<<G4endl; |
---|
| 317 | #endif |
---|
| 318 | return -2; |
---|
| 319 | } |
---|
| 320 | else if (PDGC>80000000&&PDGC<100000000) // Try to convert the NUCCoding to PDGCoding |
---|
| 321 | { |
---|
| 322 | if(PDGC==90000000) return 6; |
---|
| 323 | ConvertPDGToZNS(PDGC, z, n, s); |
---|
| 324 | G4int b=n+z+s; // Baryon number |
---|
| 325 | #ifdef debug |
---|
| 326 | G4cout<<"***G4QPDGCode::Z="<<z<<",N="<<n<<",S="<<s<<G4endl; |
---|
| 327 | #endif |
---|
| 328 | if(b<0) // ---> Baryons & Fragments |
---|
| 329 | { |
---|
| 330 | b=-b; |
---|
| 331 | n=-n; |
---|
| 332 | z=-z; |
---|
| 333 | s=-s; |
---|
| 334 | PDGC=90000000+s*1000000+z*1000+n; // New PDGC for anti-baryons |
---|
| 335 | } |
---|
| 336 | else if(!b) // --> Mesons |
---|
| 337 | { |
---|
| 338 | //G4bool anti=false; // For the PDG conversion |
---|
| 339 | if(z<0) // --> Mesons conversion |
---|
| 340 | { |
---|
| 341 | n=-n; |
---|
| 342 | z=-z; |
---|
| 343 | s=-s; |
---|
| 344 | //anti=true; // For the PDG conversion |
---|
| 345 | } |
---|
| 346 | if(!z) |
---|
| 347 | { |
---|
| 348 | if(s>0) |
---|
| 349 | { |
---|
| 350 | n=-n; |
---|
| 351 | s=-s; |
---|
| 352 | //anti=true; // For the PDG conversion |
---|
| 353 | } |
---|
| 354 | if (s==-1) return 17; // K0 |
---|
| 355 | else if(s==-2) return -1; // K0+K0 chipolino |
---|
| 356 | else return -2; // Not supported by Q Code |
---|
| 357 | } |
---|
| 358 | else // --> z>0 |
---|
| 359 | { |
---|
| 360 | if(z==1) |
---|
| 361 | { |
---|
| 362 | if (s==-1) return 18; // K+ |
---|
| 363 | else return 15; // pi+ |
---|
| 364 | } |
---|
| 365 | else if(z==2) return -1; // Chipolino |
---|
| 366 | else return -2; // Not supported by Q Code |
---|
| 367 | } |
---|
| 368 | } // End of meson case |
---|
| 369 | if(b>0) // --> Baryon case |
---|
| 370 | { |
---|
| 371 | if(b==1) |
---|
| 372 | { |
---|
| 373 | if(!s) // --> Baryons |
---|
| 374 | { |
---|
| 375 | if(z==-1) return 34; // Delta- |
---|
| 376 | else if(!z) return 91; // neutron |
---|
| 377 | else if(z==1)return 91; // proton |
---|
| 378 | else if(z==2)return 37; // Delta++ |
---|
| 379 | else if(z==3||z==-2)return -1; // Delta+pi Chipolino |
---|
| 380 | else return -2; // Not supported by Q Code |
---|
| 381 | } |
---|
| 382 | else if(s==1) // --> Hyperons |
---|
| 383 | { |
---|
| 384 | if(z==-1) return 93; // Sigma- |
---|
| 385 | else if(!z) return 92; // Lambda (@@ 24->Sigma0) |
---|
| 386 | else if(z==1)return 94; // Sigma+ |
---|
| 387 | else if(z==2||z==-2) return -1; // Sigma+pi Chipolino |
---|
| 388 | else return -2; // Not supported by Q Code |
---|
| 389 | } |
---|
| 390 | else if(s==2) // --> Xi Hyperons |
---|
| 391 | { |
---|
| 392 | if(z==-1) return 95; // Xi- |
---|
| 393 | else if(!z) return 96; // Xi0 |
---|
| 394 | else if(z==1||z==-2)return -1; // Xi+pi Chipolino |
---|
| 395 | else return -2; // Not supported by Q Code |
---|
| 396 | } |
---|
| 397 | else if(s==3) // --> Xi Hyperons |
---|
| 398 | { |
---|
| 399 | if(z==-1) return 97; // Omega- |
---|
| 400 | else if(!z||z==-2) return -1; // Omega+pi Chipolino |
---|
| 401 | else return -2; // Not supported by Q Code |
---|
| 402 | } |
---|
| 403 | } |
---|
| 404 | else |
---|
| 405 | { |
---|
| 406 | if(b==2) |
---|
| 407 | { |
---|
| 408 | if (PDGC==90002999) return 82; // p DEL++ |
---|
| 409 | else if(PDGC==89999003) return 83; // n DEL- |
---|
| 410 | else if(PDGC==90003998) return 84; // DEL++ DEL++ |
---|
| 411 | else if(PDGC==89998004) return 85; // DEL- DEL- |
---|
| 412 | else if(PDGC==90999002) return 104; // n Sigma- |
---|
| 413 | else if(PDGC==91001999) return 105; // p Sigma+ |
---|
| 414 | } |
---|
| 415 | if(b==3) |
---|
| 416 | { |
---|
| 417 | if (PDGC==90003999) return 86; // p p DEL++ |
---|
| 418 | else if(PDGC==89999004) return 87; // n n DEL- |
---|
| 419 | else if(PDGC==90004998) return 88; // p DEL++ DEL++ |
---|
| 420 | else if(PDGC==89998005) return 89; // n DEL- DEL- |
---|
| 421 | else if(PDGC==90999003) return 113; // n n Sigma- |
---|
| 422 | } |
---|
| 423 | } |
---|
| 424 | } |
---|
| 425 | } |
---|
| 426 | if (PDGC<80000000) // ----> Direct Baryons & Mesons |
---|
| 427 | { |
---|
| 428 | if (PDGC<100) // => Leptons and field bosons |
---|
| 429 | { |
---|
| 430 | if (PDGC==10) return -1; // Chipolino |
---|
| 431 | else if(PDGC==11) return 0; // e- |
---|
| 432 | else if(PDGC==12) return 1; // nu_e |
---|
| 433 | else if(PDGC==13) return 2; // mu- |
---|
| 434 | else if(PDGC==14) return 3; // nu_mu |
---|
| 435 | else if(PDGC==15) return 4; // tau- |
---|
| 436 | else if(PDGC==16) return 5; // nu_tau |
---|
| 437 | else if(PDGC==22) return 6; // Photon |
---|
| 438 | else if(PDGC==23) return 7; // Z0 boson |
---|
| 439 | else if(PDGC==24) return 8; // W- boson |
---|
| 440 | else if(PDGC==25) return 9; // H0 (neutral Higs boson) |
---|
| 441 | else if(PDGC==37) return 10; // H- (charged Higs boson) |
---|
| 442 | } |
---|
| 443 | G4int r=PDGC%10; // 2s+1 |
---|
| 444 | G4int Q= 0; |
---|
| 445 | if (!r) |
---|
| 446 | { |
---|
| 447 | // Internal CHIPS codes for the wide f_0 states must be 9000221, 9010221, 10221 |
---|
| 448 | if (PDGC==110) return 11; // Low R-P: Sigma (pi,pi S-wave) |
---|
| 449 | else if(PDGC==220) return 12; // Midle Regeon-Pomeron |
---|
| 450 | else if(PDGC==330) return 13; // High Regeon-Pomeron |
---|
| 451 | #ifdef pdebug |
---|
| 452 | G4cout<<"***G4QPDGCode::MakeQCode: (0) Unknown in Q-System code: "<<PDGCode<<G4endl; |
---|
| 453 | #endif |
---|
| 454 | return -2; |
---|
| 455 | } |
---|
| 456 | else Q=qr[r]; |
---|
| 457 | G4int p=PDGC/10; // Quark Content |
---|
| 458 | if(r%2) // (2s+1 is odd) Mesons are all the same |
---|
| 459 | { |
---|
| 460 | if (p==11) return Q+=1; |
---|
| 461 | else if(p==21) return Q+=2; |
---|
| 462 | else if(p==22) return Q+=3; |
---|
| 463 | else if(p==31) return Q+=4; |
---|
| 464 | else if(p==32) return Q+=5; |
---|
| 465 | else if(p==33) return Q+=6; |
---|
| 466 | else |
---|
| 467 | { |
---|
| 468 | #ifdef pdebug |
---|
| 469 | G4cout<<"***G4QPDGCode::MakeQCode:(1) Unknown in Q-System code: "<<PDGCode<<G4endl; |
---|
| 470 | #endif |
---|
| 471 | return -2; |
---|
| 472 | } |
---|
| 473 | } |
---|
| 474 | else // (2s+1 is even) Baryons |
---|
| 475 | { |
---|
| 476 | G4int s=r/2; |
---|
| 477 | if(s%2) // ((2s+1)/2 is odd) N Family |
---|
| 478 | { |
---|
| 479 | if (p==211) return Q+=1; |
---|
| 480 | else if(p==221) return Q+=2; |
---|
| 481 | else if(p==312) return Q+=3; |
---|
| 482 | else if(p==311) return Q+=4; |
---|
| 483 | else if(p==321) return Q+=5; |
---|
| 484 | else if(p==322) return Q+=6; |
---|
| 485 | else if(p==331) return Q+=7; |
---|
| 486 | else if(p==332) return Q+=8; |
---|
| 487 | else |
---|
| 488 | { |
---|
| 489 | #ifdef pdebug |
---|
| 490 | G4cout<<"**G4QPDGCode::MakeQCode:(2) Unknown in Q-System code:"<<PDGCode<<G4endl; |
---|
| 491 | #endif |
---|
| 492 | return -2; |
---|
| 493 | } |
---|
| 494 | } |
---|
| 495 | else // ((2s+1)/2 is odd) Delta Family |
---|
| 496 | { |
---|
| 497 | if (p==111) return Q+= 1; |
---|
| 498 | else if(p==211) return Q+= 2; |
---|
| 499 | else if(p==221) return Q+= 3; |
---|
| 500 | else if(p==222) return Q+= 4; |
---|
| 501 | else if(p==312) return Q+= 5; |
---|
| 502 | else if(p==311) return Q+= 6; |
---|
| 503 | else if(p==321) return Q+= 7; |
---|
| 504 | else if(p==322) return Q+= 8; |
---|
| 505 | else if(p==331) return Q+= 9; |
---|
| 506 | else if(p==332) return Q+=10; |
---|
| 507 | else if(p==333) return Q+=11; |
---|
| 508 | else |
---|
| 509 | { |
---|
| 510 | #ifdef pdebug |
---|
| 511 | G4cout<<"**G4QPDGCode::MakeQCode:(3) Unknown in Q-System code:"<<PDGCode<<G4endl; |
---|
| 512 | #endif |
---|
| 513 | return -2; |
---|
| 514 | } |
---|
| 515 | } |
---|
| 516 | } |
---|
| 517 | } |
---|
| 518 | else // Nuclear Fragments |
---|
| 519 | { |
---|
| 520 | G4int d=n+n+z+s; // a#of d quarks |
---|
| 521 | G4int u=n+z+z+s; // a#of u quarks |
---|
| 522 | G4int t=d+u+s; // tot#of quarks |
---|
| 523 | if(t%3 || abs(t)<3) // b=0 are in mesons |
---|
| 524 | { |
---|
| 525 | #ifdef pdebug |
---|
| 526 | G4cout<<"***G4QPDGCode::MakeQCode: Unknown PDGCode="<<PDGCode<<", t="<<t<<G4endl; |
---|
| 527 | #endif |
---|
| 528 | return -2; |
---|
| 529 | } |
---|
| 530 | else |
---|
| 531 | { |
---|
| 532 | G4int b=t/3; // baryon number |
---|
| 533 | if(b==1) // baryons |
---|
| 534 | { |
---|
| 535 | if (s==0&&u==1&&d==2) return 90; |
---|
| 536 | else if(s==0&&u==2&&d==1) return 91; |
---|
| 537 | else if(s==1&&u==1&&d==1) return 92; |
---|
| 538 | else |
---|
| 539 | { |
---|
| 540 | #ifdef pdebug |
---|
| 541 | G4cout<<"**G4QPDGCode::MakeQCode:(5) Unknown in Q-System code:"<<PDGCode<<G4endl; |
---|
| 542 | #endif |
---|
| 543 | return -2; |
---|
| 544 | } |
---|
| 545 | } |
---|
| 546 | else if(b==2) // di-baryons |
---|
| 547 | { |
---|
| 548 | if (s==0&&u==2&&d==4) return 98; |
---|
| 549 | else if(s==0&&u==3&&d==3) return 99; |
---|
| 550 | else if(s==0&&u==4&&d==2) return 100; |
---|
| 551 | else if(s==1&&u==2&&d==3) return 101; |
---|
| 552 | else if(s==1&&u==3&&d==2) return 102; |
---|
| 553 | else if(s==2&&u==2&&d==2) return 103; |
---|
| 554 | else |
---|
| 555 | { |
---|
| 556 | #ifdef pdebug |
---|
| 557 | G4cout<<"**G4QPDGCode::MakeQCode:(6) Unknown in Q-System code:"<<PDGCode<<G4endl; |
---|
| 558 | #endif |
---|
| 559 | return -2; |
---|
| 560 | } |
---|
| 561 | } |
---|
| 562 | else if(b==3) // tri-baryons |
---|
| 563 | { |
---|
| 564 | if (s==0&&u==4&&d==5) return 106; // pnn |
---|
| 565 | else if(s==0&&u==5&&d==4) return 107; // npp |
---|
| 566 | else if(s==1&&u==3&&d==5) return 108; // Lnn |
---|
| 567 | else if(s==1&&u==4&&d==4) return 109; // Lnp |
---|
| 568 | else if(s==1&&u==5&&d==3) return 110; // Lpp |
---|
| 569 | else if(s==2&&u==3&&d==4) return 111; // LLn |
---|
| 570 | else if(s==2&&u==4&&d==3) return 112; // LLp |
---|
| 571 | else if(s==1&&u==2&&d==6) return 113; // SIG-nn |
---|
| 572 | else |
---|
| 573 | { |
---|
| 574 | #ifdef pdebug |
---|
| 575 | G4cout<<"**G4QPDGCode::MakeQCode:(7) Unknown in Q-System code:"<<PDGCode<<G4endl; |
---|
| 576 | #endif |
---|
| 577 | return -2; |
---|
| 578 | } |
---|
| 579 | } |
---|
| 580 | G4int c=b/2; |
---|
| 581 | G4int g=b%2; |
---|
| 582 | G4int h=c*3; |
---|
| 583 | //G4int Q=57+c*15; |
---|
| 584 | //G4int Q=65+c*15; // "IsoNuclei" |
---|
| 585 | G4int Q=83+c*15; // "Leptons/Hyperons" |
---|
| 586 | u-=h; |
---|
| 587 | d-=h; |
---|
| 588 | if(g) |
---|
| 589 | { |
---|
| 590 | if (s==0&&u==1&&d==2) return Q+= 9; |
---|
| 591 | else if(s==0&&u==2&&d==1) return Q+=10; |
---|
| 592 | else if(s==1&&u==0&&d==2) return Q+=11; |
---|
| 593 | else if(s==1&&u==1&&d==1) return Q+=12; |
---|
| 594 | else if(s==1&&u==2&&d==0) return Q+=13; |
---|
| 595 | else if(s==2&&u==0&&d==1) return Q+=14; |
---|
| 596 | else if(s==2&&u==1&&d==0) return Q+=15; |
---|
| 597 | else |
---|
| 598 | { |
---|
| 599 | #ifdef debug |
---|
| 600 | G4cout<<"**G4QPDGCode::MakeQCode:(8) Unknown in Q-System code:"<<PDGCode<<G4endl; |
---|
| 601 | #endif |
---|
| 602 | return -2; |
---|
| 603 | } |
---|
| 604 | } |
---|
| 605 | else |
---|
| 606 | { |
---|
| 607 | if (s==0&&u==-1&&d== 1) return Q+=1; |
---|
| 608 | else if(s==0&&u== 0&&d== 0) return Q+=2; |
---|
| 609 | else if(s==0&&u== 1&&d==-1) return Q+=3; |
---|
| 610 | else if(s==1&&u==-1&&d== 0) return Q+=4; |
---|
| 611 | else if(s==1&&u== 0&&d==-1) return Q+=5; |
---|
| 612 | else if(s==2&&u==-2&&d== 0) return Q+=6; |
---|
| 613 | else if(s==2&&u==-1&&d==-1) return Q+=7; |
---|
| 614 | else if(s==2&&u== 0&&d==-2) return Q+=8; |
---|
| 615 | else |
---|
| 616 | { |
---|
| 617 | #ifdef debug |
---|
| 618 | G4cout<<"**G4QPDGCode::MakeQCode:(9) Unknown in Q-System code:"<<PDGCode<<G4endl; |
---|
| 619 | #endif |
---|
| 620 | return -2; |
---|
| 621 | } |
---|
| 622 | } |
---|
| 623 | } |
---|
| 624 | } |
---|
| 625 | #ifdef pdebug |
---|
| 626 | G4cout<<"***G4QPDGCode::MakeQCode: () Unknown in Q-System code: "<<PDGCode<<G4endl; |
---|
| 627 | #endif |
---|
| 628 | // return -2; not reachable statement |
---|
| 629 | } |
---|
| 630 | |
---|
| 631 | // Get the mean mass value for the PDG |
---|
| 632 | G4double G4QPDGCode::GetMass() |
---|
| 633 | {// ===================== |
---|
| 634 | G4int ab=theQCode; |
---|
| 635 | #ifdef debug |
---|
| 636 | G4cout<<"G4QPDGCode::GetMass: Mass for Q="<<ab<<",PDG="<<thePDGCode<<",N="<<nQHM<<G4endl; |
---|
| 637 | #endif |
---|
| 638 | if ( (ab < 0 && thePDGCode < 80000000) || !thePDGCode) { |
---|
| 639 | #ifdef debug |
---|
| 640 | if(thePDGCode!=10) |
---|
| 641 | G4cout<<"**G4QPDGCode::GetMass:m=100000.,QC="<<theQCode<<",PDG="<<thePDGCode<<G4endl; |
---|
| 642 | #endif |
---|
| 643 | return 100000.; |
---|
| 644 | } |
---|
| 645 | else if(ab>-1 && ab<nQHM) |
---|
| 646 | { |
---|
| 647 | #ifdef debug |
---|
| 648 | G4cout<<"G4QPDGCode::GetMa:m="<<QHaM(ab)<<",Q="<<theQCode<<",PDG="<<thePDGCode<<G4endl; |
---|
| 649 | #endif |
---|
| 650 | return QHaM(ab); // Get mass from the table |
---|
| 651 | } |
---|
| 652 | //if(szn==0) return m[15]; // @@ mPi0 @@ MK ? |
---|
| 653 | if(thePDGCode==90000000) |
---|
| 654 | { |
---|
| 655 | #ifdef debug |
---|
| 656 | G4cout<<"G4QPDGCode::GetMass:***m=0, QC="<<theQCode<<",PDG="<<thePDGCode<<G4endl; |
---|
| 657 | #endif |
---|
| 658 | return 0.; |
---|
| 659 | } |
---|
| 660 | G4int s=0; |
---|
| 661 | G4int z=0; |
---|
| 662 | G4int n=0; |
---|
| 663 | ConvertPDGToZNS(thePDGCode, z, n, s); |
---|
| 664 | G4double m=GetNuclMass(z,n,s); |
---|
| 665 | #ifdef debug |
---|
| 666 | G4cout<<"G4QPDG::GetM:PDG="<<thePDGCode<<"=>Z="<<z<<",N="<<n<<",S="<<s<<",M="<<m<<G4endl; |
---|
| 667 | #endif |
---|
| 668 | return m; |
---|
| 669 | } |
---|
| 670 | |
---|
| 671 | // Get the width value for the PDG |
---|
| 672 | G4double G4QPDGCode::GetWidth() |
---|
| 673 | // ===================== |
---|
| 674 | { |
---|
| 675 | //static const int nW = 72; |
---|
| 676 | //static const int nW = 80; // "Isobars" |
---|
| 677 | static const int nW = 90; // "Leptons/Hypernuclei" |
---|
| 678 | static G4double width[nW] = {0.,0.,0.,0.,0.,0.,0.,2.495,2.118,10. |
---|
| 679 | , 10., 800., 75., 350., 0., 0., .00118, 0., 0., .203 |
---|
| 680 | , 0., 0., 0., 0., 0., 0., 0., 0., 160., 160. |
---|
| 681 | , 8.41, 50.5, 50.8, 4.43, 120., 120., 120., 120., 15.6, 39. |
---|
| 682 | , 36., 35.8, 10., 9., 0., 107., 107., 185.5, 109., 98.5 |
---|
| 683 | , 76., 130., 130., 80., 120., 120., 120., 20., 20., 160. |
---|
| 684 | , 160., 168., 159., 159., 87., 300., 300., 300., 300., 200. |
---|
| 685 | , 180., 180., 180., 99., 99., 55., 387., 387., 208., 198. |
---|
| 686 | , 198., 149., 120., 120., 170., 170., 120., 120., 170., 170.}; |
---|
| 687 | G4int ab=abs(theQCode); |
---|
| 688 | if(ab<nW) return width[ab]; |
---|
| 689 | return 0.; // @@ May be real width should be implemented for nuclei e.g. pp |
---|
| 690 | } |
---|
| 691 | |
---|
| 692 | // Get a nuclear mass for Z (a#of protons), N (a#of neutrons), & S (a#of lambdas) |
---|
| 693 | G4double G4QPDGCode::GetNuclMass(G4int z, G4int n, G4int s) |
---|
| 694 | // ==================================================== |
---|
| 695 | { |
---|
| 696 | static const G4double anb = .01; // Antibinding for Ksi-n,Sig-n,Sig+p,Sig-nn, |
---|
| 697 | static const G4double mNeut= QHaM(20); |
---|
| 698 | static const G4double mProt= QHaM(21); |
---|
| 699 | static const G4double mLamb= QHaM(22); |
---|
| 700 | static const G4double mPiC = QHaM(15); |
---|
| 701 | static const G4double mKZ = QHaM(17); |
---|
| 702 | static const G4double mKM = QHaM(18); |
---|
| 703 | static const G4double mSiM = QHaM(23); |
---|
| 704 | static const G4double mSiP = QHaM(25); |
---|
| 705 | static const G4double mKsZ = QHaM(27); |
---|
| 706 | static const G4double mKsM = QHaM(26); |
---|
| 707 | static const G4double mOmM = QHaM(44); |
---|
| 708 | static const G4double mKZa = mKZ +anb; |
---|
| 709 | static const G4double mKMa = mKM +anb; |
---|
| 710 | static const G4double mSigM= mSiM+anb; |
---|
| 711 | static const G4double mSigP= mSiP+anb; |
---|
| 712 | static const G4double mKsiZ= mKsZ+anb; |
---|
| 713 | static const G4double mKsiM= mKsM+anb; |
---|
| 714 | static const G4double mOmeg= mOmM+anb; |
---|
| 715 | static const G4double mDiPi= mPiC+mPiC+anb; |
---|
| 716 | static const G4double mDiKZ= mKZa+mKZ; |
---|
| 717 | static const G4double mDiKM= mKMa+mKM; |
---|
| 718 | static const G4double mDiPr= mProt+mProt; |
---|
| 719 | static const G4double mDiNt= mNeut+mNeut; |
---|
| 720 | static const G4double mSmPi= mSiM+mDiPi; |
---|
| 721 | static const G4double mSpPi= mSiP+mDiPi; |
---|
| 722 | static const G4double mOmN = mOmeg+mNeut; |
---|
| 723 | static const G4double mSpP = mSigP+mProt; |
---|
| 724 | static const G4double mSpPP= mSpP +mProt; |
---|
| 725 | static const G4double mSmN = mSigM+mNeut; |
---|
| 726 | static const G4double mSmNN= mSmN +mNeut; |
---|
| 727 | // -------------- DAM Arrays ---------------------- |
---|
| 728 | static const G4int iNR=76; // Neutron maximum range for each Z |
---|
| 729 | static const G4int nEl = 105; // Maximum Z of the associative memory is "nEl-1=104" |
---|
| 730 | static const G4int iNF[nEl]={0,0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 14 |
---|
| 731 | 1 , 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, // 29 |
---|
| 732 | 16 , 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, // 44 |
---|
| 733 | 31 , 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 53, 54, 55, // 59 |
---|
| 734 | 56 , 56, 57, 57, 58, 60, 61, 63, 64, 65, 66, 67, 68, 69, 70, // 74 |
---|
| 735 | 71 , 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, // 89 |
---|
| 736 | 86 , 87, 88, 89, 91, 94, 98,103,109,115,122,128,134,140,146};//104 |
---|
| 737 | #ifdef qdebug |
---|
| 738 | static G4int iNmin[nEl]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 14 |
---|
| 739 | 1 , 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, // 29 |
---|
| 740 | 16 , 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, // 44 |
---|
| 741 | 31 , 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 53, 54, 55, // 59 |
---|
| 742 | 56 , 56, 57, 57, 58, 60, 61, 63, 64, 65, 66, 67, 68, 69, 70, // 74 |
---|
| 743 | 71 , 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, // 89 |
---|
| 744 | 86 , 87, 88, 89, 91, 94, 98,103,109,115,122,128,134,140,146};//104 |
---|
| 745 | static G4int iNmax=iNR; |
---|
| 746 | static G4int iNran[nEl]={19,20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, // 14 |
---|
| 747 | 34 , 35, 36, 37, 38, 39, 40, 48, 48, 48, 48, 50, 50, 50, 52, // 29 |
---|
| 748 | 53 , 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, // 44 |
---|
| 749 | 68 , 69, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, // 59 |
---|
| 750 | 73 , 73, 73, 73, 74, 74, 74, 74, 74, 74, 74, 74, 74, 75, 76, // 74 |
---|
| 751 | 76 , 76, 76, 76, 76, 75, 74, 73, 72, 71, 70, 70, 69, 69, 69, // 89 |
---|
| 752 | 68 , 68, 68, 67, 63, 59, 55, 51, 47, 43, 39, 35, 31, 27, 23};//104 |
---|
| 753 | #endif |
---|
| 754 | static const G4int iNL[nEl]={19,20,21,22,23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, // 14 |
---|
| 755 | 34 , 35, 36, 37, 38, 39, 40, 48, 48, 48, 48, 50, 50, 50, 52, // 29 |
---|
| 756 | 53 , 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, // 44 |
---|
| 757 | 68 , 69, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, // 59 |
---|
| 758 | 73 , 73, 73, 73, 74, 74, 74, 74, 74, 74, 74, 74, 74, 75, 76, // 74 |
---|
| 759 | 76 , 76, 76, 76, 76, 75, 74, 73, 72, 71, 70, 70, 69, 69, 69, // 89 |
---|
| 760 | 68 , 68, 68, 67, 63, 59, 55, 51, 47, 43, 39, 35, 31, 27, 23};//104 |
---|
| 761 | // ********* S=-4 vectors ************* |
---|
| 762 | static G4bool iNin6[nEl]={false,false,false,false,false,false,false, |
---|
| 763 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 764 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 765 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 766 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 767 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 768 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 769 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 770 | static G4double VZ6[nEl][iNR]; |
---|
| 771 | //********* S=-3 vectors ************* |
---|
| 772 | static G4bool iNin7[nEl]={false,false,false,false,false,false,false, |
---|
| 773 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 774 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 775 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 776 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 777 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 778 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 779 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 780 | static G4double VZ7[nEl][iNR]; |
---|
| 781 | // ********* S=-2 vectors ************* |
---|
| 782 | static G4bool iNin8[nEl]={false,false,false,false,false,false,false, |
---|
| 783 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 784 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 785 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 786 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 787 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 788 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 789 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 790 | static G4double VZ8[nEl][iNR]; |
---|
| 791 | // ********* S=-1 vectors ************* |
---|
| 792 | static G4bool iNin9[nEl]={false,false,false,false,false,false,false, |
---|
| 793 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 794 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 795 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 796 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 797 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 798 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 799 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 800 | static G4double VZ9[nEl][iNR]; |
---|
| 801 | // ********* S=0 vectors ************* |
---|
| 802 | static G4bool iNin0[nEl]={false,false,false,false,false,false,false, |
---|
| 803 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 804 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 805 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 806 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 807 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 808 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 809 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 810 | static G4double VZ0[nEl][iNR]; |
---|
| 811 | // ********* S=1 vectors ************* |
---|
| 812 | static G4bool iNin1[nEl]={false,false,false,false,false,false,false, |
---|
| 813 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 814 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 815 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 816 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 817 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 818 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 819 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 820 | static G4double VZ1[nEl][iNR]; |
---|
| 821 | // ********* S=2 vectors ************* |
---|
| 822 | static G4bool iNin2[nEl]={false,false,false,false,false,false,false, |
---|
| 823 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 824 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 825 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 826 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 827 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 828 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 829 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 830 | static G4double VZ2[nEl][iNR]; |
---|
| 831 | // ********* S=3 vectors ************* |
---|
| 832 | static G4bool iNin3[nEl]={false,false,false,false,false,false,false, |
---|
| 833 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 834 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 835 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 836 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 837 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 838 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 839 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 840 | static G4double VZ3[nEl][iNR]; |
---|
| 841 | // ********* S=2 vectors ************* |
---|
| 842 | static G4bool iNin4[nEl]={false,false,false,false,false,false,false, |
---|
| 843 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 844 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 845 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 846 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 847 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 848 | false,false,false,false,false,false,false,false,false,false,false,false,false,false, |
---|
| 849 | false,false,false,false,false,false,false,false,false,false,false,false,false,false}; |
---|
| 850 | static G4double VZ4[nEl][iNR]; |
---|
| 851 | // |
---|
| 852 | #ifdef qdebug |
---|
| 853 | static G4int Smin=-1; // Dynamic Associated memory is implemented for S>-2 nuclei |
---|
| 854 | static G4int Smax= 2; // Dynamic Associated memory is implemented for S< 3 nuclei |
---|
| 855 | static G4int NZmin= 0; // Dynamic Associated memory is implemented for Z>-1 nuclei |
---|
| 856 | static G4int NNmin= 0; // Dynamic Associated memory is implemented for N>-1 nuclei |
---|
| 857 | static G4int NZS1max= 0; // Dynamic Associated memory is implemented for S<3, Z=-1, N<2 |
---|
| 858 | static G4int NNS1max= 0; // Dynamic Associated memory is implemented for S<3, Z=-1, N<2 |
---|
| 859 | #endif |
---|
| 860 | // ------------------------------------------------------------------------------------- |
---|
| 861 | G4double rm=0.; |
---|
| 862 | G4int nz=n+z; |
---|
| 863 | G4int zns=nz+s; |
---|
| 864 | if(nz+s<0) |
---|
| 865 | { |
---|
| 866 | z=-z; |
---|
| 867 | n=-n; |
---|
| 868 | s=-s; |
---|
| 869 | nz=-nz; |
---|
| 870 | zns=-zns; |
---|
| 871 | } |
---|
| 872 | if(z<0) |
---|
| 873 | { |
---|
| 874 | if(z==-1) |
---|
| 875 | { |
---|
| 876 | if(!s) |
---|
| 877 | { |
---|
| 878 | if(n==1) return mPiC; // pi- |
---|
| 879 | else return mPiC+(n-1)*mNeut; // Delta- + (N-1)*n |
---|
| 880 | } |
---|
| 881 | else if(s==1) // Strange negative hadron |
---|
| 882 | { |
---|
| 883 | if(!n) return mKM; // K- |
---|
| 884 | else if(n==1) return mSiM; // Sigma- |
---|
| 885 | else if(n==2) return mSmN ; // Sigma- + n DiBaryon |
---|
| 886 | else if(n==3) return mSmNN; // Sigma- +2n TriBaryon |
---|
| 887 | else return mSigM+mNeut*(n-1); // Sigma- + (N-1)*n |
---|
| 888 | } |
---|
| 889 | else if(s==2) // --> Double-strange negative hadrons |
---|
| 890 | { |
---|
| 891 | if(!n) return mKsM; // Ksi- |
---|
| 892 | else if(n==1) return mKsiM+mNeut; // Ksi- + n |
---|
| 893 | else if(n==2) return mKsiM+mNeut+mNeut; // Ksi- + 2n |
---|
| 894 | else return mKsiM+mNeut*n; // Ksi- + Z*n |
---|
| 895 | } |
---|
| 896 | else if(s==-2) |
---|
| 897 | { |
---|
| 898 | if (nz==2) return mDiKZ+mPiC; // 2K0 + Pi- |
---|
| 899 | else return mDiKZ+mPiC+(nz-2)*mProt; |
---|
| 900 | } |
---|
| 901 | else if(s==3) // --> Triple-strange negative hadrons |
---|
| 902 | { |
---|
| 903 | if (n==-1) return mOmM; // Triple-strange Omega minus |
---|
| 904 | else if(!n ) return mOmN; // Triple-strange (Omega-) + n DiBaryon |
---|
| 905 | else if(n==-2) return mDiKZ+mKM; // Triple-strange K- + 2*K0 |
---|
| 906 | else return mOmeg+mNeut*(n+2); |
---|
| 907 | } |
---|
| 908 | else if(s==4) |
---|
| 909 | { |
---|
| 910 | if(n==-2) return mOmeg+mKM; // Omega- + K- |
---|
| 911 | else if(n==-1) return mOmeg+mLamb;// Omega- + Lambda |
---|
| 912 | else return mOmeg+mLamb+(n+1)*mNeut; // Omega- + Lambda + (n+1)*Neutrons |
---|
| 913 | } |
---|
| 914 | else if(!n) return mOmeg+(s-2)*mLamb; // Multy-Lambda + Omega minus |
---|
| 915 | else |
---|
| 916 | { |
---|
| 917 | #ifdef qdebug |
---|
| 918 | if(s>NZS1max) |
---|
| 919 | { |
---|
| 920 | NZS1max=s; |
---|
| 921 | G4cout<<">>>>>>>>>>>>>G4QPDGCode::GetMass: Z=-1, S="<<s<<">2 with N="<<n<<G4endl; |
---|
| 922 | } |
---|
| 923 | #endif |
---|
| 924 | return CalculateNuclMass(z,n,s); |
---|
| 925 | } |
---|
| 926 | } |
---|
| 927 | else if(!s) |
---|
| 928 | { |
---|
| 929 | if(!nz) |
---|
| 930 | { |
---|
| 931 | if(n==2) return mDiPi; |
---|
| 932 | else return mDiPi+(n-2)*mPiC; |
---|
| 933 | } |
---|
| 934 | else return mNeut*nz-z*mPiC+anb; |
---|
| 935 | } |
---|
| 936 | else if(!zns) // !!! s=0 is treated above !! |
---|
| 937 | { |
---|
| 938 | if(s>0) return anb+s*mKM+n*mPiC; // s*K- + n*Pi- |
---|
| 939 | else return anb-s*mKZ-z*mPiC; // (-s)*aK0 + (-z)*Pi- |
---|
| 940 | } |
---|
| 941 | else if(s==1) |
---|
| 942 | { |
---|
| 943 | if(!nz) |
---|
| 944 | { |
---|
| 945 | if(n==2) return mSmPi; |
---|
| 946 | else return mSmPi+(n-2)*mPiC; |
---|
| 947 | } |
---|
| 948 | else return mSigM+nz*mNeut-(z+1)*mPiC; |
---|
| 949 | } |
---|
| 950 | else if(s==-1) return mKZa-z*mPiC+(nz-1)*mNeut; // aK0+(nz-1)n+(-z)*Pi- |
---|
| 951 | else if(s==2) |
---|
| 952 | { |
---|
| 953 | if (nz==-1) return mKsiM+n*mPiC; |
---|
| 954 | else if(!nz) return mKsiM+mNeut-(z+1)*mPiC; |
---|
| 955 | else return mKsiM+(nz+1)*mNeut-(z+1)*mPiC; |
---|
| 956 | } |
---|
| 957 | else if(s==-2) return mDiKZ-z*mPiC+(nz-2)*mNeut; |
---|
| 958 | else if(s==3) |
---|
| 959 | { |
---|
| 960 | if (nz==-2) return mOmeg+(n+1)*mPiC; // Omega- + (n+1)* Pi- |
---|
| 961 | else if(nz==-1) return mOmeg+mNeut+n*mPiC; // Omega- + n * Pi- |
---|
| 962 | else if(!nz) return mOmeg+mDiNt+(n-1)*mPiC; // Omega- + 2N + (n-1)*Pi- |
---|
| 963 | else return mOmeg+(nz+2)*mProt-(z+1)*mPiC; |
---|
| 964 | } |
---|
| 965 | else if(s<-2) return anb-s*mKZ-z*mPiC+(nz+s)*mNeut; |
---|
| 966 | else if(s==4) |
---|
| 967 | { |
---|
| 968 | if (nz==-3) return mOmeg+mKM+(n+1)*mPiC; // Om- + K- + (n+1)*Pi- |
---|
| 969 | else if(nz==-2) return mOmeg+mSigM+n*mPiC; // Om- + Sig- + n*Pi- |
---|
| 970 | else if(nz==-1) return mOmeg+mSigM+mNeut+(n-1)*mPiC;//Om- + Sig- +N +(n-1)*Pi- |
---|
| 971 | else if(!nz) return mOmeg+mSigM+mDiNt+(n-2)*mPiC;//Om- +Sig- +2N +(n-2)*Pi- |
---|
| 972 | else return mOmeg+mSigM+(nz+2)*mDiNt-(z+2)*mPiC;//+(nz+2)N-(z+2)Pi- |
---|
| 973 | } |
---|
| 974 | // s=5: 2*K-, Ksi-; s=6: 3*K-, Omega-; s>6 adds K- and Sigma- instead of Protons |
---|
| 975 | else // !!All s<0 are done and s>4 can be easy extended if appear!! |
---|
| 976 | { |
---|
| 977 | #ifdef qdebug |
---|
| 978 | if(z<NZmin) |
---|
| 979 | { |
---|
| 980 | NZmin=z; |
---|
| 981 | G4cout<<">>>>>>>>>G4QPDGCode::GetMass: Z="<<z<<"<-1 with N="<<n<<", S="<<s<<G4endl; |
---|
| 982 | } |
---|
| 983 | #endif |
---|
| 984 | return CalculateNuclMass(z,n,s); |
---|
| 985 | } |
---|
| 986 | } |
---|
| 987 | else if(n<0) |
---|
| 988 | { |
---|
| 989 | if(n==-1) |
---|
| 990 | { |
---|
| 991 | if(!s) |
---|
| 992 | { |
---|
| 993 | if(z==1) return mPiC; // pi+ |
---|
| 994 | else return mPiC+(z-1)*mProt; // Delta++ + (Z-1)*p |
---|
| 995 | } |
---|
| 996 | else if(s==1) // --> Strange neutral hadrons |
---|
| 997 | { |
---|
| 998 | if(!z) return mKZ; // K0 |
---|
| 999 | else if(z==1) return mSiP; // Sigma+ |
---|
| 1000 | else if(z==2) return mSpP ; // Sigma+ + p DiBaryon |
---|
| 1001 | else if(z==3) return mSpPP; // Sigma+ +2p TriBaryon |
---|
| 1002 | else return mSigP+mProt*(z-1); // Sigma+ + (Z-1)*p |
---|
| 1003 | } |
---|
| 1004 | else if(s==2) // --> Double-strange negative hadrons |
---|
| 1005 | { |
---|
| 1006 | if(!z) return mKsZ; // Ksi0 |
---|
| 1007 | else if(z==1) return mKsiZ+mProt; // Ksi- + p |
---|
| 1008 | else if(z==2) return mKsiZ+mProt+mProt; // Ksi- + 2p |
---|
| 1009 | else return mKsiZ+mProt*z; // Ksi- + Z*p |
---|
| 1010 | } |
---|
| 1011 | else if(s==-2) |
---|
| 1012 | { |
---|
| 1013 | if (nz==2) return mDiKM+mPiC; // 2K+ + Pi+ |
---|
| 1014 | else return mDiKM+mPiC+(nz-2)*mProt; |
---|
| 1015 | } |
---|
| 1016 | else if(s==3) |
---|
| 1017 | { |
---|
| 1018 | if(z==1) return mOmeg+mDiPr; |
---|
| 1019 | else return mOmeg+(z+1)*mProt; |
---|
| 1020 | } |
---|
| 1021 | else if(s==4) return mOmeg+mLamb+(z+1)*mProt; |
---|
| 1022 | else if(!z) return mKZa+(s-1)*mLamb; // Multy-Lambda + K0 |
---|
| 1023 | else |
---|
| 1024 | { |
---|
| 1025 | #ifdef qdebug |
---|
| 1026 | if(s>NNS1max) |
---|
| 1027 | { |
---|
| 1028 | NNS1max=s; |
---|
| 1029 | G4cout<<">>>>>>>>>>>>>G4QPDGCode::GetMass: N=-1, S="<<s<<">2 with Z="<<z<<G4endl; |
---|
| 1030 | } |
---|
| 1031 | #endif |
---|
| 1032 | return CalculateNuclMass(z,n,s); |
---|
| 1033 | } |
---|
| 1034 | } |
---|
| 1035 | else if(!s) |
---|
| 1036 | { |
---|
| 1037 | if(!nz) |
---|
| 1038 | { |
---|
| 1039 | if(z==2) return mDiPi; |
---|
| 1040 | else return mDiPi+(z-2)*mPiC; |
---|
| 1041 | } |
---|
| 1042 | else return mProt*nz-n*mPiC+anb; |
---|
| 1043 | } |
---|
| 1044 | else if(!zns) // !!! s=0 is treated above !! |
---|
| 1045 | { |
---|
| 1046 | if(s>0) return anb+s*mKZ+z*mPiC; // s*K0 + n*Pi+ |
---|
| 1047 | else return anb-s*mKM-n*mPiC; // (-s)*aK+ + (-n)*Pi+ |
---|
| 1048 | } |
---|
| 1049 | else if(s==1) |
---|
| 1050 | { |
---|
| 1051 | if(!nz) |
---|
| 1052 | { |
---|
| 1053 | if(z==2) return mSpPi; |
---|
| 1054 | else return mSpPi+(z-2)*mPiC; |
---|
| 1055 | } |
---|
| 1056 | else return mSigP+nz*mProt-(n+1)*mPiC; |
---|
| 1057 | } |
---|
| 1058 | else if(s==-1) return mKMa-n*mPiC+(nz-1)*mProt; // K+ + (nz-1)*P + (-n)*Pi+ |
---|
| 1059 | else if(s==2) |
---|
| 1060 | { |
---|
| 1061 | if (nz==-1) return mKsiZ+z*mPiC; |
---|
| 1062 | else if(!nz) return mKsiZ+mProt-(n+1)*mPiC; |
---|
| 1063 | else return mKsiZ+(nz+1)*mProt-(n+1)*mPiC; |
---|
| 1064 | } |
---|
| 1065 | else if(s==-2) return mDiKM-n*mPiC+(nz-2)*mProt; |
---|
| 1066 | else if(s==3) |
---|
| 1067 | { |
---|
| 1068 | if (nz==-2) return mOmeg+(z+1)*mPiC; // Omega + (z+1)*Pi+ |
---|
| 1069 | else if(nz==-1) return mOmeg+mProt+z*mPiC; // Omega- + P +z*Pi+ |
---|
| 1070 | else if(!nz) return mOmeg+mDiPr+(z-1)*mPiC; // Omega- + 2P + (z-1)* Pi+ |
---|
| 1071 | else return mOmeg+(nz+2)*mProt-(n+1)*mPiC; |
---|
| 1072 | } |
---|
| 1073 | else if(s<-2) return anb-s*mKM-n*mPiC+(nz+s)*mProt; |
---|
| 1074 | else if(s==4) |
---|
| 1075 | { |
---|
| 1076 | if (nz==-3) return mOmeg+mKZ+(z+1)*mPiC; // Om- + K0 + (z+1)*Pi+ |
---|
| 1077 | else if(nz==-2) return mOmeg+mSigP+z*mPiC; // Om- + Sig+ + z*Pi+ |
---|
| 1078 | else if(nz==-1) return mOmeg+mSigP+mProt+(z-1)*mPiC;// Om- +Sig+ +P +(z-1)*Pi+ |
---|
| 1079 | else if(!nz) return mOmeg+mSigP+mDiPr+(z-2)*mPiC;//Om- +Sig+ +2P +(z-2)*Pi+ |
---|
| 1080 | else return mOmeg+mSigP+(nz+2)*mProt-(n+2)*mPiC;//+(nz+2)P-(n+2)Pi+ |
---|
| 1081 | } |
---|
| 1082 | // s=5: 2*KZ, Ksi0; s=6: 3*KZ, Omega-; s>6 adds K0 and Sigma+ instead of Protons |
---|
| 1083 | else // !!All s<0 are done and s>4 can be easy extended if appear!! |
---|
| 1084 | { |
---|
| 1085 | #ifdef qdebug |
---|
| 1086 | if(n<NNmin) |
---|
| 1087 | { |
---|
| 1088 | NNmin=n; |
---|
| 1089 | G4cout<<">>>>>>>>>G4QPDGCode::GetMass: N="<<n<<"<-1 with Z="<<z<<", S="<<s<<G4endl; |
---|
| 1090 | } |
---|
| 1091 | #endif |
---|
| 1092 | return CalculateNuclMass(z,n,s); |
---|
| 1093 | } |
---|
| 1094 | } |
---|
| 1095 | //return CalculateNuclMass(z,n,s); // @@ This is just to compare the calculation time @@ |
---|
| 1096 | if(!s) // **************> S=0 nucleus |
---|
| 1097 | { |
---|
| 1098 | if(nz==256) return 256000.; |
---|
| 1099 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1100 | if(!iNin0[z]) // ====> This element is already initialized |
---|
| 1101 | { |
---|
| 1102 | #ifdef idebug |
---|
| 1103 | G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=0 is initialized. F="<<iNin0[z]<<G4endl; |
---|
| 1104 | #endif |
---|
| 1105 | G4int iNfin=iNL[z]; |
---|
| 1106 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1107 | for (G4int in=0; in<iNfin; in++) VZ0[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1108 | iNin0[z]=true; |
---|
| 1109 | } |
---|
| 1110 | G4int dNn=n-Nmin; |
---|
| 1111 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1112 | { |
---|
| 1113 | #ifdef qdebug |
---|
| 1114 | if(n<iNmin[z]) |
---|
| 1115 | { |
---|
| 1116 | G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=0 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1117 | iNmin[z]=n; |
---|
| 1118 | } |
---|
| 1119 | #endif |
---|
| 1120 | return CalculateNuclMass(z,n,s); |
---|
| 1121 | } |
---|
| 1122 | else if(dNn<iNL[z]) return VZ0[z][dNn]; // Found in DAM |
---|
| 1123 | else // --> The maximum N must be increased |
---|
| 1124 | { |
---|
| 1125 | #ifdef qdebug |
---|
| 1126 | if(dNn>iNmax) |
---|
| 1127 | { |
---|
| 1128 | G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1129 | iNmax=dNn; |
---|
| 1130 | } |
---|
| 1131 | if(dNn>iNran[z]) |
---|
| 1132 | { |
---|
| 1133 | G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1134 | iNran[z]=dNn; |
---|
| 1135 | } |
---|
| 1136 | #endif |
---|
| 1137 | return CalculateNuclMass(z,n,s); |
---|
| 1138 | } |
---|
| 1139 | } |
---|
| 1140 | else if(s==1) // ******************> S=1 nucleus |
---|
| 1141 | { |
---|
| 1142 | |
---|
| 1143 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1144 | if(!iNin1[z]) // ====> This element is already initialized |
---|
| 1145 | { |
---|
| 1146 | #ifdef idebug |
---|
| 1147 | G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=1 is initialized. F="<<iNin1[z]<<G4endl; |
---|
| 1148 | #endif |
---|
| 1149 | G4int iNfin=iNL[z]; |
---|
| 1150 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1151 | for (G4int in=0; in<iNfin; in++) VZ1[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1152 | iNin1[z]=true; |
---|
| 1153 | } |
---|
| 1154 | G4int dNn=n-Nmin; |
---|
| 1155 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1156 | { |
---|
| 1157 | #ifdef qdebug |
---|
| 1158 | if(n<iNmin[z]) |
---|
| 1159 | { |
---|
| 1160 | G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=1 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1161 | iNmin[z]=n; |
---|
| 1162 | } |
---|
| 1163 | #endif |
---|
| 1164 | return CalculateNuclMass(z,n,s); |
---|
| 1165 | } |
---|
| 1166 | else if(dNn<iNL[z]) return VZ1[z][dNn]; // Found in DAM |
---|
| 1167 | else // --> The maximum N must be increased |
---|
| 1168 | { |
---|
| 1169 | #ifdef qdebug |
---|
| 1170 | if(dNn>iNmax) |
---|
| 1171 | { |
---|
| 1172 | G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1173 | iNmax=dNn; |
---|
| 1174 | } |
---|
| 1175 | if(dNn>iNran[z]) |
---|
| 1176 | { |
---|
| 1177 | G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1178 | iNran[z]=dNn; |
---|
| 1179 | } |
---|
| 1180 | #endif |
---|
| 1181 | return CalculateNuclMass(z,n,s); |
---|
| 1182 | } |
---|
| 1183 | } |
---|
| 1184 | else if(s==-1) // ******************> S=-1 nucleus |
---|
| 1185 | { |
---|
| 1186 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1187 | if(!iNin9[z]) // ====> This element is already initialized |
---|
| 1188 | { |
---|
| 1189 | #ifdef idebug |
---|
| 1190 | G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-1 is initialized. F="<<iNin9[z]<<G4endl; |
---|
| 1191 | #endif |
---|
| 1192 | G4int iNfin=iNL[z]; |
---|
| 1193 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1194 | for (G4int in=0; in<iNfin; in++) VZ9[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1195 | iNin9[z]=true; |
---|
| 1196 | } |
---|
| 1197 | G4int dNn=n-Nmin; |
---|
| 1198 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1199 | { |
---|
| 1200 | #ifdef qdebug |
---|
| 1201 | if(n<iNmin[z]) |
---|
| 1202 | { |
---|
| 1203 | G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<" ,S=-1 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1204 | iNmin[z]=n; |
---|
| 1205 | } |
---|
| 1206 | #endif |
---|
| 1207 | return CalculateNuclMass(z,n,s); |
---|
| 1208 | } |
---|
| 1209 | else if(dNn<iNL[z]) return VZ9[z][dNn]; // Found in DAM |
---|
| 1210 | else // --> The maximum N must be increased |
---|
| 1211 | { |
---|
| 1212 | #ifdef qdebug |
---|
| 1213 | if(dNn>iNmax) |
---|
| 1214 | { |
---|
| 1215 | G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1216 | iNmax=dNn; |
---|
| 1217 | } |
---|
| 1218 | if(dNn>iNran[z]) |
---|
| 1219 | { |
---|
| 1220 | G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1221 | iNran[z]=dNn; |
---|
| 1222 | } |
---|
| 1223 | #endif |
---|
| 1224 | return CalculateNuclMass(z,n,s); |
---|
| 1225 | } |
---|
| 1226 | } |
---|
| 1227 | else if(s==2) // ******************> S=2 nucleus |
---|
| 1228 | { |
---|
| 1229 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1230 | if(!iNin2[z]) // ====> This element is already initialized |
---|
| 1231 | { |
---|
| 1232 | #ifdef idebug |
---|
| 1233 | G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=2 is initialized. F="<<iNin2[z]<<G4endl; |
---|
| 1234 | #endif |
---|
| 1235 | G4int iNfin=iNL[z]; |
---|
| 1236 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1237 | for (G4int in=0; in<iNfin; in++) VZ2[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1238 | iNin2[z]=true; |
---|
| 1239 | } |
---|
| 1240 | G4int dNn=n-Nmin; |
---|
| 1241 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1242 | { |
---|
| 1243 | #ifdef qdebug |
---|
| 1244 | if(n<iNmin[z]) |
---|
| 1245 | { |
---|
| 1246 | G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=2 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1247 | iNmin[z]=n; |
---|
| 1248 | } |
---|
| 1249 | #endif |
---|
| 1250 | return CalculateNuclMass(z,n,s); |
---|
| 1251 | } |
---|
| 1252 | else if(dNn<iNL[z]) return VZ2[z][dNn]; // Found in DAM |
---|
| 1253 | else // --> The maximum N must be increased |
---|
| 1254 | { |
---|
| 1255 | #ifdef qdebug |
---|
| 1256 | if(dNn>iNmax) |
---|
| 1257 | { |
---|
| 1258 | G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1259 | iNmax=dNn; |
---|
| 1260 | } |
---|
| 1261 | if(dNn>iNran[z]) |
---|
| 1262 | { |
---|
| 1263 | G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1264 | iNran[z]=dNn; |
---|
| 1265 | } |
---|
| 1266 | #endif |
---|
| 1267 | return CalculateNuclMass(z,n,s); |
---|
| 1268 | } |
---|
| 1269 | } |
---|
| 1270 | else if(s==-2) // ******************> S=-2 nucleus |
---|
| 1271 | { |
---|
| 1272 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1273 | if(!iNin8[z]) // ====> This element is already initialized |
---|
| 1274 | { |
---|
| 1275 | #ifdef idebug |
---|
| 1276 | G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-2 is initialized. F="<<iNin8[z]<<G4endl; |
---|
| 1277 | #endif |
---|
| 1278 | G4int iNfin=iNL[z]; |
---|
| 1279 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1280 | for (G4int in=0; in<iNfin; in++) VZ8[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1281 | iNin8[z]=true; |
---|
| 1282 | } |
---|
| 1283 | G4int dNn=n-Nmin; |
---|
| 1284 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1285 | { |
---|
| 1286 | #ifdef qdebug |
---|
| 1287 | if(n<iNmin[z]) |
---|
| 1288 | { |
---|
| 1289 | G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-2 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1290 | iNmin[z]=n; |
---|
| 1291 | } |
---|
| 1292 | #endif |
---|
| 1293 | return CalculateNuclMass(z,n,s); |
---|
| 1294 | } |
---|
| 1295 | else if(dNn<iNL[z]) return VZ8[z][dNn]; // Found in DAM |
---|
| 1296 | else // --> The maximum N must be increased |
---|
| 1297 | { |
---|
| 1298 | #ifdef qdebug |
---|
| 1299 | if(dNn>iNmax) |
---|
| 1300 | { |
---|
| 1301 | G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1302 | iNmax=dNn; |
---|
| 1303 | } |
---|
| 1304 | if(dNn>iNran[z]) |
---|
| 1305 | { |
---|
| 1306 | G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1307 | iNran[z]=dNn; |
---|
| 1308 | } |
---|
| 1309 | #endif |
---|
| 1310 | return CalculateNuclMass(z,n,s); |
---|
| 1311 | } |
---|
| 1312 | } |
---|
| 1313 | else if(s==-3) // ******************> S=-3 nucleus |
---|
| 1314 | { |
---|
| 1315 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1316 | if(!iNin7[z]) // ====> This element is already initialized |
---|
| 1317 | { |
---|
| 1318 | #ifdef idebug |
---|
| 1319 | G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-3 is initialized. F="<<iNin7[z]<<G4endl; |
---|
| 1320 | #endif |
---|
| 1321 | G4int iNfin=iNL[z]; |
---|
| 1322 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1323 | for (G4int in=0; in<iNfin; in++) VZ7[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1324 | iNin7[z]=true; |
---|
| 1325 | } |
---|
| 1326 | G4int dNn=n-Nmin; |
---|
| 1327 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1328 | { |
---|
| 1329 | #ifdef qdebug |
---|
| 1330 | if(n<iNmin[z]) |
---|
| 1331 | { |
---|
| 1332 | G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-3 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1333 | iNmin[z]=n; |
---|
| 1334 | } |
---|
| 1335 | #endif |
---|
| 1336 | return CalculateNuclMass(z,n,s); |
---|
| 1337 | } |
---|
| 1338 | else if(dNn<iNL[z]) return VZ7[z][dNn]; // Found in DAM |
---|
| 1339 | else // --> The maximum N must be increased |
---|
| 1340 | { |
---|
| 1341 | #ifdef qdebug |
---|
| 1342 | if(dNn>iNmax) |
---|
| 1343 | { |
---|
| 1344 | G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1345 | iNmax=dNn; |
---|
| 1346 | } |
---|
| 1347 | if(dNn>iNran[z]) |
---|
| 1348 | { |
---|
| 1349 | G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1350 | iNran[z]=dNn; |
---|
| 1351 | } |
---|
| 1352 | #endif |
---|
| 1353 | return CalculateNuclMass(z,n,s); |
---|
| 1354 | } |
---|
| 1355 | } |
---|
| 1356 | else if(s==3) // ******************> S=3 nucleus |
---|
| 1357 | { |
---|
| 1358 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1359 | if(!iNin3[z]) // ====> This element is already initialized |
---|
| 1360 | { |
---|
| 1361 | #ifdef idebug |
---|
| 1362 | G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=3 is initialized. F="<<iNin3[z]<<G4endl; |
---|
| 1363 | #endif |
---|
| 1364 | G4int iNfin=iNL[z]; |
---|
| 1365 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1366 | for (G4int in=0; in<iNfin; in++) VZ3[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1367 | iNin3[z]=true; |
---|
| 1368 | } |
---|
| 1369 | G4int dNn=n-Nmin; |
---|
| 1370 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1371 | { |
---|
| 1372 | #ifdef qdebug |
---|
| 1373 | if(n<iNmin[z]) |
---|
| 1374 | { |
---|
| 1375 | G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=3 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1376 | iNmin[z]=n; |
---|
| 1377 | } |
---|
| 1378 | #endif |
---|
| 1379 | return CalculateNuclMass(z,n,s); |
---|
| 1380 | } |
---|
| 1381 | else if(dNn<iNL[z]) return VZ3[z][dNn]; // Found in DAM |
---|
| 1382 | else // --> The maximum N must be increased |
---|
| 1383 | { |
---|
| 1384 | #ifdef qdebug |
---|
| 1385 | if(dNn>iNmax) |
---|
| 1386 | { |
---|
| 1387 | G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1388 | iNmax=dNn; |
---|
| 1389 | } |
---|
| 1390 | if(dNn>iNran[z]) |
---|
| 1391 | { |
---|
| 1392 | G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1393 | iNran[z]=dNn; |
---|
| 1394 | } |
---|
| 1395 | #endif |
---|
| 1396 | return CalculateNuclMass(z,n,s); |
---|
| 1397 | } |
---|
| 1398 | } |
---|
| 1399 | else if(s==-4) // ******************> S=-4 nucleus |
---|
| 1400 | { |
---|
| 1401 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1402 | if(!iNin6[z]) // ====> This element is already initialized |
---|
| 1403 | { |
---|
| 1404 | #ifdef idebug |
---|
| 1405 | G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-4 is initialized. F="<<iNin6[z]<<G4endl; |
---|
| 1406 | #endif |
---|
| 1407 | G4int iNfin=iNL[z]; |
---|
| 1408 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1409 | for (G4int in=0; in<iNfin; in++) VZ6[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1410 | iNin6[z]=true; |
---|
| 1411 | } |
---|
| 1412 | G4int dNn=n-Nmin; |
---|
| 1413 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1414 | { |
---|
| 1415 | #ifdef qdebug |
---|
| 1416 | if(n<iNmin[z]) |
---|
| 1417 | { |
---|
| 1418 | G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-4 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1419 | iNmin[z]=n; |
---|
| 1420 | } |
---|
| 1421 | #endif |
---|
| 1422 | return CalculateNuclMass(z,n,s); |
---|
| 1423 | } |
---|
| 1424 | else if(dNn<iNL[z]) return VZ6[z][dNn]; // Found in DAM |
---|
| 1425 | else // --> The maximum N must be increased |
---|
| 1426 | { |
---|
| 1427 | #ifdef qdebug |
---|
| 1428 | if(dNn>iNmax) |
---|
| 1429 | { |
---|
| 1430 | G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1431 | iNmax=dNn; |
---|
| 1432 | } |
---|
| 1433 | if(dNn>iNran[z]) |
---|
| 1434 | { |
---|
| 1435 | G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1436 | iNran[z]=dNn; |
---|
| 1437 | } |
---|
| 1438 | #endif |
---|
| 1439 | return CalculateNuclMass(z,n,s); |
---|
| 1440 | } |
---|
| 1441 | } |
---|
| 1442 | else if(s==4) // ******************> S=4 nucleus |
---|
| 1443 | { |
---|
| 1444 | G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM) |
---|
| 1445 | if(!iNin4[z]) // ====> This element is already initialized |
---|
| 1446 | { |
---|
| 1447 | #ifdef idebug |
---|
| 1448 | G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=4 is initialized. F="<<iNin4[z]<<G4endl; |
---|
| 1449 | #endif |
---|
| 1450 | G4int iNfin=iNL[z]; |
---|
| 1451 | if(iNfin>iNR) iNfin=iNR; |
---|
| 1452 | for (G4int in=0; in<iNfin; in++) VZ4[z][in] = CalculateNuclMass(z,in+Nmin,s); |
---|
| 1453 | iNin4[z]=true; |
---|
| 1454 | } |
---|
| 1455 | G4int dNn=n-Nmin; |
---|
| 1456 | if(dNn<0) // --> The minimum N must be reduced |
---|
| 1457 | { |
---|
| 1458 | #ifdef qdebug |
---|
| 1459 | if(n<iNmin[z]) |
---|
| 1460 | { |
---|
| 1461 | G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=4 with N="<<n<<"<"<<iNmin[z]<<G4endl; |
---|
| 1462 | iNmin[z]=n; |
---|
| 1463 | } |
---|
| 1464 | #endif |
---|
| 1465 | return CalculateNuclMass(z,n,s); |
---|
| 1466 | } |
---|
| 1467 | else if(dNn<iNL[z]) return VZ4[z][dNn]; // Found in DAM |
---|
| 1468 | else // --> The maximum N must be increased |
---|
| 1469 | { |
---|
| 1470 | #ifdef qdebug |
---|
| 1471 | if(dNn>iNmax) |
---|
| 1472 | { |
---|
| 1473 | G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNmax<<G4endl; |
---|
| 1474 | iNmax=dNn; |
---|
| 1475 | } |
---|
| 1476 | if(dNn>iNran[z]) |
---|
| 1477 | { |
---|
| 1478 | G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNran[z]<<G4endl; |
---|
| 1479 | iNran[z]=dNn; |
---|
| 1480 | } |
---|
| 1481 | #endif |
---|
| 1482 | return CalculateNuclMass(z,n,s); |
---|
| 1483 | } |
---|
| 1484 | } |
---|
| 1485 | else |
---|
| 1486 | { |
---|
| 1487 | #ifdef qdebug |
---|
| 1488 | if(s<Smin || s>Smax) |
---|
| 1489 | { |
---|
| 1490 | if(s<Smin) Smin=s; |
---|
| 1491 | if(s>Smax) Smax=s; |
---|
| 1492 | G4cout<<">>G4QPDGCode::GetMass:Z="<<z<<" with S="<<s<<",N="<<n<<" (Improve)"<<G4endl; |
---|
| 1493 | } |
---|
| 1494 | #endif |
---|
| 1495 | rm=CalculateNuclMass(z,n,s); |
---|
| 1496 | } |
---|
| 1497 | #ifdef debug |
---|
| 1498 | G4cout<<"G4QPDGCode::GetMass:GetNuclMass="<<rm<<",Z="<<z<<",N="<<n<<",S="<<s<<G4endl; |
---|
| 1499 | #endif |
---|
| 1500 | return rm; |
---|
| 1501 | } |
---|
| 1502 | |
---|
| 1503 | // Calculate a nuclear mass for Z (a#of protons), N (a#of neutrons), & S (a#of lambdas) |
---|
| 1504 | G4double G4QPDGCode::CalculateNuclMass(G4int z, G4int n, G4int s) |
---|
| 1505 | // ========================================================= |
---|
| 1506 | { |
---|
| 1507 | static const G4double mP = QHaM(21); // Proton |
---|
| 1508 | static const G4double mN = QHaM(20); // Neutron |
---|
| 1509 | static const G4double mL = QHaM(22); // Lambda |
---|
| 1510 | static const G4double mD = G4Deuteron::Deuteron()->GetPDGMass(); // Deuteron (H-2) |
---|
| 1511 | static const G4double mT = G4Triton::Triton()->GetPDGMass(); // Triton (H-3) |
---|
| 1512 | static const G4double mHe3= G4He3::He3()->GetPDGMass(); // Hetrium (He-3) |
---|
| 1513 | static const G4double mAl = G4Alpha::Alpha()->GetPDGMass(); // Alpha (He-4) |
---|
| 1514 | static const G4double dmP = mP+mP; // DiProton |
---|
| 1515 | static const G4double dmN = mN+mN; // DiNeutron |
---|
| 1516 | static const G4double dmL = mL+mL; // DiLambda |
---|
| 1517 | static const G4double dLN = mL+mN; // LambdaNeutron |
---|
| 1518 | static const G4double dLP = mL+mP; // LambdaProton |
---|
| 1519 | static const G4double mSm = QHaM(23); // Sigma- |
---|
| 1520 | static const G4double mSp = QHaM(25); // Sigma+ |
---|
| 1521 | static const G4double dSP = mSp+mP; // ProtonSigma+ |
---|
| 1522 | static const G4double dSN = mSm+mN; // NeutronSigma- |
---|
| 1523 | static const G4double dnS = dSN+mN; // 2NeutronsSigma- |
---|
| 1524 | static const G4double dpS = dSP+mP; // 2ProtonsSigma+ |
---|
| 1525 | static const G4double mXm = QHaM(26); // Ksi- |
---|
| 1526 | static const G4double mXz = QHaM(27); // Ksi0 |
---|
| 1527 | static const G4double mOm = QHaM(44); // Omega- |
---|
| 1528 | static const G4double dXN = mXm+mN; // NeutronKsi- |
---|
| 1529 | static const G4double dXP = mXz+mP; // ProtonKsi0 |
---|
| 1530 | static const G4double dOP = mOm+mP; // ProtonOmega- |
---|
| 1531 | static const G4double dON = mOm+mN; // NeutronOmega- |
---|
| 1532 | static const G4double mK = QHaM(18); // Charged Kaon |
---|
| 1533 | static const G4double mK0 = QHaM(17); // Neutral Kaon |
---|
| 1534 | static const G4double mPi = QHaM(15); // Charged Pion |
---|
| 1535 | //////////static const G4double mPi0= QHaM(14); |
---|
| 1536 | //static const G4int nSh = 164; |
---|
| 1537 | //static G4double sh[nSh] = {0., // Fake element for C++ N=Z=0 |
---|
| 1538 | // -4.315548, 2.435504, -1.170501, 3.950887, 5.425238, |
---|
| 1539 | // 13.342524, 15.547586, 22.583129, 23.983480, 30.561036, |
---|
| 1540 | // 33.761971, 41.471027, 45.532156, 53.835880, 58.495514, |
---|
| 1541 | // 65.693445, 69.903344, 76.899581, 81.329361, 88.979438, |
---|
| 1542 | // 92.908703, 100.316636, 105.013393, 113.081686, 118.622601, |
---|
| 1543 | // 126.979113, 132.714435, 141.413182, 146.433488, 153.746754, |
---|
| 1544 | // 158.665225, 165.988967, 170.952395, 178.473011, 183.471531, |
---|
| 1545 | // 191.231310, 196.504414, 204.617158, 210.251108, 218.373984, |
---|
| 1546 | // 223.969281, 232.168660, 237.925619, 246.400505, 252.392471, |
---|
| 1547 | // 260.938644, 267.191321, 276.107788, 282.722682, 291.881502, |
---|
| 1548 | // 296.998590, 304.236025, 309.562296, 316.928655, 322.240263, |
---|
| 1549 | // 329.927236, 335.480630, 343.233705, 348.923475, 356.911659, |
---|
| 1550 | // 362.785757, 370.920926, 376.929998, 385.130316, 391.197741, |
---|
| 1551 | // 399.451554, 405.679971, 414.101869, 420.346260, 428.832412, |
---|
| 1552 | // 435.067445, 443.526983, 449.880034, 458.348602, 464.822352, |
---|
| 1553 | // 473.313779, 479.744524, 488.320887, 495.025069, 503.841579, |
---|
| 1554 | // 510.716379, 519.451976, 525.036156, 532.388151, 537.899017, |
---|
| 1555 | // 545.252264, 550.802469, 558.402181, 564.101100, 571.963429, |
---|
| 1556 | // 577.980340, 586.063802, 592.451334, 600.518525, 606.832027, |
---|
| 1557 | // 614.831626, 621.205330, 629.237413, 635.489106, 643.434167, |
---|
| 1558 | // 649.691284, 657.516479, 663.812101, 671.715021, 678.061128, |
---|
| 1559 | // 686.002970, 692.343712, 700.360477, 706.624091, 714.617848, |
---|
| 1560 | // 721.100390, 729.294717, 735.887170, 744.216084, 751.017094, |
---|
| 1561 | // 759.551764, 766.377807, 775.080204, 781.965673, 790.552795, |
---|
| 1562 | // 797.572494, 806.088030, 813.158751, 821.655631, 828.867137, |
---|
| 1563 | // 836.860955, 842.183292, 849.195302, 854.731798, 861.898839, |
---|
| 1564 | // 867.783606, 875.313342, 881.443441, 889.189065, 895.680189, |
---|
| 1565 | // 903.679729, 910.368085, 918.579876, 925.543547, 933.790028, |
---|
| 1566 | // 940.811396, 949.122548, 956.170201, 964.466810, 971.516490, |
---|
| 1567 | // 979.766905, 986.844659, 995.113552,1002.212760,1010.418770, |
---|
| 1568 | // 1018.302560,1025.781870,1033.263560,1040.747880,1048.234460, |
---|
| 1569 | // 1055.723430,1063.214780,1070.708750,1078.204870,1085.703370, |
---|
| 1570 | // 1093.204260,1100.707530,1108.213070}; |
---|
| 1571 | //static const G4double b1=8.09748564; // MeV |
---|
| 1572 | //static const G4double b2=-0.76277387; |
---|
| 1573 | //static const G4double b3=83.487332; // MeV |
---|
| 1574 | //static const G4double b4=0.090578206;// 2*b4 |
---|
| 1575 | //static const G4double b5=0.676377211;// MeV |
---|
| 1576 | //static const G4double b6=5.55231981; // MeV |
---|
| 1577 | static const G4double b7=25.; // MeV (Lambda binding energy predexponent) |
---|
| 1578 | // --- even-odd difference is 3.7(MeV)/X |
---|
| 1579 | // --- S(X>151)=-57.56-5.542*X**1.05 |
---|
| 1580 | static const G4double b8=10.5; // (Lambda binding energy exponent) |
---|
| 1581 | //static const G4double b9=-1./3.; |
---|
| 1582 | static const G4double a2=0.13; // LambdaBindingEnergy for deutron+LambdaSystem(MeV) |
---|
| 1583 | static const G4double a3=2.2; // LambdaBindingEnergy for (t/He3)+LambdaSystem(MeV) |
---|
| 1584 | static const G4double um=931.49432; // Umified atomic mass unit (MeV) |
---|
| 1585 | //static const G4double me =0.511; // electron mass (MeV) :: n:8.071, p:7.289 |
---|
| 1586 | static const G4double eps =0.0001; // security addition for multybaryons |
---|
| 1587 | //static G4double c[9][9]={// z=1 =2 =3 =4 =5 =6 =7 =8 =9 |
---|
| 1588 | // {13.136,14.931,25.320,38.000,45.000,55.000,65.000,75.000,85.000},//n=1 |
---|
| 1589 | // {14.950, 2.425,11.680,18.374,27.870,35.094,48.000,60.000,72.000}, //n=2 |
---|
| 1590 | // {25.930,11.390,14.086,15.770,22.921,28.914,39.700,49.000,60.000}, //n=3 |
---|
| 1591 | // {36.830,17.594,14.908, 4.942,12.416,15.699,24.960,32.048,45.000}, //n=4 |
---|
| 1592 | // {41.860,26.110,20.946,11.348,12.051,10.650,17.338,23.111,33.610}, //n=5 |
---|
| 1593 | // {45.000,31.598,24.954,12.607, 8.668, 0.000, 5.345, 8.006,16.780}, //n=6 |
---|
| 1594 | // {50.000,40.820,33.050,20.174,13.369, 3.125, 2.863, 2.855,10.680}, //n=7 |
---|
| 1595 | // {55.000,48.810,40.796,25.076,16.562, 3.020, 0.101,-4.737,1.9520}, //n=8 |
---|
| 1596 | // {60.000,55.000,50.100,33.660,23.664, 9.873, 5.683,-0.809,0.8730}}; //n=9 |
---|
| 1597 | if(z>107) |
---|
| 1598 | { |
---|
| 1599 | #ifdef debug |
---|
| 1600 | G4cout<<"***G4QPDGCode::CalcNuclMass: Z="<<z<<">107, N="<<n<<", S="<<s<<G4endl; |
---|
| 1601 | #endif |
---|
| 1602 | return 256000.; |
---|
| 1603 | } |
---|
| 1604 | G4int Z=z; |
---|
| 1605 | G4int N=n; |
---|
| 1606 | G4int S=s; |
---|
| 1607 | #ifdef debug |
---|
| 1608 | G4cout<<"G4QPDGCode::CalcNuclMass called with Z="<<Z<<",N="<<N<<", S="<<S<<G4endl; |
---|
| 1609 | #endif |
---|
| 1610 | if(!N&&!Z&&!S) |
---|
| 1611 | { |
---|
| 1612 | #ifdef debug |
---|
| 1613 | //G4cout<<"G4QPDGCode::GetNuclMass(0,0,0)="<<mPi0<<"#0"<<G4endl; |
---|
| 1614 | #endif |
---|
| 1615 | //return mPi0; // Pi0 mass (dangerose) |
---|
| 1616 | return 0.; // Photon mass |
---|
| 1617 | } |
---|
| 1618 | else if(!N&&!Z&&S>1) return mL*S+eps; |
---|
| 1619 | else if(!N&&Z>1&&!S) return mP*Z+eps; |
---|
| 1620 | else if(N>1&&!Z&&!S) return mN*N+eps; |
---|
| 1621 | G4int A=Z+N; |
---|
| 1622 | G4int Bn=A+S; |
---|
| 1623 | //if((Z<0||N<0)&&!Bn) |
---|
| 1624 | //{ |
---|
| 1625 | // if (N<0) return Bn*mL-Z*mK - N*mK0+eps* S ; |
---|
| 1626 | // else return Bn*mL+N*mPi-A*mK +eps*(N+S); |
---|
| 1627 | //} |
---|
| 1628 | //if(A<0&&Bn>=0) // Bn*LAMBDA's + (-(N+Z))*antiKaons |
---|
| 1629 | //{ |
---|
| 1630 | // if (N<0&&Z<0) return Bn*mL-Z*mK -N*mK0+eps* S ; |
---|
| 1631 | // else if(N<0) return Bn*mL+Z*mPi-A*mK0+eps*(Z+S); |
---|
| 1632 | // else return Bn*mL+N*mPi-A*mK +eps*(N+S); |
---|
| 1633 | //} |
---|
| 1634 | if(!Bn) // => "GS Mesons - octet" case (without eta&pi0) |
---|
| 1635 | { |
---|
| 1636 | if (!S&&Z<0) return mPi*N; |
---|
| 1637 | else if(!S&&N<0) return mPi*Z; |
---|
| 1638 | else if ( (N == 1 && S == -1) || (N == -1 && S == 1) ) |
---|
| 1639 | return mK0; // Simple decision |
---|
| 1640 | else if ( (S == 1 && Z == -1) || (S == -1 && Z == 1) ) |
---|
| 1641 | return mK; // Simple decision |
---|
| 1642 | else if(S>0) // General decision |
---|
| 1643 | { |
---|
| 1644 | if (-Z>S) return S*mK-(S+Z)*mPi+eps; |
---|
| 1645 | else if(Z>=0) return S*mK0+Z*mPi+eps; |
---|
| 1646 | else return (S+Z)*mK0-Z*mK+eps; |
---|
| 1647 | } |
---|
| 1648 | else if(S<0) // General decision |
---|
| 1649 | { |
---|
| 1650 | if (Z>-S) return -S*mK+(S+Z)*mPi+eps; |
---|
| 1651 | else if(Z<=0) return -S*mK0-Z*mPi+eps; |
---|
| 1652 | else return -(S+Z)*mK0+Z*mK+eps; |
---|
| 1653 | } |
---|
| 1654 | } |
---|
| 1655 | else if(Bn==1) // => "GS Baryons - octet" case (withoutSigma0) |
---|
| 1656 | { |
---|
| 1657 | if (Z== 1 && N== 0 && S== 0) return mP; |
---|
| 1658 | else if(Z== 0 && N== 1 && S== 0) return mN; |
---|
| 1659 | else if(Z== 0 && N== 0 && S== 1) return mL; |
---|
| 1660 | else if(Z== 1 && N==-1 && S== 1) return mSp; // Lower than Sigma+ (Simp.Decis) |
---|
| 1661 | else if(Z==-1 && N== 1 && S== 1) return mSm; // Lower than Sigma- (Simp.Decis) |
---|
| 1662 | else if(Z== 0 && N==-1 && S== 2) return mXz; // Lower than Xi0 (Simp.Decis) |
---|
| 1663 | else if(Z==-1 && N== 0 && S== 2) return mXm; // Lower than Xi- (Simp.Decis) |
---|
| 1664 | else if(Z==-1 && N==-1 && S== 3) return mOm; // Lower than Omega- (Simp.Decis) |
---|
| 1665 | else if(!S&&Z<0) return mN-mPi*Z+eps; // Negative Isonuclei |
---|
| 1666 | else if(!S&&N<0) return mP-mPi*N+eps; // Positive Isonuclei |
---|
| 1667 | else if(S==1) // --> General decision |
---|
| 1668 | { |
---|
| 1669 | if (N>1) return mSm+(N-1)*mPi+eps; // (Sigma-)+(n*PI-) |
---|
| 1670 | else if(Z>1) return mSp+(Z-1)*mPi+eps; // (Sigma+)+(n*PI+) |
---|
| 1671 | } |
---|
| 1672 | else if(S==2) // --> General decision |
---|
| 1673 | { |
---|
| 1674 | if (N>0) return mXm+N*mPi+eps; // (Xi-)+(n*PI-) |
---|
| 1675 | else if(Z>0) return mXz+Z*mPi+eps; // (Xi0)+(n*PI+) |
---|
| 1676 | } |
---|
| 1677 | else if(S==3) // --> General decision |
---|
| 1678 | { |
---|
| 1679 | if (N>-1) return mOm+(N+1)*mPi+eps; // (Omega-)+(n*PI-) |
---|
| 1680 | else if(Z>-1) return mOm+(Z+1)*mPi+eps; // (Omega-)+(n*PI+) |
---|
| 1681 | } |
---|
| 1682 | else if(S>3) // --> General Omega- decision |
---|
| 1683 | { |
---|
| 1684 | if (-Z>S-2) return mOm+(S-3)*mK +(2-Z-S)*mPi+eps; |
---|
| 1685 | else if(Z>-1) return mOm+(S-3)*mK0+(Z+1)+mPi+eps; |
---|
| 1686 | else return mOm+(S+Z-2)*mK0-(Z+1)*mK+eps; |
---|
| 1687 | } |
---|
| 1688 | } |
---|
| 1689 | else if(Bn==2) // => "GS Baryons - decuplet" case (NP,LP, and LN are made below) |
---|
| 1690 | { |
---|
| 1691 | if (Z== 2 && N== 0 && S== 0) return dmP; |
---|
| 1692 | //else if(Z== 1 && N== 1 && S== 0) return 1875.61282; // Exact deuteron PDG Mass |
---|
| 1693 | else if(Z== 1 && N== 1 && S== 0) return mD; // Exact deuteron PDG Mass |
---|
| 1694 | else if(Z== 0 && N== 2 && S== 0) return dmN; |
---|
| 1695 | else if(Z== 2 && N==-1 && S== 1) return dSP; |
---|
| 1696 | else if(Z== 1 && N== 0 && S== 1) return dLP; |
---|
| 1697 | else if(Z== 0 && N== 1 && S== 1) return dLN; |
---|
| 1698 | else if(Z==-1 && N== 2 && S== 1) return dSN; |
---|
| 1699 | else if(Z== 1 && N==-1 && S== 2) return dXP; |
---|
| 1700 | else if(Z== 0 && N== 0 && S== 2) return dmL; |
---|
| 1701 | else if(Z==-1 && N== 1 && S== 2) return dXN; |
---|
| 1702 | else if(Z== 0 && N==-1 && S== 3) return dOP; |
---|
| 1703 | else if(Z==-1 && N== 0 && S== 3) return dON; |
---|
| 1704 | else if(!S&&Z<0) return dmN-mPi*Z+eps; // Negative Isonuclei |
---|
| 1705 | else if(!S&&N<0) return dmP-mPi*N+eps; // Positive Isonuclei |
---|
| 1706 | else if(S==1) // --> General decision |
---|
| 1707 | { |
---|
| 1708 | if (N>2) return dSP+(N-2)*mPi+eps; // (nSigma-)+(n*PI-) |
---|
| 1709 | else if(Z>2) return dSN+(Z-1)*mPi+eps; // (pSigma+)+(n*PI+) |
---|
| 1710 | } |
---|
| 1711 | else if(S==2) // --> General decision |
---|
| 1712 | { |
---|
| 1713 | if (N>1) return dXN+(N-1)*mPi+eps; // (nXi-)+(n*PI-) |
---|
| 1714 | else if(Z>1) return dXP+(Z-1)*mPi+eps; // (pXi0)+(n*PI+) |
---|
| 1715 | } |
---|
| 1716 | else if(S==3) // --> General decision |
---|
| 1717 | { |
---|
| 1718 | if (N>0) return dON+N*mPi+eps; // (nOmega-)+(n*PI-) |
---|
| 1719 | else if(Z>0) return dOP+Z*mPi+eps; // (pOmega-)+(n*PI+) |
---|
| 1720 | } |
---|
| 1721 | else if(S>3) // --> General Omega- decision |
---|
| 1722 | { |
---|
| 1723 | if (-Z>S-2) return dON+(S-3)*mK +(2-Z-S)*mPi+eps; |
---|
| 1724 | else if(Z>0) return dOP+(S-3)*mK0+Z+mPi+eps; |
---|
| 1725 | else return dOP+(S+Z-3)*mK0-Z*mK+eps; |
---|
| 1726 | } |
---|
| 1727 | //else if(S>0) // @@ Implement General Decision |
---|
| 1728 | //{ |
---|
| 1729 | // //#ifdef debug |
---|
| 1730 | // G4cout<<"***G4QPDGCode::GetNuclMass:B=2, Z="<<Z<<",N="<<N<<",S="<<S<<G4endl; |
---|
| 1731 | // //#endif |
---|
| 1732 | // return bigM; // Exotic dibaryons (?) |
---|
| 1733 | //} |
---|
| 1734 | } |
---|
| 1735 | else if(Bn==3) |
---|
| 1736 | { |
---|
| 1737 | if(!S) // Bn=3 |
---|
| 1738 | { |
---|
| 1739 | if (Z==1 && N== 2) return mT; // tritium |
---|
| 1740 | else if(Z==2 && N== 1) return mHe3; // hetrium |
---|
| 1741 | } |
---|
| 1742 | if(S== 1 && Z==-1 && N== 3) // nnSig- |
---|
| 1743 | { |
---|
| 1744 | if (Z==-1 && N== 3) return dnS; // nnSig- |
---|
| 1745 | else if(Z== 3 && N==-1) return dpS; // ppSig+ |
---|
| 1746 | } |
---|
| 1747 | } |
---|
| 1748 | else if(!S) |
---|
| 1749 | { |
---|
| 1750 | if(Z==2 && N==2) return mAl; // Alpha |
---|
| 1751 | else if(Z<0) return A*mN-Z*mPi+eps; // Multybaryon Negative Isonuclei |
---|
| 1752 | else if(Z>A) return A*mP+(Z-A)*mPi+eps; // Multybaryon Positive Isonuclei |
---|
| 1753 | } |
---|
| 1754 | // === Start mesonic extraction === |
---|
| 1755 | G4double km=0.; // Mass Sum of K mesons (G4E::DecayAntiStrang.) |
---|
| 1756 | G4int Zm=Z; |
---|
| 1757 | G4int Nm=N; |
---|
| 1758 | G4int Sm=S; |
---|
| 1759 | if(S<0&&Bn>0) // NEW: the free mass minimum |
---|
| 1760 | { |
---|
| 1761 | if(Zm>=-S) // Enough charge for K+'s |
---|
| 1762 | { |
---|
| 1763 | km=-S*mK; // Anti-Lambdas are compensated by protons |
---|
| 1764 | Zm+=S; |
---|
| 1765 | } |
---|
| 1766 | else if(Zm>0) |
---|
| 1767 | { |
---|
| 1768 | km=Zm*mK-(S+Zm)*mK0; // Anti-Lambdas are partially compensated by neutrons |
---|
| 1769 | Zm=0; |
---|
| 1770 | Nm+=S+Zm; |
---|
| 1771 | } |
---|
| 1772 | } |
---|
| 1773 | else Sm=0; // No alternative calculations |
---|
| 1774 | // Old algorithm |
---|
| 1775 | G4double k=0.; // Mass Sum of K mesons |
---|
| 1776 | if(S<0&&Bn>0) // OLD @@ Can be improved by K0/K+ search of minimum |
---|
| 1777 | { |
---|
| 1778 | G4int sH=(-S)/2; // SmallHalfS || The algorithm must be the same |
---|
| 1779 | G4int bH=-S-sH; // BigHalhS || as in G4QE::DecayAntiStrange |
---|
| 1780 | if(Z>0 && Z>N) |
---|
| 1781 | { |
---|
| 1782 | if(Z>=bH) // => "Enough protons in nucleus" case |
---|
| 1783 | { |
---|
| 1784 | if(N>=sH) |
---|
| 1785 | { |
---|
| 1786 | k=bH*mK+sH*mK0; |
---|
| 1787 | Z-=bH; |
---|
| 1788 | N-=sH; |
---|
| 1789 | } |
---|
| 1790 | else |
---|
| 1791 | { |
---|
| 1792 | G4int dN=Z-N; |
---|
| 1793 | if(dN>=-S) |
---|
| 1794 | { |
---|
| 1795 | k=-S*mK; |
---|
| 1796 | Z+=S; |
---|
| 1797 | } |
---|
| 1798 | else |
---|
| 1799 | { |
---|
| 1800 | G4int sS=(-S-dN)/2; |
---|
| 1801 | G4int bS=-S-dN-sS; |
---|
| 1802 | sS+=dN; |
---|
| 1803 | if(Z>=sS&&N>=bS) |
---|
| 1804 | { |
---|
| 1805 | k=sS*mK+bS*mK0; |
---|
| 1806 | Z-=sS; |
---|
| 1807 | N-=bS; |
---|
| 1808 | } |
---|
| 1809 | else if(Z<sS) |
---|
| 1810 | { |
---|
| 1811 | G4int dS=-S-Z; |
---|
| 1812 | k=Z*mK+dS*mK0; |
---|
| 1813 | N-=dS; |
---|
| 1814 | Z=0; |
---|
| 1815 | } |
---|
| 1816 | else |
---|
| 1817 | { |
---|
| 1818 | G4int dS=-S-N; |
---|
| 1819 | k=dS*mK+N*mK0; |
---|
| 1820 | Z-=dS; |
---|
| 1821 | N=0; |
---|
| 1822 | } |
---|
| 1823 | } |
---|
| 1824 | } |
---|
| 1825 | } |
---|
| 1826 | else // Must not be here |
---|
| 1827 | { |
---|
| 1828 | #ifdef debug |
---|
| 1829 | G4cout<<"***G4QPDGC::CalcNuclMass:Antimatter? Z="<<Z<<",N="<<N<<",S="<<S<<G4endl; |
---|
| 1830 | #endif |
---|
| 1831 | return 0.; // @@ Antiparticles aren't implemented @@ |
---|
| 1832 | } |
---|
| 1833 | } |
---|
| 1834 | else if(N>=bH) |
---|
| 1835 | { |
---|
| 1836 | if(Z>=sH) |
---|
| 1837 | { |
---|
| 1838 | k=sH*mK+bH*mK0; |
---|
| 1839 | Z-=sH; |
---|
| 1840 | N-=bH; |
---|
| 1841 | } |
---|
| 1842 | else |
---|
| 1843 | { |
---|
| 1844 | G4int dN=N-Z; |
---|
| 1845 | if(dN>=-S) |
---|
| 1846 | { |
---|
| 1847 | k=-S*mK0; |
---|
| 1848 | N+=S; |
---|
| 1849 | } |
---|
| 1850 | else |
---|
| 1851 | { |
---|
| 1852 | G4int sS=(-S-dN)/2; |
---|
| 1853 | G4int bS=-S-dN-sS; |
---|
| 1854 | bS+=dN; |
---|
| 1855 | if(N>=bS&&Z>=sS) |
---|
| 1856 | { |
---|
| 1857 | k=sS*mK+bS*mK0; |
---|
| 1858 | Z-=sS; |
---|
| 1859 | N-=bS; |
---|
| 1860 | } |
---|
| 1861 | else if(N<bS) |
---|
| 1862 | { |
---|
| 1863 | G4int dS=-S-N; |
---|
| 1864 | k=dS*mK+N*mK0; |
---|
| 1865 | Z-=dS; |
---|
| 1866 | N=0; |
---|
| 1867 | } |
---|
| 1868 | else |
---|
| 1869 | { |
---|
| 1870 | G4int dS=-S-Z; |
---|
| 1871 | k=Z*mK+dS*mK0; |
---|
| 1872 | N-=dS; |
---|
| 1873 | Z=0; |
---|
| 1874 | } |
---|
| 1875 | } |
---|
| 1876 | } |
---|
| 1877 | } |
---|
| 1878 | else // Must not be here |
---|
| 1879 | { |
---|
| 1880 | return 0.; // @@ Antiparticles aren't implemented @@ |
---|
| 1881 | #ifdef debug |
---|
| 1882 | G4cout<<"***G4QPDGC::CalcNuclMass:Antimatter? N="<<N<<",Z="<<Z<<",S="<<S<<G4endl; |
---|
| 1883 | #endif |
---|
| 1884 | } |
---|
| 1885 | S=0; |
---|
| 1886 | } |
---|
| 1887 | if(N<0) |
---|
| 1888 | { |
---|
| 1889 | k+=-N*mPi; |
---|
| 1890 | Z+=N; |
---|
| 1891 | N=0; |
---|
| 1892 | } |
---|
| 1893 | if(Z<0) |
---|
| 1894 | { |
---|
| 1895 | k+=-Z*mPi; |
---|
| 1896 | N+=Z; |
---|
| 1897 | Z=0; |
---|
| 1898 | } |
---|
| 1899 | A=Z+N; |
---|
| 1900 | if (!A) return k+S*mL+S*eps; // @@ multy LAMBDA states are not implemented |
---|
| 1901 | G4double m=k+A*um; // Expected mass in atomic units |
---|
| 1902 | //G4double D=N-Z; // Isotopic shift of the nucleus |
---|
| 1903 | if ( (A+S < 1 && k==0.) || Z < 0 || N < 0 ) |
---|
| 1904 | { // @@ Can be generalized to anti-nuclei |
---|
| 1905 | #ifdef debug |
---|
| 1906 | G4cout<<"**G4QPDGCode::CalcNuclMass:A="<<A<<"<1 || Z="<<Z<<"<0 || N="<<N<<"<0"<<G4endl; |
---|
| 1907 | //@@throw G4QException("***G4QPDGCode::GetNuclMass: Impossible nucleus"); |
---|
| 1908 | #endif |
---|
| 1909 | return 0.; // @@ Temporary |
---|
| 1910 | } |
---|
| 1911 | if (!Z) return k+N*(mN+.1)+S*(mL+.1); // @@ n+LAMBDA states are not implemented |
---|
| 1912 | else if(!N) return k+Z*(mP+1.)+S*(mL+.1); // @@ p+LAMBDA states are not implemented |
---|
| 1913 | //else if(N<=9&&Z<=9) m+=1.433e-5*pow(double(Z),2.39)-Z*me+c[N-1][Z-1];// Geant4 Comp.now |
---|
| 1914 | else |
---|
| 1915 | { |
---|
| 1916 | //G4double fA=A; |
---|
| 1917 | // IsInTable is inside G4NucleiProperties::GetNuclearMass |
---|
| 1918 | // G4ParticleTable::GetParticleTable()->FindIon(Zm,Am,0,Zm) creates new Ion! |
---|
| 1919 | if(A==256 && Z==128) m=256000.; |
---|
| 1920 | else m=k+G4NucleiProperties::GetNuclearMass(A,Z); |
---|
| 1921 | //if(G4NucleiPropertiesTable::IsInTable(Z,A)) |
---|
| 1922 | // m=k+G4NucleiProperties::GetNuclearMass(A,Z); |
---|
| 1923 | //else if(A==256 && Z==128) m=256000.; |
---|
| 1924 | //else |
---|
| 1925 | // m=k+G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z)->GetPDGMass(); |
---|
| 1926 | //m+=-sh[Z]-sh[N]+b1*D*D*pow(fA,b2)+b3*(1.-2./(1.+exp(b4*D)))+Z*Z*(b5*pow(fA,b9)+b6/fA); |
---|
| 1927 | } |
---|
| 1928 | //@@//G4double maxM= k+Z*mP+N*mN+S*mL+eps; // @@ eps -- Wings of the Mass parabola |
---|
| 1929 | //@@//if(m>maxM) m=maxM; |
---|
| 1930 | G4double mm=m; |
---|
| 1931 | if(Sm<0) // For the new algorithm of calculation |
---|
| 1932 | { |
---|
| 1933 | if(Nm<0) |
---|
| 1934 | { |
---|
| 1935 | km+=-Nm*mPi; |
---|
| 1936 | Zm+=Nm; |
---|
| 1937 | Nm=0; |
---|
| 1938 | } |
---|
| 1939 | if(Zm<0) |
---|
| 1940 | { |
---|
| 1941 | km+=-Zm*mPi; |
---|
| 1942 | Nm+=Zm; |
---|
| 1943 | Zm=0; |
---|
| 1944 | } |
---|
| 1945 | G4int Am=Zm+Nm; |
---|
| 1946 | if(!Am) return km+eps; |
---|
| 1947 | mm=km+Am*um; // Expected mass in atomic units |
---|
| 1948 | //G4double Dm=Nm-Zm; // Isotopic shift of the nucleus |
---|
| 1949 | if ( (Am < 1 && km==0.) || Zm < 0 || Nm < 0 ) |
---|
| 1950 | { // @@ Can be generalized to anti-nuclei |
---|
| 1951 | #ifdef debug |
---|
| 1952 | G4cerr<<"**G4QPDGCode::CalcNucM:A="<<Am<<"<1 || Z="<<Zm<<"<0 || N="<<Nm<<"<0"<<G4endl; |
---|
| 1953 | #endif |
---|
| 1954 | } |
---|
| 1955 | if (!Zm) return km+Nm*(mN+.1); |
---|
| 1956 | else if(!Nm) return km+Zm*(mP+1.); |
---|
| 1957 | //else if(Nm<=9&&Zm<=9) mm+=1.433e-5*pow(double(Zm),2.39)-Zm*me+c[Nm-1][Zm-1];// Geant4 |
---|
| 1958 | else |
---|
| 1959 | { |
---|
| 1960 | // IsInTable is inside G4NucleiProperties::GetNuclearMass |
---|
| 1961 | // G4ParticleTable::GetParticleTable()->FindIon(Zm,Am,0,Zm) creates new Ion! |
---|
| 1962 | mm=km+G4NucleiProperties::GetNuclearMass(Am,Zm); |
---|
| 1963 | //G4double fA=Am; |
---|
| 1964 | //if(G4NucleiPropertiesTable::IsInTable(Zm,Am)) |
---|
| 1965 | // mm=km+G4NucleiProperties::GetNuclearMass(Am,Zm); |
---|
| 1966 | //else |
---|
| 1967 | // mm=km+G4ParticleTable::GetParticleTable()->FindIon(Zm,Am,0,Zm)->GetPDGMass(); |
---|
| 1968 | // //mm+=-sh[Zm]-sh[Nm]+b1*Dm*Dm*pow(fA,b2)+b3*(1.-2./(1.+exp(b4*Dm))) |
---|
| 1969 | // // +Zm*Zm*(b5*pow(fA,b9)+b6/Am); |
---|
| 1970 | } |
---|
| 1971 | //@@//G4double mM= km+Zm*mP+Nm*mN+eps; |
---|
| 1972 | //@@//if(mm>mM) mm=mM; |
---|
| 1973 | } |
---|
| 1974 | if(m>mm) m=mm; |
---|
| 1975 | if(S>0) |
---|
| 1976 | { |
---|
| 1977 | G4double bs=0.; |
---|
| 1978 | if (A==2) bs=a2; |
---|
| 1979 | else if(A==3) bs=a3; |
---|
| 1980 | else if(A>3) bs=b7*exp(-b8/(A+1.)); |
---|
| 1981 | m+=S*(mL-bs); |
---|
| 1982 | } |
---|
| 1983 | #ifdef debug |
---|
| 1984 | G4cout<<"G4QPDGCode::CalcNuclMass: >>>OUT<<< m="<<m<<G4endl; |
---|
| 1985 | #endif |
---|
| 1986 | //if(z==8&&n==9&&!s) G4cout<<"G4QPDGC::GetNucM:m="<<m<<",mm="<<mm<<G4endl; |
---|
| 1987 | return m; |
---|
| 1988 | } |
---|
| 1989 | |
---|
| 1990 | // Get a Quark Content {d,u,s,ad,au,as} for the PDG |
---|
| 1991 | G4QContent G4QPDGCode::GetQuarkContent() const |
---|
| 1992 | // ====================================== |
---|
| 1993 | { |
---|
| 1994 | G4int a=0; // particle |
---|
| 1995 | if(thePDGCode<0) a=1; // anti-particle |
---|
| 1996 | G4int d=0; |
---|
| 1997 | G4int u=0; |
---|
| 1998 | G4int s=0; |
---|
| 1999 | G4int ad=0; |
---|
| 2000 | G4int au=0; |
---|
| 2001 | G4int as=0; |
---|
| 2002 | G4int ab=abs(thePDGCode); |
---|
| 2003 | if (!ab) |
---|
| 2004 | { |
---|
| 2005 | G4cerr<<"***G4QPDGCode::GetQuarkContent: PDG=0, return (0,0,0,0,0,0)"<<G4endl; |
---|
| 2006 | return G4QContent(0,0,0,0,0,0); |
---|
| 2007 | } |
---|
| 2008 | else if(ab<4) // @@ for SU(6) : ab<7 |
---|
| 2009 | { |
---|
| 2010 | if (thePDGCode== 1) return G4QContent(1,0,0,0,0,0); |
---|
| 2011 | else if(thePDGCode== 2) return G4QContent(0,1,0,0,0,0); |
---|
| 2012 | else if(thePDGCode== 3) return G4QContent(0,0,1,0,0,0); |
---|
| 2013 | else if(thePDGCode==-1) return G4QContent(0,0,0,1,0,0); |
---|
| 2014 | else if(thePDGCode==-2) return G4QContent(0,0,0,0,1,0); |
---|
| 2015 | else if(thePDGCode==-3) return G4QContent(0,0,0,0,0,1); |
---|
| 2016 | } |
---|
| 2017 | else if(ab<99) |
---|
| 2018 | { |
---|
| 2019 | #ifdef debug |
---|
| 2020 | if (ab==22) G4cout<<"-W-G4QPDGC::GetQuarkCont: For the Photon? - Return 0"<<G4endl; |
---|
| 2021 | else if(ab==10) G4cout<<"-W-G4QPDGC::GetQuarkCont: For Chipolino? - Return 0"<<G4endl; |
---|
| 2022 | else G4cout<<"-W-G4QPDGCode::GetQuarkCont: For PDG="<<thePDGCode<<" Return 0"<<G4endl; |
---|
| 2023 | #endif |
---|
| 2024 | return G4QContent(0,0,0,0,0,0); // Photon, bosons, leptons |
---|
| 2025 | } |
---|
| 2026 | else if(ab<80000000) // Baryons & Mesons |
---|
| 2027 | { |
---|
| 2028 | G4int c=ab/10; // temporary (quarks) |
---|
| 2029 | G4int f=c%10; // (1) low quark |
---|
| 2030 | G4int v=c/10; // (2,3) temporary(B) or (2) final(M) (high quarks, high quark) |
---|
| 2031 | G4int t=0; // (3)prototype of highest quark (B) |
---|
| 2032 | #ifdef sdebug |
---|
| 2033 | G4cout<<"G4QPDGCode::GetQuarkContent: a="<<ab<<", c="<<c<<", f="<<f<<", v="<<v<<G4endl; |
---|
| 2034 | #endif |
---|
| 2035 | if(v>10) // Baryons |
---|
| 2036 | { |
---|
| 2037 | t=v/10; // (3) highest quark |
---|
| 2038 | v%=10; // (2) high quark |
---|
| 2039 | if (f==1) |
---|
| 2040 | { |
---|
| 2041 | if(a) ad++; |
---|
| 2042 | else d++; |
---|
| 2043 | } |
---|
| 2044 | else if(f==2) |
---|
| 2045 | { |
---|
| 2046 | if(a) au++; |
---|
| 2047 | else u++; |
---|
| 2048 | } |
---|
| 2049 | else if(f==3) |
---|
| 2050 | { |
---|
| 2051 | if(a) as++; |
---|
| 2052 | else s++; |
---|
| 2053 | } |
---|
| 2054 | else G4cerr<<"***G4QPDGCode::GetQContent:1 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl; |
---|
| 2055 | if (v==1) |
---|
| 2056 | { |
---|
| 2057 | if(a) ad++; |
---|
| 2058 | else d++; |
---|
| 2059 | } |
---|
| 2060 | else if(v==2) |
---|
| 2061 | { |
---|
| 2062 | if(a) au++; |
---|
| 2063 | else u++; |
---|
| 2064 | } |
---|
| 2065 | else if(v==3) |
---|
| 2066 | { |
---|
| 2067 | if(a) as++; |
---|
| 2068 | else s++; |
---|
| 2069 | } |
---|
| 2070 | else G4cerr<<"***G4QPDGCode::GetQContent:2 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl; |
---|
| 2071 | if (t==1) |
---|
| 2072 | { |
---|
| 2073 | if(a) ad++; |
---|
| 2074 | else d++; |
---|
| 2075 | } |
---|
| 2076 | else if(t==2) |
---|
| 2077 | { |
---|
| 2078 | if(a) au++; |
---|
| 2079 | else u++; |
---|
| 2080 | } |
---|
| 2081 | else if(t==3) |
---|
| 2082 | { |
---|
| 2083 | if(a) as++; |
---|
| 2084 | else s++; |
---|
| 2085 | } |
---|
| 2086 | else G4cerr<<"***G4QPDGCode::GetQCont:3 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl; |
---|
| 2087 | return G4QContent(d,u,s,ad,au,as); |
---|
| 2088 | } |
---|
| 2089 | else // Mesons |
---|
| 2090 | { |
---|
| 2091 | if(f==v) |
---|
| 2092 | { |
---|
| 2093 | if (f==1) return G4QContent(1,0,0,1,0,0); |
---|
| 2094 | else if(f==2) return G4QContent(0,1,0,0,1,0); |
---|
| 2095 | else if(f==3) return G4QContent(0,0,1,0,0,1); |
---|
| 2096 | else G4cerr<<"***G4QPDGCode::GetQCont:4 PDG="<<thePDGCode<<",i="<<f<<","<<v<<","<<t<<G4endl; |
---|
| 2097 | } |
---|
| 2098 | else |
---|
| 2099 | { |
---|
| 2100 | if (f==1 && v==2) |
---|
| 2101 | { |
---|
| 2102 | if(a)return G4QContent(1,0,0,0,1,0); |
---|
| 2103 | else return G4QContent(0,1,0,1,0,0); |
---|
| 2104 | } |
---|
| 2105 | else if(f==1 && v==3) |
---|
| 2106 | { |
---|
| 2107 | if(a)return G4QContent(0,0,1,1,0,0); |
---|
| 2108 | else return G4QContent(1,0,0,0,0,1); |
---|
| 2109 | } |
---|
| 2110 | else if(f==2 && v==3) |
---|
| 2111 | { |
---|
| 2112 | if(a)return G4QContent(0,0,1,0,1,0); |
---|
| 2113 | else return G4QContent(0,1,0,0,0,1); |
---|
| 2114 | } |
---|
| 2115 | else G4cerr<<"***G4QPDGCode::GetQCont:5 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl; |
---|
| 2116 | } |
---|
| 2117 | } |
---|
| 2118 | } |
---|
| 2119 | else |
---|
| 2120 | { |
---|
| 2121 | G4int szn=ab-90000000; |
---|
| 2122 | G4int ds=0; |
---|
| 2123 | G4int dz=0; |
---|
| 2124 | G4int dn=0; |
---|
| 2125 | if(szn<-100000) |
---|
| 2126 | { |
---|
| 2127 | G4int ns=(-szn)/1000000+1; |
---|
| 2128 | szn+=ns*1000000; |
---|
| 2129 | ds+=ns; |
---|
| 2130 | } |
---|
| 2131 | else if(szn<-100) |
---|
| 2132 | { |
---|
| 2133 | G4int nz=(-szn)/1000+1; |
---|
| 2134 | szn+=nz*1000; |
---|
| 2135 | dz+=nz; |
---|
| 2136 | } |
---|
| 2137 | else if(szn<0) |
---|
| 2138 | { |
---|
| 2139 | G4int nn=-szn; |
---|
| 2140 | szn=0; |
---|
| 2141 | dn+=nn; |
---|
| 2142 | } |
---|
| 2143 | G4int sz =szn/1000; |
---|
| 2144 | G4int n =szn%1000; |
---|
| 2145 | if(n>700) |
---|
| 2146 | { |
---|
| 2147 | n-=1000; |
---|
| 2148 | dz--; |
---|
| 2149 | } |
---|
| 2150 | G4int z =sz%1000-dz; |
---|
| 2151 | if(z>700) |
---|
| 2152 | { |
---|
| 2153 | z-=1000; |
---|
| 2154 | ds--; |
---|
| 2155 | } |
---|
| 2156 | s =sz/1000-ds; |
---|
| 2157 | G4int b=z+n+s; |
---|
| 2158 | d=n+b; |
---|
| 2159 | u=z+b; |
---|
| 2160 | if (d<0&&u<0&&s<0) return G4QContent(0,0,0,-d,-u,-s); |
---|
| 2161 | else if (u<0&&s<0) return G4QContent(d,0,0,0,-u,-s); |
---|
| 2162 | else if (d<0&&s<0) return G4QContent(0,u,0,-d,0,-s); |
---|
| 2163 | else if (d<0&&u<0) return G4QContent(0,0,s,-d,-u,0); |
---|
| 2164 | else if (u<0) return G4QContent(d,0,s,0,-u,0); |
---|
| 2165 | else if (s<0) return G4QContent(d,u,0,0,0,-s); |
---|
| 2166 | else if (d<0) return G4QContent(0,u,s,-d,0,0); |
---|
| 2167 | else return G4QContent(d,u,s,0,0,0); |
---|
| 2168 | } |
---|
| 2169 | return G4QContent(0,0,0,0,0,0); |
---|
| 2170 | } |
---|
| 2171 | |
---|
| 2172 | // Quark Content of pseudo-meson for q_i->q_o Quark Exchange |
---|
| 2173 | G4QContent G4QPDGCode::GetExQContent(G4int i, G4int o) const |
---|
| 2174 | {// ================================================== |
---|
| 2175 | G4QContent cQC(0,0,0,0,0,0); |
---|
| 2176 | if (!i) cQC.IncAD(); |
---|
| 2177 | else if(i==1) cQC.IncAU(); |
---|
| 2178 | else if(i==2) cQC.IncAS(); |
---|
| 2179 | else G4cerr<<"***G4QPDGCode::GetExQContent: strange entering quark i="<<i<<G4endl; |
---|
| 2180 | if (!o) cQC.IncD(); |
---|
| 2181 | else if(o==1) cQC.IncU(); |
---|
| 2182 | else if(o==2) cQC.IncS(); |
---|
| 2183 | else G4cerr<<"***G4QPDGCode::GetExQContent: strange exiting quark o="<<o<<G4endl; |
---|
| 2184 | return cQC; |
---|
| 2185 | } |
---|
| 2186 | |
---|
| 2187 | |
---|
| 2188 | // Relative Cross Index for q_i->q_o (q=0(d),1(u),2(s), 7 means there is no such cluster) |
---|
| 2189 | G4int G4QPDGCode::GetRelCrossIndex(G4int i, G4int o) const |
---|
| 2190 | {// ===================================================== |
---|
| 2191 | // pD,nD,DD,DD,ppDnnDpDDnDD n, p, L,S-,S+,X-,X0,O- |
---|
| 2192 | static const G4int b01[16]={ 7,15, 7, 7, 7, 7, 7, 7, 1, 7, 7,-1, 7, 7, 7, 7}; |
---|
| 2193 | static const G4int b02[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 2, 7, 7, 7, 7, 7, 7, 7}; |
---|
| 2194 | static const G4int b10[16]={18, 7, 7, 7, 7, 7, 7, 7, 7,-1, 7, 7,-2, 7, 7, 7}; // Iso+Bary |
---|
| 2195 | static const G4int b12[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 1, 7, 7, 7, 7, 7, 7}; |
---|
| 2196 | static const G4int b20[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,-2, 7,-3, 7,-4, 7}; |
---|
| 2197 | static const G4int b21[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,-1,-3, 7,-3, 7, 7}; |
---|
| 2198 | // nn,np,pp,Ln,Lp,LL,nS,pS,nnpnppLnnLpnLppLLnLLpnnS |
---|
| 2199 | static const G4int d01[16]={ 1, 1, 7, 1, 7, 7,-3, 7, 1, 7, 1, 1, 7, 1, 7,-5}; |
---|
| 2200 | static const G4int d02[16]={ 3, 3, 7, 2, 7, 7, 7, 7, 3, 3, 3, 3, 7, 7, 7, 7}; |
---|
| 2201 | static const G4int d10[16]={ 7,-1,-1, 7,-1, 7, 7,-4, 7,-1, 7,-1,-1, 7,-1, 7}; //B=2,3 |
---|
| 2202 | static const G4int d12[16]={ 7, 2, 2, 7, 1, 7, 7, 7, 2, 2, 7, 2, 2, 7, 7, 7}; |
---|
| 2203 | static const G4int d20[16]={ 7, 7, 7,-3,-3,-2, 7,-5, 7, 7, 7,-3,-3,-3,-3, 7}; |
---|
| 2204 | static const G4int d21[16]={ 7, 7, 7,-2,-2,-1,-6, 7, 7, 7,-2,-2, 7,-2,-2, 7}; |
---|
| 2205 | // nn,np,pp,Ln,Lp,nn,np,pp! n, p,nn,np,pp,Ln,Lp |
---|
| 2206 | static const G4int m01[15]={ 1, 1, 7, 1, 7, 1, 1, 7, 1, 7, 1, 1, 7, 1, 7};// Multibaryons |
---|
| 2207 | static const G4int m02[15]={ 3, 3, 7, 3, 3, 7, 7, 7, 3, 3, 3, 3, 7, 7, 7};// @@ Regular ! |
---|
| 2208 | static const G4int m10[15]={ 7,-1,-1, 7,-1, 7,-1,-1, 7,-1, 7,-1,-1, 7,-1};// 01->1,10->-1 |
---|
| 2209 | static const G4int m12[15]={ 7, 2, 2, 2, 2, 7, 7, 7, 2, 2, 7, 2, 2, 7, 7};// 12->2,21->-2 |
---|
| 2210 | static const G4int m20[15]={ 7, 7, 7,-3,-3, 7,-3,-3, 7, 7, 7,-3,-3,-3,-3};// 02->3,20->-3 |
---|
| 2211 | static const G4int m21[15]={ 7, 7, 7,-2,-2,-2,-2, 7, 7, 7,-2,-2, 7,-2,-2}; |
---|
| 2212 | static const G4int fragmStart = 82; // "Isonuclei are added && Leptons are added" |
---|
| 2213 | |
---|
| 2214 | if(theQCode<fragmStart) return 7; |
---|
| 2215 | G4int sub=theQCode-fragmStart; |
---|
| 2216 | if ( (sub > 1 && sub < 8) || sub == 15) return 7; //@@Why they are in clusters?-Residuals(?) |
---|
| 2217 | G4int rel=sub; // case of nuclear baryons and isonuclei |
---|
| 2218 | if (sub>31)rel =(sub-32)%15; // case of heavy fragments (BaryNum>3) |
---|
| 2219 | else if(sub>15)rel = sub-16; // case of nuclear di-baryon & tri-baryons |
---|
| 2220 | #ifdef debug |
---|
| 2221 | G4cout<<"G4QPDGCode::RelGetCrossIndex:i="<<i<<",o="<<o<<",su="<<sub<<",re="<<rel<<G4endl; |
---|
| 2222 | #endif |
---|
| 2223 | if (!i) // ==> input quark = 0(d) (d=-1/3) |
---|
| 2224 | { |
---|
| 2225 | if (!o) return 0; // -> output quark = 0(d) => 0 = the same cluster |
---|
| 2226 | else if(o==1) // -> output quark = 1(u) (ch--) |
---|
| 2227 | { |
---|
| 2228 | if (sub<16) return b01[rel]; |
---|
| 2229 | else if(sub<32) return d01[rel]; |
---|
| 2230 | else return m01[rel]; |
---|
| 2231 | } |
---|
| 2232 | else if(o==2) |
---|
| 2233 | { |
---|
| 2234 | if (sub<16) return b02[rel]; |
---|
| 2235 | else if(sub<32) return d02[rel]; |
---|
| 2236 | else return m02[rel]; |
---|
| 2237 | } |
---|
| 2238 | } |
---|
| 2239 | else if(i==1) // ==> input quark = 1(u) (u=2/3) |
---|
| 2240 | { |
---|
| 2241 | if (!o) // -> output quark = 0(d) (ch++) |
---|
| 2242 | { |
---|
| 2243 | if (sub<16) return b10[rel]; |
---|
| 2244 | else if(sub<32) return d10[rel]; |
---|
| 2245 | else return m10[rel]; |
---|
| 2246 | } |
---|
| 2247 | else if(o==1) return 0; // -> output quark = 1(u) => 0 = the same cluster |
---|
| 2248 | else if(o==2) |
---|
| 2249 | { |
---|
| 2250 | if (sub<16) return b12[rel]; |
---|
| 2251 | else if(sub<32) return d12[rel]; |
---|
| 2252 | else return m12[rel]; |
---|
| 2253 | } |
---|
| 2254 | } |
---|
| 2255 | else if(i==2) |
---|
| 2256 | { |
---|
| 2257 | if (!o) |
---|
| 2258 | { |
---|
| 2259 | if (sub<16) return b20[rel]; |
---|
| 2260 | else if(sub<32) return d20[rel]; |
---|
| 2261 | else return m20[rel]; |
---|
| 2262 | } |
---|
| 2263 | else if(o==1) |
---|
| 2264 | { |
---|
| 2265 | if (sub<16) return b21[rel]; |
---|
| 2266 | else if(sub<32) return d21[rel]; |
---|
| 2267 | else return m21[rel]; |
---|
| 2268 | } |
---|
| 2269 | else if(o==2) return 0; |
---|
| 2270 | } |
---|
| 2271 | return 7; |
---|
| 2272 | } |
---|
| 2273 | |
---|
| 2274 | // Get number of Combinations for q_i->q_o |
---|
| 2275 | G4int G4QPDGCode::GetNumOfComb(G4int i, G4int o) const |
---|
| 2276 | {// ================================================ |
---|
| 2277 | if(i>-1&&i<3) |
---|
| 2278 | { |
---|
| 2279 | G4int shiftQ=GetRelCrossIndex(i, o); |
---|
| 2280 | G4int sQCode=theQCode; // QCode of the parent cluster |
---|
| 2281 | if (shiftQ==7) return 0; // A parent cluster doesn't exist |
---|
| 2282 | else if(!shiftQ) sQCode+=shiftQ; // Shift QCode of son to QCode of his parent |
---|
| 2283 | G4QPDGCode parent; // Create a temporary (fake) parent cluster |
---|
| 2284 | parent.InitByQCode(sQCode); // Initialize it by Q-Code |
---|
| 2285 | G4QContent parentQC=parent.GetQuarkContent(); // Quark Content of the parent cluster |
---|
| 2286 | if (!o) return parentQC.GetD(); |
---|
| 2287 | else if(o==1) return parentQC.GetU(); |
---|
| 2288 | else if(o==2) return parentQC.GetS(); |
---|
| 2289 | else G4cerr<<"***G4QPDGCode:::GetNumOfComb: strange exiting quark o="<<o<<G4endl; |
---|
| 2290 | } |
---|
| 2291 | else G4cerr<<"***G4QPDGCode:::GetNumOfComb: strange entering quark i="<<i<<G4endl; |
---|
| 2292 | return 0; |
---|
| 2293 | } |
---|
| 2294 | |
---|
| 2295 | // Get a total number of Combinations for q_i |
---|
| 2296 | G4int G4QPDGCode::GetTotNumOfComb(G4int i) const |
---|
| 2297 | {// ========================================== |
---|
| 2298 | G4int tot=0; |
---|
| 2299 | if(i>-1&&i<3) for(int j=0; j<3; j++) tot+=GetNumOfComb(i, j); |
---|
| 2300 | else G4cerr<<"***G4QPDGCode:::GetTotNumOfComb: strange entering quark i="<<i<<G4endl; |
---|
| 2301 | return tot; |
---|
| 2302 | } |
---|
| 2303 | |
---|
| 2304 | // Converts nuclear PDGCode to Z(#of protons), N(#of neutrons), S(#of lambdas) values |
---|
| 2305 | void G4QPDGCode::ConvertPDGToZNS(G4int nucPDG, G4int& z, G4int& n, G4int& s) |
---|
| 2306 | {// ======================================================================= |
---|
| 2307 | if(nucPDG>80000000&&nucPDG<100000000) // Condition of conversion |
---|
| 2308 | { |
---|
| 2309 | z=0; |
---|
| 2310 | n=0; |
---|
| 2311 | s=0; |
---|
| 2312 | G4int r=nucPDG; |
---|
| 2313 | if(r==90000000) return; |
---|
| 2314 | G4int cn =r%1000; // candidate to #of neutrons |
---|
| 2315 | if(cn) |
---|
| 2316 | { |
---|
| 2317 | if(cn>500) cn-=1000; // AntiNeutrons |
---|
| 2318 | n=cn; // Increment neutrons |
---|
| 2319 | r-=cn; // Subtract them from the residual |
---|
| 2320 | if(r==90000000) return; |
---|
| 2321 | } |
---|
| 2322 | G4int cz =r%1000000; // candidate to #of neutrons |
---|
| 2323 | if(cz) |
---|
| 2324 | { |
---|
| 2325 | if(cz>500000) cz-=1000000; // AntiProtons |
---|
| 2326 | z=cz/1000; // Number of protons |
---|
| 2327 | r-=cz; // Subtract them from the residual |
---|
| 2328 | if(r==90000000) return; |
---|
| 2329 | } |
---|
| 2330 | G4int cs =r%10000000; // candidate to #of neutrons |
---|
| 2331 | if(cs) |
---|
| 2332 | { |
---|
| 2333 | if(cs>5000000) cs-=10000000; // AntiLambda |
---|
| 2334 | s=cs/1000000; // Number of Lambdas |
---|
| 2335 | } |
---|
| 2336 | } |
---|
| 2337 | return; |
---|
| 2338 | } |
---|