//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: G4QHadron.cc,v 1.64 2009/09/02 15:45:19 mkossov Exp $
// GEANT4 tag $Name: hadr-chips-V09-03-08 $
//
//      ---------------- G4QHadron ----------------
//             by Mikhail Kossov, Sept 1999.
//      class for Quasmon initiated Hadrons generated by CHIPS Model
// -------------------------------------------------------------------
// Short description: In CHIPS all particles are G4QHadrons, while they
// can be leptons, gammas or nuclei. The G4QPDGCode makes the difference.
// In addition the 4-momentum is a basic value, so the mass can be
// different from the GS mass (e.g. for the virtual gamma).
// -------------------------------------------------------------------
//
//#define debug
//#define edebug
//#define pdebug
//#define sdebug
//#define ppdebug

#include "G4QHadron.hh"
#include <cmath>
using namespace std;

G4double G4QHadron::StrangeSuppress = 0.48;         // ? M.K.
G4double G4QHadron::sigmaPt = 1.7*GeV;              // Can be 0 ?
G4double G4QHadron::widthOfPtSquare = 0.01*GeV*GeV; // ? M.K.

G4QHadron::G4QHadron(): theMomentum(0.,0.,0.,0.), theQPDG(0), valQ(0,0,0,0,0,0), nFragm(0),
  thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
  Color(), AntiColor(), bindE(0.), formTime(0.) {}

G4QHadron::G4QHadron(G4LorentzVector p): theMomentum(p), theQPDG(0), valQ(0,0,0,0,0,0),
  nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
  Color(), AntiColor(), bindE(0.), formTime(0.) {}

// For Chipolino or Quasmon doesn't make any sense
G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p): theMomentum(p), theQPDG(PDGCode),
  nFragm(0),thePosition(0.,0.,0.),theCollisionCount(0),isSplit(false),Direction(true),
  Color(), AntiColor(), bindE(0.), formTime(0.)
{
#ifdef debug
  G4cout<<"G4QHadron must be created with PDG="<<PDGCode<<", 4M="<<p<<G4endl;
#endif
  if(GetQCode()>-1)
  {
    if(theMomentum.e()==0.) theMomentum.setE(theQPDG.GetMass());
    valQ=theQPDG.GetQuarkContent();
  }
  else if(PDGCode>80000000) DefineQC(PDGCode);
  else G4cerr<<"***G4QHadron:(P) PDG="<<PDGCode<<", use other constructor"<<G4endl;
#ifdef debug
  G4cout<<"G4QHadron is created with QCode="<<GetQCode()<<", QC="<<valQ<<G4endl;
#endif
}

// For Chipolino or Quasmon doesn't make any sense
G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p): theMomentum(p), theQPDG(QPDG),
  nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
  Color(), AntiColor(), bindE(0.), formTime(0.)
{
  if(theQPDG.GetQCode()>-1)
  {
    if(theMomentum.e()==0.) theMomentum.setE(theQPDG.GetMass());
    valQ=theQPDG.GetQuarkContent();
  }
  else
  {
    G4int cPDG=theQPDG.GetPDGCode();
    if(cPDG>80000000) DefineQC(cPDG);
    else G4cerr<<"***G4QHadr:(QP) PDG="<<cPDG<<" use other constructor"<<G4endl;
  }
}

// Make sense Chipolino or Quasmon
G4QHadron::G4QHadron(G4QContent QC, G4LorentzVector p): theMomentum(p),theQPDG(0),valQ(QC),
  nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
  Color(), AntiColor(), bindE(0.), formTime(0.)
{
  G4int curPDG=valQ.GetSPDGCode();
  if(curPDG==10&&valQ.GetBaryonNumber()>0) curPDG=valQ.GetZNSPDGCode();
  if(curPDG&&curPDG!=10) theQPDG.SetPDGCode(curPDG);
  else theQPDG.InitByQCont(QC);
}

G4QHadron::G4QHadron(G4int PDGCode, G4double aMass, G4QContent QC) :
  theMomentum(0.,0.,0.,aMass), theQPDG(PDGCode), valQ(QC), nFragm(0),thePosition(0.,0.,0.),
  theCollisionCount(0), isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.),
  formTime(0.)
{}

G4QHadron::G4QHadron(G4QPDGCode QPDG, G4double aMass, G4QContent QC) :
  theMomentum(0.,0.,0.,aMass), theQPDG(QPDG), valQ(QC), nFragm(0), thePosition(0.,0.,0.),
  theCollisionCount(0), isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.),
  formTime(0.)
{}

G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p, G4QContent QC) : theMomentum(p),
  theQPDG(PDGCode), valQ(QC), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
  isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
{}

G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p, G4QContent QC) : theMomentum(p),
  theQPDG(QPDG), valQ(QC), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
  isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
{}

G4QHadron::G4QHadron(G4QParticle* pPart, G4double maxM) : theMomentum(0.,0.,0.,0.),
  theQPDG(pPart->GetQPDG()), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
  isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
{
#ifdef debug
  G4cout<<"G4QHadron is created & randomized with maxM="<<maxM<<G4endl;
#endif
  G4int PDGCode = theQPDG.GetPDGCode();
  if(PDGCode<2)G4cerr<<"***G4QHadron:(M) PDGC="<<PDGCode<<" use other constructor"<<G4endl;
  valQ=theQPDG.GetQuarkContent();
  theMomentum.setE(RandomizeMass(pPart, maxM));
}

G4QHadron::G4QHadron(const G4QHadron& right)
{
  theMomentum         = right.theMomentum;
  theQPDG             = right.theQPDG;
  valQ                = right.valQ;
  nFragm              = right.nFragm;
  thePosition         = right.thePosition;      
  theCollisionCount   = 0;
  isSplit             = false;
  Direction           = right.Direction;
  bindE               = right.bindE;
  formTime            = right.formTime;
}

G4QHadron::G4QHadron(const G4QHadron* right)
{
  theMomentum         = right->theMomentum;
  theQPDG             = right->theQPDG;
  valQ                = right->valQ;
  nFragm              = right->nFragm;
  thePosition         = right->thePosition;      
  theCollisionCount   = 0;
  isSplit             = false;
  Direction           = right->Direction;
  bindE               = right->bindE;
  formTime            = right->formTime;
}

G4QHadron::G4QHadron(const G4QHadron* right, G4int C, G4ThreeVector P, G4LorentzVector M)
{
  theMomentum         = M;
  theQPDG             = right->theQPDG;
  valQ                = right->valQ;
  nFragm              = right->nFragm;
  thePosition         = P;      
  theCollisionCount   = C;
  isSplit             = false;
  Direction           = right->Direction;
  bindE               = right->bindE;
  formTime            = right->formTime;
}

const G4QHadron& G4QHadron::operator=(const G4QHadron &right)
{
  if(this != &right)                          // Beware of self assignment
  {
    theMomentum         = right.theMomentum;
    theQPDG             = right.theQPDG;
    valQ                = right.valQ;
    nFragm              = right.nFragm;
    thePosition         = right.thePosition;      
    theCollisionCount   = 0;
    isSplit             = false;
    Direction           = right.Direction;
    bindE               = right.bindE;
  }
  return *this;
}

G4QHadron::~G4QHadron()
{
  std::list<G4QParton*>::iterator ipos = Color.begin();
  std::list<G4QParton*>::iterator epos = Color.end();
  for( ; ipos != epos; ipos++) {delete [] *ipos;}
  Color.clear();

  ipos = AntiColor.begin();
  epos = AntiColor.end();
  for( ; ipos != epos; ipos++) {delete [] *ipos;}
  AntiColor.clear();
}

// Define quark content of the particle with a particular PDG Code
void G4QHadron::DefineQC(G4int PDGCode)
{
  //G4cout<<"G4QHadron::DefineQC is called with PDGCode="<<PDGCode<<G4endl;
  G4int szn=PDGCode-90000000;
  G4int ds=0;
  G4int dz=0;
  G4int dn=0;
  if(szn<-100000)
  {
    G4int ns=(-szn)/1000000+1;
    szn+=ns*1000000;
    ds+=ns;
  }
  else if(szn<-100)
  {
    G4int nz=(-szn)/1000+1;
    szn+=nz*1000;
    dz+=nz;
  }
  else if(szn<0)
  {
    G4int nn=-szn;
    szn=0;
    dn+=nn;
  }
  G4int sz =szn/1000;
  G4int n  =szn%1000;
  if(n>700)
  {
    n-=1000;
    dz--;
  }
  G4int z  =sz%1000-dz;
  if(z>700)
  {
    z-=1000;
    ds--;
  }
  G4int Sq =sz/1000-ds;
  G4int zns=z+n+Sq;
  G4int Dq=n+zns;
  G4int Uq=z+zns;
  if      (Dq<0&&Uq<0&&Sq<0)valQ=G4QContent(0 ,0 ,0 ,-Dq,-Uq,-Sq);
  else if (Uq<0&&Sq<0)      valQ=G4QContent(Dq,0 ,0 ,0  ,-Uq,-Sq);
  else if (Dq<0&&Sq<0)      valQ=G4QContent(0 ,Uq,0 ,-Dq,0  ,-Sq);
  else if (Dq<0&&Uq<0)      valQ=G4QContent(0 ,0 ,Sq,-Dq,-Uq,0  );
  else if (Uq<0)            valQ=G4QContent(Dq,0 ,Sq,0  ,-Uq,0  );
  else if (Sq<0)            valQ=G4QContent(Dq,Uq,0 ,0  ,0  ,-Sq);
  else if (Dq<0)            valQ=G4QContent(0 ,Uq,Sq,-Dq,0  ,0  );
  else                      valQ=G4QContent(Dq,Uq,Sq,0  ,0  ,0  );
}

