#ifndef __JCOMPAREHISTOGRAMS__JTESTKOLMOGOROV_T__
#define __JCOMPAREHISTOGRAMS__JTESTKOLMOGOROV_T__

#include <istream>
#include <ostream>

#include "JCompareHistograms/JTest_t.hh"
#include "JCompareHistograms/JResultTitle.hh"


/**
 * \author rgruiz
 */
namespace JCOMPAREHISTOGRAMS {

  using JGIZMO::JRootObjectID;
  
  /**
   * Implementation of the different Kolmogorov-related tests.
   */
  class JTestKolmogorov_t 
  { 
  public:  
    /**
     * Default constructor.
     */    
    JTestKolmogorov_t(){}

    /**
     * Kolmogorov test for 1D histograms.\n
     * The input parameter threshold, is used to evaluate whether the test is passed or failed.\n
     * The evaluation is done by comparing the threshold value with the result of the Kolmogorov test. The output of a Kolmogorov test is a p-value.\n 
     * The parameter threshold should therefore be a real value between 0 and 1.
     * 
     * \param h1                First histogram
     * \param h2                Second historgram
     * \param threshold         p-value
     * \param parameterName     Name of the parameter used to test the histograms
     * \param testName          Name of the test used to compare the histograms
     *
     * \return Test result 
     */
    JTestResult JKolmogorovTest(TH1* h1, TH1* h2, double threshold, std::string testName, std::string parameterName) {    

      using namespace std;
      using namespace JPP;
      
      int n1 = h1 -> GetNbinsX();
      int n2 = h2 -> GetNbinsX();

      if(n1 != n2)
	ERROR("Histograms with different bining. The objects: " << h1 -> GetName() << " can not be compared." << endl);
     
      TH2D* h3 = (TH2D*)h1->Clone(h1->GetName()==h2->GetName() ?
				  MAKE_CSTRING(to_string(h1->GetName())) :
				  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));

      bool afunc1 = false;
      bool afunc2 = false;
      double s1 = 1./h1->Integral();
      double s2 = 1./h2->Integral();

      double ew1, ew2, w1 = 0, w2 = 0;

      for (int bin = 1; bin <= n1; ++bin){
	ew1 = h1->GetBinError(bin);
	ew2 = h2->GetBinError(bin);
	w1 += ew1*ew1;
	w2 += ew2*ew2;
      }
    
      double esum1 = 0, esum2 = 0;
      if (w1 > 0)
	esum1 = 1./s1/s1/w1;
      else
	afunc1 = true;
      if (w2 > 0)
	esum2 = 1./s2/s2/w2;
      else
	afunc2 = true;
      if (afunc2 && afunc1) {
	ERROR("Errors are zero for both histograms\n");
      }
 
      double c1 = 0, c2 = 0;
 
      double dmax = 0;
      for (int bin=1 ; bin<=n1 ; ++bin){
	c1 += s1*h1->GetBinContent(bin);
	c2 += s2*h2->GetBinContent(bin);
	double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
	double p = TMath::KolmogorovProb(d);
	h3->SetBinContent(bin,p);
	dmax = TMath::Max(dmax,TMath::Abs(c1-c2));
      }
 
      double z;
    
      if  (afunc1)
	z = dmax*TMath::Sqrt(esum2);
      else if (afunc2)
	z = dmax*TMath::Sqrt(esum1);
      else
	z = dmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
    
      double pValue = TMath::KolmogorovProb(z);

      bool passed;
    
      (pValue < threshold ? passed = false : passed = true);
    
      JResultTitle title(testName, parameterName, passed , pValue);
    
      h3->SetTitle(title.getTitle().c_str());
    
      JTestResult r (testName,
		     JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
		     JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
		     parameterName, pValue, threshold, h3, passed);

      return r;
    };


