#define CutAndCount_cxx
// The class definition in MyAnalysis.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("MyAnalysis.C")
// Root > T->Process("MyAnalysis.C","some options")
// Root > T->Process("MyAnalysis.C+")
//

#include "CutAndCount.h"
#include <iostream>
#include <TH1F.h>
#include <TLatex.h>

using namespace std;

void CutAndCount::BuildEvent() {
   
  Muons.clear();
  for (int i = 0; i < NMuon; ++i) {
     MyMuon muon(Muon_Px[i], Muon_Py[i], Muon_Pz[i], Muon_E[i]);
      muon.SetIsolation(Muon_Iso[i]);
      muon.SetCharge(Muon_Charge[i]);
      Muons.push_back(muon);
  }
   
   Electrons.clear();
   for (int i = 0; i < NElectron; ++i) {
      MyElectron electron(Electron_Px[i], Electron_Py[i], Electron_Pz[i], Electron_E[i]);
      electron.SetIsolation(Electron_Iso[i]);
      electron.SetCharge(Electron_Charge[i]);
      Electrons.push_back(electron);
   }
   
   Photons.clear();
   for (int i = 0; i < NPhoton; ++i) {
      MyPhoton photon(Photon_Px[i], Photon_Py[i], Photon_Pz[i], Photon_E[i]);
      photon.SetIsolation(Photon_Iso[i]);
      Photons.push_back(photon);
   }
   
   Jets.clear();
   for (int i = 0; i < NJet; ++i) {
      MyJet jet(Jet_Px[i], Jet_Py[i], Jet_Pz[i], Jet_E[i]);
      jet.SetBTagDiscriminator(Jet_btag[i]);
      jet.SetJetID(Jet_ID[i]);
      Jets.push_back(jet);
   }
   
   hadB.SetXYZM(MChadronicBottom_px, MChadronicBottom_py, MChadronicBottom_pz, 4.8);
   lepB.SetXYZM(MCleptonicBottom_px, MCleptonicBottom_py, MCleptonicBottom_pz, 4.8);
   hadWq.SetXYZM(MChadronicWDecayQuark_px, MChadronicWDecayQuark_py, MChadronicWDecayQuark_pz, 0.0);
   hadWqb.SetXYZM(MChadronicWDecayQuarkBar_px, MChadronicWDecayQuarkBar_py, MChadronicWDecayQuarkBar_pz, 0.0);
   lepWl.SetXYZM(MClepton_px, MClepton_py, MClepton_pz, 0.0);
   lepWn.SetXYZM(MCneutrino_px, MCneutrino_py, MCneutrino_pz, 0.0);
   met.SetXYZM(MET_px, MET_py, 0., 0.);
   
   EventWeight *= weight_factor;
   for(int j=0;j<EtaLow.size();j++)NW[j]=0;
   
}

void CutAndCount::Begin(TTree * /*tree*/) {
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).
   
   TString option = GetOption();
   
}

void CutAndCount::SetEtaBounds(vector<float> etabounds){

  EtaLow.clear();
  EtaHigh.clear();
  
  for (int i=0;i<etabounds.size()-1;i++){
    EtaLow.push_back(etabounds[i]);
    EtaHigh.push_back(etabounds[i+1]);
  }
    
}

void CutAndCount::SlaveBegin(TTree * /*tree*/) {

   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).
   
   TString option = GetOption();
   
   Asymetry.clear();
   AsymetryError.clear();
   
   EtaLow.clear();
   EtaHigh.clear();

   // default eta bounds
   vector<float> etabounds;
   etabounds.push_back(0);
   etabounds.push_back(0.4);
   etabounds.push_back(0.8);
   etabounds.push_back(1.2);
   etabounds.push_back(1.5);
   etabounds.push_back(1.8);
   etabounds.push_back(2.1);
   
   SetEtaBounds(etabounds);

}