// Redefine a Hadron with a new PDGCode
void G4QHadron::SetQPDG(const G4QPDGCode& newQPDG)
{
  theQPDG  = newQPDG;
  G4int PDG= newQPDG.GetPDGCode();
  G4int Q  = newQPDG.GetQCode();
#ifdef debug
  G4cout<<"G4QHadron::SetQPDG is called with PDGCode="<<PDG<<", QCode="<<Q<<G4endl;
#endif
  if     (Q>-1) valQ=theQPDG.GetQuarkContent();
  else if(PDG>80000000) DefineQC(PDG);
  else
  {
    G4cerr<<"***G4QHadron::SetQPDG: QPDG="<<newQPDG<<G4endl;
    throw G4QException("***G4QHadron::SetQPDG: Impossible QPDG Probably a Chipolino");
  }
}

// Decay of Hadron In2Particles f&s, f is in respect to the direction of HadronMomentumDir
G4bool G4QHadron::RelDecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
       G4LorentzVector& dir, G4double maxCost, G4double minCost)
{//    ===================================================================
  G4double fM2 = f4Mom.m2();
  G4double fM  = sqrt(fM2);              // Mass of the 1st Hadron
  G4double sM2 = s4Mom.m2();
  G4double sM  = sqrt(sM2);              // Mass of the 2nd Hadron
  G4double iM2 = theMomentum.m2();
  G4double iM  = sqrt(iM2);              // Mass of the decaying hadron
  G4double vP  = theMomentum.rho();      // Momentum of the decaying hadron
  G4double dE  = theMomentum.e();        // Energy of the decaying hadron
  if(dE<vP)
  {
    G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
    G4double accuracy=.000001*vP;
    G4double emodif=std::fabs(dE-vP);
    //if(emodif<accuracy)
    //{
      G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
      theMomentum.setE(vP+emodif+.01*accuracy);
    //}
  }
  G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans.
  G4ThreeVector ltf = -ltb;              // Boost vector for forward Lorentz Trans.
  G4LorentzVector cdir = dir;            // A copy to make a transformation to CMS
#ifdef ppdebug
  if(cdir.e()+.001<cdir.rho()) G4cerr<<"*G4QH::RDIn2:*Boost* cd4M="<<cdir<<",e-p="
                                     <<cdir.e()-cdir.rho()<<G4endl;
#endif
  cdir.boost(ltf);                       // Direction transpormed to CMS of the Momentum
  G4ThreeVector vdir = cdir.vect();      // 3-Vector of the direction-particle
#ifdef ppdebug
  G4cout<<"G4QHad::RelDI2:dir="<<dir<<",ltf="<<ltf<<",cdir="<<cdir<<",vdir="<<vdir<<G4endl;
#endif
  G4ThreeVector vx(0.,0.,1.);            // Ort in the direction of the reference particle
  G4ThreeVector vy(0.,1.,0.);            // First ort orthogonal to the direction
  G4ThreeVector vz(1.,0.,0.);            // Second ort orthoganal to the direction
  if(vdir.mag2() > 0.)                   // the refference particle isn't at rest in CMS
  {
    vx = vdir.unit();                    // Ort in the direction of the reference particle
#ifdef ppdebug
    G4cout<<"G4QH::RelDecIn2:Vx="<<vx<<",M="<<theMomentum<<",d="<<dir<<",c="<<cdir<<G4endl;
#endif
    G4ThreeVector vv= vx.orthogonal();   // Not normed orthogonal vector (!)
    vy = vv.unit();                      // First ort orthogonal to the direction
    vz = vx.cross(vy);                   // Second ort orthoganal to the direction
  }
#ifdef ppdebug
  G4cout<<"G4QHad::RelDecIn2:iM="<<iM<<"=>fM="<<fM<<"+sM="<<sM<<",ob="<<vx<<vy<<vz<<G4endl;
#endif
  if(maxCost> 1.) maxCost= 1.;
  if(minCost<-1.) minCost=-1.;
  if(maxCost<-1.) maxCost=-1.;
  if(minCost> 1.) minCost= 1.;
  if(minCost> maxCost) minCost=maxCost;
  if(fabs(iM-fM-sM)<.00000001)
  {
    G4double fR=fM/iM;
    G4double sR=sM/iM;
    f4Mom=fR*theMomentum;
    s4Mom=sR*theMomentum;
    return true;
  }
  else if (iM+.001<fM+sM || iM==0.)
  {//@@ Later on make a quark content check for the decay
    G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
    return false;
  }
  G4double d2 = iM2-fM2-sM2;
  G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;    // Decay momentum(^2) in CMS of Quasmon
  if(p2<0.)
  {
#ifdef ppdebug
    G4cout<<"**G4QH:RDIn2:p2="<<p2<<"<0,d2^2="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
#endif
    p2=0.;
  }
  G4double p  = sqrt(p2);
  G4double ct = maxCost;
  if(maxCost>minCost)
  {
    G4double dcost=maxCost-minCost;
    ct = minCost+dcost*G4UniformRand();
  }
  G4double phi= twopi*G4UniformRand();  // @@ Change 360.*deg to M_TWOPI (?)
  G4double ps=0.;
  if(fabs(ct)<1.) ps = p * sqrt(1.-ct*ct);
  else
  {
#ifdef ppdebug
    G4cout<<"**G4QH::RDIn2:ct="<<ct<<",mac="<<maxCost<<",mic="<<minCost<<G4endl;
    //throw G4QException("***G4QHadron::RDIn2: bad cos(theta)");
#endif
    if(ct>1.) ct=1.;
    if(ct<-1.) ct=-1.;
  }
  G4ThreeVector pVect=(ps*sin(phi))*vz+(ps*cos(phi))*vy+p*ct*vx;
#ifdef ppdebug
  G4cout<<"G4QH::RelDIn2:ct="<<ct<<",p="<<p<<",ps="<<ps<<",ph="<<phi<<",v="<<pVect<<G4endl;
#endif

  f4Mom.setVect(pVect);
  f4Mom.setE(sqrt(fM2+p2));
  s4Mom.setVect((-1)*pVect);
  s4Mom.setE(sqrt(sM2+p2));
  
#ifdef ppdebug
  G4cout<<"G4QHadr::RelDecIn2:p2="<<p2<<",v="<<ltb<<",f4M="<<f4Mom<<" + s4M="<<s4Mom<<" = "
        <<f4Mom+s4Mom<<", M="<<iM<<G4endl;
#endif
  if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
                                      <<f4Mom.e()-f4Mom.rho()<<G4endl;
  f4Mom.boost(ltb);                        // Lor.Trans. of 1st hadron back to LS
  if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
                                      <<s4Mom.e()-s4Mom.rho()<<G4endl;
  s4Mom.boost(ltb);                        // Lor.Trans. of 2nd hadron back to LS
#ifdef ppdebug
  G4cout<<"G4QHadron::RelDecayIn2:Output, f4Mom="<<f4Mom<<" + s4Mom="<<s4Mom<<" = "
        <<f4Mom+s4Mom<<", d4M="<<theMomentum-f4Mom-s4Mom<<G4endl;
#endif
  return true;
} // End of "RelDecayIn2"

