#include "SigmaSgluonPairs.h"

Sigma2qqbar2sgluonsgluonBar::Sigma2qqbar2sgluonsgluonBar(){
  std::cout << "ERROR YOU SHOULD NOT HAVE LANDED HERE:  Sigma2qqbar2sgluonsgluonBar::Sigma2qqbar2sgluonsgluonBar() " << std::endl;
}
Sigma2qqbar2sgluonsgluonBar::Sigma2qqbar2sgluonsgluonBar(int idResonanceIn) : Sigma2SgluonSgluonBase(idResonanceIn) {
}
Sigma2qqbar2sgluonsgluonBar::~Sigma2qqbar2sgluonsgluonBar(){
}
// Initialize process. 
  
void Sigma2qqbar2sgluonsgluonBar::initProc() {

  // Set pointer to particle properties and decay table.
  m_particlePtr = particleDataPtr->particleDataEntryPtr(m_idsgluon);
  
} 

//--------------------------------------------------------------------------

double Sigma2qqbar2sgluonsgluonBar::sigmaHat() {
  return m_sigma;
}
// Info on the subprocess.
string Sigma2qqbar2sgluonsgluonBar::name()       const {
  return "q qbar -> sgluon sgluonBar";
}
int    Sigma2qqbar2sgluonsgluonBar::code()       const {
  return 10400;
}
string Sigma2qqbar2sgluonsgluonBar::inFlux()     const {
  return "qqbarSame";
}
int    Sigma2qqbar2sgluonsgluonBar::resonanceA() const {
  return m_idsgluon;
}

int Sigma2qqbar2sgluonsgluonBar::id3Mass() const {
  return abs(m_idsgluon);
}
int Sigma2qqbar2sgluonsgluonBar::id4Mass() const {
  return abs(m_idsgluonBar);
}

// Evaluate sigmaHat(sHat); first step when inflavours unknown. 

void Sigma2qqbar2sgluonsgluonBar::sigmaKin() { 

  // here we calculate the cross section independent of flavour
  m_sigma = 0.;

  // check that sHat > Q2=4*mSgluon**2
  double sgluonMass = m_particlePtr->m0();

  if ( sH <= 4.*sgluonMass*sgluonMass ) {
    return;
  }

  // convert tH to costheta
  double ctp = 2*tH-2*sgluonMass*sgluonMass+sH;
  ctp = ctp/(sqrtf(sH*(sH-4*sgluonMass*sgluonMass)));
  if (fabs(ctp)>1.e0 ) return;

  // Unit is GeV-2
  m_sigma = dsigmadcostheta(ctp);

  // translate from dsigma/dcostheta to dsigma/dtheta:
  // dt = s/2 * sqrt(1-4m^2/s) * dcostheta
  double factor = convFactordsigmadcostheta2dsigmadt(sH,sgluonMass);
  m_sigma = m_sigma/factor;
}

double Sigma2qqbar2sgluonsgluonBar::dsigmadcostheta(double ctp) {
  
  double sgluonMass = m_particlePtr->m0();
  double Q2     = sgluonMass*sgluonMass;
  double alphaS = couplingsPtr->alphaS(Q2);

  double beta=sqrtf(1.e0-4.e0*sgluonMass*sgluonMass/sH);

  double crossSection = (4.*M_PI*alphaS*alphaS)/(9.*sH)*beta*beta*beta;
  crossSection = crossSection*0.75e0*(1.e0-ctp*ctp);

  return crossSection;
}

// Select identity, colour and anticolour.

void Sigma2qqbar2sgluonsgluonBar::setIdColAcol() {

  // Flavours trivial.
  setId( id1, id2, m_idsgluon, m_idsgluonBar);

  // Colour flow topologies. Swap when antiquarks.
  setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
  if (id1 < 0) swapColAcol();

}

