Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:
trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QEnvironment.cc

    r1228 r1315  
    2828//
    2929//
    30 // $Id: G4QEnvironment.cc,v 1.160 2009/12/03 18:09:07 mkossov Exp $
    31 // GEANT4 tag $Name: geant4-09-03 $
     30// $Id: G4QEnvironment.cc,v 1.166 2010/06/10 16:09:49 mkossov Exp $
     31// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    3232//
    3333//      ---------------- G4QEnvironment ----------------
    … …  
    103103        G4cout<<"*G4QE::Const:iH#"<<ih<<","<<curQH->GetQC()<<curQH->Get4Momentum()<<G4endl;
    104104#endif
     105        if(curQH->GetPDGCode() == 10)        // Chipolino is found in the input -> Decay
     106        {
     107          G4QContent   chQC=curQH->GetQC();  // Quark content of the Hadron-Chipolino
     108          G4QChipolino QCh(chQC);            // Define a Chipolino instance for the Hadron
     109          G4LorentzVector ch4M=curQH->Get4Momentum(); // 4Mom of the Hadron-Chipolino
     110          G4QPDGCode h1QPDG=QCh.GetQPDG1();  // QPDG of the first hadron
     111          G4double   h1M   =h1QPDG.GetMass();// Mass of the first hadron
     112          G4QPDGCode h2QPDG=QCh.GetQPDG2();  // QPDG of the second hadron
     113          G4double   h2M   =h2QPDG.GetMass();// Mass of the second hadron
     114          G4double   chM2  =ch4M.m2();       // Squared Mass of the Chipolino
     115          if( sqr(h1M+h2M) < chM2 )          // Decay is possible
     116          {
     117            G4LorentzVector h14M(0.,0.,0.,h1M);
     118            G4LorentzVector h24M(0.,0.,0.,h2M);
     119            if(!G4QHadron(ch4M).DecayIn2(h14M,h24M))
     120            {
     121              G4cerr<<"***G4QE::Constr: CM="<<std::sqrt(chM2)<<" -> h1="<<h1QPDG<<"("<<h1M
     122                    <<") + h2="<<h1QPDG<<"("<<h2M<<") = "<<h1M+h2M<<" **Failed**"<<G4endl;
     123              throw G4QException("**G4QEnvironment::Constr:QChip DecIn2 error");
     124            }
     125            delete curQH;                    // Kill the primary Chipolino
     126            G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
     127            theQHadrons.push_back(h1H);      // (delete equivalent)
     128            curQH    = new G4QHadron(h1H);   // ... just to remember independently
     129            theProjectiles.push_back(curQH); // Remember it for the error message
     130
     131#ifdef debug
     132            G4cout<<"G4QE::Constr: QChipolino -> H1="<<h1QPDG<<h14M<<G4endl;
     133#endif
     134            curQH = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
     135            theQHadrons.push_back(curQH);    // (delete equivalent)
     136#ifdef debug
     137            G4cout<<"G4QE::Constr: QChipolino -> H2="<<h2QPDG<<h24M<<G4endl;
     138#endif
     139          }
     140          else
     141          {
     142            G4cerr<<"***G4QE::Const:iH#"<<ih<<","<<curQH->GetQC()<<curQH->Get4Momentum()
     143                  <<", chipoM="<<std::sqrt(chM2)<<" < m1="<<h1M<<"("<<h1QPDG<<") + m2="
     144                  <<h2M<<"("<<h2QPDG<<") = "<<h1M+h2M<<G4endl;
     145            throw G4QException("G4QEnvironment::HadronizeQEnv: LowMassChipolino in Input");
     146          }
     147        }
    105148        theQHadrons.push_back(curQH);        // (delete equivalent)
    106         curQH    = new G4QHadron(projHadrons[ih]); // ... to remember
    107         theProjectiles.push_back(curQH);      // Remenber it for the error message
     149        curQH    = new G4QHadron(curQH);     // ... just to remember independently
     150        theProjectiles.push_back(curQH);     // Remenber it for the error message
    108151      }
    109152    }
    … …  
    116159      theQHadrons.push_back(curQH);          // (delete equivalent)
    117160    }
    118     if (nHadrons<0) G4cout<<"**G4QE::Const:NH="<<nHadrons<<" < 0 !"<<G4endl;
     161    if (nHadrons<0) G4cout<<"***Warning****G4QE::Const:NH="<<nHadrons<<" < 0 !"<<G4endl;
    119162    return;
    120163  }
    … …  
    13901433  static const G4QPDGCode s0QPDG(3122);
    13911434  static const G4double mPi0 = G4QPDGCode(111).GetMass();
     1435  //static const G4double dPi0 = mPi0+mPi0;
     1436  static const G4double fPi0 = 4*mPi0;
    13921437  static const G4double mPi  = G4QPDGCode(211).GetMass();
    13931438  static const G4double mK   = G4QPDGCode(321).GetMass();
    … …  
    17221767      totInC            += pQ->GetQC().GetCharge();
    17231768    } // End of TotInitial4Momentum summation LOOP over Quasmons
    1724     G4int nsHadr  = theQHadrons.size();        // Update the value of OUTPUT entries
     1769    G4int nsHadr  = theQHadrons.size();        // Update the value of #of OUTPUT entries
    17251770    if(nsHadr) for(G4int jso=0; jso<nsHadr; jso++)// LOOP over output hadrons
    17261771    {
    … …  
    17341779    }
    17351780#endif
    1736     //G4int   c3Max = 27;
     1781    // @@ Experimental calculations
     1782    G4QContent totInQC=theEnvironment.GetQCZNS();
     1783    G4LorentzVector totIn4M=theEnvironment.Get4Momentum();
     1784    for (G4int is=0; is<nQuasmons; is++) // Sum 4mom's of Quasmons for comparison
     1785    {
     1786      G4Quasmon*      pQ = theQuasmons[is];
     1787      totIn4M           += pQ->Get4Momentum();
     1788      totInQC           += pQ->GetQC();
     1789    } // End of TotInitial4Momentum/QC summation LOOP over Quasmons
     1790    G4double totMass=totIn4M.m();
     1791    G4QNucleus totN(totInQC);
     1792    G4double totM=totN.GetMZNS();
     1793    G4double excE = totMass-totM;
     1794    // @@ End of experimental calculations
     1795    G4int   envA=theEnvironment.GetA();
     1796    //G4int   c3Max = 27;                   // Big number (and any #0) slowes dow a lot
    17371797    //G4int   c3Max = 9;                    // Max#of "no hadrons" steps (reduced below?)
    1738     //G4int   c3Max = 3;
     1798    G4int   c3Max = 3;
     1799    if(excE > fPi0) c3Max=(G4int)(excE/mPi0); // Try more for big excess
    17391800    //G4int   c3Max = 1;
    1740     G4int   c3Max = 0;
     1801    //if(excE > dPi0) c3Max=(G4int)(excE/mPi0); // Try more for big excess
     1802    //G4int   c3Max = 0;
     1803    //if(excE > mPi0) c3Max=(G4int)(excE/mPi0); // Try more for big excess
     1804    //G4int   c3Max = 0;                    // It closes the force decay of Quasmon at all!
     1805    //
    17411806    //G4int   premC = 27;
    17421807    //G4int   premC = 3;
    17431808    G4int   premC = 1;
    1744     G4int   envA=theEnvironment.GetA();
    17451809    //if(envA>1&&envA<13) premC = 24/envA;
    1746     //if(envA>1&&envA<19) premC = 36/envA;
    1747     //if(envA>1&&envA<25) premC = 48/envA;
    1748     //if(envA>1&&envA<41) premC = 80/envA;
     1810    if(envA>1&&envA<19) premC = 36/envA;
     1811    //if(envA>1&&envA<31) premC = 60/envA;
     1812    //if(envA>1&&envA<61) premC = 120/envA;
    17491813    G4int  sumstat= 2;                     // Sum of statuses of all Quasmons
    17501814    G4bool force  = false;                 // Prototype of the Force Major Flag
    17511815    G4int cbR     =0;                      // Counter of the "Stoped by Coulomb Barrier"
    17521816    //
    1753     G4int cbRM    =0;                      //** MaxCounter of "StopedByCoulombBarrier" **
     1817    //G4int cbRM    =0;                      //** MaxCounter of "StopedByCoulombBarrier" **
    17541818    //G4int cbRM    =1;                      // MaxCounter of the "StopedByCoulombBarrier"
    17551819    //G4int cbRM    =3;                      // MaxCounter of the "StopedByCoulombBarrier"
    17561820    //G4int cbRM    =9;                      // MaxCounter of the "Stoped byCoulombBarrier"
     1821    G4int cbRM    = c3Max;                 // MaxCounter of the "StopedByCoulombBarrier"
    17571822    G4int totC    = 0;                     // Counter to break the "infinit" loop
    17581823    //G4int totCM   = 227;                   // Limit for the "infinit" loop counter
    17591824    G4int totCM   = envA;                   // Limit for the "infinit" loop counter
    17601825    //G4int totCM   = 27;                    // Limit for this counter
    1761     G4int nCnMax = 1;                      // MaxCounterOfHadrFolts for shortCutSolutions
     1826    //G4int nCnMax = 1;                      // MaxCounterOfHadrFolts for shortCutSolutions
    17621827    //G4int nCnMax = 3;                      // MaxCounterOfHadrFolts for shortCutSolutions
    17631828    //G4int nCnMax = 9;                      // MaxCounterOfHadrFolts for shortCutSolutions
     1829    G4int nCnMax = c3Max;                  // MaxCounterOfHadrFolts for shortCutSolutions
    17641830    G4bool first=true;                     // Flag of the first interaction (only NucMedia)
    17651831    G4int cAN=0;                           // Counter of the nucleon absorptions
    … …  
    23892455                        {                     //                                          ^
    23902456                          G4cerr<<"***G4QE::HQE:tM="<<ttM<<"->h1="<<h1QPDG<<"("<<h1M //   ^
    2391                                 <<")+h2="<<h1QPDG<<"("<<h2M<<"="<<h1M+h2M<<G4endl;   //   ^
     2457                                <<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;  //   ^
    23922458                          throw G4QException("**G4QE::HadrQEnv:QChip (1) DecIn2 error");//^
    23932459                        }                     //                                          ^
    … …  
    40534119
    40544120//Evaporate Residual Nucleus
    4055 void G4QEnvironment::EvaporateResidual(G4QHadron* qH)
    4056 {//  ================================================
     4121void G4QEnvironment::EvaporateResidual(G4QHadron* qH,  G4bool fCGS)
     4122{//  ==============================================================
    40574123  static const G4double mAlph = G4QPDGCode(2112).GetNuclMass(2,2,0);
    40584124  static const G4double mDeut = G4QPDGCode(2112).GetNuclMass(1,1,0);
    … …  
    41024168  if(!thePDG) thePDG = theQC.GetSPDGCode();  // If there is no PDG code, get it from QC
    41034169  if(thePDG==10&&theQC.GetBaryonNumber()>0) thePDG=theQC.GetZNSPDGCode();
    4104   if(theS>0) thePDG-=theS*999999;            // @@ May hide hypernuclear problems (G4)
     4170  if(theS>0) thePDG-=theS*999999;            // @@ May hide hypernuclear problems (G4) ! @@
    41054171  G4double totGSM = G4QNucleus(thePDG).GetGSMass();// TheGroundStMass of theTotalResNucleus
    41064172  if(theBN==2)
    … …  
    41364202#endif
    41374203    G4Quasmon* quasH = new G4Quasmon(theQC,q4M);
    4138     if(!CheckGroundState(quasH,true))
     4204    if(fCGS && !CheckGroundState(quasH, true) )
    41394205    {
    41404206#ifdef debug
    … …  
    43004366  if(ExCount>=MaxExCnt)
    43014367  {
    4302     G4cerr<<"*G4QEnv::Fragment:Exception.Target="<<theTargetPDG<<". Projectiles:"<<G4endl;
    43034368    G4int nProj=theProjectiles.size();
     4369    G4cerr<<"*G4QEnv::Fragment:Exception.Target="<<theTargetPDG<<". #Proj="<<nProj<<G4endl;
    43044370    if(nProj) for(G4int ipr=0; ipr<nProj; ipr++)
    43054371    {
    … …  
    43084374    }
    43094375    throw G4QException("G4QEnvironment::Fragment:This reaction caused the CHIPSException");
    4310     //G4Exception("G4QEnvironment::Fragment","027",FatalException,"GeneralCHIPSException");
    43114376  }
    43124377  // Put the postponed hadrons in the begining of theFragments and clean them up
    … …  
    43694434  static const G4double mXiZ = G4QPDGCode(3322).GetMass();
    43704435  static const G4double mXiM = G4QPDGCode(3312).GetMass();
     4436  static const G4double mOmM = G4QPDGCode(3334).GetMass();
    43714437#ifdef debug
    43724438  static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);
    … …  
    70067072      G4int nSM=0;                             // A#0f unavoidable Sigma-
    70077073      G4int nSP=0;                             // A#0f unavoidable Sigma+
    7008       G4int nXZ=0;                             // A#0f unavoidable Xi-
    7009       G4int nXM=0;                             // A#0f unavoidable Xi0
    7010       //G4int nOM=0;                             // A#0f unavoidable Omega-
     7074      G4int nXZ=0;                             // A#0f unavoidable Xi0
     7075      G4int nXM=0;                             // A#0f unavoidable Xi-
     7076      G4int nOM=0;                             // A#0f unavoidable Omega-
    70117077      // @@ Now it does not include more complicated clusters as Omega-,Xi- (Improve)
    70127078      if(nC < 0)                               // Negative hypernucleus
    … …  
    70347100              nL-=nX-nC;
    70357101            }
     7102            // @@ Don't close this warning as it costs days to find out this growing point!
    70367103            else G4cout<<"-Warning-G4QEn::FSI:*Improve*,-nC<=nL,nL>nB, PDG="<<hPDG<<G4endl;
    70377104          }
    … …  
    70687135            }
    70697136          }
    7070 #ifdef pdebug
     7137          else if(nL > nB && ((nL-nB)/2)*2 == nL-nB && nC+nL <= nB+nB)// Omega- can be used
     7138          {
     7139            nOM=(nL-nB)/2;                     // This is how many Omega- can be made
     7140            nL=nB+nB-nL-nC;                    // a#of Lambdas should be 0 !
     7141            nSP=nOM-nC;
     7142          }
     7143          // @@ Do not close this warning as it costs days to find out this growing point !
    70717144          else G4cout<<"-Warning-G4QEn::FSI:*Improve*,nC>nB-nL,nC<=nB, PDG="<<hPDG<<G4endl;
    7072 #endif
    70737145        }
    70747146        else
    … …  
    70817153      G4LorentzVector r4M=curHadr->Get4Momentum(); // Real 4-momentum of the hypernucleus
    70827154      G4double reM=r4M.m();                    // Real mass of the hypernucleus
    7083       // Subtract Lamb/Sig/Xi from the Nucleus and decay
     7155      // Subtract Lamb/Sig/Xi/Omega from the Nucleus and decay
    70847156      G4int SS=nL+nSP+nSM+nXZ+nXM;
    7085       if(SS<nB)                                // Should not be Xi's in the nucleus
     7157      if(SS<nB && !nOM)                        // Should not be Xi's or Omeg in the nucleus
    70867158      {
    70877159        G4int rlPDG = hPDG-nL*1000000-nSP*1000999-nSM*999001;
    … …  
    70917163        G4cout<<"G4QEnvironment::FSI:*G4* nS+="<<nSP<<", nS-="<<nSM<<", nL="<<nL<<G4endl;
    70927164#endif
    7093         if(nSP||nSM)
     7165        if(nSP || nSM)
    70947166        {
    70957167          if(nL)
    … …  
    72757347        else // If this Error shows up (lowProbable appearance) => now it is left as is
    72767348        {
    7277           G4double d=rlM+MLa-reM;
    7278           G4cerr<<"G4QE::F:R="<<rlM<<",S+="<<nSP<<",S-="<<nSM<<",L="<<nL<<",d="<<d<<G4endl;
    7279           d=rnM+mPi0-reM;
    7280           G4cerr<<"-W-G4QE::FSI:HypN="<<hPDG<<", M="<<reM<<"<"<<rnM+mPi0<<",d="<<d<<G4endl;
     7349          G4cout<<"-Warning-G4QEnv::F:R="<<rlM<<",S+="<<nSP<<",S-="<<nSM<<",L="<<nL<<",d1="
     7350                <<rlM+MLa-reM<<G4endl;
     7351          G4cout<<"-Warning-G4QEnv::FSI:HyperN="<<hPDG<<", M="<<reM<<"<"<<rnM+mPi0<<", d2="
     7352                <<rnM+mPi0-reM<<G4endl;
    72817353          //throw G4QException("G4QEnvironment::FSInteract: Hypernuclear conversion");
    72827354        }
    72837355#endif
    72847356      }
    7285       else if(SS==nB)
     7357      else if(SS==nB && !nOM) // Decay with Xi, but without residual nucleus
    72867358      {
    72877359#ifdef pdebug
    … …  
    74447516            fM*=fS;
    74457517            sM+=sS;
    7446             if(reM>fM+sM-eps)                 // can be split or decayed
     7518            G4double mm=fM+sM;
     7519            G4double MM=reM+eps;
     7520            if(MM<=mm && fPDG==3122 && sPDG==3222) // Lamba,Sigma+ => Xi0,p (can be 50%)
     7521            {
     7522              fPDG= 2212;
     7523              fM  = mProt;
     7524              sPDG= 3322;
     7525              sM  = mXiZ;
     7526              mm  = fM+sM;
     7527            }
     7528            if(MM>mm)                         // can be split or decayed
    74477529            {
    74487530              G4LorentzVector f4M(0.,0.,0.,fM);
    74497531              G4LorentzVector s4M(0.,0.,0.,sM);
    74507532              G4double sum=fM+sM;
    7451               if(fabs(reM-sum)<eps)           // splitting
     7533              if(fabs(reM-sum)<=eps)           // splitting
    74527534              {
    74537535                f4M=r4M*(fM/sum);
    … …  
    74817563              }
    74827564            }
    7483             else G4cout<<"-Warning-G4QE::FSI:*Improve* S2, PDG="<<hPDG<<",M="<<reM<<G4endl;
    7484             //else if(reM>rnM+mPi0-eps&&!nSP&&!nSM)  // Lambda->N only if Sigmas are absent
    7485             //{
    7486             //  G4int nPi=static_cast<G4int>((reM-rnM)/mPi0);
    7487             //  if (nPi>nL) nPi=nL;
    7488             //  G4double npiM=nPi*mPi0;
    7489             //  G4LorentzVector n4M(0.,0.,0.,rnM);
    7490             //  G4LorentzVector h4M(0.,0.,0.,npiM);
    7491             //  G4double sum=rnM+npiM;
    7492             //  if(fabs(reM-sum)<eps)
    7493             //  {
    7494             //    n4M=r4M*(rnM/sum);
    7495             //    h4M=r4M*(npiM/sum);
    7496             //  }
    7497             //  else if(reM<sum || !G4QHadron(r4M).DecayIn2(n4M,h4M))
    7498             //  {
    7499             //    G4cerr<<"**G4QE::FSI:HypN,M="<<reM<<"<A+n*Pi0="<<sum<<",d="<<sum-reM<<G4endl;
    7500             //    throw G4QException("***G4QEnvironment::FSInter: Hypernuclear decay error");
    7501             //  }
    7502             //  curHadr->Set4Momentum(n4M);
    7503             //  curHadr->SetQPDG(G4QPDGCode(rnPDG)); // converted hypernucleus
     7565            // @@ Don't close this warning as it costs days to find out this growing point!
     7566            else G4cout<<"-Warning-G4QE::FSI:*Improve* S2, PDG="<<hPDG<<",M="<<reM<<",fS="
     7567                       <<fS<<",sS="<<sS<<",fM="<<fM<<",sM="<<sM<<G4endl;
     7568          }
     7569          // @@ Do not close this warning as it costs days to find out this growing point !
     7570          else G4cout<<"-Warning-G4QE::FSI:*Improve* 3Body's necessary,PDG="<<hPDG<<G4endl;
     7571        }
     7572        // @@ Do not close this warning as it costs days to find out this growing point !
     7573        else G4cout<<"-Warning-G4QE::FSI:*Improve* MultyD, SP="<<SP<<",PDG="<<hPDG<<G4endl;
     7574      }
     7575      else if(nOM) // Decay with Omega-, but without residual nucleus
     7576      {
    75047577#ifdef pdebug
    7505             //  G4cout<<"*G4QE::FSI:HN="<<r4M<<"->A="<<rnPDG<<n4M<<",n*Pi0="<<nPi<<h4M
    7506             //        <<G4endl;
    7507 #endif
    7508             //  if(nPi>1) h4M=h4M/nPi;               // split the 4-mom if necessary
    7509             //  for(G4int ihn=0; ihn<nPi; ihn++)     // loop over additional pions
    7510             //  {
    7511             //    G4QHadron* thePion = new G4QHadron(111,h4M);// Make a New Hadr for the pion
    7512             //    theQHadrons.push_back(thePion);    // (del.eq.- user is responsible for del)
    7513             //    //theFragments->push_back(thePion);  // (delete equivalent for the pion)
    7514             //  }
    7515             //  if(rnPDG==90000002)                  // Additional action with curHadr change
    7516             //  {
    7517             //    G4LorentzVector newLV=n4M/2.;
    7518             //    curHadr->Set4Momentum(newLV);
    7519             //    curHadr->SetQPDG(G4QPDGCode(90000001));
    7520             //    G4QHadron* secHadr = new G4QHadron(curHadr);
    7521             //    theQHadrons.push_back(secHadr);    // (del.eq.- user is responsible for del)
    7522             //  }
    7523             //  else if(rnPDG==90002000)             // Additional action with curHadr change
    7524             //  {
    7525             //    G4LorentzVector newLV=n4M/2.;
    7526             //    curHadr->Set4Momentum(newLV);
    7527             //    curHadr->SetQPDG(G4QPDGCode(90001000));
    7528             //    G4QHadron* secHadr = new G4QHadron(curHadr);
    7529             //    theQHadrons.push_back(secHadr);    // (del.eq.- user is responsible for del)
    7530             //  }
    7531             //  // @@ Add multybaryon decays if necessary
    7532             //}
    7533             //else if(reM>rnM-eps)    // Decay in nonstange and gamma
    7534             //{
    7535             //  G4LorentzVector n4M(0.,0.,0.,rnM);
    7536             //  G4LorentzVector h4M(0.,0.,0.,0.);
    7537             //  G4double sum=rnM;
    7538             //  if(fabs(reM-sum)<eps) n4M=r4M;
    7539             //  else if(reM<sum || !G4QHadron(r4M).DecayIn2(n4M,h4M))
    7540             //  {
    7541             //    G4cerr<<"**G4QE::FSI:HypN,M="<<reM<<"<A+n*Pi0="<<sum<<",d="<<sum-reM<<G4endl;
    7542             //    throw G4QException("***G4QEnvironment::FSInt:Hypernuclear GammaDecay error");
    7543             //  }
    7544             //  curHadr->Set4Momentum(n4M);
    7545             //  curHadr->SetQPDG(G4QPDGCode(rnPDG)); // converted hypernucleus
     7578        G4cout<<"G4QEnvir::FSI:Omega-,B="<<nB<<",SP="<<nSP<<",OM="<<nOM<<",L="<<nL<<G4endl;
     7579#endif
     7580        G4int SP=0;
     7581        if(nL) ++SP;
     7582        if(nSP) ++SP;
     7583        if(nOM) ++SP;
     7584        if(!nL)
     7585        {
     7586          G4int    fPDG= 3334;
     7587          G4double fM  = mOmM;
     7588          G4int    fS  = nOM;
     7589          G4int    sPDG= 3222;
     7590          G4double sM  = mSigP;
     7591          G4int    sS  = nSP;
    75467592#ifdef pdebug
    7547             //  G4cout<<"*G4QE::FSI:HyN="<<r4M<<"->A="<<rnPDG<<n4M<<",gamma="<<h4M<<G4endl;
    7548 #endif
    7549             //  G4QHadron* theGamma = new G4QHadron(22,h4M);// Make a New Hadr for the gamma
    7550             //  theQHadrons.push_back(theGamma);     // (del.eq.- user is responsible for del)
    7551             //  if(rnPDG==90000002)                  // Additional action with curHadr change
    7552             //  {
    7553             //    G4LorentzVector newLV=n4M/2.;
    7554             //    curHadr->Set4Momentum(newLV);
    7555             //    curHadr->SetQPDG(G4QPDGCode(90000001));
    7556             //    G4QHadron* secHadr = new G4QHadron(curHadr);
    7557             //    theQHadrons.push_back(secHadr);    // (del.eq.- user is responsible for del)
    7558             //  }
    7559             //  else if(rnPDG==90002000)             // Additional action with curHadr change
    7560             //  {
    7561             //    G4LorentzVector newLV=n4M/2.;
    7562             //    curHadr->Set4Momentum(newLV);
    7563             //    curHadr->SetQPDG(G4QPDGCode(90001000));
    7564             //    G4QHadron* secHadr = new G4QHadron(curHadr);
    7565             //    theQHadrons.push_back(secHadr);    // (del.eq.- user is responsible for del)
    7566             //  }
    7567             //  // @@ Add multybaryon decays if necessary
    7568             //}
    7569             //else // If this Error shows up (lowProbable appearance) => now it is left as is
    7570             //{
    7571             //  G4double d=rlM+MLa-reM;
    7572             //G4cerr<<"G4QE::F:R="<<rlM<<",S+="<<nSP<<",S-="<<nSM<<",L="<<nL<<",d="<<d<<G4endl;
    7573             //  d=rnM+mPi0-reM;
    7574             //  G4cerr<<"-W-G4QE::FSI:HyN="<<hPDG<<",M="<<reM<<"<"<<rnM+mPi0<<",d="<<d<<G4endl;
    7575             //  //throw G4QException("G4QEnvironment::FSInteract: Hypernuclear conversion");
    7576             //}
    7577           } // End of decay in 2
    7578         }
    7579       }
     7593          G4cout<<"G4QEnvironment::FSI:2HypD,fS="<<fS<<",fPDG="<<sPDG<<",fM="<<fM<<",sS="
     7594                <<sS<<",sPDG="<<sPDG<<",sM="<<sM<<",rM="<<reM<<G4endl;
     7595#endif
     7596          fM*=fS;
     7597          sM+=sS;
     7598          if(reM>fM+sM-eps)                 // can be split or decayed
     7599          {
     7600            G4LorentzVector f4M(0.,0.,0.,fM);
     7601            G4LorentzVector s4M(0.,0.,0.,sM);
     7602            G4double sum=fM+sM;
     7603            if(fabs(reM-sum)<eps)           // splitting
     7604            {
     7605              f4M=r4M*(fM/sum);
     7606              s4M=r4M*(sM/sum);
     7607            }
     7608            else if(reM<sum || !G4QHadron(r4M).DecayIn2(f4M,s4M))
     7609            {
     7610              G4cerr<<"***G4QE::FSI:Hyp,M="<<reM<<"<OM+SP="<<sum<<",d="<<sum-reM<<G4endl;
     7611              throw G4QException("***G4QEnvironment::FSInter: HypernucOnlyStran3 error");
     7612            }
    75807613#ifdef pdebug
     7614            G4cout<<"G4QEnv::FSI:OHNuc="<<r4M<<hPDG<<"->First="<<fS<<f4M<<fPDG<<",Second="
     7615                  <<sS<<s4M<<sPDG<<G4endl;
     7616#endif
     7617            if(fS>1)
     7618            {
     7619              f4M/=fS;                           // split the Hyperons 4-mom if necessary
     7620              for(G4int il=1; il<fS; il++)       // loop over excessive Hyperons
     7621              {
     7622                G4QHadron* theHyp = new G4QHadron(fPDG,f4M);// Make NewHadr for theHyper
     7623                theQHadrons.push_back(theHyp);    // (del.eq.- user is respons for del)
     7624              }
     7625            }
     7626            curHadr->Set4Momentum(f4M);
     7627            curHadr->SetQPDG(G4QPDGCode(fPDG));  // converted hypernucleus
     7628            if(sS>1) s4M/=sS;                    // split the Hyperon 4-mom if necessary
     7629            for(G4int il=0; il<sS; il++)         // loop over excessive
     7630            {
     7631              G4QHadron* theHyp = new G4QHadron(sPDG,s4M);// Make NewHadr for theHyperon
     7632              theQHadrons.push_back(theHyp);     // (del.eq.- user is respons for del)
     7633            }
     7634          }
     7635          // @@ Don't close this warning as it costs days to find out this growing point!
     7636          else G4cout<<"-Warning-G4QE::FSI:*Improve* OMSP, PDG="<<hPDG<<",M="<<reM<<",nSP="
     7637                     <<nSP<<",nOM="<<nOM<<",mOM="<<fM<<",mSP="<<sM<<G4endl;
     7638        }
     7639        // @@ Do not close this warning as it costs days to find out this growing point !
     7640        else G4cout<<"-Warning-G4QE::FSI:*Improve* nLamb="<<nL<<" # 0, PDG="<<hPDG<<G4endl;
     7641      }
     7642      // @@ Do not close this warning as it costs days to find out this growing point !
    75817643      else G4cout<<"-Warning-G4QE::FSI:*Improve*,S="<<SS<<">B="<<nB<<",PDG="<<hPDG<<G4endl;
    7582 #endif
    75837644    }
    7584     //unsigned nHd=theQHadrons.size()-1;
    7585     //if(hd==nHd && !fOK)
    7586     //{
    7587     //  G4cout<<"---Warning---G4QE::FSI: The Last #"<<hd<<" hadron must be filled"<<G4endl;
    7588     //  fOK=true;
    7589     //  curHadr = new G4QHadron(theQHadrons[hd]);
    7590     //}
    7591     //if(fOK) theFragments->push_back(curHadr);  //(del eq. - user is responsible for del)
    7592     //else hd--;
    75937645    if(!fOK) hd--;
    75947646  } //
    … …  
    86058657          G4cout<<"G4QE::CGS: Fill Quasm "<<valQ<<quas4M<<" in any form"<<G4endl;
    86068658#endif
    8607           EvaporateResidual(quasH);     // Try to evaporate residual (del.eq.)
     8659          EvaporateResidual(quasH, false);   // Try to evaporate residual quasmon (del.eq.)
    86088660#ifdef debug
    86098661          G4cout<<"G4QE::CGS: Fill envir "<<theEQC<<enva4M<<" in any form"<<G4endl;
    … …  
    86118663          // @@ Substitute by EvaporateResidual (if it is not used in the evaporateResid)
    86128664          envaH->Set4Momentum(enva4M);
    8613           EvaporateResidual(envaH);     // Try to evaporate residual (del.eq.)
     8665          EvaporateResidual(envaH);     // Try to evaporate residual environment(del.eq.)
    86148666          // Kill environment as it is already included in the decay
    86158667          theEnvironment=G4QNucleus(G4QContent(0,0,0,0,0,0), G4LorentzVector(0.,0.,0.,0.));
    … …  
    86388690          G4cout<<"G4QE::CGS: fill newQH "<<valQ<<quas4M<<quas4M.m()<<" inAnyForm"<<G4endl;
    86398691#endif
    8640           EvaporateResidual(quasH);     // Try to evaporate residual (del.eq.)
     8692          EvaporateResidual(quasH, false);   // Try to evaporate residual quasmon (del.eq.)
    86418693        }
    86428694        else
    … …  
    86898741            G4cout<<"G4QE::CGS: Fill nucleus "<<reTQC<<nuc4M<<" in any form"<<G4endl;
    86908742#endif
    8691             EvaporateResidual(nucH);     // Try to evaporate residual (del.eq.)
     8743            EvaporateResidual(nucH, false);   // Try to evaporate residual nucleus(del.eq.)
    86928744          }
    86938745        }
    … …  
    87498801                G4cout<<"G4QE::CGS: Fill Resid "<<reTQC<<quas4M<<" in any form"<<G4endl;
    87508802#endif
    8751                 EvaporateResidual(rqH);     // Try to evaporate residual (del.eq.)
     8803                EvaporateResidual(rqH, false);//Try to evaporate residual quasmon (del.eq.)
    87528804                if(envPDG!=90000000) theEnvironment=G4QNucleus(G4QContent(0,0,0,0,0,0),
    87538805                                                          G4LorentzVector(0.,0.,0.,0.));
    … …  
    89719023                              <<ipiQC<<ipi4M<<G4endl;
    89729024#endif
    8973                       EvaporateResidual(rqH); // Try to evaporate residual (del.eq.)
     9025                      EvaporateResidual(rqH, false); // Try to evaporate residual (del.eq.)
    89749026                      if(envPDG!=90000000)theEnvironment=G4QNucleus(G4QContent(0,0,0,0,0,0)
    89759027                                                        ,G4LorentzVector(0.,0.,0.,0.));
    … …  
    90599111                          G4cout<<"G4QE::CGS:FilFr "<<nnQPDG<<quas4M<<" inAnyForm"<<G4endl;
    90609112#endif
    9061                           EvaporateResidual(rqH); // Try to evaporate residual (del.eq.)
     9113                          EvaporateResidual(rqH, false); // Try to evaporate resQ (del.eq.)
    90629114                          if(envPDG!=90000000) theEnvironment=
    90639115                          G4QNucleus(G4QContent(0,0,0,0,0,0),G4LorentzVector(0.,0.,0.,0.));
    … …  
    91179169                            G4cout<<"G4QE::CGS:FilTC "<<tcQPDG<<tc4M<<" inAnyForm"<<G4endl;
    91189170#endif
    9119                             EvaporateResidual(tcH); // Try to evaporate(d.e.)
     9171                            EvaporateResidual(tcH, false);// Try to evaporate hadron (d.e.)
    91209172                          }
    91219173                        }
    … …  
    91589210                      G4cout<<"G4QE::CGS:Fill GSRes "<<reTQC<<quas4M<<" inAnyForm"<<G4endl;
    91599211#endif
    9160                       EvaporateResidual(rqH); // Try to evaporate(d.e.)
     9212                      EvaporateResidual(rqH, false); // Try to evaporate residHadron (d.e.)
    91619213                      if(envPDG!=90000000)theEnvironment=G4QNucleus(G4QContent(0,0,0,0,0,0)
    91629214                                                            ,G4LorentzVector(0.,0.,0.,0.));

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QIsotope.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4QIsotope.cc,v 1.13 2009/08/28 14:49:10 mkossov Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4QIsotope.cc,v 1.16 2010/05/28 15:03:46 mkossov Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//      ---------------- G4QIsotope class ----------------
    … …  
    8686#endif
    8787  // If an error is found in this initialization, please, correct the simple tree below
    88   vector<pair<G4int,G4double> >*a1=new vector<pair<G4int,G4double> >;
     88  vector<pair<G4int,G4double> >*a1=new vector<pair<G4int,G4double> >; // 1-H
    8989  a1->push_back(make_pair(0,.99985));
    9090  a1->push_back(make_pair(1,1.));
    9191  natEl.push_back(a1);
    92   vector<pair<G4int,G4double> >*a2=new vector<pair<G4int,G4double> >;
     92  vector<pair<G4int,G4double> >*a2=new vector<pair<G4int,G4double> >; // 2-He
    9393  a2->push_back(make_pair(2,.999999863));
    9494  a2->push_back(make_pair(1,1.));
    9595  natEl.push_back(a2);
    96   vector<pair<G4int,G4double> >*a3=new vector<pair<G4int,G4double> >;
     96  vector<pair<G4int,G4double> >*a3=new vector<pair<G4int,G4double> >; // 3-Li
    9797  a3->push_back(make_pair(4,.925));
    9898  a3->push_back(make_pair(3,1.));
    9999  natEl.push_back(a3);
    100   vector<pair<G4int,G4double> >*a4=new vector<pair<G4int,G4double> >;
     100  vector<pair<G4int,G4double> >*a4=new vector<pair<G4int,G4double> >; // 4-Be
    101101  a4->push_back(make_pair(5,1.));
    102102  natEl.push_back(a4);
    103   vector<pair<G4int,G4double> >*a5=new vector<pair<G4int,G4double> >;
     103  vector<pair<G4int,G4double> >*a5=new vector<pair<G4int,G4double> >; // 5-B
    104104  a5->push_back(make_pair(6,.801));
    105105  a5->push_back(make_pair(5,1.));
    106106  natEl.push_back(a5);
    107   vector<pair<G4int,G4double> >*a6=new vector<pair<G4int,G4double> >;
     107  vector<pair<G4int,G4double> >*a6=new vector<pair<G4int,G4double> >; // 6-C
    108108  a6->push_back(make_pair(6,.989));
    109109  a6->push_back(make_pair(7,1.));
    110110  natEl.push_back(a6);
    111   vector<pair<G4int,G4double> >*a7=new vector<pair<G4int,G4double> >;
     111  vector<pair<G4int,G4double> >*a7=new vector<pair<G4int,G4double> >; // 7-N
    112112  a7->push_back(make_pair(7,.9963));
    113113  a7->push_back(make_pair(8,1.));
    114114  natEl.push_back(a7);
    115   vector<pair<G4int,G4double> >*a8=new vector<pair<G4int,G4double> >;
     115  vector<pair<G4int,G4double> >*a8=new vector<pair<G4int,G4double> >; // 8-O
    116116  a8->push_back(make_pair(8,.9976));
    117117  a8->push_back(make_pair(10,.9996));
    118118  a8->push_back(make_pair(9,1.));
    119119  natEl.push_back(a8);
    120   vector<pair<G4int,G4double> >*a9=new vector<pair<G4int,G4double> >;
     120  vector<pair<G4int,G4double> >*a9=new vector<pair<G4int,G4double> >; // 9-F
    121121  a9->push_back(make_pair(10,1.));
    122122  natEl.push_back(a9);
    123   vector<pair<G4int,G4double> >*b0=new vector<pair<G4int,G4double> >;
     123  vector<pair<G4int,G4double> >*b0=new vector<pair<G4int,G4double> >; // 10-Ne
    124124  b0->push_back(make_pair(10,.9948));
    125125  b0->push_back(make_pair(11,.9975));
    126126  b0->push_back(make_pair(12,1.));
    127127  natEl.push_back(b0);
    128   vector<pair<G4int,G4double> >*b1=new vector<pair<G4int,G4double> >;
     128  vector<pair<G4int,G4double> >*b1=new vector<pair<G4int,G4double> >; // 11-Na
    129129  b1->push_back(make_pair(12,1.));
    130130  natEl.push_back(b1);
    131   vector<pair<G4int,G4double> >*b2=new vector<pair<G4int,G4double> >;
     131  vector<pair<G4int,G4double> >*b2=new vector<pair<G4int,G4double> >; // 12-Mg
    132132  b2->push_back(make_pair(12,.7899));
    133133  b2->push_back(make_pair(13,.8899));
    134134  b2->push_back(make_pair(14,1.));
    135135  natEl.push_back(b2);
    136   vector<pair<G4int,G4double> >*b3=new vector<pair<G4int,G4double> >;
     136  vector<pair<G4int,G4double> >*b3=new vector<pair<G4int,G4double> >; // 13-Al
    137137  b3->push_back(make_pair(14,1.));
    138138  natEl.push_back(b3);
    139   vector<pair<G4int,G4double> >*b4=new vector<pair<G4int,G4double> >;
     139  vector<pair<G4int,G4double> >*b4=new vector<pair<G4int,G4double> >; // 14-Si
    140140  b4->push_back(make_pair(14,.9223));
    141141  b4->push_back(make_pair(15,.969));
    142142  b4->push_back(make_pair(16,1.));
    143143  natEl.push_back(b4);
    144   vector<pair<G4int,G4double> >*b5=new vector<pair<G4int,G4double> >;
     144  vector<pair<G4int,G4double> >*b5=new vector<pair<G4int,G4double> >; // 15-P
    145145  b5->push_back(make_pair(16,1.));
    146146  natEl.push_back(b5);
    147   vector<pair<G4int,G4double> >*b6=new vector<pair<G4int,G4double> >;
     147  vector<pair<G4int,G4double> >*b6=new vector<pair<G4int,G4double> >; // 16-S
    148148  b6->push_back(make_pair(16,.9502));
    149149  b6->push_back(make_pair(18,.9923));
    … …  
    151151  b6->push_back(make_pair(20,1.));
    152152  natEl.push_back(b6);
    153   vector<pair<G4int,G4double> >*b7=new vector<pair<G4int,G4double> >;
     153  vector<pair<G4int,G4double> >*b7=new vector<pair<G4int,G4double> >; // 17-Cl
    154154  b7->push_back(make_pair(18,.7577));
    155155  b7->push_back(make_pair(20,1.));
    156156  natEl.push_back(b7);
    157   vector<pair<G4int,G4double> >*b8=new vector<pair<G4int,G4double> >;
     157  vector<pair<G4int,G4double> >*b8=new vector<pair<G4int,G4double> >; // 18-Ar
    158158  b8->push_back(make_pair(22,.996));
    159159  b8->push_back(make_pair(18,.99937));
    160160  b8->push_back(make_pair(20,1.));
    161161  natEl.push_back(b8);
    162   vector<pair<G4int,G4double> >*b9=new vector<pair<G4int,G4double> >;
     162  vector<pair<G4int,G4double> >*b9=new vector<pair<G4int,G4double> >; // 19-K
    163163  b9->push_back(make_pair(20,.932581));
    164164  b9->push_back(make_pair(22,.999883));
    165165  b9->push_back(make_pair(21,1.));
    166166  natEl.push_back(b9);
    167   vector<pair<G4int,G4double> >*c0=new vector<pair<G4int,G4double> >;
     167  vector<pair<G4int,G4double> >*c0=new vector<pair<G4int,G4double> >; // 20-Ca
    168168  c0->push_back(make_pair(20,.96941));
    169169  c0->push_back(make_pair(24,.99027));
    … …  
    173173  c0->push_back(make_pair(26,1.));
    174174  natEl.push_back(c0);
    175   vector<pair<G4int,G4double> >*c1=new vector<pair<G4int,G4double> >;
     175  vector<pair<G4int,G4double> >*c1=new vector<pair<G4int,G4double> >; // 21-Sc
    176176  c1->push_back(make_pair(24,1.));
    177177  natEl.push_back(c1);
    178   vector<pair<G4int,G4double> >*c2=new vector<pair<G4int,G4double> >;
    179   c2->push_back(make_pair(1,.738));
    180   c2->push_back(make_pair(1,.818));
    181   c2->push_back(make_pair(1,.891));
    182   c2->push_back(make_pair(1,.946));
    183   c2->push_back(make_pair(1,1.));
     178  vector<pair<G4int,G4double> >*c2=new vector<pair<G4int,G4double> >; // 22-Ti
     179  c2->push_back(make_pair(26,.738));
     180  c2->push_back(make_pair(24,.818));
     181  c2->push_back(make_pair(25,.891));
     182  c2->push_back(make_pair(27,.946));
     183  c2->push_back(make_pair(28,1.));
    184184  natEl.push_back(c2);
    185   vector<pair<G4int,G4double> >*c3=new vector<pair<G4int,G4double> >;
     185  vector<pair<G4int,G4double> >*c3=new vector<pair<G4int,G4double> >; // 23-V
    186186  c3->push_back(make_pair(28,.9975));
    187187  c3->push_back(make_pair(27,1.));
    188188  natEl.push_back(c3);
    189   vector<pair<G4int,G4double> >*c4=new vector<pair<G4int,G4double> >;
     189  vector<pair<G4int,G4double> >*c4=new vector<pair<G4int,G4double> >; // 24-Cr
    190190  c4->push_back(make_pair(28,.8379));
    191191  c4->push_back(make_pair(29,.9329));
    … …  
    193193  c4->push_back(make_pair(30,1.));
    194194  natEl.push_back(c4);
    195   vector<pair<G4int,G4double> >*c5=new vector<pair<G4int,G4double> >;
     195  vector<pair<G4int,G4double> >*c5=new vector<pair<G4int,G4double> >; // 25-Mn
    196196  c5->push_back(make_pair(30,1.));
    197197  natEl.push_back(c5);
    198   vector<pair<G4int,G4double> >*c6=new vector<pair<G4int,G4double> >;
     198  vector<pair<G4int,G4double> >*c6=new vector<pair<G4int,G4double> >; // 26-Fe
    199199  c6->push_back(make_pair(30,.9172));
    200200  c6->push_back(make_pair(28,.9762));
    … …  
    202202  c6->push_back(make_pair(32,1.));
    203203  natEl.push_back(c6);
    204   vector<pair<G4int,G4double> >*c7=new vector<pair<G4int,G4double> >;
     204  vector<pair<G4int,G4double> >*c7=new vector<pair<G4int,G4double> >; // 27-Co
    205205  c7->push_back(make_pair(32,1.));
    206206  natEl.push_back(c7);
    207   vector<pair<G4int,G4double> >*c8=new vector<pair<G4int,G4double> >;
     207  vector<pair<G4int,G4double> >*c8=new vector<pair<G4int,G4double> >; // 28-Ni
    208208  c8->push_back(make_pair(30,.68077));
    209209  c8->push_back(make_pair(32,.943));
    … …  
    212212  c8->push_back(make_pair(36,1.));
    213213  natEl.push_back(c8);
    214   vector<pair<G4int,G4double> >*c9=new vector<pair<G4int,G4double> >;
     214  vector<pair<G4int,G4double> >*c9=new vector<pair<G4int,G4double> >; // 29-Cu
    215215  c9->push_back(make_pair(34,.6917));
    216216  c9->push_back(make_pair(36,1.));
    217217  natEl.push_back(c9);
    218   vector<pair<G4int,G4double> >*d0=new vector<pair<G4int,G4double> >;
     218  vector<pair<G4int,G4double> >*d0=new vector<pair<G4int,G4double> >; // 30-Zn
    219219  d0->push_back(make_pair(34,.486));
    220220  d0->push_back(make_pair(36,.765));
    … …  
    223223  d0->push_back(make_pair(40,1.));
    224224  natEl.push_back(d0);
    225   vector<pair<G4int,G4double> >*d1=new vector<pair<G4int,G4double> >;
     225  vector<pair<G4int,G4double> >*d1=new vector<pair<G4int,G4double> >; // 31-Ga
    226226  d1->push_back(make_pair(38,.60108));
    227227  d1->push_back(make_pair(40,1.));
    228228  natEl.push_back(d1);
    229   vector<pair<G4int,G4double> >*d2=new vector<pair<G4int,G4double> >;
     229  vector<pair<G4int,G4double> >*d2=new vector<pair<G4int,G4double> >; // 32-Ge
    230230  d2->push_back(make_pair(42,.3594));
    231231  d2->push_back(make_pair(40,.6360));
    … …  
    234234  d2->push_back(make_pair(44,1.));
    235235  natEl.push_back(d2);
    236   vector<pair<G4int,G4double> >*d3=new vector<pair<G4int,G4double> >;
     236  vector<pair<G4int,G4double> >*d3=new vector<pair<G4int,G4double> >; // 33-As
    237237  d3->push_back(make_pair(42,1.));
    238238  natEl.push_back(d3);
    239   vector<pair<G4int,G4double> >*d4=new vector<pair<G4int,G4double> >;
     239  vector<pair<G4int,G4double> >*d4=new vector<pair<G4int,G4double> >; // 34-Se
    240240  d4->push_back(make_pair(46,.4961));
    241241  d4->push_back(make_pair(44,.7378));
    … …  
    245245  d4->push_back(make_pair(40,1.));
    246246  natEl.push_back(d4);
    247   vector<pair<G4int,G4double> >*d5=new vector<pair<G4int,G4double> >;
     247  vector<pair<G4int,G4double> >*d5=new vector<pair<G4int,G4double> >; // 35-Br
    248248  d5->push_back(make_pair(44,.5069));
    249249  d5->push_back(make_pair(46,1.));
    250250  natEl.push_back(d5);
    251   vector<pair<G4int,G4double> >*d6=new vector<pair<G4int,G4double> >;
     251  vector<pair<G4int,G4double> >*d6=new vector<pair<G4int,G4double> >; // 36-Kr
    252252  d6->push_back(make_pair(48,.57));
    253253  d6->push_back(make_pair(50,.743));
    … …  
    257257  d6->push_back(make_pair(42,1.));
    258258  natEl.push_back(d6);
    259   vector<pair<G4int,G4double> >*d7=new vector<pair<G4int,G4double> >;
     259  vector<pair<G4int,G4double> >*d7=new vector<pair<G4int,G4double> >; // 37-Rb
    260260  d7->push_back(make_pair(48,.7217));
    261261  d7->push_back(make_pair(50,1.));
    262262  natEl.push_back(d7);
    263   vector<pair<G4int,G4double> >*d8=new vector<pair<G4int,G4double> >;
     263  vector<pair<G4int,G4double> >*d8=new vector<pair<G4int,G4double> >; // 38-sr
    264264  d8->push_back(make_pair(50,.8258));
    265265  d8->push_back(make_pair(48,.9244));
    … …  
    267267  d8->push_back(make_pair(46,1.));
    268268  natEl.push_back(d8);
    269   vector<pair<G4int,G4double> >*d9=new vector<pair<G4int,G4double> >;
     269  vector<pair<G4int,G4double> >*d9=new vector<pair<G4int,G4double> >; // 39-Y
    270270  d9->push_back(make_pair(50,1.));
    271271  natEl.push_back(d9);
    272   vector<pair<G4int,G4double> >*e0=new vector<pair<G4int,G4double> >;
     272  vector<pair<G4int,G4double> >*e0=new vector<pair<G4int,G4double> >; // 40-Zr
    273273  e0->push_back(make_pair(50,.5145));
    274274  e0->push_back(make_pair(54,.6883));
    275   e0->push_back(make_pair(53,.8598));
     275  e0->push_back(make_pair(52,.8598));
    276276  e0->push_back(make_pair(51,.972));
    277277  e0->push_back(make_pair(56,1.));
    278278  natEl.push_back(e0);
    279   vector<pair<G4int,G4double> >*e1=new vector<pair<G4int,G4double> >;
     279  vector<pair<G4int,G4double> >*e1=new vector<pair<G4int,G4double> >; // 41-Nb
    280280  e1->push_back(make_pair(52,1.));
    281281  natEl.push_back(e1);
    282   vector<pair<G4int,G4double> >*e2=new vector<pair<G4int,G4double> >;
     282  vector<pair<G4int,G4double> >*e2=new vector<pair<G4int,G4double> >; // 42-Mo
    283283  e2->push_back(make_pair(56,.2413));
    284284  e2->push_back(make_pair(54,.4081));
    … …  
    289289  e2->push_back(make_pair(52,1.));
    290290  natEl.push_back(e2);
    291   vector<pair<G4int,G4double> >*e3=new vector<pair<G4int,G4double> >;
     291  vector<pair<G4int,G4double> >*e3=new vector<pair<G4int,G4double> >; // 43-Tc
    292292  e3->push_back(make_pair(55,1.));
    293293  natEl.push_back(e3);
    294   vector<pair<G4int,G4double> >*e4=new vector<pair<G4int,G4double> >;
     294  vector<pair<G4int,G4double> >*e4=new vector<pair<G4int,G4double> >; // 44-Ru
    295295  e4->push_back(make_pair(58,.316));
    296296  e4->push_back(make_pair(60,.502));
    … …  
    301301  e4->push_back(make_pair(54,1.));
    302302  natEl.push_back(e4);
    303   vector<pair<G4int,G4double> >*e5=new vector<pair<G4int,G4double> >;
     303  vector<pair<G4int,G4double> >*e5=new vector<pair<G4int,G4double> >; // 45-Rh
    304304  e5->push_back(make_pair(58,1.));
    305305  natEl.push_back(e5);
    306   vector<pair<G4int,G4double> >*e6=new vector<pair<G4int,G4double> >;
     306  vector<pair<G4int,G4double> >*e6=new vector<pair<G4int,G4double> >; // 46-Pd
    307307  e6->push_back(make_pair(60,.2733));
    308308  e6->push_back(make_pair(62,.5379));
    … …  
    312312  e6->push_back(make_pair(56,1.));
    313313  natEl.push_back(e6);
    314   vector<pair<G4int,G4double> >*e7=new vector<pair<G4int,G4double> >;
     314  vector<pair<G4int,G4double> >*e7=new vector<pair<G4int,G4double> >; // 47-Ag
    315315  e7->push_back(make_pair(60,.51839));
    316316  e7->push_back(make_pair(62,1.));
    317317  natEl.push_back(e7);
    318   vector<pair<G4int,G4double> >*e8=new vector<pair<G4int,G4double> >;
     318  vector<pair<G4int,G4double> >*e8=new vector<pair<G4int,G4double> >; // 48-Cd
    319319  e8->push_back(make_pair(66,.2873));
    320320  e8->push_back(make_pair(64,.5286));
    … …  
    326326  e8->push_back(make_pair(60,1.));
    327327  natEl.push_back(e8);
    328   vector<pair<G4int,G4double> >*e9=new vector<pair<G4int,G4double> >;
     328  vector<pair<G4int,G4double> >*e9=new vector<pair<G4int,G4double> >; // 49-In
    329329  e9->push_back(make_pair(66,.9577));
    330330  e9->push_back(make_pair(64,1.));
    331331  natEl.push_back(e9);
    332   vector<pair<G4int,G4double> >*f0=new vector<pair<G4int,G4double> >;
     332  vector<pair<G4int,G4double> >*f0=new vector<pair<G4int,G4double> >; // 50-Sn
    333333  f0->push_back(make_pair(70,.3259));
    334334  f0->push_back(make_pair(68,.5681));
    … …  
    341341  f0->push_back(make_pair(64,1.));
    342342  //f0->push_back(make_pair(64,.9964));
    343   //f0->push_back(make_pair(65,1.)); // Nine isotopes is the maximum, so Sn155 is out
     343  //f0->push_back(make_pair(65,1.)); // Nine isotopes is the maximum, so Sn115 is out
    344344  natEl.push_back(f0);
    345   vector<pair<G4int,G4double> >*f1=new vector<pair<G4int,G4double> >;
     345  vector<pair<G4int,G4double> >*f1=new vector<pair<G4int,G4double> >; // 51-Sb
    346346  f1->push_back(make_pair(70,.5736));
    347347  f1->push_back(make_pair(72,1.));
    348348  natEl.push_back(f1);
    349   vector<pair<G4int,G4double> >*f2=new vector<pair<G4int,G4double> >;
     349  vector<pair<G4int,G4double> >*f2=new vector<pair<G4int,G4double> >; // 52-Te
    350350  f2->push_back(make_pair(78,.3387));
    351351  f2->push_back(make_pair(76,.6557));
    … …  
    357357  f2->push_back(make_pair(68,1.));
    358358  natEl.push_back(f2);
    359   vector<pair<G4int,G4double> >*f3=new vector<pair<G4int,G4double> >;
     359  vector<pair<G4int,G4double> >*f3=new vector<pair<G4int,G4double> >; // 53-I
    360360  f3->push_back(make_pair(74,1.));
    361361  natEl.push_back(f3);
    362   vector<pair<G4int,G4double> >*f4=new vector<pair<G4int,G4double> >;
     362  vector<pair<G4int,G4double> >*f4=new vector<pair<G4int,G4double> >; // 54-Xe
    363363  f4->push_back(make_pair(78,.269));
    364364  f4->push_back(make_pair(75,.533));
    … …  
    371371  f4->push_back(make_pair(72,1.));
    372372  natEl.push_back(f4);
    373   vector<pair<G4int,G4double> >*f5=new vector<pair<G4int,G4double> >;
     373  vector<pair<G4int,G4double> >*f5=new vector<pair<G4int,G4double> >; // 55-Cs
    374374  f5->push_back(make_pair(78,1.));
    375375  natEl.push_back(f5);
    376   vector<pair<G4int,G4double> >*f6=new vector<pair<G4int,G4double> >;
     376  vector<pair<G4int,G4double> >*f6=new vector<pair<G4int,G4double> >; // 56-Ba
    377377  f6->push_back(make_pair(82,.717));
    378378  f6->push_back(make_pair(81,.8293));
    … …  
    383383  f6->push_back(make_pair(76,1.));
    384384  natEl.push_back(f6);
    385   vector<pair<G4int,G4double> >*f7=new vector<pair<G4int,G4double> >;
     385  vector<pair<G4int,G4double> >*f7=new vector<pair<G4int,G4double> >; // 57-La
    386386  f7->push_back(make_pair(82,.999098));
    387387  f7->push_back(make_pair(81,1.));
    388388  natEl.push_back(f7);
    389   vector<pair<G4int,G4double> >*f8=new vector<pair<G4int,G4double> >;
     389  vector<pair<G4int,G4double> >*f8=new vector<pair<G4int,G4double> >; // 58-Ce
    390390  f8->push_back(make_pair(82,.8843));
    391391  f8->push_back(make_pair(84,.9956));
    … …  
    393393  f8->push_back(make_pair(78,1.));
    394394  natEl.push_back(f8);
    395   vector<pair<G4int,G4double> >*f9=new vector<pair<G4int,G4double> >;
     395  vector<pair<G4int,G4double> >*f9=new vector<pair<G4int,G4double> >; // 59-Pr
    396396  f9->push_back(make_pair(82,1.));
    397397  natEl.push_back(f9);
    398   vector<pair<G4int,G4double> >*g0=new vector<pair<G4int,G4double> >;
     398  vector<pair<G4int,G4double> >*g0=new vector<pair<G4int,G4double> >; // 60-Nd
    399399  g0->push_back(make_pair(82,.2713));
    400400  g0->push_back(make_pair(84,.5093));
    … …  
    405405  g0->push_back(make_pair(90,1.));
    406406  natEl.push_back(g0);
    407   vector<pair<G4int,G4double> >*g1=new vector<pair<G4int,G4double> >;
     407  vector<pair<G4int,G4double> >*g1=new vector<pair<G4int,G4double> >; // 61-Pm
    408408  g1->push_back(make_pair(85,1.));
    409409  natEl.push_back(g1);
    410   vector<pair<G4int,G4double> >*g2=new vector<pair<G4int,G4double> >;
     410  vector<pair<G4int,G4double> >*g2=new vector<pair<G4int,G4double> >; // 62-Sm
    411411  g2->push_back(make_pair(90,.267));
    412412  g2->push_back(make_pair(92,.494));
    … …  
    417417  g2->push_back(make_pair(82,1.));
    418418  natEl.push_back(g2);
    419   vector<pair<G4int,G4double> >*g3=new vector<pair<G4int,G4double> >;
     419  vector<pair<G4int,G4double> >*g3=new vector<pair<G4int,G4double> >; // 63-Eu
    420420  g3->push_back(make_pair(90,.522));
    421421  g3->push_back(make_pair(89,1.));
    422422  natEl.push_back(g3);
    423   vector<pair<G4int,G4double> >*g4=new vector<pair<G4int,G4double> >;
     423  vector<pair<G4int,G4double> >*g4=new vector<pair<G4int,G4double> >; // 64-Gd
    424424  g4->push_back(make_pair(94,.2484));
    425425  g4->push_back(make_pair(96,.4670));
    … …  
    430430  g4->push_back(make_pair(88,1.));
    431431  natEl.push_back(g4);
    432   vector<pair<G4int,G4double> >*g5=new vector<pair<G4int,G4double> >;
     432  vector<pair<G4int,G4double> >*g5=new vector<pair<G4int,G4double> >; // 65-Tb
    433433  g5->push_back(make_pair(94,1.));
    434434  natEl.push_back(g5);
    435   vector<pair<G4int,G4double> >*g6=new vector<pair<G4int,G4double> >;
     435  vector<pair<G4int,G4double> >*g6=new vector<pair<G4int,G4double> >; // 66-Dy
    436436  g6->push_back(make_pair(98,.282));
    437437  g6->push_back(make_pair(96,.537));
    … …  
    442442  g6->push_back(make_pair(90,1.));
    443443  natEl.push_back(g6);
    444   vector<pair<G4int,G4double> >*g7=new vector<pair<G4int,G4double> >;
     444  vector<pair<G4int,G4double> >*g7=new vector<pair<G4int,G4double> >; // 67-Ho
    445445  g7->push_back(make_pair(98,1.));
    446446  natEl.push_back(g7);
    447   vector<pair<G4int,G4double> >*g8=new vector<pair<G4int,G4double> >;
     447  vector<pair<G4int,G4double> >*g8=new vector<pair<G4int,G4double> >; // 68-Er
    448448  g8->push_back(make_pair( 98,.3360));
    449449  g8->push_back(make_pair(100,.6040));
    … …  
    453453  g8->push_back(make_pair( 94,1.));
    454454  natEl.push_back(g8);
    455   vector<pair<G4int,G4double> >*g9=new vector<pair<G4int,G4double> >;
     455  vector<pair<G4int,G4double> >*g9=new vector<pair<G4int,G4double> >; // 69-Tm
    456456  g9->push_back(make_pair(100,1.));
    457457  natEl.push_back(g9);
    458   vector<pair<G4int,G4double> >*h0=new vector<pair<G4int,G4double> >;
     458  vector<pair<G4int,G4double> >*h0=new vector<pair<G4int,G4double> >; // 70-Yb
    459459  h0->push_back(make_pair(104,.3180));
    460460  h0->push_back(make_pair(102,.5370));
    … …  
    465465  h0->push_back(make_pair( 98,1.));
    466466  natEl.push_back(h0);
    467   vector<pair<G4int,G4double> >*h1=new vector<pair<G4int,G4double> >;
     467  vector<pair<G4int,G4double> >*h1=new vector<pair<G4int,G4double> >; // 71-Lu
    468468  h1->push_back(make_pair(104,.9741));
    469469  h1->push_back(make_pair(105,1.));
    470470  natEl.push_back(h1);
    471   vector<pair<G4int,G4double> >*h2=new vector<pair<G4int,G4double> >;
     471  vector<pair<G4int,G4double> >*h2=new vector<pair<G4int,G4double> >; // 72-Hf
    472472  h2->push_back(make_pair(108,.35100));
    473473  h2->push_back(make_pair(106,.62397));
    … …  
    477477  h2->push_back(make_pair(102,1.));
    478478  natEl.push_back(h2);
    479   vector<pair<G4int,G4double> >*h3=new vector<pair<G4int,G4double> >;
     479  vector<pair<G4int,G4double> >*h3=new vector<pair<G4int,G4double> >; // 73-Ta
    480480  h3->push_back(make_pair(108,.99988));
    481481  h3->push_back(make_pair(107,1.));
    482482  natEl.push_back(h3);
    483   vector<pair<G4int,G4double> >*h4=new vector<pair<G4int,G4double> >;
     483  vector<pair<G4int,G4double> >*h4=new vector<pair<G4int,G4double> >; // 74-W
    484484  h4->push_back(make_pair(110,.307));
    485485  h4->push_back(make_pair(112,.593));
    … …  
    488488  h4->push_back(make_pair(106,1.));
    489489  natEl.push_back(h4);
    490   vector<pair<G4int,G4double> >*h5=new vector<pair<G4int,G4double> >;
     490  vector<pair<G4int,G4double> >*h5=new vector<pair<G4int,G4double> >; // 75-Re
    491491  h5->push_back(make_pair(112,.626));
    492492  h5->push_back(make_pair(110,1.));
    493493  natEl.push_back(h5);
    494   vector<pair<G4int,G4double> >*h6=new vector<pair<G4int,G4double> >;
     494  vector<pair<G4int,G4double> >*h6=new vector<pair<G4int,G4double> >; // 78-Os
    495495  h6->push_back(make_pair(116,.410));
    496496  h6->push_back(make_pair(114,.674));
    … …  
    501501  h6->push_back(make_pair(108,1.));
    502502  natEl.push_back(h6);
    503   vector<pair<G4int,G4double> >*h7=new vector<pair<G4int,G4double> >;
     503  vector<pair<G4int,G4double> >*h7=new vector<pair<G4int,G4double> >; // 77-Ir
    504504  h7->push_back(make_pair(116,.627));
    505505  h7->push_back(make_pair(114,1.));
    506506  natEl.push_back(h7);
    507   vector<pair<G4int,G4double> >*h8=new vector<pair<G4int,G4double> >;
     507  vector<pair<G4int,G4double> >*h8=new vector<pair<G4int,G4double> >; // 78-Pt
    508508  h8->push_back(make_pair(117,.338));
    509509  h8->push_back(make_pair(116,.667));
    … …  
    513513  h8->push_back(make_pair(112,1.));
    514514  natEl.push_back(h8);
    515   vector<pair<G4int,G4double> >*h9=new vector<pair<G4int,G4double> >;
     515  vector<pair<G4int,G4double> >*h9=new vector<pair<G4int,G4double> >; // 79-Au
    516516  h9->push_back(make_pair(118,1.));
    517517  natEl.push_back(h9);
    518   vector<pair<G4int,G4double> >*i0=new vector<pair<G4int,G4double> >;
     518  vector<pair<G4int,G4double> >*i0=new vector<pair<G4int,G4double> >; // 80-Hg
    519519  i0->push_back(make_pair(122,.2986));
    520520  i0->push_back(make_pair(120,.5296));
    … …  
    525525  i0->push_back(make_pair(116,1.));
    526526  natEl.push_back(i0);
    527   vector<pair<G4int,G4double> >*i1=new vector<pair<G4int,G4double> >;
     527  vector<pair<G4int,G4double> >*i1=new vector<pair<G4int,G4double> >; // 81-Tl
    528528  i1->push_back(make_pair(124,.70476));
    529529  i1->push_back(make_pair(122,1.));
    530530  natEl.push_back(i1);
    531   vector<pair<G4int,G4double> >*i2=new vector<pair<G4int,G4double> >;
     531  vector<pair<G4int,G4double> >*i2=new vector<pair<G4int,G4double> >; // 82-Pb
    532532  i2->push_back(make_pair(126,.524));
    533533  i2->push_back(make_pair(124,.765));
    … …  
    535535  i2->push_back(make_pair(122,1.));
    536536  natEl.push_back(i2);
    537   vector<pair<G4int,G4double> >*i3=new vector<pair<G4int,G4double> >;
     537  vector<pair<G4int,G4double> >*i3=new vector<pair<G4int,G4double> >; // 83-Bi
    538538  i3->push_back(make_pair(126,1.));
    539539  natEl.push_back(i3);
    540   vector<pair<G4int,G4double> >*i4=new vector<pair<G4int,G4double> >;
     540  vector<pair<G4int,G4double> >*i4=new vector<pair<G4int,G4double> >; // 84-Po
    541541  i4->push_back(make_pair(125,1.));
    542542  natEl.push_back(i4);
    543   vector<pair<G4int,G4double> >*i5=new vector<pair<G4int,G4double> >;
     543  vector<pair<G4int,G4double> >*i5=new vector<pair<G4int,G4double> >; // 85-At
    544544  i5->push_back(make_pair(136,1.));
    545545  natEl.push_back(i5);
    546   vector<pair<G4int,G4double> >*i6=new vector<pair<G4int,G4double> >;
     546  vector<pair<G4int,G4double> >*i6=new vector<pair<G4int,G4double> >; // 86-Ru
    547547  i6->push_back(make_pair(136,1.));
    548548  natEl.push_back(i6);
    549   vector<pair<G4int,G4double> >*i7=new vector<pair<G4int,G4double> >;
     549  vector<pair<G4int,G4double> >*i7=new vector<pair<G4int,G4double> >; // 87-Fr
    550550  i7->push_back(make_pair(138,1.));
    551551  natEl.push_back(i7);
    552   vector<pair<G4int,G4double> >*i8=new vector<pair<G4int,G4double> >;
     552  vector<pair<G4int,G4double> >*i8=new vector<pair<G4int,G4double> >; // 88-Ra
    553553  i8->push_back(make_pair(138,1.));
    554554  natEl.push_back(i8);
    555   vector<pair<G4int,G4double> >*i9=new vector<pair<G4int,G4double> >;
     555  vector<pair<G4int,G4double> >*i9=new vector<pair<G4int,G4double> >; // 89-Ac
    556556  i9->push_back(make_pair(142,1.));
    557557  natEl.push_back(i9);
    558   vector<pair<G4int,G4double> >*j0=new vector<pair<G4int,G4double> >;
     558  vector<pair<G4int,G4double> >*j0=new vector<pair<G4int,G4double> >; // 90-Th
    559559  j0->push_back(make_pair(142,1.));
    560560  natEl.push_back(j0);
    561   vector<pair<G4int,G4double> >*j1=new vector<pair<G4int,G4double> >;
     561  vector<pair<G4int,G4double> >*j1=new vector<pair<G4int,G4double> >; // 91-Pa
    562562  j1->push_back(make_pair(140,1.));
    563563  natEl.push_back(j1);
    564   vector<pair<G4int,G4double> >*j2=new vector<pair<G4int,G4double> >;
     564  vector<pair<G4int,G4double> >*j2=new vector<pair<G4int,G4double> >; // 92-U
    565565  j2->push_back(make_pair(146,.992745));
    566566  j2->push_back(make_pair(143,.999945));
    567567  j2->push_back(make_pair(142,1.));
    568568  natEl.push_back(j2);
    569   vector<pair<G4int,G4double> >*j3=new vector<pair<G4int,G4double> >;
     569  vector<pair<G4int,G4double> >*j3=new vector<pair<G4int,G4double> >; // 93-Np
    570570  j3->push_back(make_pair(144,1.));
    571571  natEl.push_back(j3);
    572   vector<pair<G4int,G4double> >*j4=new vector<pair<G4int,G4double> >;
     572  vector<pair<G4int,G4double> >*j4=new vector<pair<G4int,G4double> >; // 94-Pu
    573573  j4->push_back(make_pair(150,1.));
    574574  natEl.push_back(j4);
    575   vector<pair<G4int,G4double> >*j5=new vector<pair<G4int,G4double> >;
     575  vector<pair<G4int,G4double> >*j5=new vector<pair<G4int,G4double> >; // 95-Am
    576576  j5->push_back(make_pair(148,1.));
    577577  natEl.push_back(j5);
    578   vector<pair<G4int,G4double> >*j6=new vector<pair<G4int,G4double> >;
     578  vector<pair<G4int,G4double> >*j6=new vector<pair<G4int,G4double> >; // 96-Cm
    579579  j6->push_back(make_pair(151,1.));
    580580  natEl.push_back(j6);
    581   vector<pair<G4int,G4double> >*j7=new vector<pair<G4int,G4double> >;
     581  vector<pair<G4int,G4double> >*j7=new vector<pair<G4int,G4double> >; // 97-Bk
    582582  j7->push_back(make_pair(150,1.));
    583583  natEl.push_back(j7);
    584   vector<pair<G4int,G4double> >*j8=new vector<pair<G4int,G4double> >;
     584  vector<pair<G4int,G4double> >*j8=new vector<pair<G4int,G4double> >; // 98-Cf
    585585  j8->push_back(make_pair(153,1.));
    586586  natEl.push_back(j8);
    587   vector<pair<G4int,G4double> >*j9=new vector<pair<G4int,G4double> >;
     587  vector<pair<G4int,G4double> >*j9=new vector<pair<G4int,G4double> >; // 99-Es
    588588  j9->push_back(make_pair(157,1.));
    589589  natEl.push_back(j9);
    590   vector<pair<G4int,G4double> >*k0=new vector<pair<G4int,G4double> >;
     590  vector<pair<G4int,G4double> >*k0=new vector<pair<G4int,G4double> >; // 100-Fm
    591591  k0->push_back(make_pair(157,1.));
    592592  natEl.push_back(k0);
    593   vector<pair<G4int,G4double> >*k1=new vector<pair<G4int,G4double> >;
     593  vector<pair<G4int,G4double> >*k1=new vector<pair<G4int,G4double> >; // 101-Md
    594594  k1->push_back(make_pair(157,1.));
    595595  natEl.push_back(k1);
    596   vector<pair<G4int,G4double> >*k2=new vector<pair<G4int,G4double> >;
     596  vector<pair<G4int,G4double> >*k2=new vector<pair<G4int,G4double> >; // 102-No
    597597  k2->push_back(make_pair(157,1.));
    598598  natEl.push_back(k2);
    599   vector<pair<G4int,G4double> >*k3=new vector<pair<G4int,G4double> >;
     599  vector<pair<G4int,G4double> >*k3=new vector<pair<G4int,G4double> >; // 103-Lr
    600600  k3->push_back(make_pair(157,1.));
    601601  natEl.push_back(k3);
    602   vector<pair<G4int,G4double> >*k4=new vector<pair<G4int,G4double> >;
     602  vector<pair<G4int,G4double> >*k4=new vector<pair<G4int,G4double> >; // 104-Rf
    603603  k4->push_back(make_pair(157,1.));
    604604  natEl.push_back(k4);
    605   vector<pair<G4int,G4double> >*k5=new vector<pair<G4int,G4double> >;
     605  vector<pair<G4int,G4double> >*k5=new vector<pair<G4int,G4double> >; // 105-Db
    606606  k5->push_back(make_pair(157,1.));
    607607  natEl.push_back(k5);
    608   vector<pair<G4int,G4double> >*k6=new vector<pair<G4int,G4double> >;
     608  vector<pair<G4int,G4double> >*k6=new vector<pair<G4int,G4double> >; // 106-Sg
    609609  k6->push_back(make_pair(157,1.));
    610610  natEl.push_back(k6);
    611   vector<pair<G4int,G4double> >*k7=new vector<pair<G4int,G4double> >;
     611  vector<pair<G4int,G4double> >*k7=new vector<pair<G4int,G4double> >; // 107-Bh
    612612  k7->push_back(make_pair(155,1.));
    613613  natEl.push_back(k7);
    614   vector<pair<G4int,G4double> >*k8=new vector<pair<G4int,G4double> >;
     614  vector<pair<G4int,G4double> >*k8=new vector<pair<G4int,G4double> >; // 108-Hs
    615615  k8->push_back(make_pair(157,1.));
    616616  natEl.push_back(k8);
    617   vector<pair<G4int,G4double> >*k9=new vector<pair<G4int,G4double> >;
     617  vector<pair<G4int,G4double> >*k9=new vector<pair<G4int,G4double> >; // 109-Mt
    618618  k9->push_back(make_pair(157,1.));
    619619  natEl.push_back(k9);

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QNucleus.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4QNucleus.cc,v 1.116 2009/12/14 16:41:52 mkossov Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4QNucleus.cc,v 1.118 2010/06/10 08:37:27 mkossov Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//      ---------------- G4QNucleus ----------------
    … …  
    958958  static const G4int      nPDG = 2112;            // PDGCode of neutron
    959959  static const G4QPDGCode nQPDG(nPDG);            // QPDGCode of neutron
     960  static const G4QPDGCode anQPDG(-nPDG);          // QPDGCode of anti-neutron
    960961  static const G4int      pPDG = 2212;            // PDGCode of proton
    961962  static const G4QPDGCode pQPDG(pPDG);            // QPDGCode of proton
     963  static const G4QPDGCode apQPDG(-pPDG);          // QPDGCode of anti-proton
    962964  static const G4int      lPDG = 3122;            // PDGCode of Lambda
    963965  static const G4QPDGCode lQPDG(lPDG);            // QPDGCode of Lambda
     966  static const G4QPDGCode aDppQPDG(-2224);        // QPDGCode of anti-Delta++
     967  static const G4QPDGCode aDmQPDG(-1114);         // QPDGCode of anti-Delta-
     968  static const G4QPDGCode alQPDG(-lPDG);          // QPDGCode of anti-Lambda
    964969  static const G4int      dPDG = 90001001;        // PDGCode of deutron
    965970  static const G4int      aPDG = 90002002;        // PDGCode of ALPHA
    … …  
    10311036  G4cout<<"G4QN::EB:pB="<<PBarr<<",aB="<<ABarr<<",ppB="<<PPBarr<<",paB="<<PABarr<<G4endl;
    10321037#endif
    1033   if(a==2)
     1038  if(a==-2)
     1039  {
     1040    if(Z==1 || N==1)
     1041    {
     1042      G4int  nucPDG  = -2112;
     1043      G4int  piPDG   =  211;
     1044      G4double nucM  = mNeut;
     1045      G4QPDGCode del = aDmQPDG;
     1046      G4QPDGCode nuc = anQPDG;
     1047      if(N>0)
     1048      {
     1049        nucPDG = -2212;
     1050        piPDG  = -211;
     1051        nucM   = mProt;
     1052        del    = aDppQPDG;
     1053        nuc    = apQPDG;
     1054      }
     1055      if(totMass > mPi+nucM+nucM)
     1056      {
     1057        G4LorentzVector n14M(0.,0.,0.,nucM);
     1058        G4LorentzVector n24M(0.,0.,0.,nucM);
     1059        G4LorentzVector pi4M(0.,0.,0.,mPi);
     1060        if(!DecayIn3(n14M, n24M, pi4M))
     1061        {
     1062          G4cerr<<"***G4QNucl::EvapBary: (anti) tM="<<totMass<<"-> 2N="<<nucPDG<<"(M="
     1063                <<nucM<<") + pi="<<piPDG<<"(M="<<mPi<<")"<<G4endl;
     1064          //throw G4QException("G4QNucl::EvapBary:ISO-dibaryon DecayIn3 did not succeed");
     1065          return false;
     1066        }
     1067        n14M+=pi4M;
     1068        h1->SetQPDG(del);
     1069        h2->SetQPDG(nuc);
     1070        h1->Set4Momentum(n14M);
     1071        h2->Set4Momentum(n24M);
     1072        return true;
     1073      }
     1074      else
     1075      {
     1076        G4cerr<<"***G4QNucleus::EvaporateBaryon: M="<<totMass
     1077              <<", M="<<totMass<<" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<G4endl;
     1078        //throw G4QException("***G4QNucl::EvaporateBaryon: ISO-dibaryon under Mass Shell");
     1079        return false;
     1080      }     
     1081    }
     1082    else if(Z==2 || N==2)
     1083    {
     1084      G4int  nucPDG  = -2112;
     1085      G4int  piPDG   =  211;
     1086      G4double nucM  = mNeut;
     1087      G4QPDGCode del = aDmQPDG;
     1088      if(N==2)
     1089      {
     1090        nucPDG = -2212;
     1091        piPDG  = -211;
     1092        nucM   = mProt;
     1093        del    = aDppQPDG;
     1094      }
     1095      if(totMass > mPi+mPi+nucM+nucM)
     1096      {
     1097        G4LorentzVector n14M(0.,0.,0.,nucM);
     1098        G4LorentzVector n24M(0.,0.,0.,nucM);
     1099        G4LorentzVector pi4M(0.,0.,0.,mPi+mPi);
     1100        if(!DecayIn3(n14M, n24M, pi4M))
     1101        {
     1102          G4cerr<<"***G4QNucl::EvapBary: (anti) tM="<<totMass<<"-> 2N="<<nucPDG<<"(M="
     1103                <<nucM<<") + 2pi="<<piPDG<<"(M="<<mPi<<")"<<G4endl;
     1104          //throw G4QException("G4QNucl::EvapBary:ISO-dibaryon DecayIn3 did not succeed");
     1105          return false;
     1106        }
     1107        G4LorentzVector hpi4M=pi4M/2.;
     1108        n14M+=hpi4M;
     1109        n24M+=hpi4M;
     1110        h1->SetQPDG(del);
     1111        h2->SetQPDG(del);
     1112        h1->Set4Momentum(n14M);
     1113        h2->Set4Momentum(n24M);
     1114        return true;
     1115      }
     1116      else
     1117      {
     1118        G4cerr<<"***G4QNucleus::EvaporateBaryon: M="<<totMass
     1119              <<", M="<<totMass<<" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<G4endl;
     1120        //throw G4QException("***G4QNucl::EvaporateBaryon: ISO-dibaryon under Mass Shell");
     1121        return false;
     1122      }     
     1123    }
     1124    else if(Z==-2)
     1125    {
     1126      h1mom=G4LorentzVector(0.,0.,0.,mProt);
     1127      h2mom=h1mom;
     1128      h1->SetQPDG(apQPDG);
     1129      h2->SetQPDG(apQPDG);
     1130      if(!DecayIn2(h1mom,h2mom)) return false;
     1131    }
     1132    else if(N==-2)
     1133    {
     1134      h1mom=G4LorentzVector(0.,0.,0.,mNeut);
     1135      h2mom=h1mom;
     1136      h1->SetQPDG(anQPDG);
     1137      h2->SetQPDG(anQPDG);
     1138      if(!DecayIn2(h1mom,h2mom)) return false;
     1139    }
     1140    else if(N==-1 && Z==-1)
     1141    {
     1142      h1mom=G4LorentzVector(0.,0.,0.,mProt);
     1143      h2mom=G4LorentzVector(0.,0.,0.,mNeut);
     1144      h1->SetQPDG(apQPDG);
     1145      h2->SetQPDG(anQPDG);
     1146      if(!DecayIn2(h1mom,h2mom)) return false;
     1147    }
     1148    else if(Z==-1 && S==-1)
     1149    {
     1150      h1mom=G4LorentzVector(0.,0.,0.,mProt);
     1151      h2mom=G4LorentzVector(0.,0.,0.,mLamb);
     1152      h1->SetQPDG(apQPDG);
     1153      h2->SetQPDG(alQPDG);
     1154      if(!DecayIn2(h1mom,h2mom)) return false;
     1155    }
     1156    else
     1157    {
     1158      h1mom=G4LorentzVector(0.,0.,0.,mNeut);
     1159      h2mom=G4LorentzVector(0.,0.,0.,mLamb);
     1160      h1->SetQPDG(anQPDG);
     1161      h2->SetQPDG(alQPDG);
     1162      if(!DecayIn2(h1mom,h2mom)) return false;
     1163    }
     1164    h1->Set4Momentum(h1mom);
     1165    h2->Set4Momentum(h2mom);
     1166    return true;
     1167  }
     1168  else if(a==2)
    10341169  {
    10351170    if(Z<0||N<0)
    … …  
    15651700      if(evalph&&aFlag&&aMin<minE) minE=aMin;
    15661701
    1567 #ifdef ppdebug
     1702#ifdef pdebug
    15681703      G4cout<<"G4QNucleus::EvapBaryon: nE="<<nExcess<<">"<<nMin<<",pE="<<pExcess<<">"<<pMin
    15691704            <<",sE="<<lExcess<<">"<<lMin<<",E="<<aExcess<<">"<<aMin<<",miE="<<minE<<"<maE="
    … …  
    20172152          else if(PDG==lPDG&&(lnCond&&lpCond&&llCond&&laCond)) // l+b+RN decay can't happen
    20182153          { //@@ Take into account Coulomb Barier Penetration Probability
    2019 #ifdef ppdebug
     2154#ifdef pdebug
    20202155            G4cout<<"G4QN::EB:*l*: n="<<lnCond<<",p="<<lpCond<<",l="<<llCond<<",a="<<laCond
    20212156                  <<G4endl;
    … …  
    22742409        if(!DecayIn3(h1mom,h2mom,h3mom))
    22752410        {
    2276 #ifdef ppdebug
     2411#ifdef pdebug
    22772412          G4cout<<"*G4QNucl::EvaporateBaryon:Decay M="<<totMass<<",b="<<eMass<<bQPDG
    22782413          <<",f="<<fMass<<fQPDG<<",r="<<rMass<<rQPDG<<G4endl;
    … …  
    26892824          rQPDG=LQPDG;
    26902825          rMass=GSResNl;
    2691 #ifdef ppdebug
     2826#ifdef pdebug
    26922827          G4cout<<"G4QNucleus::EvaporateBaryon: Decay in L + rA="<<GSResNl+mLamb<<G4endl;
    26932828#endif
    … …  
    27992934    return true;
    28002935  }
    2801   if(a<1) G4cerr<<"***G4QNucleus::EvaporateBaryon: A="<<a<<G4endl;
     2936  if(a<-2) G4cerr<<"***G4QNucleus::EvaporateBaryon: A="<<a<<G4endl;
    28022937  return false;
    28032938}
    … …  
    33873522  G4ThreeVector  delta(0.,0.,0.);             // Prototype of the distance between nucleons
    33883523  G4ThreeVector* places= new G4ThreeVector[theA]; // Vector of 3D positions
    3389   G4bool      freeplace= false;               // flag of free space available
     3524  G4bool         freeplace= false;            // flag of free space available
    33903525  G4double nucDist2= nucleonDistance*nucleonDistance; // @@ can be a common static
    33913526  G4double maxR= GetRadius(0.01);             // there are cond no nucleons at this density
    3392   G4ThreeVector sumPos(0.,0.,0.);             // Vector of the current 3D sum
    3393   while(i<theA)                               // LOOP over all nucleons
     3527  G4ThreeVector  sumPos(0.,0.,0.);            // Vector of the current 3D sum
     3528  G4ThreeVector  minPos(0.,0.,0.);            // Minimum nucleon 3D position
     3529  G4double       mirPos=maxR;                 // Proto MinimumRadius of the nucleonPosition
     3530  G4int failCNT=0;                            // Counter of bad attempts
     3531  G4int maxCNT=27;                            // Limit for the bad attempts
     3532  while( i < theA && maxR > 0.)               // LOOP over all nucleons
    33943533  {
    33953534    rPos = maxR*pow(G4UniformRand(),third);   // Get random radius of the nucleon position
    … …  
    34043543      // @@ Gaussian oscilator distribution is good only up to He4 (s-wave). Above: p-wave
    34053544      // (1+k*(r^2/R^2)]*exp[-r^2/R^2]. A=s+p=4+3*4=16 (M.K.) So Li,Be,C,N,O are wrong
    3406       if(i==lastN) aPos=-rPos*sumPos.unit();  // TheLast tries to compensate CenterOfGravity
     3545      if(i==lastN) aPos=-rPos*sumPos.unit();  // TheLast tries toCompensate CenterOfGravity
    34073546      else     aPos=rPos*G4RandomDirection(); // It uses the standard G4 function
    34083547      freeplace = true;                       // Imply that there is a free space
    … …  
    34263565        if (eFermi <= CoulombBarrier()) freeplace=false;
    34273566      }
    3428       if(freeplace)
    3429       {
    3430 #ifdef debug
    3431         G4cout<<"G4QNucl::ChoosePositions: fill nucleon i="<<i<<", V="<<aPos<<G4endl;
     3567      if(rPos<mirPos)
     3568      {
     3569        mirPos=rPos;
     3570        minPos=aPos;
     3571      }
     3572      if( freeplace || failCNT > maxCNT )
     3573      {
     3574        if( failCNT > maxCNT ) aPos=minPos;
     3575#ifdef debug
     3576        G4cout<<"G4QNuc::ChoosePos:->> fill N["<<i<<"], R="<<aPos<<", f="<<failCNT<<G4endl;
    34323577#endif     
    34333578        places[i]=aPos;
    34343579        sumPos+=aPos;
    34353580        ++i;
    3436       }
     3581        failCNT=0;
     3582        mirPos=maxR;
     3583      }
     3584      else ++failCNT;
    34373585    }
    34383586  }
    … …  
    39614109  if((thePDG || thePDG==10) && theQC.GetBaryonNumber()>0) thePDG=theQC.GetZNSPDGCode();
    39624110  G4LorentzVector q4M = qH->Get4Momentum(); // Get 4-momentum of theTotalNucleus
    3963   if(theBN<1 || thePDG<80000000 || thePDG==90000000) // Hadron, anti-nucleous, or vacuum
     4111  if(!theBN || thePDG<80000000 || thePDG==90000000) // Hadron, anti-nucleous, or vacuum
    39644112  {
    39654113#ifdef debug
    … …  
    39694117    {
    39704118#ifdef qdebug
    3971       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (1) qH="<<G4endl;
    3972       else
     4119      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (1) qH=0"<<G4endl;
    39734120#endif
    39744121      delete qH;
    3975 #ifdef qdebug
    3976       qH=0;                                           // @Not necessary@(can't be checked)
    3977 #endif
    39784122      G4cout<<"-Warning-G4QNucleus::EvapNuc:vacuum,4Mom="<<q4M<<G4endl;
    39794123    }
    39804124    else   // Put input to the output (delete equivalent)
    39814125    {
     4126      G4cout<<"-Warning-G4QNucleus::EvapNuc:vacuum, Called for Meson PDG="<<thePDG<<G4endl;
    39824127      evaHV->push_back(qH);
    3983 #ifdef qdebug
    3984       qH=0;                                           // @Not necessary@(can't be checked)
    3985 #endif
    39864128    }
    39874129    return;
    … …  
    40054147#endif
    40064148    evaHV->push_back(qH);     // TheFundamentalParticles must be FilledAsTheyAre (del.eq)
    4007 #ifdef qdebug
    4008       qH=0;                                           // @Not necessary@(can't be checked)
    4009 #endif
    40104149    return;
    40114150  }
    … …  
    40194158        <<theBN<<", Z="<<theC<<", N="<<theN<<", S="<<theS<<G4endl;
    40204159#endif
    4021   if(thePDG==91000000||thePDG==90001000||thePDG==90000001)
     4160  if(theBN<-2)
     4161  {
     4162    G4cout<<"-Warning-G4QNuc::EvapNuc: Evapor of anti-nuclei is not implemented yet PDG="
     4163          <<thePDG<<G4endl;
     4164    evaHV->push_back(qH);
     4165    return;
     4166  }
     4167  else if(thePDG==91000000||thePDG==90001000||thePDG==90000001) // Excited Lambda* or N*
    40224168  //else if(2>3)// One can easily close this decay as it will be done later (time of calc?)
    40234169  {
    … …  
    40314177#endif
    40324178      evaHV->push_back(qH); // (delete equivalent)
     4179      return;
     4180    }
     4181    else if(totMass<gsM)
     4182    {
    40334183#ifdef qdebug
    4034       qH=0;                                           // @Not necessary@(can't be checked)
    4035 #endif
    4036       return;
    4037     }
    4038     else if(totMass<gsM)
    4039     {
    4040 #ifdef qdebug
    4041       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (2) qH="<<G4endl;
    4042       else
     4184      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (2) qH=0"<<G4endl;
    40434185#endif
    40444186      delete qH;
    4045 #ifdef qdebug
    4046       qH=0;                                              // @Not necessary@
    4047 #endif
    40484187      G4cerr<<"***G4QN::EvaNuc:Baryon "<<thePDG<<" is belowMassShell, M="<<totMass<<G4endl;
    40494188      throw G4QException("G4QNucleus::EvaporateNucleus: Baryon is below Mass Shell");
    … …  
    40594198      if(d>142.)                           // @@ to avoid more precise calculations
    40604199      {
    4061         if(thePDG==90001000)               // D+ -> n + pi+
     4200        if(thePDG==90001000)               // p* -> n + pi+
    40624201        {
    40634202          gsM=mNeut;
    … …  
    40664205          decM=mPi;
    40674206        }
    4068         else if(thePDG==90000001)          // D+ -> n + pi+
     4207        else if(thePDG==90000001)          // n* -> p + pi-
    40694208        {
    40704209          gsM=mProt;
    … …  
    40734212          decM=mPi;
    40744213        }
    4075         else
     4214        else                               // decay in Pi0
    40764215        {
    40774216          decPDG=111;
    … …  
    40844223      {
    40854224#ifdef qdebug
    4086         if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (3) qH="<<G4endl;
    4087         else
     4225        if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (3) qH=0"<<G4endl;
    40884226#endif
    40894227        delete qH;
    4090 #ifdef qdebug
    4091         qH=0;                                              // @Not necessary@
    4092 #endif
    40934228        G4cerr<<"***G4QNuc::EvaNuc:h="<<thePDG<<"(GSM="<<gsM<<")+gam>tM="<<totMass<<G4endl;
    40944229        throw G4QException("G4QNucleus::EvaporateNucleus:BaryonDecay In Baryon+Gam Error");
    … …  
    41094244      evaHV->push_back(curG);         // Fill gamma/pion (delete equivalent)
    41104245#ifdef qdebug
    4111       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (4) qH="<<G4endl;
    4112       else
     4246      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (4) qH=0"<<G4endl;
    41134247#endif
    41144248      delete qH;
     4249    }
     4250  }
     4251  else if(thePDG==89000000||thePDG==89999000||thePDG==89999999) // anti-Lambda* or anti-N*
     4252  //else if(2>3)// One can easily close this decay as it will be done later (time of calc?)
     4253  {
     4254    G4double gsM=mNeut;
     4255    if(thePDG==89999000)      gsM=mProt;
     4256    else if(thePDG==89000000) gsM=mLamb;
     4257    if(fabs(totMass-gsM)<.001)
     4258    {
     4259#ifdef debug
     4260      G4cout<<"G4QNu::EvaNucl:(aB*),GSM="<<gsM<<", H="<<thePDG<<qH->Get4Momentum()<<G4endl;
     4261#endif
     4262      evaHV->push_back(qH); // (delete equivalent)
     4263      return;
     4264    }
     4265    else if(totMass<gsM)
     4266    {
    41154267#ifdef qdebug
    4116       qH=0;                                              // @Not necessary@
    4117 #endif
     4268      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (2a) qH=0"<<G4endl;
     4269#endif
     4270      delete qH;
     4271      G4cerr<<"*G4QN::EvaNuc:antiBaryon="<<thePDG<<" below MassShell, M="<<totMass<<G4endl;
     4272      throw G4QException("G4QNucleus::EvaporateNucleus: anti-Baryon is below Mass Shell");
     4273    }
     4274    else                                 // Decay in gamma or charged pion (@@ neutral)
     4275    {
     4276      G4double d=totMass-gsM;
     4277#ifdef debug
     4278      G4cout<<"G4QN::EvaporNucl: PDG="<<thePDG<<",M="<<totMass<<">"<<gsM<<",d="<<d<<G4endl;
     4279#endif
     4280      G4int decPDG=22;
     4281      G4double decM=0.;
     4282      if(d>142.)                           // @@ to avoid more precise calculations
     4283      {
     4284        if(thePDG==89999000)               // anti (p* -> n + pi+)
     4285        {
     4286          gsM=mNeut;
     4287          thePDG=89999999;
     4288          decPDG=-211;
     4289          decM=mPi;
     4290        }
     4291        else if(thePDG==89999999)          // anti (n* -> p + pi-)
     4292        {
     4293          gsM=mProt;
     4294          thePDG=89999000;
     4295          decPDG=211;
     4296          decM=mPi;
     4297        }
     4298        else                               // decay in Pi0
     4299        {
     4300          decPDG=111;
     4301          decM=mPi0;
     4302        }
     4303      }
     4304      G4LorentzVector h4Mom(0.,0.,0.,gsM); // GSMass must be updated in previous while-LOOP
     4305      G4LorentzVector g4Mom(0.,0.,0.,decM);
     4306      if(!G4QHadron(q4M).DecayIn2(h4Mom, g4Mom))
     4307      {
     4308#ifdef qdebug
     4309        if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (3a) qH=0"<<G4endl;
     4310#endif
     4311        delete qH;
     4312        G4cerr<<"**G4QNuc::EvaNuc:ah="<<thePDG<<"(GSM="<<gsM<<")+gam>tM="<<totMass<<G4endl;
     4313        throw G4QException("G4QNucleus::EvaporateNucleus:BaryonDecay In Baryon+Gam Error");
     4314      }
     4315#ifdef debug
     4316      G4cout<<"G4QNucl::EvaNuc:aM="<<totMass<<q4M<<"->"<<thePDG<<h4Mom<<"+g="<<g4Mom<<",n="
     4317            <<evaHV->size()<<G4endl;
     4318#endif
     4319      G4QHadron* curH = new G4QHadron(thePDG, h4Mom);
     4320#ifdef debug
     4321      G4cout<<"G4QNucleus::EvaporateNucleus: antiHadr="<<thePDG<<h4Mom<<G4endl;
     4322#endif
     4323      evaHV->push_back(curH);         // Fill Baryon (delete equiv.)
     4324      G4QHadron* curMes = new G4QHadron(decPDG, g4Mom);
     4325#ifdef debug
     4326      G4cout<<"G4QNucleus::EvaporateNucleus: (anti) Gamma(pion)4M="<<g4Mom<<G4endl;
     4327#endif
     4328      evaHV->push_back(curMes);         // Fill gamma/pion (delete equivalent)
     4329#ifdef qdebug
     4330      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (4a) qH=0"<<G4endl;
     4331#endif
     4332      delete qH;
    41184333    }
    41194334  }
    … …  
    41424357      evaHV->push_back(curM); // (delete equivalent)
    41434358#ifdef qdebug
    4144       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (5) qH="<<G4endl;
    4145       else
    4146 #endif
    4147       delete qH;                                         // *** New correction ***
     4359      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (5) qH=0"<<G4endl;
     4360#endif
     4361      delete qH;
     4362      return;
     4363    }
     4364    else if(totMass<gsM+mPi)
     4365    {
    41484366#ifdef qdebug
    4149       qH=0;                                              // @Not necessary@
    4150 #endif
    4151       return;
    4152     }
    4153     else if(totMass<gsM+mPi)
    4154     {
    4155 #ifdef qdebug
    4156       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (6) qH="<<G4endl;
    4157       else
     4367      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (6) qH=0"<<G4endl;
    41584368#endif
    41594369      delete qH;
    4160 #ifdef qdebug
    4161       qH=0;                                              // @Not necessary@
    4162 #endif
    41634370      G4cerr<<"***G4QN::EvaNuc:Delta "<<thePDG<<" is belowMassShell, M="<<totMass<<G4endl;
    41644371      throw G4QException("G4QNucleus::EvaporateNucleus: Delta is below Mass Shell");
    … …  
    41714378      {
    41724379#ifdef qdebug
    4173         if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (7) qH="<<G4endl;
    4174         else
     4380        if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (7) qH=0"<<G4endl;
    41754381#endif
    41764382        delete qH;
    4177 #ifdef qdebug
    4178         qH=0;                                              // @Not necessary@
    4179 #endif
    41804383        G4cerr<<"***G4QNuc::EvaNuc:Dh="<<thePDG<<"N+pi="<<gsM+mPi<<">tM="<<totMass<<G4endl;
    41814384        throw G4QException("G4QNucleus::EvaporateNucleus: DeltaDecInBaryon+Pi Error");
    … …  
    41964399      evaHV->push_back(curG);         // Fill the pion (delete equivalent)
    41974400#ifdef qdebug
    4198       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (8) qH="<<G4endl;
    4199       else
     4401      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (8) qH=0"<<G4endl;
    42004402#endif
    42014403      delete qH;
     4404    }
     4405  }
     4406  else if((thePDG==89998001||thePDG==90000998) && totMass>1080.) // @@ to avoid threshold
     4407  //else if(2>3)// One can easily close this decay as it will be done later (time of calc?)
     4408  {
     4409    G4double gsM=mNeut;
     4410    G4int barPDG=-2112;
     4411    G4int mesPDG=211;
     4412    if(thePDG==89998001)
     4413    {
     4414      gsM=mProt;
     4415      barPDG=-2212;
     4416      mesPDG=-211;
     4417    }
     4418    if(fabs(totMass-gsM-mPi)<.001)
     4419    {
     4420#ifdef debug
     4421      G4cout<<"G4QN::EvaporateNuc:(A)GSM="<<gsM<<",H="<<thePDG<<qH->Get4Momentum()<<G4endl;
     4422#endif
     4423      G4LorentzVector h4Mom=q4M*(gsM/totMass);           // At rest in CM
     4424      G4QHadron* curB = new G4QHadron(barPDG,h4Mom);
     4425      evaHV->push_back(curB); // (delete equivalent)
     4426      G4LorentzVector g4Mom=q4M*(mPi/totMass);
     4427      G4QHadron* curM = new G4QHadron(mesPDG,g4Mom);
     4428      evaHV->push_back(curM); // (delete equivalent)
    42024429#ifdef qdebug
    4203       qH=0;
    4204 #endif
    4205     }
    4206   }
    4207   else if(theBN>0&&theS<0)
    4208   {
    4209     DecayAntiStrange(qH,evaHV);
     4430      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (5a) qH=0"<<G4endl;
     4431#endif
     4432      delete qH;
     4433      return;
     4434    }
     4435    else if(totMass<gsM+mPi)
     4436    {
    42104437#ifdef qdebug
    4211     qH=0;                                              // @Not necessary@
    4212 #endif
    4213   }// "AntyStrangeNucleus" (del.eq.)
    4214   else if(theBN>0&&(theC<0||theC>theBN-theS))
    4215   {
    4216     DecayIsonucleus(qH,evaHV);
     4438      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (6a) qH=0"<<G4endl;
     4439#endif
     4440      delete qH;
     4441      G4cerr<<"***G4QN::EvaNuc:aDelta "<<thePDG<<" is belowMassShell, M="<<totMass<<G4endl;
     4442      throw G4QException("G4QNucleus::EvaporateNucleus: anti-Delta is below Mass Shell");
     4443    }
     4444    else                                 // Decay in gamma or charged pion (@@ neutral)
     4445    {
     4446      G4LorentzVector h4Mom(0.,0.,0.,gsM); // GSMass must be updated in previous while-LOOP
     4447      G4LorentzVector g4Mom(0.,0.,0.,mPi);
     4448      if(!G4QHadron(q4M).DecayIn2(h4Mom, g4Mom))
     4449      {
    42174450#ifdef qdebug
    4218       qH=0;                                              // @Not necessary@
    4219 #endif
    4220   }//"UnavoidIsonucl"
    4221   else if(theBN==2)
    4222   {
    4223     DecayDibaryon(qH,evaHV);
     4451        if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (7a) qH=0"<<G4endl;
     4452#endif
     4453        delete qH;
     4454        G4cerr<<"***G4QNuc::EvaNuc:aD="<<thePDG<<"N+pi="<<gsM+mPi<<">tM="<<totMass<<G4endl;
     4455        throw G4QException("G4QNucleus::EvaporateNucleus:AntiDeltaDecayInBaryon+Pi Error");
     4456      }
     4457#ifdef debug
     4458      G4cout<<"G4QNuc::EvaNuc:(aD) "<<totMass<<q4M<<"->"<<thePDG<<h4Mom<<" + pi="<<g4Mom
     4459            <<", nH="<<evaHV->size()<<G4endl;
     4460#endif
     4461      G4QHadron* curH = new G4QHadron(barPDG,h4Mom);
     4462#ifdef debug
     4463      G4cout<<"G4QNucleus::EvaporateNucl: Nucleon="<<thePDG<<h4Mom<<G4endl;
     4464#endif
     4465      evaHV->push_back(curH);         // Fill the nucleon (delete equiv.)
     4466      G4QHadron* curMes = new G4QHadron(mesPDG,g4Mom);
     4467#ifdef debug
     4468      G4cout<<"G4QE::EvaporateR: Pion="<<g4Mom<<G4endl;
     4469#endif
     4470      evaHV->push_back(curMes);         // Fill the pion (delete equivalent)
    42244471#ifdef qdebug
    4225       qH=0;                                              // @Not necessary@
    4226 #endif
    4227   }    //=> "Dibaryon" case (del eq.)
    4228   else if(thePDG==89999003||thePDG==90002999)   //=> "ISO-dibarion" (analyse and decay)
     4472      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (8a) qH=0"<<G4endl;
     4473#endif
     4474      delete qH;
     4475    }
     4476  }
     4477  else if(theBN>0&&theS<0) DecayAntiStrange(qH,evaHV); // "AntyStrangeNucleus" (del.eq.)
     4478  else if(theBN>0&&(theC<0||theC>theBN-theS))DecayIsonucleus(qH,evaHV);//"Isonucleus"(d.e.)
     4479  else if((thePDG==89999003 || thePDG==90002999) && totMass>2020.) //=> "ISO-dibarion"
    42294480  {
    42304481    // @@ Check that it never comes here !!
    … …  
    42464497      {
    42474498#ifdef qdebug
    4248         if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (9) qH="<<G4endl;
    4249         else
     4499        if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (9) qH=0"<<G4endl;
    42504500#endif
    42514501        delete qH;
    4252 #ifdef qdebug
    4253         qH=0;                                              // @Not necessary@
    4254 #endif
    42554502        G4cerr<<"***G4QNucleus::EvaporateNucleus: tM="<<totMass<<"-> 2N="<<nucPDG<<"(M="
    42564503              <<nucM<<") + pi="<<piPDG<<"(M="<<mPi<<")"<<G4endl;
    … …  
    42584505      }
    42594506#ifdef qdebug
    4260       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (10) qH="<<G4endl;
    4261       else
     4507      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (10) qH=0"<<G4endl;
    42624508#endif
    42634509      delete qH;
    4264 #ifdef qdebug
    4265       qH=0;                                              // @Not necessary@
    4266 #endif
    42674510      G4QHadron* h1H = new G4QHadron(nucPDG,n14M);
    42684511#ifdef debug
    … …  
    42844527    {
    42854528#ifdef qdebug
    4286       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (11) qH="<<G4endl;
    4287       else
     4529      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (11) qH=0"<<G4endl;
    42884530#endif
    42894531      delete qH;
    4290 #ifdef qdebug
    4291       qH=0;                                              // @Not necessary@
    4292 #endif
    42934532      G4cerr<<"***G4QNucleus::EvapNucleus: IdPDG="<<thePDG<<", q4M="<<q4M<<", M="<<totMass
    42944533            <<" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<G4endl;
    … …  
    42964535    }
    42974536  }
    4298   else if((thePDG==90000002||thePDG==90001001||thePDG==90002000)&&totMass>2020.) //=> IsoDB
     4537  else if((thePDG==90000002||thePDG==90001001||thePDG==90002000)&&totMass>2020.) //=> IsoBP
    42994538  {
    43004539    // @@ Pi0 can be taken into account !
    … …  
    43334572      }
    43344573#ifdef qdebug
    4335       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (12) qH="<<G4endl;
    4336       else
     4574      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (12) qH=0"<<G4endl;
    43374575#endif
    43384576      delete qH;
    4339 #ifdef qdebug
    4340       qH=0;                                              // @Not necessary@
    4341 #endif
    43424577      G4QHadron* h1H = new G4QHadron(n1PDG,n14M);
    43434578#ifdef debug
    … …  
    43594594    {
    43604595#ifdef qdebug
    4361       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (13) qH="<<G4endl;
    4362       else
     4596      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (13) qH=0"<<G4endl;
    43634597#endif
    43644598      delete qH;
    4365 #ifdef qdebug
    4366       qH=0;                                              // @Not necessary@
    4367 #endif
    43684599      G4cerr<<"***G4QNuc::EvaporateNucleus: IdPDG="<<thePDG<<", q4M="<<q4M<<", M="<<totMass
    43694600            <<" < M1+M2+Pi, d="<<totMass-n1M-n2M-mPi<<G4endl;
    … …  
    43714602    }
    43724603  }
     4604  else if(theBN==2) DecayDibaryon(qH, evaHV);       //=> "Dibaryon" case (del eq.)
     4605  else if((thePDG==90000997 || thePDG==89997001) && totMass>2020.) //=> "anti-ISO-dibarion"
     4606  {
     4607    // @@ Check that it never comes here !!
     4608    G4int  nucPDG = -2112;
     4609    G4double nucM = mNeut;
     4610    G4int   piPDG = 211;
     4611    if(thePDG==90002999)
     4612    {
     4613      nucPDG = -2212;
     4614      nucM   = mProt;
     4615      piPDG  = -211;
     4616    }
     4617    if(totMass>mPi+nucM+nucM)
     4618    {
     4619      G4LorentzVector n14M(0.,0.,0.,nucM);
     4620      G4LorentzVector n24M(0.,0.,0.,nucM);
     4621      G4LorentzVector pi4M(0.,0.,0.,mPi);
     4622      if(!G4QHadron(q4M).DecayIn3(n14M,n24M,pi4M))
     4623      {
     4624#ifdef qdebug
     4625        if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (9a) qH=0"<<G4endl;
     4626#endif
     4627        delete qH;
     4628        G4cerr<<"***G4QNucleus::EvaporateNucleus:antiM="<<totMass<<"-> 2aN="<<nucPDG<<"(M="
     4629              <<nucM<<") + pi="<<piPDG<<"(M="<<mPi<<")"<<G4endl;
     4630        throw G4QException("G4QNucleus::EvaporateNucleus:Anti-ISO-DibaryonDecayIn3 error");
     4631      }
     4632#ifdef qdebug
     4633      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (10a) qH=0"<<G4endl;
     4634#endif
     4635      delete qH;
     4636      G4QHadron* h1H = new G4QHadron(nucPDG,n14M);
     4637#ifdef debug
     4638      G4cout<<"G4QNucleus::EvaporateNucleus:(I) antiBar1="<<nucPDG<<n14M<<G4endl;
     4639#endif
     4640      evaHV->push_back(h1H);                // (delete equivalent)
     4641      G4QHadron* h2H = new G4QHadron(nucPDG,n24M);
     4642#ifdef debug
     4643      G4cout<<"G4QNucleus::EvaporateNucleus:(I) antiBar2="<<nucPDG<<n24M<<G4endl;
     4644#endif
     4645      evaHV->push_back(h2H);                // (delete equivalent)
     4646      G4QHadron* piH = new G4QHadron(piPDG,pi4M);
     4647#ifdef debug
     4648      G4cout<<"G4QNucleus::EvaporateNucleus:(I) (anti) Pi="<<piPDG<<pi4M<<G4endl;
     4649#endif
     4650      evaHV->push_back(piH);                // (delete equivalent)
     4651    }
     4652    else
     4653    {
     4654#ifdef qdebug
     4655      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (11a) qH=0"<<G4endl;
     4656#endif
     4657      delete qH;
     4658      G4cerr<<"***G4QNucleus::EvapNucleus: aIdPDG="<<thePDG<<", q4M="<<q4M<<", M="<<totMass
     4659            <<" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<G4endl;
     4660      throw G4QException("G4QNucleus::EvaporateNucleus:AntiISODiBaryon is underMassShell");
     4661    }
     4662  }
     4663  else if((thePDG==89999998||thePDG==89998999||thePDG==89998000)&&totMass>2020.)//=>AnIsoBP
     4664  {
     4665    // @@ Pi0 can be taken into account !
     4666    G4int  n1PDG = -2212;
     4667    G4int  n2PDG = -2112;
     4668    G4int  piPDG = 211;                           // dummy initialization
     4669    G4double n1M = mProt;
     4670    G4double n2M = mNeut;
     4671    if     (thePDG==89998000) piPDG  = -211;      // anti ( pp -> np + pi- )
     4672    else if(thePDG==89999998) piPDG  =  211;      // anti ( nn -> np + pi- )
     4673    else                                          // anti ( np -> 50%(nnpi+) 50%(pppi-) )
     4674    {
     4675      if(G4UniformRand()>.5)
     4676      {
     4677        n1PDG=-2112;
     4678        n1M=mNeut;
     4679        piPDG  = -211;
     4680      }
     4681      else
     4682      {
     4683        n2PDG=-2212;
     4684        n2M=mProt;
     4685        piPDG  =  211;
     4686      }
     4687    }
     4688    if(totMass>mPi+n1M+n2M)
     4689    {
     4690      G4LorentzVector n14M(0.,0.,0.,n1M);
     4691      G4LorentzVector n24M(0.,0.,0.,n2M);
     4692      G4LorentzVector pi4M(0.,0.,0.,mPi);
     4693      if(!G4QHadron(q4M).DecayIn3(n14M,n24M,pi4M))
     4694      {
     4695        G4cerr<<"**G4QNucl::EvapNucleus:IsoDN,antiM="<<totMass<<"-> N1="<<n1PDG<<"(M="<<n1M
     4696              <<") + N2="<<n2PDG<<"(M="<<n2M<<") + pi="<<piPDG<<" (Mpi="<<mPi<<")"<<G4endl;
     4697        throw G4QException("G4QNucl::EvaporateNucleus:AntiExcitedDibaryon DecayIn3 error");
     4698      }
     4699#ifdef qdebug
     4700      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (12a) qH=0"<<G4endl;
     4701#endif
     4702      delete qH;
     4703      G4QHadron* h1H = new G4QHadron(n1PDG,n14M);
     4704#ifdef debug
     4705      G4cout<<"G4QNucleus::EvaporateNucleus: antiBar1="<<n1PDG<<n14M<<G4endl;
     4706#endif
     4707      evaHV->push_back(h1H);                // (delete equivalent)
     4708      G4QHadron* h2H = new G4QHadron(n2PDG,n24M);
     4709#ifdef debug
     4710      G4cout<<"G4QNucleus::EvaporateNucleus: antiBar2="<<n2PDG<<n24M<<G4endl;
     4711#endif
     4712      evaHV->push_back(h2H);                // (delete equivalent)
     4713      G4QHadron* piH = new G4QHadron(piPDG,pi4M);
     4714#ifdef debug
     4715      G4cout<<"G4QNucleus::EvaporateNucleus: (anti)Pi="<<piPDG<<pi4M<<G4endl;
     4716#endif
     4717      evaHV->push_back(piH);                // (delete equivalent)
     4718    }
     4719    else
     4720    {
     4721#ifdef qdebug
     4722      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (13a) qH=0"<<G4endl;
     4723#endif
     4724      delete qH;
     4725      G4cerr<<"***G4QNuc::EvaporateNucleus:andPDG="<<thePDG<<", q4M="<<q4M<<", M="<<totMass
     4726            <<" < M1+M2+Pi, d="<<totMass-n1M-n2M-mPi<<G4endl;
     4727      throw G4QException("G4QNucleus::EvaporateNucleus:AntiDiBarState is under MassShell");
     4728    }
     4729  }
     4730  else if(theBN==-2) DecayAntiDibaryon(qH,evaHV);       //=> "Anti-Dibaryon" case (del eq.)
    43734731  else if(fabs(totMass-totGSM)<.001)  // Fill as it is or decay Be8, He5, Li5 (@@ add more)
    43744732  {
    … …  
    43794737    {
    43804738      DecayAlphaAlpha(qH,evaHV);
    4381 #ifdef qdebug
    4382       qH=0;                                              // @Not necessary@
    4383 #endif
    43844739    } // "Alpha+Alpha Decay" case (del eq.)
    43854740    else if(thePDG==90004002)
    43864741    {
    43874742      DecayAlphaDiN(qH,evaHV);
    4388 #ifdef qdebug
    4389       qH=0;                                              // @Not necessary@
    4390 #endif
    43914743    } // Decay alpha+2p(alpha+2n is stable)
    43924744    else if((theC==theBN||theN==theBN||theS==theBN)&&theBN>1)
    43934745    {
    43944746      DecayMultyBaryon(qH,evaHV);
    4395 #ifdef qdebug
    4396       qH=0;                                              // @Not necessary@
    4397 #endif
    43984747    }
    43994748    else if(theBN==5)
    44004749    {
    44014750      DecayAlphaBar(qH,evaHV);
    4402 #ifdef qdebug
    4403       qH=0;                                              // @Not necessary@
    4404 #endif
    44054751    }   // Decay unstable A5 system (del eq.)
    44064752    else
    44074753    {
    44084754      evaHV->push_back(qH);
    4409 #ifdef qdebug
    4410       qH=0;
    4411 #endif
    44124755    }      // Fill as it is (del eq.)
    44134756  }
    4414   else if(theBN>1&&thePDG>88000000&&thePDG<89000000)//==> 2antiK in the nucleus (!Comment!)
     4757  else if(theBN>1 && thePDG>88000000 && thePDG<89000000) //==> 2antiK in the nucleus
    44154758  {
    44164759    G4cout<<"---Warning---G4QNucl::EvaNuc:MustNotBeHere.PDG="<<thePDG<<",S="<<theS<<G4endl;
    … …  
    44444787    {
    44454788#ifdef qdebug
    4446       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (14) qH="<<G4endl;
    4447       else
     4789      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (14) qH=0"<<G4endl;
    44484790#endif
    44494791      delete qH;
    4450 #ifdef qdebug
    4451       qH=0;                                              // @Not necessary@
    4452 #endif
    44534792      G4cout<<"***G4QNucleus::EvaNuc:tM="<<totMass<<"-> N="<<nucPDG<<"(M="<<nucM<<") + k1="
    44544793            <<k1PDG<<"(M="<<mK1<<") + k2="<<k2PDG<<"(M="<<mK2<<")"<<G4endl;
    … …  
    44564795    }
    44574796#ifdef qdebug
    4458     if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (15) qH="<<G4endl;
    4459     else
     4797    if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (15) qH=0"<<G4endl;
    44604798#endif
    44614799    delete qH;
    4462 #ifdef qdebug
    4463     qH=0;                                              // @Not necessary@
    4464 #endif
    44654800    G4QHadron* k1H = new G4QHadron(k1PDG,k14M);
    44664801#ifdef debug
    … …  
    44924827    G4double GSMass =qNuc.GetGSMass();         // GSMass of the Total Residual Nucleus
    44934828    G4QContent totQC=qNuc.GetQCZNS();          // QuarkCont of theTotalResidNucleus (theQC)
    4494     G4int    bA     =qNuc.GetA();              // A#of baryons in theTotal Residual Nucleus
    4495     G4int    bZ     =qNuc.GetZ();              // A#of protons in theTotal Residual Nucleus
    4496     G4int    bN     =qNuc.GetN();              // A#of neutrons in theTotal ResidualNucleus
    4497     G4int    bS     =qNuc.GetS();              // A#of lambdas in theTotal Residual Nucleus
     4829    G4int    bA     =qNuc.GetA();              // A#of baryons in Total Residual Nucleus
     4830    G4int    bZ     =qNuc.GetZ();              // A#of protons in the Total ResidualNucleus
     4831    G4int    bN     =qNuc.GetN();              // A#of neutrons in the TotalResidualNucleus
     4832#ifdef ppdebug
     4833    G4cout<<"G4QN::EvaNuc: theBN="<<theBN<<", bA="<<bA<<", bZ="<<bZ<<", bN="<<bN<<G4endl;
     4834#endif
     4835    G4int    bS     =qNuc.GetS();              // A#of lambdas in the Total ResidualNucleus
    44984836#ifdef debug
    44994837    if(bZ==2&&bN==5)G4cout<<"G4QNucleus::EvaNucl: (2,5) GSM="<<GSMass<<" > "
    … …  
    45164854#endif
    45174855      evaHV->push_back(qH);
    4518 #ifdef qdebug
    4519       qH=0;                                     // @Not necessary@ (Can't be checked)
    4520 #endif
    45214856      return;
    45224857    }
    … …  
    45864921              delete theLast; //When kill, DON'T forget to delete lastQHadron asAnInstance!
    45874922#ifdef qdebug
    4588               if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (16) qH="<<G4endl;
    4589               else
     4923              if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (16) qH=0"<<G4endl;
    45904924#endif
    45914925              delete qH;
    4592 #ifdef qdebug
    4593               qH=0;                                              // @Not necessary@
    4594 #endif
    45954926#ifdef debug
    45964927              G4cout<<"G4QNucleus::EvaporateNucl: EVH "<<totPDG<<q4M<<" fill AsIs"<<G4endl;
    … …  
    46074938                delete theLast; //When kill,DON'T forget to delete lastQHadron asAnInstance
    46084939#ifdef qdebug
    4609                 if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (17) qH="<<G4endl;
    4610                 else
     4940                if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (17) qH=0"<<G4endl;
    46114941#endif
    46124942                delete qH;
    4613 #ifdef qdebug
    4614                 qH=0;                                              // @Not necessary@
    4615 #endif
    46164943#ifdef debug
    46174944                G4cout<<"***G4QNucleus::EvaNucl: EVH "<<totPDG<<q4M<<" fill AsIs"<<G4endl;
    … …  
    46264953                delete evH;
    46274954#ifdef qdebug
    4628                 if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (18) qH="<<G4endl;
    4629                 else
     4955                if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (18) qH=0"<<G4endl;
    46304956#endif
    46314957                delete qH;
    4632 #ifdef qdebug
    4633                 qH=0;                                              // @Not necessary@
    4634 #endif
    46354958                theLast->Set4Momentum(last4M);// Already exists:don't create&fill,->set4Mom
    46364959                G4QHadron* nucH = new G4QHadron(thePDG,r4Mom); // Create QHadron for qH-nuc
    … …  
    46614984        {
    46624985#ifdef qdebug
    4663         if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (19) qH="<<G4endl;
    4664         else
     4986          if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (19) qH=0"<<G4endl;
    46654987#endif
    46664988          delete qH;
    4667 #ifdef qdebug
    4668           qH=0;                                              // @Not necessary@
    4669 #endif
    46704989          G4cerr<<"**G4QN::EvaNuc:h="<<thePDG<<"(GSM="<<GSMass<<")+g>tM="<<totMass<<G4endl;
    46714990          throw G4QException("G4QNucleus::EvaporateNucleus: Decay in Gamma failed");
    … …  
    46865005        evaHV->push_back(curG);       // Fill the gamma (delete equivalent)
    46875006#ifdef qdebug
    4688         if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (20) qH="<<G4endl;
    4689         else
     5007        if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (20) qH=0"<<G4endl;
    46905008#endif
    46915009        delete qH;
    4692 #ifdef qdebug
    4693         qH=0;                                              // @Not necessary@
    4694 #endif
    4695       }
    4696     }
    4697     else if(bA==2)
    4698     {
    4699       DecayDibaryon(qH,evaHV);
    4700 #ifdef qdebug
    4701       qH=0;                                              // @Not necessary@
    4702 #endif
    4703     }// Decay the residual dibaryon (del.equivalent)
    4704     else if(bA>0&&bS<0)
    4705     {
    4706       DecayAntiStrange(qH,evaHV);
    4707 #ifdef qdebug
    4708       qH=0;                                              // @Not necessary@
    4709 #endif
    4710     }// Decay nucleus with antistrangeness
     5010      }
     5011    }
     5012    else if(bA>0&&bS<0) DecayAntiStrange(qH,evaHV);// Decay nucleus with antistrangeness
     5013    else if(bA==2) DecayDibaryon(qH,evaHV); // Decay the residual dibaryon (del.equivalent)
     5014    else if(bA==-2) DecayAntiDibaryon(qH,evaHV);   // Decay residual anti-dibaryon (del.eq)
    47115015    else if(totMass<GSMass+.003&&(bsCond||dbsCond))//==>" M<GSM but decay is possible" case
    47125016    {
    4713 #ifdef debug
     5017#ifdef pdebug
    47145018      G4cout<<"G4QN::EvN:2B="<<dbsCond<<",B="<<bsCond<<",M="<<totMass<<"<"<<GSMass<<G4endl;
    47155019#endif
    … …  
    48585162            <<",p="<<totMass-pResM-mProt<<",l="<<totMass-lResM-mLamb<<G4endl;
    48595163#endif
    4860 
    4861       if ( thePDG == 90004004 ||
    4862 
    4863           (thePDG == 90002004 && totMass > mHel6+.003) ||
    4864 
    4865           (bA > 4 && bsCond && bN > 1 && bZ > 1 && totMass > aResM+mAlph) ||
    4866 
     5164      if ( thePDG == 90004004                                                 ||
     5165          (thePDG == 90002004 && totMass > mHel6+.003)                        ||
     5166          (bA > 4 && bsCond && bN > 1 && bZ > 1 && totMass > aResM+mAlph)     ||
    48675167          (bA > 1 && bsCond && ( (bN > 0 && totMass > nResM+mNeut) ||
    48685168                                 (bZ > 0 && totMass > pResM+mProt) ||
    4869                                  (bS > 0 && totMass > lResM+mLamb) ) )    ||
     5169                                 (bS > 0 && totMass > lResM+mLamb) ) )        ||
    48705170          (bA > 2 &&
    48715171           (( bN > 0 && bZ > 0 &&
    … …  
    50015301          {
    50025302            evaHV->push_back(qH);     // Fill as it is (delete equivalent)
    5003 #ifdef qdebug
    5004             qH=0;                           // @Not necessary@ (Can't be checked)
    5005 #endif
    50065303            G4cout<<"---Warning---G4QNucleus::EvaNuc:rP="<<pResPDG<<",rN="<<nResPDG<<",rL="
    50075304                  <<lResPDG<<",N="<<bN<<",Z="<<bZ<<",L="<<bS<<",totM="<<totMass<<",n="
    … …  
    50145311          {
    50155312#ifdef qdebug
    5016             if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (21) qH="<<G4endl;
    5017             else
     5313            if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (21) qH=0"<<G4endl;
    50185314#endif
    50195315            delete qH;
    5020 #ifdef qdebug
    5021             qH=0;
    5022 #endif
    50235316            G4QHadron* HadrB = new G4QHadron(barPDG,a4Mom);
    50245317#ifdef debug
    … …  
    50415334          {
    50425335            evaHV->push_back(qH);    // Fill as it is (delete equivalent)
    5043 #ifdef qdebug
    5044             qH=0;                                // @Not necessary@ (Can't be checked)
    5045 #endif
    50465336            G4cout<<"---Warning---G4QN::EvN:rNN="<<nnResPDG<<",rNP="<<npResPDG<<",rPP="
    50475337                  <<ppResPDG<<",N="<<bN<<",Z="<<bZ<<",L="<<bS<<",tM="<<totMass<<",nn="
    … …  
    50545344          {
    50555345#ifdef qdebug
    5056             if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (22) qH="<<G4endl;
    5057             else
     5346            if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (22) qH=0"<<G4endl;
    50585347#endif
    50595348            delete qH;
    5060 #ifdef qdebug
    5061             qH=0;                                              // @Not necessary@
    5062 #endif
    50635349            G4QHadron* HadrB = new G4QHadron(barPDG,a4Mom);
    50645350#ifdef debug
    … …  
    50875373#endif
    50885374        evaHV->push_back(qH);  // FillAsItIs (del.eq.)
    5089 #ifdef qdebug
    5090         qH=0;                                         // @Not necessary@ (Can't be checked)
    5091 #endif
    50925375        return;
    50935376      }
    … …  
    50995382#endif
    51005383        evaHV->push_back(qH);                   // Correct or fill as it is
    5101 #ifdef qdebug
    5102         qH=0;                                   // @Not necessary@ (can't be checked)
    5103 #endif
    51045384        return;
    51055385      }
    51065386    }
    5107     else                                        // ===> Evaporation of excited system
    5108     {
    5109 #ifdef debug
     5387    else                                        // ===> Evaporation of the excited system
     5388    {
     5389#ifdef pdebug
    51105390      G4cout<<"G4QN::EvaNuc:***EVA***tPDG="<<thePDG<<",M="<<totMass<<">GSM="<<GSMass<<",d="
    51115391            <<totMass-GSMass<<", N="<<qNuc.Get4Momentum()<<qNuc.Get4Momentum().m()<<G4endl;
    … …  
    51415421#endif
    51425422          evaHV->push_back(qH);               // fill AsItIs
    5143 #ifdef qdebug
    5144           qH=0;                               // @Not necessary@ (Can't be checked)
    5145 #endif
    51465423          return;
    51475424        }
    … …  
    51705447#endif
    51715448#ifdef qdebug
    5172       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (23) qH="<<G4endl;
    5173       else
     5449      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (23) qH=0"<<G4endl;
    51745450#endif
    51755451      delete qH;
    5176 #ifdef qdebug
    5177       qH=0;                                              // @Not necessary@
    5178 #endif
    51795452      if(bB<2) evaHV->push_back(bHadron);         // Fill EvaporatedBaryon (del.equivalent)
    51805453      else if(bB==2) DecayDibaryon(bHadron,evaHV);// => "Dibaryon" case needs decay
    … …  
    52325505        {
    52335506#ifdef qdebug
    5234           if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (24) qH="<<G4endl;
    5235           else
     5507          if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (24) qH=0"<<G4endl;
    52365508#endif
    52375509          delete qH;                            // Chipolino should not be in a sequence
    5238 #ifdef qdebug
    5239           qH=0;                                              // @Not necessary@
    5240 #endif
    52415510          G4LorentzVector fq4M(0.,0.,0.,m1);
    52425511          G4LorentzVector qe4M(0.,0.,0.,m2);
    … …  
    52605529        {
    52615530#ifdef qdebug
    5262           if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (25) qH="<<G4endl;
    5263           else
     5531          if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (25) qH=0"<<G4endl;
    52645532#endif
    52655533          delete qH;
    5266 #ifdef qdebug
    5267           qH=0;                                              // @Not necessary@
    5268 #endif
    52695534          G4cerr<<"**G4QN::EN:M="<<totMass<<"<"<<m1<<"+"<<m2<<",d="<<m1+m2-totMass<<G4endl;
    52705535          throw G4QException("G4QNucleus::EvaporateNucleus: Chipolino is under MassShell");
    … …  
    52805545#endif
    52815546          evaHV->push_back(qH);
     5547        }
     5548        else if ((thePDG==221||thePDG==331)&&totMass>mPi+mPi) // "Decay in pipi" case
     5549        {
    52825550#ifdef qdebug
    5283           qH=0;
    5284 #endif
    5285         }
    5286         else if ((thePDG==221||thePDG==331)&&totMass>mPi+mPi) // "Decay in pipi" case
    5287         {
    5288 #ifdef qdebug
    5289           if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (26) qH="<<G4endl;
    5290           else
     5551          if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (26) qH=0"<<G4endl;
    52915552#endif
    52925553          delete qH;
    5293 #ifdef qdebug
    5294           qH=0;                                              // @Not necessary@
    5295 #endif
    52965554          G4LorentzVector fq4M(0.,0.,0.,mPi);
    52975555          G4LorentzVector qe4M(0.,0.,0.,mPi);
    … …  
    53155573        {
    53165574#ifdef qdebug
    5317           if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (27) qH="<<G4endl;
    5318           else
     5575          if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (27) qH=0"<<G4endl;
    53195576#endif
    53205577          delete qH;
    5321 #ifdef qdebug
    5322           qH=0;                                              // @Not necessary@
    5323 #endif
    53245578          G4LorentzVector fq4M(0.,0.,0.,mPi0);
    53255579          G4LorentzVector qe4M(0.,0.,0.,mPi0);
    … …  
    53435597        {
    53445598#ifdef qdebug
    5345           if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (28) qH="<<G4endl;
    5346           else
     5599          if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (28) qH=0"<<G4endl;
    53475600#endif
    53485601          delete qH;
    5349 #ifdef qdebug
    5350           qH=0;                                              // @Not necessary@
    5351 #endif
    53525602          G4LorentzVector fq4M(0.,0.,0.,0.);
    53535603          G4LorentzVector qe4M(0.,0.,0.,totM);
    … …  
    53715621        {
    53725622#ifdef qdebug
    5373           if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (29) qH="<<G4endl;
    5374           else
     5623          if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (29) qH=0"<<G4endl;
    53755624#endif
    53765625          delete qH;
    5377 #ifdef qdebug
    5378           qH=0;                                              // @Not necessary@
    5379 #endif
    53805626          G4LorentzVector fq4M(0.,0.,0.,0.);
    53815627          G4LorentzVector qe4M(0.,0.,0.,0.);
    … …  
    53995645        {
    54005646#ifdef qdebug
    5401           if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (30) qH="<<G4endl;
    5402           else
     5647          if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (30) qH=0"<<G4endl;
    54035648#endif
    54045649          delete qH;
    5405 #ifdef qdebug
    5406           qH=0;                                              // @Not necessary@
    5407 #endif
    54085650          G4cerr<<"***G4QNucl::EvaNuc: Nuc="<<thePDG<<theQC<<", q4M="<<q4M<<", M="<<totMass
    54095651                <<" < GSM="<<totM<<", 2Pi="<<mPi+mPi<<", 2Pi0="<<mPi0+mPi0<<G4endl;
    … …  
    54155657    {
    54165658#ifdef qdebug
    5417       if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (31) qH="<<G4endl;
    5418       else
     5659      if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (31) qH=0"<<G4endl;
    54195660#endif
    54205661      delete qH;
    5421 #ifdef qdebug
    5422       qH=0;                                              // @Not necessary@
    5423 #endif
    54245662      G4cerr<<"**G4QNuc::EvaNuc:RN="<<thePDG<<theQC<<",q4M="<<q4M<<",qM="<<totMass<<G4endl;
    54255663      throw G4QException("G4QNucleus::EvaporateNucleus: This is not aNucleus nor aHadron");
    … …  
    54315669    G4cout<<"G4QNucleus::EvaporateNucleus: deletedAtEnd, PDG="<<qH->GetPDGCode()<<G4endl;
    54325670    if(!qH) G4cout<<"G4QNucleus::EvaporateNucleus: (20) qH="<<G4endl;
    5433     else
    5434       delete qH;
     5671    else delete qH;
    54355672  }
    54365673#endif
    … …  
    54655702    G4cout<<"--Warning(Upgrade)--G4QNuc::DecIsonuc:FillAsIs,4M="<<q4M<<",QC="<<qQC<<G4endl;
    54665703    evaHV->push_back(qH);                        // fill as it is (delete equivalent)
    5467 #ifdef qdebug
    5468     qH=0;
    5469 #endif
    54705704    return;
    54715705  }
    … …  
    57525986#endif
    57535987    evaHV->push_back(qH);                  // fill as it is (delete equivalent)
    5754 #ifdef qdebug
    5755     qH=0;
    5756 #endif
    57575988    return;
    57585989  }
    … …  
    57645995#endif
    57655996    evaHV->push_back(qH);                  // fill as it is (delete equivalent)
    5766 #ifdef qdebug
    5767     qH=0;
    5768 #endif
    57695997    return;
    57705998  }
    … …  
    57746002#endif
    57756003  delete qH;
    5776 #ifdef qdebug
    5777     qH=0;
    5778 #endif
    57796004  if(qBN)
    57806005  {
    … …  
    58166041//Decay of the excited dibaryon in two baryons
    58176042void G4QNucleus::DecayDibaryon(G4QHadron* qH, G4QHadronVector* evaHV)
    5818 {//  ============================================
     6043{//  ================================================================
    58196044  static const G4double mPi  = G4QPDGCode(211).GetMass();
    58206045  static const G4double mNeut= G4QPDGCode(2112).GetMass();
    … …  
    59776202    {
    59786203      DecayAntiStrange(qH,evaHV);
    5979 #ifdef qdebug
    5980       qH=0;
    5981 #endif
    59826204      return;
    59836205    }
    … …  
    59856207    {
    59866208      delete qH;
    5987 #ifdef qdebug
    5988       qH=0;
    5989 #endif
    59906209      G4cerr<<"***G4QN::DecDiBar: badPDG="<<qPDG<<" or smallM="<<qM<<",2mP="<<dProt
    59916210            <<",2mN="<<dNeut<<G4endl;
    … …  
    60126231      //throw G4QException("***G4QNucleus::DecayDibaryon: DiBaryon DecayIn2 error");
    60136232      evaHV->push_back(qH);
    6014 #ifdef qdebug
    6015       qH=0;
    6016 #endif
    60176233      return;
    60186234    }
    … …  
    60226238#endif
    60236239    delete qH;
    6024 #ifdef qdebug
    6025     qH=0;
    6026 #endif
    60276240    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom); // Create a Hadron for the 1-st baryon
    60286241    evaHV->push_back(H1);                 // Fill "H1" (delete equivalent)
    … …  
    60476260      //throw G4QException("***G4QNucleus::DecDibaryon: Dibaryon DecayIn2 error");
    60486261      evaHV->push_back(qH);
    6049 #ifdef qdebug
    6050       qH=0;
    6051 #endif
    60526262      return;
    60536263    }
    … …  
    60706280      // Should not be here as sum was already compared with qM above for the first delta
    60716281      delete qH;
    6072 #ifdef qdebug
    6073       qH=0;
    6074 #endif
    60756282      G4cerr<<"***G4QNucl::DecDibar:fPDG="<<fPDG<<"(fM="<<fMass<<") + sPDG="<<sPDG<<"(sM="
    60766283            <<sMass<<")="<<sum<<" >? (DD2,Can't be here) TotM="<<q4M.m()<<q4M<<G4endl;
    … …  
    60866293    evaHV->push_back(H4);                 // Fill "H2" (delete equivalent)
    60876294    delete qH;
    6088 #ifdef qdebug
    6089     qH=0;
    6090 #endif
    60916295  }
    60926296  else
    … …  
    61066310      //throw G4QException("G4QNucleus::DecayDibaryon: diBar DecayIn3 error");
    61076311      evaHV->push_back(qH);
    6108 #ifdef qdebug
    6109       qH=0;
    6110 #endif
    61116312      return;
    61126313    }
    … …  
    61196320    // Instead
    61206321    delete qH;
    6121 #ifdef qdebug
    6122     qH=0;
    6123 #endif
    61246322    //
    61256323    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom); // Create a Hadron for the 1-st baryon
    … …  
    61396337#endif
    61406338} // End of DecayDibaryon
     6339
     6340//Decay of the excited anti-dibaryon in two anti-baryons
     6341void G4QNucleus::DecayAntiDibaryon(G4QHadron* qH, G4QHadronVector* evaHV)
     6342{//  ====================================================================
     6343  static const G4double mPi  = G4QPDGCode(211).GetMass();
     6344  static const G4double mNeut= G4QPDGCode(2112).GetMass();
     6345  static const G4double mProt= G4QPDGCode(2212).GetMass();
     6346  static const G4double mSigM= G4QPDGCode(3112).GetMass();
     6347  static const G4double mLamb= G4QPDGCode(3122).GetMass();
     6348  static const G4double mSigP= G4QPDGCode(3222).GetMass();
     6349  static const G4double mKsiM= G4QPDGCode(3312).GetMass();
     6350  static const G4double mKsiZ= G4QPDGCode(3322).GetMass();
     6351  static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);
     6352  static const G4double mPiN = mPi+mNeut;
     6353  static const G4double mPiP = mPi+mProt;
     6354  static const G4double dmPiN= mPiN+mPiN;
     6355  static const G4double dmPiP= mPiP+mPiP;
     6356  static const G4double nnPi = mNeut+mPiN;
     6357  static const G4double ppPi = mProt+mPiP;
     6358  static const G4double lnPi = mLamb+mPiN;
     6359  static const G4double lpPi = mLamb+mPiP;
     6360  static const G4double dNeut= mNeut+mNeut;
     6361  static const G4double dProt= mProt+mProt;
     6362  static const G4double dLamb= mLamb+mLamb;
     6363  static const G4double dLaNe= mLamb+mNeut;
     6364  static const G4double dLaPr= mLamb+mProt;
     6365  static const G4double dSiPr= mSigP+mProt;
     6366  static const G4double dSiNe= mSigM+mNeut;
     6367  static const G4double dKsPr= mKsiZ+mProt;
     6368  static const G4double dKsNe= mKsiM+mNeut;
     6369  static const G4double eps  = 0.003;
     6370  static const G4QNucleus vacuum(90000000);
     6371  G4bool four=false;                           // defFALSE for 4-particle decay of diDelta
     6372  G4LorentzVector q4M = qH->Get4Momentum();    // Get 4-momentum of the Dibaryon
     6373  G4int          qPDG = qH->GetPDGCode();      // PDG Code of the decaying dybaryon
     6374  G4double         qM = q4M.m();               // Mass of the decaying anti-di-baryon
     6375  G4double         rM = qM+eps;                // Just to avoid the computer accuracy
     6376#ifdef debug
     6377  G4cout<<"G4QNucl::DecayAntiDibar:*Called* PDG="<<qPDG<<",4Mom="<<q4M<<", M="<<qM<<G4endl;
     6378#endif
     6379  // Select a chanel of the dibaryon decay (including Delta+Delta-> 4 particle decay
     6380  G4int          fPDG = -2212;                 // Prototype for anti-pp case
     6381  G4int          sPDG = -2212;
     6382  G4int          tPDG = 0;                     // Zero prototype to separate 3 from 2
     6383  G4double       fMass= mProt;
     6384  G4double       sMass= mProt;
     6385  G4double       tMass= mPi;
     6386  if     (qPDG==89996002 && rM>=dmPiP)         // "anti-diDelta++" case
     6387  {
     6388    sPDG = -211;
     6389    sMass= mPi;
     6390    four = true;
     6391  }
     6392  else if(qPDG==90001996 && rM>=dmPiN)         // "diDelta--" case
     6393  {
     6394    sPDG = 211;
     6395    fPDG = -2112;
     6396    sMass= mPi;
     6397    fMass= mNeut;
     6398    four = true;
     6399  }
     6400  else if(qPDG==89999998 && rM>=dNeut)         // "dineutron" case
     6401  {
     6402    fPDG = -2112;
     6403    sPDG = -2112;
     6404    fMass= mNeut;
     6405    sMass= mNeut;   
     6406  }
     6407  else if(qPDG==89998999 && rM>=mDeut)         // "exited deutron" case
     6408  {
     6409    if(fabs(qM-mDeut)<eps)
     6410    {
     6411      evaHV->push_back(qH);                    // Fill as it is (delete equivalent)
     6412      return;
     6413    }
     6414    else if(mProt+mNeut<rM)
     6415    {
     6416      fPDG = -2112;
     6417      fMass= mNeut;   
     6418    }
     6419    else
     6420    {
     6421      fPDG = 22;
     6422      sPDG = 89998999;                         // Anti-deuteron
     6423      fMass= 0.;
     6424      sMass= mDeut;   
     6425      G4cout<<"--Warning--G4QNucl::DecayAntiDibar:ANTI-DEUTERON is created M="<<rM<<G4endl;
     6426    }
     6427  }
     6428  else if(qPDG==88999999 && rM>=dLaNe)         // "Lambda-neutron" case
     6429  {
     6430    fPDG = -2112;
     6431    sPDG = -3122;
     6432    fMass= mNeut;
     6433    sMass= mLamb;   
     6434  }
     6435  else if(qPDG==88999999 && rM>=dLaPr)         // "Lambda-proton" case
     6436  {
     6437    sPDG = -3122;
     6438    sMass= mLamb;   
     6439  }
     6440  else if(qPDG==90000997 && rM>=nnPi)         // "neutron/Delta-" case
     6441  {
     6442    fPDG = -2112;
     6443    sPDG = -2112;
     6444    tPDG = 211;
     6445    fMass= mNeut;
     6446    sMass= mNeut;   
     6447  }
     6448  else if(qPDG==89997001 && rM>=ppPi)         // "proton/Delta++" case
     6449  {
     6450    tPDG = -211;
     6451  }
     6452  else if(qPDG==89000998 && rM>=lnPi)         // "lambda/Delta-" case
     6453  {
     6454    fPDG = -2112;
     6455    sPDG = -3122;
     6456    tPDG = 211;
     6457    fMass= mNeut;
     6458    sMass= mLamb;   
     6459  }
     6460  else if(qPDG==889998001 && rM>=lpPi)         // "lambda/Delta+" case
     6461  {
     6462    sPDG = -3122;
     6463    tPDG = -211;
     6464    sMass= mLamb;   
     6465  }
     6466  else if(qPDG==89000998 && rM>=dSiNe)         // "Sigma-/neutron" case
     6467  {
     6468    fPDG = -2112;
     6469    sPDG = -3112;
     6470    fMass= mNeut;
     6471    sMass= mSigM;   
     6472  }
     6473  else if(qPDG==88998001 && rM>=dSiPr)         // "Sigma+/proton" case
     6474  {
     6475    sPDG = -3222;
     6476    sMass= mSigP;   
     6477  }
     6478  else if(qPDG==88000000 && rM>=dLamb)         // "diLambda" case
     6479  {
     6480    fPDG = -3122;
     6481    sPDG = -3122;
     6482    fMass= mLamb;
     6483    sMass= mLamb;   
     6484  }
     6485  else if(qPDG==88000999 && rM>=dKsNe)         // "Ksi-/neutron" case
     6486  {
     6487    fPDG = -2112;
     6488    sPDG = -3312;
     6489    fMass= mNeut;
     6490    sMass= mKsiM;   
     6491  }
     6492  else if(qPDG==87999001 && rM>=dKsPr)         // "Ksi0/proton" case
     6493  {
     6494    sPDG = -3322;
     6495    sMass= mKsiZ;   
     6496  }
     6497  else if(qPDG!=89998000|| rM<dProt)           // Other possibilities (if not a default)
     6498  {
     6499    G4int qS = qH->GetStrangeness();
     6500    G4int qB = qH->GetBaryonNumber();
     6501    if(qB>0&&qS<0)                             // Antistrange diBarion
     6502    {
     6503      DecayAntiStrange(qH,evaHV);
     6504      return;
     6505    }
     6506    else
     6507    {
     6508      delete qH;
     6509      G4cerr<<"**G4QNuc::DecayAntiDiBar: badPDG="<<qPDG<<" or smallM="<<qM<<", 2mP="<<dProt
     6510            <<", 2mN="<<dNeut<<G4endl;
     6511      // @@ Nothing to do. Just 2 GeV disappears... Very rare! Just to avoid the exception.
     6512      //throw G4QException("G4QNucleus::DecayDibar: Unknown PDG code or small Mass of DB");
     6513    }
     6514  }
     6515  G4LorentzVector f4Mom(0.,0.,0.,fMass);
     6516  G4LorentzVector s4Mom(0.,0.,0.,sMass);
     6517  G4LorentzVector t4Mom(0.,0.,0.,tMass);
     6518  if(!tPDG&&!four)
     6519  {
     6520    G4double sum=fMass+sMass;
     6521    if(fabs(qM-sum)<eps)
     6522    {
     6523      f4Mom=q4M*(fMass/sum);
     6524      s4Mom=q4M*(sMass/sum);
     6525    }
     6526    else if(qM<sum || !G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
     6527    {
     6528      G4cout<<"---Warning---G4QN::DecAntiDib:fPDG="<<fPDG<<"(M="<<fMass<<")+sPDG="<<sPDG
     6529            <<"(M="<<sMass<<")="<<sum<<" >? TotM="<<q4M.m()<<q4M<<G4endl;
     6530      //G4cerr<<"***G4QNucl::DecayDiBar: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
     6531      //throw G4QException("***G4QNucleus::DecayDibaryon: DiBaryon DecayIn2 error");
     6532      evaHV->push_back(qH);
     6533      return;
     6534    }
     6535#ifdef debug
     6536    G4cout<<"G4QNucleus::DecayAntiDibaryon:(2) *DONE* f4M="<<f4Mom<<",fPDG="<<fPDG
     6537          <<", s4M="<<s4Mom<<",sPDG="<<sPDG<<G4endl;
     6538#endif
     6539    delete qH;
     6540    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom); // Create a Hadron for the 1-st baryon
     6541    evaHV->push_back(H1);                 // Fill "H1" (delete equivalent)
     6542    G4QHadron* H2 = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the 2-nd baryon
     6543    evaHV->push_back(H2);                 // Fill "H2" (delete equivalent)
     6544  }
     6545  else if(four)
     6546  {
     6547    q4M=q4M/2.;                                // Divided in 2 !!!
     6548    qM/=2.;                                    // Divide the mass in 2 !
     6549    G4double sum=fMass+sMass;
     6550    if(fabs(qM-sum)<eps)
     6551    {
     6552      f4Mom=q4M*(fMass/sum);
     6553      s4Mom=q4M*(sMass/sum);
     6554    }
     6555    else if(qM<sum || !G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
     6556    {
     6557      G4cout<<"---Warning---G4QN::DecAntiDib:fPDG="<<fPDG<<"(M="<<fMass<<")+sPDG="<<sPDG
     6558            <<"(M="<<sMass<<")"<<"="<<sum<<">tM="<<q4M.m()<<q4M<<G4endl;
     6559      //G4cerr<<"***G4QN::DecayDibaryon: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
     6560      //throw G4QException("***G4QNucleus::DecDibaryon: Dibaryon DecayIn2 error");
     6561      evaHV->push_back(qH);
     6562      return;
     6563    }
     6564#ifdef debug
     6565    G4cout<<"G4QNucleus::DecayAntiDibaryon:(3) *DONE* f4M="<<f4Mom<<",fPDG="<<fPDG
     6566          <<", s4M="<<s4Mom<<",sPDG="<<sPDG<<G4endl;
     6567#endif
     6568    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom); // Create a Hadron for the 1-st baryon
     6569    evaHV->push_back(H1);                      // Fill "H1" (delete equivalent)
     6570    G4QHadron* H2 = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the 2-nd baryon
     6571    evaHV->push_back(H2);                      // Fill "H2" (delete equivalent)
     6572    // Now the second pair mus be decayed
     6573    if(fabs(qM-sum)<eps)
     6574    {
     6575      f4Mom=q4M*(fMass/sum);
     6576      s4Mom=q4M*(sMass/sum);
     6577    }
     6578    else if(!G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
     6579    {
     6580      // Should not be here as sum was already compared with qM above for the first delta
     6581      delete qH;
     6582      G4cerr<<"**G4QNucl::DecAntiDibar:fPDG="<<fPDG<<"(fM="<<fMass<<")+sPDG="<<sPDG<<"(sM="
     6583            <<sMass<<")="<<sum<<" >? (DD2,Can't be here) TotM="<<q4M.m()<<q4M<<G4endl;
     6584      throw G4QException("G4QNucleus::DecayAntiDibaryon: General DecayIn2 error");
     6585    }
     6586#ifdef debug
     6587    G4cout<<"G4QNucl::DecayAntiDibaryon:(4) *DONE* f4M="<<f4Mom<<",fPDG="<<fPDG
     6588          <<", s4M="<<s4Mom<<",sPDG="<<sPDG<<G4endl;
     6589#endif
     6590    G4QHadron* H3 = new G4QHadron(fPDG,f4Mom); // Create a Hadron for the 1-st baryon
     6591    evaHV->push_back(H3);                      // Fill "H1" (delete equivalent)
     6592    G4QHadron* H4 = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the 2-nd baryon
     6593    evaHV->push_back(H4);                      // Fill "H2" (delete equivalent)
     6594    delete qH;
     6595  }
     6596  else
     6597  {
     6598    G4double sum=fMass+sMass+tMass;
     6599    if(fabs(qM-sum)<eps)
     6600    {
     6601      f4Mom=q4M*(fMass/sum);
     6602      s4Mom=q4M*(sMass/sum);
     6603      t4Mom=q4M*(tMass/sum);
     6604    }
     6605    else if(qM<sum || !G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom))
     6606    {
     6607      G4cout<<"-Warning-G4QN::DecAntiDib:fPDG="<<fPDG<<"(M="<<fMass<<")+sPDG="<<sPDG<<"(M="
     6608            <<sMass<<")+tPDG="<<tPDG<<"(tM="<<tMass<<")="<<sum<<">TotM="<<q4M.m()<<G4endl;
     6609      //G4cerr<<"***G4QNuc::DecayDibaryon:qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
     6610      //throw G4QException("G4QNucleus::DecayDibaryon: diBar DecayIn3 error");
     6611      evaHV->push_back(qH);
     6612      return;
     6613    }
     6614#ifdef debug
     6615    G4cout<<"G4QNuc::DecayAbtiDibaryon:(5) *DONE* f4M="<<f4Mom<<",fPDG="<<fPDG<<", s4M="
     6616          <<s4Mom<<",sPDG="<<sPDG<<", t4M="<<t4Mom<<",tPDG="<<tPDG<<G4endl;
     6617#endif
     6618    //qH->SetNFragments(2);                    // Fill a#of fragments to decaying Dibaryon
     6619    //evaHV->push_back(qH);               // Fill hadron with nf=2 (delete equivalent)
     6620    // Instead
     6621    delete qH;
     6622    //
     6623    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom); // Create a Hadron for the 1-st baryon
     6624    evaHV->push_back(H1);                 // Fill "H1" (delete equivalent)
     6625    G4QHadron* H2 = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the 2-nd baryon
     6626    evaHV->push_back(H2);                 // Fill "H2" (delete equivalent)
     6627    G4QHadron* H3 = new G4QHadron(tPDG,t4Mom); // Create a Hadron for the meson
     6628    evaHV->push_back(H3);                 // Fill "H3" (delete equivalent)
     6629  }
     6630#ifdef qdebug
     6631  if (qH)
     6632  {
     6633    G4cout<<"G4QNucleus::DecayDiBaryon: deleted at end - PDG="<<qH->GetPDGCode()<<G4endl;
     6634    delete qH;
     6635  }
     6636#endif
     6637} // End of DecayAntiDibaryon
    61416638
    61426639//Decay of the nuclear states with antistrangeness (K:/K0)
    … …  
    61596656  {
    61606657    delete qH;
    6161 #ifdef qdebug
    6162     qH=0;
    6163 #endif
    61646658    G4cerr<<"G4QNuc::DecayAntiStrange:QC="<<qQC<<",S="<<qS<<",B="<<qB<<",4M="<<q4M<<G4endl;
    61656659    throw G4QException("G4QNucleus::DecayAntiStrange: not an Anti Strange Nucleus");
    … …  
    63626856#endif
    63636857    delete qH;
    6364 #ifdef qdebug
    6365     qH=0;
    6366 #endif
    63676858    //
    63686859    f4Mom/=n1;
    … …  
    64056896#endif
    64066897    delete qH;
    6407 #ifdef qdebug
    6408     qH=0;
    6409 #endif
    64106898    //
    64116899    f4Mom/=n1;
    … …  
    64716959  {
    64726960    delete qH;
    6473 #ifdef qdebug
    6474     qH=0;
    6475 #endif
    64766961    G4cerr<<"***G4QNuc::DecayMultyBaryon: PDG="<<qPDG<<G4endl;
    64776962    throw G4QException("***G4QNuc::DecayMultyBaryon: Can not decay this PDG Code");
    … …  
    64816966  {
    64826967    delete qH;
    6483 #ifdef qdebug
    6484     qH=0;
    6485 #endif
    64866968    G4cerr<<"**G4QNucleus::DecayMultyBaryon: PDG="<<qPDG<<G4endl;
    64876969    throw G4QException("***G4QNuc::DecayMultyBaryon: Unknown PDG code of the MultiBaryon");
    … …  
    65126994#endif
    65136995    delete qH;
    6514 #ifdef qdebug
    6515     qH=0;
    6516 #endif
    65176996    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom);   // Create a Hadron for the 1-st baryon
    65186997    evaHV->push_back(H1);                   // Fill "H1" (delete equivalent)
    … …  
    65467025#endif
    65477026    delete qH;
    6548 #ifdef qdebug
    6549     qH=0;
    6550 #endif
    65517027    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom);   // Create a Hadron for the 1-st baryon
    65527028    evaHV->push_back(H1);                   // Fill "H1" (delete equivalent)
    … …  
    65677043#endif
    65687044    delete qH;
    6569 #ifdef qdebug
    6570     qH=0;
    6571 #endif
    65727045    for(G4int h=0; h<totBN; h++)
    65737046    {
    … …  
    66197092    {
    66207093      delete qH;
    6621 #ifdef qdebug
    6622       qH=0;
    6623 #endif
    66247094      G4cerr<<"***G4QNu::DecAlDiN:M(He6="<<mHel6<<")="<<qM<<"<"<<mNeut+mNeut+mAlph<<G4endl;
    66257095      throw G4QException("G4QNuc::DecayAlphaDiN: Cannot decay excited He6 with this mass");
    … …  
    66297099  {
    66307100    delete qH;
    6631 #ifdef qdebug
    6632     qH=0;
    6633 #endif
    66347101    G4cerr<<"***G4QNuc::DecayAlphaDiN: PDG="<<qPDG<<G4endl;
    66357102    throw G4QException("G4QNuc::DecayAlphaDiN: Can not decay this PDG Code");
    … …  
    66587125#endif
    66597126  delete qH;
    6660 #ifdef qdebug
    6661   qH=0;
    6662 #endif
    66637127  G4QHadron* H1 = new G4QHadron(fPDG,f4Mom);    // Create a Hadron for the 1-st baryon
    66647128  evaHV->push_back(H1);                    // Fill "H1" (delete equivalent)
    … …  
    67367200      //throw G4QException("G4QNucleus::DecayAlphaBar: DecayIn2 didn't succeed for 3/2");
    67377201      evaHV->push_back(qH);
    6738 #ifdef qdebug
    6739       qH=0;
    6740 #endif
    67417202      return;
    67427203    }
    … …  
    67457206#endif
    67467207    delete qH;
    6747 #ifdef qdebug
    6748     qH=0;
    6749 #endif
    67507208    G4LorentzVector rf4Mom=f4Mom/2;
    67517209    G4QHadron* H1 = new G4QHadron(fPDG,rf4Mom); // Create a Hadron for the 1-st baryon
    … …  
    67887246      //throw G4QException("G4QNucleus::DecayAlphaBar: t/nn,He3/pp DecayIn3 didn't");
    67897247      evaHV->push_back(qH);
    6790 #ifdef qdebug
    6791       qH=0;
    6792 #endif
    67937248      return;
    67947249    }
    … …  
    67977252#endif
    67987253    delete qH;
    6799 #ifdef qdebug
    6800     qH=0;
    6801 #endif
    68027254    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom);   // Create a Hadron for the 1-st baryon
    68037255    evaHV->push_back(H1);                   // Fill "H1" (delete equivalent)
    … …  
    68477299      //throw G4QException("G4QNucl::DecayAlphaBar:QuintBaryon DecayIn2 didn't succeed");
    68487300      evaHV->push_back(qH);
    6849 #ifdef qdebug
    6850       qH=0;
    6851 #endif
    68527301      return;
    68537302    }
    … …  
    68567305#endif
    68577306    delete qH;
    6858 #ifdef qdebug
    6859     qH=0;
    6860 #endif
    68617307    G4LorentzVector rf4Mom=f4Mom/4;
    68627308    G4QHadron* H1 = new G4QHadron(fPDG,rf4Mom); // Create a Hadron for the 1-st baryon
    … …  
    68877333    {
    68887334      evaHV->push_back(qH);                     // Fill hadron as it is (delete equivalent)
    6889 #ifdef qdebug
    6890       qH=0;
    6891 #endif
    68927335      //EvaporateNucleus(qH, evaHV);            // Evaporate Nucleus (delete equivivalent)
    68937336      return;
    … …  
    69197362      //throw G4QException("***G4QNucl::DecayAlphaBar:Alpha+Baryon DecIn2 didn't succeed");
    69207363      evaHV->push_back(qH);
    6921 #ifdef qdebug
    6922       qH=0;
    6923 #endif
    69247364      return;
    69257365    }
    … …  
    69287368#endif
    69297369    delete qH;
    6930 #ifdef qdebug
    6931     qH=0;
    6932 #endif
    69337370    G4QHadron* H1 = new G4QHadron(fPDG,f4Mom);      // Create a Hadron for the alpha
    69347371    evaHV->push_back(H1);                      // Fill "H1" (delete equivalent)
    … …  
    69577394  {
    69587395    delete qH;
    6959 #ifdef qdebug
    6960     qH=0;
    6961 #endif
    69627396    G4cerr<<"***G4QNucleus::DecayAlphaAlpha: qPDG="<<qPDG<<G4endl;
    69637397    throw G4QException("***G4QNucleus::DecayAlphaAlpha: Not Be8 state decais in 2 alphas");
    … …  
    70247458#endif
    70257459  delete qH;
    7026 #ifdef qdebug
    7027   qH=0;
    7028 #endif
    70297460  G4QHadron* H1 = new G4QHadron(fPDG,f4Mom);      // Create a Hadron for the 1-st alpha
    70307461  evaHV->push_back(H1);                      // Fill "H1" (delete equivalent)

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QPDGCode.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4QPDGCode.cc,v 1.63 2009/11/03 16:13:37 mkossov Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4QPDGCode.cc,v 1.65 2010/06/10 08:37:27 mkossov Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//      ---------------- G4QPDGCode ----------------
    … …  
    6767  }
    6868#ifdef debug
    69   G4cout<<"G4QPDGCode:Constructer(PDG) PDG="<<PDGCode<<", QCode="<<theQCode<<G4endl; 
     69  if(PDGCode==3222)G4cout<<"G4QPDGCd:Con(PDG) PDG="<<PDGCode<<", QCode="<<theQCode<<G4endl;
    7070#endif
    7171}
    … …  
    327327    return -2;
    328328  }
    329   else if (PDGC>80000000&&PDGC<100000000) // Try to convert the NUCCoding to PDGCoding
     329  else if (PDGC>80000000 && PDGC<100000000) // Try to convert the NUCCoding to PDGCoding
    330330  {
    331331    if(PDGC==90000000) return 6;
    … …  
    376376      }
    377377    } // End of meson case
    378     if(b>0)                                        // --> Baryon case
    379     {
    380       if(b==1)
     378    if(b>0)                                        // --> Baryoniums case
     379    {
     380      if(b==1)                                     // --> Baryons+Hyperons
    381381      {
    382382        if(!s)                                     // --> Baryons
    383383        {
    384384          if(z==-1)    return 34;                  // Delta-
    385           else if(!z)  return 91;                  // neutron
    386           else if(z==1)return 91;                  // proton
     385          else if(!z)  return 20;                  // neutron
     386          else if(z==1)return 21;                  // proton
    387387          else if(z==2)return 37;                  // Delta++
    388388          else if(z==3||z==-2)return -1;           // Delta+pi Chipolino
    … …  
    391391        else if(s==1)                              // --> Hyperons
    392392        {
    393           if(z==-1)    return 93;                  // Sigma-
    394           else if(!z)  return 92;                  // Lambda (@@ 24->Sigma0)
    395           else if(z==1)return 94;                  // Sigma+
     393          if(z==-1)    return 23;                  // Sigma-
     394          else if(!z)  return 22;                  // Lambda (@@ 24->Sigma0)
     395          else if(z==1)return 25;                  // Sigma+
    396396          else if(z==2||z==-2) return -1;          // Sigma+pi Chipolino
    397397          else         return -2;                  // Not supported by Q Code
    … …  
    399399        else if(s==2)                              // --> Xi Hyperons
    400400        {
    401           if(z==-1)    return 95;                  // Xi-
    402           else if(!z)  return 96;                  // Xi0
     401          if(z==-1)    return 26;                  // Xi-
     402          else if(!z)  return 27;                  // Xi0
    403403          else if(z==1||z==-2)return -1;           // Xi+pi Chipolino
    404404          else         return -2;                  // Not supported by Q Code
    … …  
    406406        else if(s==3)                              // --> Xi Hyperons
    407407        {
    408           if(z==-1)    return 97;                  // Omega-
     408          if(z==-1)    return 44;                  // Omega-
    409409          else if(!z||z==-2)  return -1;           // Omega+pi Chipolino
    410410          else         return -2;                  // Not supported by Q Code
    … …  
    433433    }
    434434  }
    435   if (PDGC<80000000)              // ----> Direct Baryons & Mesons
    436   {
    437     if     (PDGC<100)             // => Leptons and field bosons
     435  if (PDGC<80000000)                // ----> Direct Baryons & Mesons
     436  {
     437    if     (PDGC<100)               // => Leptons and field bosons
    438438    {
    439439      if     (PDGC==10)  return -1; // Chipolino
    … …  
    458458      else if(PDGC==220) return 12; // Midle Regeon-Pomeron
    459459      else if(PDGC==330) return 13; // High Regeon-Pomeron
    460 #ifdef pdebug
     460#ifdef debug
    461461      G4cout<<"***G4QPDGCode::MakeQCode: (0) Unknown in Q-System code: "<<PDGCode<<G4endl;
    462462#endif
    … …  
    465465    else Q=qr[r];
    466466    G4int p=PDGC/10;                // Quark Content
    467     if(r%2)                         // (2s+1 is odd) Mesons are all the same
     467    if(r%2)                         // (2s+1 is odd) Mesons
    468468    {
    469469      if     (p==11) return Q+=1;
    … …  
    475475      else
    476476      {
    477 #ifdef pdebug
    478         G4cout<<"***G4QPDGCode::MakeQCode:(1) Unknown in Q-System code: "<<PDGCode<<G4endl;
     477#ifdef debug
     478        G4cout<<"*Warning*G4QPDGCode::MakeQCode:(1)UnknownQCode for PDG="<<PDGCode<<G4endl;
    479479#endif
    480480        return -2;
    … …  
    496496        else
    497497        {
    498 #ifdef pdebug
    499           G4cout<<"**G4QPDGCode::MakeQCode:(2) Unknown in Q-System code:"<<PDGCode<<G4endl;
     498#ifdef debug
     499          G4cout<<"*Warning*G4QPDGCode::MakeQCode:(2) UnknownQCode, PDG="<<PDGCode<<G4endl;
    500500#endif
    501501          return -2;
    … …  
    517517        else
    518518        {
    519 #ifdef pdebug
     519#ifdef debug
    520520          G4cout<<"**G4QPDGCode::MakeQCode:(3) Unknown in Q-System code:"<<PDGCode<<G4endl;
    521521#endif
    … …  
    532532    if(t%3 || abs(t)<3)       // b=0 are in mesons
    533533    {
    534 #ifdef pdebug
     534#ifdef debug
    535535      G4cout<<"***G4QPDGCode::MakeQCode: Unknown PDGCode="<<PDGCode<<", t="<<t<<G4endl;
    536536#endif
    … …  
    547547        else
    548548        {
    549 #ifdef pdebug
     549#ifdef debug
    550550          G4cout<<"**G4QPDGCode::MakeQCode:(5) Unknown in Q-System code:"<<PDGCode<<G4endl;
    551551#endif
    … …  
    563563        else
    564564        {
    565 #ifdef pdebug
     565#ifdef debug
    566566          G4cout<<"**G4QPDGCode::MakeQCode:(6) Unknown in Q-System code:"<<PDGCode<<G4endl;
    567567#endif
    … …  
    581581        else
    582582        {
    583 #ifdef pdebug
     583#ifdef debug
    584584          G4cout<<"**G4QPDGCode::MakeQCode:(7) Unknown in Q-System code:"<<PDGCode<<G4endl;
    585585#endif
    … …  
    632632    }
    633633  }
    634 #ifdef pdebug
    635   G4cout<<"***G4QPDGCode::MakeQCode: () Unknown in Q-System code: "<<PDGCode<<G4endl;
    636 #endif
    637 // return -2; not reachable statement 
     634  G4cout<<"*Warning*G4QPDGCode::MakeQCode:() Unknown Q Code for PDG = "<<PDGCode<<G4endl;
     635  return -2; //not reachable statement (fake, only for compiler) 
    638636}
    639637
    … …  
    642640{//      =====================
    643641  G4int ab=theQCode;
    644 #ifdef debug
     642#ifdef pdebug
     643  G4bool pPrint = thePDGCode == 3222 || ab == 25;
     644  if(pPrint)
    645645  G4cout<<"G4QPDGCode::GetMass: Mass for Q="<<ab<<",PDG="<<thePDGCode<<",N="<<nQHM<<G4endl;
    646646#endif
    647647  if ( (ab < 0 && thePDGCode < 80000000) || !thePDGCode) {
    648 #ifdef debug
    649     if(thePDGCode!=10)
     648#ifdef pdebug
     649    if(thePDGCode!=10 && pPrint)
    650650      G4cout<<"**G4QPDGCode::GetMass:m=100000.,QC="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
    651651#endif
    … …  
    654654  else if(ab>-1 && ab<nQHM)
    655655  {
    656 #ifdef debug
    657     G4cout<<"G4QPDGCode::GetMa:m="<<QHaM(ab)<<",Q="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
    658 #endif
    659     return QHaM(ab);            // Get mass from the table
     656    G4double mass = QHaM(ab);
     657#ifdef pdebug
     658    //if(pPrint)
     659    if(thePDGCode == 3222 || ab == 25)
     660    G4cout<<"G4QPDGCode::GetMa:m="<<mass<<",Q="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
     661#endif
     662    return mass;                                // Get mass from the table
    660663  }
    661664  //if(szn==0) return m[15];                    // @@ mPi0   @@ MK ?
    662665  if(thePDGCode==90000000)
    663666  {
    664 #ifdef debug
     667#ifdef pdebug
     668    if(pPrint)
    665669    G4cout<<"G4QPDGCode::GetMass:***m=0, QC="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
    666670#endif
    … …  
    672676  ConvertPDGToZNS(thePDGCode, z, n, s);
    673677  G4double m=GetNuclMass(z,n,s);
    674 #ifdef debug
     678#ifdef pdebug
     679  if(pPrint)
    675680  G4cout<<"G4QPDG::GetM:PDG="<<thePDGCode<<"=>Z="<<z<<",N="<<n<<",S="<<s<<",M="<<m<<G4endl;
    676681#endif
    … …  
    684689  //static const int nW = 72;
    685690  //static const int nW = 80; // "Isobars"
    686   static const int nW = 90; // "Leptons/Hypernuclei"
    687   static G4double width[nW] = {0.,0.,0.,0.,0.,0.,0.,2.495,2.118,10.
     691  static G4double width[nQHM] = {0.,0.,0.,0.,0.,0.,0.,2.495,2.118,10.
    688692    ,  10., 800.,  75., 350.,   0.,   0., .00118,  0.,   0., .203
    689693    ,   0.,   0.,   0.,   0.,   0.,   0.,   0.,    0., 160., 160.
    … …  
    695699    , 198., 149., 120., 120., 170., 170., 120.,  120., 170., 170.};
    696700  G4int ab=abs(theQCode);
    697   if(ab<nW) return width[ab];
     701  if(ab<nQHM) return width[ab];
    698702  return 0.;             // @@ May be real width should be implemented for nuclei e.g. pp
    699703}
    … …  
    928932        {
    929933          NZS1max=s;
    930           G4cout<<">>>>>>>>>>>>>G4QPDGCode::GetMass: Z=-1, S="<<s<<">2 with N="<<n<<G4endl;
     934          G4cout<<">>>>>>>>>>G4QPDGCode::GetNucMass: Z=-1, S="<<s<<">2 with N="<<n<<G4endl;
    931935        }
    932936#endif
    … …  
    988992      {
    989993        NZmin=z;
    990         G4cout<<">>>>>>>>>G4QPDGCode::GetMass: Z="<<z<<"<-1 with N="<<n<<", S="<<s<<G4endl;
     994        G4cout<<">>>>>>G4QPDGCode::GetNucMass: Z="<<z<<"<-1 with N="<<n<<", S="<<s<<G4endl;
    991995      }
    992996#endif
    … …  
    10361040        {
    10371041          NNS1max=s;
    1038           G4cout<<">>>>>>>>>>>>>G4QPDGCode::GetMass: N=-1, S="<<s<<">2 with Z="<<z<<G4endl;
     1042          G4cout<<">>>>>>>>>>G4QPDGCode::GetNucMass: N=-1, S="<<s<<">2 with Z="<<z<<G4endl;
    10391043        }
    10401044#endif
    … …  
    10961100      {
    10971101        NNmin=n;
    1098         G4cout<<">>>>>>>>>G4QPDGCode::GetMass: N="<<n<<"<-1 with Z="<<z<<", S="<<s<<G4endl;
     1102        G4cout<<">>>>>>G4QPDGCode::GetNucMass: N="<<n<<"<-1 with Z="<<z<<", S="<<s<<G4endl;
    10991103      }
    11001104#endif
    … …  
    11031107  }
    11041108  //return CalculateNuclMass(z,n,s); // @@ This is just to compare the calculation time @@
     1109  if(nz >= 256 || z >= nEl) return 256000.;
    11051110  if(!s)                   // **************> S=0 nucleus
    11061111  {
    1107     if(nz==256) return 256000.;
    11081112    G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
    11091113    if(!iNin0[z])          // ====> This element is already initialized
    11101114    {
    11111115#ifdef idebug
    1112       G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=0 is initialized. F="<<iNin0[z]<<G4endl;
     1116      G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=0 is initialized. F="<<iNin0[z]<<G4endl;
    11131117#endif
    11141118      G4int iNfin=iNL[z];
    … …  
    11231127      if(n<iNmin[z])
    11241128      {
    1125         G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=0 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1129        G4cout<<">>>>G4QPDGCode::GetNucM:Z="<<z<<", S=0 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    11261130        iNmin[z]=n;
    11271131      }
    … …  
    11351139      if(dNn>iNmax)
    11361140      {
    1137         G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1141        G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNmax<<G4endl;
    11381142        iNmax=dNn;
    11391143      }
    11401144      if(dNn>iNran[z])
    11411145      {
    1142         G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1146        G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    11431147        iNran[z]=dNn;
    11441148      }
    … …  
    11491153  else if(s==1)            // ******************> S=1 nucleus
    11501154  {
    1151 
    11521155    G4int Nmin=iNF[z];     // Minimun N(Z) for the Dynamic Associative Memory (DAM)
     1156#ifdef pdebug
     1157    G4bool pPrint = !z && !n;
     1158    if(pPrint)
     1159      G4cout<<"G4QPDGC::GetNucM:Nmin="<<Nmin<<",iNin1="<<iNin1[0]<<",iNL="<<iNL[0]<<G4endl;
     1160#endif
    11531161    if(!iNin1[z])          // ====> This element is already initialized
    11541162    {
    1155 #ifdef idebug
    1156       G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=1 is initialized. F="<<iNin1[z]<<G4endl;
     1163#ifdef pdebug
     1164      if(pPrint)
     1165      G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=1 is initialized. F="<<iNin1[z]<<G4endl;
    11571166#endif
    11581167      G4int iNfin=iNL[z];
    … …  
    11671176      if(n<iNmin[z])
    11681177      {
    1169         G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=1 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1178        G4cout<<">>>>G4QPDGCode::GetNucM:Z="<<z<<", S=1 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    11701179        iNmin[z]=n;
    11711180      }
    … …  
    11791188      if(dNn>iNmax)
    11801189      {
    1181         G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1190        G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNmax<<G4endl;
    11821191        iNmax=dNn;
    11831192      }
    11841193      if(dNn>iNran[z])
    11851194      {
    1186         G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1195        G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    11871196        iNran[z]=dNn;
    11881197      }
    … …  
    11971206    {
    11981207#ifdef idebug
    1199       G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-1 is initialized. F="<<iNin9[z]<<G4endl;
     1208      G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=-1 is initialized. F="<<iNin9[z]<<G4endl;
    12001209#endif
    12011210      G4int iNfin=iNL[z];
    … …  
    12101219      if(n<iNmin[z])
    12111220      {
    1212         G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<" ,S=-1 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1221        G4cout<<">>>G4QPDGCode::GetNucM:Z="<<z<<" ,S=-1 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    12131222        iNmin[z]=n;
    12141223      }
    … …  
    12221231      if(dNn>iNmax)
    12231232      {
    1224         G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1233        G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNmax<<G4endl;
    12251234        iNmax=dNn;
    12261235      }
    12271236      if(dNn>iNran[z])
    12281237      {
    1229         G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1238        G4cout<<"G4QPDGCode::GetNucM:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    12301239        iNran[z]=dNn;
    12311240      }
    … …  
    12401249    {
    12411250#ifdef idebug
    1242       G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=2 is initialized. F="<<iNin2[z]<<G4endl;
     1251      G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=2 is initialized. F="<<iNin2[z]<<G4endl;
    12431252#endif
    12441253      G4int iNfin=iNL[z];
    … …  
    12531262      if(n<iNmin[z])
    12541263      {
    1255         G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=2 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1264        G4cout<<">>>>G4QPDGCode::GetNucM:Z="<<z<<", S=2 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    12561265        iNmin[z]=n;
    12571266      }
    … …  
    12651274      if(dNn>iNmax)
    12661275      {
    1267         G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1276        G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNmax<<G4endl;
    12681277        iNmax=dNn;
    12691278      }
    12701279      if(dNn>iNran[z])
    12711280      {
    1272         G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1281        G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    12731282        iNran[z]=dNn;
    12741283      }
    … …  
    12831292    {
    12841293#ifdef idebug
    1285       G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-2 is initialized. F="<<iNin8[z]<<G4endl;
     1294      G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=-2 is initialized. F="<<iNin8[z]<<G4endl;
    12861295#endif
    12871296      G4int iNfin=iNL[z];
    … …  
    12961305      if(n<iNmin[z])
    12971306      {
    1298         G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-2 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1307        G4cout<<">>>G4QPDGCode::GetNucM:Z="<<z<<", S=-2 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    12991308        iNmin[z]=n;
    13001309      }
    … …  
    13081317      if(dNn>iNmax)
    13091318      {
    1310         G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1319        G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNmax<<G4endl;
    13111320        iNmax=dNn;
    13121321      }
    13131322      if(dNn>iNran[z])
    13141323      {
    1315         G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1324        G4cout<<"G4QPDGCode::GetNucM:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    13161325        iNran[z]=dNn;
    13171326      }
    … …  
    13261335    {
    13271336#ifdef idebug
    1328       G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-3 is initialized. F="<<iNin7[z]<<G4endl;
     1337      G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=-3 is initialized. F="<<iNin7[z]<<G4endl;
    13291338#endif
    13301339      G4int iNfin=iNL[z];
    … …  
    13391348      if(n<iNmin[z])
    13401349      {
    1341         G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-3 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1350        G4cout<<">>>G4QPDGCode::GetNucM:Z="<<z<<", S=-3 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    13421351        iNmin[z]=n;
    13431352      }
    … …  
    13511360      if(dNn>iNmax)
    13521361      {
    1353         G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1362        G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNmax<<G4endl;
    13541363        iNmax=dNn;
    13551364      }
    13561365      if(dNn>iNran[z])
    13571366      {
    1358         G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1367        G4cout<<"G4QPDGCode::GetNucM:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    13591368        iNran[z]=dNn;
    13601369      }
    … …  
    13691378    {
    13701379#ifdef idebug
    1371       G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=3 is initialized. F="<<iNin3[z]<<G4endl;
     1380      G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=3 is initialized. F="<<iNin3[z]<<G4endl;
    13721381#endif
    13731382      G4int iNfin=iNL[z];
    … …  
    13821391      if(n<iNmin[z])
    13831392      {
    1384         G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=3 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1393        G4cout<<">>>>G4QPDGCode::GetNucM:Z="<<z<<", S=3 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    13851394        iNmin[z]=n;
    13861395      }
    … …  
    13941403      if(dNn>iNmax)
    13951404      {
    1396         G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1405        G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNmax<<G4endl;
    13971406        iNmax=dNn;
    13981407      }
    13991408      if(dNn>iNran[z])
    14001409      {
    1401         G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1410        G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    14021411        iNran[z]=dNn;
    14031412      }
    … …  
    14121421    {
    14131422#ifdef idebug
    1414       G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-4 is initialized. F="<<iNin6[z]<<G4endl;
     1423      G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=-4 is initialized. F="<<iNin6[z]<<G4endl;
    14151424#endif
    14161425      G4int iNfin=iNL[z];
    … …  
    14251434      if(n<iNmin[z])
    14261435      {
    1427         G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-4 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1436        G4cout<<">>>G4QPDGCode::GetNucM:Z="<<z<<", S=-4 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    14281437        iNmin[z]=n;
    14291438      }
    … …  
    14371446      if(dNn>iNmax)
    14381447      {
    1439         G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1448        G4cout<<"**>G4QPDGCode::GetNucM:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNmax<<G4endl;
    14401449        iNmax=dNn;
    14411450      }
    14421451      if(dNn>iNran[z])
    14431452      {
    1444         G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1453        G4cout<<"G4QPDGCode::GetNucM:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    14451454        iNran[z]=dNn;
    14461455      }
    … …  
    14551464    {
    14561465#ifdef idebug
    1457       G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=4 is initialized. F="<<iNin4[z]<<G4endl;
     1466      G4cout<<"*>G4QPDGCode::GetNucM:Z="<<z<<", S=4 is initialized. F="<<iNin4[z]<<G4endl;
    14581467#endif
    14591468      G4int iNfin=iNL[z];
    … …  
    14681477      if(n<iNmin[z])
    14691478      {
    1470         G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=4 with N="<<n<<"<"<<iNmin[z]<<G4endl;
     1479        G4cout<<">>>>G4QPDGCode::GetNucM:Z="<<z<<", S=4 with N="<<n<<"<"<<iNmin[z]<<G4endl;
    14711480        iNmin[z]=n;
    14721481      }
    … …  
    14801489      if(dNn>iNmax)
    14811490      {
    1482         G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNmax<<G4endl;
     1491        G4cout<<"**>>G4QPDGCode::GetNucM:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNmax<<G4endl;
    14831492        iNmax=dNn;
    14841493      }
    14851494      if(dNn>iNran[z])
    14861495      {
    1487         G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
     1496        G4cout<<">G4QPDGCode::GetNucM:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
    14881497        iNran[z]=dNn;
    14891498      }
    … …  
    14991508      if(s<Smin) Smin=s;
    15001509      if(s>Smax) Smax=s;
    1501       G4cout<<">>G4QPDGCode::GetMass:Z="<<z<<" with S="<<s<<",N="<<n<<" (Improve)"<<G4endl;
     1510      G4cout<<">>G4QPDGCode::GetNucM:Z="<<z<<" with S="<<s<<",N="<<n<<" (Improve)"<<G4endl;
    15021511    }
    15031512#endif
    … …  
    15081517#endif
    15091518  return rm;
    1510 }
     1519} // End of GetNuclearMass
    15111520
    15121521// Calculate a nuclear mass for Z (a#of protons), N (a#of neutrons), & S (a#of lambdas)
    … …  
    16251634    return 0.;     // Photon mass
    16261635  }
     1636  else if(!N&&!Z&&S==1) return mL;
     1637  else if(!N&&Z==1&&!S) return mP;
     1638  else if(N==1&&!Z&&!S) return mN;
    16271639  else if(!N&&!Z&&S>1) return mL*S+eps;
    16281640  else if(!N&&Z>1&&!S) return mP*Z+eps;

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4Quasmon.cc

    r1228 r1315  
    2828//
    2929//
    30 // $Id: G4Quasmon.cc,v 1.121 2009/12/03 18:09:07 mkossov Exp $
    31 // GEANT4 tag $Name: geant4-09-03 $
     30// $Id: G4Quasmon.cc,v 1.125 2010/06/10 08:37:27 mkossov Exp $
     31// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    3232//
    3333//      ---------------- G4Quasmon ----------------
    … …  
    596596    G4LorentzVector cr4Mom=zeroLV;        // 4-momentum prototype for the ColResQuasmon
    597597    G4int kCount =0;                      // Counter of attempts of k for hadronization
     598    //
     599    //G4int qCountMax=27;                   // Try different q to come over CoulBar or SepE
     600    //G4int qCountMax=12;                   // Try different q to come over CoulBar or SepE
     601    //G4int qCountMax=9;                    // Try different q to come over CoulBar or SepE
     602    //G4int qCountMax=3;                    // Try different q to come over CoulBar or SepE
     603    G4int qCountMax=1;                    // Try different q to come over CoulBar or SepE
     604    if(excE > diPiM) qCountMax=(G4int)(excE/mPi0); // Try more for big excess
     605    //
    598606    //G4int kCountMax=27;
    599607    //G4int kCountMax=9;
    600608    //G4int kCountMax=3;                    // Try different k if they are below minK
    601609    G4int kCountMax=1;                    // "No reson to increase it"
     610    //G4int kCountMax=qCountMax;            // "No reson to increase it"
    602611    //G4int kCountMax=0;                    //@@ *** Close search for the minimum k ***
    603     //
    604     //G4int qCountMax=27;                   // Try different q to come over CoulBar or SepE
    605     G4int qCountMax=12;                   // Try different q to come over CoulBar or SepE
    606     //G4int qCountMax=9;                    // Try different q to come over CoulBar or SepE
    607     //G4int qCountMax=3;                    // Try different q to come over CoulBar or SepE
    608     //G4int qCountMax=1;                    // Try different q to come over CoulBar or SepE
    609612    //
    610613    //G4int pCountMax=27;                   //Try differentHadrons(Parents) forBetterRecoil
    … …  
    613616    G4int pCountMax=1;                    //Try differentHadrons(Parents) forBetterRecoil
    614617    //if(envA>0) pCountMax=3;
    615     //if(envA>0&&envA<31)
    616     if(envA>0&&envA<61)
    617     {
    618       //pCountMax=60/envA;
    619       pCountMax=120/envA;
    620       //kCountMax=pCountMax;                // See above !
    621     }
     618    if(envA>0&&envA<19) pCountMax=36/envA;
     619    //if(envA>0&&envA<31) pCountMax=60/envA;
     620    //if(envA>0&&envA<61) pCountMax=120/envA;
    622621    G4bool gintFlag=false;                // Flag of gamma interaction with one quark
    623622    while(kCount<kCountMax&&kCond)
    … …  
    11091108          ClearQuasmon();                   // This Quasmon is done
    11101109        }
    1111 #ifdef pdebug
     1110#ifdef debug
    11121111        G4cout<<"***G4Q::HQ:dE="<<excE<<">minK="<<minK<<",Env="<<theEnvironment<<",k="<<kLS
    11131112              <<",Q="<<valQ<<quasM<<", nQ="<<nQuasms<<G4endl;
    … …  
    11771176        if(!sPDG&&theEnvironment.GetPDG()!=NUCPDG&&totBN>1&&totMass>totM&&totS>=0)
    11781177        {
    1179 #ifdef ppdebug
     1178#ifdef pdebug
    11801179          G4cout<<"G4Quas::HQ:NEED-EVAP-1:Q="<<q4Mom<<valQ<<",En="<<theEnvironment<<G4endl;
    11811180          throw G4QException("G4Quasmon::HadronizeQuasmon: Why Fail? (1)"); //@@ TMP
    … …  
    12031202            if(totBN>1&&totMass>totM&&totS>=0)
    12041203            {
    1205 #ifdef ppdebug
     1204#ifdef pdebug
    12061205              G4cout<<"G4Q::HQ:NEED-EVAP-2:Q="<<q4Mom<<valQ<<",E="<<theEnvironment<<G4endl;
    12071206              throw G4QException("G4Quasmon::HadronizeQuasmon: Why Fail? (2)"); //@@ TMP
    … …  
    12201219        else if(nQuasms>1)
    12211220        {
    1222 #ifdef ppdebug
     1221#ifdef pdebug
    12231222          G4cout<<"G4Quas::HQ:NEED-EVAP-0:Q="<<q4Mom<<valQ<<",En="<<theEnvironment<<G4endl;
    12241223          throw G4QException("G4Quasmon::HadronizeQuasmon: Why Fail? (0)"); //@@ TMP
    … …  
    12421241          else if(rPDG!=NUCPDG&&totBN>1&&totMass>totM&&totS>=0)  // ===> "Evaporation" case
    12431242          {
    1244 #ifdef ppdebug
     1243#ifdef pdebug
    12451244            G4QContent nTotQC=totQC-neutQC;
    12461245            G4QNucleus nTotN(nTotQC);         // PseudoNucleus for TotalResidual to neutron
    … …  
    15611560            qCond=true;
    15621561            G4int qCount=0;
    1563 #ifdef ppdebug
     1562#ifdef pdebug
    15641563            G4double dM=0.;
    15651564            if(sPDG==90001001) dM=2.25;       // Binding energy of the deuteron ???
    … …  
    15731572            if(resTNM && totMass<resTNM+sMass)// Probably it never takes place
    15741573            {
    1575 #ifdef ppdebug
     1574#ifdef pdebug
    15761575              G4cout<<"***G4Quasmon::HadronizeQuasmon:***PANIC#1***TotalDE="<<excE<<"< bE="
    15771576                    <<sMass-pMass-dM<<", dM="<<dM<<", sM="<<sMass<<", bM="<<pMass<<G4endl;
    … …  
    16441643              //**PSLtest**//G4double z= pow((loli+hili)/2,pw);// ***TMP Mid q-RANDOMIZE***
    16451644              G4double ctkk=1.-(dex/(1.-z)-pMass)/kLS;// cos(theta_k,kappa) in LS
    1646 #ifdef ppdebug
     1645#ifdef pdebug
    16471646              if(qCount) G4cout<<"G4Q::HQ:qC="<<qCount<<",ct="<<ctkk<<",M="<<pMass<<",z="
    16481647                               <<z<<",zl="<<pow(loli,pw)<<",zh="<<pow(hili,pw)<<",dE="
    … …  
    19221921        {
    19231922          hsflag=true;                          // Decay in ThisHadron/Fragment+(SHadron)
    1924 #ifdef pdebug
     1923#ifdef debug
    19251924          G4cout<<"G4Q::HQ:NO,hsfl=1,kt="<<kt<<"<"<<minSqB<<" or M="<<tM<<"<"<<rtM<<G4endl;
    19261925#endif
    … …  
    21382137            PMEMhsfl=hsflag;
    21392138            PMEMnucf=nucflag;
    2140 #ifdef pdebug
     2139#ifdef debug
    21412140            G4cout<<"G4Q::HQ:RemTheBest rPDG="<<rPDG<<",sPDG="<<sPDG<<",kt="<<kt<<G4endl;
    21422141#endif
    … …  
    21692168            PMEMhsfl=hsflag;
    21702169            PMEMnucf=nucflag;
    2171 #ifdef pdebug
     2170#ifdef debug
    21722171            G4cout<<"G4Q::HQ:RemTheFirst rPDG="<<rPDG<<",sPDG="<<sPDG<<",kt="<<kt<<G4endl;
    21732172#endif
    … …  
    22052204      pCount++;
    22062205    } // End of the WHILE of the parent choice
    2207 #ifdef pdebug
     2206#ifdef debug
    22082207    G4cout<<"G4Q::HQ:>rPDG="<<rPDG<<curQ<<",sPDG="<<sPDG<<",kt="<<kt<<",F="<<fprob
    22092208          <<",totQC="<<totQC<<",sQC="<<sQC<<G4endl;
    … …  
    22342233       ) aMass=0.;         // No Pi0 cond.(eg in NucE)
    22352234
    2236 #ifdef pdebug
     2235#ifdef debug
    22372236    G4cout <<"G4Q::HQ:Is hsfl="<<hsflag<<" or fdul="<<fdul<<" or [rM="<<rMass<<"<"<<reMass
    22382237           <<" + "<<aMass<<" or rM2="<<reTNM2<<" < miM2="<<tmpTM2<<" and ePDG="<<envPDG
    … …  
    22642263        if(rQ4Mom==zeroLV)
    22652264        {
    2266 #ifdef ppdebug
     2265#ifdef pdebug
    22672266          G4cout<<"G4Q::HQ:NEEDS-EVAP-5,Q="<<q4Mom<<valQ<<",QEnv="<<theEnvironment<<G4endl;
    22682267          throw G4QException("G4Quasmon::HadronizeQuasmon: Why Fail? (5)"); //@@ TMP
    … …  
    23762375        //else if(2>3) // ********** Forced Evaporation is Closed ************
    23772376        {
    2378 #ifdef ppdebug
     2377#ifdef pdebug
    23792378          //@@ May be recalculate hadronization ??
    23802379          G4double fraM=fr4Mom.m();
    … …  
    26782677            if(CheckGroundState()) ClearQuasmon();// This Quasmon is done
    26792678            //if(CheckGroundState(true)) KillQuasmon();// This Quasmon is done
    2680 #ifdef ppdebug
     2679#ifdef pdebug
    26812680            G4cerr<<"G4Q::HQ:NeedsEvap7:s="<<sPDG<<",Q="<<q4Mom<<valQ<<",r="<<rPDG<<G4endl;
    26822681            throw G4QException("G4Quasmon::HadronizeQuasmon: Why Fail? (7)"); //@@ TMP
    … …  
    27252724      //if(2>3)                                      // This Evaporation channel is closed
    27262725      {
    2727 #ifdef ppdebug
     2726#ifdef pdebug
    27282727        G4cerr<<"***G4Q::HQ:E-8:RQM+SM="<<reMass+sMass<<">QM="<<quasM<<", sCB="<<sCB
    27292728              <<" + rCB="<<rCB<<" + rM="<<reMass<<" + sMass="<<sMass<<" + eM="<<envM<<" = "
    … …  
    28022801      {
    28032802        G4QContent resNQC=totQC-sQC;          // Quark Content of the totNucleus-Fragment
    2804 #ifdef pdebug
     2803#ifdef debug
    28052804        G4cerr<<"---Warning---G4Q::HQ:M="<<tmM<<"=>rPDG="<<rPDG<<"(rM="<<reMass<<")+sPDG="
    28062805              <<sPDG<<"(sM="<<sMass<<")="<<sum<<",resNQC="<<resNQC<<G4endl;
    … …  
    30263025          G4QNucleus nucZ(lanQC-Pi0QC);            // New neucleus for the residual for Pi-
    30273026          G4double nreZ=nucZ.GetGSMass();          // Mass of the residual for Pi-
    3028 #ifdef pdebug
     3027#ifdef debug
    30293028          G4cout<<"G4Q::HQ:LsPDG="<<sPDG<<",rPDG="<<rPDG<<",Z="<<nucZ<<",M="<<nucM<<G4endl;
    30303029#endif
    … …  
    32203219      if(sKE<sCB)                             // => "KinEn is below CB, try once more" case
    32213220      {
    3222 #ifdef ppdebug
     3221#ifdef pdebug
    32233222        G4cout<<"****G4Q::HQ:E-9: sKE="<<sKE<<"<sCB="<<sCB<<G4endl;
    32243223        throw G4QException("G4Quasmon::HadronizeQuasmon: Why Fail? (9)"); //@@ TMP
    … …  
    32533252        {
    32543253          //if(totBN>1&&totMass>totM&&totS>=0)         //@@ ??
    3255 #ifdef ppdebug
     3254#ifdef pdebug
    32563255          G4cout<<"G4Q::HQ:*NEEDS-EVAP-10* M="<<resTotN4Mom.m()<<"<miM="<<resTotNM<<G4endl;
    32573256          throw G4QException("G4Quasmon::HadronizeQuasmon: Why Fail? (10)"); //@@ TMP
    … …  
    33293328               rK0PDG == NUCPDG )
    33303329          {
    3331 #ifdef ppdebug
     3330#ifdef pdebug
    33323331            G4cout<<"G4Q::HQ:***PANIC#2***tM="<<totMass<<"<KM="<<mK<<","<<mK0<<",rM="<<rKPM
    33333332                  <<","<<rK0M<<",d="<<mK+rKPM-totMass<<","<<mK0+rK0M-totMass<<G4endl;
    … …  
    34453444                <<",d="<<rMass+sMass-q4Mom.m()<<G4endl;
    34463445#endif
    3447 #ifdef ppdebug
     3446#ifdef pdebug
    34483447          throw G4QException("G4Quasmon::HadronizeQuasmon: Why PANIC? (3)"); //@@ TMP
    34493448#endif
    … …  
    36363635  G4QHadronVector* theFragments = new G4QHadronVector; // user is responsible to delete!
    36373636  G4int nHadrs=theQHadrons.size();
    3638 #ifdef pdebug
     3637#ifdef debug
    36393638  G4cout<<"G4Q::DecayQuasmon:After decay (FillHadronVector byItself) nH="<<nHadrs<<G4endl;
    36403639#endif
    … …  
    36443643    theFragments->push_back(curHadr);                  // (user must delete)
    36453644  }
    3646 #ifdef ppdebug
     3645#ifdef pdebug
    36473646  else G4cerr<<"*******G4Quasmon::DecayQuasmon: *** Nothing is in the output ***"<<G4endl;
    36483647#endif
    … …  
    37203719    else if(thePDG==89999002) thePDG=1114;  // Delta-(resonance) === bary-DECUPLET/OCTET
    37213720    else if(thePDG==90001999) thePDG=2224;  // Delta++(res)(Delta0&Delta+ a covered by n&p)
     3721    else if(thePDG==91000000) thePDG=3122;  // Lambda
    37223722    else if(thePDG==90999001) thePDG=3112;  // Sigma-
    37233723    else if(thePDG==91000999) thePDG=3222;  // Sigma+ (Sigma0 iz covered by Lambda)
    … …  
    37693769    //G4double fragMas=qH->GetMass();          // GrStMass of the nuclear fragment (wrong?)
    37703770    G4QNucleus qNuc(t,thePDG);               // Make a Nucleus out of the Hadron
    3771     G4double GSMass =qNuc.GetGSMass();       // GrState Mass of the nuclear fragment
     3771    // @@ Probably, when nucleus is initialized, the mass is not initialized ? Was OK! Why?
     3772    //G4double GSMass =qNuc.GetGSMass();       // GrState Mass of the nuclear fragment (?)
     3773    G4double GSMass = G4QPDGCode(thePDG).GetMass(); // More robust definition
    37723774    G4QContent totQC=qNuc.GetQCZNS();        // Total Quark Content of Residual Nucleus
    37733775    G4int    nN     =qNuc.GetN();            // A#of neutrons in the Nucleus
    … …  
    37783780    G4cout<<"G4Quasm::FillHadrVect:Nucl="<<qNuc<<",nPDG="<<thePDG<<",GSM="<<GSMass<<G4endl;
    37793781#endif
    3780     if((nN<0||nZ<0||nS<0)&&bA>0)             // => "Anti-strangeness or ISOBAR" case
     3782    if((nN<0 || nZ<0 || nS<0) && bA>0)       // => "Anti-strangeness or ISOBAR" case
    37813783    {
    37823784      G4double m1=mPi;                       // Prototypes for the nZ<0 case
    … …  
    38623864      else
    38633865      {
    3864 #ifdef pdebug
     3866#ifdef debug
    38653867        G4cerr<<"-Warning-G4Q::FillHVec:PDG="<<thePDG<<"("<<t.m()<<","<<fragMas<<") < Mes="
    38663868              <<PDG1<<"("<<m1<<") + ResA="<<PDG2<<"("<<m2<<"), d="<<fragMas-m1-m2<<G4endl;
    … …  
    39113913        else lResM  =resN.GetMZNS();           // min mass of the Residual Nucleus
    39123914      }
    3913 #ifdef pdebug
     3915#ifdef debug
    39143916      G4cout<<"G4Quasm::FillHadrVec:rP="<<pResPDG<<",rN="<<nResPDG<<",rL="<<lResPDG<<",nN="
    39153917            <<nN<<",nZ="<<nZ<<",nL="<<nS<<",totM="<<fragMas<<",n="<<fragMas-nResM-mNeut
    … …  
    39193921           (bA > 1 && ( (nN > 0 && fragMas > nResM+mNeut) ||
    39203922                        (nZ > 0 && fragMas > pResM+mProt) ||
    3921    (nS > 0 && fragMas > lResM+mLamb) ) ) )
     3923                        (nS > 0 && fragMas > lResM+mLamb) ) ) )
    39223924      {
    39233925        G4int barPDG = 90002002;
    … …  
    39263928        G4double resM= mAlph;
    39273929
    3928  if (fragMas > nResM+mNeut) {  // Can radiate a neutron (priority 1)
     3930        if (fragMas > nResM+mNeut) {  // Can radiate a neutron (priority 1)
    39293931          barPDG = 90000001;
    39303932          resPDG = nResPDG;
    39313933          barM= mNeut;
    39323934          resM= nResM;
    3933  }
    3934  else if(fragMas>pResM+mProt)  // Can radiate a proton (priority 2)
    3935  {
     3935        }
     3936        else if(fragMas>pResM+mProt)  // Can radiate a proton (priority 2)
     3937        {
    39363938          barPDG=90001000;
    39373939          resPDG=pResPDG;
    39383940          barM  =mProt;
    39393941          resM  =pResM;
    3940  }
    3941  else if(fragMas>lResM+mLamb)  // Can radiate a Lambda (priority 3)
    3942      {
     3942        }
     3943        else if(fragMas>lResM+mLamb)  // Can radiate a Lambda (priority 3)
     3944        {
    39433945          barPDG=91000000;
    39443946          resPDG=lResPDG;
    39453947          barM  =mLamb;
    39463948          resM  =lResM;
    3947  }
    3948         else if(thePDG!=90004004&&fragMas>GSMass)// If it's not Be8 decay in gamma
    3949  {
     3949        }
     3950        else if(thePDG!=90004004 && fragMas>GSMass)// If it's not Be8 decay in gamma
     3951        {
    39503952          barPDG=22;
    39513953          resPDG=thePDG;
    39523954          barM  =0.;
    39533955          resM  =pResM;
    3954  }
     3956        }
    39553957        else if(thePDG!=90004004)
    3956  {
     3958        {
    39573959          G4cerr<<"***G4Q::FillHadV:PDG="<<thePDG<<",M="<<fragMas<<"<GSM="<<GSMass<<G4endl;
    39583960          throw G4QException("***G4Quasmon::FillHadronVector: Below GSM but cann't decay");
    3959  }
    3960 
     3961        }
    39613962        G4LorentzVector a4Mom(0.,0.,0.,barM);
    39623963        G4LorentzVector b4Mom(0.,0.,0.,resM);
    … …  
    39813982      else
    39823983      {
    3983 #ifdef pdebug
     3984#ifdef debug
    39843985        G4cout<<"G4Quasm::FillHadrVect: Leave as it is"<<G4endl;
    39853986#endif
    … …  
    39873988      }
    39883989    }
    3989     else if (fragMas<GSMass)
     3990    else if (fragMas < GSMass)                 // Approximate equality was already checked
    39903991    {
    3991       G4cerr<<"***G4Q::FillHV:M="<<fragMas<<">GSM="<<GSMass<<",d="<<fragMas-GSMass<<G4endl;
     3992      G4cerr<<"***G4Quasmon::FillHV:M="<<fragMas<<">GSM="<<GSMass<<"(PDG="<<thePDG<<"),d="
     3993            <<fragMas-GSMass<<", NZS="<<nN<<","<<nZ<<","<<nS<<G4endl;
    39923994      //throw G4QException("*G4Quasmon::FillHadronVector:Mass is below theGroundStateVal");
    39933995      G4cout<<"***>>>G4Quasm::FillHadrVect: Leave as it is Instead of Exception"<<G4endl;
    39943996      theQHadrons.push_back(qH);              // Fill As Is  (delete equivalent)
    39953997    }
    3996     else if (bA==1&&fragMas>GSMass)
     3998    else if (bA==1 && fragMas>GSMass)
    39973999    {
    39984000      G4int gamPDG=22;
    … …  
    40204022    else                                   // ===> Evaporation of excited system
    40214023    {
    4022 #ifdef pdebug
     4024#ifdef ppdebug
    40234025      G4cout<<"G4Quasm::FillHadrVect:Evaporate "<<thePDG<<",tM="<<fragMas<<" > GS="<<GSMass
    40244026            <<qNuc.Get4Momentum()<<", m="<<qNuc.Get4Momentum().m()<<G4endl;
    … …  
    40354037      else
    40364038      {
    4037 #ifdef pdebug
     4039#ifdef debug
    40384040        G4cout<<"G4Q::FlHV:Done,b="<<bHadron->GetQPDG()<<",r="<<rHadron->GetQPDG()<<G4endl;
    40394041#endif
    40404042        delete qH;
    4041         FillHadronVector(bHadron);          // Fill Evapor. Baryon (delete equivalent)
    4042         FillHadronVector(rHadron);          // Fill Residual Nucl. (delete equivalent)
     4043        if(bHadron->GetPDGCode() < -1111)
     4044        {
     4045          G4QHadronVector* tmpQHadVec=DecayQHadron(bHadron); // (delete equivalent)
     4046          G4int tmpS=tmpQHadVec->size();
     4047          theQHadrons.resize(tmpS+theQHadrons.size());       // Resize theQHadrons length
     4048          copy( tmpQHadVec->begin(), tmpQHadVec->end(), theQHadrons.end()-tmpS);
     4049          tmpQHadVec->clear();
     4050          delete tmpQHadVec;  // Who calls DecayQHadron is responsible for clear & delete
     4051        }
     4052        else FillHadronVector(bHadron);          // Fill Evapor. Baryon (delete equivalent)
     4053        if(rHadron->GetPDGCode() < -1111)
     4054        {
     4055          G4QHadronVector* tmpQHadVec=DecayQHadron(rHadron); // (delete equivalent)
     4056          G4int tmpS=tmpQHadVec->size();
     4057          theQHadrons.resize(tmpS+theQHadrons.size());       // Resize theQHadrons length
     4058          copy( tmpQHadVec->begin(), tmpQHadVec->end(), theQHadrons.end()-tmpS);
     4059          tmpQHadVec->clear();
     4060          delete tmpQHadVec;  // Who calls DecayQHadron is responsible for clear & delete
     4061        }
     4062        else FillHadronVector(rHadron);          // Fill Residual Nucl. (delete equivalent)
    40434063      }
    40444064    }
    … …  
    40744094{
    40754095  //gives k>kMin QParton Momentum for the current Quasmon
    4076 #ifdef pdebug
     4096#ifdef debug
    40774097  //G4cout<<"G4Quasmon::GetQPartonMom:**called**mR="<<sqrt(mR2)<<",mC="<<sqrt(mC2)<<G4endl;
    40784098  G4cout<<"G4Quas::GetQPartonMomentum:***called*** kMax="<<kMax<<",mC="<<sqrt(mC2)<<G4endl;
    … …  
    41014121    throw G4QException("G4Quasmon::GetQPartonMomentum: Can not generate quark-parton");   
    41024122  }
    4103 #ifdef pdebug
     4123#ifdef debug
    41044124  G4cout<<"G4Q::GetQPM: kLim="<<kLim<<",kMin="<<kMin<<",kMax="<<kMax<<",nQ="<<nOfQ<<G4endl;
    41054125#endif
    … …  
    41544174    if(f>27.)
    41554175    {
    4156 #ifdef pdebug
     4176#ifdef debug
    41574177      G4cout<<"G4Q::GetQPMom: f="<<f<<" is changed to 99"<<G4endl;
    41584178#endif
    … …  
    41674187    G4double c  = p*r*(1.+nx); // vRndm=(1-x)**n*(1+n*x) is the equation too be solved
    41684188    G4double od = c - vRndm;   // the old Residual Difference
    4169 #ifdef pdebug
     4189#ifdef debug
    41704190    G4cout<<"G4Q::GetQPMom:>>>First x="<<x<<", n="<<n<<", f="<<f<<", d/R(first)="
    41714191            <<od/vRndm<<G4endl;
    … …  
    42044224        if (f>=27.) f=27.;
    42054225      }
    4206 #ifdef pdebug
     4226#ifdef debug
    42074227      G4cout<<"G4Q::GetQPMom: Iter#"<<it<<": (c="<<c<<" - R="<<vRndm<<")/R ="<<d/vRndm
    42084228            <<", x="<<x<<", f="<<f<<G4endl;
    … …  
    42194239      it++;
    42204240    }
    4221 #ifdef pdebug
     4241#ifdef debug
    42224242    if(it>nitMax) G4cout<<"G4Q::GetQPMom: a#of iterations > nitMax="<<nitMax<<G4endl;
    42234243#endif
    … …  
    42994319  else if ( ev && nOfQ<3) nOfQ=3; // #of valence quarks is odd
    43004320  //
    4301 #ifdef pdebug
     4321#ifdef debug
    43024322  G4cout<<"G4Q::Calc#ofQP:QM="<<q4Mom<<qMass<<",T="<<Temperature<<",QC="<<valQ<<",n="<<nOfQ
    43034323         <<G4endl;
    … …  
    43144334  G4int nMaxStrangeSea=static_cast<int>((qMass-stran*mK0)/(mK0+mK0));//KK is min for s-sea
    43154335  if (absb) nMaxStrangeSea=static_cast<int>((qMass-absb)/672.); //LambdaK is min for s-sea
    4316 #ifdef pdebug
     4336#ifdef debug
    43174337  G4cout<<"G4Q::Calc#ofQP:"<<valQ<<",INtot="<<valc<<",nOfQ="<<nOfQ<<",SeaPairs="<<nSeaPairs
    43184338        <<G4endl;
    … …  
    43234343    if(nSeaPairs>0)valQ.IncQAQ(nSeaPairs,SSin2Gluons);
    43244344    else    morDec=valQ.DecQAQ(-nSeaPairs);
    4325 #ifdef pdebug
     4345#ifdef debug
    43264346    if(morDec) G4cout<<"G4Q::Calc#ofQP: "<<morDec<<" pairs can be reduced more"<<G4endl;
    43274347#endif
    … …  
    43314351    if(sSea>nMaxStrangeSea) // @@@@@@@ Too many strange sea ??????
    43324352    {
    4333 #ifdef pdebug
     4353#ifdef debug
    43344354      G4cout<<"G4Q::Calc#ofQP:**Reduce** S="<<sSea<<",aS="<<asSea<<",maxS="<<nMaxStrangeSea
    43354355            <<G4endl;
    … …  
    43464366  //if(nOfQ<nmin) nOfQ=nmin;
    43474367  // --- End of Temporary
    4348 #ifdef pdebug
     4368#ifdef debug
    43494369  G4cout<<"G4Quasmon::Calc#ofQP: *** RESULT IN*** nQ="<<nOfQ<<", FinalQC="<<valQ<<G4endl;
    43504370#endif
    … …  
    44714491  G4QContent envQC=theEnvironment.GetQCZNS(); // QuarkContent of the CurrentNuclearEnviron.
    44724492  //////G4double addPhoton=phot4M.e();              // PhotonEn for capture by quark-parton
    4473 #ifdef pdebug
     4493#ifdef debug
    44744494  G4int absb = abs(qBar);                     // Abs BaryNum of Quasmon
    44754495  G4int maxB = theEnvironment.GetMaxClust();  // Maximum BaryNum for clusters
    … …  
    44794499  // ================= Calculate probabilities for candidates
    44804500  unsigned nHC=theQCandidates.size();
    4481 #ifdef pdebug
     4501#ifdef debug
    44824502  G4cout<<"G4Q::CHP: *** nHC="<<nHC<<G4endl;
    44834503#endif
    … …  
    45094529      //G4bool pos=curCand->GetPossibility();
    45104530#ifdef pdebug
    4511       G4bool pPrint = abs(cPDG)%10<3 && cPDG<80000000 ||cPDG==90001000||cPDG==90000001||
    4512         cPDG==90000002||cPDG==90001001||cPDG==90001002||cPDG==90002001||cPDG==90002002;
     4531      //G4bool pPrint = abs(cPDG)%10<3 && cPDG<80000000 ||cPDG==90001000||cPDG==90000001||
     4532      //  cPDG==90000002||cPDG==90001001||cPDG==90001002||cPDG==90002001||cPDG==90002002;
     4533      //G4bool pPrint = cPDG==2212 || cPDG==2112 ||cPDG==90001000||cPDG==90000001;
     4534      G4bool pPrint = false;
    45134535      if(pPrint) G4cout<<"G4Q::CHP:==****==>>>c="<<cPDG<<",dUD="<<dUD<<",pos="<<pos<<",eA="
    45144536                       <<envA<<",tM="<<totMass<<" > tmpTM+frM="<<tmpTM+frM<<G4endl;
    … …  
    45314553        //G4int cNQ= candQC.GetTot()+3*baryn-4; // A#of q-partons+b+2*(b-1)-2 (choc q-link)
    45324554        G4double resM=0.;                       // Prototype for minMass of residual hadron
    4533 #ifdef pdebug
     4555#ifdef debug
    45344556        if(pPrint)G4cout<<"G4Q::CHP:B="<<baryn<<",C="<<cC<<",CB="<<CB<<",#q="<<cNQ<<G4endl;
    45354557#endif
    … …  
    45544576            iQmax=1;
    45554577          }
    4556 #ifdef pdebug
     4578#ifdef debug
    45574579          if(pPrint)G4cout<<"G4Q::CHP:***!!!***>>F="<<cPDG<<",mF="<<frM<<",iq:"<<iQmin<<","
    45584580                         <<iQmax<<",kLS="<<kLS<<",kQCM="<<kVal<<",eA="<<envA<<G4endl;
    … …  
    45654587            {
    45664588              qFact=4.;
    4567 #ifdef pdebug
     4589#ifdef debug
    45684590              if(pPrint) G4cout<<"G4Q::CHP:photon cap(gaF) is enhanced for Uquark"<<G4endl;
    45694591#endif
    … …  
    46294651#ifdef pdebug
    46304652                  if(pPrint) G4cout<<"G4Q::CHP:parentPossibility="<<possib<<",pZ="<<charge
    4631                                    <<",pN="<<barot-charge<<", cPDG="<<cPDG<<G4endl;
     4653                                   <<" <= envZ="<<envZ<<", pN="<<barot-charge<<" <= envN="
     4654                                   <<envN<<", cPDG="<<cPDG<<G4endl;
    46324655#endif
    46334656                  if(possib && charge<=envZ && barot-charge<=envN)
    … …  
    46944717                    if ( (!piF && first && baryn < 3) ||
    46954718                         (!piF && !first && baryn < 5 ) ||
    4696                          (piF && abs(dS) < 3) ) // Isotopic focusing for AtRest Reactions
     4719                         ( piF && abs(dS) < 3) ) // Isotope Focusing for AtRest Reactions
    46974720                    //if(!qIso&&!dC||qIso>0&&dC<0||qIso<0&&dC>0)//MediumIsoFocusingForAll
    46984721                    //if(abs(dS)<3) // Universal IsotopeFocusing(<3) (Best for pi-capture)
    … …  
    47054728                      G4double boundM=parCand->GetEBMass();//EnvironBoundMass ofParentClust
    47064729                      G4double nucBM =parCand->GetNBMass();//TotNuclBoundMass ofParentClust
    4707 #ifdef pdebug
     4730#ifdef debug
    47084731                      if(pPrint) G4cout<<"G4Q::CHP:c="<<cPDG<<",p="<<parPDG<<",bM="<<boundM
    47094732                                       <<",i="<<is<<",adr="<<parCand<<",pPP="<<pPP<<G4endl;
    … …  
    47264749                      if(resPDG && minM>0.) // Kinematical analysis of hadronization
    47274750                      {
    4728 #ifdef pdebug
     4751#ifdef debug
    47294752                        if(pPrint) G4cout<<"G4Q::CHP:fM="<<frM<<",bM="<<boundM<<",rM="
    47304753                                         <<tmpTM<<",tM="<<totMass<<G4endl;
    … …  
    47374760                        G4double dked=ked+ked;
    47384761                        //////G4double dkedC=dked-CB-CB;
    4739                         //G4double kd =kLS-nDelta; //For TotalNucleus (includingQuasmon)
    4740                         //G4double dkd=kd+kd;
    4741                         //////////G4double dkdC=dkd-CB-CB;
     4762                        G4double kd =kLS-nDelta; //For TotalNucleus (includingQuasmon)
     4763                        G4double dkd=kd+kd;
     4764                        //G4double dkdC=dkd-CB-CB;
    47424765                        G4double dkLS=kLS+kLS;
    47434766                        //G4double Em=(E-eDelta)*(1.-frM/totMass);
    47444767                        //G4double Em=(E-nDelta)*(1.-frM/totMass);
    4745                         //G4double Em=(E-nDelta-CB)*(1.-frM/totMass);
     4768                        G4double Em=(E-nDelta-CB)*(1.-frM/totMass);
    47464769                        // *** START LIMITS ***
    47474770                        G4double ne=1.-dked/(boundM+dkLS);// qmin=DEFOULT=bM*(k-de)/(bM+2k)
    47484771                        G4double kf=1.;
    47494772                        if(ne>0.&&ne<1.)kf=pow(ne,cNQ);
    4750 #ifdef pdebug
     4773#ifdef debug
    47514774                        if(pPrint)G4cout<<"G4Q::CHP:<qi_DEF>="<<ne<<",k="<<kf<<",dk="<<dked
    47524775                                        <<",dkLS="<<dkLS<<",M="<<boundM<<",C="<<CB<<G4endl;
    … …  
    47674790                        G4double rtE=rtEP-rMo;        // Energy of RealTotColouredResidSyst
    47684791                        ////////////G4double rtEMP=rtE-rMo;
    4769 #ifdef pdebug
     4792#ifdef debug
    47704793                        G4QContent tmpEQ=envQC-parQC;//QuarkContent for ResidualEnvironment
    47714794                        if(pPrint) G4cout<<"G4Q::CHP:RN="<<tmpEQ<<"="<<envM-boundM<<"=eM="
    … …  
    47884811                        G4double minBM2=minBM*minBM+.1;
    47894812                        G4double minM2=minM*minM+.1;
    4790 #ifdef pdebug
     4813#ifdef debug
    47914814                        G4double ph=kf;                     // Just for printing
    47924815                        if(pPrint) G4cout<<"G4Q::CHP:M2="<<minM2<<",R="<<rQ2<<",m="<<mintM2
    … …  
    47994822                        G4double nc=1.-(dkLS-E-E+CB+CB)/boundM;   // q_min=k-E+CB
    48004823                        G4double newl=0.;
    4801 #ifdef pdebug
     4824#ifdef debug
    48024825                        if(pPrint) G4cout<<"G4Q::CHP:qi_k-E="<<nc<<",k="<<kLS<<",E="<<E
    48034826                                         <<",M="<<boundM<<G4endl;
    … …  
    48114834                        else if(nc <= 0.) kf=0.;
    48124835
    4813                         //G4double nk=1.-(dkd-Em-Em)/boundM;  // q_min=(k-delta)-E*(M-m)/M
    4814 #ifdef pdebug
    4815                         //if(pPrint) G4cout<<"G4Q::CHP:qi_R="<<nk<<",kd="<<kd<<",E="<<Em
    4816                         //               <<",M="<<totMass<<G4endl;
    4817 #endif
    4818                         //if(nk > 0. && nk < 1. && nk < ne)
    4819                         //{
    4820                         //  ne=nk;
    4821                         //  newh=pow(nk,cNQ);
    4822                         //  if(newh<kf) kf=newh;
    4823                         //}
    4824                         //else if(nk <= 0.) kf=0.;
     4836                        G4double nk=1.-(dkd-Em-Em)/boundM;  // q_min=(k-delta)-E*(M-m)/M
     4837#ifdef debug
     4838                        if(pPrint) G4cout<<"G4Q::CHP:qi_R="<<nk<<",kd="<<kd<<",E="<<Em
     4839                                       <<",M="<<totMass<<G4endl;
     4840#endif
     4841                        if(nk > 0. && nk < 1. && nk < ne)
     4842                        {
     4843                          ne=nk;
     4844                          newh=pow(nk,cNQ);
     4845                          if(newh<kf) kf=newh;
     4846                        }
     4847                        else if(nk <= 0.) kf=0.;
    48254848
    48264849                        //G4double mex=frM+Em;
    48274850                        //G4double sr=sqrt(mex*mex-frM2);//qmin=k-sqrt((m+E*(M-m)/M)^2-m^2)
    48284851                        //G4double np=1.-(dkLS-sr-sr)/boundM;
    4829 #ifdef pdebug
     4852#ifdef debug
    48304853                        //if(pPrint)G4cout<<"G4Q::CHP:qi_k-sr="<<np<<",sr="<<sr<<",m="<<mex
    48314854                        //                <<",M="<<frM<<G4endl;
    … …  
    48464869                        if(mix > frM) st=sqrt(mix*mix-frM2);
    48474870                        G4double nq=1.-(dkLS-st-st)/boundM;//qi=k-sq((m+E*(M-m)/M)^2-m^2)
    4848 #ifdef pdebug
     4871#ifdef debug
    48494872                        if(pPrint) G4cout<<"G4Q::CHP:qi_k-st="<<nq<<",st="<<st<<",m="
    48504873                                         <<mix<<",M="<<frM<<G4endl;
    … …  
    48734896                          //if(atrest) nz=1.-(mintM2-rtQ2+pmk*dkd)/(nucBM*(rtEP+pmk));
    48744897                          //else       nz=1.-(mintM2-rtQ2)/(nucBM*rtEP);
    4875 #ifdef pdebug
     4898#ifdef debug
    48764899                          if(pPrint) G4cout<<"G4Q::CHP:q="<<nz<<",a="<<atrest<<",M2="
    48774900                                           <<mintM2<<">"<<rtQ2<<G4endl;
    … …  
    49234946                          //if(atrest) nz=1.-(minBM2-rQ2+pmk*dkd)/(nucBM*(rEP+pmk));
    49244947                          //else       nz=1.-(minBM2-rQ2)/(nucBM*rEP);
    4925 #ifdef pdebug   
     4948#ifdef debug   
    49264949                          if(pPrint) G4cout<<"G4Q::CHP:q="<<nz<<",a="<<atrest<<",QM2="
    49274950                                           <<minM2<<">"<<rQ2<<G4endl;
    … …  
    49674990                          //if(atrest) nz=1.-(minM2-rQ2+pmk*dkd)/(nucBM*(rEP+pmk));
    49684991                          //else       nz=1.-(minM2-rQ2)/(nucBM*rEP);
    4969 #ifdef pdebug
     4992#ifdef debug
    49704993                          if(pPrint) G4cout<<"G4Q::CHP:q="<<nz<<",a="<<atrest<<",QM2="
    49714994                                           <<minM2<<">"<<rQ2<<G4endl;
    … …  
    49825005                        if(kf>1.)kf=1.;
    49835006                        G4double high = kf;                 // after this kf can be changed
    4984 #ifdef pdebug
     5007#ifdef debug
    49855008                        if(pPrint) G4cout<<"G4Q::CHP:"<<kf<<",minM2="<<minM2<<",rQ2="<<rQ2
    49865009                                         <<G4endl;
    … …  
    49945017                        if(lz>0.&&lz<1.)low=pow(lz,cNQ);
    49955018                        else if(lz>=1.)low=1.;
    4996 #ifdef pdebug
     5019#ifdef debug
    49975020                        G4double pl=low;                      // Just for printing
    49985021                        if(pPrint) G4cout<<"G4Q::CHP:<qa_DEF>="<<lz<<", eDel="<<eDelta
    … …  
    50015024                        // == (@@) Historical additional cuts for q_max ===
    50025025                        //G4double tms=kLS+nDelta+Em;
    5003                         ////G4double tms=kLS+eDelta+Em;        // The same don't change ***
    5004                         //G4double le=1.-(tms+tms)/boundM;     // q_max=k+delta+E*(M-m)/M
    5005 #ifdef pdebug
    5006                         //if(pPrint) G4cout<<"G4Q::CHP:qa_t="<<le<<",k="<<kLS<<",E="<<Em
    5007                         //                 <<",bM="<<boundM<<G4endl;
    5008 #endif
    5009                         //if(le>0.&&le<1.&&le>lz)
    5010                         //{
    5011                         //  lz=le;
    5012                         //  newl=pow(le,cNQ);
    5013                         //  if(newl>low) low=newl;
    5014                         //}
    5015                         //else if(le>=1.)low=1.;
     5026                        G4double tms=kLS+eDelta+Em;        // The same don't change ***
     5027                        G4double le=1.-(tms+tms)/boundM;     // q_max=k+delta+E*(M-m)/M
     5028#ifdef debug
     5029                        if(pPrint) G4cout<<"G4Q::CHP:qa_t="<<le<<",k="<<kLS<<",E="<<Em
     5030                                         <<",bM="<<boundM<<G4endl;
     5031#endif
     5032                        if(le>0.&&le<1.&&le>lz)
     5033                        {
     5034                          lz=le;
     5035                          newl=pow(le,cNQ);
     5036                          if(newl>low) low=newl;
     5037                        }
     5038                        else if(le>=1.)low=1.;
    50165039                        // === End of historical cuts
    50175040
    50185041                        //G4double lk=1.-(dkLS+E+E)/boundM;    // q_max=k+E
    50195042                        G4double lk=1.-(dkLS+E+E-CB-CB)/boundM;//qmax=k+E-CB(surfaceCond)
    5020 #ifdef pdebug
     5043#ifdef debug
    50215044                        if(pPrint) G4cout<<"G4Q::CHP:qa_k+E="<<lk<<",k="<<kLS<<",E="<<E
    50225045                                         <<",M="<<boundM<<G4endl;
    … …  
    50335056                        // === Instead one can try this ===
    50345057                        //G4double lq=1.-(dkLS+st+st)/boundM;//qm=k+sqrt((E*(M-m)/M)^2-m^2)
    5035 #ifdef pdebug
     5058#ifdef debug
    50365059                        //if(pPrint)G4cout<<"G4Q::CHP:qa_k+st="<<lq<<",st="<<st<<",m="<<mix
    50375060                        //                <<",M="<<frM<<G4endl;
    … …  
    50465069
    50475070                        // === The same as previous but "sr" instead of "st" ===
    5048                         //G4double lp=1.-(dkLS+sr+sr)/boundM;//qm=k+sqrt((E*(M-m)/M)^2-m^2)
    5049 #ifdef pdebug
    5050                         //if(pPrint) G4cout<<"G4Q::CHP:qa_k+sr="<<lp<<",sr="<<sr<<",m="
    5051                         //                 <<mex<<",M="<<frM<<G4endl;
    5052 #endif
    5053                         //if(lp>0.&&lp<1.&&lp>lz)
    5054                         //{
    5055                         //  lz=lp;
    5056                         //  newl=pow(lp,cNQ);
    5057                         //  if(newl>low) low=newl;
    5058                         //}
    5059                         //else if(lp>=1.)low=1.;
     5071                        G4double lp=1.-(dkLS+sr+sr)/boundM;//qm=k+sqrt((E*(M-m)/M)^2-m^2)
     5072#ifdef debug
     5073                        if(pPrint) G4cout<<"G4Q::CHP:qa_k+sr="<<lp<<",sr="<<sr<<",m="
     5074                                         <<mex<<",M="<<frM<<G4endl;
     5075#endif
     5076                        if(lp>0.&&lp<1.&&lp>lz)
     5077                        {
     5078                          lz=lp;
     5079                          newl=pow(lp,cNQ);
     5080                          if(newl>low) low=newl;
     5081                        }
     5082                        else if(lp>=1.)low=1.;
    50605083                        // ............................................................
    50615084                        //It's SpecificCoulombBarrierLimit forChargedParticles(canBeSkiped)
    5062                         if(totZ>cC)                         // ==> Check of Coulomb Barrier
     5085                        if(totZ>cC)                         // ==> Check CoulombBarrier
     5086                        //if(2>3)
    50635087                        {
    50645088                          G4double qmaCB=boundM-qMax;
    50655089                          //G4double qmaCB=nucBM-qMax;
    50665090                          G4double nz=1.-(qmaCB+qmaCB)/boundM;//q=Mb-Mf-CB+kLS,qM=Mf+CB-kLS
    5067 #ifdef pdebug
     5091#ifdef debug
    50685092                          if(pPrint) G4cout<<"G4Q::CHP:<qa_CB>="<<nz<<",m="<<qmaCB<<",CB="
    50695093                                           <<CB<<G4endl;
    … …  
    50785102                        // ***** End of restrictions *****
    50795103                        kf-=low;
    5080 #ifdef pdebug
     5104#ifdef debug
    50815105                        if(pPrint) G4cout<<"G4Q::CHP:>>"<<cPDG<<",l="<<low<<",h="<<high
    50825106                                         <<",ol="<<pl<<",oh="<<ph<<",nl="<<newl<<",nh="
    … …  
    51065130                            else          probab*=2;           // One S + three P
    51075131                          }
    5108 #ifdef pdebug
     5132#ifdef debug
    51095133                          if(pPrint)G4cout<<"G4Q::CHP:prob="<<probab<<",qF="<<qFact<<",iq="
    51105134                                          <<iq<<",oq="<<oq<<",Pho4M="<<phot4M<<",pUD="<<pUD
    … …  
    51895213            if(cPDG==223) comb*=2.; //@@
    51905214            if(cPDG==111&&npiz>0) comb*=(npiz+1); // Bose multyplication
    5191 #ifdef pdebug
     5215#ifdef debug
    51925216            if(abs(cPDG)<3) G4cout<<"G4Q::CHP:comb="<<comb<<",cmd="<<cmd<<",cmuu="<<cmu
    51935217                                  <<",cms="<<cms<<G4endl;
    … …  
    51985222          G4QContent resTQC = curQ+envQC;    // Total nuclear Residual Quark Content
    51995223          G4double resTM=G4QPDGCode(resTQC.GetSPDGCode()).GetMass();
    5200 #ifdef pdebug
     5224#ifdef debug
    52015225          G4bool priCon = aPDG < 10000 && aPDG%10 < 3;
    52025226          if(priCon) G4cout<<"G4Q::CHP:***>>cPDG="<<cPDG<<",cQC="<<candQC<<",comb="<<comb
    … …  
    52085232            resTM=mPi0;
    52095233          }
    5210 #ifdef pdebug
     5234#ifdef debug
    52115235          if(priCon) G4cout<<"G4Q::CHP:cPDG="<<cPDG<<",c="<<comb<<",rPDG/QC="<<resPDG<<curQ
    52125236                           <<",tM="<<totMass<<">"<<frM-CB+resTM<<"=fM="<<frM<<"+rM="<<resTM
    … …  
    52165240             ((resPDG > 80000000 && resPDG != 90000000) || resPDG<10000) )
    52175241          {
    5218 #ifdef pdebug
     5242#ifdef debug
    52195243            if(priCon) G4cout<<"G4Q::CHP:ind="<<index<<",qQC="<<valQ<<mQ<<",cPDG="<<cPDG
    52205244                             <<",rPDG="<<resPDG<<curQ<<G4endl;
    … …  
    52235247            else resM=G4QChipolino(curQ).GetMass(); // Chipolino mass for theResidualHadron
    52245248            G4int resQCode=G4QPDGCode(curQ).GetQCode();
    5225 #ifdef pdebug
     5249#ifdef debug
    52265250            if(priCon) G4cout<<"G4Q::CHP:rM/QC="<<resM<<curQ<<",E="<<envPDGC<<",rQC="
    52275251                             <<resQCode<<G4endl;
    … …  
    52355259              G4double rtM =rtN.GetMZNS();          // Min Mass of total residual Nucleus
    52365260              G4double bnRQ=rtM-envM;               // Bound mass of residual Quasmon
    5237 #ifdef pdebug
     5261#ifdef debug
    52385262              if(priCon) G4cout<<"G4Q::CHP: **Rec**,RQMass="<<bnRQ<<",envM="<<envM<<",rtM="
    52395263                               <<rtM<<G4endl;
    … …  
    52425266              if(bnRQ<resM) resM=bnRQ;
    52435267            }
    5244 #ifdef pdebug
     5268#ifdef debug
    52455269            if(aPDG<10000&&aPDG%10<3)
    52465270            //if(aPDG<10000&&aPDG%10<5)
    … …  
    52515275              G4double limM=mQ-resM;
    52525276              G4double rndM=GetRandomMass(cPDG,limM);// Candidate's Mass randomization
    5253 #ifdef pdebug
     5277#ifdef debug
    52545278              G4double cMass=G4QPDGCode(cPDG).GetMass();
    52555279              if(aPDG<10000&&aPDG%10<3)
    … …  
    52725296                if(qBar) zMin= mR2/mQ/(mQ-dk);       // z_min for Quasmon-Baryon @@ ?? @@
    52735297                G4double possibility=zMax-zMin;
    5274 #ifdef pdebug
     5298#ifdef debug
    52755299                if(priCon) G4cout<<"G4Q::CHP:M="<<rndM<<",ps="<<possibility<<",zMax="<<zMax
    52765300                                 <<",rHk="<<rHk<<",mQ="<<mQ<<",dk="<<dk<<",zMin="<<zMin
    … …  
    52855309                {
    52865310                  probability = vaf*(pow(zMax, vap)-pow(zMin, vap));
    5287 #ifdef pdebug
     5311#ifdef debug
    52885312                  if(priCon) G4cout<<"G4Q::CHP:#"<<index<<",mH2="<<mH2<<",nQ="<<nOfQ<<",p="
    52895313                                   <<probability<<",vf="<<vaf<<",vp="<<vap<<",zMax="<<zMax
    … …  
    53155339              else
    53165340              {
    5317 #ifdef pdebug
     5341#ifdef debug
    53185342                if(priCon) G4cout<<"G4Q::CHP:cM=0[cPDG"<<cPDG<<"],mQ/QC="<<mQ<<valQ<<",rM="
    53195343                                 <<resM<<curQ<<G4endl;
    … …  
    53235347            else
    53245348            {
    5325 #ifdef pdebug
     5349#ifdef debug
    53265350              if(priCon) G4cout<<"***G4Q::CHP: M=0, #"<<index<<valQ<<",cPDH="<<cPDG<<"+rP="
    53275351                               <<resPDG<<curQ<<G4endl;
    … …  
    53325356          {
    53335357            probability=0.;
    5334 #ifdef pdebug
     5358#ifdef debug
    53355359            if(priCon) G4cout<<"G4Q::CHP:"<<index<<valQ<<",PDG="<<cPDG<<"+r="<<resPDG<<curQ
    53365360                             <<":c=0(!) || tM="<<totMass<<"<"<<frM-CB+resTM<<" = fM="<<frM
    … …  
    53415365        }   // ---> End of Hadronic Case of fragmentation
    53425366        else probability=0.;
    5343 #ifdef pdebug
     5367#ifdef debug
    53445368        G4int aPDG = abs(cPDG);
    53455369        if(cPDG>90000000&&baryn<5||aPDG<10000&&aPDG%10<3) G4cout<<"G4Q::CHP:^^^cPDG="<<cPDG
    … …  
    53825406  G4LorentzVector reTLV=q4Mom;           // Prototyoe of the 4-Mom of the Residual Nucleus
    53835407  G4double resSMa=resQMa;                // Prototype of MinSplitMass of ResidualNucleus
    5384 #ifdef pdebug
     5408#ifdef debug
    53855409  G4cout<<"G4Q::CheckGS: EnvPDG="<<theEnvironment.GetPDG()<<",Quasmon="<<resQPDG<<G4endl;
    53865410#endif
    … …  
    53885412  {
    53895413    resEMa=theEnvironment.GetMZNS();     // GSMass of the Residual Environment
    5390 #ifdef pdebug
     5414#ifdef debug
    53915415    G4cout<<"G4Q::CheckGS: Environment Mass="<<resEMa<<G4endl;
    53925416#endif
    … …  
    54005424    //G4QNucleus tmpQN(valQ);              // TemporaryNucleus for the VacuumQuasmon
    54015425    bsCond = tmpQN.SplitBaryon();        // Possibility to split Fragment from the VacuumQ
    5402 #ifdef pdebug
     5426#ifdef debug
    54035427    G4cout<<"G4Q::CheckGS: No environment, theOnlyQ="<<tmpQN<<",bsCond="<<bsCond<<G4endl;
    54045428#endif
    … …  
    54535477  //if(resTMa>resSMa && (resEMa || bsCond)) return true;// Why not ?? @@ (see G4E the same)
    54545478  G4int nOfOUT = theQHadrons.size();      // Total #of QHadrons at this point
    5455 #ifdef pdebug
     5479#ifdef debug
    54565480  G4cout<<"G4Q::CheckGS: (totM="<<resTMa<<" < rQM+rEM="<<resSMa<<" || rEM="<<resEMa
    54575481        <<"=0 && "<<bsCond<<"=0) && n="<<nOfOUT<<" >0"<<G4endl;
    … …  
    54665490      G4double  hadrMa=hadr4M.m();        // Mass of the Last hadron (==GSMass)
    54675491      G4LorentzVector tmpTLV=reTLV+hadr4M;// Tot (ResidNucl+LastHadron) 4-Mom
    5468 #ifdef pdebug
     5492#ifdef debug
    54695493      G4cout<<"G4Q::CheckGS:YES,T="<<tmpTLV<<tmpTLV.m()<<">rM+hM="<<resSMa+hadrMa<<G4endl;
    54705494#endif
    … …  
    55245548          G4double   totQMa=G4QPDGCode(totPDG).GetMass(); // GS Mass of the ResidualNucleus
    55255549          G4double   totNMa=tmpTLV.m();                   // RealMass of TotResidualNucleus
    5526 #ifdef pdebug
     5550#ifdef debug
    55275551          G4cout<<"G4Q::CheckGS:NO, M="<<tmpTLV<<totNMa<<">"<<totQMa+hadrMa+prevMa<<G4endl;
    55285552#endif
    … …  
    55435567              theLast->Set4Momentum(hadr4M);
    55445568              thePrev->Set4Momentum(prev4M);
    5545 #ifdef pdebug
     5569#ifdef debug
    55465570           G4cout<<"G4Q::CheckGS: Yes, Check D4M="<<tmpTLV-hadr4M-prev4M-nuc4M<<G4endl;
    55475571#endif
    … …  
    55575581            G4double resPrevM=G4QPDGCode(resPPDG).GetMass();//GSM of ResidNucl for thePrevH
    55585582            //////G4bool    which = true;          // LastH is absorbed, PrevH is radiated
    5559 #ifdef pdebug
     5583#ifdef debug
    55605584            G4cout<<"G4Quasm::CheckGS: NO, tM="<<totNMa<<" > rp+l="<<resLastM+hadrMa
    55615585                  <<" || > rl+p="<<resPrevM+prevMa<<G4endl;
    … …  
    55915615              theEnvironment = vacuum;
    55925616              thePrev->Set4Momentum(prev4M);
    5593 #ifdef pdebug
     5617#ifdef debug
    55945618              G4cout<<"G4Q::CheckGS:Yes, Check D4M="<<tmpTLV-prev4M-nuc4M<<G4endl;
    55955619#endif
    … …  
    56145638  G4int        pap = 0;                          // --- particle
    56155639  if(thePDG<0) pap = 1;                          // --- anti-particle
    5616   G4LorentzVector t = qH->Get4Momentum();   // Get 4-momentum of decaying hadron
     5640  G4LorentzVector t = qH->Get4Momentum();        // Get 4-momentum of decaying hadron
    56175641  G4double m = t.m();                            // Get the mass value of decaying Hadron
    56185642  // --- Randomize a channel of decay
    56195643  G4QDecayChanVector decV = theWorld->GetQParticle(theQPDG)->GetDecayVector();
    56205644  G4int nChan = decV.size();
    5621 #ifdef pdebug
     5645#ifdef debug
    56225646  G4cout<<"G4Quasm::DecQHadron: PDG="<<thePDG<<",m="<<m<<",("<<nChan<<" channels)"<<G4endl;
    56235647#endif
    … …  
    56315655      {
    56325656        G4QDecayChan* dC = decV[i];              // The i-th Decay Channel
    5633 #ifdef pdebug
     5657#ifdef debug
    56345658        G4cout<<"G4Quasmon::DecaQHadr:i="<<i<<",r="<<rnd<<"<dl="<<dC->GetDecayChanLimit()
    56355659              <<", mm="<<dC->GetMinMass()<<G4endl;
    … …  
    56415665    G4QPDGCodeVector cV=decV[i]->GetVecOfSecHadrons();// PDGVector of theSelectedDecChannel
    56425666    G4int nPart=cV.size();                         // A#of particles to decay in
    5643 #ifdef pdebug
     5667#ifdef debug
    56445668    G4cout<<"G4Quasmon::DecayQHadron: resi="<<i<<",nP="<<nPart<<":"<<cV[0]->GetPDGCode()
    56455669          <<","<<cV[1]->GetPDGCode();
    … …  
    56535677      return theFragments;
    56545678    }
    5655 #ifdef pdebug
     5679#ifdef debug
    56565680    G4cout<<"G4Q::DecQH:Decay("<<ElMaDecays<<") PDG="<<thePDG<<t<<m<<",nP="<<nPart<<G4endl;
    56575681#endif
    … …  
    56645688      // Radiative decays In2 (eta, eta', Sigma0) are closed if the ElMaDecays=false
    56655689      if ( (fPDG != 22 && sPDG != 22) || ElMaDecays) {
    5666 #ifdef pdebug
     5690#ifdef debug
    56675691        G4cout<<"G4Q::DecQH:Yes2,fPDG="<<fPDG<<",sPDG="<<sPDG<<",EMF="<<ElMaDecays<<G4endl;
    56685692#endif
    … …  
    56915715          sHadr = new G4QHadron(sPart,sdm);               // the 2-nd Hadron is initialized
    56925716          if(sPDG<0) sHadr->MakeAntiHadron();
    5693 #ifdef pdebug
     5717#ifdef debug
    56945718          G4cout<<"G4Q::DQH:M="<<m<<",mi="<<mim<<",fd="<<fdm<<",fm="<<fm<<",sd="<<sdm
    56955719                <<",sm="<<sHadr->GetMass()<<G4endl;
    56965720#endif
    56975721        }
    5698 #ifdef pdebug
     5722#ifdef debug
    56995723        G4cout<<"G4Q::DQH:(DecayIn2)1="<<fHadr->GetMass()<<",2="<<sHadr->GetMass()<<G4endl;
    57005724#endif
    … …  
    57105734          delete fHadr;                                   // Delete "new fHadr"
    57115735          delete sHadr;                                   // Delete "new sHadr"
    5712 #ifdef pdebug
     5736#ifdef debug
    57135737          G4cerr<<"---Warning---G4Q::DecayQHadron:in2,PDGC="<<thePDG<<", ch#"<<i<<": 4M="
    57145738                <<qH->Get4Momentum()<<"("<<qH->GetMass()<<")->"<<f4Mom<<"+"<<s4Mom<<G4endl;
    … …  
    57285752          G4QHadronVector* theTmpQHV=DecayQHadron(fHadr); // Try to decay
    57295753          G4int nProd=theTmpQHV->size();
    5730 #ifdef pdebug
     5754#ifdef debug
    57315755          G4cout<<"G4Q::DecayQHadr:(DecayIn2) nOfProdForQH1="<<nProd<<G4endl;
    57325756#endif
    … …  
    57425766              copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
    57435767            }
    5744 #ifdef pdebug
     5768#ifdef debug
    57455769            G4cout<<"G4Q::DecayQHadr:(DecayIn2) Copy Sec11 nProd="<<tmpS<<G4endl;
    57465770#endif
    … …  
    57535777          theTmpQHV=DecayQHadron(sHadr);          // Try to decay
    57545778          nProd=theTmpQHV->size();
    5755 #ifdef pdebug
     5779#ifdef debug
    57565780          G4cout<<"G4Q::DecayQHadr:(DecayIn2) nOfProdForQH2="<<nProd<<G4endl;
    57575781#endif
    … …  
    57675791              copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
    57685792            }
    5769 #ifdef pdebug
     5793#ifdef debug
    57705794            G4cout<<"G4Q::DecayQHadr:(DecayIn2) Copy Sec12 nProd="<<tmpS<<G4endl;
    57715795#endif
    … …  
    57765800          delete theTmpQHV;
    57775801        }
    5778 #ifdef pdebug
     5802#ifdef debug
    57795803        G4cout<<"G4Q::DecQHadr: DecayIn2 is made with nH="<<theFragments->size()<<G4endl;
    57805804#endif
    … …  
    57825806      else
    57835807      {
    5784 #ifdef pdebug
     5808#ifdef debug
    57855809        if(thePDG==89999003||thePDG==90002999)G4cerr<<"*G4Q::DQH:8999003/90002999"<<G4endl;
    57865810#endif
    … …  
    57995823      if ( (fPDG != 22 && sPDG != 22 && tPDG != 22) || ElMaDecays)
    58005824      {
    5801 #ifdef pdebug
     5825#ifdef debug
    58025826        G4cout<<"G4Q::DQH:Y,f="<<fPDG<<",s="<<sPDG<<",t="<<tPDG<<",F="<<ElMaDecays<<G4endl;
    58035827#endif
    … …  
    58495873          if(tPDG<0) tHadr->MakeAntiHadron();
    58505874        }     
    5851 #ifdef pdebug
     5875#ifdef debug
    58525876        G4cout<<"G4Quasmon::DecayQHadron:3Dec. m1="<<fHadr->GetMass()
    58535877              <<",m2="<<sHadr->GetMass()<<",m3="<<tHadr->GetMass()<<G4endl;
    … …  
    58815905          G4QHadronVector* theTmpQHV=DecayQHadron(fHadr); // Try to decay
    58825906          G4int nProd=theTmpQHV->size();
    5883 #ifdef pdebug
     5907#ifdef debug
    58845908          G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH1="<<nProd<<G4endl;
    58855909#endif
    … …  
    58955919              copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
    58965920            }
    5897 #ifdef pdebug
     5921#ifdef debug
    58985922            G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec11 nProd="<<tmpS<<G4endl;
    58995923#endif
    … …  
    59075931          theTmpQHV=DecayQHadron(sHadr);          // Try to decay
    59085932          nProd=theTmpQHV->size();
    5909 #ifdef pdebug
     5933#ifdef debug
    59105934          G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH2="<<nProd<<G4endl;
    59115935#endif
    … …  
    59215945              copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
    59225946            }
    5923 #ifdef pdebug
     5947#ifdef debug
    59245948            G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec12 nProd="<<tmpS<<G4endl;
    59255949#endif
    … …  
    59335957          theTmpQHV=DecayQHadron(tHadr);          // Try to decay
    59345958          nProd=theTmpQHV->size();
    5935 #ifdef pdebug
     5959#ifdef debug
    59365960          G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH3="<<nProd<<G4endl;
    59375961#endif
    … …  
    59475971              copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
    59485972            }
    5949 #ifdef pdebug
     5973#ifdef debug
    59505974            G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec13 nProd="<<tmpS<<G4endl;
    59515975#endif
    … …  
    59575981
    59585982        }
    5959 #ifdef pdebug
     5983#ifdef debug
    59605984        G4cout<<"G4Q::DecQHadr: DecayIn3 is made with nH="<<theFragments->size()<<G4endl;
    59615985#endif
    … …  
    59665990  else
    59675991  {
    5968 #ifdef pdebug
     5992#ifdef debug
    59695993    G4cout<<"G4Quas::DecQHadr:Fill PDG= "<<thePDG<<t<<m<<" as it is ***0***>>>>"<<G4endl;
    59705994#endif
    … …  
    59725996    theFragments->push_back(qH);                   // Fill as it is (delete equivalent)
    59735997  }
    5974 #ifdef pdebug
     5998#ifdef debug
    59755999  G4cout<<"G4Q::DecQHadr:====HADRON IS DECAYED==== with nH="<<theFragments->size()<<G4endl;
    59766000#endif
    … …  
    60096033G4QHadronVector* G4Quasmon::Fragment(G4QNucleus& nucEnviron, G4int nQ)
    60106034{//              =====================================================
    6011 #ifdef pdebug
     6035#ifdef debug
    60126036  G4cout<<"G4Quasmon::Fragment called E="<<nucEnviron<<nucEnviron.GetProbability()<<G4endl;
    60136037#endif
    … …  
    60156039  HadronizeQuasmon(nucEnviron,nQs);
    60166040  G4int nHadrs=theQHadrons.size();
    6017 #ifdef pdebug
     6041#ifdef debug
    60186042  G4cout<<"G4Quasmon::Fragment: after HadronizeQuasmon nH="<<nHadrs<<G4endl;
    60196043#endif
    … …  
    60246048    theFragments->push_back(curHadr);         // (delete equivalent - user)
    60256049  }
    6026 #ifdef ppdebug
     6050#ifdef pdebug
    60276051  else G4cerr<<"*******G4Quasmon::Fragment *** Nothing is in the output ***"<<G4endl;
    60286052#endif
Note: See TracChangeset for help on using the changeset viewer.