// Decay of Hadron In2Particles f&s, f w/r/to dN/dO [cp>0: ~cost^cp, cp<0: ~(1-cost)^(-cp)]
G4bool G4QHadron::CopDecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
                              G4LorentzVector& dir, G4double cosp)
{//    ===================================================================
  G4double fM2 = f4Mom.m2();
  G4double fM  = sqrt(fM2);              // Mass of the 1st Hadron
  G4double sM2 = s4Mom.m2();
  G4double sM  = sqrt(sM2);              // Mass of the 2nd Hadron
  G4double iM2 = theMomentum.m2();
  G4double iM  = sqrt(iM2);              // Mass of the decaying hadron
  G4double vP  = theMomentum.rho();      // Momentum of the decaying hadron
  G4double dE  = theMomentum.e();        // Energy of the decaying hadron
  G4bool neg=false;                // Negative (backward) distribution of t
  if(cosp<0)
  {
    cosp=-cosp;
    neg=true;
  }
  if(dE<vP)
  {
    G4cerr<<"***G4QHad::CopDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
    G4double accuracy=.000001*vP;
    G4double emodif=std::fabs(dE-vP);
    //if(emodif<accuracy)
    //{
      G4cerr<<"G4QHadron::CopDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
      theMomentum.setE(vP+emodif+.01*accuracy);
    //}
  }
  G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans.
  G4ThreeVector ltf = -ltb;              // Boost vector for forward Lorentz Trans.
  G4LorentzVector cdir = dir;            // A copy to make a transformation to CMS
#ifdef ppdebug
  if(cdir.e()+.001<cdir.rho()) G4cerr<<"*G4QH::RDIn2:*Boost* cd4M="<<cdir<<",e-p="
                                     <<cdir.e()-cdir.rho()<<G4endl;
#endif
  cdir.boost(ltf);                       // Direction transpormed to CMS of the Momentum
  G4ThreeVector vdir = cdir.vect();      // 3-Vector of the direction-particle
#ifdef ppdebug
  G4cout<<"G4QHad::CopDI2:dir="<<dir<<",ltf="<<ltf<<",cdir="<<cdir<<",vdir="<<vdir<<G4endl;
#endif
  G4ThreeVector vx(0.,0.,1.);            // Ort in the direction of the reference particle
  G4ThreeVector vy(0.,1.,0.);            // First ort orthogonal to the direction
  G4ThreeVector vz(1.,0.,0.);            // Second ort orthoganal to the direction
  if(vdir.mag2() > 0.)                   // the refference particle isn't at rest in CMS
  {
    vx = vdir.unit();                    // Ort in the direction of the reference particle
#ifdef ppdebug
    G4cout<<"G4QH::CopDecIn2:Vx="<<vx<<",M="<<theMomentum<<",d="<<dir<<",c="<<cdir<<G4endl;
#endif
    G4ThreeVector vv= vx.orthogonal();   // Not normed orthogonal vector (!)
    vy = vv.unit();                      // First ort orthogonal to the direction
    vz = vx.cross(vy);                   // Second ort orthoganal to the direction
  }
#ifdef ppdebug
  G4cout<<"G4QHad::CopDecIn2:iM="<<iM<<"=>fM="<<fM<<"+sM="<<sM<<",ob="<<vx<<vy<<vz<<G4endl;
#endif
  if(fabs(iM-fM-sM)<.00000001)
  {
    G4double fR=fM/iM;
    G4double sR=sM/iM;
    f4Mom=fR*theMomentum;
    s4Mom=sR*theMomentum;
    return true;
  }
  else if (iM+.001<fM+sM || iM==0.)
  {//@@ Later on make a quark content check for the decay
    G4cerr<<"***G4QH::CopDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
    return false;
  }
  G4double d2 = iM2-fM2-sM2;
  G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;    // Decay momentum(^2) in CMS of Quasmon
  if(p2<0.)
  {
#ifdef ppdebug
    G4cout<<"*G4QH:CopDI2:p2="<<p2<<"<0,d4/4="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
#endif
    p2=0.;
  }
  G4double p  = sqrt(p2);
  G4double ct = 0;
  G4double rn = pow(G4UniformRand(),cosp+1.);
  if(neg)  ct = rn+rn-1.;                  // More backward than forward
  else     ct = 1.-rn-rn;                  // More forward than backward
  //
  G4double phi= twopi*G4UniformRand();  // @@ Change 360.*deg to M_TWOPI (?)
  G4double ps=0.;
  if(fabs(ct)<1.) ps = p * sqrt(1.-ct*ct);
  else
  {
#ifdef ppdebug
    G4cout<<"**G4QH::CopDecayIn2:ct="<<ct<<",mac="<<maxCost<<",mic="<<minCost<<G4endl;
    //throw G4QException("***G4QHadron::RDIn2: bad cos(theta)");
#endif
    if(ct>1.) ct=1.;
    if(ct<-1.) ct=-1.;
  }
  G4ThreeVector pVect=(ps*sin(phi))*vz+(ps*cos(phi))*vy+p*ct*vx;
#ifdef ppdebug
  G4cout<<"G4QH::CopDIn2:ct="<<ct<<",p="<<p<<",ps="<<ps<<",ph="<<phi<<",v="<<pVect<<G4endl;
#endif

  f4Mom.setVect(pVect);
  f4Mom.setE(sqrt(fM2+p2));
  s4Mom.setVect((-1)*pVect);
  s4Mom.setE(sqrt(sM2+p2));
  
#ifdef ppdebug
  G4cout<<"G4QHadr::CopDecIn2:p2="<<p2<<",v="<<ltb<<",f4M="<<f4Mom<<" + s4M="<<s4Mom<<" = "
        <<f4Mom+s4Mom<<", M="<<iM<<G4endl;
#endif
  if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
                                      <<f4Mom.e()-f4Mom.rho()<<G4endl;
  f4Mom.boost(ltb);                        // Lor.Trans. of 1st hadron back to LS
  if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
                                      <<s4Mom.e()-s4Mom.rho()<<G4endl;
  s4Mom.boost(ltb);                        // Lor.Trans. of 2nd hadron back to LS
#ifdef ppdebug
  G4cout<<"G4QHadron::CopDecayIn2:Output, f4Mom="<<f4Mom<<" + s4Mom="<<s4Mom<<" = "
        <<f4Mom+s4Mom<<", d4M="<<theMomentum-f4Mom-s4Mom<<G4endl;
#endif
  return true;
} // End of "CopDecayIn2"

// Decay of the Hadron in 2 particles (f + s)
G4bool G4QHadron::DecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom)
{//    ===================================================================
  G4double fM2 = f4Mom.m2();
  if(fM2<0.) fM2=0.;
  G4double fM  = sqrt(fM2);              // Mass of the 1st Hadron
  G4double sM2 = s4Mom.m2();
  if(sM2<0.) sM2=0.;
  G4double sM  = sqrt(sM2);              // Mass of the 2nd Hadron
  G4double iM2 = theMomentum.m2();
  if(iM2<0.) iM2=0.;
  G4double iM  = sqrt(iM2);              // Mass of the decaying hadron
#ifdef debug
  G4cout<<"G4QHadron::DecIn2: iM="<<iM<<" => fM="<<fM<<" + sM="<<sM<<" = "<<fM+sM<<G4endl;
#endif
  //@@ Later on make a quark content check for the decay
  if (fabs(iM-fM-sM)<.0000001)
  {
    G4double fR=fM/iM;
    G4double sR=sM/iM;
    f4Mom=fR*theMomentum;
    s4Mom=sR*theMomentum;
    return true;
  }
  else if (iM+.001<fM+sM || iM==0.)
  {
#ifdef debug
    G4cerr<<"***G4QHadron::DecayIn2*** fM="<<fM<<" + sM="<<sM<<"="<<fM+sM<<" > iM="<<iM
          <<", d="<<iM-fM-sM<<G4endl;
#endif
    return false;
  }

  G4double d2 = iM2-fM2-sM2;
  G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;    // Decay momentum(^2) in CMS of Quasmon
  if (p2<0.)
  {
#ifdef debug
    G4cerr<<"***G4QH::DI2:p2="<<p2<<"<0,d2^2="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
#endif
    p2=0.;
  }
  G4double p  = sqrt(p2);
  G4double ct = 1.-2*G4UniformRand();
#ifdef debug
  G4cout<<"G4QHadron::DecayIn2: ct="<<ct<<", p="<<p<<G4endl;
#endif
  G4double phi= twopi*G4UniformRand();  // @@ Change 360.*deg to M_TWOPI (?)
  G4double ps = p * sqrt(1.-ct*ct);
  G4ThreeVector pVect(ps*sin(phi),ps*cos(phi),p*ct);

  f4Mom.setVect(pVect);
  f4Mom.setE(sqrt(fM2+p2));
  s4Mom.setVect((-1)*pVect);
  s4Mom.setE(sqrt(sM2+p2));

  if(theMomentum.e()<theMomentum.rho())
  {
    G4cerr<<"*G4QH::DecIn2:*Boost* 4M="<<theMomentum<<",e-p="
          <<theMomentum.e()-theMomentum.rho()<<G4endl;
    //throw G4QException("G4QHadron::DecayIn2: Decay of particle with zero mass")
    theMomentum.setE(1.0000001*theMomentum.rho());
  }
  G4double vP  = theMomentum.rho();      // Momentum of the decaying hadron
  G4double dE  = theMomentum.e();        // Energy of the decaying hadron
  if(dE<vP)
  {
    G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
    G4double accuracy=.000001*vP;
    G4double emodif=std::fabs(dE-vP);
    if(emodif<accuracy)
    {
      G4cerr<<"G4QHadron::DecayIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
      theMomentum.setE(vP+emodif+.01*accuracy);
    }
  }
  G4ThreeVector ltb = theMomentum.boostVector(); // Boost vector for backward Lor.Trans.
#ifdef debug
  G4cout<<"G4QHadron::DecIn2:LorTrans v="<<ltb<<",f4Mom="<<f4Mom<<",s4Mom="<<s4Mom<<G4endl;
#endif
  if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* f4M="<<f4Mom<<G4endl;
  f4Mom.boost(ltb);                        // Lor.Trans. of 1st hadron back to LS
  if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* s4M="<<s4Mom<<G4endl; 
  s4Mom.boost(ltb);                        // Lor.Trans. of 2nd hadron back to LS
#ifdef debug
  G4cout<<"G4QHadron::DecayIn2: ROOT OUTPUT f4Mom="<<f4Mom<<", s4Mom="<<s4Mom<<G4endl;
#endif
  return true;
} // End of "DecayIn2"

