#include <iostream>
#include <sstream>
#include <iomanip>
#include "TCanvas.h"
#include "TH1F.h"


#define nrT 100

void readMark(bool saving=false)
{

  std::string filen="../data/MarkSystNumbers.txt";
  
  
  Double_t X0[nrT], X0stat[nrT], JESup[nrT], JESdw[nrT], JER[nrT], bTAGup[nrT], bTAGdw[nrT];
  int m_debug;
    //Int_t nr = 0;
    if (m_debug>1) std::cout<<" file to be read = <"<<filen<<">"<<std::endl;
    
    int nr=0;
    std::string STRING;
    stringstream ss;
    ifstream infile(filen.c_str());
    //Double_t* X1,X2,X3,X4;
    if (infile.is_open())                     //if the file is open
    {
	while(!infile.eof()) // To get you all the lines.
	{
	    ss.clear();
	    getline(infile,STRING); // Saves the line in STRING.
	    if (STRING.substr(0,1)=="*" || STRING.size()==0) {
	      if (STRING.substr(0,1)=="*") std::cout<<" ----------------OFFSET "<<nr<<std::endl;  
	      continue; // skip comments
	    }
	    if (m_debug==0) std::cout<<" now read --- "<<STRING<<std::endl; // Prints our STRING.
	    ss << STRING;
	    //	    std::cout<<"size of line read ="<<STRING.size()<<std::endl;
	    // 
	    Double_t ax1,ax2,ax3,ax4,ax5,ax6,ax7;
	    ss>>ax1>>ax2>>ax3>>ax4>>ax5>>ax6>>ax7;
	    //if (m_debug==0) std::cout<<"==== at index "<<nr<<" aX1-7 are "<<ax1<<" "<<ax2<<" "<<ax3<<" "<<ax4<<" "<<ax5<<" "<<ax6<<" "<<ax7<<std::endl;	
	    //Double_t X0[nr], X0stat[nr], Double_t JERup[nr], Double_t JERdw[nr], Double_t JES[nr], Double_t bTAGup[nr], Double_t bTAGdw[nr];
	    X0[nr]=ax1;
	    X0stat[nr]=       ax2/ax1+1;
	    JESup[nr]=   abs(ax3)/ax1+1.;
	    JESdw[nr]=   1.-abs(ax4)/ax1;
	    JER[nr]=     abs(ax5)/ax1+1.;
	    bTAGup[nr]=  abs(ax6)/ax1+1.;
	    bTAGdw[nr]=  1.-abs(ax7)/ax1;
	    ++nr;
	    if (nr==100) break;
	}
    }
    else 
    {
	std::cout<<"Cannot open file "<<std::endl;
	return;	
    }
    
    infile.close();

    std::cout<<" number of lines "<<std::endl;

    TFile* newFile ;
    if (saving) newFile  = new TFile("MarkSyst.root","recreate");

    

    TH1F *h_Data0stat       ;
    TH1F *h_Data0JERup      ;
    TH1F *h_Data0JERdw      ;
    TH1F *h_Data0JES        ;
    TH1F *h_Data0BTAGup     ;
    TH1F *h_Data0BTAGdw     ;

    TH1F *h_Data1stat       ;
    TH1F *h_Data1JERup      ;
    TH1F *h_Data1JERdw      ;
    TH1F *h_Data1JES        ;
    TH1F *h_Data1BTAGup     ;
    TH1F *h_Data1BTAGdw     ;

    TH1F *h_Data2stat       ;
    TH1F *h_Data2JERup      ;
    TH1F *h_Data2JERdw      ;
    TH1F *h_Data2JES        ;
    TH1F *h_Data2BTAGup     ;
    TH1F *h_Data2BTAGdw     ;

    TH1F *h_Data3stat       ;
    TH1F *h_Data3JERup      ;
    TH1F *h_Data3JERdw      ;
    TH1F *h_Data3JES        ;
    TH1F *h_Data3BTAGup     ;
    TH1F *h_Data3BTAGdw     ;

    TH1F *h_Data4stat       ;
    TH1F *h_Data4JERup      ;
    TH1F *h_Data4JERdw      ;
    TH1F *h_Data4JES        ;
    TH1F *h_Data4BTAGup     ;
    TH1F *h_Data4BTAGdw     ;


    TH1F *h_Data5stat       ;
    TH1F *h_Data5JERup      ;
    TH1F *h_Data5JERdw      ;
    TH1F *h_Data5JES        ;
    TH1F *h_Data5BTAGup     ;
    TH1F *h_Data5BTAGdw     ;


    for (int jj=0; jj<6; ++jj){
      int analysisBin = jj;
      std::cout<<" analysisBin = "<<jj<<std::endl;
      if(analysisBin==0)
    {
      const int myn = 7 ;
      Double_t mynewBinsPlot0   [myn] = {20,30,50,75,110,200,500};

      // dummy[0]= 0 ;
      // dummy[1]= 0 ;
      // dummy[2]= 0 ;
      // dummy[3]= 0 ;
      // dummy[4]= 0 ;
      // dummy[5]= 0 ;
      h_Data0stat        = new TH1F("h_Data0stat",""        ,myn-1,mynewBinsPlot0);
      h_Data0JESup       = new TH1F("h_Data0JERup",""       ,myn-1,mynewBinsPlot0);
      h_Data0JESdw       = new TH1F("h_Data0JERdw",""       ,myn-1,mynewBinsPlot0);
      h_Data0JER         = new TH1F("h_Data0JES",""         ,myn-1,mynewBinsPlot0);
      h_Data0BTAGup      = new TH1F("h_Data0BTAGup",""      ,myn-1,mynewBinsPlot0);
      h_Data0BTAGdw      = new TH1F("h_Data0BTAGdw",""      ,myn-1,mynewBinsPlot0);
      int hoff = 0;
      for (int j=1; j<myn; ++j) h_Data0stat ->SetBinContent(j,X0stat[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data0JESup->SetBinContent(j,JESup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data0JESdw->SetBinContent(j,JESdw[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data0JER  ->SetBinContent(j,JER[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data0BTAGup->SetBinContent(j,bTAGup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data0BTAGdw->SetBinContent(j,bTAGdw[hoff+j-1]);
    }
  if(analysisBin==1)
    {	  
      const int myn = 9 ;
      Double_t mynewBinsPlot1   [myn] = {0,0.2,0.4,0.6,0.8,1.2,1.6,2.0,2.4};
	
      // dummy[0]= 0 ;
      // dummy[1]= 0 ;
      // dummy[2]= 0 ;
      // dummy[3]= 0 ;
      // dummy[4]= 0 ;
      // dummy[5]= 0 ;
      // dummy[6]= 0 ;
      // dummy[7]= 0 ;
      h_Data1stat        = new TH1F("h_Data1stat",""        ,myn-1,mynewBinsPlot1);
      h_Data1JESup       = new TH1F("h_Data1JERup",""       ,myn-1,mynewBinsPlot1);
      h_Data1JESdw       = new TH1F("h_Data1JERdw",""       ,myn-1,mynewBinsPlot1);
      h_Data1JER         = new TH1F("h_Data1JES",""         ,myn-1,mynewBinsPlot1);
      h_Data1BTAGup      = new TH1F("h_Data1BTAGup",""      ,myn-1,mynewBinsPlot1);
      h_Data1BTAGdw      = new TH1F("h_Data1BTAGdw",""      ,myn-1,mynewBinsPlot1);
      int hoff = 6;
      for (int j=1; j<myn; ++j) h_Data1stat ->SetBinContent(j,X0stat[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data1JESup->SetBinContent(j,JESup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data1JESdw->SetBinContent(j,JESdw[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data1JER  ->SetBinContent(j,JER[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data1BTAGup->SetBinContent(j,bTAGup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data1BTAGdw->SetBinContent(j,bTAGdw[hoff+j-1]);

    }
  if(analysisBin==2)		  	 
    {
      const int myn = 9 ;
      Double_t mynewBinsPlot2   [myn] = {1,20,30,40,60,80,110,200,500};
	
      // dummy[0]= 0 ;
      // dummy[1]= 0 ;
      // dummy[2]= 0 ;
      // dummy[3]= 0 ;
      // dummy[4]= 0 ;
      // dummy[5]= 0 ;
      // dummy[6]= 0 ;
      // dummy[7]= 0 ;
      h_Data2stat        = new TH1F("h_Data2stat",""        ,myn-1,mynewBinsPlot2);
      h_Data2JESup       = new TH1F("h_Data2JERup",""       ,myn-1,mynewBinsPlot2);
      h_Data2JESdw       = new TH1F("h_Data2JERdw",""       ,myn-1,mynewBinsPlot2);
      h_Data2JER         = new TH1F("h_Data2JES",""         ,myn-1,mynewBinsPlot2);
      h_Data2BTAGup      = new TH1F("h_Data2BTAGup",""      ,myn-1,mynewBinsPlot2);
      h_Data2BTAGdw      = new TH1F("h_Data2BTAGdw",""      ,myn-1,mynewBinsPlot2);
      int hoff = 31;
      for (int j=1; j<myn; ++j) h_Data2stat ->SetBinContent(j,X0stat[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data2JESup->SetBinContent(j,JESup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data2JESdw->SetBinContent(j,JESdw[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data2JER  ->SetBinContent(j,JER[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data2BTAGup->SetBinContent(j,bTAGup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data2BTAGdw->SetBinContent(j,bTAGdw[hoff+j-1]);
    }
  if(analysisBin==3)		  	 
    {
      const int myn = 9 ;
      Double_t mynewBinsPlot3   [myn] = {0.0,0.2,0.4,0.6,0.8,1.2,1.6,2.0,2.5};
	
      // dummy[0]= 0 ;
      // dummy[1]= 0 ;
      // dummy[2]= 0 ;
      // dummy[3]= 0 ;
      // dummy[4]= 0 ;
      // dummy[5]= 0 ;
      // dummy[6]= 0 ;
      // dummy[7]= 0 ;
      h_Data3stat        = new TH1F("h_Data3stat",""        ,myn-1,mynewBinsPlot3);
      h_Data3JESup       = new TH1F("h_Data3JERup",""       ,myn-1,mynewBinsPlot3);
      h_Data3JESdw       = new TH1F("h_Data3JERdw",""       ,myn-1,mynewBinsPlot3);
      h_Data3JER         = new TH1F("h_Data3JES",""         ,myn-1,mynewBinsPlot3);
      h_Data3BTAGup      = new TH1F("h_Data3BTAGup",""      ,myn-1,mynewBinsPlot3);
      h_Data3BTAGdw      = new TH1F("h_Data3BTAGdw",""      ,myn-1,mynewBinsPlot3);
      int hoff = 39;
      for (int j=1; j<myn; ++j) h_Data3stat ->SetBinContent(j,X0stat[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data3JESup->SetBinContent(j,JESup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data3JESdw->SetBinContent(j,JESdw[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data3JER  ->SetBinContent(j,JER[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data3BTAGup->SetBinContent(j,bTAGup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data3BTAGdw->SetBinContent(j,bTAGdw[hoff+j-1]);
    }
  if(analysisBin==4)
    {
      const int myn = 9 ;
      Double_t mynewBinsPlot4   [myn] = {0,0.5,1,1.5,2.0,2.4,2.8,3,TMath::Pi()};
	
      // dummy[0]= 0 ;
      // dummy[1]= 0 ;
      // dummy[2]= 0 ;
      // dummy[3]= 0 ;
      // dummy[4]= 0 ;
      // dummy[5]= 0 ;
      // dummy[6]= 0 ;
      // dummy[7]= 0 ;
      // dummy[8]= 0 ;
      h_Data4stat        = new TH1F("h_Data4stat",""        ,myn-1,mynewBinsPlot4);
      h_Data4JESup       = new TH1F("h_Data4JERup",""       ,myn-1,mynewBinsPlot4);
      h_Data4JESdw       = new TH1F("h_Data4JERdw",""       ,myn-1,mynewBinsPlot4);
      h_Data4JER         = new TH1F("h_Data4JES",""         ,myn-1,mynewBinsPlot4);
      h_Data4BTAGup      = new TH1F("h_Data4BTAGup",""      ,myn-1,mynewBinsPlot4);
      h_Data4BTAGdw      = new TH1F("h_Data4BTAGdw",""      ,myn-1,mynewBinsPlot4);
      int hoff = 14;
      for (int j=1; j<myn; ++j) h_Data4stat ->SetBinContent(j,X0stat[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data4JESup->SetBinContent(j,JESup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data4JESdw->SetBinContent(j,JESdw[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data4JER  ->SetBinContent(j,JER[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data4BTAGup->SetBinContent(j,bTAGup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data4BTAGdw->SetBinContent(j,bTAGdw[hoff+j-1]);
    }
  if(analysisBin==5)
    {				  	 
      const int myn = 10 ;
      Double_t mynewBinsPlot5   [myn] = {0,1,1.5,2,2.5,3,3.5,4,4.5,6};   
	
      // dummy[0]= 0 ;
      // dummy[1]= 0 ;
      // dummy[2]= 0 ;
      // dummy[3]= 0 ;
      // dummy[4]= 0 ;
      // dummy[5]= 0 ;
      // dummy[6]= 0 ;
      // dummy[7]= 0 ;
      h_Data5stat        = new TH1F("h_Data5stat",""        ,myn-1,mynewBinsPlot5);
      h_Data5JESup       = new TH1F("h_Data5JERup",""       ,myn-1,mynewBinsPlot5);
      h_Data5JESdw       = new TH1F("h_Data5JERdw",""       ,myn-1,mynewBinsPlot5);
      h_Data5JER         = new TH1F("h_Data5JES",""         ,myn-1,mynewBinsPlot5);
      h_Data5BTAGup      = new TH1F("h_Data5BTAGup",""      ,myn-1,mynewBinsPlot5);
      h_Data5BTAGdw      = new TH1F("h_Data5BTAGdw",""      ,myn-1,mynewBinsPlot5);
      int hoff = 22;
      for (int j=1; j<myn; ++j) h_Data5stat ->SetBinContent(j,X0stat[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data5JESup->SetBinContent(j,JESup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data5JESdw->SetBinContent(j,JESdw[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data5JER  ->SetBinContent(j,JER[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data5BTAGup->SetBinContent(j,bTAGup[hoff+j-1]);
      for (int j=1; j<myn; ++j) h_Data5BTAGdw->SetBinContent(j,bTAGdw[hoff+j-1]);
    }
    }

  TCanvas* c0 = new TCanvas("slide 0","0");
  h_Data0stat->GetYaxis()->SetRangeUser(0.5,1.5);
  h_Data0stat->Draw();
  h_Data0JESup->Draw("SAME");
  h_Data0JESdw->Draw("SAME");
  h_Data0JER->Draw("SAME");
  h_Data0BTAGup->Draw("SAME");
  h_Data0BTAGdw->Draw("SAME");

  TCanvas* c1 = new TCanvas("slide 1","1");
  h_Data1stat->GetYaxis()->SetRangeUser(0.5,1.5);
  h_Data1stat->Draw();
  h_Data1JESup->Draw("SAME");
  h_Data1JESdw->Draw("SAME");
  h_Data1JER->Draw("SAME");
  h_Data1BTAGup->Draw("SAME");
  h_Data1BTAGdw->Draw("SAME");

  TCanvas* c2 = new TCanvas("slide 2","2");
  h_Data2stat->GetYaxis()->SetRangeUser(0.5,1.5);
  h_Data2stat->Draw();
  h_Data2JESup->Draw("SAME");
  h_Data2JESdw->Draw("SAME");
  h_Data2JER->Draw("SAME");
  h_Data2BTAGup->Draw("SAME");
  h_Data2BTAGdw->Draw("SAME");

  TCanvas* c3 = new TCanvas("slide 3","3");
  h_Data3stat->GetYaxis()->SetRangeUser(0.5,1.5);
  h_Data3stat->Draw();
  h_Data3JESup->Draw("SAME");
  h_Data3JESdw->Draw("SAME");
  h_Data3JER->Draw("SAME");
  h_Data3BTAGup->Draw("SAME");
  h_Data3BTAGdw->Draw("SAME");

  TCanvas* c4 = new TCanvas("slide 4","4");
  h_Data4stat->GetYaxis()->SetRangeUser(0.5,1.5);
  h_Data4stat->Draw();
  h_Data4JESup->Draw("SAME");
  h_Data4JESdw->Draw("SAME");
  h_Data4JER->Draw("SAME");
  h_Data4BTAGup->Draw("SAME");
  h_Data4BTAGdw->Draw("SAME");

  TCanvas* c5 = new TCanvas("slide 5","5");
  h_Data5stat->GetYaxis()->SetRangeUser(0.5,1.5); 
  h_Data5stat->Draw();
  //h_Data5JESup->Draw("SAME");
  //h_Data5JESdw->Draw("SAME");
  h_Data5JER->Draw("SAME");
  //h_Data5BTAGup->Draw("SAME");
  //h_Data5BTAGdw->Draw("SAME");



  if (saving){
    newFile->Write();
    newFile->Close();
  }
    return;
}