    /**
     * Kolmogorov test for sliced 2D histograms.\n   
     * The histograms are sliced along the axis specified by the slice parameter. A slice per bin is made.\n
     * For each of the slices, the input parameter threshold is used to evaluate whether the test is passed or failed.\n
     * The evaluation is done by comparing the threshold value with the result of the Kolmogorov test. The output of the Kolmogorov test is a p-value.\n 
     * The parameter threshold should therefore be a real value between 0 and 1.
     * The fraction of failed tests is compared to the input parameter failuresThreshold. If this fraction is larger than failuresThreshold, the test fails. 
     *
     * \param h1 First histogram
     * \param h2 Second histogram
     * \param threshold Threshold value for the test result
     * \param failuresThreshold Threshold value for the test result
     * \param parameterName Name of the parameter used to test the histograms
     * \param testName Name of the test used to compare the histograms
     * \param slice The axis along which the histogram is sliced.
     *
     * \return Test result.
     */
    JTestResult JKolmogorovTestSlice(TH2* h1, TH2* h2, double threshold, double failuresThreshold, std::string testName, std::string parameterName, char slice) {

      using namespace std;
      using namespace JPP;

      int nFailures = 0;

      JTestResult r;
    
      if(slice == 'x' || slice == 'X'){
    
	int nSlices1 = h1->GetNbinsX();
	int nSlices2 = h2->GetNbinsX();

	TH1* h3 = h1->ProjectionX(MAKE_CSTRING(h1->GetName() << "_VS_" <<
					       h2->GetName() << "_KolmogorovTestSliceX"));
      
	if(nSlices1 != nSlices2) {
	  ERROR("Histograms with different binning. The objects: " << h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
	}

	for (int i=1 ; i<=nSlices1 ; ++i){

	  std::string sliceName = MAKE_STRING(h3->GetName() << "_" << to_string(i));
	     
	  TH1D* s1 = h1->ProjectionY (sliceName.c_str(),i,i);
	  TH1D* s2 = h2->ProjectionY (sliceName.c_str(),i,i);

	  if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; }	  
	
	  double p = s1 -> KolmogorovTest (s2);

	  bool passed = (p < threshold ? false : true);

	  if (!passed) nFailures++;
	
	  h3->SetBinContent(i,p);
	
	}

	bool passed = (nFailures/nSlices1 > failuresThreshold ? false : true);
      
	JResultTitle title(testName, parameterName, passed, nFailures);
      
	h3->SetTitle(title.getTitle().c_str());
      
	r = JTestResult  (testName,
			  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
			  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
			  parameterName, nFailures, failuresThreshold, h3, passed);
            
      } else if (slice == 'y' || slice == 'Y') {
	   
	int nSlices1 = h1->GetNbinsY();
	int nSlices2 = h2->GetNbinsY();

	TH1* h3 = h1->ProjectionY(MAKE_CSTRING(h1->GetName() << "_VS_" <<
					       h2->GetName() << "_KolmogorovTestSliceY"));
      
	if(nSlices1 != nSlices2) {
	  ERROR("Histograms with different binning. The objects: " << h1->GetName() << " can not be compared." << endl);
	}

	for (int i=1 ; i<=nSlices1 ; ++i){
	     
	  std::string sliceName = MAKE_STRING(h3->GetName() << "_" << to_string(i));

	  TH1D* s1 = h1->ProjectionX (sliceName.c_str(),i,i);
	  TH1D* s2 = h2->ProjectionX (sliceName.c_str(),i,i);

	  if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; }

	  double p = s1 -> KolmogorovTest (s2);

	  bool passed = (p < threshold ? false : true);

	  if (!passed) nFailures++;
	
	  h3->SetBinContent(i,p);
	}
      
	bool passed = (nFailures/nSlices1 > failuresThreshold ? false : true);
      
	JResultTitle title(testName, parameterName, passed, nFailures);
      
	h3->SetTitle(title.getTitle().c_str());
      
	r = JTestResult (testName,
			 JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
			 JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
			 parameterName, nFailures, failuresThreshold, h3, passed);
      
      }