// Correction for the Hadron + fr decay in case of the new corM mass of the Hadron
G4bool G4QHadron::CorMDecayIn2(G4double corM, G4LorentzVector& fr4Mom)
{//    ===============================================================
  G4double fM  = fr4Mom.m();                // Mass of the Fragment
  G4LorentzVector comp=theMomentum+fr4Mom;  // 4Mom of the decaying compound system
  G4double iM  = comp.m();                  // mass of the decaying compound system
#ifdef debug
  G4cout<<"G4QH::CMDIn2: iM="<<iM<<comp<<"=>fM="<<fM<<"+corM="<<corM<<"="<<fM+corM<<G4endl;
#endif
  G4double dE=iM-fM-corM;
  //@@ Later on make a quark content check for the decay
  if (fabs(dE)<.001)
  {
    G4double fR=fM/iM;
    G4double cR=corM/iM;
    fr4Mom=fR*comp;
    theMomentum=cR*comp;
    return true;
  }
  else if (dE<-.001 || iM==0.)
  {
    G4cerr<<"***G4QH::CorMDIn2***fM="<<fM<<" + cM="<<corM<<" > iM="<<iM<<",d="<<dE<<G4endl;
    return false;
  }
  G4double corM2= corM*corM;
  G4double fM2 = fM*fM;
  G4double iM2 = iM*iM;
  G4double d2 = iM2-fM2-corM2;
  G4double p2 = (d2*d2/4.-fM2*corM2)/iM2;    // Decay momentum(^2) in CMS of Quasmon
  if (p2<0.)
  {
#ifdef debug
    G4cerr<<"**G4QH::CMDI2:p2="<<p2<<"<0,d="<<d2*d2/4.<<"<4*fM2*hM2="<<4*fM2*corM2<<G4endl;
#endif
    p2=0.;
  }
  G4double p  = sqrt(p2);
  if(comp.e()<comp.rho())G4cerr<<"*G4QH::CorMDecayIn2:*Boost* comp4M="<<comp<<",e-p="
                               <<comp.e()-comp.rho()<<G4endl;
  G4ThreeVector ltb = comp.boostVector();      // Boost vector for backward Lor.Trans.
  G4ThreeVector ltf = -ltb;                    // Boost vector for forward Lorentz Trans.
  G4LorentzVector cm4Mom=fr4Mom;               // Copy of fragment 4Mom to transform to CMS
  if(cm4Mom.e()<cm4Mom.rho())
  {
    G4cerr<<"*G4QH::CorMDecIn2:*Boost* c4M="<<cm4Mom<<G4endl;
    cm4Mom.setE(1.0000001*cm4Mom.rho());
  }
  cm4Mom.boost(ltf);                           // Now it is in CMS (Forward Lor.Trans.)
  G4double pfx= cm4Mom.px();
  G4double pfy= cm4Mom.py();
  G4double pfz= cm4Mom.pz();
  G4double pt2= pfx*pfx+pfy*pfy;
  G4double tx=0.;
  G4double ty=0.;
  if(pt2<=0.)
  {
    G4double phi= 360.*deg*G4UniformRand();  // @@ Change 360.*deg to M_TWOPI (?)
    tx=sin(phi);
    ty=cos(phi);
  }
  else
  {
    G4double pt=sqrt(pt2);
    tx=pfx/pt;
    ty=pfy/pt;
  }
  G4double pc2=pt2+pfz*pfz;
  G4double ct=0.;
  if(pc2<=0.)
  {
    G4double rnd= G4UniformRand();
    ct=1.-rnd-rnd;
  }
  else
  {
    G4double pc=sqrt(pc2);
    ct=pfz/pc;
  }
#ifdef debug
  G4cout<<"G4QHadron::CorMDecayIn2: ct="<<ct<<", p="<<p<<G4endl;
#endif
  G4double ps = p * sqrt(1.-ct*ct);
  G4ThreeVector pVect(ps*tx,ps*ty,p*ct);
  fr4Mom.setVect(pVect);
  fr4Mom.setE(sqrt(fM2+p2));
  theMomentum.setVect((-1)*pVect);
  theMomentum.setE(sqrt(corM2+p2));
#ifdef debug
  G4LorentzVector dif2=comp-fr4Mom-theMomentum;
  G4cout<<"G4QH::CorMDIn2:c="<<comp<<"-f="<<fr4Mom<<"-4M="<<theMomentum<<"="<<dif2<<G4endl;
#endif
  if(fr4Mom.e()+.001<fr4Mom.rho())G4cerr<<"*G4QH::CorMDecIn2:*Boost*fr4M="<<fr4Mom<<G4endl;
  fr4Mom.boost(ltb);                        // Lor.Trans. of the Fragment back to LS
  if(theMomentum.e()<theMomentum.rho())
  {
    G4cerr<<"*G4QH::CMDI2:4="<<theMomentum<<G4endl;
    theMomentum.setE(1.0000001*theMomentum.rho());
  }
  theMomentum.boost(ltb);                  // Lor.Trans. of the Hadron back to LS
#ifdef debug
  G4LorentzVector dif3=comp-fr4Mom-theMomentum;
  G4cout<<"G4QH::CorMDecIn2:OUTPUT:f4M="<<fr4Mom<<",h4M="<<theMomentum<<"d="<<dif3<<G4endl;
#endif
  return true;
} // End of "CorMDecayIn2"


// Fragment fr4Mom louse energy corE and transfer it to This Hadron 
G4bool G4QHadron::CorEDecayIn2(G4double corE, G4LorentzVector& fr4Mom)
{//    ===============================================================
  G4double fE  = fr4Mom.m();                // Energy of the Fragment
#ifdef debug
  G4cout<<"G4QH::CorEDecIn2:fE="<<fE<<fr4Mom<<">corE="<<corE<<",h4M="<<theMomentum<<G4endl;
#endif
  if (fE+.001<=corE)
  {
#ifdef debug
    G4cerr<<"***G4QHadron::CorEDecIn2*** fE="<<fE<<"<corE="<<corE<<", d="<<corE-fE<<G4endl;
#endif
    return false;
  }
  G4double fM2=fr4Mom.m2();                 // Squared Mass of the Fragment
  if(fM2<0.) fM2=0.;
  G4double iPx=fr4Mom.px();                 // Initial Px of the Fragment
  G4double iPy=fr4Mom.py();                 // Initial Py of the Fragment
  G4double iPz=fr4Mom.pz();                 // Initial Pz of the Fragment
  G4double fP2=iPx*iPx+iPy*iPy+iPz*iPz;     // Initial Squared 3-momentum of the Fragment
  G4double finE = fE - corE;                // Final energy of the fragment
  G4double rP = sqrt((finE*finE-fM2)/fP2);  // Reduction factor for the momentum
  G4double fPx=iPx*rP;
  G4double fPy=iPy*rP;
  G4double fPz=iPz*rP;
  fr4Mom= G4LorentzVector(fPx,fPy,fPz,finE);
  G4double Px=theMomentum.px()+iPx-fPx;
  G4double Py=theMomentum.py()+iPy-fPy;
  G4double Pz=theMomentum.pz()+iPz-fPz;
  G4double mE=theMomentum.e();
  ///////////G4double mM2=theMomentum.m2();
  theMomentum= G4LorentzVector(Px,Py,Pz,mE+corE);
#ifdef debug
  G4double difF=fr4Mom.m2()-fM2;
  G4cout<<"G4QH::CorEDecIn2: dF="<<difF<<",out:"<<theMomentum<<fr4Mom<<G4endl;
#endif
  return true;
} // End of "CorEDecayIn2"

