/**
   @file ToolsLib.h
   @author Benoit Revenu
   @brief this file defines some functions used in the analysis
   (fits, distances...)
*/

#ifndef _TOOLSLIB_H_
#define _TOOLSLIB_H_

#include <vector>
#include <iostream>
#include <fstream>
#include <unistd.h>
#include <string>
#include <map>

#include "TCanvas.h"
#include "TROOT.h"
#include "TGraph.h"
#include "TFile.h"
#include "TH1.h"

#include "TSmSp.h"

using namespace std;

static const double kDref = 1000;
static const double kDmin = 300;

#define sqr(x) ((x)*(x))
#define maxval(a,b) ( (a)>(b) ? (a) : (b) )
#define minval(a,b) ( (a)<(b) ? (a) : (b) )

template<class T>
ostream& operator<<(ostream& s, const vector<T>& vc);
//template<class T> ostream& operator<<(ostream& s, const valarray<T>& vc);
template<class T>
ostream& operator<<(ostream& s, const map<T,int>& m);

/**
   @brief this function converts u, v into theta (rad), phi (rad)
*/
void uv2thetaphi(double u, double v, double& theta, double& phi);

/**
   @brief this function converts theta (rad), phi (rad) into u, v
*/
void thetaphi2uv(double theta, double phi, double& u, double& v);

/**
 * @brief this function computes the direction parameters u and v
 * @param listcoords size = number of stations, 
 *  listcoords[i].size = 3 (x,y,z of each station)
 * @param times size = number of stations and times[i] is the trigger time for Tank i
 * @param t0 is the time of the impact point on the ground
 * @param u is sin(theta) cos(phi)
 * @param v is sin(theta) sin(phi)
 */
void compute_uv(const vector< vector<double> > & listcoords,
                const vector<double> & times,
                double * t0, double *u, double *v);

/**
 * @brief this function computes the position of the core of the shower
 */
void compute_barycentre(const vector< vector<double> > & listcoords,
                        const vector<double> & sumpmt, const vector<double> listtimes,
                        double *xav, double *yav, double *zav, double *tav);

/**
 * @brief this function computes the distances of all the stations in the event to the shower axis
 */
void compute_distances(double theta, double phi, double ximp, double yimp, double zimp,
                       const vector< vector<double> > & listcoords,
                       vector<double> & listdist);

/**
 * @brief this function computes the distances of all the stations in the event to the shower axis
 */
void compute_distancesUV(double u, double v, double ximp, double yimp, double zimp,
                         const vector< vector<double> > & listcoords,
                         vector<double> & listdist);

/**
 * @brief this function compute the maximum distance between all pairs of stations in the event
 * @return the value of this maximum distance
 */
double compute_dmax(const vector< vector<double> >& listcoords);

/**
 * @brief this function returns the area of the triangle
 * must be used only with 3 tanks events !
 * @return the area of the triangle
 */
double compute_surface(const vector< vector<double> >& listcoords);

double compute_angleXY(double x1, double y1, double x2, double y2);

double max_angle_in_triangle(const vector<double>& a,
                             const vector<double>& b,
                             const vector<double>& c);

/**
 * @brief this function returns the distance between the shower axis
 * (intersecting the ground at (0,0,0) and the point 
 * (ximp,yimp,zimp)).
 * @param theta is the zenithal angle of the shower
 * @param phi is the azimuth of the shower
 * @param ximp is the x distance of the considered point to the core
 * @param yimp is the y distance of the considered point to the core
 * @param zimp is the z distance of the considered point to the core
 * @return the distance between the shower axis and the considered point
 */
double compute_dist(double theta, double phi, double ximp, double yimp, double zimp);

/**
 * @brief this function returns the distance between the shower axis
 * (intersecting the ground at (0,0,0) and the point 
 * (ximp,yimp,zimp)).
 * @param u is sin(theta)cos(phi) angle of the shower
 * @param v is sin(theta)sin(phi) of the shower
 * @param ximp is the x distance of the considered point to the core
 * @param yimp is the y distance of the considered point to the core
 * @param zimp is the z distance of the considered point to the core
 * @return the distance between the shower axis and the considered point
 */
double compute_distUV(double u, double v, double ximp, double yimp, double zimp);

double compute_LDFHPFunc(double dist, double theta);

double compute_LDFLOGFunc(double dist, double beta, double gamma);

void GetBetaGammaFromTheta(int primary, double theta, double& beta, double& gamma);

double compute_LDFPhotonFunc(double dist, double theta, double energy);
double compute_LDF2EXPOFunc(double dist, double theta, double energy, int primary=0);

double compute_energy(string ldftype, double SRef, double theta, int primary=0, double ds1000 = 0, double * de = 0x0);

double compute_S1000(string ldftype, double E, double theta, TSMSPEstim *apSmspEstim=0x0,int primary=0);

void compute_error_angle(double u,double du,double v,double dv,
                         double *dtheta,double *dphi);

void compute_error_energy(string ldftype, double SRef,double dSRef,double theta,double dtheta,
                          double *de);

void compute_dE_dS1000_2EXPO(double E, double theta, double SRef, double *de, double *ds1000);

void fit_ellipse(unsigned int npts, double *x, double *y, double *w, double& xm, double& ym,
                 double &sa, double &sb, double &theta, string method_orientation="distances",
                 string method_area="convexity");

void GetPointsOrientation(int npts, double * x, double * y, double * w, double & xm, double & ym,
                          double & a, double & b, double & sa, double & sb, double & theta,
                          string method="distances");

double ComputePolygonArea(int npts, double * x, double * y);

vector<int> get_edges(const vector<int>& lsids, const vector< vector<double> >& coords);

vector<int> get_edges(int npts, double * x, double * y);

// compute_asymetry : x and y are actually dx=x-xcore, dy=y-ycore
double compute_asymetry(double x, double y, double dist, double u, double v, int method=0, double e=0);

void intp2dim(double x, double y, double x1, double x2, double y1, double y2, double z11, double z12, double z21, double z22, double &z);

double lnStirling(int n);

double ProbPoisson(int k, double mu);

float a25(float r);
float a36(float r);
float a45(float r);
float a53(float r);
float a60(float r);
float a66(float r);
float a72(float r);

double distance2Core(double dx, double dy, double u, double v);
float interpol(float x1, float x2, float y1, float y2, float x);
int matinv4(double *a, double *b);

#endif