#include "PHASIC++/Channels/Channel_Elements_KK.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Math/Vector.H"
#include "ATOOLS/Math/MathTools.H"
#include "ATOOLS/Math/Random.H"
#include "ATOOLS/Org/My_MPI.H"
#include "MODEL/Main/Model_Base.H"

using namespace PHASIC;
using namespace ATOOLS;
using namespace std;

Channel_Elements_KK PHASIC::CEKK;

Channel_Elements_KK::Channel_Elements_KK()
{
  m_ed    = 0;
  m_rs    = 0;
}

Channel_Elements_KK::~Channel_Elements_KK()
{
}

void Channel_Elements_KK::Init(int nin,int nout,ATOOLS::Flavour* fl)
{
  if (m_rs>0) return;
  m_nin = nin;
  m_nout = nout;
  m_kkp=-1;m_mpss=1.;
  int mode = MODEL::s_model->ScalarNumber(std::string("KK_mode"));
  for (int i=m_nin;i<m_nin+m_nout;i++) {
    if(fl[i].IsKK() && (mode==1 || mode==2 || mode==5)){
      if(ATOOLS::IsZero(ATOOLS::sqr(fl[i].Mass()))){
	msg_Error()<<"Error in Channel_Elements_KK: "<<endl
		   <<"   Please initialize with nonzero particle mass ("<<fl[i]<<") !"<<std::endl;
	Abort();
      }
      m_kkp=i;
      
      m_ed  = MODEL::s_model->ScalarNumber(std::string("ED"));
      m_r2  = sqr(MODEL::s_model->ScalarConstant(std::string("Radius")));
      m_gn  = MODEL::s_model->ScalarConstant(std::string("G_Newton"));

      //Calculation of Gamma(ed/2)
      if(m_ed%2==0) m_gam=1.;
      else m_gam=sqrt(M_PI);
      for(int k=2-m_ed%2;k<m_ed;k+=2)m_gam*=0.5*k;

      m_prevET = 0.;
      break;
    }
  }
  m_rs = 1;
}

void Channel_Elements_KK::SetKKmass(double *ms, double ET, Cut_Data* cuts,double ran)
{
  if (!IsEqual(ET,m_prevET) && m_kkp>-1) {
    double mm = m_prevET = ET;
    for(int j=m_nin;j<m_nin+m_nout;j++)
      if(j!=m_kkp) mm -= Max(sqrt(ms[j]),cuts->etmin[j]);
    if (m_nout==2) mm = Min(mm,sqrt(sqr(ET)-2*cuts->etmin[5-m_kkp]*ET));
    m_maxm2 = sqr(mm);
    m_maxn  = sqrt(m_maxm2*m_r2);
    m_mpss  = 2.*pow(sqrt(M_PI)*m_maxn,double(m_ed))/m_gam;
  }
  m_sran = ran;
  double ms2=sqr(ran)*m_maxm2;
  m_weight=pow(ran,double(m_ed-1));

  ms[m_kkp]=ms2;
}

double Channel_Elements_KK::GetWeightKK(double& ran){  
  ran = m_sran;
  return m_mpss*m_weight;
}