// Decay of the hadron in 3 particles i=>r+s+t
G4bool G4QHadron::DecayIn3
                   (G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, G4LorentzVector& t4Mom)
{//    ====================================================================================
#ifdef debug
  G4cout<<"G4QH::DIn3:"<<theMomentum<<"=>pf="<<f4Mom<<"+ps="<<s4Mom<<"+pt="<<t4Mom<<G4endl;
#endif
  G4double iM  = theMomentum.m();  // Mass of the decaying hadron
  G4double fM  = f4Mom.m();        // Mass of the 1st hadron
  G4double sM  = s4Mom.m();        // Mass of the 2nd hadron
  G4double tM  = t4Mom.m();        // Mass of the 3rd hadron
  G4double eps = 0.001;            // Accuracy of the split condition
  if (fabs(iM-fM-sM-tM)<=eps)
  {
    G4double fR=fM/iM;
    G4double sR=sM/iM;
    G4double tR=tM/iM;
    f4Mom=fR*theMomentum;
    s4Mom=sR*theMomentum;
    t4Mom=tR*theMomentum;
    return true;
  }
  if (iM+eps<fM+sM+tM)
  {
    G4cout<<"***G4QHadron::DecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
          <<",d="<<iM-fM-sM-tM<<G4endl;
    return false;
  }
  G4double fM2 = fM*fM;
  G4double sM2 = sM*sM;
  G4double tM2 = tM*tM;
  G4double iM2 = iM*iM;
  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
  G4double m12sMin =(fM+sM)*(fM+sM);
  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
  G4double rR = 0.;
  G4double rnd= 1.;
#ifdef debug
  G4int    tr = 0;                 //@@ Comment if "cout" below is skiped @@
#endif
  G4double m12s = 0.;              // Fake definition before the Loop
  while (rnd > rR)
  {
    m12s = m12sMin + m12sBase*G4UniformRand();
    G4double e1=m12s+fM2-sM2;
    G4double e2=iM2-m12s-tM2;
    G4double four12=4*m12s;
    G4double m13sRange=0.;
    G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
    if(dif<0.)
    {
#ifdef debug
      if(dif<-.01) G4cerr<<"*G4QHadron::DecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
                         <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
#endif
    }
    else m13sRange=sqrt(dif)/m12s;
    rR = m13sRange/m13sBase;
    rnd= G4UniformRand();
#ifdef debug
    G4cout<<"G4QHadron::DecayIn3: try to decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
#endif
  }
  G4double m12 = sqrt(m12s);       // Mass of the H1+H2 system
  G4LorentzVector dh4Mom(0.,0.,0.,m12);
  
  if(!DecayIn2(t4Mom,dh4Mom))
  {
    G4cerr<<"***G4QHadron::DecayIn3: Exception1"<<G4endl;
    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
    return false;
  }
#ifdef debug
  G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
#endif
  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
  {
    G4cerr<<"***G4QHadron::DecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
    return false;
  }
  return true;
} // End of DecayIn3

// Relative Decay of the hadron in 3 particles i=>f+s+t (t is with respect to minC<ct<maxC)
G4bool G4QHadron::RelDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
                              G4LorentzVector& t4Mom, G4LorentzVector& dir,
                              G4double maxCost, G4double minCost)
{//    ====================================================================================
#ifdef debug
  G4cout<<"G4QH::RelDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl;
#endif
  G4double iM  = theMomentum.m();  // Mass of the decaying hadron
  G4double fM  = f4Mom.m();        // Mass of the 1st hadron
  G4double sM  = s4Mom.m();        // Mass of the 2nd hadron
  G4double tM  = t4Mom.m();        // Mass of the 3rd hadron
  G4double eps = 0.001;            // Accuracy of the split condition
  if (fabs(iM-fM-sM-tM)<=eps)
  {
    G4double fR=fM/iM;
    G4double sR=sM/iM;
    G4double tR=tM/iM;
    f4Mom=fR*theMomentum;
    s4Mom=sR*theMomentum;
    t4Mom=tR*theMomentum;
    return true;
  }
  if (iM+eps<fM+sM+tM)
  {
    G4cout<<"***G4QHadron::RelDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
          <<",d="<<iM-fM-sM-tM<<G4endl;
    return false;
  }
  G4double fM2 = fM*fM;
  G4double sM2 = sM*sM;
  G4double tM2 = tM*tM;
  G4double iM2 = iM*iM;
  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
  G4double m12sMin =(fM+sM)*(fM+sM);
  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
  G4double rR = 0.;
  G4double rnd= 1.;
#ifdef debug
  G4int    tr = 0;                 //@@ Comment if "cout" below is skiped @@
#endif
  G4double m12s = 0.;              // Fake definition before the Loop
  while (rnd > rR)
  {
    m12s = m12sMin + m12sBase*G4UniformRand();
    G4double e1=m12s+fM2-sM2;
    G4double e2=iM2-m12s-tM2;
    G4double four12=4*m12s;
    G4double m13sRange=0.;
    G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
    if(dif<0.)
    {
#ifdef debug
      if(dif<-.01) G4cerr<<"G4QHadron::RelDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
                         <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
#endif
    }
    else m13sRange=sqrt(dif)/m12s;
    rR = m13sRange/m13sBase;
    rnd= G4UniformRand();
#ifdef debug
    G4cout<<"G4QHadron::RelDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
#endif
  }
  G4double m12 = sqrt(m12s);       // Mass of the H1+H2 system
  G4LorentzVector dh4Mom(0.,0.,0.,m12);
  
  if(!RelDecayIn2(t4Mom,dh4Mom,dir,maxCost,minCost))
  {
    G4cerr<<"***G4QHadron::RelDecayIn3: Exception1"<<G4endl;
    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
    return false;
  }
#ifdef debug
  G4cout<<"G4QHadron::RelDecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
#endif
  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
  {
    G4cerr<<"***G4QHadron::RelDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
    return false;
  }
  return true;
} // End of RelDecayIn3

// Relative Decay of hadron in 3: i=>f+s+t.  dN/dO [cp>0:~cost^cp, cp<0:~(1-cost)^(-cp)]
G4bool G4QHadron::CopDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
                              G4LorentzVector& t4Mom, G4LorentzVector& dir, G4double cosp)
{//    ====================================================================================
#ifdef debug
  G4cout<<"G4QH::CopDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl;
#endif
  G4double iM  = theMomentum.m();  // Mass of the decaying hadron
  G4double fM  = f4Mom.m();        // Mass of the 1st hadron
  G4double sM  = s4Mom.m();        // Mass of the 2nd hadron
  G4double tM  = t4Mom.m();        // Mass of the 3rd hadron
  G4double eps = 0.001;            // Accuracy of the split condition
  if (fabs(iM-fM-sM-tM)<=eps)
  {
    G4double fR=fM/iM;
    G4double sR=sM/iM;
    G4double tR=tM/iM;
    f4Mom=fR*theMomentum;
    s4Mom=sR*theMomentum;
    t4Mom=tR*theMomentum;
    return true;
  }
  if (iM+eps<fM+sM+tM)
  {
    G4cout<<"***G4QHadron::CopDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
          <<",d="<<iM-fM-sM-tM<<G4endl;
    return false;
  }
  G4double fM2 = fM*fM;
  G4double sM2 = sM*sM;
  G4double tM2 = tM*tM;
  G4double iM2 = iM*iM;
  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
  G4double m12sMin =(fM+sM)*(fM+sM);
  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
  G4double rR = 0.;
  G4double rnd= 1.;
#ifdef debug
  G4int    tr = 0;                 //@@ Comment if "cout" below is skiped @@
#endif
  G4double m12s = 0.;              // Fake definition before the Loop
  while (rnd > rR)
  {
    m12s = m12sMin + m12sBase*G4UniformRand();
    G4double e1=m12s+fM2-sM2;
    G4double e2=iM2-m12s-tM2;
    G4double four12=4*m12s;
    G4double m13sRange=0.;
    G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
    if(dif<0.)
    {
#ifdef debug
      if(dif<-.01) G4cerr<<"G4QHadron::CopDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
                         <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
#endif
    }
    else m13sRange=sqrt(dif)/m12s;
    rR = m13sRange/m13sBase;
    rnd= G4UniformRand();
#ifdef debug
    G4cout<<"G4QHadron::CopDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
#endif
  }
  G4double m12 = sqrt(m12s);       // Mass of the H1+H2 system
  G4LorentzVector dh4Mom(0.,0.,0.,m12);
  
  if(!CopDecayIn2(t4Mom,dh4Mom,dir,cosp))
  {
    G4cerr<<"***G4QHadron::CopDecayIn3: Exception1"<<G4endl;
    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
    return false;
  }
#ifdef debug
  G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
#endif
  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
  {
    G4cerr<<"***G4QHadron::CopDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
    return false;
  }
  return true;
} // End of CopDecayIn3

