#include "PHOTONS++/PhaseSpace/Avarage_Photon_Number.H"
#include "ATOOLS/Math/Random.H"
#include "ATOOLS/Math/Vector.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Phys/Particle.H"
#include "PHOTONS++/Main/Photons.H"
#include "PHOTONS++/Main/Dipole_Type.H"

using namespace PHOTONS;
using namespace ATOOLS;
using namespace std;

Avarage_Photon_Number::Avarage_Photon_Number
(const Particle_Vector& dip, const double& wmax, const double& wmin) :
  m_omegaMax(wmax), m_omegaMin(wmin), m_dipole(dip), m_nbar(0.)
{
  DEBUG_FUNC(dip.size());
  CalculateAvaragePhotonNumber();
  msg_Debugging()<<"all dipoles:\n";
  for (unsigned int i=0; i<m_nbars.size(); i++) {
    msg_Debugging()<<m_nbars[i]
                   <<": ["<<m_dipole[m_nbars[i].ij.i]->Flav()<<","
                   <<m_dipole[m_nbars[i].ij.j]->Flav()<<"]"<<std::endl;
  }
  msg_Debugging()<<"-> "<<m_nbar<<endl;
}

Avarage_Photon_Number::~Avarage_Photon_Number() {
}

void Avarage_Photon_Number::CalculateAvaragePhotonNumber() {
  double sum      = 0.;
  for(unsigned int j=0; j<m_dipole.size(); j++) {
    for(unsigned int i=0; i<j; i++) {
      double Zi       = m_dipole[i]->Flav().Charge();
      double Zj       = m_dipole[j]->Flav().Charge();
      double titj     = TiTj(i,j);
      double betai    = CalculateBeta(m_dipole[i]->Momentum());
      double betaj    = CalculateBeta(m_dipole[j]->Momentum());
      Vec3D  pi       = m_dipole[i]->Momentum();
      Vec3D  pj       = m_dipole[j]->Momentum();
      double alpha    = 0.;
      if (!IsZero(pi.Abs()) && !IsZero(pj.Abs())) {
        double arg = pi*pj/(pi.Abs()*pj.Abs());
        if (IsEqual(arg,-1.))     alpha = 0.;
        else if (IsEqual(arg,1.)) alpha = M_PI/2.;
        else                      alpha = (M_PI - acos(arg))/2.;
      }
      double ai       = betai*sin(alpha);
      double aj       = betaj*sin(alpha);
      double bi       = betai*cos(alpha);
      double bj       = betaj*cos(alpha);
      double intterm  = 0.;
      if (IsZero(alpha,2.e-4)) {
        intterm = log(((1.+betai)*(1.+betaj))/((1.-betai)*(1.-betaj)))
                  /(betai+betaj);
        msg_Debugging()<<"back-to-back pair discovered: "<<intterm
                       <<" <-> "<<InterferenceTerm(ai,aj,bi,bj)<<std::endl;
      }
      else {
        intterm = InterferenceTerm(ai,aj,bi,bj);
      }
      double dipoleij = Photons::s_alpha/M_PI*Zi*Zj*titj
			  *log(m_omegaMax/m_omegaMin)
			  *(2.-(1.-ai*aj+bi*bj)*intterm);
#ifdef PHOTONS_DEBUG
      msg_Debugging()<<"ana: "<<dipoleij<<endl;
      msg_Debugging()<<"alpha/pi*ZiZjtitj: "
                     <<Photons::s_alpha/M_PI*Zi*Zj*titj<<endl;
      msg_Debugging()<<"log(Emax/Emin): "<<log(m_omegaMax/m_omegaMin)<<endl;
      msg_Debugging()<<"int: "
                     <<(1-ai*aj+bi*bj)*intterm<<endl;
      msg_Debugging()<<"Emax,Emin: "<<m_omegaMax<<" ,  "<<m_omegaMin<<endl;
#endif
      // check we are entering something sensible, nan's will also fail
      if (intterm!=0.) m_nbars.push_back(IdPairNbar(IdPair(i,j),dipoleij));
      else             m_nbars.push_back(IdPairNbar(IdPair(i,j),0.));
      sum += m_nbars.back().nbar;
    }
  }
  m_nbar = sum;
}

