/////////////////////////////////////////////////////////////////
//      Code for documenting all reweighting factors          ///
//                 Junjie Zhu (July 4, 2011)                  ///
/////////////////////////////////////////////////////////////////

#include "ReweightFactorsManager.h"
#include <iostream>
#include "TFile.h"
#include "TROOT.h"

using namespace std;

/// Global definitions
ReweightFactorsManager* ReweightFactorsManager::_instance=(ReweightFactorsManager*)0;

/// Singleton.  If no instance exists, create one.  Get pointer to instance.
ReweightFactorsManager* ReweightFactorsManager::get_instance() {
  if (_instance == 0) _instance = new ReweightFactorsManager();
  return _instance;
}

/// initialization
void ReweightFactorsManager::initialize() {
  // open root file
  TFile *file_Muon_Trig_Eff = new TFile("SF_for_EPS_Muon.root");
  
  // read histograms
  _MuonTrigEff_EtaPhi_Barrel_Data = (TH2D *)(((TH2D *)(gROOT->FindObject("Eff_Barrel_Data_StacoCB_mu18MG_mu40MSonly_barrel"))))->Clone();
  _MuonTrigEff_EtaPhi_Barrel_MC = (TH2D *)(((TH2D *)(gROOT->FindObject("Eff_Barrel_MC_StacoCB_mu18MG_mu40MSonly_barrel"))))->Clone();
  _MuonTrigEff_Pt_Endcap_Data = (TH1D *)(((TH1D *)(gROOT->FindObject("Eff_Endcap_Data_StacoCB_mu18MG_mu40MSonly_barrel"))))->Clone();
  _MuonTrigEff_Pt_Endcap_MC = (TH1D *)(((TH1D *)(gROOT->FindObject("Eff_Endcap_MC_StacoCB_mu18MG_mu40MSonly_barrel"))))->Clone();

  /// open root file
  TFile *file_Elec_Trig_Eff = new TFile("TrigEff_for_EPS_Elec.root");
  _ElecTrigEff_Eta_MC = (TH1D *)(((TH1D *)(gROOT->FindObject("Trig_Eff_Elec_Eta"))))->Clone();

  /// close file
  ///  file_Muon_Trig_Eff->Close();
  ///  file_Elec_Trig_Eff->Close();
}

/// Constructor initializes pointers to zero
ReweightFactorsManager::ReweightFactorsManager() {
  initialize();
}

/// Destructor 
ReweightFactorsManager::~ReweightFactorsManager() { 
  delete _MuonTrigEff_EtaPhi_Barrel_Data;
  delete _MuonTrigEff_EtaPhi_Barrel_MC;
  delete _MuonTrigEff_Pt_Endcap_Data;
  delete _MuonTrigEff_Pt_Endcap_MC;
  delete _ElecTrigEff_Eta_MC;
}

/// get muon trigger efficiency from data or MC
void ReweightFactorsManager::getMuonTrigEff(bool isData, bool useGeV, double pT, double eta, double phi, double& eff, double& err) {
  Int_t bin = -1;
  if(!useGeV) pT = pT/1000.;
  
  double phi_min = -15./16 * TMath::Pi();
  double phi_max = 17./16 * TMath::Pi();

  if(TMath::Abs(eta) < 1.05) {  // central region
    if(isData) {  
      bin = _MuonTrigEff_EtaPhi_Barrel_Data->FindBin(eta, phi<phi_min?phi+2*TMath::Pi():phi);
      eff = _MuonTrigEff_EtaPhi_Barrel_Data->GetBinContent(bin);
      err = _MuonTrigEff_EtaPhi_Barrel_Data->GetBinError(bin);
    } else {
      bin = _MuonTrigEff_EtaPhi_Barrel_MC->FindBin(eta, phi<phi_min?phi+2*TMath::Pi():phi);
      eff = _MuonTrigEff_EtaPhi_Barrel_MC->GetBinContent(bin);
      err = _MuonTrigEff_EtaPhi_Barrel_MC->GetBinError(bin);
    } // separate for data and MC
  } else {                    // endcap region
    if(isData) {
      if(pT > 120.) pT = 120.;   // use GeV unit
      bin = _MuonTrigEff_Pt_Endcap_Data->FindBin(pT);
      eff = _MuonTrigEff_Pt_Endcap_Data->GetBinContent(bin);
      err = _MuonTrigEff_Pt_Endcap_Data->GetBinError(bin);
    } else {
      bin = _MuonTrigEff_Pt_Endcap_MC->FindBin(pT);
      eff = _MuonTrigEff_Pt_Endcap_MC->GetBinContent(bin);
      err = _MuonTrigEff_Pt_Endcap_MC->GetBinError(bin);
    } // separate for data and MC
  }

  if(fabs(eta) > 2.4) {
    eff = 0.;
    err = 0.;
    cout<<"Muon with eta>2.4, set trigger efficiency to 0"<<endl;
  }     
}

