// Based on the macro examples/Example1.C provided by Delphes
// 21 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 macro2.cc

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

// Jetset, JADE and Durham tools at ftp://ftp.slac.stanford.edu/groups/lcd/Physics_tools/ (M. Iwasaki)
#include "JetFinder.cxx" 
#include "EventShape.cxx" 

void macro2() {

  gSystem->Load("libDelphes");

  EventShape* eshape=new EventShape();
  JetFinder* jetfinder=new JetFinder(0.004); 

  // 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);

  TH1F *histThrust=new TH1F("Jetset Thrust", "Jetset Thrust", 100, 0.0, 1.0);
  TH1F *histDurhamY23=new TH1F("Durham y_{23}", "Durham y_{23}", 100, 0.0, 0.1);
  TH1F *histProcessID=new TH1F("Process ID", "Process ID", 40000, 0, 40000);
  TH1F *histMissingP3=new TH1F("Missing 3D Momentum", "Missing 3D Momentum", 250, 0, cme);
  TH1F *histJetMultiplicity=new TH1F("Durham y_{cut}=0.004 Jet Multiplicity", "Durham y_{cut}=0.004 Jet Multiplicity", 10, 0, 10);
  TH1F *histTrackMultiplicity=new TH1F("Track Multiplicity", "Track Multiplicity", 20, 0.0, 20.);

  TChain chain("Delphes");
  chain.Add(inputFile);

  ExRootTreeReader *treeReader=new ExRootTreeReader(&chain);
  TClonesArray *branchEvent=treeReader->UseBranch("Event");
  TClonesArray *branchTrack=treeReader->UseBranch("Track");
  TClonesArray *branchElectron=treeReader->UseBranch("Electron");
  TClonesArray *branchPhoton=treeReader->UseBranch("Photon");
  TClonesArray *branchChargedHadron=treeReader->UseBranch("ChargedHadron");
  TClonesArray *branchNeutralHadron=treeReader->UseBranch("NeutralHadron");
  TClonesArray *branchMuon=treeReader->UseBranch("Muon");
  TClonesArray *branchMissingET=treeReader->UseBranch("MissingET");

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

    if(entry%100==0) cout << "Event " << entry << endl;
    treeReader->ReadEntry(entry);

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

    MissingET* MET=(MissingET*) branchMissingET->At(0);
    Double_t met=MET->MET;
    Double_t eta=MET->Eta;
    TVector3* missingTV3=new TVector3(0,0,0);
    missingTV3->SetPtEtaPhi(MET->MET, MET->Eta, MET->Phi);

    TObjArray* jetpartTLV=new TObjArray();

    TObjArray* tracksTLV=new TObjArray(); 
    TObjArray* tracksTV3=new TObjArray(); 
    for(Int_t iTrack=0; iTrack<branchTrack->GetEntries(); iTrack++) {
      Track* object=(Track*) branchTrack->At(iTrack);
      TLorentzVector* object4V=new TLorentzVector(object->P4());
      TVector3* object3V=new TVector3(object4V->Vect());
      tracksTV3->Add(object3V);
      tracksTLV->Add(object4V);
      jetpartTLV->Add(object4V);
    }

    eshape->setPartList(tracksTV3); 
    Double_t thrust=float(eshape->thrust().X());

    TObjArray* eplusTLV=new TObjArray(); 
    TObjArray* eminusTLV=new TObjArray(); 
    for(Int_t iElectron=0; iElectron<branchElectron->GetEntries(); iElectron++) {
      Electron* object=(Electron*) branchElectron->At(iElectron);
      TLorentzVector* object4V=new TLorentzVector(object->P4());
      if(object->Charge>0) eplusTLV->Add(object4V);
      if(object->Charge<0) eminusTLV->Add(object4V);
    }

    TObjArray* photonsTLV=new TObjArray(); 
    for(Int_t iPhoton=0; iPhoton<branchPhoton->GetEntries(); iPhoton++) {
      Photon* object=(Photon*) branchPhoton->At(iPhoton);
      TLorentzVector* object4V=new TLorentzVector(object->P4());
      photonsTLV->Add(object4V);
      jetpartTLV->Add(object4V);
    }

    TObjArray* hplusTLV=new TObjArray(); 
    TObjArray* hminusTLV=new TObjArray(); 
    for(Int_t iChargedHadron=0; iChargedHadron<branchChargedHadron->GetEntries(); iChargedHadron++) {
      Track* object=(Track*) branchChargedHadron->At(iChargedHadron);
      TLorentzVector* object4V=new TLorentzVector(object->P4());
      if(object->Charge>0) hplusTLV->Add(object4V);
      if(object->Charge<0) hminusTLV->Add(object4V);
    }

    TObjArray* nhadronsTLV=new TObjArray(); 
    for(Int_t iNeutralHadron=0; iNeutralHadron<branchNeutralHadron->GetEntries(); iNeutralHadron++) {
      Tower* object=(Tower*) branchNeutralHadron->At(iNeutralHadron);
      TLorentzVector* object4V=new TLorentzVector(object->P4());
      nhadronsTLV->Add(object4V);
      jetpartTLV->Add(object4V);
    }

    TObjArray* muplusTLV=new TObjArray(); 
    TObjArray* muminusTLV=new TObjArray(); 
    for(Int_t iMuon=0; iMuon<branchMuon->GetEntries(); iMuon++) {
      Muon* object=(Muon*) branchMuon->At(iMuon);
      TLorentzVector* object4V=new TLorentzVector(object->P4());
      if(object->Charge>0) muplusTLV->Add(object4V);
      if(object->Charge<0) muminusTLV->Add(object4V);
    }

    TObjArray* durhamjetTLV=new TObjArray(); 
    jetfinder->setDURHAM();
    jetfinder->setYCut(0.004);
    jetfinder->setPartList(jetpartTLV);
    jetfinder->doFindJets();
    Int_t jetMultiplicity=jetfinder->njets();
    for(Int_t iJet=0; iJet<jetMultiplicity; iJet++) durhamjetTLV->Add(jetfinder->jet4vec(iJet));
        

    Int_t next, last=0;
    Double_t ycut=0.100;
    Double_t y23=0;
    while(ycut>0&&y23==0.) {
      jetfinder->setYCut(double(ycut));
      jetfinder->doFindJets();
      next=jetfinder->njets();
      if(last==2&&next==3) y23=ycut;
      ycut-=0.001;
      last=next;
    }

    histThrust->Fill(thrust, eventWeight);
    histDurhamY23->Fill(y23, eventWeight);
    histProcessID->Fill(processID, eventWeight);
    histMissingP3->Fill(missingTV3->Mag(), eventWeight);
    histJetMultiplicity->Fill(jetMultiplicity, eventWeight);
    histTrackMultiplicity->Fill(tracksTLV->GetEntries(), eventWeight);

  }

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

  c0->cd(1);
  histProcessID->Draw("HIST");
  c0->cd(2);
  histDurhamY23->Draw("HIST");
  c0->cd(3);
  histThrust->Draw("HIST");
  c0->cd(4);
  histMissingP3->Draw("HIST");
  c0->cd(5);
  histJetMultiplicity->Draw("HIST");
  c0->cd(6);
  histTrackMultiplicity->Draw("HIST");

  c0->SaveAs("macro2.eps");
}