double Avarage_Photon_Number::CalculateBeta(const Vec4D& p) {
  return Vec3D(p).Abs()/p[0];
}

double Avarage_Photon_Number::InterferenceTerm
(const double& ai, const double& aj, const double& bi, const double& bj) {
  if (IsZero(ai) && IsZero(aj) && IsZero(bi) && !IsZero(bj))
    return 1./bj*log((1.+bj)/(1.-bj));
  if (IsZero(ai) && IsZero(aj) && !IsZero(bi) && IsZero(bj))
    return 1./bi*log((1.+bi)/(1.-bi));
  double A  = bi-bj;
  double B  = 2.*bi*bj;
  double Ci = 1.-ai*ai;
  double Cj = 1.-aj*aj;
  double Di = -2.*bi;
  double Dj = 2.*bj;
  double Ei = ai*ai+bi*bi;
  double Ej = aj*aj+bj*bj;
  double rooti = sqrt(B*B*Ci-A*B*Di+A*A*Ei);
  double rootj = sqrt(B*B*Cj-A*B*Dj+A*A*Ej);
  // avoid divergences in A+-B when alpha = 0 and beta1=1 and beta2=0.3333
  if (ai/bi < 1E-4) // ai/bi = tan(alpha)
    return 1./(bi+bj)*log(((1.+bi)*(1.+bj))/((1.-bi)*(1.-bj)));
  // avoid nan's (brute force sofar)
  else if (!(2.*sqrt(Ci-Di+Ei)+(B*(2.*Ci-Di)-A*(Di-2.*Ei))/rooti>0) ||
           !(2.*sqrt(Cj-Dj+Ej)+(B*(2.*Cj-Dj)-A*(Dj-2.*Ej))/rootj>0) ||
           !(2.*sqrt(Ci+Di+Ei)+(B*(2.*Ci+Di)-A*(Di+2.*Ei))/rooti>0) ||
           !(2.*sqrt(Cj+Dj+Ej)+(B*(2.*Cj+Dj)-A*(Dj+2.*Ej))/rootj>0)) return 0.;
  else
    return bi*bj*(1./(bj*rooti)
                  *log(abs((A+B)/(A-B))
                       *(2.*sqrt(Ci-Di+Ei)+(B*(2.*Ci-Di)-A*(Di-2.*Ei))/rooti)
                       /(2.*sqrt(Ci+Di+Ei)+(B*(2.*Ci+Di)-A*(Di+2.*Ei))/rooti))
                 -1./(bi*rootj)
                  *log(abs((A+B)/(A-B))
                       *(2.*sqrt(Cj-Dj+Ej)+(B*(2.*Cj-Dj)-A*(Dj-2.*Ej))/rootj)
                       /(2.*sqrt(Cj+Dj+Ej)+(B*(2.*Cj+Dj)-A*(Dj+2.*Ej))/rootj)));
}

double Avarage_Photon_Number::TiTj(const size_t& i, const size_t& j)
{
  if      ((m_dipole[i]->ProductionBlob() == m_dipole[j]->ProductionBlob()) &&
           (m_dipole[i]->ProductionBlob() != NULL))  return  1.;
  else if ((m_dipole[i]->DecayBlob() == m_dipole[j]->ProductionBlob()) &&
           (m_dipole[i]->DecayBlob() != NULL))       return -1.;
  else if ((m_dipole[i]->ProductionBlob() == m_dipole[j]->DecayBlob()) &&
           (m_dipole[i]->ProductionBlob() != NULL))  return -1.;
  else if ((m_dipole[i]->DecayBlob() == m_dipole[j]->DecayBlob()) &&
           (m_dipole[i]->DecayBlob() != NULL))       return  1.;
  else                                               return  0.;
}