#include <iostream>
#include <fstream>
#include "libstat.h"
#include "binning2.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TString.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"
#include "THStack.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TLatex.h"
#include "TMath.h"
#include "TGraphAsymmErrors.h"
#include "TLine.h"
#include "TString.h"
#include "TKey.h"
#include "TColor.h"
#include <sys/stat.h>
#include <TError.h>
#include <map>
#include <string>
#include <sstream>
//#include <fstream>
//#include <ifstream>
//#include <vector>
using namespace std;

//base class for the background
/**Base class for the background **/
class background : public TObject
{
 public:
  background(){};
  /** Standard constructor for the background class**/
  background(TString Name, Color_t Color, TFile* File){
    bkgName = Name; 
    bkgFile = File;
    bkgColor = Color;
  };
  ~background(){};
  TString bkgName; /**< Name of background **/
  TFile* bkgFile; /**< Pointer to the backgroundfile **/
  Color_t bkgColor; /**< bkgColor is no longer used*/
  TH1F* bkgHisto; /**< HIstogram containing the backgroundhistogram taken from TFile* bkgFile **/
/**Class definition, standard**/
  ClassDef(background, 1); 
};


//base class for systematics, only a name is associated
/**Base class for systematics**/
class systematic : public TObject
{
 public:
  systematic(){};
  systematic(TString Name){
    sysName = Name;
  };
  virtual ~systematic(){};
  TString sysName; /**< String containing the systematics name **/
  TH1F* hSysTotal; /**<Histogram containint the total systematics **/
  TH1F* hSysToy; /**<Histogram containing the systematics specific for this toy **/
/**Class definition, standard**/
  ClassDef(systematic, 1);
};

//class for the systematic that is applied to all bins of the distribution
//Needs only a value, since it applies to all
/** Class for a systematic that is applied to all bins of the distribution**/
class OverallSys : public systematic
{
 public:
 OverallSys(): systematic() { };
 /** Default constructor for OveralSys class**/
 OverallSys(TString Name, double Value): systematic() {
    sysName = Name;
    sysValue = Value;
  };
  virtual ~OverallSys(){};
 double sysValue; /**< Scaling factor for the systematic**/
     /** set Scaling factor which accounts how much of the systematic is added to the toy ***/
 void setAlpha(double alpha){sysAlpha = alpha;}; 
 double sysAlpha; /**<Addiional scaling factor for the toy deviations **/
/** CLass definition **/
 ClassDef(OverallSys, 1);
};
//class for detector systematics. Taken from variaton of cp tools
//applied to all backgrounds, gets filled from histograms in 
//void getSystematicHistogram(TString plotName, ExperimentalSys* myExpSys);
/**Class for detector systematics **/
class ExperimentalSys : public systematic
{
 public:
 ExperimentalSys(): systematic() { };
 /** Default constructor for class ExperimentalSys**/
 ExperimentalSys(TString Name, TString Type, TString Suffix): systematic() {
    sysName = Name;
    sysType = Type;
    sysSuffix = Suffix;
    sysAlpha = 0;
  };
  virtual ~ExperimentalSys(){};
    /** set Scaling factor which accounts how much of the systematic is added to the toy ***/
  void setAlpha(double alpha){sysAlpha = alpha;}; 
  TString sysType; /**< String giving the systematic type **/
  TString sysSuffix; /**< String given the suffix of the given systematic **/
  double sysAlpha; /**< Scaling factor for the systematic **/
  TH1F* hSysUp; /**< Systematics for the upper range of the bins **/
  TH1F* hSysDown; /**< Systematics for the lower range of the bins **/
/** Class definition **/
  ClassDef(ExperimentalSys, 1);
};

//class for systematics that is applied only on one background
//and is correlated in all bins and distributions
/**Class for systematic that is applied on only one background
 * and is correlated in all bins and distributions**/
