#include "MET_Cleaning_Utils.h"
#include <iostream>

// Constructor
MET_Cleaning_Utils::MET_Cleaning_Utils():
  GeV(1000)
  ,jet_THRESHOLD(20*GeV)
  ,jet_THRESHOLD_1sigma(4*GeV)
  ,jet_iso(0.1)
{
  eta_min.push_back(0.0);eta_max.push_back(1.45);phi_min.push_back(-0.78847);phi_max.push_back(-0.59213);
}

//destructor
MET_Cleaning_Utils::~MET_Cleaning_Utils(){

}

// Common funtions
float MET_Cleaning_Utils::DeltaPhi(float phi1, float phi2){
  float dphi=fabs(phi1-phi2);
  if (dphi>TMath::Pi())   dphi=2*TMath::Pi()-dphi;
  return dphi;
}
float MET_Cleaning_Utils::DeltaEta(float eta1, float eta2){
  float deta=fabs(eta1-eta2);
  return deta;
}
float MET_Cleaning_Utils::dr(float phi1, float eta1,float phi2, float eta2){
  return sqrt(pow(DeltaEta(eta1,eta2),2.0)+DeltaPhi(phi1,phi2));
}

float MET_Cleaning_Utils::dphi_MET_MPT(float &met_phi,std::vector<float> *trk_pt,std::vector<float> *trk_eta,
				       std::vector<float> *trk_d0_wrtPV,std::vector<float> *trk_z0_wrtPV,
				       std::vector<float> *trk_theta_wrtPV,std::vector<float> *trk_phi_wrtPV,
				       std::vector<int> *trk_nPixHits,std::vector<int> *trk_nSCTHits){
  float sign=1;
  float mpt_phi=MPT_Loose_phi(trk_pt,trk_eta,trk_d0_wrtPV,trk_z0_wrtPV,trk_theta_wrtPV,trk_phi_wrtPV,trk_nPixHits,trk_nSCTHits);
  if(met_phi-mpt_phi<0) sign=-1;
  return DeltaPhi(met_phi,mpt_phi)*sign;
}

float MET_Cleaning_Utils::MPT_Loose_phi(std::vector<float> *trk_pt,std::vector<float> *trk_eta,std::vector<float> *trk_d0_wrtPV,
					std::vector<float> *trk_z0_wrtPV,std::vector<float> *trk_theta_wrtPV,
					std::vector<float> *trk_phi_wrtPV,std::vector<int> *trk_nPixHits,
					std::vector<int> *trk_nSCTHits ){
  TLorentzVector P(0,0,0,0);
  for(int i=0; i<trk_pt->size();i++){
    if(trk_pt->at(i)>500 &&  fabs(trk_eta->at(i))<2.5 && fabs(trk_d0_wrtPV->at(i))<1.5 &&
       fabs(trk_z0_wrtPV->at(i)*trk_theta_wrtPV->at(i))<1.5 && trk_nPixHits->at(i)>=1 &&
       trk_nSCTHits->at(i)>=6 ){
      TLorentzVector LT(0,0,0,0);
      LT.SetPtEtaPhiM(trk_pt->at(i),0,trk_phi_wrtPV->at(i),0);
      P-=LT;
    }
  }
  return P.Phi();
}



bool MET_Cleaning_Utils::iswithinDeadFEB(float eta, float phi, float objCone, float objIso,bool dosquare){
  for(int i=0;i<eta_min.size();i++){
    //std::cout << "iswithinDeadFEB" << std::endl;
    // compute if the jet cone overlaps with a rectangular dead feb region
    //step 1 compute rectangular cut
    float etamid=(eta_max.at(i)+eta_min.at(i))/2;
    float phimid=(phi_max.at(i)+phi_min.at(i))/2;
    float etadiff=fabs(eta_max.at(i)-eta_min.at(i))/2+objIso;
    float phidiff=fabs(DeltaPhi(phi_max.at(i),phi_min.at(i)))/2+objIso;
    float dphi=fabs(DeltaPhi(phi,phimid));
    float deta=fabs(eta-etamid);
    if(fabs(etadiff+objCone)<fabs(deta) || fabs(phidiff+objCone)<fabs(dphi))break;
    if(dosquare) return true;
    // step to now look at corners to take into accont that you have a square and a circle
    if(fabs(etadiff+objCone)>fabs(deta) && fabs(phidiff)>fabs(dphi)) return true;
    if(fabs(etadiff)>deta && fabs(phidiff+objCone)>fabs(dphi)) return true;
    if(sqrt(pow(deta-etadiff,2.0)+pow(dphi-phidiff,2.0)) < objCone)return true;
  }
  return false;
}


