#ifndef STATUTILINCLUDED
#define STATUTILINCLUDED

#include "Math/ProbFunc.h" // gaussian cdf 
#include "Math/QuantFuncMathCore.h" // gaussian quantile
#include "TGraph.h"

/*! obtain p value from Z in std dev
 * slide 3 https://idpasc.lip.pt/uploads/talk/file/192/Stats-Tutorial2-2015.pdf
 * */
inline double pvalue_from_Z(double Z, bool onesided = true)
{
	if (onesided)
		return ROOT::Math::gaussian_cdf(-Z, 1, 0); // sigma 1, mean 0
	else
		return ROOT::Math::gaussian_cdf(-Z, 1, 0) * 2;
}

/*! obtain Z in std dev from p value */
inline double Z_from_pvalue(double p_value, bool onesided = true)
{
	if (onesided)
		return ROOT::Math::gaussian_quantile(1 - p_value, 1); // sigma 1, mean 0
	else
        return ROOT::Math::gaussian_quantile(1 - p_value / 2, 1);
}



/*! helper for dealing with which fraction of a bin falls in a range */

struct Interval
{
	double x1, x2;
	Interval(double a, double b) : x1(a), x2(b) {}

	/*! return the fraction of the range between x1 and x2 that is below cut */
	double fraction_below(double cut)
	{
		if (x2 <= x1) return -999; // error
		if (x2 <= cut) return 1; // whole range below cut
		if (x1 >= cut) return 0; // whole range above cut
		// else x1<cut<x2
		return (cut - x1) / (x2 - x1);
	}

	/*! return fraction of interval that is above cut */
	double fraction_above(double cut)
	{
		return 1 - fraction_below( cut);
	}

	/* return the fraction of the interval that is between cut1, and cut2 */
	double fraction_in_range(double cut1, double cut2)
	{
		return fraction_below(cut2) - fraction_below(cut1);
	}
};


/*! Get Median of a histogram */
inline double median( TH1& h )
{
	double prob[1] = {0.5};
	double q_median[1] = {0.0};
	h.GetQuantiles(1, q_median, prob);
	return q_median[0];
}


/*! Get Median of a histogram */
inline double fraction_below( TH1& h , double x , bool with_overflow = true )
{
	double c(0),d(0);

	for (int i = with_overflow?0:1 ; i<=h.GetNbinsX() + with_overflow?1:0 ; i++)
		{
			Interval ival( h.GetXaxis()->GetBinLowEdge(i), h.GetXaxis()->GetBinUpEdge(i) );
			c += ival.fraction_below( x ) * h.GetBinContent(i);
			d +=                            h.GetBinContent(i);
		}

	if (d==0) return 0;
	return c/d;
}

/*! return the first x-value where the graph intersects level. 
    graph must have a least 2 points. 
	linear interpolation. 
*/
inline double graph_intersect( TGraph& g , double level, double sentinel = 0 )
{
	for (int i = 1; i < g.GetN() ; i++ )
	{
		double x0 = g.GetX()[i-1];
		double y0=  g.GetY()[i-1];
		double x1 = g.GetX()[i];
		double y1=  g.GetY()[i];

		if ( y0 >= level && level > y1 ) // line going down, y0>level>y1
		{
			return x0 + (y0-level)/(y0-y1)*(x1-x0);
		}
		if ( y0 < level && level <= y1 ) // line going up, y0 < level < y1
		{
			return x0 + ( level-y0 )/ ( y1-y0 )*(x1-x0);
		}
	}
	
	// no intersection.
	return sentinel;

}

#endif