G4QPDGCode.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 //      ---------------- G4QPDGCode ----------------
00030 //             by Mikhail Kossov, Sept 1999.
00031 //      class for Hadron definitions in CHIPS Model
00032 // -------------------------------------------------------------------
00033 // Short description: The PDG Code is made on the basis of the Quark
00034 // Content (G4QuarkContent) of the hadronic state (including nuclear
00035 // fragments). The PDG code of the ground state (e.g. pi, N, etc.) is
00036 // calculated. It includes a complicated algortithm of the G.S. mass
00037 // calculation for nuclear fragments (now it is synchronised with the
00038 // G4 nuclear massess).
00039 // -------------------------------------------------------------------
00040 
00041 //#define debug
00042 //#define pdebug
00043 //#define qdebug
00044 //#define idebug
00045 //#define sdebug
00046 
00047 #include <cmath>
00048 #include <cstdlib>
00049 
00050 #include "G4QPDGCodeVector.hh"
00051 
00052 using namespace std;
00053 
00054 G4QPDGCode::G4QPDGCode(G4int PDGCode): thePDGCode(PDGCode)
00055 {
00056 #ifdef sdebug
00057   G4cout<<"G4QPDGCode:Constructer is called with PDGCode="<<PDGCode<<G4endl;  
00058 #endif
00059   if(PDGCode==130) PDGCode= 311; // Safety. Should not happen.
00060   if(PDGCode==310) PDGCode=-311; // Safety. Should not happen.
00061   if(PDGCode==90000000)
00062   {
00063     thePDGCode=22;
00064     theQCode=6;
00065   }
00066   else if(PDGCode) theQCode=MakeQCode(PDGCode);
00067   else        
00068   {
00069 #ifdef sdebug
00070     G4cout<<"***G4QPDGCode: Constructed with PDGCode=0, QCode=-2"<<G4endl;  
00071 #endif
00072     theQCode=-2;
00073   }
00074 #ifdef debug
00075   if(PDGCode==3222)G4cout<<"G4QPDGCd:Con(PDG) PDG="<<PDGCode<<", QCode="<<theQCode<<G4endl;
00076 #endif
00077 }
00078 
00079 G4QPDGCode::G4QPDGCode(G4bool f, G4int QCode): theQCode(QCode)
00080 {
00081   if(f&&QCode<0)G4cerr<<"***G4QPDGCode::Constr. QCode="<<QCode<<G4endl;
00082   thePDGCode = MakePDGCode(QCode);
00083 #ifdef debug
00084   G4cout<<"G4QPDGCode::Constr: PDGCode="<<thePDGCode<<G4endl;
00085 #endif
00086 }
00087 
00088 G4QPDGCode::G4QPDGCode(G4QContent QCont)
00089 {
00090   InitByQCont(QCont);
00091 }
00092 
00093 G4QPDGCode::G4QPDGCode(const G4QPDGCode& rhs)
00094 {
00095   thePDGCode =rhs.thePDGCode;
00096   theQCode   =rhs.theQCode;
00097 }
00098 
00099 G4QPDGCode::G4QPDGCode(G4QPDGCode* rhs)
00100 {
00101   thePDGCode =rhs->thePDGCode;
00102   theQCode   =rhs->theQCode;
00103 }
00104 
00105 const G4QPDGCode& G4QPDGCode::operator=(const G4QPDGCode& rhs)
00106 {
00107   if(this != &rhs)                          // Beware of self assignment
00108   {
00109     thePDGCode =rhs.thePDGCode;
00110     theQCode   =rhs.theQCode;
00111   }
00112   return *this;
00113 }
00114 
00115 G4QPDGCode::~G4QPDGCode() {}
00116 
00117 // Standard output for QPDGCode
00118 ostream& operator<<(ostream& lhs, G4QPDGCode& rhs)
00119 {
00120   lhs << "[ PDG=" << rhs.GetPDGCode() << ", Q=" << rhs.GetQCode() << "]";
00121   return lhs;
00122 }
00123 
00124 // Standard output for const QPDGCode
00125 ostream& operator<<(ostream& lhs, const G4QPDGCode& rhs)
00126 {
00127   lhs << "[ PDG=" << rhs.GetPDGCode() << ", Q=" << rhs.GetQCode() << "]";
00128   return lhs;
00129 }
00130 
00131 // Overloading of QPDGCode addition
00132 G4int operator+(const G4QPDGCode& lhs, const G4QPDGCode& rhs)
00133 {
00134   G4int  s_value  = lhs.GetPDGCode();
00135   return s_value += rhs.GetPDGCode();
00136 }
00137 G4int operator+(const G4QPDGCode& lhs, const G4int& rhs)
00138 {
00139   G4int  s_value  = lhs.GetPDGCode();
00140   return s_value += rhs;
00141 }
00142 G4int operator+(const G4int& lhs, const G4QPDGCode& rhs)
00143 {
00144   G4int  s_value  = lhs;
00145   return s_value += rhs.GetPDGCode();
00146 }
00147 
00148 // Overloading of QPDGCode subtraction
00149 G4int operator-(const G4QPDGCode& lhs, const G4QPDGCode& rhs)
00150 {
00151   G4int  s_value  = lhs.GetPDGCode();
00152   return s_value -= rhs.GetPDGCode();
00153 }
00154 G4int operator-(const G4QPDGCode& lhs, const G4int& rhs)
00155 {
00156   G4int  s_value  = lhs.GetPDGCode();
00157   return s_value -= rhs;
00158 }
00159 G4int operator-(const G4int& lhs, const G4QPDGCode& rhs)
00160 {
00161   G4int  s_value  = lhs;
00162   return s_value -= rhs.GetPDGCode();
00163 }
00164 
00165 // Overloading of QPDGCode multiplication
00166 G4int operator*(const G4QPDGCode& lhs, const G4QPDGCode& rhs)
00167 {
00168   G4int  s_value  = lhs.GetPDGCode();
00169   return s_value *= rhs.GetPDGCode();
00170 }
00171 
00172 G4int operator*(const G4QPDGCode& lhs, const G4int& rhs)
00173 {
00174   G4int  s_value  = lhs.GetPDGCode();
00175   return s_value *= rhs;
00176 }
00177 
00178 G4int operator*(const G4int& lhs, const G4QPDGCode& rhs)
00179 {
00180   G4int  s_value  = lhs;
00181   return s_value *= rhs.GetPDGCode();
00182 }
00183 
00184 // Overloading of QPDGCode division
00185 G4int operator/(const G4QPDGCode& lhs, const G4QPDGCode& rhs)
00186 {
00187   G4int  s_value  = lhs.GetPDGCode();
00188   return s_value /= rhs.GetPDGCode();
00189 }
00190 
00191 G4int operator/(const G4QPDGCode& lhs, const G4int& rhs)
00192 {
00193   G4int  s_value  = lhs.GetPDGCode();
00194   return s_value /= rhs;
00195 }
00196 
00197 G4int operator/(const G4int& lhs, const G4QPDGCode& rhs)
00198 {
00199   G4int  s_value  = lhs;
00200   return s_value /= rhs.GetPDGCode();
00201 }
00202 
00203 // Overloading of QPDGCode residual
00204 G4int operator%(const G4QPDGCode& lhs, const G4int& rhs)
00205 {
00206   G4int  s_value  = lhs.GetPDGCode();
00207   return s_value %= rhs;
00208 }
00209 
00210 // TRUE if it is not RealNeutral (111,221,331 etc), FALSE if it is.
00211 G4bool G4QPDGCode::TestRealNeutral(const G4int& PDGCode)
00212 {
00213   if(PDGCode>0 && PDGCode<999)    // RealNeutral are always positive && mesons
00214   {
00215     if(PDGCode==22) return false; // Photon
00216     G4int p=PDGCode/10;
00217     if(p/10==p%10) return false; // This is a RealNeutral
00218   }
00219   return true;
00220 }
00221 
00222 // Make a Q Code out of the PDG Code
00223 G4int G4QPDGCode::MakePDGCode(const G4int& QCode)
00224 {
00225   //static const G4int modi = 81;  // Q Codes for more than di-baryon nuclei
00226   //static const G4int modi = 89;  // Q Codes for more than di-baryon nuclei "IsoNuclei"
00227   //static const G4int modi = 122; // Q Codes: more than quarta-baryon nuclei "Lept/Hyper"
00228   static const G4int modi = 85;    // Reduced Q Codes: > quarta-baryon nuclei "Lept/Hyper"
00229   //static G4int qC[modi]={ 11,   12,   13,   14,   15,   16,   22,   23,   24,   25, // 10
00230   //                        37,  110,  220,  330,  111,  211,  221,  311,  321,  331, // 20
00231   //                      2112, 2212, 3122, 3112, 3212, 3222, 3312, 3322,  113,  213, // 30
00232   //                       223,  313,  323,  333, 1114, 2114, 2214, 2224, 3124, 3114, // 40
00233   //                      3214, 3224, 3314, 3324, 3334,  115,  215,  225,  315,  325, // 50
00234   //                       335, 2116, 2216, 3126, 3116, 3216, 3226, 3316, 3326,  117, // 60
00235   //                       217,  227,  317,  327,  337, 1118, 2118, 2218, 2228, 3128, // 70
00236   //                      3118, 3218, 3228, 3318, 3328, 3338,  119,  219,  229,  319, // 80
00237   //                       329,  339, 90002999 , 89999003 , 90003998 ,
00238   //                       89998004 , 90003999 , 89999004 , 90004998 , 89998005 ,     // 90
00239   //                       90000001 , 90001000 , 91000000 , 90999001 , 91000999 ,
00240   //                       91999000 , 91999999 , 92999000 , 90000002 , 90001001 ,     //100
00241   //                       90002000 , 91000001 , 91001000 , 92000000 , 90999002 ,
00242   //                       91001999 , 90001002 , 90002001 , 91000002 , 91001001 ,     //110
00243   //                       91002000 , 92000001 , 92001000 , 90999003 , 90001003 ,
00244   //                       90002002 , 90003001 , 91001002 , 91002001 , 92000002 ,     //120
00245   //                       92001001 , 92002000};                                      //122
00246   static G4int qC[modi] ={  11,   12,   13,   14,   15,   16,   22,   23,   24,   25, // 10
00247                             37,  110,  220,  330,  111,  211,  221,  311,  321,  331, // 20
00248                           2112, 2212, 3122, 3112, 3212, 3222, 3312, 3322,  113,  213, // 30
00249                            223,  313,  323,  333, 1114, 2114, 2214, 2224, 3124, 3114, // 40
00250                                   3214, 3224, 3314, 3324, 3334,                               // 45
00251                                     90002999 , 89999003 , 90003998 , 89998004 , 90003999 ,    // 50
00252                                     89999004 , 90004998 , 89998005 , 90000001 , 90001000 ,    // 55
00253                                     91000000 , 90999001 , 91000999 , 91999000 , 91999999 ,    // 60
00254                                     92999000 , 90000002 , 90001001 , 90002000 , 91000001 ,    // 65
00255                                     91001000 , 92000000 , 90999002 , 91001999 , 90001002 ,    // 70
00256                                     90002001 , 91000002 , 91001001 , 91002000 , 92000001 ,    // 75
00257                                     92001000 , 90999003 , 90001003 , 90002002 , 90003001 ,    // 80
00258                                     91001002 , 91002001 , 92000002 , 92001001 , 92002000};    // 85
00259   static G4int aC[15] = {1,1000,999001,1000000,1000999,1999000,1999999,        // sum 1
00260                          2,1001,2000,1000001,1001000,1999001,2000000,2000999}; // sum 2
00261   if      (QCode<0)
00262   {
00263     G4cerr<<"***G4QPDGCode::MakePDGCode: negative Q Code ="<<QCode<<G4endl;
00264     return 0;
00265   }
00266   else if (QCode>=modi)
00267   {
00268     //G4int q=QCode-modi;              // Starting BarNum=3
00269     //G4int a=q/15+1;                  // BarNum/2
00270     //G4int b=q%15;
00271     G4int q=QCode-modi;              // Starting BarNum=5
00272     G4int a=q/15+2;                  // BarNum/2
00273     G4int b=q%15;
00274     return 90000000+a*1001+aC[b];
00275   }
00276   return qC[theQCode];
00277 }
00278 
00279 // Hadronic masses synhronized with the Geant4 hadronic masses
00280 G4double G4QPDGCode:: QHaM(G4int nQ)
00281 {
00282   static G4bool iniFlag=true;
00283   //static G4double mass[nQHM]={.511,0.,105.65837, 0., 1777., 0.,   0., 91188., 80403., 140.00
00284   //,120.000,    800.,     980.,    1370.,  134.98,  139.57, 547.51, 497.65, 493.68, 957.78
00285   //,939.5654,938.272, 1115.683,  1197.45, 1192.64, 1189.37,1321.31,1314.83,  775.5,  775.5
00286   //, 782.65,   896.0,   891.66,  1019.46,   1232.,   1232.,  1232.,  1232., 1519.5, 1387.2
00287   //, 1383.7,  1382.8,    1535.,   1531.8, 1672.45,  1318.3, 1318.3, 1275.4, 1432.4, 1425.6
00288   //,  1525.,   1680.,    1680.,    1820.,   1915.,   1915.,  1915.,  2025.,  2025.,  1691.
00289   //,  1691.,   1667.,    1776.,    1776.,   1854.,   1950.,  1950.,  1950.,  1950.,  2100.
00290   //,  2030.,   2030.,    2030.,    2127.,   2127.,   2252.,  2020.,  2020.,  2044.,  2045.
00291   //, 2045., 2297., 2170.272, 2171.565, 2464., 2464., 3108.544, 3111.13,3402.272,3403.565};
00292   // Reduced:
00293   static G4double mass[nQHM]={.511, 0., 105.65837, 0., 1777., 0.,   0., 91188., 80403., 140.00
00294     ,120.000,    800.,     980.,    1370.,  134.98,  139.57, 547.51, 497.65, 493.68, 957.78
00295     ,939.5654,938.272, 1115.683,  1197.45, 1192.64, 1189.37,1321.31,1314.83,  775.5,  775.5
00296     , 782.65,   896.0,   891.66,  1019.46,   1232.,   1232.,  1232.,  1232., 1519.5, 1387.2
00297     , 1383.7,  1382.8,    1535.,   1531.8, 1672.45,2170.272,2171.565, 2464., 2464.,3108.544
00298     ,3111.13,3402.272, 3403.565};
00299   if(iniFlag) // Initialization of the Geant4 hadronic masses
00300   {
00301     mass[ 0]=      G4Electron::Electron()->GetPDGMass();
00302     mass[ 1]=    G4NeutrinoE::NeutrinoE()->GetPDGMass();
00303     mass[ 2]=    G4MuonMinus::MuonMinus()->GetPDGMass();
00304     mass[ 3]=  G4NeutrinoMu::NeutrinoMu()->GetPDGMass();
00305     mass[ 4]=      G4TauMinus::TauMinus()->GetPDGMass();
00306     mass[ 5]=G4NeutrinoTau::NeutrinoTau()->GetPDGMass();
00307     mass[14]=      G4PionZero::PionZero()->GetPDGMass();
00308     mass[15]=    G4PionMinus::PionMinus()->GetPDGMass();
00309     mass[16]=                G4Eta::Eta()->GetPDGMass();
00310     mass[17]=      G4KaonZero::KaonZero()->GetPDGMass();
00311     mass[18]=    G4KaonMinus::KaonMinus()->GetPDGMass();
00312     mass[19]=      G4EtaPrime::EtaPrime()->GetPDGMass();
00313     mass[20]=        G4Neutron::Neutron()->GetPDGMass();
00314     mass[21]=          G4Proton::Proton()->GetPDGMass();
00315     mass[22]=          G4Lambda::Lambda()->GetPDGMass();
00316     mass[23]=  G4SigmaMinus::SigmaMinus()->GetPDGMass();
00317     mass[24]=    G4SigmaZero::SigmaZero()->GetPDGMass();
00318     mass[25]=    G4SigmaPlus::SigmaPlus()->GetPDGMass();
00319     mass[26]=        G4XiMinus::XiMinus()->GetPDGMass();
00320     mass[27]=          G4XiZero::XiZero()->GetPDGMass();
00321     mass[44]=  G4OmegaMinus::OmegaMinus()->GetPDGMass();
00322     iniFlag=false;
00323   }
00324   if(nQ<0 || nQ>=nQHM)
00325   {
00326     G4cout<<"***G4QPDGCode::QHaM: negative Q-code or Q="<<nQ<<" >= nQmax = "<<nQHM<<G4endl;
00327     return 0.;
00328   }
00329   return mass[nQ];
00330 }
00331 
00332 // Make a Q Code out of the PDG Code
00333 G4int G4QPDGCode::MakeQCode(const G4int& PDGCode)
00334 {
00335   static const G4int qr[10]={0,13,19,27,33,44,50,58,64,75};
00336   G4int PDGC=abs(PDGCode);        // Qcode is always not negative
00337   G4int s_value=0;
00338   G4int z=0;
00339   G4int n=0;
00340   if (PDGC>100000000)             // Not supported
00341   {
00342 #ifdef debug
00343     G4cout<<"***G4QPDGCode::MakeQCode: Unknown in Q-System code: "<<PDGCode<<G4endl;
00344 #endif
00345     return -2;
00346   }
00347   else if (PDGC>80000000 && PDGC<100000000) // Try to convert the NUCCoding to PDGCoding
00348   {
00349     //if(PDGC==90000000) return 6;            // @@ already done in the constructor
00350     ConvertPDGToZNS(PDGC, z, n, s_value);
00351     G4int b=n+z+s_value;                           // Baryon number
00352 #ifdef debug
00353     G4cout<<"***G4QPDGCode::Z="<<z<<",N="<<n<<",S="<<s_value<<G4endl;
00354 #endif
00355     if(b<0)                                        // ---> Baryons & Fragments
00356     {
00357       b=-b;
00358       n=-n;
00359       z=-z;
00360       s_value=-s_value;
00361       PDGC=90000000+s_value*1000000+z*1000+n;      // New PDGC for anti-baryons
00362     }
00363     else if(!b)                                    // --> Mesons
00364     {
00365       //G4bool anti=false;                           // For the PDG conversion
00366       if(z<0)                                      // --> Mesons conversion
00367       {
00368         n=-n;
00369         z=-z;
00370         s_value=-s_value;
00371         //anti=true;                                 // For the PDG conversion
00372       }
00373       if(!z)
00374       {
00375         if(s_value>0)
00376         {
00377           n=-n;
00378           s_value=-s_value;
00379           //anti=true;                               // For the PDG conversion
00380         }
00381         if     (s_value==-1) return 17;            // K0
00382         else if(s_value==-2) return -1;            // K0+K0 chipolino
00383         else           return -2;                  // Not supported by Q Code
00384       }
00385       else                                         // --> z>0
00386       {
00387         if(z==1)
00388         {
00389           if   (s_value==-1) return 18;            // K+
00390           else         return 15;                  // pi+
00391         }
00392         else if(z==2)  return -1;                  // Chipolino
00393         else           return -2;                  // Not supported by Q Code
00394       }
00395     } // End of meson case
00396     if(b>0)                                        // --> Baryoniums case
00397     {
00398       if(b==1)                                     // --> Baryons+Hyperons
00399       {
00400         if(PDGC>80000000)
00401         {
00402           if(!s_value)                               // --> Baryons
00403           {
00404             if     (!z)  return 53;                  // neutron
00405             else if(z==1)return 54;                  // proton
00406             else         return -2;                  // Not supported by Q Code
00407           }
00408           else if(s_value==1)                        // --> Hyperons
00409           {
00410             if(z==-1)    return 56;                  // Sigma-
00411             else if(!z)  return 55;                  // Lambda
00412             else if(z==1)return 57;                  // Sigma+
00413             else         return -2;                  // Not supported by Q Code
00414           }
00415           else if(s_value==2)                        // --> Xi Hyperons
00416           {
00417             if(z==-1)    return 58;                  // Xi-
00418             else if(!z)  return 59;                  // Xi0
00419             else         return -2;                  // Not supported by Q Code
00420           }
00421           else if(s_value==3)                        // --> Xi Hyperons
00422           {
00423             if(z==-1)    return 60;                  // Omega-
00424             else         return -2;                  // Not supported by Q Code
00425           }
00426         }
00427         else
00428         {
00429           if(!s_value)                               // --> Baryons
00430           {
00431             if(z==-1)    return 34;                  // Delta-
00432             else if(!z)  return 20;                  // neutron
00433             else if(z==1)return 21;                  // proton
00434             else if(z==2)return 37;                  // Delta++
00435             else if(z==3||z==-2)return -1;           // Delta+pi Chipolino
00436             else         return -2;                  // Not supported by Q Code
00437           }
00438           else if(s_value==1)                        // --> Hyperons
00439           {
00440             if(z==-1)    return 23;                  // Sigma-
00441             else if(!z)  return 22;                  // Lambda (@@ 24->Sigma0)
00442             else if(z==1)return 25;                  // Sigma+
00443             else if(z==2||z==-2) return -1;          // Sigma+pi Chipolino
00444             else         return -2;                  // Not supported by Q Code
00445           }
00446           else if(s_value==2)                        // --> Xi Hyperons
00447           {
00448             if(z==-1)    return 26;                  // Xi-
00449             else if(!z)  return 27;                  // Xi0
00450             else if(z==1||z==-2)return -1;           // Xi+pi Chipolino
00451             else         return -2;                  // Not supported by Q Code
00452           }
00453           else if(s_value==3)                        // --> Xi Hyperons
00454           {
00455             if(z==-1)    return 44;                  // Omega-
00456             else if(!z||z==-2)  return -1;           // Omega+pi Chipolino
00457             else         return -2;                  // Not supported by Q Code
00458           }
00459         }
00460       }
00461       else
00462       {
00463         if(b==2)
00464         {
00465           if     (PDGC==90002999) return 45;       // p DEL++
00466           else if(PDGC==89999003) return 46;       // n DEL-
00467           else if(PDGC==90003998) return 47;       // DEL++ DEL++
00468           else if(PDGC==89998004) return 48;       // DEL-  DEL-
00469           else if(PDGC==90999002) return 67;      // n Sigma-
00470           else if(PDGC==91001999) return 68;      // p Sigma+
00471         }
00472         if(b==3)
00473         {
00474           if     (PDGC==90003999) return 49;       // p p DEL++
00475           else if(PDGC==89999004) return 50;       // n n DEL-
00476           else if(PDGC==90004998) return 51;       // p DEL++ DEL++
00477           else if(PDGC==89998005) return 52;       // n DEL-  DEL-
00478           else if(PDGC==90999003) return 76;      // n n Sigma-
00479         }
00480       }
00481     }
00482   }
00483   if (PDGC<80000000)                // ----> Direct Baryons & Mesons
00484   {
00485     if     (PDGC<100)               // => Leptons and field bosons
00486     {
00487       if     (PDGC==10)  return -1; // Chipolino
00488       else if(PDGC==11)  return  0; // e-
00489       else if(PDGC==12)  return  1; // nu_e
00490       else if(PDGC==13)  return  2; // mu-
00491       else if(PDGC==14)  return  3; // nu_mu
00492       else if(PDGC==15)  return  4; // tau-
00493       else if(PDGC==16)  return  5; // nu_tau
00494       else if(PDGC==22)  return  6; // Photon
00495       else if(PDGC==23)  return  7; // Z0 boson
00496       else if(PDGC==24)  return  8; // W- boson
00497       else if(PDGC==25)  return  9; // H0 (neutral Higs boson)
00498       else if(PDGC==37)  return 10; // H- (charged Higs boson)
00499     }
00500     G4int r=PDGC%10;                // 2s+1
00501     G4int         Q= 0;
00502     if     (!r)
00503     {
00504       // Internal CHIPS codes for the wide f_0 states must be 9000221, 9010221, 10221
00505       if     (PDGC==110) return 11; // Low R-P: Sigma (pi,pi S-wave)
00506       else if(PDGC==220) return 12; // Midle Regeon-Pomeron
00507       else if(PDGC==330) return 13; // High Regeon-Pomeron
00508 #ifdef debug
00509       G4cout<<"***G4QPDGCode::MakeQCode: (0) Unknown in Q-System code: "<<PDGCode<<G4endl;
00510 #endif
00511       return -2;
00512     }
00513     else Q=qr[r];
00514     G4int p=PDGC/10;                // Quark Content
00515     if(r%2)                         // (2s+1 is odd) Mesons
00516     {
00517       if     (p==11) return Q+=1;
00518       else if(p==21) return Q+=2;
00519       else if(p==22) return Q+=3;
00520       else if(p==31) return Q+=4;
00521       else if(p==32) return Q+=5;
00522       else if(p==33) return Q+=6;
00523       else
00524       {
00525 #ifdef debug
00526         G4cout<<"*Warning*G4QPDGCode::MakeQCode:(1)UnknownQCode for PDG="<<PDGCode<<G4endl;
00527 #endif
00528         return -2;
00529       }
00530     }
00531     else                    // (2s+1 is even) Baryons
00532     {
00533       s_value=r/2;
00534       if(s_value%2)         // ((2s+1)/2 is odd) N Family
00535       {
00536         if     (p==211) return Q+=1;
00537         else if(p==221) return Q+=2;
00538         else if(p==312) return Q+=3;
00539         else if(p==311) return Q+=4;
00540         else if(p==321) return Q+=5;
00541         else if(p==322) return Q+=6;
00542         else if(p==331) return Q+=7;
00543         else if(p==332) return Q+=8;
00544         else
00545         {
00546 #ifdef debug
00547           G4cout<<"*Warning*G4QPDGCode::MakeQCode:(2) UnknownQCode, PDG="<<PDGCode<<G4endl;
00548 #endif
00549           return -2;
00550         }
00551       }
00552       else                  // ((2s+1)/2 is odd) Delta Family
00553       {
00554         if     (p==111) return Q+= 1;
00555         else if(p==211) return Q+= 2;
00556         else if(p==221) return Q+= 3;
00557         else if(p==222) return Q+= 4;
00558         else if(p==312) return Q+= 5;
00559         else if(p==311) return Q+= 6;
00560         else if(p==321) return Q+= 7;
00561         else if(p==322) return Q+= 8;
00562         else if(p==331) return Q+= 9;
00563         else if(p==332) return Q+=10;
00564         else if(p==333) return Q+=11;
00565         else
00566         {
00567 #ifdef debug
00568           G4cout<<"**G4QPDGCode::MakeQCode:(3) Unknown in Q-System code:"<<PDGCode<<G4endl;
00569 #endif
00570           return -2;
00571         }
00572       }
00573     }
00574   }
00575   else                        // Nuclear Fragments
00576   {
00577     G4int d=n+n+z+s_value;    // a#of d quarks
00578     G4int u=n+z+z+s_value;    // a#of u quarks
00579     G4int t=d+u+s_value;      // tot#of quarks
00580     if(t%3 || abs(t)<3)       // b=0 are in mesons
00581     {
00582 #ifdef debug
00583       G4cout<<"***G4QPDGCode::MakeQCode: Unknown PDGCode="<<PDGCode<<", t="<<t<<G4endl;
00584 #endif
00585       return -2;
00586     }
00587     else
00588     {
00589       G4int b=t/3;            // baryon number
00590       if(b==1)                // baryons
00591       {
00592         if     (s_value==0&&u==1&&d==2) return 53; // n
00593         else if(s_value==0&&u==2&&d==1) return 54; // p
00594         else if(s_value==1&&u==1&&d==1) return 55; // Lambda
00595         else if(s_value==1&&u==0&&d==2) return 56; // Sigma-
00596         else if(s_value==1&&u==2&&d==0) return 57; // Sigma+
00597         else if(s_value==2&&u==0&&d==1) return 58; // Xi-
00598         else if(s_value==2&&u==1&&d==0) return 59; // Xi0
00599         else if(s_value==3&&u==0&&d==0) return 60; // Omega-
00600         else
00601         {
00602 #ifdef debug
00603           G4cout<<"**G4QPDGCode::MakeQCode:(5) Unknown in Q-System code:"<<PDGCode<<G4endl;
00604 #endif
00605           return -2;
00606         }
00607       }
00608       else if(b==2)           // di-baryons
00609       {
00610         if     (s_value==0&&u==2&&d==4) return 61; // nn
00611         else if(s_value==0&&u==3&&d==3) return 62; // np
00612         else if(s_value==0&&u==4&&d==2) return 63; // pp
00613         else if(s_value==1&&u==2&&d==3) return 64; // nLambda
00614         else if(s_value==1&&u==3&&d==2) return 65; // pLambda
00615         else if(s_value==2&&u==2&&d==2) return 66; // LambdaLambda
00616         else
00617         {
00618 #ifdef debug
00619           G4cout<<"**G4QPDGCode::MakeQCode:(6) Unknown in Q-System code:"<<PDGCode<<G4endl;
00620 #endif
00621           return -2;
00622         }
00623       }
00624       else if(b==3)           // tri-baryons
00625       {
00626         if     (s_value==0&&u==4&&d==5) return 69; // pnn
00627         else if(s_value==0&&u==5&&d==4) return 70; // npp
00628         else if(s_value==1&&u==3&&d==5) return 71; // Lnn
00629         else if(s_value==1&&u==4&&d==4) return 72; // Lnp
00630         else if(s_value==1&&u==5&&d==3) return 73; // Lpp
00631         else if(s_value==2&&u==3&&d==4) return 74; // LLn
00632         else if(s_value==2&&u==4&&d==3) return 75; // LLp
00633         else if(s_value==1&&u==2&&d==6) return 76; // nnSigma-
00634         else
00635         {
00636 #ifdef debug
00637           G4cout<<"**G4QPDGCode::MakeQCode:(7) Unknown in Q-System code:"<<PDGCode<<G4endl;
00638 #endif
00639           return -2;
00640         }
00641       }
00642       G4int c=b/2;              // From here b>3: (4,5):c=2,g=0,1; (6,7):c=3,g=0,1; ...
00643       G4int g_value=b%2;
00644       G4int h=c*3;
00645       //G4int Q=57+c*15;
00646       //G4int Q=65+c*15;           // "IsoNuclei"
00647       //G4int Q=83+c*15;           // "Leptons/Hyperons"
00648       G4int Q=46+c*15;           // Reduced "Leptons/Hyperons"
00649       u-=h;
00650       d-=h;
00651       if(g_value)
00652       {
00653         if     (s_value==0&&u==1&&d==2) return Q+= 9;
00654         else if(s_value==0&&u==2&&d==1) return Q+=10;
00655         else if(s_value==1&&u==0&&d==2) return Q+=11;
00656         else if(s_value==1&&u==1&&d==1) return Q+=12;
00657         else if(s_value==1&&u==2&&d==0) return Q+=13;
00658         else if(s_value==2&&u==0&&d==1) return Q+=14;
00659         else if(s_value==2&&u==1&&d==0) return Q+=15;
00660         else
00661         {
00662 #ifdef debug
00663           G4cout<<"**G4QPDGCode::MakeQCode:(8) Unknown in Q-System code:"<<PDGCode<<G4endl;
00664 #endif
00665           return -2;
00666         }
00667       }
00668       else
00669       {
00670         if     (s_value==0&&u==-1&&d== 1) return Q+=1;
00671         else if(s_value==0&&u== 0&&d== 0) return Q+=2;
00672         else if(s_value==0&&u== 1&&d==-1) return Q+=3;
00673         else if(s_value==1&&u==-1&&d== 0) return Q+=4;
00674         else if(s_value==1&&u== 0&&d==-1) return Q+=5;
00675         else if(s_value==2&&u==-2&&d== 0) return Q+=6;
00676         else if(s_value==2&&u==-1&&d==-1) return Q+=7;
00677         else if(s_value==2&&u== 0&&d==-2) return Q+=8;
00678         else
00679         {
00680 #ifdef debug
00681           G4cout<<"**G4QPDGCode::MakeQCode:(9) Unknown in Q-System code:"<<PDGCode<<G4endl;
00682 #endif
00683           return -2;
00684         }
00685       }
00686     }
00687   }
00688   G4cout<<"*Warning*G4QPDGCode::MakeQCode:() Unknown Q Code for PDG = "<<PDGCode<<G4endl;
00689   return -2; //not reachable statement (fake, only for compiler)  
00690 }
00691 
00692 // Get the mean mass value for the PDG
00693 G4double G4QPDGCode::GetMass()
00694 {
00695   G4int ab=theQCode;
00696 #ifdef pdebug
00697   G4bool pPrint = thePDGCode == 3222 || ab == 25;
00698   if(pPrint)
00699   G4cout<<"G4QPDGCode::GetMass: Mass for Q="<<ab<<",PDG="<<thePDGCode<<",N="<<nQHM<<G4endl;
00700 #endif
00701   if ( (ab < 0 && thePDGCode < 80000000) || !thePDGCode) {
00702 #ifdef pdebug
00703     if(thePDGCode!=10 && pPrint)
00704       G4cout<<"**G4QPDGCode::GetMass:m=100000.,QC="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
00705 #endif
00706     return 100000.;
00707   }
00708   else if(ab>-1 && ab<nQHM)
00709   {
00710     G4double mass = QHaM(ab);
00711 #ifdef pdebug
00712     //if(pPrint)
00713     if(thePDGCode == 3222 || ab == 25)
00714     G4cout<<"G4QPDGCode::GetMa:m="<<mass<<",Q="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
00715 #endif
00716     return mass;                                // Get mass from the table
00717   }
00718   //if(szn==0) return m[15];                    // @@ mPi0   @@ MK ?
00719   if(thePDGCode==90000000)
00720   {
00721 #ifdef pdebug
00722     if(pPrint)
00723     G4cout<<"G4QPDGCode::GetMass:***m=0, QC="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
00724 #endif
00725     return 0.;
00726   }
00727   G4int s_value=0;
00728   G4int z=0;
00729   G4int n=0;
00730   ConvertPDGToZNS(thePDGCode, z, n, s_value);
00731   G4double m_value=GetNuclMass(z,n,s_value);
00732 #ifdef pdebug
00733   if(pPrint)
00734   G4cout<<"G4QPDG::GetM:PDG="<<thePDGCode<<"=>Z="<<z<<",N="<<n<<",S="<<s_value<<",M="<<m_value<<G4endl;
00735 #endif
00736   return m_value;
00737 }
00738 
00739 // Get the width value for the PDG
00740 G4double G4QPDGCode::GetWidth()
00741 {
00742   //static const int nW = 72;
00743   //static const int nW = 80; // "Isobars"
00744   //static G4double width[nQHM] = {0.,0.,0.,0.,0.,0.,0.,2.495,2.118,10.
00745   //  ,  10., 800.,  75., 350.,   0.,   0., .00118,  0.,   0., .203
00746   //  ,   0.,   0.,   0.,   0.,   0.,   0.,   0.,    0., 160., 160.
00747   //  , 8.41, 50.5, 50.8, 4.43, 120., 120., 120.,  120., 15.6,  39.
00748   //  ,  36., 35.8,  10.,   9.,   0., 107., 107., 185.5, 109., 98.5
00749   //  ,  76., 130., 130.,  80., 120., 120., 120.,   20.,  20., 160.
00750   //  , 160., 168., 159., 159.,  87., 300., 300.,  300., 300., 200.
00751   //  , 180., 180., 180.,  99.,  99.,  55., 387.,  387., 208., 198.
00752   //  , 198., 149., 120., 120., 170., 170., 120.,  120., 170., 170.};
00753   // Reduced:
00754   static G4double width[nQHM] = {0.,0.,0.,0.,0.,0.,0.,2.495,2.118,10.
00755     ,  10., 800.,  75., 350.,   0.,   0., .00118,  0.,   0., .203
00756     ,   0.,   0.,   0.,   0.,   0.,   0.,   0.,    0., 160., 160.
00757     , 8.41, 50.5, 50.8, 4.43, 120., 120., 120.,  120., 15.6,  39.
00758     ,  36., 35.8,  10.,   9.,   0., 120., 120.,  170., 170., 120.
00759     , 120., 170., 170.};
00760   G4int ab=abs(theQCode);
00761   if(ab<nQHM) return width[ab];
00762   return 0.;             // @@ May be real width should be implemented for nuclei e.g. pp
00763 }
00764 
00765 // Get a nuclear mass for Z (a#of protons), N (a#of neutrons), & S (a#of lambdas) 
00766 G4double G4QPDGCode::GetNuclMass(G4int z, G4int n, G4int s_value)
00767 {
00768   static const G4double anb = .01; // Antibinding for Ksi-n,Sig-n,Sig+p,Sig-nn,
00769   static const G4double mNeut= QHaM(20);
00770   static const G4double mProt= QHaM(21);
00771   static const G4double mLamb= QHaM(22);
00772   static const G4double mPiC = QHaM(15);
00773   static const G4double mKZ  = QHaM(17);
00774   static const G4double mKM  = QHaM(18);
00775   static const G4double mSiM = QHaM(23);
00776   static const G4double mSiP = QHaM(25);
00777   static const G4double mKsZ = QHaM(27);
00778   static const G4double mKsM = QHaM(26);
00779   static const G4double mOmM = QHaM(44);
00780   static const G4double mKZa = mKZ +anb;
00781   static const G4double mKMa = mKM +anb;
00782   static const G4double mSigM= mSiM+anb;
00783   static const G4double mSigP= mSiP+anb;
00784   static const G4double mKsiZ= mKsZ+anb;
00785   static const G4double mKsiM= mKsM+anb;
00786   static const G4double mOmeg= mOmM+anb;
00787   static const G4double mDiPi= mPiC+mPiC+anb;
00788   static const G4double mDiKZ= mKZa+mKZ;
00789   static const G4double mDiKM= mKMa+mKM;
00790   static const G4double mDiPr= mProt+mProt;
00791   static const G4double mDiNt= mNeut+mNeut;
00792   static const G4double mSmPi= mSiM+mDiPi;
00793   static const G4double mSpPi= mSiP+mDiPi;
00794   static const G4double mOmN = mOmeg+mNeut;
00795   static const G4double mSpP = mSigP+mProt;
00796   static const G4double mSpPP= mSpP +mProt;
00797   static const G4double mSmN = mSigM+mNeut;
00798   static const G4double mSmNN= mSmN +mNeut;
00799   // -------------- DAM Arrays ----------------------
00800   static const G4int iNR=76;    // Neutron maximum range for each Z
00801   static const G4int nEl = 105; // Maximum Z of the associative memory is "nEl-1=104"
00802   static const G4int iNF[nEl]={0,0,0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, // 14
00803                          1  ,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, // 29
00804                          16 , 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, // 44
00805                          31 , 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 53, 54, 55, // 59
00806                          56 , 56, 57, 57, 58, 60, 61, 63, 64, 65, 66, 67, 68, 69, 70, // 74
00807                          71 , 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, // 89
00808                          86 , 87, 88, 89, 91, 94, 98,103,109,115,122,128,134,140,146};//104
00809 #ifdef qdebug
00810   static G4int iNmin[nEl]={0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, // 14
00811                          1  ,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, // 29
00812                          16 , 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, // 44
00813                          31 , 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 53, 54, 55, // 59
00814                          56 , 56, 57, 57, 58, 60, 61, 63, 64, 65, 66, 67, 68, 69, 70, // 74
00815                          71 , 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, // 89
00816                          86 , 87, 88, 89, 91, 94, 98,103,109,115,122,128,134,140,146};//104
00817   static G4int iNmax=iNR;
00818   static G4int iNran[nEl]={19,20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, // 14
00819                          34 , 35, 36, 37, 38, 39, 40, 48, 48, 48, 48, 50, 50, 50, 52, // 29
00820                          53 , 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, // 44
00821                          68 , 69, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, // 59
00822                          73 , 73, 73, 73, 74, 74, 74, 74, 74, 74, 74, 74, 74, 75, 76, // 74
00823                          76 , 76, 76, 76, 76, 75, 74, 73, 72, 71, 70, 70, 69, 69, 69, // 89
00824                          68 , 68, 68, 67, 63, 59, 55, 51, 47, 43, 39, 35, 31, 27, 23};//104
00825 #endif
00826   static const G4int iNL[nEl]={19,20,21,22,23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, // 14
00827                          34 , 35, 36, 37, 38, 39, 40, 48, 48, 48, 48, 50, 50, 50, 52, // 29
00828                          53 , 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, // 44
00829                          68 , 69, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, // 59
00830                          73 , 73, 73, 73, 74, 74, 74, 74, 74, 74, 74, 74, 74, 75, 76, // 74
00831                          76 , 76, 76, 76, 76, 75, 74, 73, 72, 71, 70, 70, 69, 69, 69, // 89
00832                          68 , 68, 68, 67, 63, 59, 55, 51, 47, 43, 39, 35, 31, 27, 23};//104
00833    // ********* S=-4 vectors *************
00834   static G4bool iNin6[nEl]={false,false,false,false,false,false,false,
00835     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00836     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00837     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00838     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00839     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00840     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00841     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00842   static G4double VZ6[nEl][iNR];
00843   //********* S=-3 vectors *************
00844   static G4bool iNin7[nEl]={false,false,false,false,false,false,false,
00845     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00846     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00847     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00848     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00849     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00850     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00851     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00852   static G4double VZ7[nEl][iNR];
00853   // ********* S=-2 vectors *************
00854   static G4bool iNin8[nEl]={false,false,false,false,false,false,false,
00855     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00856     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00857     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00858     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00859     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00860     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00861     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00862   static G4double VZ8[nEl][iNR];
00863   // ********* S=-1 vectors *************
00864   static G4bool iNin9[nEl]={false,false,false,false,false,false,false,
00865     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00866     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00867     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00868     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00869     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00870     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00871     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00872   static G4double VZ9[nEl][iNR];
00873   // ********* S=0 vectors *************
00874   static G4bool iNin0[nEl]={false,false,false,false,false,false,false,
00875     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00876     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00877     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00878     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00879     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00880     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00881     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00882  static G4double VZ0[nEl][iNR];
00883   // ********* S=1 vectors *************
00884   static G4bool iNin1[nEl]={false,false,false,false,false,false,false,
00885     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00886     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00887     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00888     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00889     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00890     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00891     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00892   static G4double VZ1[nEl][iNR];
00893   // ********* S=2 vectors *************
00894   static G4bool iNin2[nEl]={false,false,false,false,false,false,false,
00895     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00896     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00897     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00898     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00899     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00900     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00901     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00902   static G4double VZ2[nEl][iNR];
00903   // ********* S=3 vectors *************
00904   static G4bool iNin3[nEl]={false,false,false,false,false,false,false,
00905     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00906     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00907     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00908     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00909     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00910     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00911     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00912   static G4double VZ3[nEl][iNR];
00913   // ********* S=2 vectors *************
00914   static G4bool iNin4[nEl]={false,false,false,false,false,false,false,
00915     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00916     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00917     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00918     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00919     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00920     false,false,false,false,false,false,false,false,false,false,false,false,false,false,
00921     false,false,false,false,false,false,false,false,false,false,false,false,false,false};
00922   static G4double VZ4[nEl][iNR];
00923   //
00924 #ifdef qdebug
00925   static G4int Smin=-1; // Dynamic Associated memory is implemented for S>-2 nuclei
00926   static G4int Smax= 2; // Dynamic Associated memory is implemented for S< 3 nuclei
00927   static G4int NZmin= 0; // Dynamic Associated memory is implemented for Z>-1 nuclei
00928   static G4int NNmin= 0; // Dynamic Associated memory is implemented for N>-1 nuclei
00929   static G4int NZS1max= 0; // Dynamic Associated memory is implemented for S<3, Z=-1, N<2
00930   static G4int NNS1max= 0; // Dynamic Associated memory is implemented for S<3, Z=-1, N<2 
00931 #endif
00932   // -------------------------------------------------------------------------------------
00933   G4double rm=0.;
00934   G4int nz=n+z;
00935   G4int zns=nz+s_value;
00936   if(nz+s_value<0)
00937   {
00938     z=-z;
00939     n=-n;
00940     s_value=-s_value;
00941     nz=-nz;
00942     zns=-zns;
00943   }
00944   if(z<0)
00945   {
00946     if(z==-1)
00947     {
00948       if(!s_value)
00949       {
00950         if(n==1)      return mPiC;              // pi-
00951         else          return mPiC+(n-1)*mNeut;  // Delta- + (N-1)*n
00952       }
00953       else if(s_value==1)                       // Strange negative hadron
00954       {
00955         if(!n)        return mKM;               // K-
00956         else if(n==1) return mSiM;              // Sigma-
00957         else if(n==2) return mSmN ;             // Sigma- + n DiBaryon
00958         else if(n==3) return mSmNN;             // Sigma- +2n TriBaryon
00959         else          return mSigM+mNeut*(n-1); // Sigma- + (N-1)*n
00960       }
00961       else if(s_value==2)                       // --> Double-strange negative hadrons
00962       {
00963         if(!n)        return mKsM;              // Ksi-
00964         else if(n==1) return mKsiM+mNeut;       // Ksi- + n
00965         else if(n==2) return mKsiM+mNeut+mNeut; // Ksi- + 2n
00966         else          return mKsiM+mNeut*n;     // Ksi- + Z*n
00967       }
00968       else if(s_value==-2)
00969       {
00970         if     (nz==2)         return mDiKZ+mPiC; // 2K0 + Pi-
00971         else                   return mDiKZ+mPiC+(nz-2)*mProt;
00972       }
00973       else if(s_value==3)                       // --> Triple-strange negative hadrons
00974       {
00975         if     (n==-1) return mOmM;       // Triple-strange Omega minus
00976         else if(!n   ) return mOmN;       // Triple-strange (Omega-) + n DiBaryon
00977         else if(n==-2) return mDiKZ+mKM;  // Triple-strange K- + 2*K0
00978         else           return mOmeg+mNeut*(n+2);
00979       }
00980       else if(s_value==4)
00981       {
00982         if(n==-2)      return mOmeg+mKM;  // Omega- + K-
00983         else if(n==-1) return mOmeg+mLamb;// Omega- + Lambda
00984         else           return mOmeg+mLamb+(n+1)*mNeut; // Omega- + Lambda + (n+1)*Neutrons
00985       }
00986       else if(!n)            return mOmeg+(s_value-2)*mLamb; // Multy-Lambda + Omega minus
00987       else
00988       {
00989 #ifdef qdebug
00990         if(s_value>NZS1max)
00991         {
00992           NZS1max=s_value;
00993           G4cout<<"-------->>G4QPDGCode::GetNucMass: Z=-1, S="<<s_value<<">2 with N="<<n<<G4endl;
00994         }
00995 #endif
00996           return CalculateNuclMass(z,n,s_value);
00997       }
00998     }
00999     else if(!s_value)
01000     {
01001       if(!nz)
01002       {
01003         if(n==2)             return mDiPi;
01004         else                 return mDiPi+(n-2)*mPiC;
01005       }
01006       else                   return mNeut*nz-z*mPiC+anb;
01007     }
01008     else if(!zns)           // !!! s=0 is treated above !!
01009     {
01010       if(s_value>0)          return anb+s_value*mKM+n*mPiC; // s*K- + n*Pi-
01011       else                   return anb-s_value*mKZ-z*mPiC; // (-s)*aK0 + (-z)*Pi-
01012     }
01013     else if(s_value==1)
01014     {
01015       if(!nz)
01016       {
01017         if(n==2)             return mSmPi;
01018         else                 return mSmPi+(n-2)*mPiC;
01019       }
01020       else                   return mSigM+nz*mNeut-(z+1)*mPiC;
01021     }
01022     else if(s_value==-1)     return mKZa-z*mPiC+(nz-1)*mNeut; // aK0+(nz-1)n+(-z)*Pi-
01023     else if(s_value==2)
01024     {
01025       if     (nz==-1)        return mKsiM+n*mPiC;
01026       else if(!nz)           return mKsiM+mNeut-(z+1)*mPiC;
01027       else                   return mKsiM+(nz+1)*mNeut-(z+1)*mPiC;
01028     }
01029     else if(s_value==-2)     return mDiKZ-z*mPiC+(nz-2)*mNeut;
01030     else if(s_value==3)
01031     {
01032       if     (nz==-2)        return mOmeg+(n+1)*mPiC; // Omega- + (n+1)* Pi-
01033       else if(nz==-1)        return mOmeg+mNeut+n*mPiC; // Omega- + n * Pi-
01034       else if(!nz)           return mOmeg+mDiNt+(n-1)*mPiC; // Omega- + 2N + (n-1)*Pi-
01035       else                   return mOmeg+(nz+2)*mProt-(z+1)*mPiC;
01036     }
01037     else if(s_value<-2)      return anb-s_value*mKZ-z*mPiC+(nz+s_value)*mNeut;
01038     else if(s_value==4)
01039     {
01040       if     (nz==-3)        return mOmeg+mKM+(n+1)*mPiC;        // Om- + K- + (n+1)*Pi-
01041       else if(nz==-2)        return mOmeg+mSigM+n*mPiC;          // Om- + Sig- + n*Pi-
01042       else if(nz==-1)        return mOmeg+mSigM+mNeut+(n-1)*mPiC;//Om- + Sig- +N +(n-1)*Pi-
01043       else if(!nz)           return mOmeg+mSigM+mDiNt+(n-2)*mPiC;//Om- +Sig- +2N +(n-2)*Pi-
01044       else                   return mOmeg+mSigM+(nz+2)*mDiNt-(z+2)*mPiC;//+(nz+2)N-(z+2)Pi-
01045     }
01046     // s=5: 2*K-, Ksi-; s=6: 3*K-, Omega-; s>6 adds K- and Sigma- instead of Protons
01047     else                     // !!All s<0 are done and s>4 can be easy extended if appear!!
01048     {
01049 #ifdef qdebug
01050       if(z<NZmin)
01051       {
01052         NZmin=z;
01053         G4cout<<"---->>G4QPDGCode::GetNucMass: Z="<<z<<"<-1 with N="<<n<<", S="<<s_value<<G4endl;
01054       }
01055 #endif
01056       return CalculateNuclMass(z,n,s_value);
01057     }
01058   }
01059   else if(n<0)
01060   {
01061     if(n==-1)
01062     {
01063       if(!s_value)
01064       {
01065         if(z==1)      return mPiC;              // pi+
01066         else          return mPiC+(z-1)*mProt;  // Delta++ + (Z-1)*p
01067       }
01068       else if(s_value==1)                       // --> Strange neutral hadrons
01069       {
01070         if(!z)        return mKZ;               // K0
01071         else if(z==1) return mSiP;              // Sigma+
01072         else if(z==2) return mSpP ;             // Sigma+ + p DiBaryon
01073         else if(z==3) return mSpPP;             // Sigma+ +2p TriBaryon
01074         else          return mSigP+mProt*(z-1); // Sigma+ + (Z-1)*p
01075       }
01076       else if(s_value==2)                       // --> Double-strange negative hadrons
01077       {
01078         if(!z)        return mKsZ;              // Ksi0
01079         else if(z==1) return mKsiZ+mProt;       // Ksi- + p
01080         else if(z==2) return mKsiZ+mProt+mProt; // Ksi- + 2p
01081         else          return mKsiZ+mProt*z;     // Ksi- + Z*p
01082       }
01083       else if(s_value==-2)
01084       {
01085         if     (nz==2)         return mDiKM+mPiC; // 2K+ + Pi+
01086         else                   return mDiKM+mPiC+(nz-2)*mProt;
01087       }
01088       else if(s_value==3)
01089       {
01090         if(z==1) return mOmeg+mDiPr;
01091         else     return mOmeg+(z+1)*mProt;
01092       }
01093       else if(s_value==4) return mOmeg+mLamb+(z+1)*mProt;
01094       else if(!z) return mKZa+(s_value-1)*mLamb;      // Multy-Lambda + K0
01095       else
01096       {
01097 #ifdef qdebug
01098         if(s_value>NNS1max)
01099         {
01100           NNS1max=s_value;
01101           G4cout<<"-------->>G4QPDGCode::GetNucMass: N=-1, S="<<s_value<<">2 with Z="<<z<<G4endl;
01102         }
01103 #endif
01104           return CalculateNuclMass(z,n,s_value);
01105       }
01106     }
01107     else if(!s_value)
01108     {
01109       if(!nz)
01110       {
01111         if(z==2)             return mDiPi;
01112         else                 return mDiPi+(z-2)*mPiC;
01113       }
01114       else                   return mProt*nz-n*mPiC+anb;
01115     }
01116     else if(!zns)           // !!! s=0 is treated above !!
01117     {
01118       if(s_value>0)          return anb+s_value*mKZ+z*mPiC; // s*K0 + n*Pi+
01119       else                   return anb-s_value*mKM-n*mPiC; // (-s)*aK+ + (-n)*Pi+
01120     }
01121     else if(s_value==1)
01122     {
01123       if(!nz)
01124       {
01125         if(z==2)             return mSpPi;
01126         else                 return mSpPi+(z-2)*mPiC;
01127       }
01128       else                   return mSigP+nz*mProt-(n+1)*mPiC;
01129     }
01130     else if(s_value==-1)     return mKMa-n*mPiC+(nz-1)*mProt; // K+ + (nz-1)*P + (-n)*Pi+
01131     else if(s_value==2)
01132     {
01133       if     (nz==-1)        return mKsiZ+z*mPiC;
01134       else if(!nz)           return mKsiZ+mProt-(n+1)*mPiC;
01135       else                   return mKsiZ+(nz+1)*mProt-(n+1)*mPiC;
01136     }
01137     else if(s_value==-2)     return mDiKM-n*mPiC+(nz-2)*mProt;
01138     else if(s_value==3)
01139     {
01140       if     (nz==-2)        return mOmeg+(z+1)*mPiC;       // Omega + (z+1)*Pi+
01141       else if(nz==-1)        return mOmeg+mProt+z*mPiC;     // Omega- + P +z*Pi+
01142       else if(!nz)           return mOmeg+mDiPr+(z-1)*mPiC; // Omega- + 2P + (z-1)* Pi+
01143       else                   return mOmeg+(nz+2)*mProt-(n+1)*mPiC;
01144     }
01145     else if(s_value<-2)      return anb-s_value*mKM-n*mPiC+(nz+s_value)*mProt;
01146     else if(s_value==4)
01147     {
01148       if     (nz==-3)        return mOmeg+mKZ+(z+1)*mPiC;        // Om- + K0 + (z+1)*Pi+
01149       else if(nz==-2)        return mOmeg+mSigP+z*mPiC;          // Om- + Sig+ + z*Pi+
01150       else if(nz==-1)        return mOmeg+mSigP+mProt+(z-1)*mPiC;// Om- +Sig+ +P +(z-1)*Pi+
01151       else if(!nz)           return mOmeg+mSigP+mDiPr+(z-2)*mPiC;//Om- +Sig+ +2P +(z-2)*Pi+
01152       else                   return mOmeg+mSigP+(nz+2)*mProt-(n+2)*mPiC;//+(nz+2)P-(n+2)Pi+
01153     }
01154     // s=5: 2*KZ, Ksi0; s=6: 3*KZ, Omega-; s>6 adds K0 and Sigma+ instead of Protons
01155     else                     // !!All s<0 are done and s>4 can be easy extended if appear!!
01156     {
01157 #ifdef qdebug
01158       if(n<NNmin)
01159       {
01160         NNmin=n;
01161         G4cout<<"---->>G4QPDGCode::GetNucMass: N="<<n<<"<-1 with Z="<<z<<", S="<<s_value<<G4endl;
01162       }
01163 #endif
01164       return CalculateNuclMass(z,n,s_value);
01165     }
01166   }
01167   //return CalculateNuclMass(z,n,s_value); // @@ This is just to compare the calculation time @@
01168   if(nz >= 256 || z >= nEl) return 256000.;
01169   if(!s_value)             // **************> S=0 nucleus
01170   {
01171     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01172     if(!iNin0[z])          // =--=> This element is already initialized
01173     {
01174 #ifdef idebug
01175       G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=0 is initialized. F="<<iNin0[z]<<G4endl;
01176 #endif
01177       G4int iNfin=iNL[z];
01178       if(iNfin>iNR) iNfin=iNR;
01179       for (G4int in=0; in<iNfin; in++) VZ0[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01180       iNin0[z]=true;
01181     }
01182     G4int dNn=n-Nmin;
01183     if(dNn<0)              // --> The minimum N must be reduced 
01184     {
01185 #ifdef qdebug
01186       if(n<iNmin[z])
01187       {
01188         G4cout<<"-->>G4QPDGCode::GetNucM:Z="<<z<<", S=0 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01189         iNmin[z]=n;
01190       }
01191 #endif
01192       return CalculateNuclMass(z,n,s_value);
01193     }
01194     else if(dNn<iNL[z]) return VZ0[z][dNn]; // Found in DAM
01195     else                   // --> The maximum N must be increased
01196     {
01197 #ifdef qdebug
01198       if(dNn>iNmax)
01199       {
01200         G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNmax<<G4endl;
01201         iNmax=dNn;
01202       }
01203       if(dNn>iNran[z])
01204       {
01205         G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01206         iNran[z]=dNn;
01207       }
01208 #endif
01209       return CalculateNuclMass(z,n,s_value);
01210     }
01211   }
01212   else if(s_value==1)      // ******************> S=1 nucleus
01213   {
01214     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01215 #ifdef pdebug
01216     G4bool pPrint = !z && !n;
01217     if(pPrint)
01218       G4cout<<"G4QPDGC::GetNucM:Nmin="<<Nmin<<",iNin1="<<iNin1[0]<<",iNL="<<iNL[0]<<G4endl;
01219 #endif
01220     if(!iNin1[z])          // =--=> This element is already initialized
01221     {
01222 #ifdef pdebug
01223       if(pPrint)
01224       G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=1 is initialized. F="<<iNin1[z]<<G4endl;
01225 #endif
01226       G4int iNfin=iNL[z];
01227       if(iNfin>iNR) iNfin=iNR;
01228       for (G4int in=0; in<iNfin; in++) VZ1[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01229       iNin1[z]=true;
01230     }
01231     G4int dNn=n-Nmin;
01232     if(dNn<0)              // --> The minimum N must be reduced 
01233     {
01234 #ifdef qdebug
01235       if(n<iNmin[z])
01236       {
01237         G4cout<<"-->>G4QPDGCode::GetNucM:Z="<<z<<", S=1 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01238         iNmin[z]=n;
01239       }
01240 #endif
01241       return CalculateNuclMass(z,n,s_value);
01242     }
01243     else if(dNn<iNL[z]) return VZ1[z][dNn]; // Found in DAM
01244     else                   // --> The maximum N must be increased
01245     {
01246 #ifdef qdebug
01247       if(dNn>iNmax)
01248       {
01249         G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNmax<<G4endl;
01250         iNmax=dNn;
01251       }
01252       if(dNn>iNran[z])
01253       {
01254         G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01255         iNran[z]=dNn;
01256       }
01257 #endif
01258       return CalculateNuclMass(z,n,s_value);
01259     }
01260   }
01261   else if(s_value==-1)     // ******************> S=-1 nucleus
01262   {
01263     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01264     if(!iNin9[z])          // =--=> This element is already initialized
01265     {
01266 #ifdef idebug
01267       G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=-1 is initialized. F="<<iNin9[z]<<G4endl;
01268 #endif
01269       G4int iNfin=iNL[z];
01270       if(iNfin>iNR) iNfin=iNR;
01271       for (G4int in=0; in<iNfin; in++) VZ9[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01272       iNin9[z]=true;
01273     }
01274     G4int dNn=n-Nmin;
01275     if(dNn<0)              // --> The minimum N must be reduced 
01276     {
01277 #ifdef qdebug
01278       if(n<iNmin[z])
01279       {
01280         G4cout<<"->>G4QPDGCode::GetNucM:Z="<<z<<" ,S=-1 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01281         iNmin[z]=n;
01282       }
01283 #endif
01284       return CalculateNuclMass(z,n,s_value);
01285     }
01286     else if(dNn<iNL[z]) return VZ9[z][dNn]; // Found in DAM
01287     else                   // --> The maximum N must be increased
01288     {
01289 #ifdef qdebug
01290       if(dNn>iNmax)
01291       {
01292         G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNmax<<G4endl;
01293         iNmax=dNn;
01294       }
01295       if(dNn>iNran[z])
01296       {
01297         G4cout<<"G4QPDGCode::GetNucM:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01298         iNran[z]=dNn;
01299       }
01300 #endif
01301       return CalculateNuclMass(z,n,s_value);
01302     }
01303   }
01304   else if(s_value==2)      // ******************> S=2 nucleus
01305   {
01306     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01307     if(!iNin2[z])          // =--=> This element is already initialized
01308     {
01309 #ifdef idebug
01310       G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=2 is initialized. F="<<iNin2[z]<<G4endl;
01311 #endif
01312       G4int iNfin=iNL[z];
01313       if(iNfin>iNR) iNfin=iNR;
01314       for (G4int in=0; in<iNfin; in++) VZ2[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01315       iNin2[z]=true;
01316     }
01317     G4int dNn=n-Nmin;
01318     if(dNn<0)              // --> The minimum N must be reduced 
01319     {
01320 #ifdef qdebug
01321       if(n<iNmin[z])
01322       {
01323         G4cout<<"-->>G4QPDGCode::GetNucM:Z="<<z<<", S=2 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01324         iNmin[z]=n;
01325       }
01326 #endif
01327       return CalculateNuclMass(z,n,s_value);
01328     }
01329     else if(dNn<iNL[z]) return VZ2[z][dNn]; // Found in DAM
01330     else                   // --> The maximum N must be increased
01331     {
01332 #ifdef qdebug
01333       if(dNn>iNmax)
01334       {
01335         G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNmax<<G4endl;
01336         iNmax=dNn;
01337       }
01338       if(dNn>iNran[z])
01339       {
01340         G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01341         iNran[z]=dNn;
01342       }
01343 #endif
01344       return CalculateNuclMass(z,n,s_value);
01345     }
01346   }
01347   else if(s_value==-2)     // ******************> S=-2 nucleus
01348   {
01349     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01350     if(!iNin8[z])          // =--=> This element is already initialized
01351     {
01352 #ifdef idebug
01353       G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=-2 is initialized. F="<<iNin8[z]<<G4endl;
01354 #endif
01355       G4int iNfin=iNL[z];
01356       if(iNfin>iNR) iNfin=iNR;
01357       for (G4int in=0; in<iNfin; in++) VZ8[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01358       iNin8[z]=true;
01359     }
01360     G4int dNn=n-Nmin;
01361     if(dNn<0)              // --> The minimum N must be reduced 
01362     {
01363 #ifdef qdebug
01364       if(n<iNmin[z])
01365       {
01366         G4cout<<"->>G4QPDGCode::GetNucM:Z="<<z<<", S=-2 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01367         iNmin[z]=n;
01368       }
01369 #endif
01370       return CalculateNuclMass(z,n,s_value);
01371     }
01372     else if(dNn<iNL[z]) return VZ8[z][dNn]; // Found in DAM
01373     else                   // --> The maximum N must be increased
01374     {
01375 #ifdef qdebug
01376       if(dNn>iNmax)
01377       {
01378         G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNmax<<G4endl;
01379         iNmax=dNn;
01380       }
01381       if(dNn>iNran[z])
01382       {
01383         G4cout<<"G4QPDGCode::GetNucM:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01384         iNran[z]=dNn;
01385       }
01386 #endif
01387       return CalculateNuclMass(z,n,s_value);
01388     }
01389   }
01390   else if(s_value==-3)     // ******************> S=-3 nucleus
01391   {
01392     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01393     if(!iNin7[z])          // =--=> This element is already initialized
01394     {
01395 #ifdef idebug
01396       G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=-3 is initialized. F="<<iNin7[z]<<G4endl;
01397 #endif
01398       G4int iNfin=iNL[z];
01399       if(iNfin>iNR) iNfin=iNR;
01400       for (G4int in=0; in<iNfin; in++) VZ7[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01401       iNin7[z]=true;
01402     }
01403     G4int dNn=n-Nmin;
01404     if(dNn<0)              // --> The minimum N must be reduced 
01405     {
01406 #ifdef qdebug
01407       if(n<iNmin[z])
01408       {
01409         G4cout<<"->>G4QPDGCode::GetNucM:Z="<<z<<", S=-3 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01410         iNmin[z]=n;
01411       }
01412 #endif
01413       return CalculateNuclMass(z,n,s_value);
01414     }
01415     else if(dNn<iNL[z]) return VZ7[z][dNn]; // Found in DAM
01416     else                   // --> The maximum N must be increased
01417     {
01418 #ifdef qdebug
01419       if(dNn>iNmax)
01420       {
01421         G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNmax<<G4endl;
01422         iNmax=dNn;
01423       }
01424       if(dNn>iNran[z])
01425       {
01426         G4cout<<"G4QPDGCode::GetNucM:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01427         iNran[z]=dNn;
01428       }
01429 #endif
01430       return CalculateNuclMass(z,n,s_value);
01431     }
01432   }
01433   else if(s_value==3)      // ******************> S=3 nucleus
01434   {
01435     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01436     if(!iNin3[z])          // =--=> This element is already initialized
01437     {
01438 #ifdef idebug
01439       G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=3 is initialized. F="<<iNin3[z]<<G4endl;
01440 #endif
01441       G4int iNfin=iNL[z];
01442       if(iNfin>iNR) iNfin=iNR;
01443       for (G4int in=0; in<iNfin; in++) VZ3[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01444       iNin3[z]=true;
01445     }
01446     G4int dNn=n-Nmin;
01447     if(dNn<0)              // --> The minimum N must be reduced 
01448     {
01449 #ifdef qdebug
01450       if(n<iNmin[z])
01451       {
01452         G4cout<<"-->>G4QPDGCode::GetNucM:Z="<<z<<", S=3 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01453         iNmin[z]=n;
01454       }
01455 #endif
01456       return CalculateNuclMass(z,n,s_value);
01457     }
01458     else if(dNn<iNL[z]) return VZ3[z][dNn]; // Found in DAM
01459     else                   // --> The maximum N must be increased
01460     {
01461 #ifdef qdebug
01462       if(dNn>iNmax)
01463       {
01464         G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNmax<<G4endl;
01465         iNmax=dNn;
01466       }
01467       if(dNn>iNran[z])
01468       {
01469         G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01470         iNran[z]=dNn;
01471       }
01472 #endif
01473       return CalculateNuclMass(z,n,s_value);
01474     }
01475   }
01476   else if(s_value==-4)     // ******************> S=-4 nucleus
01477   {
01478     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01479     if(!iNin6[z])          // =--=> This element is already initialized
01480     {
01481 #ifdef idebug
01482       G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=-4 is initialized. F="<<iNin6[z]<<G4endl;
01483 #endif
01484       G4int iNfin=iNL[z];
01485       if(iNfin>iNR) iNfin=iNR;
01486       for (G4int in=0; in<iNfin; in++) VZ6[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01487       iNin6[z]=true;
01488     }
01489     G4int dNn=n-Nmin;
01490     if(dNn<0)              // --> The minimum N must be reduced 
01491     {
01492 #ifdef qdebug
01493       if(n<iNmin[z])
01494       {
01495         G4cout<<"->>G4QPDGCode::GetNucM:Z="<<z<<", S=-4 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01496         iNmin[z]=n;
01497       }
01498 #endif
01499       return CalculateNuclMass(z,n,s_value);
01500     }
01501     else if(dNn<iNL[z]) return VZ6[z][dNn]; // Found in DAM
01502     else                   // --> The maximum N must be increased
01503     {
01504 #ifdef qdebug
01505       if(dNn>iNmax)
01506       {
01507         G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNmax<<G4endl;
01508         iNmax=dNn;
01509       }
01510       if(dNn>iNran[z])
01511       {
01512         G4cout<<"G4QPDGCode::GetNucM:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01513         iNran[z]=dNn;
01514       }
01515 #endif
01516       return CalculateNuclMass(z,n,s_value);
01517     }
01518   }
01519   else if(s_value==4)      // ******************> S=4 nucleus
01520   {
01521     G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
01522     if(!iNin4[z])          // =--=> This element is already initialized
01523     {
01524 #ifdef idebug
01525       G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=4 is initialized. F="<<iNin4[z]<<G4endl;
01526 #endif
01527       G4int iNfin=iNL[z];
01528       if(iNfin>iNR) iNfin=iNR;
01529       for (G4int in=0; in<iNfin; in++) VZ4[z][in] = CalculateNuclMass(z,in+Nmin,s_value);
01530       iNin4[z]=true;
01531     }
01532     G4int dNn=n-Nmin;
01533     if(dNn<0)              // --> The minimum N must be reduced 
01534     {
01535 #ifdef qdebug
01536       if(n<iNmin[z])
01537       {
01538         G4cout<<"-->>G4QPDGCode::GetNucM:Z="<<z<<", S=4 with N="<<n<<"<"<<iNmin[z]<<G4endl;
01539         iNmin[z]=n;
01540       }
01541 #endif
01542       return CalculateNuclMass(z,n,s_value);
01543     }
01544     else if(dNn<iNL[z]) return VZ4[z][dNn]; // Found in DAM
01545     else                   // --> The maximum N must be increased
01546     {
01547 #ifdef qdebug
01548       if(dNn>iNmax)
01549       {
01550         G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNmax<<G4endl;
01551         iNmax=dNn;
01552       }
01553       if(dNn>iNran[z])
01554       {
01555         G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
01556         iNran[z]=dNn;
01557       }
01558 #endif
01559       return CalculateNuclMass(z,n,s_value);
01560     }
01561   }
01562   else
01563   {
01564 #ifdef qdebug
01565     if(s_value<Smin || s_value>Smax)
01566     {
01567       if(s_value<Smin) Smin=s_value;
01568       if(s_value>Smax) Smax=s_value;
01569       G4cout<<">>G4QPDGCode::GetNucM:Z="<<z<<" with S="<<s_value<<",N="<<n<<" (Improve)"<<G4endl;
01570     }
01571 #endif
01572     rm=CalculateNuclMass(z,n,s_value);
01573   }
01574 #ifdef debug
01575   G4cout<<"G4QPDGCode::GetMass:GetNuclMass="<<rm<<",Z="<<z<<",N="<<n<<",S="<<s_value<<G4endl;
01576 #endif
01577   return rm;
01578 } // End of GetNuclearMass
01579 
01580 // Calculate a nuclear mass for Z (a#of protons), N (a#of neutrons), & S (a#of lambdas) 
01581 G4double G4QPDGCode::CalculateNuclMass(G4int z, G4int n, G4int s_value)
01582 {
01583   static const G4double mP  = QHaM(21); // Proton
01584   static const G4double mN  = QHaM(20); // Neutron
01585   static const G4double mL  = QHaM(22); // Lambda
01586   static const G4double mD  = G4Deuteron::Deuteron()->GetPDGMass(); // Deuteron (H-2)
01587   static const G4double mT  = G4Triton::Triton()->GetPDGMass();     // Triton (H-3)
01588   static const G4double mHe3= G4He3::He3()->GetPDGMass();           // Hetrium (He-3)
01589   static const G4double mAl = G4Alpha::Alpha()->GetPDGMass();       // Alpha (He-4)
01590   static const G4double dmP = mP+mP;    // DiProton
01591   static const G4double dmN = mN+mN;    // DiNeutron
01592   static const G4double dmL = mL+mL;    // DiLambda
01593   static const G4double dLN = mL+mN;    // LambdaNeutron
01594   static const G4double dLP = mL+mP;    // LambdaProton
01595   static const G4double mSm = QHaM(23); // Sigma-
01596   static const G4double mSp = QHaM(25); // Sigma+
01597   static const G4double dSP = mSp+mP;   // ProtonSigma+
01598   static const G4double dSN = mSm+mN;   // NeutronSigma-
01599   static const G4double dnS = dSN+mN;   // 2NeutronsSigma-
01600   static const G4double dpS = dSP+mP;   // 2ProtonsSigma+
01601   static const G4double mXm = QHaM(26); // Ksi-
01602   static const G4double mXz = QHaM(27); // Ksi0
01603   static const G4double mOm = QHaM(44); // Omega-
01604   static const G4double dXN = mXm+mN;   // NeutronKsi-
01605   static const G4double dXP = mXz+mP;   // ProtonKsi0
01606   static const G4double dOP = mOm+mP;   // ProtonOmega-
01607   static const G4double dON = mOm+mN;   // NeutronOmega-
01608   static const G4double mK  = QHaM(18); // Charged Kaon
01609   static const G4double mK0 = QHaM(17); // Neutral Kaon
01610   static const G4double mPi = QHaM(15); // Charged Pion
01612   //static const G4int    nSh = 164;
01613   //static G4double sh[nSh] = {0.,                        // Fake element for C++ N=Z=0
01614   //                             -4.315548,   2.435504,  -1.170501,   3.950887,   5.425238,
01615   //                             13.342524,  15.547586,  22.583129,  23.983480,  30.561036,
01616   //                             33.761971,  41.471027,  45.532156,  53.835880,  58.495514,
01617   //                             65.693445,  69.903344,  76.899581,  81.329361,  88.979438,
01618   //                             92.908703, 100.316636, 105.013393, 113.081686, 118.622601,
01619   //                            126.979113, 132.714435, 141.413182, 146.433488, 153.746754,
01620   //                            158.665225, 165.988967, 170.952395, 178.473011, 183.471531,
01621   //                            191.231310, 196.504414, 204.617158, 210.251108, 218.373984,
01622   //                            223.969281, 232.168660, 237.925619, 246.400505, 252.392471,
01623   //                            260.938644, 267.191321, 276.107788, 282.722682, 291.881502,
01624   //                            296.998590, 304.236025, 309.562296, 316.928655, 322.240263,
01625   //                            329.927236, 335.480630, 343.233705, 348.923475, 356.911659,
01626   //                            362.785757, 370.920926, 376.929998, 385.130316, 391.197741,
01627   //                            399.451554, 405.679971, 414.101869, 420.346260, 428.832412,
01628   //                            435.067445, 443.526983, 449.880034, 458.348602, 464.822352,
01629   //                            473.313779, 479.744524, 488.320887, 495.025069, 503.841579,
01630   //                            510.716379, 519.451976, 525.036156, 532.388151, 537.899017,
01631   //                            545.252264, 550.802469, 558.402181, 564.101100, 571.963429,
01632   //                            577.980340, 586.063802, 592.451334, 600.518525, 606.832027,
01633   //                            614.831626, 621.205330, 629.237413, 635.489106, 643.434167,
01634   //                            649.691284, 657.516479, 663.812101, 671.715021, 678.061128,
01635   //                            686.002970, 692.343712, 700.360477, 706.624091, 714.617848,
01636   //                            721.100390, 729.294717, 735.887170, 744.216084, 751.017094,
01637   //                            759.551764, 766.377807, 775.080204, 781.965673, 790.552795,
01638   //                            797.572494, 806.088030, 813.158751, 821.655631, 828.867137,
01639   //                            836.860955, 842.183292, 849.195302, 854.731798, 861.898839,
01640   //                            867.783606, 875.313342, 881.443441, 889.189065, 895.680189,
01641   //                            903.679729, 910.368085, 918.579876, 925.543547, 933.790028,
01642   //                            940.811396, 949.122548, 956.170201, 964.466810, 971.516490,
01643   //                            979.766905, 986.844659, 995.113552,1002.212760,1010.418770,
01644   //                           1018.302560,1025.781870,1033.263560,1040.747880,1048.234460,
01645   //                           1055.723430,1063.214780,1070.708750,1078.204870,1085.703370,
01646   //                           1093.204260,1100.707530,1108.213070};
01647   //static const G4double b1=8.09748564; // MeV
01648   //static const G4double b2=-0.76277387;
01649   //static const G4double b3=83.487332;  // MeV
01650   //static const G4double b4=0.090578206;// 2*b4
01651   //static const G4double b5=0.676377211;// MeV
01652   //static const G4double b6=5.55231981; // MeV
01653   static const G4double b7=25.;        // MeV (Lambda binding energy predexponent)
01654   // --- even-odd difference is 3.7(MeV)/X
01655   // --- S(X>151)=-57.56-5.542*X**1.05
01656   static const G4double b8=10.5;       // (Lambda binding energy exponent)
01657   //static const G4double b9=-1./3.;
01658   static const G4double a2=0.13;       // LambdaBindingEnergy for deutron+LambdaSystem(MeV)
01659   static const G4double a3=2.2;        // LambdaBindingEnergy for (t/He3)+LambdaSystem(MeV)
01660   static const G4double um_value=931.49432;  // Umified atomic mass unit (MeV)
01661   //static const G4double me =0.511;     // electron mass (MeV) :: n:8.071, p:7.289
01662   static const G4double eps =0.0001;   // security addition for multybaryons
01663   //static G4double c[9][9]={// z=1     =2     =3     =4     =5     =6     =7     =8     =9
01664   //                 {13.136,14.931,25.320,38.000,45.000,55.000,65.000,75.000,85.000},//n=1
01665   //     {14.950, 2.425,11.680,18.374,27.870,35.094,48.000,60.000,72.000},  //n=2
01666   //     {25.930,11.390,14.086,15.770,22.921,28.914,39.700,49.000,60.000},  //n=3
01667   //     {36.830,17.594,14.908, 4.942,12.416,15.699,24.960,32.048,45.000},  //n=4
01668   //     {41.860,26.110,20.946,11.348,12.051,10.650,17.338,23.111,33.610},  //n=5
01669   //     {45.000,31.598,24.954,12.607, 8.668, 0.000, 5.345, 8.006,16.780},  //n=6
01670   //     {50.000,40.820,33.050,20.174,13.369, 3.125, 2.863, 2.855,10.680},  //n=7
01671   //     {55.000,48.810,40.796,25.076,16.562, 3.020, 0.101,-4.737,1.9520},  //n=8
01672   //     {60.000,55.000,50.100,33.660,23.664, 9.873, 5.683,-0.809,0.8730}}; //n=9
01673   if(z>107)
01674   {
01675 #ifdef debug
01676     G4cout<<"***G4QPDGCode::CalcNuclMass: Z="<<z<<">107, N="<<n<<", S="<<s_value<<G4endl;
01677 #endif
01678     return 256000.;
01679   }
01680   G4int Z=z;
01681   G4int N=n;
01682   G4int S=s_value;
01683 #ifdef debug
01684   G4cout<<"G4QPDGCode::CalcNuclMass called with Z="<<Z<<",N="<<N<<", S="<<S<<G4endl;
01685 #endif
01686   if(!N&&!Z&&!S)
01687   {
01688 #ifdef debug
01689     //G4cout<<"G4QPDGCode::GetNuclMass(0,0,0)="<<mPi0<<"#0"<<G4endl;
01690 #endif
01691     //return mPi0; // Pi0 mass (dangerose)
01692     return 0.;     // Photon mass
01693   }
01694   else if(!N&&!Z&&S==1) return mL;
01695   else if(!N&&Z==1&&!S) return mP;
01696   else if(N==1&&!Z&&!S) return mN;
01697   else if(!N&&!Z&&S>1) return mL*S+eps;
01698   else if(!N&&Z>1&&!S) return mP*Z+eps;
01699   else if(N>1&&!Z&&!S) return mN*N+eps;
01700   G4int A=Z+N;
01701   G4int Bn=A+S;
01702   //if((Z<0||N<0)&&!Bn)
01703   //{
01704   //  if     (N<0) return Bn*mL-Z*mK - N*mK0+eps*   S ;
01705   //  else              return Bn*mL+N*mPi-A*mK +eps*(N+S);
01706   //}
01707   //if(A<0&&Bn>=0)                      // Bn*LAMBDA's + (-(N+Z))*antiKaons
01708   //{
01709   //  if     (N<0&&Z<0) return Bn*mL-Z*mK -N*mK0+eps*   S ;
01710   //  else if(N<0)      return Bn*mL+Z*mPi-A*mK0+eps*(Z+S);
01711   //  else              return Bn*mL+N*mPi-A*mK +eps*(N+S);
01712   //}
01713   if(!Bn)                     // => "GS Mesons - octet" case (without eta&pi0)
01714   {
01715     if     (!S&&Z<0) return mPi*N;
01716     else if(!S&&N<0) return mPi*Z;
01717     else if ( (N == 1 && S == -1) || (N == -1 && S == 1) ) 
01718       return mK0; // Simple decision
01719     else if ( (S == 1 && Z == -1) || (S == -1 && Z == 1) ) 
01720       return mK;  // Simple decision
01721     else if(S>0)                                  // General decision
01722     {
01723       if     (-Z>S) return S*mK-(S+Z)*mPi+eps;
01724       else if(Z>=0) return S*mK0+Z*mPi+eps;
01725       else          return (S+Z)*mK0-Z*mK+eps;
01726     }
01727     else if(S<0)                                  // General decision
01728     {
01729       if     (Z>-S) return -S*mK+(S+Z)*mPi+eps;
01730       else if(Z<=0) return -S*mK0-Z*mPi+eps;
01731       else          return -(S+Z)*mK0+Z*mK+eps;
01732     }
01733   }
01734   else if(Bn==1)                    // => "GS Baryons - octet" case (withoutSigma0)
01735   {
01736     if     (Z== 1 && N== 0 && S== 0) return mP;
01737     else if(Z== 0 && N== 1 && S== 0) return mN;
01738     else if(Z== 0 && N== 0 && S== 1) return mL;
01739     else if(Z== 1 && N==-1 && S== 1) return mSp;  // Lower than Sigma+ (Simp.Decis)
01740     else if(Z==-1 && N== 1 && S== 1) return mSm;  // Lower than Sigma- (Simp.Decis)
01741     else if(Z== 0 && N==-1 && S== 2) return mXz;  // Lower than Xi0    (Simp.Decis)
01742     else if(Z==-1 && N== 0 && S== 2) return mXm;  // Lower than Xi-    (Simp.Decis)
01743     else if(Z==-1 && N==-1 && S== 3) return mOm;  // Lower than Omega- (Simp.Decis)
01744     else if(!S&&Z<0) return mN-mPi*Z+eps;         // Negative Isonuclei
01745     else if(!S&&N<0) return mP-mPi*N+eps;         // Positive Isonuclei
01746     else if(S==1)                                 // --> General decision
01747     {
01748       if     (N>1)   return mSm+(N-1)*mPi+eps;    // (Sigma-)+(n*PI-)
01749       else if(Z>1)   return mSp+(Z-1)*mPi+eps;    // (Sigma+)+(n*PI+)
01750     }
01751     else if(S==2)                                 // --> General decision
01752     {
01753       if     (N>0)   return mXm+N*mPi+eps;        // (Xi-)+(n*PI-)
01754       else if(Z>0)   return mXz+Z*mPi+eps;        // (Xi0)+(n*PI+)
01755     }
01756     else if(S==3)                                 // --> General decision
01757     {
01758       if     (N>-1)  return mOm+(N+1)*mPi+eps;    // (Omega-)+(n*PI-)
01759       else if(Z>-1)  return mOm+(Z+1)*mPi+eps;    // (Omega-)+(n*PI+)
01760     }
01761     else if(S>3)                                  // --> General Omega- decision
01762     {
01763       if   (-Z>S-2)  return mOm+(S-3)*mK +(2-Z-S)*mPi+eps;
01764       else if(Z>-1)  return mOm+(S-3)*mK0+(Z+1)+mPi+eps;
01765       else           return mOm+(S+Z-2)*mK0-(Z+1)*mK+eps;
01766     }
01767   }
01768   else if(Bn==2) // => "GS Baryons - decuplet" case (NP,LP, and LN are made below)
01769   {
01770     if     (Z== 2 && N== 0 && S== 0) return dmP;
01771     //else if(Z== 1 && N== 1 && S== 0) return 1875.61282; // Exact deuteron PDG Mass
01772     else if(Z== 1 && N== 1 && S== 0) return mD;   // Exact deuteron PDG Mass
01773     else if(Z== 0 && N== 2 && S== 0) return dmN;
01774     else if(Z== 2 && N==-1 && S== 1) return dSP;
01775     else if(Z== 1 && N== 0 && S== 1) return dLP;
01776     else if(Z== 0 && N== 1 && S== 1) return dLN;
01777     else if(Z==-1 && N== 2 && S== 1) return dSN;
01778     else if(Z== 1 && N==-1 && S== 2) return dXP;
01779     else if(Z== 0 && N== 0 && S== 2) return dmL;
01780     else if(Z==-1 && N== 1 && S== 2) return dXN;
01781     else if(Z== 0 && N==-1 && S== 3) return dOP;
01782     else if(Z==-1 && N== 0 && S== 3) return dON;
01783     else if(!S&&Z<0) return dmN-mPi*Z+eps;        // Negative Isonuclei
01784     else if(!S&&N<0) return dmP-mPi*N+eps;        // Positive Isonuclei
01785     else if(S==1)                                 // --> General decision
01786     {
01787       if     (N>2)   return dSP+(N-2)*mPi+eps;    // (nSigma-)+(n*PI-)
01788       else if(Z>2)   return dSN+(Z-1)*mPi+eps;    // (pSigma+)+(n*PI+)
01789     }
01790     else if(S==2)                                 // --> General decision
01791     {
01792       if     (N>1)   return dXN+(N-1)*mPi+eps;    // (nXi-)+(n*PI-)
01793       else if(Z>1)   return dXP+(Z-1)*mPi+eps;    // (pXi0)+(n*PI+)
01794     }
01795     else if(S==3)                                 // --> General decision
01796     {
01797       if     (N>0)   return dON+N*mPi+eps;        // (nOmega-)+(n*PI-)
01798       else if(Z>0)   return dOP+Z*mPi+eps;        // (pOmega-)+(n*PI+)
01799     }
01800     else if(S>3)                                  // --> General Omega- decision
01801     {
01802       if   (-Z>S-2)  return dON+(S-3)*mK +(2-Z-S)*mPi+eps;
01803       else if(Z>0)   return dOP+(S-3)*mK0+Z+mPi+eps;
01804       else           return dOP+(S+Z-3)*mK0-Z*mK+eps;
01805     }
01806     //else if(S>0)                                // @@ Implement General Decision
01807     //{
01808     //  //#ifdef debug
01809     //  G4cout<<"***G4QPDGCode::GetNuclMass:B=2, Z="<<Z<<",N="<<N<<",S="<<S<<G4endl;
01810     //  //#endif
01811     //  return bigM;                              // Exotic dibaryons (?)
01812     //}
01813   }
01814   else if(Bn==3)
01815   {
01816     if(!S)                                      // Bn=3
01817     {
01818       if     (Z==1 && N== 2) return mT;         // tritium
01819       else if(Z==2 && N== 1) return mHe3;       // hetrium
01820     }
01821     if(S== 1 && Z==-1 && N== 3)                 // nnSig-
01822     {
01823       if     (Z==-1 && N== 3) return dnS;       // nnSig-
01824       else if(Z== 3 && N==-1) return dpS;       // ppSig+
01825     }
01826   }
01827   else if(!S)
01828   {
01829     if(Z==2 && N==2) return mAl;                // Alpha
01830     else if(Z<0) return A*mN-Z*mPi+eps;         // Multybaryon Negative Isonuclei
01831     else if(Z>A) return A*mP+(Z-A)*mPi+eps;     // Multybaryon Positive Isonuclei
01832   }
01833   // === Start mesonic extraction ===
01834   G4double km_value=0.;               // Mass Sum of K mesons (G4E::DecayAntiStrang.)
01835   G4int Zm=Z;
01836   G4int Nm=N;
01837   G4int Sm=S;
01838   if(S<0&&Bn>0)                       // NEW: the free mass minimum
01839   {
01840     if(Zm>=-S)                        // Enough charge for K+'s
01841     {
01842       km_value=-S*mK;                 // Anti-Lambdas are compensated by protons
01843       Zm+=S;
01844     }
01845     else if(Zm>0)
01846     {
01847       km_value=Zm*mK-(S+Zm)*mK0;      // Anti-Lambdas are partially compensated by neutrons
01848       Zm=0;
01849       Nm+=S+Zm;
01850     }
01851   }
01852   else Sm=0;                          // No alternative calculations
01853   // Old algorithm
01854   G4double k=0.;                      // Mass Sum of K mesons
01855   if(S<0&&Bn>0)                       // OLD @@ Can be improved by K0/K+ search of minimum
01856   {
01857     G4int sH=(-S)/2;                  // SmallHalfS || The algorithm must be the same
01858     G4int bH=-S-sH;                   // BigHalhS   || as in G4QE::DecayAntiStrange
01859     if(Z>0 && Z>N)
01860     {
01861       if(Z>=bH)                       // => "Enough protons in nucleus" case
01862       {
01863         if(N>=sH)
01864         {
01865           k=bH*mK+sH*mK0;
01866           Z-=bH;
01867           N-=sH;
01868         }
01869         else
01870         {
01871           G4int dN=Z-N;
01872           if(dN>=-S)
01873           {
01874             k=-S*mK;
01875             Z+=S;
01876           }
01877           else
01878           {
01879             G4int sS=(-S-dN)/2;
01880             G4int bS=-S-dN-sS;
01881             sS+=dN;
01882             if(Z>=sS&&N>=bS)
01883             {
01884               k=sS*mK+bS*mK0;
01885               Z-=sS;
01886               N-=bS;
01887             }
01888             else if(Z<sS)
01889             {
01890               G4int dS=-S-Z;
01891               k=Z*mK+dS*mK0;
01892               N-=dS;
01893               Z=0;
01894             }
01895             else
01896             {
01897               G4int dS=-S-N;
01898               k=dS*mK+N*mK0;
01899               Z-=dS;
01900               N=0;
01901             }
01902           }
01903         }
01904       }
01905       else // Must not be here
01906       {
01907 #ifdef debug
01908         G4cout<<"***G4QPDGC::CalcNuclMass:Antimatter? Z="<<Z<<",N="<<N<<",S="<<S<<G4endl;
01909 #endif
01910         return 0.;                          // @@ Antiparticles aren't implemented @@
01911       }
01912     }
01913     else if(N>=bH)
01914     {
01915       if(Z>=sH)
01916       {
01917         k=sH*mK+bH*mK0;
01918         Z-=sH;
01919         N-=bH;
01920       }
01921       else
01922       {
01923         G4int dN=N-Z;
01924         if(dN>=-S)
01925         {
01926           k=-S*mK0;
01927           N+=S;
01928         }
01929         else
01930         {
01931           G4int sS=(-S-dN)/2;
01932           G4int bS=-S-dN-sS;
01933           bS+=dN;
01934           if(N>=bS&&Z>=sS)
01935           {
01936             k=sS*mK+bS*mK0;
01937             Z-=sS;
01938             N-=bS;
01939           }
01940           else if(N<bS)
01941           {
01942             G4int dS=-S-N;
01943             k=dS*mK+N*mK0;
01944             Z-=dS;
01945             N=0;
01946           }
01947           else
01948           {
01949             G4int dS=-S-Z;
01950             k=Z*mK+dS*mK0;
01951             N-=dS;
01952             Z=0;
01953           }
01954         }
01955       }
01956     }
01957     else // Must not be here
01958     {
01959       return 0.;                            // @@ Antiparticles aren't implemented @@
01960 #ifdef debug
01961         G4cout<<"***G4QPDGC::CalcNuclMass:Antimatter? N="<<N<<",Z="<<Z<<",S="<<S<<G4endl;
01962 #endif
01963     }
01964     S=0;
01965   }
01966   if(N<0)
01967   {
01968     k+=-N*mPi;
01969     Z+=N;
01970     N=0;
01971   }
01972   if(Z<0)
01973   {
01974     k+=-Z*mPi;
01975     N+=Z;
01976     Z=0;
01977   }
01978   A=Z+N;
01979   if (!A) return k+S*mL+S*eps;              // @@ multy LAMBDA states are not implemented
01980   G4double m_value=k+A*um_value;            // Expected mass in atomic units
01981   //G4double D=N-Z;                         // Isotopic shift of the nucleus
01982   if ( (A+S < 1 && k==0.) || Z < 0 || N < 0 ) 
01983   {  // @@ Can be generalized to anti-nuclei
01984 #ifdef debug
01985     G4cout<<"**G4QPDGCode::CalcNuclMass:A="<<A<<"<1 || Z="<<Z<<"<0 || N="<<N<<"<0"<<G4endl;
01986     //@@throw G4QException("***G4QPDGCode::GetNuclMass: Impossible nucleus");
01987 #endif
01988     return 0.;                              // @@ Temporary
01989   }
01990   if     (!Z) return k+N*(mN+.1)+S*(mL+.1); // @@ n+LAMBDA states are not implemented
01991   else if(!N) return k+Z*(mP+1.)+S*(mL+.1); // @@ p+LAMBDA states are not implemented
01992   //else if(N<=9&&Z<=9) m_value+=1.433e-5*pow(double(Z),2.39)-Z*me+c[N-1][Z-1];// Geant4 Comp.now
01993   else 
01994   {
01995     //G4double fA=A;
01996     // IsInTable is inside G4NucleiProperties::GetNuclearMass
01997     // G4ParticleTable::GetParticleTable()->FindIon(Zm,Am,0,Zm) creates new Ion!
01998     if(A==256 && Z==128) m_value=256000.;
01999     else                 m_value=k+G4NucleiProperties::GetNuclearMass(A,Z);
02000   }
02001   //@@//G4double maxM= k+Z*mP+N*mN+S*mL+eps;      // @@ eps -- Wings of the Mass parabola
02002   //@@//if(m_value>maxM) m_value=maxM;
02003   G4double mm_value=m_value;
02004   if(Sm<0)                                  // For the new algorithm of calculation 
02005   {
02006     if(Nm<0)
02007     {
02008       km_value+=-Nm*mPi;
02009       Zm+=Nm;
02010       Nm=0;
02011     }
02012     if(Zm<0)
02013     {
02014       km_value+=-Zm*mPi;
02015       Nm+=Zm;
02016       Zm=0;
02017     }
02018     G4int Am=Zm+Nm;
02019     if(!Am) return km_value+eps;
02020     mm_value=km_value+Am*um_value;          // Expected mass in atomic units
02021     //G4double Dm=Nm-Zm;                    // Isotopic shift of the nucleus
02022     if ( (Am < 1 && km_value==0.) || Zm < 0 || Nm < 0 ) 
02023     {   // @@ Can be generalized to anti-nuclei
02024 #ifdef debug
02025       G4cerr<<"*G4QPDGCode::CalcNucM:A="<<Am<<"<1 || Z="<<Zm<<"<0 || N="<<Nm<<"<0"<<G4endl;
02026 #endif
02027     }
02028     if     (!Zm) return km_value+Nm*(mN+.1);
02029     else if(!Nm) return km_value+Zm*(mP+1.);
02030     //else if(Nm<=9&&Zm<=9) mm_value+=1.433e-5*pow(double(Zm),2.39)-Zm*me+c[Nm-1][Zm-1];// Geant4
02031     else 
02032     {
02033       // IsInTable is inside G4NucleiProperties::GetNuclearMass
02034       // G4ParticleTable::GetParticleTable()->FindIon(Zm,Am,0,Zm) creates new Ion!
02035       mm_value=km_value+G4NucleiProperties::GetNuclearMass(Am,Zm);
02036     }
02037     //@@//G4double mM= km_value+Zm*mP+Nm*mN+eps;
02038     //@@//if(mm_value>mM) mm_value=mM;
02039   }
02040   if(m_value>mm_value) m_value=mm_value;
02041   if(S>0)
02042   {
02043     G4double bs=0.;
02044     if     (A==2) bs=a2;
02045     else if(A==3) bs=a3;
02046     else if(A>3)  bs=b7*exp(-b8/(A+1.));
02047     m_value+=S*(mL-bs);
02048   }  
02049 #ifdef debug
02050   G4cout<<"G4QPDGCode::CalcNuclMass: >->-> OUT <-<-< m="<<m_value<<G4endl;
02051 #endif
02052   //if(z==8&&n==9&&!s_value) G4cout<<"G4QPDGC::GetNucM:m="<<m_value<<",mm="<<mm_value<<G4endl;
02053   return m_value;
02054 }
02055 
02056 // Get a Quark Content {d,u,s,ad,au,as} for the PDG
02057 G4QContent G4QPDGCode::GetQuarkContent() const
02058 {
02059   G4int            a=0; // particle
02060   if(thePDGCode<0) a=1; // anti-particle
02061   G4int d=0;
02062   G4int u=0;
02063   G4int s_value=0;
02064   G4int ad=0;
02065   G4int au=0;
02066   G4int as=0;
02067   G4int ab=abs(thePDGCode);
02068   if     (!ab)
02069   {
02070     G4cerr<<"***G4QPDGCode::GetQuarkContent: PDG=0, return (0,0,0,0,0,0)"<<G4endl;
02071     return G4QContent(0,0,0,0,0,0);
02072   }
02073   else if(ab<4) // @@ for SU(6) : ab<7
02074   {
02075     if     (thePDGCode== 1) return G4QContent(1,0,0,0,0,0);
02076     else if(thePDGCode== 2) return G4QContent(0,1,0,0,0,0);
02077     else if(thePDGCode== 3) return G4QContent(0,0,1,0,0,0);
02078     else if(thePDGCode==-1) return G4QContent(0,0,0,1,0,0);
02079     else if(thePDGCode==-2) return G4QContent(0,0,0,0,1,0);
02080     else if(thePDGCode==-3) return G4QContent(0,0,0,0,0,1);
02081   }
02082   else if(ab<99)
02083   {
02084     if     (thePDGCode== 11) return G4QContent(1,0,0,0,1,0);
02085     else if(thePDGCode==-11) return G4QContent(0,1,0,1,0,0);
02086     else if(thePDGCode== 13) return G4QContent(1,0,0,0,1,0);
02087     else if(thePDGCode==-13) return G4QContent(0,1,0,1,0,0);
02088     else if(thePDGCode== 15) return G4QContent(1,0,0,0,1,0);
02089     else if(thePDGCode==-15) return G4QContent(0,1,0,1,0,0);
02090 #ifdef debug
02091     if     (ab==22) G4cout<<"-W-G4QPDGC::GetQuarkCont: For the Photon? - Return 0"<<G4endl;
02092     else if(ab==10) G4cout<<"-W-G4QPDGC::GetQuarkCont: For Chipolino? - Return 0"<<G4endl;
02093     else G4cout<<"-W-G4QPDGCode::GetQuarkCont: For PDG="<<thePDGCode<<" Return 0"<<G4endl;
02094 #endif
02095     return G4QContent(0,0,0,0,0,0); // Photon, bosons, leptons
02096   }
02097   else if(ab<80000000) // Baryons & Mesons
02098   {
02099     G4int c=ab/10;     // temporary (quarks)
02100     G4int f=c%10;      // (1) low quark
02101     G4int v=c/10;      // (2,3) temporary(B) or (2) final(M) (high quarks, high quark)
02102     G4int t=0;         // (3)prototype of highest quark (B)
02103 #ifdef sdebug
02104     G4cout<<"G4QPDGCode::GetQuarkContent: a="<<ab<<", c="<<c<<", f="<<f<<", v="<<v<<G4endl;
02105 #endif
02106     if(v>10)           // Baryons
02107     {
02108       t=v/10;          // (3) highest quark
02109       v%=10;           // (2) high quark
02110       if     (f==1)
02111       {
02112         if(a) ad++;
02113         else   d++;
02114       }
02115       else if(f==2)
02116       {
02117         if(a) au++;
02118         else   u++;
02119       }
02120       else if(f==3)
02121       {
02122         if(a) as++;
02123         else   s_value++;
02124       }
02125       else G4cerr<<"*G4QPDGC::GetQCont:1 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
02126       if     (v==1)
02127       {
02128         if(a) ad++;
02129         else   d++;
02130       }
02131       else if(v==2)
02132       {
02133         if(a) au++;
02134         else   u++;
02135       }
02136       else if(v==3)
02137       {
02138         if(a) as++;
02139         else   s_value++;
02140       }
02141       else G4cerr<<"*G4QPDGC::GetQCont:2 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
02142       if     (t==1)
02143       {
02144         if(a) ad++;
02145         else   d++;
02146       }
02147       else if(t==2)
02148       {
02149         if(a) au++;
02150         else   u++;
02151       }
02152       else if(t==3)
02153       {
02154         if(a) as++;
02155         else   s_value++;
02156       }
02157       else G4cerr<<"*G4QPDGC::GetQCont:3 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
02158       return G4QContent(d,u,s_value,ad,au,as);
02159     }
02160     else        // Mesons
02161     {
02162       if(f==v)
02163       {
02164         if     (f==1) return G4QContent(1,0,0,1,0,0);
02165         else if(f==2) return G4QContent(0,1,0,0,1,0);
02166         else if(f==3) return G4QContent(0,0,1,0,0,1);
02167         else G4cerr<<"*G4QPDGC::GetQC:4 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
02168       }
02169       else
02170       {
02171         if     (f==1 && v==2)
02172         {
02173           if(a)return G4QContent(1,0,0,0,1,0);
02174           else return G4QContent(0,1,0,1,0,0);
02175         }
02176         else if(f==1 && v==3)
02177         {
02178           if(a)return G4QContent(0,0,1,1,0,0);
02179           else return G4QContent(1,0,0,0,0,1);
02180         }
02181         else if(f==2 && v==3)
02182         {
02183           if(a)return G4QContent(0,0,1,0,1,0);
02184           else return G4QContent(0,1,0,0,0,1);
02185         }
02186         else G4cerr<<"*G4QPDGC::GetQC:5 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
02187       }
02188     }
02189   }
02190   else                    
02191   {
02192     G4int szn=ab-90000000;
02193     G4int ds=0;
02194     G4int dz=0;
02195     G4int dn=0;
02196     if(szn<-100000)
02197     {
02198       G4int ns_value=(-szn)/1000000+1;
02199       szn+=ns_value*1000000;
02200       ds+=ns_value;
02201     }
02202     else if(szn<-100)
02203     {
02204       G4int nz=(-szn)/1000+1;
02205       szn+=nz*1000;
02206       dz+=nz;
02207     }
02208     else if(szn<0)
02209     {
02210       G4int nn=-szn;
02211       szn=0;
02212       dn+=nn;
02213     }
02214     G4int sz =szn/1000;
02215     G4int n  =szn%1000;
02216     if(n>700)
02217     {
02218       n-=1000;
02219       dz--;
02220     }
02221     G4int z  =sz%1000-dz;
02222     if(z>700)
02223     {
02224       z-=1000;
02225       ds--;
02226     }
02227     s_value  =sz/1000-ds;
02228     G4int b=z+n+s_value;
02229     d=n+b;
02230     u=z+b;
02231     if      (d<0&&u<0&&s_value<0) return G4QContent(0,0,0,-d,-u,-s_value);
02232     else if (u<0&&s_value<0)      return G4QContent(d,0,0,0,-u,-s_value);
02233     else if (d<0&&s_value<0)      return G4QContent(0,u,0,-d,0,-s_value);
02234     else if (d<0&&u<0)            return G4QContent(0,0,s_value,-d,-u,0);
02235     else if (u<0)                 return G4QContent(d,0,s_value,0,-u,0);
02236     else if (s_value<0)           return G4QContent(d,u,0,0,0,-s_value);
02237     else if (d<0)                 return G4QContent(0,u,s_value,-d,0,0);
02238     else                          return G4QContent(d,u,s_value,0,0,0);
02239   }
02240   return G4QContent(0,0,0,0,0,0);
02241 }
02242 
02243 // Quark Content of pseudo-meson for q_i->q_o Quark Exchange
02244 G4QContent G4QPDGCode::GetExQContent(G4int i, G4int o)  const
02245 {
02246   G4QContent cQC(0,0,0,0,0,0);
02247   if     (!i)   cQC.IncAD();
02248   else if(i==1) cQC.IncAU();
02249   else if(i==2) cQC.IncAS();
02250   else G4cerr<<"***G4QPDGCode::GetExQContent: strange entering quark i="<<i<<G4endl;
02251   if     (!o)   cQC.IncD();
02252   else if(o==1) cQC.IncU();
02253   else if(o==2) cQC.IncS();
02254   else G4cerr<<"***G4QPDGCode::GetExQContent: strange exiting quark o="<<o<<G4endl;
02255   return cQC;
02256 }
02257 
02258 // Relative Cross Index for q_i->q_o (q=0(d),1(u),2(s), 7 means there is no such cluster)
02259 G4int G4QPDGCode::GetRelCrossIndex(G4int i, G4int o)  const
02260 {
02261   //                          pD,nD,DD,DD,ppDnnDpDDnDD n, p, L,S-,S+,X-,X0,O-
02262   static const G4int b01[16]={ 7,15, 7, 7, 7, 7, 7, 7, 1, 7, 7,-1, 7, 7, 7, 7};
02263   static const G4int b02[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 2, 7, 7, 7, 7, 7, 7, 7};
02264   static const G4int b10[16]={18, 7, 7, 7, 7, 7, 7, 7, 7,-1, 7, 7,-2, 7, 7, 7}; // Iso+Bary
02265   static const G4int b12[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 1, 7, 7, 7, 7, 7, 7};
02266   static const G4int b20[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,-2, 7,-3, 7,-4, 7};
02267   static const G4int b21[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,-1,-3, 7,-3, 7, 7};
02268   //                          nn,np,pp,Ln,Lp,LL,nS,pS,nnpnppLnnLpnLppLLnLLpnnS
02269   static const G4int d01[16]={ 1, 1, 7, 1, 7, 7,-3, 7, 1, 7, 1, 1, 7, 1, 7,-5};
02270   static const G4int d02[16]={ 3, 3, 7, 2, 7, 7, 7, 7, 3, 3, 3, 3, 7, 7, 7, 7};
02271   static const G4int d10[16]={ 7,-1,-1, 7,-1, 7, 7,-4, 7,-1, 7,-1,-1, 7,-1, 7}; //B=2,3
02272   static const G4int d12[16]={ 7, 2, 2, 7, 1, 7, 7, 7, 2, 2, 7, 2, 2, 7, 7, 7};
02273   static const G4int d20[16]={ 7, 7, 7,-3,-3,-2, 7,-5, 7, 7, 7,-3,-3,-3,-3, 7};
02274   static const G4int d21[16]={ 7, 7, 7,-2,-2,-1,-6, 7, 7, 7,-2,-2, 7,-2,-2, 7};
02275   //                          nn,np,pp,Ln,Lp,nn,np,pp! n, p,nn,np,pp,Ln,Lp
02276   static const G4int m01[15]={ 1, 1, 7, 1, 7, 1, 1, 7, 1, 7, 1, 1, 7, 1, 7};// Multibaryons
02277   static const G4int m02[15]={ 3, 3, 7, 3, 3, 7, 7, 7, 3, 3, 3, 3, 7, 7, 7};// @@ Regular !
02278   static const G4int m10[15]={ 7,-1,-1, 7,-1, 7,-1,-1, 7,-1, 7,-1,-1, 7,-1};// 01->1,10->-1
02279   static const G4int m12[15]={ 7, 2, 2, 2, 2, 7, 7, 7, 2, 2, 7, 2, 2, 7, 7};// 12->2,21->-2
02280   static const G4int m20[15]={ 7, 7, 7,-3,-3, 7,-3,-3, 7, 7, 7,-3,-3,-3,-3};// 02->3,20->-3
02281   static const G4int m21[15]={ 7, 7, 7,-2,-2,-2,-2, 7, 7, 7,-2,-2, 7,-2,-2};
02282   //static const G4int fragmStart = 82;    // "Isonuclei are added & Leptons are added"
02283   static const G4int fragmStart = 45;    // "Isonuclei & Leptons are added, Reduced"
02284 
02285   if(theQCode<fragmStart) return 7;
02286   G4int sub=theQCode-fragmStart;
02287   if ( (sub > 1 && sub < 8) || (sub > 12 && sub <16)) return 7; // SuperIso & SuperStrange 
02288   G4int rel=sub;                         // case of nuclear baryons and isonuclei
02289   if     (sub>31) rel =(sub-32)%15;      // case of heavy fragments (BaryNum>3)
02290   else if(sub>15) rel = sub-16;          // case of nuclear di-baryon & tri-baryons
02291 #ifdef debug
02292   G4cout<<"G4QPDGCode::RelGetCrossIndex:i="<<i<<",o="<<o<<",su="<<sub<<",re="<<rel<<G4endl;
02293 #endif
02294   //if ( (rel > 5 && rel < 9) || (rel > 13 && rel <16)) return 7; // SuperStrange're closed
02295   if     (!i)                            // ==> input quark = 0(d) (d=-1/3)
02296   {
02297     if     (!o)       return 0;          // -> output quark = 0(d) => 0 = the same cluster
02298     else if(o==1)                        // -> output quark = 1(u) (ch--)
02299     {
02300       if     (sub<16) return b01[rel];
02301       else if(sub<32) return d01[rel];
02302       else            return m01[rel];
02303     }
02304     else if(o==2)
02305     {
02306       if     (sub<16) return b02[rel];
02307       else if(sub<32) return d02[rel];
02308       else            return m02[rel];
02309     }
02310   }
02311   else if(i==1)                          // ==> input quark = 1(u) (u=2/3)
02312   {
02313     if     (!o)                          // -> output quark = 0(d) (ch++)
02314     {
02315       if     (sub<16) return b10[rel];
02316       else if(sub<32) return d10[rel];
02317       else            return m10[rel];
02318     }
02319     else if(o==1)     return 0;         // -> output quark = 1(u) => 0 = the same cluster
02320     else if(o==2)
02321     {
02322       if     (sub<16) return b12[rel];
02323       else if(sub<32) return d12[rel];
02324       else            return m12[rel];
02325     }
02326   }
02327   else if(i==2)
02328   {
02329     if     (!o)
02330     {
02331       if     (sub<16) return b20[rel];
02332       else if(sub<32) return d20[rel];
02333       else            return m20[rel];
02334     }
02335     else if(o==1)
02336     {
02337       if     (sub<16) return b21[rel];
02338       else if(sub<32) return d21[rel];
02339       else            return m21[rel];
02340     }
02341     else if(o==2)     return 0;
02342   }
02343   return 7;
02344 }
02345 
02346 // Get number of Combinations for q_i->q_o
02347 G4int G4QPDGCode::GetNumOfComb(G4int i, G4int o) const
02348 {
02349   if(i>-1&&i<3)
02350   {
02351     G4int shiftQ=GetRelCrossIndex(i, o);
02352     G4int sQCode=theQCode;                     // QCode of the parent cluster
02353     if     (shiftQ==7) return 0;               // A parent cluster doesn't exist
02354     else if(!shiftQ) sQCode+=shiftQ;           // Shift QCode of son to QCode of his parent
02355     G4QPDGCode parent;                         // Create a temporary (fake) parent cluster
02356     parent.InitByQCode(sQCode);                // Initialize it by Q-Code
02357     G4QContent parentQC=parent.GetQuarkContent();// Quark Content of the parent cluster
02358     if     (!o)   return parentQC.GetD();
02359     else if(o==1) return parentQC.GetU();
02360     else if(o==2) return parentQC.GetS();
02361     else G4cerr<<"***G4QPDGCode:::GetNumOfComb: strange exiting quark o="<<o<<G4endl;
02362   }
02363   else G4cerr<<"***G4QPDGCode:::GetNumOfComb: strange entering quark i="<<i<<G4endl;
02364   return 0;
02365 }
02366 
02367 // Get a total number of Combinations for q_i
02368 G4int G4QPDGCode::GetTotNumOfComb(G4int i) const
02369 {
02370   G4int tot=0;
02371   if(i>-1&&i<3) for(int j=0; j<3; j++) tot+=GetNumOfComb(i, j);
02372   else G4cerr<<"***G4QPDGCode:::GetTotNumOfComb: strange entering quark i="<<i<<G4endl;
02373   return tot;
02374 }
02375 
02376 // Converts nuclear PDGCode to Z(#of protons), N(#of neutrons), S(#of lambdas) values
02377 void G4QPDGCode::ConvertPDGToZNS(G4int nucPDG, G4int& z, G4int& n, G4int& s_value)
02378 {
02379   if(nucPDG>80000000&&nucPDG<100000000)            // Condition of conversion
02380   {
02381     z=0;
02382     n=0;
02383     s_value=0;
02384     G4int r=nucPDG;
02385     if(r==90000000) return;
02386     G4int cn =r%1000;                              // candidate to #of neutrons
02387     if(cn)
02388     {
02389       if(cn>500) cn-=1000;                         // AntiNeutrons
02390       n=cn;                                        // Increment neutrons
02391       r-=cn;                                       // Subtract them from the residual
02392       if(r==90000000) return;
02393     }
02394     G4int cz =r%1000000;                           // candidate to #of neutrons
02395     if(cz)
02396     {
02397       if(cz>500000) cz-=1000000;                   // AntiProtons
02398       z=cz/1000;                                   // Number of protons
02399       r-=cz;                                       // Subtract them from the residual
02400       if(r==90000000) return;
02401     }
02402     G4int cs =r%10000000;                           // candidate to #of neutrons
02403     if(cs)
02404     {
02405       if(cs>5000000) cs-=10000000;                 // AntiLambda
02406       s_value=cs/1000000;                          // Number of Lambdas
02407     }
02408   }
02409   return;
02410 }
02411 
02412 // Only for irreducable DiQaDiQ! L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2
02413 std::pair<G4int,G4int> G4QPDGCode::MakeTwoBaryons(G4int L1, G4int L2, G4int R1, G4int R2)
02414 {
02415   G4int dl=0;
02416   G4int ul=0;
02417   //G4int sl=0;
02418   if     (L1==1) ++dl;
02419   else if(L1==2) ++ul;
02420   //else if(L1==3) ++sl;
02421   if     (L2==1) ++dl;
02422   else if(L2==2) ++ul;
02423   //else if(L2==3) ++sl;
02424   if     (R1==1) ++dl;
02425   else if(R1==2) ++ul;
02426   //else if(R1==3) ++sl;
02427   if     (R2==1) ++dl;
02428   else if(R2==2) ++ul;
02429   //else if(R2==3) ++sl;
02430   if     (dl==2 && ul==2) return make_pair(1114,2212); // @@ can be (2112,2224)
02431   else if(dl==1 && ul==2) return make_pair(3112,2212);
02432   else if(dl==0 && ul==2) return make_pair(3212,3212); // @@ can be (3312,2212)
02433   else if(dl==2 && ul==1) return make_pair(3222,2112);
02434   else if(dl==1 && ul==1) return make_pair(3312,2112); // @@ can be (3322,2212)
02435   else if(dl==2 && ul==0) return make_pair(3112,3112); // @@ can be (3322,1122)
02436   //#ifdef debug
02437   else G4cout<<"-Warning-G4QPDGCode::MakeTwoBaryons: Irreduceble? L1="<<L1<<",L2="<<L2
02438              <<",R1="<<R1<<",R2="<<R2<<G4endl;
02439   //#endif
02440   return make_pair(2212,2112);                         // @@ Just theMinimum, makeException
02441 }

Generated on Mon May 27 17:49:38 2013 for Geant4 by  doxygen 1.4.7