/// get electron trigger efficiency from data or MC
void ReweightFactorsManager::getElecTrigEff(bool isData, bool useGeV, double pT, double eta, double phi, double& eff, double& err) {
  Int_t bin = -1;
  if(!useGeV) pT = pT/1000.;
  
  if(!isData) {   // MC
    bin = _ElecTrigEff_Eta_MC->FindBin(eta);
    eff = _ElecTrigEff_Eta_MC->GetBinContent(bin);
    err = _ElecTrigEff_Eta_MC->GetBinError(bin);
  } else {  // copy the code from https://espace.cern.ch/atlas-phys-sm-ew/wz/Shared%20Documents/TriggerSF/electron_WZ_SF_EPS.h for data so far
    double ele_eta = eta;
    int etaI = -1;
    
    const double etabins[18] = {-2.47,-2.37,-2.01,-1.81,-1.37,-1.15,-0.8,-0.6,-0.1,0,0.1,0.6,0.8,1.15,1.52,1.81,2.01,2.37};
    const double Eff_data_matrix[18] = {0.9607, 0.9690, 0.9793, 0.9818, 0.9920, 0.9847,
					0.9798, 0.9880, 0.9862, 0.9878, 0.9873, 0.9921,
					0.9865, 0.9907, 0.9932, 0.9733, 0.9737, 0.9741};
    const double etamax = 2.47;
    const double etacrack[2] = {1.37,1.52};
    if ( fabs(ele_eta) > etamax || (fabs(ele_eta) > etacrack[0] && fabs(ele_eta) < etacrack[1]) )  // check forward, crack regions
      eff = 1.0;       
    else {           
      for (int i=17; i>=0; i--){    // find eta index
	if ( ele_eta > etabins[i] ) {
	  etaI = i;
	  break;
	}
      }
      eff = Eff_data_matrix[etaI];
    } 
  } // for data situation
}     

/// get event trigger effciency for events with various configurations, 
/// can have several electrons, muons or electrons/muons (such as emumu, eemu etc)
double ReweightFactorsManager::getEventTrigSF(bool useGeV, vector<TLorentzVector> electrons, vector<TLorentzVector> muons) {

  double rate_not_fired_data = 1.;
  double rate_not_fired_mc = 1.;

  for(int ielec = 0; ielec < electrons.size(); ielec ++) {
    double eff_data = 0., err_data = 0.;
    double eff_mc = 0., err_mc = 0.;

    // get efficiency from data
    getElecTrigEff(true, useGeV, electrons[ielec].Pt(), electrons[ielec].Eta(), electrons[ielec].Phi(), eff_data, err_data);

    // get efficiency from MC
    getElecTrigEff(false, useGeV, electrons[ielec].Pt(), electrons[ielec].Eta(), electrons[ielec].Phi(), eff_mc, err_mc);

    rate_not_fired_data *= (1-eff_data);
    rate_not_fired_mc *= (1-eff_mc);
  } // loop over all leptons

  for(int imuon = 0; imuon < muons.size(); imuon ++) {
    double eff_data = 0., err_data = 0.;
    double eff_mc = 0., err_mc = 0.;

    // get efficiency from data
    getMuonTrigEff(true, useGeV, muons[imuon].Pt(), muons[imuon].Eta(), muons[imuon].Phi(), eff_data, err_data);

    // get efficiency from MC
    getMuonTrigEff(false, useGeV, muons[imuon].Pt(), muons[imuon].Eta(), muons[imuon].Phi(), eff_mc, err_mc);

    rate_not_fired_data *= (1-eff_data);
    rate_not_fired_mc *= (1-eff_mc);
  } // loop over all leptons

  double event_SF = 1.;

  // prevent event swith no triggered electrons or muons
  if ((electrons.size() != 0 || muons.size() != 0) && ((1-rate_not_fired_mc) != 0))
    event_SF = (1-rate_not_fired_data)/(1-rate_not_fired_mc);

  return event_SF;
}