/// -*- C++ -*-
#include "Rivet/Rivet.hh"
#include "BroadeningGamma.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Cmp.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Math/Constants.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "math.h"
#include "Rivet/Projection.hh"
#include "Rivet/Projections/AxesDefinition.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Event.hh"
#include "Rivet/Particle.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Projections/DISLepton.hh"
#include "Rivet/Projections/Beam.hh"

namespace Rivet {
    
  int BroadeningGamma::compare(const Projection& p) const {
    return mkNamedPCmp(p, "Kinematics");
  } 

  
  void BroadeningGamma::calcbg(const vector<Particle>& fsparticles, FourMomentum pgamma) {
    vector<Vector3> threeMomenta;
    threeMomenta.reserve(fsparticles.size());
    foreach (const Particle& p, fsparticles) {
      const Vector3 p3 = p.momentum().vector3();
      threeMomenta.push_back(p3);
    }
    _calbgBroadeningGamma(threeMomenta,pgamma);
  }
  
  void BroadeningGamma::calcbg(const vector<FourMomentum>& fsmomenta, FourMomentum pgamma) {
    vector<Vector3> threeMomenta;
    threeMomenta.reserve(fsmomenta.size());
    foreach (const FourMomentum& v, fsmomenta) {
      threeMomenta.push_back(v.vector3());
    }
    _calbgBroadeningGamma(threeMomenta, pgamma);
  }
  
  void BroadeningGamma::calcbg(const vector<Vector3>& fsmomenta, FourMomentum pgamma) {
    _calbgBroadeningGamma(fsmomenta, pgamma);
  }
  
  
  /////////////////////////////////////////////////
  
  
  
  inline bool mod2Cmp(const Vector3& a, const Vector3& b) {
    return a.mod2() > b.mod2();
  }
  
  
  // Do the general case thrustgamma calbgulation
  void _calcBG(const vector<Vector3>& momenta, double& b, Vector3& taxis, FourMomentum pgamma) {
   

    vector<Vector3> p = momenta;
            
    Vector3 pgam = pgamma.vector3();
    ///make a 3vector from the fourmomentum of the photon

      
    pgam = pgam.unit();
      // Calbulate the thrustgamma value 
      b=0.;
      
      for (unsigned int k=0 ; k<p.size() ; k++){
	b += mod(pgam.cross(p[k]));    
       
      }

      // Store everything
     taxis = pgam;
    
    
 
      }
  

  
  
  
  // Do the full calbulation
  void BroadeningGamma::_calbgBroadeningGamma(const vector<Vector3>& fsmomenta, FourMomentum pgamma) {
    // Make a vector of the three-momenta in the final state
    double momentumSum(0.0);
    foreach (const Vector3& p3, fsmomenta) {
      momentumSum += mod(p3);
    }
    getLog() << Log::DEBUG << "Number of particles = " << fsmomenta.size() << endl;
    
    
    // Clear the caches
    _broadeninggammas.clear();
    _broadeninggammaAxes.clear();
    
    
    // If there are fewer than 2 visible particles, we can't do much
    if (fsmomenta.size() < 2) {
      for (int i = 0; i < 3; ++i) {
	_broadeninggammas.push_back(-1);
	_broadeninggammaAxes.push_back(Vector3(0,0,0));
      }
      return;
    }
    
    
    // Handle special case of thrustgamma = 1 if there are only 2 particles
    if (fsmomenta.size() == 2) {
      Vector3 axis(0,0,0);
      _broadeninggammas.push_back(1.0);
      _broadeninggammas.push_back(0.0);
      _broadeninggammas.push_back(0.0);
      axis = fsmomenta[0].unit();
      if (axis.z() < 0) axis = -axis;
      _broadeninggammaAxes.push_back(axis);
      /// @todo Improve this --- special directions bad...
      /// (a,b,c) _|_ 1/(a^2+b^2) (b,-a,0) etc., but which combination minimises error?
      if (axis.z() < 0.75)
	_broadeninggammaAxes.push_back( (axis.cross(Vector3(0,0,1))).unit() );
       else
         _broadeninggammaAxes.push_back( (axis.cross(Vector3(0,1,0))).unit() );
       _broadeninggammaAxes.push_back( _broadeninggammaAxes[0].cross(_broadeninggammaAxes[1]) );
      return;
    }
    
    
    
    // Temporary variables for calbs
    Vector3 axis(0,0,0);
    double val = 0.;
    
    // Get thrustgamma
    _calcBG(fsmomenta, val, axis, pgamma);
    getLog() << Log::DEBUG << "Mom sum = " << momentumSum << endl;
    _broadeninggammas.push_back(val / momentumSum);
    // Make sure that thrustgamma always points along the +ve z-axis.
    if (axis.z() < 0) axis = -axis;
    axis = axis.unit();
    getLog() << Log::DEBUG << "Axis = " << axis << endl;
    _broadeninggammaAxes.push_back(axis);
    
    // Get thrustgamma major
    vector<Vector3> threeMomenta;
    foreach (const Vector3& v, fsmomenta) {
      // Get the part of each 3-momentum which is perpendicular to the thrustgamma axis
      const Vector3 vpar = dot(v, axis.unit()) * axis.unit();
      threeMomenta.push_back(v - vpar);
    }
    _calcBG(threeMomenta, val, axis, pgamma);
    _broadeninggammas.push_back(val / momentumSum);
    if (axis.x() < 0) axis = -axis;
    axis = axis.unit();
    _broadeninggammaAxes.push_back(axis);
    
    // Get thrustgamma minor
    if (_broadeninggammaAxes[0].dot(_broadeninggammaAxes[1]) < 1e-10) {
      axis = _broadeninggammaAxes[0].cross(_broadeninggammaAxes[1]);
      _broadeninggammaAxes.push_back(axis);
      val = 0.0;
      foreach (const Vector3& v, fsmomenta) {
	val += fabs(dot(axis, v));
      }
      _broadeninggammas.push_back(val / momentumSum);
    } else {
      _broadeninggammas.push_back(-1.0);
      _broadeninggammaAxes.push_back(Vector3(0,0,0));
    }
    
  }
  
  
  
}