/*
Simple macro showing how to access branches from the delphes output root file,
loop over events, and plot simple quantities such as the jet pt and the di-electron invariant
mass.

root -l examples/Example1.C'("delphes_output.root")'
*/

#ifdef __CLING__
R__LOAD_LIBRARY(libDelphes)
#include <iostream>
#include <TClonesArray.h>
#include <TH1.h>
#include <TTree.h>
#include "classes/DelphesClasses.h"
#include "external/ExRootAnalysis/ExRootTreeReader.h"
#endif

//------------------------------------------------------------------------------

void zhbb(const char *inputFile)
{
  gSystem->Load("libDelphes");

  // Create chain of root trees
  TChain chain("Delphes");
  chain.Add(inputFile);

  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();

  // Get pointers to branches used in this analysis
  TClonesArray *branchJet = treeReader->UseBranch("Jet");
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchPhoton   = treeReader->UseBranch("Photon");
  //TClonesArray *branchMET = treeReader->UseBranch("MissingET");

  // Book histograms
  TH1 *histMass = new TH1F("mass", "M_{inv}(e_{1}, e_{2})", 100, 40.0, 140.0);
  TH1 *h3 = new TH1F("h3", "Higgs boson invariant mass", 100,  100.0, 300.0);
  


  // Loop over all events
  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);

    
    Electron *elec1, *elec2;

    // If event contains at least 2 electrons
    if(branchElectron->GetEntries() > 1)
    {
      // Take first two electrons
      elec1 = (Electron *) branchElectron->At(0);
      elec2 = (Electron *) branchElectron->At(1);

      // Plot their invariant mass
      histMass->Fill(((elec1->P4()) + (elec2->P4())).M());
    }

  }

  // define the number of B-jets in the event
  int NofBs = 0;

  // Loop over all events
  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);
 
    Jet *jet1, *jet2;

    // set the number of B-jets in the event to 0
    NofBs = 0;

    // If event contains at least 1 jet
    for(unsigned int i = 0; i < branchJet->GetEntries(); i++)
    {
      // Take first jet
      jet1 = (Jet *) branchJet->At(i);
      jet2 = (Jet *) branchJet->At(i);
     
      if (jet1->BTag == 0) NofBs++;
      if (jet2->BTag == 1) NofBs++;
     
    }

      
      h3 -> Fill(((jet1 -> P4()) + (jet2 -> P4())).M());

  }

  
  // Show resulting histograms
  TCanvas * canvas1 = new TCanvas("canvas1");
  TCanvas * canvas2 = new TCanvas("canvas2");
    
  canvas1->cd();
  histMass->Draw();
  canvas2->cd();
  h3->Draw();
  
  /*TFile *file1 = new TFile("histos.root", "UPDATE");
  histJetPT->Write();
  histMass->Write();
  //histMET->Write();
  histDiPhotMass->Write();
  file1->Close();*/
}