#ifndef __JCOMPAREHISTOGRAMS__JTESTRUNS_T__
#define __JCOMPAREHISTOGRAMS__JTESTRUNS_T__

#include <iostream>

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


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

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

    /**
     * Implements the Wald-Wolfowitx runs test: <https://en.wikipedia.org/wiki/Wald%E2%80%93Wolfowitz_runs_test>
     * In this, an expected number of runs and a standard deviation are computed from the number of bins and the number of "aboves" and "belows".
     * The test returns the difference between the observed number of runs and the expected number of runs, expressed in standard deviations. 
     * This is compared to the threshold input parameter.  
     * 
     * \param h1 First histogram
     * \param h2 Second histogram
     * \param threshold
     * \param parameterName Name of the parameter used to test the histograms
     * \param testName Name of the test used to compare the histograms
     */
    JTestResult JRunsTest(TH1* h1, TH1* h2, double threshold, std::string testName, std::string parameterName) {    

      using namespace std;
      using namespace JPP;
    
      if(h1 -> GetNbinsX() != h2 -> GetNbinsX())
	ERROR("Histograms with different bining. The objects: " << h1 -> GetName() << " can not be compared." << endl);

      double R = h1->Integral();
      double T = 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())));
    
      for (int i=1 ; i<h1->GetNbinsX() ; ++i) {
	h3->SetBinContent(i , (T/R)*h1->GetBinContent(i) - h2->GetBinContent(i)); 
      }
    
      int    n = 1;
      double p = 0;
      double q = 0;
    
      bool a = ((T/R)*h1->GetBinContent(1) - h2->GetBinContent(1)) < 0;

      (a ? p++ : q++);

      for (int i = 2 ; i<h1->GetNbinsX() ; ++i){
      
	bool b = ((T/R)*h1->GetBinContent(i) - h2->GetBinContent(i)) < 0;

	(b ? p++ : q++);

	if (b != a){
	  n++;
	  a=b;
	}
      }

      const double N = 1 + 2*p*q/(p+q) ;
      const double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) );
      const double d = (n-N)/s;
    
      const bool passed = (fabs(d) > threshold ? false : true);

      JResultTitle title(testName, parameterName , passed , fabs(d));
    
      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, fabs(d), threshold, h3, passed);
    
      return r;
    };

    /**
     * Runs 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 Runs test.\n 
     * 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 fraction of failed tests.
     * \param testName Name of the test used to compare the histograms
     * \param parameterName Name of the parameter used to test the histograms
     * \param slice The axis along which the histogram is sliced.
     *
     * \return Test result.
     */
    JTestResult JRunsTestSlice(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() << "_RunsTestSliceX"));
      
	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 R = s1->Integral();
	  double T = s2->Integral();
    
	  int    n = 1;
	  double p = 0;
	  double q = 0;
    
	  bool a = ((T/R)*s1->GetBinContent(1) - s2->GetBinContent(1)) < 0;

	  (a ? p++ : q++);

	  for (int i = 2 ; i<s1->GetNbinsX() ; ++i){
      
	    bool b = ((T/R)*s1->GetBinContent(i) - s2->GetBinContent(i)) < 0;

	    (b ? p++ : q++);

	    if (b != a){
	      n++;
	      a=b;
	    }
	  }

	  double N = 1 + 2*p*q/(p+q) ;
	  double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) );
	  double d = (n-N)/s;

	  bool passed = (fabs(d) > threshold ? false : true);

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

	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->GetNbinsX();
	int nSlices2 = h2->GetNbinsX();

	TH1* h3 = h1->ProjectionY(h1->GetName()==h2->GetName() ?
				  MAKE_CSTRING(to_string(h1->GetName())) :
				  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));
      
	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->ProjectionX (sliceName.c_str(),i,i);
	  TH1D* s2 = h2->ProjectionX (sliceName.c_str(),i,i);

	  if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; }
		
	  double R = s1->Integral();
	  double T = s2->Integral();
    
	  int    n = 1;
	  double p = 0;
	  double q = 0;
    
	  bool a = ((T/R)*s1->GetBinContent(1) - s2->GetBinContent(1)) < 0;

	  (a ? p++ : q++);

	  for (int i = 2 ; i<s1->GetNbinsX() ; ++i){
      
	    bool b = ((T/R)*s1->GetBinContent(i) - s2->GetBinContent(i)) < 0;

	    (b ? p++ : q++);

	    if (b != a){
	      n++;
	      a=b;
	    }
	  }

	  double N = 1 + 2*p*q/(p+q) ;
	  double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) );
	  double d = (n-N)/s;

	  bool passed = (fabs(d) > threshold ? false : true);

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

	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;
    };  
  };
}
  
#endif