#include <cstdlib>
#include <iostream>
#include <fstream>
#include <sstream>
#include <map>
#include <string>
#include <vector>

#include "TMath.h"
#include "TChain.h"
#include "TFile.h"
#include "TTree.h"
#include "TLatex.h"
#include "TString.h"
#include "TObjString.h"
#include <TMinuit.h>
#include <TColor.h>
#include <limits.h>
#include <float.h>
#include "TSystem.h"
#include "TROOT.h"

#include "TMVA/Factory.h"
#include "TMVA/DataLoader.h"
#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/MethodCuts.h"
#include "TMVA/TMVAGui.h"

using namespace std;
using namespace TMVA;


/*--------------------------------------------------------------------------*
 *                 Significance Function for Minuit Analysis                *
 *--------------------------------------------------------------------------*/

// fcn passes back f = chi, the function to be minimized.
void fcn(int nvars, double* deriv, double f, double par[], int flag) {
  TString sigFilename = "hhttbb.root";
  TString bckFilename[] = {"ttbar_new.root", "Ztautau.root", "ZtautauB.root", "ZtautauC.root", "ZtautauL.root"};
  TString TupleName = "CxAODTuple_Nominal";
  const int nBackgrounds = 5;

// Get the TTree objects from files
  TFile* siginputfile = new TFile(sigFilename);
  TTree* sigTree = (TTree*)siginputfile->Get(TupleName);

// Get the TTree objects from files
  TFile* bckinputfile[nBackgrounds];
  TTree* bckTree[nBackgrounds];
  for ( int i=0 ; i<nBackgrounds ; i++ ) {
    bckinputfile[i] = new TFile(bckFilename[i]);
    bckTree[i] = (TTree*)bckinputfile[i]->Get(TupleName);
  }

  TString varsName[] = {"mBB", "mMMC","mHH", "drBB", "drLepTau", "mtLepMet"};
  double sf[] = {1e-3, 1., 1e-3, 1., 1., 1e-3};

  const int nSig = sigTree->GetEntries();
  int nBck[nBackgrounds];
  for ( int i=0 ; i<nBackgrounds ; i++ ) {
    nBck[i] = bckTree[i]->GetEntries();
  }

// set address of variable to cut over
  TString SVN = "SelVars_Nominal_";
  float var[nvars], weight;
  sigTree->SetBranchAddress("EventWeight_Nominal_TotalWeight", &weight);
  for ( int i=0 ; i<nvars ; i++ ) {
    sigTree->SetBranchAddress(SVN + varsName[i], &var[i]);
  }

  for ( int i=0 ; i<nBackgrounds ; i++ ) {
    bckTree[i]->SetBranchAddress("EventWeight_Nominal_TotalWeight", &weight);
    for ( int j=0 ; j<nvars ; j++ ) {
      bckTree[i]->SetBranchAddress(SVN + varsName[j], &var[j]);
    }
  }

  double cutVals[nvars][2];
  for ( int i=0 ; i<nvars ; i++ ) {
    cutVals[i][0] = par[i];
    cutVals[i][1] = par[i+nvars];
  }

// get the signal events for each cut
  double sigCount = 0.;
  double sigTotalWeight = 0.;

  for ( int i=0 ; i<nSig ; i++ ) {
    sigTree->GetEvent(i);
    bool pass = true;
    for ( int k=0 ; k<nvars ; k++ ) {
      if ( var[k] < cutVals[k][0]/sf[k] ) { pass = false;}
      if ( var[k] > cutVals[k][1]/sf[k] ) { pass = false;}
    }

    if ( pass == true ) { sigCount += weight; }

    sigTotalWeight += weight;
    }



// get the background events for each cut
// get the background events for each cut
  double bckCount[nBackgrounds];
  for( int i=0 ; i<nBackgrounds ; i++) { bckCount[i] = 0.;}

  double bckTotalWeight[nBackgrounds];
  for ( int i=0 ; i<nBackgrounds ; i++ ) {
    bckTotalWeight[i] = 0.;
    for ( int j=0 ; j<nBck[i] ; j++ ) {
      bckTree[i]->GetEvent(j);
      bool pass = true;
      for ( int l=0 ; l<nvars ; l++ ) {
        if ( var[l] < cutVals[l][0]/sf[l] ) { pass = false; }
        if ( var[l] > cutVals[l][1]/sf[l] ) { pass = false; }
      }

      if ( pass == true ) { bckCount[i] += weight; }

      bckTotalWeight[i] += weight;
    }
  }

// find the significance value
  double s = sigCount/sigTotalWeight;
  double b = 0.;
  for ( int j=0 ; j<nBackgrounds ; j++ ) {
    b += bckCount[j]/bckTotalWeight[j];
  }

  double Z;
  if ( b!=0. ) { Z = sqrt( 2. * ( (s+b)*TMath::Log(1 + (s/b)) - s ) );;}
  else {Z=0;}



  delete sigTree;
  delete siginputfile;
  for ( int i=0 ; i<nBackgrounds ; i++ ) {
    delete bckTree[i];
    delete bckinputfile[i];
  }

  f = -1*Z;

  cout << f << endl;

}

void Debugger() {

  double f;
  int nvars = 6;
  double par[nvars*2];
  double* deriv;
  int flag = 0;

/*--------------------------------------------------------------
  Cut Below			  Cut above
  --------------------------------------------------------------*/
// starting value
  par[0] = 72.;		  par[0+nvars] = 141.;		// mBB
  par[1] = 98.;		  par[1+nvars] = 476.;		// mMMC
  par[2] = 338.;	  par[2+nvars] = 3.93921e+06;	// mHH
  par[3] = 0.76;	  par[3+nvars] = 4.7;		// drBB
  par[4] = 0.6;		  par[4+nvars] = 4.9;		// drLepTau
  par[5] = 15.;	          par[5+nvars] = 190.;		// mtLepMet

  for ( int i=0 ; i<100 ; i++) {
    par[0] = 72.+(i*0.5);
    fcn(nvars, deriv, f, par, flag);
  }

}