Sigma2gg2sgluonsgluonBar::Sigma2gg2sgluonsgluonBar(){
  std::cout << "ERROR YOU SHOULD NOT HAVE LANDED HERE: Sigma2gg2sgluonsgluonBar::Sigma2gg2sgluonsgluonBar() " << std::endl;
}
Sigma2gg2sgluonsgluonBar::Sigma2gg2sgluonsgluonBar(int idResonanceIn) : Sigma2SgluonSgluonBase(idResonanceIn) {
}
Sigma2gg2sgluonsgluonBar::~Sigma2gg2sgluonsgluonBar(){
}
// Here we have gg -> sgluon sgluonBar
double Sigma2gg2sgluonsgluonBar::sigmaHat() {
  return m_sigma;
}
string Sigma2gg2sgluonsgluonBar::name()       const {
  return "g g    -> sgluon sgluonBar";
}
int    Sigma2gg2sgluonsgluonBar::code()       const {
  return 10401;
}
string Sigma2gg2sgluonsgluonBar::inFlux()     const {
  return "gg";
}
int    Sigma2gg2sgluonsgluonBar::resonanceA() const {
  return m_idsgluon;
}
int Sigma2gg2sgluonsgluonBar::id3Mass() const {
  return abs(m_idsgluon);
}
int Sigma2gg2sgluonsgluonBar::id4Mass() const {
  return abs(m_idsgluonBar);
}


// Initialize process. 
  
void Sigma2gg2sgluonsgluonBar::initProc() {

  // Set pointer to particle properties and decay table.
  m_particlePtr = particleDataPtr->particleDataEntryPtr(m_idsgluon);
  
} 

//--------------------------------------------------------------------------

// Evaluate sigmaHat(sHat); first step when inflavours unknown. 

void Sigma2gg2sgluonsgluonBar::sigmaKin() { 

  // here we calculate the cross section independent of flavour
  m_sigma = 0.;

  // check that sHat > Q2=4*mSgluon**2
  double sgluonMass = m_particlePtr->m0();

  if ( sH <= 4.*sgluonMass*sgluonMass ) {
    return;
  }

  double ctp = 2*tH-2*sgluonMass*sgluonMass+sH;
  ctp = ctp/(sqrtf(sH*(sH-4*sgluonMass*sgluonMass)));
  if (fabs(ctp)>1.e0 ) return;

  // unit is GeV-2
  m_sigma = dsigmadcostheta(ctp);

  // translate from dsigma/dcostheta to dsigma/dtheta:
  // dt = s/2 * sqrt(1-4m^2/s) * dcostheta
  double factor = convFactordsigmadcostheta2dsigmadt(sH,sgluonMass);
  m_sigma = m_sigma/factor;

}

//--------------------------------------------------------------------------

double Sigma2gg2sgluonsgluonBar::dsigmadcostheta(double ctp) {

  // process gg -> sigma sigma differential cross section

  double sgluonMass = m_particlePtr->m0();
  double Q2 = sgluonMass*sgluonMass;

  double alphaS = couplingsPtr->alphaS(Q2);
	     
  double beta= sqrtf(1.e0-4.e0*sgluonMass*sgluonMass/sH);
  double beta2 = beta*beta;
  
  double crossSection = M_PI*(alphaS*alphaS)*beta/(6.*sH)*(9./8.)*
    (27.-17.*beta2 - 6.*(1.e0-beta2)*(3.e0+beta2)
     /(2.e0*beta)*log((1.e0+beta)/(1.e0-beta)));     
  
  double ctp2 = ctp*ctp;
  double stp4 = (1.e0-ctp2)*(1.e0-ctp2);
  
  double cr = 1.5e0*((1.e0-beta2)*(1.e0-beta2) + beta2*beta2*stp4)/
    (1.e0-beta2*ctp2)*(1.e0-beta2*ctp2) * (3.e0+beta2*ctp2)/
    ( 27.e0-17.e0*beta2-6.e0*(1.e0-beta2)*(3.e0+beta2)/
      (2.e0*beta)*log((1.e0+beta)/(1.e0-beta)));

  crossSection = crossSection*cr;

  return crossSection;
}

// Select identity, colour and anticolour.

void Sigma2gg2sgluonsgluonBar::setIdColAcol() {

  // Flavours trivial.
  setId( id1, id2, m_idsgluon, m_idsgluonBar);

  // Colour flow topologies of 4 colour octets
  setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);

}