// Randomize particle mass taking into account the width
G4double G4QHadron::RandomizeMass(G4QParticle* pPart, G4double maxM)
//       ===========================================================
{
  G4double meanM = theQPDG.GetMass();
  G4double width = theQPDG.GetWidth()/2.;
#ifdef debug
  G4cout<<"G4QHadron::RandomizeMass: meanM="<<meanM<<", halfWidth="<<width<<G4endl;
#endif
  if(maxM<meanM-3*width) 
  {
#ifdef debug
    G4cout<<"***G4QH::RandM:m=0 maxM="<<maxM<<"<meanM="<<meanM<<"-3*halfW="<<width<<G4endl;
#endif
    return 0.;
  }
  ///////////////G4double theMass  = 0.;
  if(width==0.)
  {
#ifdef debug
    if(meanM>maxM) G4cerr<<"***G4QHadron::RandM:Stable m="<<meanM<<">maxM="<<maxM<<G4endl;
#endif
    return meanM;
    //return 0.;
  }
  else if(width<0.)
  {
    G4cerr<<"***G4QHadron::RandM: width="<<width<<"<0,PDGC="<<theQPDG.GetPDGCode()<<G4endl;
    throw G4QException("G4QHadron::RandomizeMass: with the width of the Hadron < 0.");
  }
  G4double minM = pPart->MinMassOfFragm();
  if(minM>maxM)
  {
#ifdef debug
    G4cout<<"***G4QHadron::RandomizeMass:for PDG="<<theQPDG.GetPDGCode()<<" minM="<<minM
          <<" > maxM="<<maxM<<G4endl;
#endif
    return 0.;
  }
  //Now calculate the Breit-Wigner distribution with two cuts
  G4double v1=atan((minM-meanM)/width);
  G4double v2=atan((maxM-meanM)/width);
  G4double dv=v2-v1;
#ifdef debug
  G4cout<<"G4QHadr::RandM:Mi="<<minM<<",i="<<v1<<",Ma="<<maxM<<",a="<<v2<<","<<dv<<G4endl;
#endif
  return meanM+width*tan(v1+dv*G4UniformRand());
}

// Split the hadron in two partons ( quark = anti-diquark v.s. anti-quark = diquark)
void G4QHadron::SplitUp()
{  
  if (isSplit) return;
#ifdef pdebug
  G4cout<<"G4QHadron::SplitUp ***IsCalled***, before Splitting nC="<<Color.size()
        <<", SoftColCount="<<theCollisionCount<<G4endl;
#endif
  isSplit=true;                                     // PutUp isSplit flag to avoid remake
  if (!Color.empty()) return;                       // Do not split if it is already split
  if (!theCollisionCount)                           // Diffractive splitting from particDef
  {
    G4QParton* Left = 0;
    G4QParton* Right = 0;
    GetValenceQuarkFlavors(Left, Right);
    G4ThreeVector Pos=GetPosition();
    Left->SetPosition(Pos);
    Right->SetPosition(Pos);
  
    G4double theMomPlus  = theMomentum.plus();      // E+pz
    G4double theMomMinus = theMomentum.minus();     // E-pz
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Dif* possition="<<Pos<<", 4M="<<theMomentum<<G4endl;
#endif

    // momenta of string ends 
    G4double pt2 = theMomentum.perp2();
    G4double transverseMass2 = theMomPlus*theMomMinus;
    if(transverseMass2<0.) transverseMass2=0.;
    G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2));
    G4ThreeVector pt(0., 0., 0.); // Prototype
    if(maxAvailMomentum2/widthOfPtSquare > 0.01)
                                         pt=GaussianPt(widthOfPtSquare, maxAvailMomentum2);
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Dif* maxMom2="<<maxAvailMomentum2<<", pt="<<pt<<G4endl;
#endif
    G4LorentzVector LeftMom(pt, 0.);
    G4LorentzVector RightMom;
    RightMom.setPx(theMomentum.px() - pt.x());
    RightMom.setPy(theMomentum.py() - pt.y());
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Dif* right4m="<<RightMom<<", left4M="<<LeftMom<<G4endl;
#endif
    G4double Local1 = theMomMinus + (RightMom.perp2() - LeftMom.perp2()) / theMomPlus;
    G4double Local2 = std::sqrt(std::max(0., Local1*Local1 -
                                             4*RightMom.perp2()*theMomMinus / theMomPlus));
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp:Dif,L1="<<Local1<<",L2="<<Local2<<",D="<<Direction<<G4endl;
#endif
    if (Direction) Local2 = -Local2;
    G4double RightMinus   = 0.5*(Local1 + Local2);
    G4double LeftMinus = theMomentum.minus() - RightMinus;
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Dif* Rminus="<<RightMinus<<",Lminus="<<LeftMinus<<",hmm="
          <<theMomentum.minus()<<G4endl;
#endif
    G4double LeftPlus  = 0.;
    if(LeftMinus) LeftPlus  = LeftMom.perp2()/LeftMinus;
    G4double RightPlus = theMomentum.plus() - LeftPlus;
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Dif* Rplus="<<RightPlus<<", Lplus="<<LeftPlus<<G4endl;
#endif
    LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
    LeftMom.setE (0.5*(LeftPlus + LeftMinus));
    RightMom.setPz(0.5*(RightPlus - RightMinus));
    RightMom.setE (0.5*(RightPlus + RightMinus));
    //G4cout<<"DSU 6: Left4M="<<LeftMom<<", Right4M="<<RightMom<<G4endl;
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Dif* -final- R4m="<<RightMom<<", L4M="<<LeftMom<<", L+R="
          <<RightMom+LeftMom<<", D4M="<<theMomentum-RightMom+LeftMom<<G4endl;
#endif
    Left->Set4Momentum(LeftMom);
    Right->Set4Momentum(RightMom);
    Color.push_back(Left);
    AntiColor.push_back(Right);
  }
  else
  {
    // Soft hadronization splitting: sample transversal momenta for sea and valence quarks
    //G4double phi, pts;
    G4ThreeVector SumP(0.,0.,0.);                         // Remember the hadron position
    G4ThreeVector Pos = GetPosition();                    // Remember the hadron position
    G4int nSeaPair = theCollisionCount-1;                 // a#of sea-pairs
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp:*Soft* Pos="<<Pos<<", nSeaPair="<<nSeaPair<<G4endl;
#endif   
    // here the condition,to ensure viability of splitting, also in cases
    // where difractive excitation occured together with soft scattering.
    for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++) // If the sea pairs exist!
    {
      //  choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
      G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress); 

      //  BuildSeaQuark() determines quark spin, isospin and colour 
      //  via parton-constructor G4QParton(aPDGCode) 
      G4QParton* aParton = BuildSeaQuark(false, aPDGCode); // quark/anti-diquark creation

      // save colour a spin-3 for anti-quark
      G4int firstPartonColour = aParton->GetColour();
      G4double firstPartonSpinZ = aParton->GetSpinZ();
#ifdef pdebug
      G4cout<<"G4QHadron::SplitUp:*Soft* Part1 PDG="<<aPDGCode<<", Col="<<firstPartonColour
            <<", SpinZ="<<firstPartonSpinZ<<", 4M="<<aParton->Get4Momentum()<<G4endl;
#endif
      SumP+=aParton->Get4Momentum();
      Color.push_back(aParton);                           // Quark/anti-diquark is filled

      // create anti-quark
      aParton = BuildSeaQuark(true, aPDGCode); // Redefine "aParton" (working pointer)
      aParton->SetSpinZ(-firstPartonSpinZ);
      aParton->SetColour(-firstPartonColour);
#ifdef pdebug
      G4cout<<"G4QHadron::SplUp:Sft,P2="<<aParton->Get4Momentum()<<",i="<<aSeaPair<<G4endl;
#endif

      SumP+=aParton->Get4Momentum();
      AntiColor.push_back(aParton);                       // Anti-quark/diquark is filled
#ifdef pdebug
      G4cout<<"G4QHadron::SplUp:*Sft* Antiquark is filled, i="<<aSeaPair<<G4endl;
#endif
    }
    // ---- Create valence quarks/diquarks
    G4QParton* pColorParton = 0;
    G4QParton* pAntiColorParton = 0;
    GetValenceQuarkFlavors(pColorParton, pAntiColorParton);
    G4int ColorEncoding = pColorParton->GetPDGCode();
#ifdef pdebug
    G4int AntiColorEncoding = pAntiColorParton->GetPDGCode();
    G4cout<<"G4QHadron::SplUp:*Sft*,C="<<ColorEncoding<<", AC="<<AntiColorEncoding<<G4endl;
