#include <iomanip>
#define DBGFlag 1

void ReadFitResults(int analysisBin = 0, int SystVar =0)
{

  char FileName[1000]       ;  
  fstream FileToOpen        ;  

  char prefix[1000];
  char DisPath[1000];
  char Channel[1000];
  char DisType[1000];
  char SystType[1000];
  char NameDirectory[1000];
  char OutputFileName[1000];

  sprintf(prefix,"/afs/le.infn.it/user/n/niko/Zb2011macros/");
  
  if(analysisBin==0)sprintf(DisPath,"DirFinalRunZbFit_jetPt/");
  if(analysisBin==1)sprintf(DisPath,"DirFinalRunZbFit_jetY/");
  if(analysisBin==2)sprintf(DisPath,"DirFinalRunZbFit_ZPt/");
  if(analysisBin==3)sprintf(DisPath,"DirFinalRunZbFit_ZY/");
  if(analysisBin==4)sprintf(DisPath,"DirFinalRunZbFit_DPhi/");
  if(analysisBin==5)sprintf(DisPath,"DirFinalRunZbFit_DR/");

  sprintf(Channel,"FinalFitResults/CombinedChannel/");
  
  if(analysisBin==0)sprintf(DisType,"JetPt");
  if(analysisBin==1)sprintf(DisType,"JetY");
  if(analysisBin==2)sprintf(DisType,"ZPt");
  if(analysisBin==3)sprintf(DisType,"ZY");
  if(analysisBin==4)sprintf(DisType,"DPhi");
  if(analysisBin==5)sprintf(DisType,"DR");
  
  if(SystVar==0)sprintf(SystType,"Central");
  if(SystVar==1)sprintf(SystType,"BkgDown");
  if(SystVar==2)sprintf(SystType,"BkgUp");
  if(SystVar==3)sprintf(SystType,"CharmRew");
  if(SystVar==4)sprintf(SystType,"EvtGenRew");
  if(SystVar==5)sprintf(SystType,"HFOR");
  if(SystVar==6)sprintf(SystType,"JER");
  if(SystVar==7)sprintf(SystType,"JESdown");
  if(SystVar==8)sprintf(SystType,"JESup");
  if(SystVar==9)sprintf(SystType,"MPIdown");
  if(SystVar==10)sprintf(SystType,"MPIup");
  if(SystVar==11)sprintf(SystType,"StatInTemplates");
  if(SystVar==12)sprintf(SystType,"MET1");
  if(SystVar==13)sprintf(SystType,"MET2");
  if(SystVar==14)sprintf(SystType,"MET3");
  if(SystVar==15)sprintf(SystType,"MET4");
  if(SystVar==16)sprintf(SystType,"MET5");
  if(SystVar==17)sprintf(SystType,"MET6");
  if(SystVar==18)sprintf(SystType,"MET7");
  if(SystVar==19)sprintf(SystType,"MET8");
  if(SystVar==20)sprintf(SystType,"MET9");
  if(SystVar==21)sprintf(SystType,"MET10");
  if(SystVar==22)sprintf(SystType,"MET11");
  if(SystVar==23)sprintf(SystType,"MET12");
  if(SystVar==24)sprintf(SystType,"MET13");
  if(SystVar==25)sprintf(SystType,"MET14");
  if(SystVar==26)sprintf(SystType,"MET15");
  if(SystVar==27)sprintf(SystType,"BTAG10_down");
  if(SystVar==28)sprintf(SystType,"BTAG10_up");
  if(SystVar==29)sprintf(SystType,"BTAG1_down");
  if(SystVar==30)sprintf(SystType,"BTAG1_up");
  if(SystVar==31)sprintf(SystType,"BTAG2_down");
  if(SystVar==32)sprintf(SystType,"BTAG2_up");
  if(SystVar==33)sprintf(SystType,"BTAG3_down");
  if(SystVar==34)sprintf(SystType,"BTAG3_up");
  if(SystVar==35)sprintf(SystType,"BTAG4_down");
  if(SystVar==36)sprintf(SystType,"BTAG4_up");
  if(SystVar==37)sprintf(SystType,"BTAG5_down");
  if(SystVar==38)sprintf(SystType,"BTAG5_up");
  if(SystVar==39)sprintf(SystType,"BTAG6_down");
  if(SystVar==40)sprintf(SystType,"BTAG6_up");
  if(SystVar==41)sprintf(SystType,"BTAG7_down");
  if(SystVar==42)sprintf(SystType,"BTAG7_up");
  if(SystVar==43)sprintf(SystType,"BTAG8_down");
  if(SystVar==44)sprintf(SystType,"BTAG8_up");
  if(SystVar==45)sprintf(SystType,"BTAG9_down");
  if(SystVar==46)sprintf(SystType,"BTAG9_up");
  if(SystVar==47)sprintf(SystType,"CTAG_down");
  if(SystVar==48)sprintf(SystType,"CTAG_up");
  if(SystVar==49)sprintf(SystType,"LTAG_down");
  if(SystVar==50)sprintf(SystType,"LTAG_up");
  
  if(SystVar!=11)sprintf(NameDirectory,"%s%s%s%s/%s/FitResultsLogForPlot.txt",prefix,DisPath,Channel,DisType,SystType);
  else sprintf(NameDirectory,"%s%s%s%s/%s/ToyMCForStatLog.txt",prefix,DisPath,Channel,DisType,SystType);

  sprintf(FileName,"%s",NameDirectory);   FileToOpen.open(FileName,fstream::in); 

  sprintf(OutputFileName,"FitSummary_CombinedChannel_%s_%s.root",DisType,SystType);
  std::cout<<"\n Opening -> "<<NameDirectory<<"\n"<<std::endl;

  const int n = 100;   float xVar[n] ; float xVarErr[n] ;  float dummy[n] ; //const int myn ;
  
  if(analysisBin==0)
    {
      const int myn = 7 ;
      Double_t mynewBins       [myn] = {20,30,50,75,110,200,500};
      Double_t mynewBinsPlot   [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 ;
    }
  if(analysisBin==1)
    {	  
      const int myn = 9 ;
      Double_t mynewBins       [myn] = {0,0.2,0.4,0.6,0.8,1.2,1.6,2.0,2.4};
      Double_t mynewBinsPlot   [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 ;
    }
  if(analysisBin==2)		  	 
    {
      const int myn = 9 ;
      Double_t mynewBins       [myn] = {0,20,30,40,60,80,110,200,500};
      Double_t mynewBinsPlot   [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 ;
    }
  if(analysisBin==3)		  	 
    {
      const int myn = 9 ;
      Double_t mynewBins       [myn] = {0.0,0.2,0.4,0.6,0.8,1.2,1.6,2.0,2.5};
      Double_t mynewBinsPlot   [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 ;
    }
  if(analysisBin==4)
    {
      const int myn = 9 ;
      Double_t mynewBins       [myn] = {0,0.5,1,1.5,2.0,2.4,2.8,3,TMath::Pi()};
      Double_t mynewBinsPlot   [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 ;
    }
  if(analysisBin==5)
    {				  	 
      const int myn = 10 ;
      Double_t mynewBins       [myn] = {0,1,1.5,2,2.5,3,3.5,4,4.5,6};   
      Double_t mynewBinsPlot   [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 ;
    }

  //*****************************************
  // All the standard systematics
  //*****************************************
 
  std::string sFileRow;
  //To read the standard fit output
  std::string cuts , tags , ptmins , ptmaxs , ymins , ymaxs , Nls , Nlerrs , Ncs , Ncerrs , Nbs , Nberrs , lNorms , cNorms , bNorms ;
  std::string mycut, tag  ;
  float ptmin , ptmax , ymin , ymax , Nl , Nlerr , Nc , Ncerr , Nb , Nberr , lNorm , cNorm , bNorm , chi ;
  //To read the output for stat. in templates
  float errMC , errMCErr ;
  std::string s1 , s2 , s3 , s4 , s5 , s6 , s7;   
  
  float vchi[n] ;  
  float vNb[n];     float vNc[n];     float vNl[n];   
  float vNbErr[n];  float vNcErr[n];  float vNlErr[n];
  float vNbMC[n];   float vNcMC[n];   float vNlMC[n];	  	 
  
  float vStatInTemp[n];  float vStatInTempErr[n];
  
  int index1  = 0;
  
  if(DBGFlag==1)std::cout<<"\n"<<std::endl;
    
  if(SystVar!=11)
    {
      while( getline(FileToOpen,sFileRow) ) 
	{
	  istringstream ss(sFileRow);
	  ss >>cuts >> mycut  >> tags   >> tag >> ptmins >> ptmin >> ptmaxs >> ptmax >> ymins >> ymin >> ymaxs >> ymax >>
	    Nls    >>  Nl   >> Nlerrs >> Nlerr >> 
	    Ncs    >>  Nc   >> Ncerrs >> Ncerr >>
	    Nbs    >>  Nb   >> Nberrs >> Nberr >>
	    lNorms >> lNorm >>
	    cNorms >> cNorm >>
	    bNorms >> bNorm >>
	    chi ;
	  
	  if(DBGFlag==1){
	    std::cout <<
	      cuts  << " " << mycut << " "<< tags   << " "<< tag   << " "<< ptmins << " "<< ptmin << " "<< ptmaxs << " "<< ptmax << " "<< ymins << " "<< ymin << " "<< ymaxs << " "<< ymax << " "<< 
	      Nls   << " " <<  Nl   << " "<< Nlerrs << " "<< Nlerr << " "<<
	      Ncs   << " " <<  Nc   << " "<< Ncerrs << " "<< Ncerr << " "<<
	      Nbs   << " " <<  Nb   << " "<< Nberrs << " "<< Nberr << " "<<
	      lNorms<< " " << lNorm << " "<<
	      cNorms<< " " << cNorm << " "<<
	      bNorms<< " " << bNorm << " "<< 
	      chi
		    <<std::endl;
	  }

	  xVar     [index1] = (float)(mynewBins[index1+1]+mynewBins[index1])/2;
	  xVarErr  [index1] = (float)(mynewBins[index1+1]-mynewBins[index1])/2;
	  vNb      [index1] = Nb/(2*xVarErr[index1]) ;
	  vNc      [index1] = Nc/(2*xVarErr[index1]) ;
	  vNl      [index1] = Nl/(2*xVarErr[index1]) ;
	  //**************************//
	  vchi     [index1] = chi                          ;
	  vNbErr   [index1] = Nberr/(2*xVarErr[index1])    ;
	  vNcErr   [index1] = Ncerr/(2*xVarErr[index1])    ;
	  vNlErr   [index1] = Nlerr/(2*xVarErr[index1])    ;
	  //*************************//
	  vNbMC    [index1] = bNorm/(2*xVarErr[index1]);   
	  vNcMC    [index1] = cNorm/(2*xVarErr[index1]);   
	  vNlMC    [index1] = lNorm/(2*xVarErr[index1]);

	  index1 ++;
	}
  
      if(DBGFlag==1)std::cout<<"\n"<<std::endl;

      TH1F *h_Data      = new TH1F("h_Data",""       ,myn-1,mynewBinsPlot);
      TH1F *h_DataErr   = new TH1F("h_DataErr",""    ,myn-1,mynewBinsPlot);
      TH1F *h_unData    = new TH1F("h_unData",""     ,myn-1,mynewBinsPlot);
      TH1F *h_MC        = new TH1F("h_MC",""         ,myn-1,mynewBinsPlot);
  
      for(int m = 1 ; m < myn ; m++)
	{
	  //std::cout<< "m = / vNbMC[m-1] = "<<m<< " " << vNbMC[m-1] << std::endl; 
	  h_MC      ->SetBinContent(m,vNbMC[m-1]);
	  h_Data    ->SetBinContent(m,vNb[m-1]);
	  h_Data    ->SetBinError(m,vNbErr[m-1]);
	  h_DataErr ->SetBinContent(m,vNbErr[m-1]);
	  h_unData  ->SetBinContent(m,(vNb[m-1]*(2*xVarErr[m-1])));
	  h_unData  ->SetBinError(m,(vNbErr[m-1]*(2*xVarErr[m-1])));
	}
  
      TGraphAsymmErrors *gChi  = new TGraphAsymmErrors(index1 ,xVar ,vchi  ,dummy   ,dummy   , dummy  , dummy  );
      TGraphAsymmErrors *gData = new TGraphAsymmErrors(index1 ,xVar ,vNb   ,xVarErr ,xVarErr , vNbErr , vNbErr );
      TGraphAsymmErrors *gMC   = new TGraphAsymmErrors(index1 ,xVar ,vNbMC ,xVarErr ,xVarErr , dummy  , dummy  );

      if(DBGFlag>1)
	{
	  TCanvas *mycnv = new TCanvas();
	  h_Data->Draw();
	  h_MC  ->Draw("same");
	  gData ->Draw("p");
	  gMC   ->Draw("l");
	}

      TFile *f = new TFile(OutputFileName,"recreate");
      h_Data   ->Write();
      h_unData ->Write();
      h_DataErr->Write();
      h_MC     ->Write();  
      gData    ->Write("gdata");
      gMC      ->Write("gmc");
      gChi     ->Write("gchi");
      f        ->Close();
  
      delete gChi;
      delete gData; 
      delete gMC;
      delete h_Data; 
      delete h_unData; 
      delete h_DataErr; 
      delete h_MC;
    }
  else
    {
      
      while( getline(FileToOpen,sFileRow) ) 
	{
	  istringstream ss(sFileRow);
	
	  ss >> ptmin >> s1 >> s2 >> s3 >> ptmax >> ymin >> s4 >> s5 >> s6 >> ymax >> errMC >> s7 >> errMCErr ; 

	  if(DBGFlag==1) std::cout<<"errMC/errMCErr  "<<errMC<<"/"<<errMCErr <<std::endl;
	
	  vStatInTemp   [index1] = errMC    ;
	  vStatInTempErr[index1] = errMCErr ;
	
	  index1 ++;
	}
    
    
      TH1F *h_Stat = new TH1F("h_Stat",""    ,myn-1,mynewBinsPlot);
  
      for(int m = 1 ; m < myn ; m++)  h_Stat->SetBinContent(m,vStatInTemp[m-1] );
      
      TGraphAsymmErrors *gStat = new TGraphAsymmErrors(index1 ,xVar , vStatInTemp,xVarErr ,xVarErr , dummy  , dummy  );

      TFile *f = new TFile(OutputFileName,"recreate");
      h_Stat->Write();
      gStat ->Write("gdata");
      f     ->Close();
  
      delete h_Stat;
      delete gStat;
    }
  
}