#include <stdio.h>
#include <stdlib.h>
#include <sstream>
#include <iomanip>
#include <cmath>
#include <iostream>
#include <fstream>
#include <math.h>
#include <string.h>
#include <exception>
#include <cstdlib>
#include <cstdio>
#include <limits>


#include "TROOT.h"
#include "TSystem.h"
#include "TApplication.h"
#include "TGClient.h"
#include "TString.h"
#include "TClonesArray.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH1I.h"
#include "TH1D.h"
#include "TH2F.h"
#include "TH2I.h"
#include "TH2D.h"
#include "THStack.h"
#include "TLorentzVector.h"

#include "/home/sayan/HEPTOOLS/Delphes-3.4.1/external/ExRootAnalysis/ExRootTreeReader.h"
#include "/home/sayan/HEPTOOLS/Delphes-3.4.1/external/ExRootAnalysis/ExRootTreeWriter.h"
#include "/home/sayan/HEPTOOLS/Delphes-3.4.1/external/ExRootAnalysis/ExRootTreeBranch.h"
#include "/home/sayan/HEPTOOLS/Delphes-3.4.1/external/ExRootAnalysis/ExRootResult.h"
#include "/home/sayan/HEPTOOLS/Delphes-3.4.1/external/ExRootAnalysis/ExRootUtilities.h"


#include "/home/sayan/HEPTOOLS/Delphes-3.4.1/classes/DelphesClasses.h"

using namespace std;



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

int main(int argc, char *argv[])
{
   
    TString in_root;
   
   
    if ( argc < 2 ){
        cout<<"usage: "<< argv[0] <<" <Input Delphes ROOT File>  \n" ;
        exit(1);
    }
    else{
       
       
        in_root = (TString) argv[1] ;
       
       
        }

  gSystem->Load("libDelphes");

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

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

  //Get pointers to branches used in this analysis
   TClonesArray *branchElectron = treeReader->UseBranch("Electron");
   TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");
  

  //Book histograms
  TH1 *histMassT = new TH1F("M_{T}", "Transverse Mass", 200, 0.0, 150.0);

  //Print number of events
  cout<<"Number of events:"<<numberOfEntries <<endl;
  
  //Event Loop
  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
    {
      //Load selected branches with dada from specified event
      treeReader->ReadEntry(entry);

      //check if there is at least one electron
      if(branchElectron->GetEntries() > 0 && branchMissingET->GetEntries() > 0 ){


      
      Electron *electron = (Electron*)branchElectron->At(0);
      MissingET *met = (MissingET*) branchMissingET->At(0);

      double x = electron->PT;
      double y = electron->Phi;
      double z = met->Phi;
      double xyz = met->MET;


      double fn2 = 2*x*xyz*(1-cos(y-z));

      //final transverse mass variable
      double fn = sqrt(fn2);
 
      //Plot the transverse mass for the decay process
      histMassT->Fill(fn);
    
      //print the invariant mass
      cout << "Tranverse Mass:  "<<fn << endl;
      }
    }

  // Show resulting histograms
    histMassT->Draw();
    gPad->Update();
    gPad->Print();
    
    
    
  return 0;
}