/*
-----------------------------------------------
------Function for calculating sphericity------
-------------Ian Connelly 13/11/11-------------
-----------------------------------------------
This function takes as input a double array, A, defined as: pxx, pxy, pxz, pyx, pyy, pyz, pzx, pzy, pzz
This array is techinally symmetric but the function has been defined using generic matricies.
A 3x3 matrix is formed, and the eigen values are calculated, sorted and then used to return the value of sphericity.
*/

#ifndef SPHERICITY_H
#define SPHERICITY_H

#include "TMatrixD.h"
#include "TMatrixDEigen.h"
#include <algorithm>

double sphericity(double A[9]) {
  
  TMatrixD M(3,3); 
  M.SetMatrixArray(A);
  TMatrixDEigen V(M);
  TMatrixD Eval = V.GetEigenValues();
  
  double e1 = TMatrixDRow(Eval,0)(0);
  double e2 = TMatrixDRow(Eval,1)(1);
  double e3 = TMatrixDRow(Eval,2)(2);
  
  std::vector<double> evalvector;
  evalvector.push_back(e1);
  evalvector.push_back(e2);
  evalvector.push_back(e3);

  //sort eigenvalues to get the lowest two for use in sphericity (lowest to highest order)
  sort (evalvector.begin(), evalvector.end());

  //error-checking
  //this number should equal zero as the off-diagonal elements should all be zero in the eigenvalue matrix
  //returns error value of -1 
 
  double check = TMatrixDRow(Eval,0)(1) + TMatrixDRow(Eval,0)(2) + TMatrixDRow(Eval,1)(0) + TMatrixDRow(Eval,1)(2) + TMatrixDRow(Eval,2)(0) + TMatrixDRow(Eval,2)(1);
  if (check != 0.0) {double err = -1; return err;}    

  //for formula, see: http://cepa.fnal.gov/psm/simulation/mcgen/lund/pythia_manual/pythia6.3/pythia6301/node213.html

  double value = evalvector.at(0)+evalvector.at(1);
  double sphericity = 1.5*value;

  return sphericity;
}

double sphericity(double A[4], int a) { //include an integer if want transverse sphericity

  TMatrixD M(2,2);
  M.SetMatrixArray(A);
  TMatrixDEigen V(M);
  TMatrixD Eval = V.GetEigenValues();

  double e1 = TMatrixDRow(Eval,0)(0);
  double e2 = TMatrixDRow(Eval,1)(1);

  std::vector<double> evalvector;
  evalvector.push_back(e1);
  evalvector.push_back(e2);

  sort(evalvector.begin(), evalvector.end());

  //error-check
  double check = TMatrixDRow(Eval,0)(1) + TMatrixDRow(Eval,1)(0);
  if (check != 0.0){double err = -1; return err;}

  //from http://www.nikhef.nl/pub/experiments/atlas/daq/Dankers/chapters_04_05.pdf pg26
  // http://arxiv.org/pdf/1110.2278v1

  double transverse_sphericity = 2*evalvector.at(0)/(evalvector.at(0) + evalvector.at(1));

  return transverse_sphericity;
}

#endif