      return r;
    };


  
    /**
     * Kolmogorov test for 2D histograms.\n
     * The input parameter threshold, is used to evaluate whether the test is passed or failed.\n
     * The evaluation is done by comparing the threshold value with the result of the Kolmogorov test. The output of a Kolmogorov test is a p-value.\n 
     * The parameter threshold should therefore be a real value between 0 and 1.
     * 
     * \param h1 First histogram
     * \param h2 Second historgram
     * \param threshold p-value
     * \param parameterName Name of the parameter used to test the histograms
     * \param testName Name of the test used to compare the histograms
     *
     * \return Test result
     */
    JTestResult JKolmogorovTest2D(TH2* h1, TH2* h2, double threshold, std::string testName, std::string parameterName) {    

      using namespace std;
      using namespace JPP;
    
      int n1x = h1 -> GetNbinsX();
      int n2x = h2 -> GetNbinsX();
      int n1y = h1 -> GetNbinsY();
      int n2y = h2 -> GetNbinsY();

      if(n1x != n2x || n1y != n2y)
	ERROR("Histograms with different bining. The objects: " << h1 -> GetName() << " can not be compared." << endl);
     
      if(h1->Integral()==0 || h2->Integral()==0)
	ERROR("Empty histogram: " << h1 -> GetName() << " can not be compared." << endl);

      double s1 = 1./h1->Integral();
      double s2 = 1./h2->Integral();

      TH2D* h3 = (TH2D*)h1->Clone(h1->GetName()==h2->GetName() ?
				  MAKE_CSTRING(to_string(h1->GetName())) :
				  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));

      bool afunc1 = false;
      bool afunc2 = false;

      double ew1, ew2, w1 = 0, w2 = 0;

      for (int i = 1; i <= n1x; ++i){
	for (int j = 1; j <= n1y; ++j){
	  ew1 = h1->GetBinError(i,j);
	  ew2 = h2->GetBinError(i,j);
	  w1 += ew1*ew1;
	  w2 += ew2*ew2;
	}
      }
    
      double esum1 = 0, esum2 = 0;
      if (w1 > 0)
	esum1 = 1./s1/s1/w1;
      else
	afunc1 = true;
      if (w2 > 0)
	esum2 = 1./s2/s2/w2;
      else
	afunc2 = true;
      if (afunc2 && afunc1) {
	ERROR("Errors are zero for both histograms\n");
      }
 
      double c1 = 0, c2 = 0;

      double dmax1 = 0;
      for (int i=1 ; i<=n1x ; ++i){
	for (int j=1 ; j<=n1y ; ++j){
	  c1 += s1*h1->GetBinContent(i,j);
	  c2 += s2*h2->GetBinContent(i,j);
	  double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
	  h3->Fill(i,j,d);
	  dmax1 = TMath::Max(dmax1,TMath::Abs(c1-c2));
	}
      }

      c1 = 0, c2 = 0;
    
      double dmax2 = 0;
      for (int j=1 ; j<=n1y ; ++j){
	for (int i=1 ; i<=n1x ; ++i){
	  c1 += s1*h1->GetBinContent(i,j);
	  c2 += s2*h2->GetBinContent(i,j);
	  double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
	  h3->Fill(i,j,d);
	  dmax1 = TMath::Max(dmax2,TMath::Abs(c1-c2));
	}
      }
        
      double dmax = 0.5*(dmax1+dmax2);
      double z;

      if  (afunc1)
	z = dmax*TMath::Sqrt(esum2);
      else if (afunc2)
	z = dmax*TMath::Sqrt(esum1);
      else
	z = dmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
    
      double pValue = TMath::KolmogorovProb(z);

      for (int i=1 ; i<=n1x ; ++i){
	for (int j=1 ; j<=n1y ; ++j){
	  h3->SetBinContent(i,j,TMath::KolmogorovProb(0.5 * h3->GetBinContent(i,j)));
	}
      }

      bool passed;
    
      (pValue < threshold ? passed = false : passed = true);
    
      JResultTitle title(testName, parameterName, passed , pValue);
    
      h3->SetTitle(title.getTitle().c_str());
    
      JTestResult r (testName,
		     JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
		     JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
		     parameterName, pValue, threshold, h3, passed);

      return r;
    };
  };
}
  
#endif