#ifndef ASIMOVHYPOTHESIS_INCLUDED_HH
#define ASIMOVHYPOTHESIS_INCLUDED_HH

#include <vector>
#include <string>
#include <TMath.h>
#include <map>
#include <iostream>
#include <TH2D.h>
#include <TRandom.h>
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"
#include "TMinuitMinimizer.h"
#include "TStopwatch.h"

#include "statutil.hh"
#include "stringutil.hh"  
#include "rootutil.hh"
using stringutil::print;

using namespace std;

/*! Class to fit simple mu* S + B models models to datasets, generate speudo-data 
    and compute distributions of poisson-likelihood-ratio-test-statisics */

struct BinnedAnalysis2D
{
	vector<TH2D> signal_histos;
	vector<TH2D> background_histos;

	int verb = 0;
	bool sparse = true;  //<! if true, use fast likelihood evaluation for sparse histogram
	double mumin = 1e-8; //<! minimum value for use of mu (to prevel log(p=0))


	BinnedAnalysis2D() {}

	void add_histos( const TH2D& hsig, const TH2D& hbg )
	{
		signal_histos.push_back( hsig );
		background_histos.push_back( hbg );
	}
	

	/*! make a histogram with pseudo-data for the given sigstrength */

	vector<TH2D> make_pseudo_dataset( double sigstrength )
	{
		vector<TH2D> r = signal_histos; // copy
		fill_pseudo_dataset( r, sigstrength );
		return r;
	}

	/*! fill histogram with pseudo-data for the given sigstrength */

	void fill_pseudo_dataset( vector<TH2D>& data, double sigstrength );
	

	/*! Log of the poisson likelihood of the dataset, given the signal strength */

	double log_likelihood( const vector<TH2D>& data , double sigstrength );


	/*! The log likelihood ratio of the data as function of the tested sig. strength */

	double log_likelihood_ratio( const vector<TH2D>& data, double sigstrength );
	

	/*! Fit the signal-strength, optimising the log-likelihood */

	double optimize_sigstrength( const vector<TH2D>& data , 
	                             double startvalue = 1, double xmin = -3, double xmax = 1000, double stepsize = 1e-3);

	/*! Return the likelihood ratio test statistic for a data-set */

	double test_statistic ( const vector<TH2D>& data )
	{
		double s = optimize_sigstrength( data );
		double r = log_likelihood_ratio( data, s );
		if (verb) print( s, r );
		return r;
	}

	/*! Return a histogram with the distribution of the test-statistic for the given signal strength */

	TH1D test_statistic_distribution( double sig_strength , int n_pe =1000 , int nbins = 1100 , double xmin = -10 , double xmax = 100  );


	/*! produce histogram to check bias of signal-strength fit */

	TH2D check_signal_fit( int n_pe = 1000, double ss_min = 0, double ss_max = 10 ) 
	{
		TH2D r("hsigfit","hsigfit;true signal strength;fitted-true signal strength", 50, ss_min, ss_max, 31, -5, 5 );

		for (int i =0; i< n_pe; i++ )
		{
			double ss = gRandom->Uniform( ss_min, ss_max );
			auto data = make_pseudo_dataset( ss );
			double s = optimize_sigstrength( data );
			r.Fill( ss, (s-ss) );
		}

		return r;
	}

};


#endif