// -*- C++ -*-
#include "Rivet/Rivet.hh"
#include "Cparam.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Particle.hh"


namespace Rivet {
  
  
  void Cparam::calc(const FinalState& fs) {//take the final state object and use it to put a list of the particles in the event to the next calc
    calc(fs.particles());
  }
  
  void Cparam::calc(const vector<Particle>& fsparticles) {//take a list of the particles and extract the 4momentums of each in a list
    vector<FourMomentum> fourMomenta;
    fourMomenta.reserve(fsparticles.size());
    foreach (const Particle& p, fsparticles) {
      const FourMomentum p4 = p.momentum();
      fourMomenta.push_back(p4);
     }
    _calcCparam(fourMomenta);
  }
  
  void Cparam::calc(const vector<FourMomentum>& fsmomenta) {//take a list of the 4momentums of particles in an event and calculate from that the CParameter of the event
    _calcCparam(fsmomenta);
  }
   
  /////////////////////////////////////////////////
  
  
  
  inline bool mod2Cmp(const Vector3& a, const Vector3& b) {
    return a.mod2() > b.mod2();
  }
  
  
  // Do the general case cparam calculation
  void _calcP(const vector<FourMomentum>& momenta, double& m) {

    //create a vector storing 3 vector objects describing all the momementums of the final state particles in a given event:
    vector<Vector3> p;
    p.reserve(momenta.size());
    foreach (const FourMomentum& v, momenta) {
      p.push_back(v.vector3());
    }

    m=0.;
    double pn=0.;
    double pma=0.;
    //Here I loop over the set of particles twice to calculate the CParameter
    for (unsigned int k=0 ; k<p.size() ; k++){ 
      for (unsigned int j=0 ; j<p.size() ; j++){ 
	pma+=(mod((p[k])))*(mod((p[j])))*(sin(p[k].angle(p[j])))*(sin(p[k].angle(p[j])));
      }
      pn+=mod((p[k]));
    }
    // finally setting the value of m here equal to the C-Parameter to be put to a array in _calcparam
    m=(3*pma)/(2*(pn*pn));
  }       

  
  
  
  
  // Do the full calculation
  void Cparam::_calcCparam(const vector<FourMomentum>& fsmomenta) {
    
    // Clear the caches
    _cparams.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) {
	_cparams.push_back(-1);
      }
      return;
    }          
    

    // Temporary variables for calcps

    double val = 0.;
    
    // Get cparam
    _calcP(fsmomenta, val);

    //store the cparameter to be called by cparam():
    _cparams.push_back(val);

    
  }
  
  
  
}