source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QHadronBuilder.cc @ 836

//
// ********************************************************************
// * 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: G4QHadronBuilder.cc,v 1.2 2006/12/12 11:02:22 mkossov Exp $
// GEANT4 tag $Name:  $
//
// -----------------------------------------------------------------------------
//      GEANT4 class implementation file
//
//      Created by Mikhail Kosov, October 2006 
//      Simple algorithm making a hadron out of two partons
//      For comparison mirror member functions are taken from G4 class:
//      G4HadronBuilder
// -----------------------------------------------------------------------------

//#define debug

#include "G4QHadronBuilder.hh"
#include "Randomize.hh"
#include "G4ParticleTable.hh"

G4QHadronBuilder::G4QHadronBuilder()
{
         mesonSpinMix = 0.5;                        // probability to create vector meson 
         baryonSpinMix= 0.5;                     // probability to create 3/2 baryon 
  scalarMesonMixings.resize(6);
  scalarMesonMixings[0] = 0.5; 
  scalarMesonMixings[1] = 0.25; 
  scalarMesonMixings[2] = 0.5; 
  scalarMesonMixings[3] = 0.25; 
  scalarMesonMixings[4] = 1.0; 
  scalarMesonMixings[5] = 0.5; 
  vectorMesonMixings.resize(6);
  vectorMesonMixings[0] = 0.5;
  vectorMesonMixings[1] = 0.0;
  vectorMesonMixings[2] = 0.5;
  vectorMesonMixings[3] = 0.0;
  vectorMesonMixings[4] = 1.0;
  vectorMesonMixings[5] = 1.0; 
}

G4QHadron* G4QHadronBuilder::Build(G4QParton* black, G4QParton* white)
{
         if(black->GetParticleSubType()=="di_quark" || white->GetParticleSubType()=="di_quark")
  {
    // Baryon consists of quark and at least one di-quark
           Spin spin = (G4UniformRand() < baryonSpinMix) ? SpinHalf : SpinThreeHalf;
           return Baryon(black,white,spin);
         }
  else
  {     
    // Meson consists of quark and abnti-quark
           Spin spin = (G4UniformRand() < mesonSpinMix) ? SpinZero : SpinOne;
           return Meson(black,white,spin);
         }
}

G4QHadron* G4QHadronBuilder::BuildLowSpin(G4QParton* black, G4QParton* white)
{
         if(black->GetParticleSubType() == "quark" && white->GetParticleSubType() == "quark")
                      return Meson(black,white, SpinZero);
  //    returns a SpinThreeHalf Baryon if all quarks are the same
         else           return Baryon(black,white, SpinHalf);
}

G4QHadron* G4QHadronBuilder::BuildHighSpin(G4QParton* black, G4QParton* white)
{
        if(black->GetParticleSubType() == "quark" && white->GetParticleSubType() == "quark")
                     return Meson(black,white, SpinOne);
         else return Baryon(black,white,SpinThreeHalf);
}

G4QHadron* G4QHadronBuilder::Meson(G4QParton* black, G4QParton* white, Spin theSpin)
{
#ifdef debug
  //  Verify Input Charge
  G4double charge =  black->GetPDGCharge() + white->GetPDGCharge();      
  if (std::abs(charge)>2 || std::abs(3*(charge-G4int(charge*1.001))) > perCent)
           G4cerr<<"-Warning-G4QHadronBuilder::Meson: invalid TotalCharge="<<charge<<", PDG_q1=" 
          <<black->GetPDGEncoding()<< ", PDG_q2="<<white->GetPDGEncoding()<<G4endl;
#endif  
        
         G4int id1= black->GetPDGCode();
         G4int id2= white->GetPDGCode();
         if (std::abs(id1) < std::abs(id2))     // exchange black and white
         {
           G4int xchg = id1; 
           id1 = id2;  
           id2 = xchg;
         }
         if(std::abs(id1)>3) 
         {
                                G4cerr<<"***G4QHadronBuilder::Meson: q1="<<id1<<", q2="<<id2
          <<" while CHIPS is only SU(3)"<<G4endl;
    G4Exception("G4QHadronBuilder::Meson:","72",FatalException,"HeavyQuarkFound");
  }     
  G4int PDGEncoding=0;
         if(!(id1+id2))                      // annihilation case (neutral)
  {     
           G4double rmix = G4UniformRand();
           G4int    imix = 2*std::abs(id1) - 1;
           if(theSpin == SpinZero)
              PDGEncoding = 110*(1 + G4int(rmix + scalarMesonMixings[imix - 1])
                                   + G4int(rmix + scalarMesonMixings[imix]    ) ) +  theSpin;
                  else
              PDGEncoding = 110*(1 + G4int(rmix + vectorMesonMixings[imix - 1])
                                                    + G4int(rmix + vectorMesonMixings[imix]    ) ) +  theSpin;
         }
  else
  {
           PDGEncoding = 100 * std::abs(id1) + 10 * std::abs(id2) +  theSpin;  
           G4bool IsUp = (std::abs(id1)&1) == 0;        // quark 1 is up type quark (u or c?)
           G4bool IsAnti = id1 < 0;                        // quark 1 is an antiquark?
           if( (IsUp && IsAnti) || (!IsUp && !IsAnti) )  PDGEncoding = - PDGEncoding;
         }
         G4QHadron* Meson= new G4QHadron(PDGEncoding);
#ifdef debug
         if(std::abs(black->GetPDGCharge() + white->GetPDGCharge() - Meson->GetCharge()) > .001)
           G4cout<<"-Warning-G4QHadronBuilder::Meson:wrongCharge, q1="<<black->GetPDGCharge()<<"("
          <<black->->GetParticleName()<<"), q2="<<white->GetPDGCharge()<<"("
          <<white->->GetParticleName()<<"), qM="<<Meson->GetCharge()<<"/"<<PDGEncoding
          <<G4endl;
#endif
        return Meson;
}