class bkgNormSys : public systematic
{
 public:
 bkgNormSys(): systematic() { };
 /**Normal constructor for bkgNormSys class **/
 bkgNormSys(TString Name, TString bkgName, double Value): systematic() {
    sysName = Name;
    sysValue = Value;
    sysBkgName = bkgName;
  };
  virtual ~bkgNormSys(){}
  /** set Scaling factor which accounts how much of the systematic is added to the toy ***/
  void setAlpha(double alpha){sysAlpha = alpha;};
  TString sysBkgName;/**< Name of background systematic**/
  double sysValue; /**<Scaling factor for the systematic AND te toy **/
  double sysAlpha; /**< Scaling factor which accounts how much of the systematic is added to the toy **/
/** Class definition **/
  ClassDef(bkgNormSys, 1);
};
//class for systematics that is applied only on one background
//and is fully uncorrelated in all bins and distributions
/** Class for systematics that is applied on only one background and is fully uncorrelated in all bins and distributions **/
class bkgShapeSys : public systematic
{
 public:
 bkgShapeSys(): systematic() { };
 /** Default constructor for the bkgShapeSys class **/
 bkgShapeSys(TString Name, TString bkgName, double Value): systematic() {
    sysName = Name;
    sysValue = Value;
    sysBkgName = bkgName;
  };
  virtual ~bkgShapeSys(){};
  TString sysBkgName; /**< Name of background systematic**/
  double sysValue; /**< Scaling factor for the systematic**/
/** Class definition**/
  ClassDef(bkgShapeSys, 1);
};

//Same ad bkgShapeSys, but with histograms as input
/**Class for systeematics that is applied on only one background ans is fully uncorrelated in all bins and distributions
 * However, values are taken from histogram as input */
class bkgShapeHistoSys : public systematic
{
 public:
 bkgShapeHistoSys(): systematic() { };
 /**Default constructor for the bkgShapeHistoSys class **/
 bkgShapeHistoSys(TString Name, TString bkgName,TString Type, TString Suffix): systematic() {
    sysName = Name;
    sysBkgName = bkgName;
    sysType=Type;
    sysSuffix = Suffix;
  };
  virtual ~bkgShapeHistoSys(){}
  TString sysBkgName; /**< Name of background systematic**/
  TString sysType; /**< String containing the systematic type **/
  TString sysSuffix;/**< String containing the suffix of the systematic type **/
/**Class definition **/
  ClassDef(bkgShapeHistoSys, 1);
};

//class for uncertainities on a single background
//Values taken from histogram
/**Class for uncertainities on a single background, wih values taken from histogram*/
class bkgHistoSys : public systematic
{
	public:
	bkgHistoSys():systematic(){};
	 /**Normal constructor for bkgHistoSys class **/
	bkgHistoSys(TString Name, TString Type, TString bkgName,  TString Suffix):systematic(){
		sysName = Name; //Of both, uncertain if needed. Will see when implementation ismade
		sysSuffix=Suffix;
		sysBkgName =bkgName;
		sysType=Type;
		};
	virtual ~bkgHistoSys(){};
	TString sysBkgName; /**< Name of background systematic**/
	/** set Scaling factor which accounts how much of the systematic is added to the toy **/
	void setAlpha(double alpha){sysAlpha = alpha;};
	double sysAlpha; /**< Scaling factor which accounts how much of the systematic is added to the toy **/
	TString sysType; /**< String containing the systematic type **/
	TString sysSuffix; /**< String containing the suffix of the systematic type **/
	TH1F* hSysUp; /**< Histogram containing the systematics up. **/ 
	//Needed to split the systematic up errors from the systematic down errors
	TH1F* hSysDown; /**< Histogram containing the systematics down **/
	/** Class definition**/
	ClassDef(bkgHistoSys,1);
	
};
/** Class which generates toys and calculates the p_values. */
/** Look at runtoystest. C for an example application */
class toygenerator
{
 public:
  /**Initializes the toygenerator object given a certain seed.*/
  toygenerator(double seed){
    rnd = new TRandom3(seed);
    topoList = new TList();
    bkgList = new TList();
    sysList = new TList();
    sysExpList = new TList();
    sysBkgNormList = new TList();
    sysBkgShapeList = new TList();
    sysOverallList = new TList();
    sysBkgHistoList = new TList();
    //    stat = new StatUtils(seed+1345, 2, true);
    stat = new StatUtils(0, 2, true); /**<Initializaton of Statutils class**/
    hPvalues = new TH1F("hPvalues", "hPvalues", 100, 0, 10); /**< Initizaliations of Hiistogram for p values**/
    hBin4Content = new TH1F("hBin4Content", "hBin4Content", 5000,0,5000);
    countertopo=0;
    ApplySignal = 0;
    doPoissonVariation=true; /**< Without specification, poisson variation is on */
    doSelect_on_p_value=false; /**< p-value selection is standard off */
	  doBinning_restriction=false; /**< No binning restriction is done by default. */
    generateRegionfile=false; /**< No region file is generated by default */
    MeVToGeV=false;/**< No conversion to GeV range is done by default. Take care of turning it on by default**/
    allowMCdeficit=true; /**< In general, all data for toys is allowd **/
  };