#endif
    G4ThreeVector ptr = GaussianPt(sigmaPt, DBL_MAX);
    SumP += ptr;
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Sft*, ptr="<<ptr<<G4endl;
#endif

    if (ColorEncoding < 0) // use particle definition
    {
      G4LorentzVector ColorMom(-SumP, 0);
      pColorParton->Set4Momentum(ColorMom);
      G4LorentzVector AntiColorMom(ptr, 0.);
      pAntiColorParton->Set4Momentum(AntiColorMom);
    }
    else
    {
      G4LorentzVector ColorMom(ptr, 0);
      pColorParton->Set4Momentum(ColorMom);
      G4LorentzVector AntiColorMom(-SumP, 0);
      pAntiColorParton->Set4Momentum(AntiColorMom);
    }
    Color.push_back(pColorParton);
    AntiColor.push_back(pAntiColorParton);
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Soft* Col&Anticol are filled PDG="<<GetPDGCode()<<G4endl;
#endif
    // Sample X
    G4int nColor=Color.size();
    G4int nAntiColor=AntiColor.size();
    if(nColor!=nAntiColor || nColor != nSeaPair+1)
    {
      G4cerr<<"***G4QHadron::SplitUp: nA="<<nAntiColor<<",nAC="<<nColor<<",nSea="<<nSeaPair
             <<G4endl;
      G4Exception("G4QHadron::SplitUp:","72",FatalException,"Colours&AntiColours notSinc");
    }
#ifdef pdebug
    G4cout<<"G4QHad::SpUp:,nPartons="<<nColor+nColor<<<<G4endl;
#endif
    G4int dnCol=nColor+nColor;
    // From here two algorithm of splitting can be used (All(default): New, OBO: Olg, Bad)
    G4double* xs=RandomX(dnCol);        // All-Non-iterative CHIPS algorithm of splitting
    // Instead one can try one-by-one CHIPS algorithm (faster? but not exact). OBO comment.
    //G4double Xmax=1.;                   // OBO
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp:*Sft* Loop ColorX="<<ColorX<<G4endl;
#endif
    std::list<G4QParton*>::iterator icolor = Color.begin();
    std::list<G4QParton*>::iterator ecolor = Color.end();
    std::list<G4QParton*>::iterator ianticolor = AntiColor.begin();
    std::list<G4QParton*>::iterator eanticolor = AntiColor.end();
    G4int xi=-1;                        // XIndex for All-Non-interactive CHIPS algorithm
    //G4double X=0.;                      // OBO
    for ( ; icolor != ecolor && ianticolor != eanticolor; ++icolor, ++ianticolor)
    {
      (*icolor)->SetX(xs[++xi]);        // All-Non-iterative CHIPS algorithm of splitting
      //X=SampleCHIPSX(Xmax, dnCol);      // OBO
      //Xmax-=X;                          // OBO
      //--dnCol;                          // OBO
      //(*icolor)->SetX(X);               // OBO
      // ----
      (*icolor)->DefineEPz(theMomentum);
      (*ianticolor)->SetX(xs[++xi]);    // All-Non-iterative CHIPS algorithm of splitting
      //X=SampleCHIPSX(Xmax, dnCol);      // OBO
      //Xmax-=X;                          // OBO
      //--dnCol;                          // OBO
      //(*ianticolor)->SetX(X);           // OBO 
      // ----
      (*ianticolor)->DefineEPz(theMomentum);
    }
    delete[] xs;                           // The calculated array must be deleted (All)
#ifdef pdebug
    G4cout<<"G4QHadron::SplitUp: *Soft* ===> End, ColSize="<<Color.size()<<G4endl;
#endif
    return;
  }
} // End of SplitUp

// Boost hadron 4-momentum, using Boost Lorentz vector
void G4QHadron::Boost(const G4LorentzVector& boost4M)
{  
  // see CERNLIB short writeup U101 for the algorithm
  G4double bm=boost4M.mag();
  G4double factor=(theMomentum.vect()*boost4M.vect()/(boost4M.e()+bm)-theMomentum.e())/bm;
  theMomentum.setE(theMomentum.dot(boost4M)/bm);
  theMomentum.setVect(factor*boost4M.vect() + theMomentum.vect());
} // End of Boost

// Build one (?M.K.) sea-quark
G4QParton* G4QHadron::BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode)
{
  if (isAntiQuark) aPDGCode*=-1;
  G4QParton* result = new G4QParton(aPDGCode);
  result->SetPosition(GetPosition());
  G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
  G4LorentzVector a4Momentum(aPtVector, 0);
  result->Set4Momentum(a4Momentum);
  return result;
} // End of BuildSeaQuark

// Fast non-iterative CHIPS algorithm
G4double* G4QHadron::RandomX(G4int nPart)
{
  G4double* x = 0;
  if(nPart<2)
  {
    G4cout<<"-Warning-G4QHadron::RandomX: nPart="<<nPart<<" < 2"<<G4endl;
    return x;
  }
  x = new G4double[nPart];
  G4int nP1=nPart-1;
  x[0]=G4UniformRand();
  for(G4int i=1; i<nP1; ++i)
  {
    G4double r=G4UniformRand();
    G4int j=0;
    for( ; j<i; ++j) if(r < x[j])
    {
      for(G4int k=i; k>j; --k) x[k]=x[k-1];
      x[j]=r;
      break;
    }
    if(j==i) x[i]=r;
  }
  x[nP1]=1.;
  for(G4int i=nP1; i>0; --i) x[i]-=x[i-1];
  return x;
}

// Non-iterative recursive phase-space CHIPS algorthm
G4double G4QHadron::SampleCHIPSX(G4double anXtot, G4int nSea)
{
  G4double ns=nSea;
  if(nSea<1 || anXtot<=0.) G4cout<<"-Warning-G4QHad::SCX:N="<<nSea<<",tX="<<anXtot<<G4endl;
  if(nSea<2) return anXtot;
  return anXtot*(1.-std::pow(G4UniformRand(),1./ns));
}

// Get flavors for the valence quarks of this hadron
void G4QHadron::GetValenceQuarkFlavors(G4QParton* &Parton1, G4QParton* &Parton2)
{
  // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
  G4int aEnd=0;
  G4int bEnd=0;
  G4int HadronEncoding = GetPDGCode();
  if(!(GetBaryonNumber())) SplitMeson(HadronEncoding, &aEnd, &bEnd);
  else                     SplitBaryon(HadronEncoding, &aEnd, &bEnd);

  Parton1 = new G4QParton(aEnd);
  Parton1->SetPosition(GetPosition());

// G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl;
// G4cerr << "Parton 1: " 
//        << " PDGcode: "  << aEnd
//        << " - Name: "   << Parton1->GetDefinition()->GetParticleName()
//        << " - Type: "   << Parton1->GetDefinition()->GetParticleType() 
//        << " - Spin-3: " << Parton1->GetSpinZ() 
//        << " - Colour: " << Parton1->GetColour() << G4endl;

  Parton2 = new G4QParton(bEnd);
  Parton2->SetPosition(GetPosition());

// G4cerr << "Parton 2: " 
//        << " PDGcode: "  << bEnd
//        << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
//        << " - Type: "   << Parton2->GetDefinition()->GetParticleType() 
//        << " - Spin-3: " << Parton2->GetSpinZ() 
//        << " - Colour: " << Parton2->GetColour() << G4endl;
// G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl;

  // colour of parton 1 choosen at random by G4QParton(aEnd)
  // colour of parton 2 is the opposite:

  Parton2->SetColour(-(Parton1->GetColour()));

 // isospin-3 of both partons is handled by G4Parton(PDGCode)

 // spin-3 of parton 1 and 2 choosen at random by G4QParton(aEnd)
 // spin-3 of parton 2 may be constrained by spin of original particle:

  if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > GetSpin()) 
                                                 Parton2->SetSpinZ(-(Parton2->GetSpinZ()));

// G4cerr << "Parton 2: " 
//        << " PDGcode: "  << bEnd
//        << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
//        << " - Type: "   << Parton2->GetDefinition()->GetParticleType() 
//        << " - Spin-3: " << Parton2->GetSpinZ() 
//        << " - Colour: " << Parton2->GetColour() << G4endl;
// G4cerr << "------------" << G4endl;

} // End of GetValenceQuarkFlavors

G4bool G4QHadron::SplitMeson(G4int PDGcode, G4int* aEnd, G4int* bEnd)
{
  G4bool result = true;
  G4int absPDGcode = std::abs(PDGcode);
  if (absPDGcode >= 1000) return false;
  if(absPDGcode == 22 || absPDGcode == 111) // only u-ubar, d-dbar configurations
  {
    G4int it=1;
    if(G4UniformRand()<.5) it++;
    *aEnd = it;
    *bEnd = -it;
  }
  else if(absPDGcode == 130 || absPDGcode == 310) // K0-K0bar mixing
  {
    G4int it=1;
    if(G4UniformRand()<.5) it=-1;
    *aEnd = it;
    if(it>0) *bEnd = -3;
    else     *bEnd =  3;
  }
  else
  {
    G4int heavy =  absPDGcode/100;
    G4int light = (absPDGcode%100)/10;
    G4int anti  = 1 - 2*(std::max(heavy, light)%2);
    if (PDGcode < 0 ) anti = -anti;
    heavy *=  anti;
    light *= -anti;
    if ( anti < 0) G4SwapObj(&heavy, &light);
    *aEnd = heavy;
    *bEnd = light;
  }
  return result;
}

