//////////////// cut on the mass: select eg 80-100 (opposite charge muon pairs) to select Z bosons ////////////////////////////////////
//////////////////// Tnourji Abdellah /////////////////////////////////////////////////////////////////////////////////////////////////

#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TObject.h"
#include "TPaveLabel.h"
#include "TGaxis.h"
#include  "TLorentzVector.h"
#include "TText.h"
#include "THStack.h"
#include "TBranch.h"
#include "Riostream.h"
#include "TStyle.h"
#include "TString.h"
#include "TLegend.h"
#include "TRandom3.h"
#include "TMath.h"
#include "math.h"
#include "TColor.h"

#include "TCanvas.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"

#include <vector>


using namespace std;

void HistoStyle(TH1* histo, int mk, int zr){
         histo-> SetMarkerStyle(mk);
         histo-> SetMinimum(0);    
         histo-> SetLineColor(kOrange);
         histo-> SetLineWidth(2);
	 histo-> SetFillStyle(zr);
	 histo-> SetFillColor(kOrange);
}

void Wmet5() {

	 TH1F::SetDefaultSumw2();
	 gStyle->SetOptTitle(1111); gStyle->SetOptStat(1111);
         gStyle->SetOptFit(111); gStyle->SetStatBorderSize(0);
         gStyle->SetStatX(.89); gStyle->SetStatY(.89);
	

	const double mumass = 105.6583745e-3;

	TFile *hfile_0 = new TFile("/eos/cms/store/group/phys_heavyions/dileptons/echapon/metPbPb2018/tree_Data_Wmumu.root");
	TFile *hfile_1 = new TFile("/eos/cms/store/group/phys_heavyions/dileptons/echapon/metPbPb2018/tree_MC_Wjets_1.root");
	
	
	TTree* myT_0 = (TTree*) hfile_0->Get("particleTree");
	TTree* myT = (TTree*) hfile_1->Get("particleTree");

	// ParticaleFlowTree Branche
	TTree* myTPF = (TTree*) hfile_0->Get("pfTree");
	TTree* myTPF_1 = (TTree*) hfile_1->Get("pfTree");

	// HiTree Branche
	TTree* myTH = (TTree*) hfile_0->Get("hiTree");
	TTree* myTH1 = (TTree*) hfile_1->Get("hiTree");

	
	// MAKE Class 
	//myT->MakeClass("pfTree"); 

	// intialisation of  variables 
	// Declaration of leaf types
	

	Int_t            hiBin = 0;
	Int_t            nMu = 0;
        vector<float>   *muPt = 0;
        vector<float>   *muEta = 0 ;
        vector<float>   *muPhi = 0;
        vector<int>     *muCharge = 0;
        vector<float>   *muChi2NDF = 0;
        vector<float>   *muInnerD0 = 0;
        vector<float>   *muInnerDz = 0;
        vector<int>     *muTrkLayers = 0;
        vector<int>     *muPixelHits = 0;
        vector<int>     *muMuonHits = 0;
        vector<int>     *muStations = 0;

	
	// Declaration of leaf types
        Int_t           nPFpart = 0;
        vector<int>     *pfId = 0;
        vector<float>   *pfPt = 0;
        vector<float>   *pfEnergy = 0;
        vector<float>   *pfEta = 0;
        vector<float>   *pfPhi = 0;

	//Data

	myT->SetBranchStatus("nMu",1);
        myT->SetBranchStatus("muPt",1);
        myT->SetBranchStatus("muEta",1);
        myT->SetBranchStatus("muPhi",1);
        myT->SetBranchStatus("muCharge",1);
        myT->SetBranchStatus("muChi2NDF",0);
        myT->SetBranchStatus("muInnerD0",0);
        myT->SetBranchStatus("muInnerDz",0);
        myT->SetBranchStatus("muMuonHits",0);
        myT->SetBranchStatus("muStations",0);
	myT->SetBranchAddress("muTrkLayers",0);
	myT->SetBranchAddress("muPixelHits",0);
	
        myT->SetBranchAddress("nMu",&nMu);
        myT->SetBranchAddress("muPt",&muPt);
        myT->SetBranchAddress("muPhi",&muPhi);
        myT->SetBranchAddress("muEta",&muEta);
	myT->SetBranchAddress("muCharge",&muCharge);
	myT->SetBranchAddress("muChi2NDF",&muChi2NDF);
	myT->SetBranchAddress("muInnerD0",&muInnerD0);
	myT->SetBranchAddress("muInnerDz",&muInnerDz);
	myT->SetBranchAddress("muTrkLayers",&muTrkLayers);
	myT->SetBranchAddress("muPixelHits",&muPixelHits);
	myT->SetBranchAddress("muMuonHits",&muMuonHits);
	myT->SetBranchAddress("muStations",&muStations);
	
	
	
	//MC

	myT_0->SetBranchStatus("nMu",1);
        myT_0->SetBranchStatus("muPt",1);
        myT_0->SetBranchStatus("muEta",1);
        myT_0->SetBranchStatus("muPhi",1);
        myT_0->SetBranchStatus("muCharge",1);
        myT_0->SetBranchStatus("muChi2NDF",0);
        myT_0->SetBranchStatus("muInnerD0",0);
        myT_0->SetBranchStatus("muInnerDz",0);
        myT_0->SetBranchStatus("muMuonHits",0);
        myT_0->SetBranchStatus("muStations",0);
	myT_0->SetBranchAddress("muTrkLayers",0);
	myT_0->SetBranchAddress("muPixelHits",0);

	myT_0->SetBranchAddress("nMu",&nMu);
        myT_0->SetBranchAddress("muPt",&muPt);
        myT_0->SetBranchAddress("muPhi",&muPhi);
        myT_0->SetBranchAddress("muEta",&muEta);
	myT_0->SetBranchAddress("muCharge",&muCharge);
	myT_0->SetBranchAddress("muChi2NDF",&muChi2NDF);
	myT_0->SetBranchAddress("muInnerD0",&muInnerD0);
	myT_0->SetBranchAddress("muInnerDz",&muInnerDz);
	myT_0->SetBranchAddress("muTrkLayers",&muTrkLayers);
	myT_0->SetBranchAddress("muPixelHits",&muPixelHits);
	myT_0->SetBranchAddress("muMuonHits",&muMuonHits);
	myT_0->SetBranchAddress("muStations",&muStations);

	
	// SetBranchSatue and link variables for PaticleFlowTree
	// DATA
        
        myTPF->SetBranchStatus("*",0);
        myTPF->SetBranchStatus("nPFpart",1);
        myTPF->SetBranchStatus("pfId",1);
        myTPF->SetBranchStatus("pfEta",1);
        myTPF->SetBranchStatus("pfPt",1);
        myTPF->SetBranchStatus("pfPhi",1);
        myTPF->SetBranchStatus("pfEnergy",1);
	
	myTPF->SetBranchAddress("nPFpart",&nPFpart);
        myTPF->SetBranchAddress("pfId",&pfId);
        myTPF->SetBranchAddress("pfEta",&pfEta);
        myTPF->SetBranchAddress("pfPt",&pfPt);
        myTPF->SetBranchAddress("pfPhi",&pfPhi);
        myTPF->SetBranchAddress("pfEnergy",&pfEnergy);

	// MC

        myTPF_1->SetBranchStatus("*",0);
        myTPF_1->SetBranchStatus("nPFpart",1);
        myTPF_1->SetBranchStatus("pfId",1);
        myTPF_1->SetBranchStatus("pfEta",1);
        myTPF_1->SetBranchStatus("pfPt",1);
        myTPF_1->SetBranchStatus("pfPhi",1);
        myTPF_1->SetBranchStatus("pfEnergy",1);
	
	myTPF_1->SetBranchAddress("nPFpart",&nPFpart);
        myTPF_1->SetBranchAddress("pfId",&pfId);
        myTPF_1->SetBranchAddress("pfEta",&pfEta);
        myTPF_1->SetBranchAddress("pfPt",&pfPt);
        myTPF_1->SetBranchAddress("pfPhi",&pfPhi);
        myTPF_1->SetBranchAddress("pfEnergy",&pfEnergy);

	myTH->SetBranchStatus("hiBin",1);
	myTH1->SetBranchStatus("hiBin",1);

	// hiBin MC
	myTH1->SetBranchAddress("hiBin",&hiBin);
	// hiBin Data	
	myTH->SetBranchAddress("hiBin",&hiBin);
	


	// Histogram Centrality 
        TH1F *hData = new TH1F("hData","hiBin Data ;Bin ;Entries",40,0,200);
	TH1F *hMC = new TH1F("hMC","hiBin MC;Bin;Entries",40,0,200);
	
	TH1F *Weight = new TH1F("Weight","",40,0,200);
	
	TH1F *hMC_W= new TH1F("hMC_W","",40,0,200);
	
	TH1F *h0 = new TH1F("h0","W mass Data;E^{miss}(GeV); Numbre of Events",100,0,1000);
	TH1F *h1 = new TH1F("h1","W mass MC;E^{miss}(GeV); Numbre of Events",100,0,1000);



	Long64_t nentriess = myT_0->GetEntries();
	Long64_t nentries_PF = myTPF->GetEntries();
	Long64_t nentries_H = myTH->GetEntries();
	
        for(Int_t i=0; i < nentriess; i++){  	
        myT_0->GetEntry(i);
	myTPF->GetEntry(i);
	myTH->GetEntry(i);

	
	if (muPt->size() < 1) continue;
	if (fabs(muEta->at(0)) > 2.4) continue;			
	if (muPt->at(0) < 50 ) continue;
        TLorentzVector mu0; mu0.SetPtEtaPhiM(muPt->at(0), muEta->at(0), muPhi->at(0), mumass);
	hData->Fill(hiBin);
	 }
	   	
	
	// MC hiBIN
	Long64_t nentriess1 = myT->GetEntries();
	Long64_t nentries_PF1 = myTPF_1->GetEntries();
	Long64_t nentries_H1 = myTH1->GetEntries();
	

        for(Int_t i=0; i < nentriess1; i++){  	
        myT->GetEntry(i);
	myTPF_1->GetEntry(i);
	myTH1->GetEntry(i);
	
	if (muPt->size() < 1) continue;
	if (fabs(muEta->at(0)) > 2.4) continue;			
	if (muPt->at(0) < 50 ) continue;
        TLorentzVector mu0; mu0.SetPtEtaPhiM(muPt->at(0), muEta->at(0), muPhi->at(0), mumass);
	hMC->Fill(hiBin);
	 }
	
	// Defined Weight 
	Weight->Divide(hData, hMC); 
	 

	//Data
	cout << "tnourji" <<endl;
	TLorentzVector mu0;
        for(Int_t i=0; i < nentriess; i++){  
     		
        myT_0->GetEntry(i);
	myTPF->GetEntry(i);
	myTH->GetEntry(i);
 		
	if (muPt->size()< 1) continue;	
	if (fabs(muEta->at(0)) > 2.4) continue;			
	if (muPt->at(0) < 50 ) continue;
	mu0.SetPtEtaPhiM(fabs(muPt->at(0)), muEta->at(0), muPhi->at(0), mumass);
	float Px (0), Py(0), Phi(0);
	float MET_Data(0);
	float  MT(0);
	TLorentzVector metvector;
	for( Int_t c = 0; c < nPFpart; c++){
	        if ( !(pfId->at(c)>=1 && pfId->at(c)<=5 ) ) continue; 
		Px +=(pfPt->at(c))*cos(pfPhi->at(c));
		Py +=(pfPt->at(c))*sin(pfPhi->at(c));
        }
	MET_Data = sqrt(pow(Px,2)+pow(Py,2));
	h0->Fill(MET_Data);
	
	}

	/////////////////////////////////////////////////////////
	
	// MC
	cout << "tnourji1" <<endl;
	TLorentzVector mu1;
        for(Int_t i=0; i < nentriess1; i++){  
        myT->GetEntry(i);
	myTPF_1->GetEntry(i);
	myTH1->GetEntry(i);	
	if (muPt->size()<1) continue;
	if (fabs(muEta->at(0)) > 2.4) continue;			
	if (muPt->at(0) < 50 ) continue;
	
	mu1.SetPtEtaPhiM(fabs(muPt->at(0)), muEta->at(0), muPhi->at(0), mumass);

	float Px (0), Py(0), Phi(0);
	float MET_MC(0);
	float  MT_MC(0);

	
	TLorentzVector metvector;
	
	for( Int_t c = 0; c < nPFpart; c++){

	        if ( !(pfId->at(c)>=1 && pfId->at(c)<=5 ) ) continue; 
		Px +=(pfPt->at(c))*cos(pfPhi->at(c));
		Py +=(pfPt->at(c))*sin(pfPhi->at(c));
        }

	
	MET_MC = sqrt(pow(Px,2)+pow(Py,2));
	h1->Fill(MET_MC, Weight->GetBinContent(Weight->FindBin(hiBin)));
	
	}
 

	TCanvas *c2 = new TCanvas("c2","MET Data/MC(Weight) ID{1,5} |eta|< 2.4 mupt > 50 Gev  ",1900, 1600);
	

	 //Histogram Normalieation 
	Double_t norm = h0->Integral();
	Double_t scale = norm/(h1->Integral());
	h1->Scale(scale);

	
	HistoStyle (h1, 8, 3001);
	h0-> SetMarkerStyle(8);

	/*
	THStack *a = new THStack ("a", "MET Data/MC(Weight) ID{1,5} |eta|< 2.4 mupt > 50 Gev ;#slash{E}_{T}(GeV); Numbre of events");
	  
	a->Add(h1, "HIST");
	a->Add(h0, "PLC PMC");
	a->Draw("nostack");
	*/

	h1->Draw ("SAME"); 
	h0->Draw ("ALP"); 
	
	TLegend *legend = new TLegend(.75,.80,.95,.95);
	legend->AddEntry(h1,"MC");
	legend->AddEntry(h0,"Data");
	legend->Draw("F");

	

		
	c2->SetLogy();
	c2->Print ("meeeee.eps");

}