  toygenerator(){
    rnd = new TRandom3(0);
    topoList = new TList();
    bkgList = new TList();
    sysList = new TList();
    sysExpList = new TList();
    sysBkgNormList = new TList();
    sysBkgShapeList = new TList();
    sysOverallList = new TList();
    sysBkgHistoList = new TList();
    stat = new StatUtils(0, 2, true);
    hPvalues = new TH1F("hPvalues", "hPvalues", 100, 0, 10);
    hBin4Content = new TH1F("hBin4Content", "hBin4Content", 5000,0,5000);
    countertopo=0;
    ApplySignal = 0;
    doPoissonVariation=true;
    doSelect_on_p_value=false;
    doBinning_restriction=false;
    doDynamic_binning=false;
	generateRegionfile=false;
	MeVToGeV=false;
	allowMCdeficit=true;
  };
  virtual ~toygenerator(){}

  void GenerateToys( int ntoys, TString filesuffix);
  void GenerateToyExperiment(TString plotname, TString filesuffix);

  void setExpSystematic(TString Name, TString Type, TString Suffix);
  void setThSystematic(TString Name, TString Type, TString bkgName, TString Suffix);
//  void setBkgNormSystematic(TString Name, TString bkg, double value);
//  void setBkgShapeSystematic(TString Name, TString bkg, double value);
  void setOverallSystematic(TString Name, double value);

  void setLumiData(double lumi=20800);
  void setTotalMcFile(TFile* totalMcFile);
  void setSignal(TString sigName, Color_t sigColor, TFile* sigFile);
  void setBackground(TString bkgName, Color_t bkgColor, TFile* bkgFile); //Ads a background. Should ask Sascha what it means phyisicaly
  void setBackground(TString bkgName, Color_t bkgColor, TFile* bkgFile, double NormSysValue);
  void setBackground(TString bkgName, Color_t bkgColor, TFile* bkgFile, double NormSysValue, double ShapeSysValue);
  void setResolutionBinning(bool dobinning);
  void setDirectory(TString inDirectory = "/lustre/users/sa1010/runs/rel17/runFinal/runMuon" );
  void setOutputDirectory(TString outDirectory = "/lustre/users/sa1010/runs/rel17/runFinal/runMuon" );
  void setAnalysisName(TString myAnalysis = "SD3PDSpecial_GeneralSearch" );
  void setOutputFileName(TString myFileName);
  void setDataFile(TFile* dataFile);
  void findDataTopologies(TString plotname, TString filename);
  void setStream(TString myStream = "Muon");
  void setPoissonVariation(bool PoissonVariation);
  void p_value_select(double max_p_value);
  void get_bin_restrictions(TString filename);
  void set_region_file(TString set_file_name);
  void	setDynamicBinning(bool dynamic);
  void setallowMCdeficit(bool deficit);
   
 private:
  TRandom3* rnd; /**< Our random number generator**/
  TH1F* emptyHistogram(TH1F* h); 
  void getSystematicHistogram(TString plotName, ExperimentalSys* myExpSys); //function gets systematic error from histogram
  void getSystematicHistogram(TString plotname, bkgHistoSys* mybkgHistoSys);// Overloading the function, since the first name isn't that practical 

