/// -*- C++ -*-
#include "Rivet/Rivet.hh"
#include "ThrustGamma.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 ThrustGamma::compare(const Projection& p) const {
    return mkNamedPCmp(p, "Kinematics");
  } 

  
  void ThrustGamma::calctg(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);
    }

    _calcThrustGamma(threeMomenta,pgamma);
  }
  
  void ThrustGamma::calctg(const vector<FourMomentum>& fsmomenta, FourMomentum pgamma) {
    vector<Vector3> threeMomenta;
    threeMomenta.reserve(fsmomenta.size());

    foreach (const FourMomentum& v, fsmomenta) {
      threeMomenta.push_back(v.vector3());
    }
    _calcThrustGamma(threeMomenta, pgamma);
  }
  
  void ThrustGamma::calctg(const vector<Vector3>& fsmomenta, FourMomentum pgamma) {
    _calcThrustGamma(fsmomenta, pgamma);
  }
  
  
  /////////////////////////////////////////////////
  
  
  
  inline bool mod2Cmp(const Vector3& a, const Vector3& b) {
    return a.mod2() > b.mod2();
  }
  
  
  // Do the general case thrustgamma calculation
  void _calcTG(const vector<Vector3>& momenta, double& t, Vector3& taxis, FourMomentum pgamma) {
   

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

      
    pgam = pgam.unit();


      // Calculate the thrustgamma value 
      t=0.;
      double foo;

      for (unsigned int k=0 ; k<p.size() ; k++){

	foo=pgam.dot(p[k]);    
      t += sqrt (foo * foo);

      }

      // Store everything
     taxis = pgam;
    
     printf("%s %f %f %f %s", "tgamma axis", pgam.unit().x(), pgam.unit().y(), pgam.unit().z(), "...\n"); 
 
      }
  

  
  
  
  // Do the full calbulation
  void ThrustGamma::_calcThrustGamma(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
    _thrustgammas.clear();
    _thrustgammaAxes.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) {

	_thrustgammas.push_back(-1);
	_thrustgammaAxes.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);
      _thrustgammas.push_back(1.0);
      _thrustgammas.push_back(0.0);
      _thrustgammas.push_back(0.0);
      axis = fsmomenta[0].unit();

      if (axis.z() < 0) axis = -axis;

      _thrustgammaAxes.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)
	_thrustgammaAxes.push_back( (axis.cross(Vector3(0,0,1))).unit() );
       else
         _thrustgammaAxes.push_back( (axis.cross(Vector3(0,1,0))).unit() );
       _thrustgammaAxes.push_back( _thrustgammaAxes[0].cross(_thrustgammaAxes[1]) );
      return;
    }
    
    
    
    // Temporary variables for calbs
    Vector3 axis(0,0,0);
    double val = 0.;
    
    // Get thrustgamma
    _calcTG(fsmomenta, val, axis, pgamma);
    getLog() << Log::DEBUG << "Mom sum = " << momentumSum << endl;
    _thrustgammas.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;
    _thrustgammaAxes.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);

    }

    _calcTG(threeMomenta, val, axis, pgamma);
    _thrustgammas.push_back(val / momentumSum);
    if (axis.x() < 0) axis = -axis;
    axis = axis.unit();
    _thrustgammaAxes.push_back(axis);
    
    // Get thrustgamma minor
    if (_thrustgammaAxes[0].dot(_thrustgammaAxes[1]) < 1e-10) {

      axis = _thrustgammaAxes[0].cross(_thrustgammaAxes[1]);
      _thrustgammaAxes.push_back(axis);
      val = 0.0;

      foreach (const Vector3& v, fsmomenta) {

	val += fabs(dot(axis, v));

      }
      _thrustgammas.push_back(val / momentumSum);

    } else {
      _thrustgammas.push_back(-1.0);
      _thrustgammaAxes.push_back(Vector3(0,0,0));
    }
    
  }
  
  
  
}