/****************************************************************************/
/**  This macro extracts variables from Delphes output to a new ROOT file  **/
/****************************************************************************/

#include <iostream>
#include <TClonesArray.h>
#include <TTree.h>
#ifdef __CLING__
R__LOAD_LIBRARY(libDelphes)
#include "classes/DelphesClasses.h"
#include "external/ExRootAnalysis/ExRootTreeReader.h"
#endif
 #define PI 3.1415927
//--------------------------------------------------------------------------------------
using namespace std;
/*******************************************************************************************/
void Zextract( const char *inputFile, const char *outputFile)
{
  gSystem->Load("libDelphes");

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

  // Creating the output file and tree
  TFile *outfile = new TFile(outputFile,"CREATE");
  TTree *outtree=new TTree("event_tree","Test");

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

//  TClonesArray *branchParticle = treeReader->UseBranch("Particle");
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchMuon = treeReader->UseBranch("Muon");
  TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");

  //Variables in new ROOT file
Double_t MET=0.0, lepPt=0.0;

//Creating Branches in new ROOT file
  outtree->Branch("lepPt",&lepPt);
  outtree->Branch("MissingET",&MET);

/******************* Cut are defined here **************************/
Double_t lep_pt_pre = 20.; //		  **************************/
Double_t lep_eta_pre = 2.45; //		  **************************/
/*******************************************************************/

   Electron *elec[50];
   Muon *mu[50];
   MissingET *misset;

   //Loop over all events
   for(Int_t entry = 0; entry < numberOfEntries; entry++)
   {
   treeReader->ReadEntry(entry);

Int_t elec_no(0), mu_no(0), miss_no(0) ;

    for(Int_t j=0;j<branchMissingET->GetEntries();j++)   //Missing energy branch
    {
     misset  = (MissingET *) branchMissingET->At(j);
    }

    for(Int_t j=0;j<branchElectron->GetEntries();j++)   //Electron branch
    {
      elec[j] = (Electron *) branchElectron->At(j);
if(elec[j]->PT > lep_pt_pre && abs(elec[j]->Eta) < lep_eta_pre){elec_no++; }
    }
    for(Int_t j=0;j<branchMuon->GetEntries();j++)   //Muon branch
    {
      mu[j] = (Muon *) branchMuon->At(j);
if((mu[j]->PT) > lep_pt_pre && abs(mu[j]->Eta) < lep_eta_pre){mu_no++; }
    }

if(elec_no >= 2 && mu_no <= 2)
{
      MET = misset->MET;
      lepPt = (elec[0]->P4() + elec[1]->P4()).Pt();
outtree->Fill();
}else if( mu_no >= 2 && elec_no < 2 ){
      MET = misset->MET;
      lepPt = (mu[0]->P4() + mu[1]->P4()).Pt();
outtree->Fill();
}
}//Event loop

  outfile->Write();
delete outtree;
delete outfile;
}