  void getHistograms(TString plotName);
  double setMinimum(TString channel);
  int resolutionRebin(TString channel);
  TH1F* makeErrorHistogram(TH1F* hMC_ERRUP, TH1F*hMC_ERRDOWN);
  bool MeVToGeV; /**< Boolian for switching the axis to GeV**/
  bool DoResolutionBinning; /**< Boolian to do resolution binning, applies only in the case of equidistant binning**/
  bool doPoissonVariation; /**< Boolian which dictates the variation of the toys **/
  bool doSelect_on_p_value; /**< Boolian that dictates if there is a p_value selection **/
  bool doBinning_restriction; /**< Boolian that dictates if there should be binning restriction **/
  bool doDynamic_binning; /**< Boolian that dictates if there is dynamic or no dynamic binning **/
  bool generateRegionfile; /**< Boolian that dictates the generation of a region file.**/
  bool allowMCdeficit; /**< Boolian that dictates if a MC higher than the data is allowed**/
  
  double lumiData; /**< Stores the integrated luminosity**/
  double p_selection_value; /**< Stores the p_value_selection threshold**/
  int rebin; /**<Stores the rebin factor for equidistant rebinning**/
  void axisMeVToGeV(TAxis& axis, TString plotname="0");
  //void axisMeVToGeV(TAxis& axis);
  void rebin_histogram(TString channel, TH1F* histogram);
  
  StatUtils* stat; /**<Object of class StatUtils**/
  TH1F* hPvalues; /**< Histogram for storing all Pvalus**/
  TH1F* hGlobalPvaluesPlot; /**<Histogram for storing the global P values, same as hPvalues**/
  TH1F* hBin4Content; /**< Stores the 4th bin content of every plot. Not used anymore**/
  int countertopo; /**< Counter for storing amount of topologies**/

  ofstream myfile; /**< Output file for storing calculations**/
  TList* sysList; /**<List with al systematics**/
  TList* sysExpList; /**<List containing all experimental systematics**/
  TList* sysOverallList; /**< List containing all overal systematics**/
  TList* sysBkgNormList; /**< List containing all bkg NOrm systematics**/
  TList* sysBkgShapeList; /**< List containing akk bkg shape systematics**/
  TList *sysBkgHistoList; /**< List containing all Backgrounds histogram systematics**/
  TList* topoList; /**< List that was used for all topo's. No longer in use **/
  TList* bkgList; /**< List containing all backgrounds **/
  vector<TString> topoVector; /**< Vector containing all topologies**/
  //The name here is really confusing, but since topoVector is already
 //taken.
  //vector<const double*> TopoArrayVector; /*<

  TH1F* hSIG; /**< Histogram for the signal**/
  TFile* fSIG; /**< File which contains the signal histograms**/ 
  bool ApplySignal; /**<Boolian that dictates if a signal is applies**/


  TH1F* hDATA; /**< HIstogram for the data of a specific channel**/
  TH1F* hMC; /**<HIstogram for the mc expectation**/
  TH1F* hMCtotal;/**< History for the mc expectation of all bkg's.**/
  TFile* fMCtotal;/**< File containing all MC data**/
  TFile* fDATA; /**<File containing all Data 2012. Barely used**/
  TFile* outputFile; /**< Pointer to file to which text output should be written**/
  ofstream outputTextFile; /**< Stream to output to outputFIle**/
  ofstream outputRegionFile; /**< Stream to output the region file**/

  TString inputDirectory; /**< String for the input directory. No longer in use**/
  TString outputDirectory; /**< String for the output directory**/
  TString outputFileName; /**< String for the filename of the output**/
  TString stream; /**<String that dictates which stream we are looking at, eg Muon, Egamma, Jettauetmiss**/
  TString analysis; /**< Name of the  analysis**/
	TString outputregionFileName; /**<Name of the ouputfile containing the regions**/

  map<TString,vector<int> > region_map; /**< Map for storing the region restrictions**/
  /** Standard class definition**/
  ClassDef(toygenerator, 1);
};