G4bool G4QHadron::SplitBaryon(G4int PDGcode, G4int* quark, G4int* diQuark)
{
  static const G4double r2=.5;
  static const G4double r3=1./3.;
  static const G4double d3=2./3.;
  static const G4double r4=1./4.;
  static const G4double r6=1./6.;
  static const G4double r12=1./12.;
  //
  std::pair<G4int,G4int> qdq[5];
  G4double               prb[5];
  G4int                  nc=0;
  G4int            aPDGcode=std::abs(PDGcode);
  if(aPDGcode==2212)        // ==> Proton
  {
    nc=3;
    qdq[0]=make_pair(2203, 1); prb[0]=r3;    // uu_1, d
    qdq[1]=make_pair(2103, 2); prb[1]=r6;    // ud_1, u
    qdq[2]=make_pair(2101, 2); prb[2]=r2;    // ud_0, u
  }
  else if(aPDGcode==2112)   // ==> Neutron
  {
    nc=3;
    qdq[0]=make_pair(2103, 1); prb[0]=r6;    // ud_1, d
    qdq[1]=make_pair(2101, 1); prb[1]=r2;    // ud_0, d
    qdq[2]=make_pair(1103, 2); prb[2]=r3;    // dd_1, u
  }
  else if(aPDGcode%10<3)    // ==> Spin 1/2 Hyperons
  {
    if(aPDGcode==3122)      // Lambda
    {
      nc=5;
      qdq[0]=make_pair(2103, 3); prb[0]=r3;  // ud_1, s
      qdq[1]=make_pair(3203, 1); prb[1]=r4;  // su_1, d
      qdq[2]=make_pair(3201, 1); prb[2]=r12; // su_0, d
      qdq[3]=make_pair(3103, 2); prb[3]=r4;  // sd_1, u
      qdq[4]=make_pair(3101, 2); prb[4]=r12; // sd_0, u
    }
    else if(aPDGcode==3222) // Sigma+
    {
      nc=3;
      qdq[0]=make_pair(2203, 3); prb[0]=r3;  // uu_1, s
      qdq[1]=make_pair(3203, 2); prb[1]=r6;  // su_1, d
      qdq[2]=make_pair(3201, 2); prb[2]=r2;  // su_0, d
    }
    else if(aPDGcode==3212) // Sigma0
    {
      nc=5;
      qdq[0]=make_pair(2103, 3); prb[0]=r3;  // ud_1, s
      qdq[1]=make_pair(3203, 1); prb[1]=r12; // su_1, d
      qdq[2]=make_pair(3201, 1); prb[2]=r4;  // su_0, d
      qdq[3]=make_pair(3103, 2); prb[3]=r12; // sd_1, u
      qdq[4]=make_pair(3101, 2); prb[4]=r4;  // sd_0, u
    }
    else if(aPDGcode==3112) // Sigma-
    {
      nc=3;
      qdq[0]=make_pair(1103, 3); prb[0]=r3;  // dd_1, s
      qdq[1]=make_pair(3103, 1); prb[1]=r6;  // sd_1, d
      qdq[2]=make_pair(3101, 1); prb[2]=r2;  // sd_0, d
    }
    else if(aPDGcode==3312) // Xi-
    {
      nc=3;
      qdq[0]=make_pair(3103, 3); prb[0]=r6;  // sd_1, s
      qdq[1]=make_pair(3101, 3); prb[1]=r2;  // sd_0, s
      qdq[2]=make_pair(3303, 1); prb[2]=r3;  // ss_1, d
    }
    else if(aPDGcode==3322) // Xi0
    {
      nc=3;
      qdq[0]=make_pair(3203, 3); prb[0]=r6;  // su_1, s
      qdq[1]=make_pair(3201, 3); prb[1]=r2;  // su_0, s
      qdq[2]=make_pair(3303, 2); prb[2]=r3;  // ss_1, u
    }
    else return false;
  }
  else                                       // ==> Spin 3/2 Baryons (Spin>3/2 is ERROR)
  {
    if(aPDGcode==3334)
    {
      nc=1;
      qdq[0]=make_pair(3303, 3); prb[0]=1.;  // ss_1, s
    }
    else if(aPDGcode==2224)
    {
      nc=1;
      qdq[0]=make_pair(2203, 2); prb[0]=1.;  // uu_1, s
    }
    else if(aPDGcode==2214)
    {
      nc=2;
      qdq[0]=make_pair(2203, 1); prb[0]=r3;  // uu_1, d
      qdq[1]=make_pair(2103, 2); prb[1]=d3;  // ud_1, u
    }
    else if(aPDGcode==2114)
    {
      nc=2;
      qdq[0]=make_pair(1103, 2); prb[0]=r3;  // dd_1, u
      qdq[1]=make_pair(2103, 1); prb[1]=d3;  // ud_1, d
    }
    else if(aPDGcode==1114)
    {
      nc=1;
      qdq[0]=make_pair(1103, 1); prb[0]=1.;  // uu_1, s
    }
    else if(aPDGcode==3224)
    {
      nc=2;
      qdq[0]=make_pair(2203, 3); prb[0]=r3;  // uu_1, s
      qdq[1]=make_pair(3203, 2); prb[1]=d3;  // su_1, u
    }
    else if(aPDGcode==3214) // @@ SU(3) is broken because of the s-quark mass
    {
      nc=3;
      qdq[0]=make_pair(2103, 3); prb[0]=r3;  // ud_1, s
      qdq[1]=make_pair(3203, 1); prb[1]=r3;  // su_1, d
      qdq[2]=make_pair(3103, 2); prb[2]=r3;  // sd_1, u
    }
    else if(aPDGcode==3114)
    {
      nc=2;
      qdq[0]=make_pair(1103, 3); prb[0]=r3;  // dd_1, s
      qdq[1]=make_pair(3103, 1); prb[1]=d3;  // sd_1, d
    }
    else if(aPDGcode==3324)
    {
      nc=2;
      qdq[0]=make_pair(3203, 3); prb[0]=r3;  // su_1, s
      qdq[1]=make_pair(3303, 2); prb[1]=d3;  // ss_1, u
    }
    else if(aPDGcode==3314)
    {
      nc=2;
      qdq[0]=make_pair(3103, 3); prb[0]=d3;  // sd_1, s
      qdq[1]=make_pair(3303, 1); prb[1]=r3;  // ss_1, d
    }
    else return false;
  }
  G4double random = G4UniformRand();
  G4double sum = 0.;
  for(G4int i=0; i<nc; i++)
  {
    sum += prb[i];
    if(sum>random)
    {
      if(PDGcode>0)
      {
        *diQuark= qdq[i].first;
        *quark  = qdq[i].second;
      }
      else
      {
        *diQuark= -qdq[i].second;
        *quark  = -qdq[i].first;
      }
      break;
    }
  }
  return true;
}

// This is not usual Gaussian, in fact it is dN/d(pt) ~ pt * exp(-pt^2/pt0^2)
G4ThreeVector G4QHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
{
  G4double R=0.;
  while ((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare){}
  R = std::sqrt(R);
  G4double phi = twopi*G4UniformRand();
  return G4ThreeVector(R*std::cos(phi), R*std::sin(phi), 0.);    
}

G4QParton* G4QHadron::GetNextParton()
{
  if(Color.size()==0) return 0;
  G4QParton* result = Color.back();
  Color.pop_back();
  return result;
}

G4QParton* G4QHadron::GetNextAntiParton()
{
  if(AntiColor.size() == 0) return 0;
  G4QParton* result = AntiColor.front();
  AntiColor.pop_front();
  return result;
}

// Random Split of the Hadron in 2 Partons (caller is responsible for G4QPartonPair delete)
G4QPartonPair* G4QHadron::SplitInTwoPartons() // If result=0: impossible to split (?)
{
  if(std::abs(GetBaryonNumber())>1) // Not Baryons or Mesons or Anti-Baryons
  {
    G4cerr<<"***G4QHadron::SplitInTwoPartons: Can not split QC="<<valQ<< G4endl;
    G4Exception("G4QFragmentation::ChooseX:","72",FatalException,"NeitherMesonNorBaryon");
  }
  std::pair<G4int,G4int> PP = valQ.MakePartonPair();
  return new G4QPartonPair(new G4QParton(PP.first), new G4QParton(PP.second));
}