G4QHadron* G4QHadronBuilder::Baryon(G4QParton* black, G4QParton* white, Spin theSpin)
{
#ifdef debug
  //  Verify Input Charge
  G4double charge =  black->GetPDGCharge() + white->GetPDGCharge();      
  if(std::abs(charge) > 2 || std::abs(3*(charge - G4int(charge*1.001))) > perCent )
           G4cerr<<"-Warning-G4QHadronBuilder::Baryon: invalid TotalCharge="<<charge<<", PDG_q1=" 
          <<black->GetPDGEncoding()<< ", PDG_q2="<<white->GetPDGEncoding()<<G4endl;
#endif  
         G4int id1= black->GetPDGCode();
         G4int id2= white->GetPDGCode();
         if(std::abs(id1) < std::abs(id2))
         {
           G4int xchg = id1; 
           id1 = id2;  
           id2 = xchg;
         }
         if(std::abs(id1)<1000 || std::abs(id2)> 3) 
         {
                                G4cerr<<"***G4QHadronBuilder::Baryon: q1="<<id1<<", q2="<<id2
          <<" can't create a Baryon"<<G4endl;
    G4Exception("G4QHadronBuilder::Baryon:","72",FatalException,"WrongQdQSequence");
  }
         G4int ifl1= std::abs(id1)/1000;
         G4int ifl2 = (std::abs(id1) - ifl1 * 1000)/100;
         G4int diquarkSpin = std::abs(id1)%10; 
         G4int ifl3 = id2;
         if (id1 < 0) {ifl1 = - ifl1; ifl2 = - ifl2;}
         //... Construct baryon, distinguish Lambda and Sigma baryons.
         G4int kfla = std::abs(ifl1);
         G4int kflb = std::abs(ifl2);
         G4int kflc = std::abs(ifl3);
         G4int kfld = std::max(kfla,kflb);
               kfld = std::max(kfld,kflc);
         G4int kflf = std::min(kfla,kflb);
               kflf = std::min(kflf,kflc);
         G4int kfle = kfla + kflb + kflc - kfld - kflf;
         //... baryon with content uuu or ddd or sss has always spin = 3/2
         if(kfla==kflb && kflb==kflc) theSpin=SpinThreeHalf;   

         G4int kfll = 0;
         if(theSpin == SpinHalf && kfld > kfle && kfle > kflf)
  { 
    // Spin J=1/2 and all three quarks different
    // Two states exist: (uds -> lambda or sigma0)
    //   -  lambda: s(ud)0 s : 3122; ie. reverse the two lighter quarks
    //   -  sigma0: s(ud)1 s : 3212
           if(diquarkSpin == 1 )
    {
              if ( kfla == kfld) kfll = 1; // heaviest quark in diquark
              else kfll = G4int(0.25 + G4UniformRand());
           }   
           if(diquarkSpin==3 && kfla!=kfld) kfll = G4int(0.75+G4UniformRand());
         }
         G4int PDGEncoding=0;
         if (kfll == 1) PDGEncoding = 1000 * kfld + 100 * kflf + 10 * kfle + theSpin;
         else           PDGEncoding = 1000 * kfld + 100 * kfle + 10 * kflf + theSpin;
        if (id1 < 0) PDGEncoding = -PDGEncoding;
         G4QHadron* Baryon= new G4QHadron(PDGEncoding);
#ifdef debug
         if(std::abs(black->GetPDGCharge() + white->GetPDGCharge()- Baryon->GetCharge()) > .001)
           G4cout<<"-Warning-G4QHadronBuilder::Baryon:wrongCharge,dq="<<black->GetPDGCharge()<<"("
          <<black->->GetParticleName()<<"), q="<<white->GetPDGCharge()<<"("
          <<white->->GetParticleName()<<"), qB="<<Baryon->GetCharge()<<"/"<<PDGEncoding
          <<G4endl;
#endif
        return Baryon;
}