{
#include <vector>
#include<TH1.h>
#include<TROOT.h>
#include<TChain.h>
#include<TFile.h>
#include<TStyle.h>
  //this is for jet graphs
  // Style setting
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetFillColor(0);
  gStyle->SetFuncColor(kRed);
  gStyle->SetFuncWidth(1.0);
  gStyle->SetOptStat(111111);
//  gStyle->SetOptFit(111111);

//open files
  // read a string via file since long string causes memory error in CINT when it is read via stdin
  std::string argStr;
  char buf[256+1];
  unsigned int delpos;
  std::ifstream ifs("input.txt");
  while (true)
    {
      ifs.read(buf,256);
      if (ifs.eof())
        {
          if (ifs.gcount() == 0) break;
          delpos = ifs.gcount()-1;
        }
      else
        {
          delpos = ifs.gcount();
        }
      buf[delpos] = 0x00;
      argStr += buf;
    }

  // split by ','
  std::vector<std::string> fileList;
  for (size_t i=0,n; i <= argStr.length(); i=n+1)
    {
      n = argStr.find_first_of(',',i);
      if (n == string::npos)
        n = argStr.length();
      string tmp = argStr.substr(i,n-i);
      fileList.push_back(tmp);
    }

  // open input files
  TChain fChain("qcd");
  for (int iFile=0; iFile<fileList.size(); ++iFile)
    {
      std::cout << "open " << fileList[iFile].c_str() << std::endl;
      fChain.Add(fileList[iFile].c_str());
    }
  const Double_t pi = TMath::Pi();
  // Make a canvas
  TCanvas c1;
  TCanvas c2; 
  TCanvas c3; 
  //TCanvas c4; 

  // ROOT-blah-blah needed for release 14
  //gSystem->Load("libCintex");
  //ROOT::Cintex::Cintex::Enable(); 

  //   Open file
  //TFile f("/scruser/3/yjf/group10.perf-jets.data10_7TeV.00155228.physics_L1Calo.recon.ESD.f259_JetEtMissDPDModifier000025.v1.EXT0._00001.qcd.root");
  //TChain * MYTREE = new TChain("JetMETD3PD");
  //MYTREE->Add("user.jyu.JetD3PD.AANT._00002.root");
  //MYTREE->Add("user.jyu.JetD3PD.AANT._00003.root");
  
  // Get the number of entries in branch
  Int_t nEntriesNtuple = (Int_t)qcd->GetEntries();
  std::cout << "Number of entries in the n-tuple is " << nEntriesNtuple << std::endl;
  
  // Declare variables to plot, remember to check the class of the variables
  std::vector<float>* m_phi;
  std::vector<float>* m_energy;
  std::vector<float>* m_pt;
  std::vector<float>* m_eta;
  Int_t m_jetNum;    

  //declare the variables for clean cuts
  std::vector<float>* m_n90;
  std::vector<float>* m_timing;
  std::vector<float>* m_quality;
  std::vector<float>* m_HECf;
  std::vector<float>* m_emf;
  std::vector<int>* m_isGood;
  std::vector<float>* m_BCH_CORR_CELL;
  std::vector<float>* m_TileGap3Frac;

  // Assign variable to the n-tuple
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_phi",&m_phi);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_E",&m_energy);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_eta",&m_eta);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_pt",&m_pt); 
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_n",&m_jetNum);
  //assign the quality variables to the n-tuple
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_n90",&m_n90);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_timing",&m_timing);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_quality",&m_quality);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_hecf",&m_HECf);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_emfrac",&m_emf);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_isGood",&m_isGood);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_tgap3f",&m_TileGap3Frac);
  qcd->SetBranchAddress("jet_AntiKt4H1TopoNew_BCH_CORR_CELL",&m_BCH_CORR_CELL);

  // Declare histograms
  TH1* deltaPhiGoodHist = new TH1F("#Delta#phi", "#Delta#phi between the Leading Jet and the Second Leading Jet with Good Cut", 100, 0, pi); 
  deltaPhiGoodHist->GetXaxis()->SetTitle("#Delta#phi"); 
  deltaPhiGoodHist->GetYaxis()->SetTitle("Number of Events");
  TH1* deltaPhiUglyHist = new TH1F("#Delta#phi", "#Delta#phi between the Leading Jet and the Second Leading Jet with Ugly Cut", 100, 0, pi); 
  deltaPhiUglyHist->GetXaxis()->SetTitle("#Delta#phi"); 
  deltaPhiUglyHist->GetYaxis()->SetTitle("Number of Events");      
  TH1* deltaPhiBadHist = new TH1F("#Delta#phi", "#Delta#phi between the Leading Jet and the Second Leading Jet with Bad Cut", 100, 0, pi); 
  deltaPhiBadHist->GetXaxis()->SetTitle("#Delta#phi"); 
  deltaPhiBadHist->GetYaxis()->SetTitle("Number of Events");   
// Fill histograms with numbers
  for(Int_t k=0; k<nEntriesNtuple; ++k) {
    qcd->GetEntry(k); // Fills all assigned variables in one go, for ith event

    for (Int_t j = 0; j<m_jetNum; ++j) { 
	//scaling of some variables
	double energy = (*m_energy)[j]/1000;
	double pt = (*m_pt)[j]/1000;
	//the quality variables declaration
  
	int n90 = (*m_n90)[j];
	double HECf = (*m_HECf)[j];
	double time = (*m_timing)[j];
	double quality = (*m_quality)[j];
	double emf = (*m_emf)[j];
	double TileGap3Frac = (*m_TileGap3Frac)[j];
	double BCH_CORR_CELL = (*m_BCH_CORR_CELL)[j];
	int isGood = (*m_isGood)[j];

	//fill the histograms
//	For Good Cuts
	if(isGood==2 && m_jetNum>=2&&pt>15){
	double deltaPhiGood = abs((*m_phi)[0]-(*m_phi)[1]);
	if(deltaPhiGood>pi){
	deltaPhiGood = deltaPhiGood-pi; 
	}
	deltaPhiGoodHist->Fill(deltaPhiGood); 
	}
//	For Ugly Cuts
//	Note if both bad and ugly, here count as ugly. 
	if(isGood==1 && m_jetNum>=2&&pt>15){
	double deltaPhiUgly = abs((*m_phi)[0]-(*m_phi)[1]);
	if(deltaPhiUgly>pi){
	deltaPhiUgly = deltaPhiUgly-pi; 
	}
	deltaPhiUglyHist->Fill(deltaPhiUgly); 
	}

//	For Bad Cuts
	if(n90<=5&&HECf>0.8)||(fabs(time)>50)||(abs(quality)<0.8&&emf<0.95)&&(TileGap3Frac<=0.5&&BCH_CORR_CELL<=0.5&&m_jetNum>=2&&pt>15)){
 // 	if(isGood==0 && m_jetNum>=2 && pt>=15){
	double deltaPhiBad = abs((*m_phi)[0]-(*m_phi)[1]);
	if(deltaPhiBad>pi){
	deltaPhiBad = deltaPhiBad-pi; 
	} 
	deltaPhiBadHist->Fill(deltaPhiBad); 
	}

    }
	
     
  }

  // Draw histograms
 c1.cd();
deltaPhiGoodHist->Draw();
c2.cd();
deltaPhiUglyHist->Draw();
c3.cd();
deltaPhiBadHist->Draw();

//write to file
TFile saveFile ("saving.root", "RECREATE"); 
deltaPhiGoodHist->Write(); 
deltaPhiUglyHist->Write();
deltaPhiBadHist->Write(); 
saveFile.Close(); 
}