bool MET_Cleaning_Utils::isJetAffected(std::vector<TLorentzVector> &jet, float iso, float JCone, float minpt ){
  for(int i=0;i<jet.size();i++){
    if(isJetAffected(jet.at(i),iso,JCone,minpt))return true;
  }
  return false;
}

bool MET_Cleaning_Utils::isJetAffected(TLorentzVector &jet, float iso, float JCone, float minpt ){
  if(jet.Pt()>minpt){
    if(iswithinDeadFEB(jet.Eta(),jet.Phi(),JCone,iso)){
      return true;
    }
  }
  return false;
}


bool MET_Cleaning_Utils::hasUnmatchTrackJet(std::vector<TLorentzVector> &jet,float JCone,std::vector<TLorentzVector> &trackjet,
					    float TJCone,float minpt){
  for(int i=0;i<trackjet.size();i++){
    float dr_min=10;
    if(isJetAffected(trackjet.at(i),0.0,(JCone+TJCone)/4,minpt)){
      for(int j=0; j<jet.size();j++){
        float currentdr=dr(trackjet.at(i).Phi(),trackjet.at(i).Eta(),
                           jet.at(j).Phi(),jet.at(j).Eta());
        if(currentdr<dr_min && currentdr<JCone+TJCone){
          dr_min=currentdr;
        }
      }
      if((JCone+TJCone)/4<dr_min) return true;
    }
  }
  return false;
}



float MET_Cleaning_Utils::JetTrackJetRatio(std::vector<TLorentzVector> &jet,TLorentzVector &trackjet){
  float JCone=0.4;
  float TJCone=0.4;
  float dr_min=999;
  int jidx=-1;
  for(int j=0; j<jet.size();j++){
    float currentdr=dr(trackjet.Phi(),trackjet.Eta(),
		       jet.at(j).Phi(),jet.at(j).Eta());
    if(currentdr<dr_min && currentdr<0.1){
      dr_min=currentdr;
        jidx=j;
    }
  }
  if(jidx!=-1) return jet.at(jidx).Pt()/trackjet.Pt();
  return -1.0;
}


bool MET_Cleaning_Utils::isJetSafe(std::vector<TLorentzVector> &jet, std::vector<float> *BCH_CORR_JET, 
				   std::vector<float> *BCH_CORR_CELL, bool isdata, float shift){
  if(jet.size()!=BCH_CORR_CELL->size() || jet.size()!=BCH_CORR_JET->size()){
    std::cout << "FATAL: jets BCH_CORR_CELL and BCH_CORR_JET do not match! FIX YOU INPUTS  " << std::endl;
  }
  for(int i=0;i<jet.size();i++){
    float THRESHOLD=jet_THRESHOLD+jet_THRESHOLD_1sigma*shift;
    if(isdata)THRESHOLD*=(1-BCH_CORR_JET->at(i))/(1-BCH_CORR_CELL->at(i));
    if(isJetAffected(jet.at(i),0.0,jet_iso,THRESHOLD))return false;
  }
  return true;
}


float MET_Cleaning_Utils::dphi_Jet_MET(std::vector<TLorentzVector> &jet, float met_phi){
  float dphi_min=10;
  for(int i=0;i<jet.size();i++){
    float Dphi= DeltaPhi(jet.at(i).Phi(),met_phi);
    if(dphi_min>Dphi){
      dphi_min=Dphi;
    }
  }
  return dphi_min;
}