#ifdef __CLING__
// This part is enable only in ROOT6.
R__LOAD_LIBRARY(libDelphes)
#include "classes/DelphesClasses.h"
#include "external/ExRootAnalysis/ExRootTreeReader.h"
#endif

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

// #define BRANCH_SIZE 4

// ## means concatenation as a symbol.
#define SUM(c, ht) for (int i = 0; i < branch##c->GetEntries(); i++) ht += ((c*) branch##c->At(i))->PT;

void meff2(const char *inputFile)
{
  gSystem->Load("~/Delphes-3.3.0/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 *branchMuon = treeReader->UseBranch("Muon");
  TClonesArray *branchMET = treeReader->UseBranch("MissingET");
  TClonesArray *branchPhoton = treeReader->UseBranch("Photon");

  /*
  TClonesArray *branches[BRANCH_SIZE];

  branches[0] = branchJet;
  branches[1] = branchElectron;
  branches[2] = branchMuon;
  branches[3] = branchPhoton;
  */

  // Book histograms
  TH2 *histMETMeff = new TH2F("MET-Meff", "MET-H_{T}", 100, 0.0, 1500.0, 100, 0.0, 3500.0);
  histMETMeff->GetXaxis()->SetTitle("Missing ET");
  histMETMeff->GetYaxis()->SetTitle("Scalar HT");

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

    Double_t ht = 0.;

    /*
    for (int i = 0; i < BRANCH_SIZE; i++) {
      TClonesArray* branch = branches[i];
      for (int j = 0; j < branch->GetEntries(); j++) {
        // This code is terrible since branch->At(j) may return
        // a non-Jet instance, like Electron or Muon.
        // Then, why does this work?
        // Roughly speaking, the reason is the memory alignment 
        // of each classes. We can find it in the header file.
        // For details, please ask me.
        // You **must not** write this kind of code unless you 
        // understand the reason.
        ht += ((Jet*) branch->At(j))->PT; 
        // a += b is equivalent to a = a + b.
        // Similarly, we can use a *= b, a /= b, and something like that.
      }
    }
    */

    // This macro should be expanded as
    // for (int i = 0; i < branchJet->GetEntries(); i++) ht += ((Jet*) branchJet->At(i)).PT;
    SUM(Jet, ht);
    // The followings are similar.
    SUM(Electron, ht);
    SUM(Muon, ht);
    SUM(Photon, ht);
    // Anyway, I do NOT want to copy and paste. That shuold be ugly.

    Double_t meff = ht;
    Double_t met = 0.;

    for (int i = 0; i < branchMET->GetEntries(); i++) {
      Double_t metv = ((MissingET*) branchMET->At(i))->MET;
      meff += metv;
      met += metv;
    }

    histMETMeff->Fill(met, meff);

  }

  TCanvas *c1=new TCanvas("c1", "test", 800, 600); 
  // Color scheme
  histMETMeff->Draw("colz");
  // Saving the output as PDF
  c1->Print("met-meff.pdf"); 
}