#include "SampleAnalyzer/User/Analyzer/cat_ued_16_013.h"
using namespace MA5;
using namespace std;

// -----------------------------------------------------------------------------
// Initialize
// function called one time at the beginning of the analysis
// -----------------------------------------------------------------------------
bool cat_ued_16_013::Initialize(const MA5::Configuration& cfg, const std::map<std::string,std::string>& parameters)
{
  cout << "BEGIN Initialization" << endl;
  // initialize variables, histos
  // Initialize signal regions
  Manager()->AddRegionSelection("200<MET<500,100<HT<500"); // SR1
  Manager()->AddRegionSelection("MET>500,100<HT<500"); // SR2
  Manager()->AddRegionSelection("200<MET<500,HT>500"); // SR3
  Manager()->AddRegionSelection("MET>500,HT>500"); // SR4


  // Baseline cuts
  Manager()->AddCut("nolepton");
  Manager()->AddCut("MET>200");
  Manager()->AddCut("jet_pT>20");
  Manager()->AddCut("HT>100");

  // Signal-region specific cuts
  // MET < 500 for signal regions 1 and 3:
  string SR_1_3[] = {"200<MET<500,100<HT<500", "200<MET<500,HT>500"};
  Manager()->AddCut("MET<500", SR_1_3);

  // HT < 500 applies to signal regions 1 and 2:
  string SR_1_2[] = {"200<MET<500,100<HT<500", "MET>500,100<HT<500"};
  Manager()->AddCut("HT<500", SR_1_2);

  // HT > 500 applies to signal regions 3 and 4:
  string SR_3_4[] = {"200<MET<500,HT>500", "MET>500,HT>500"};
  Manager()->AddCut("HT>500", SR_3_4);

  // MET > 500 for signal regions 2 and 4:
  string SR_2_4[] = {"MET>500,100<HT<500", "MET>500,HT>500"};
  Manager()->AddCut("MET>500", SR_2_4);


  cout << "END   Initialization" << endl;
  return true;
}

// -----------------------------------------------------------------------------
// Finalize
// function called one time at the end of the analysis
// -----------------------------------------------------------------------------
void cat_ued_16_013::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files)
{
  cout << "BEGIN Finalization" << endl;
  // saving histos
  cout << "END   Finalization" << endl;
}

// -----------------------------------------------------------------------------
// Execute
// function called each time one event is read
// -----------------------------------------------------------------------------
bool cat_ued_16_013::Execute(SampleFormat& sample, const EventFormat& event)
{
  // ***************************************************************************
  // Example of analysis with reconstructed particles
  // Concerned samples : ROOT
  // ***************************************************************************

  double myEventWeight;
  if(Configuration().IsNoEventWeight()) myEventWeight=1.;
  else if(event.mc()->weight()!=0.) myEventWeight=event.mc()->weight();
  else
    {
      INFO << "Found one event with a zero weight. Skipping...\n";
      return false;
    }
  Manager()->InitializeForNewEvent(myEventWeight);

  if (event.rec()!=0)
    {
    cout << "---------------NEW EVENT-------------------" << endl;
    // =====first declare the empty containers:=====¬
    vector<const RecLeptonFormat*> electron, muon;
    vector<const RecJetFormat*> jet;

    // =====fill the electrons container:=====
    for(unsigned int e=0; e<event.rec()->electrons().size(); e++)
    {
      const RecLeptonFormat *CurrentElectron = &(event.rec()->electrons()[e]);
      double pt = CurrentElectron->momentum().Pt();
      if(!(pt>10)) continue;
      electron.push_back(CurrentElectron);
    }

    // =====fill the muons container:=====
    for(unsigned int m=0; m<event.rec()->muons().size(); m++)
    {
      const RecLeptonFormat *CurrentMuon = &(event.rec()->muons()[m]);
      double pt = CurrentMuon->momentum().Pt();
      if(!(pt>10)) continue;
      muon.push_back(CurrentMuon);
    }

    // =====define HT======
    double HT = 0;
    // =====fill the jet containers and order=====
    for(unsigned int j=0; j<event.rec()->jets().size(); j++)
    {
      const RecJetFormat *CurrentJet = &(event.rec()->jets()[j]);
      double pt = CurrentJet->momentum().Pt();
      if(!(pt > 20)) continue;
      // Add pt to HT:
      HT += pt;
      jet.push_back(CurrentJet);
    }

    // =====Get the missing ET=====
    TLorentzVector pTmiss = event.rec()->MET().momentum();
    double MET = pTmiss.Pt();

    // =====Apply no-lepton cut========
    if(!Manager()->ApplyCut((electron.size()+muon.size() == 0),"nolepton")) return true;

    // =====Apply HT100 cut============
    if(!Manager()->ApplyCut((HT > 100),"HT>100")) return true;

    // =====Apply MET200 cut===========
    if(!Manager()->ApplyCut((MET > 200),"MET>200")) return true;

    // =====Apply signal region cuts=======
    Manager()->ApplyCut(MET < 500, "MET<500"); // SR1, SR3
    Manager()->ApplyCut(HT < 500, "HT<500"); // SR1, SR2
    Manager()->ApplyCut(HT >= 500, "HT>500"); // SR3, SR4
    Manager()->ApplyCut(MET >= 500, "MET>500"); // SR2, SR4

    // =====finished=====
    return true;
  }
}