Bool_t CutAndCount::Process(Long64_t entry) {
   
   ++TotalEvents;
   
   GetEntry(entry);
   
   if (TotalEvents % 10000 == 0)
      cout << "Next event -----> " << TotalEvents << endl;
   
   BuildEvent();
   
   //////////////////////////////
   // Exercise 2: Cut And Count

   //  double MuonPtCut = 25.;
   double MuonPtCut10 = 10.;
   double MuonRelIsoCut = 0.15;
   //   double MuonRelIsoCut = 0.01;
   //double METCut=30;

   if(!triggerIsoMu24) return kFALSE;

   int N_IsoMuon10 = 0;
   int N_Muon10 = 0;

   int N_IsoMuon25 = 0;
   int N_Muon25 = 0;

   MyMuon *muoniso1, *muoniso2;
   MyMuon *muon1, *muon2;

   // global distributions
   
   for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
     
     if (jt->Pt()>MuonPtCut10) {
       ++N_Muon10;     

       if (N_Muon10 == 1) muon1 = &(*jt);
       if (N_Muon10 == 2) muon2 = &(*jt);

       if (jt->IsIsolated(MuonRelIsoCut)){
	 ++N_IsoMuon10;
       
	 if (N_IsoMuon10 == 1) muoniso1 = &(*jt);
	 if (N_IsoMuon10 == 2) muoniso2 = &(*jt);
	 
       }
       
     }
     
     if (jt->Pt()>MuonPtCut) {
       
       ++N_Muon25;
       
       
       if (jt->IsIsolated(MuonRelIsoCut)) {
	 ++N_IsoMuon25;
       }
     }
   }
   // for W
   if(doW){ 
     if (N_IsoMuon25 == 1 && N_IsoMuon10<2 && met.Pt()>METCut) { 
       
       int ieta=999;
       
       for(int j=0;j<EtaHigh.size();j++){
	 if(fabs(muoniso1->Eta())<EtaHigh[j] && fabs(muoniso1->Eta())>EtaLow[j]) ieta=j;
       }	 

      
       if(muoniso1->GetCharge()==1) Nplus[ieta]+=EventWeight;
       if(muoniso1->GetCharge()==-1) Nminus[ieta]+=EventWeight;    

       // Count also the number of events for Z frac in fit:
       NW[ieta]+=EventWeight;
       
     }
   }else{
   // for Z  
     if (N_IsoMuon25==2 && (*muoniso1 + *muoniso2).M()>70 &&  (*muoniso1 + *muoniso2).M()<110 ){
     
       int ieta1=999;
       int ieta2=999;
       
       for(int j=0;j<EtaHigh.size();j++){
	 if(fabs(muoniso1->Eta())<EtaHigh[j] && fabs(muoniso1->Eta())>EtaLow[j]) ieta1=j;
	 if(fabs(muoniso2->Eta())<EtaHigh[j] && fabs(muoniso2->Eta())>EtaLow[j]) ieta2=j;
       }
       if(muoniso1->GetCharge()==1) Nplus[ieta1]++;
       if(muoniso1->GetCharge()==-1) Nminus[ieta1]++;    
       if(muoniso2->GetCharge()==1) Nplus[ieta2]++;
       if(muoniso2->GetCharge()==-1) Nminus[ieta2]++;            
     }
     

   }
   
   //////////////////////////////
   
   return kTRUE;
}

void CutAndCount::SlaveTerminate() {
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server. float Asym[10];
  float Asym[10];
  float AsymError[10];
  
  for(int j=0;j<EtaHigh.size();j++){
    Asym[j]=(Nplus[j]-Nminus[j])/(Nplus[j]+Nminus[j]);
    AsymError[j]=2.*sqrt( (Nplus[j]*Nminus[j]*Nminus[j] + Nminus[j]*Nplus[j]*Nplus[j]) / TMath::Power(Nplus[j]+Nminus[j],4));
    Asymetry.push_back(Asym[j]);
    AsymetryError.push_back(AsymError[j]);
  }
  
  for(int j=0;j<EtaHigh.size();j++){
    nWvec.push_back(NW[j]);
    nWplus.push_back(Nplus[j]);
    nWminus.push_back(Nminus[j]);
  }
   
}

void CutAndCount::Terminate() {
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
 
}