// Based on the macro examples/Example1.C provided by Delphes
// 20 July 2020, Chris Potter. Usage:
// 1) place the this macro in your Delphes installation directory
// 2) modify inputFile, normalization and cme as appropriate
// 3) bash> root -b -q macro1.cc

#ifdef __CLING__
R__LOAD_LIBRARY(libDelphes)
#include "classes/DelphesClasses.h"
#include "external/ExRootAnalysis/ExRootTreeReader.h"
#endif

void macro1() {

  gSystem->Load("libDelphes");

  // Modify the following lines as appropriate!
  TString inputFile("potter/rootFiles/zllhinv_ILC250_LR.root");
  Double_t luminosity=250.; // fb-1
  Double_t xsec=310.*0.101*0.1; // fb
  Double_t normalization=xsec*luminosity;
  //TString inputFile("potter/rootFiles/all_SM_background_ILC250_LR.root");
  //Double_t normalization=1.;
  Double_t cme=250.;
  TLorentzVector beamsP4(0,0,0,cme);

  TChain chain("Delphes");
  chain.Add(inputFile);
  ExRootTreeReader *treeReader=new ExRootTreeReader(&chain);
  TClonesArray *branchEvent=treeReader->UseBranch("Event");
  TClonesArray *branchElectron=treeReader->UseBranch("Electron");

  TH1 *histRecoilMass=new TH1F("RecoilMass", "m_{recoil}", 200, 0.0, 200.0);
  TH1 *histZMass=new TH1F("ZMass", "m_{ee}", 200, 0.0, 200.0);

  for(Int_t entry=0; entry<treeReader->GetEntries(); ++entry) {

    treeReader->ReadEntry(entry);
    LHEFEvent* event=(LHEFEvent*) branchEvent->At(0);
    Double_t eventWeight=event->Weight;
    Double_t processID=event->ProcessID;

    Electron *electron1, *electron2;
    if(branchElectron->GetEntries()>1) {
      electron1=(Electron*) branchElectron->At(0);
      electron2=(Electron*) branchElectron->At(1);
      TLorentzVector ZP4=electron1->P4()+electron2->P4();
      if(TMath::Abs(ZP4.M()-91.1876)<5.) {
        cout << "Process ID: " << processID << ", Event Weight: " << eventWeight << endl;
        histZMass->Fill(ZP4.M(), eventWeight);
        TLorentzVector recoilP4=beamsP4-ZP4;
        histRecoilMass->Fill(recoilP4.M(), eventWeight);
      }
    }
  }

  if(normalization!=1.) {
    histZMass->Scale(normalization/histZMass->GetEntries());
    histRecoilMass->Scale(normalization/histZMass->GetEntries());
  }

  TCanvas* c0=new TCanvas("c0", "c0", 500,300);
  c0->cd(1);
  c0->Divide(2,1);

  c0->cd(1);
  histZMass->Draw("HIST");
  c0->cd(2);
  histRecoilMass->Draw("HIST");

  c0->SaveAs("macro.eps");

}