#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/Delphes-3.4.1/external/ExRootAnalysis/ExRootTreeReader.h"
#include "/home/sayan/Delphes-3.4.1/external/ExRootAnalysis/ExRootTreeWriter.h"
#include "/home/sayan/Delphes-3.4.1/external/ExRootAnalysis/ExRootTreeBranch.h"
#include "/home/sayan/Delphes-3.4.1/external/ExRootAnalysis/ExRootResult.h"
#include "/home/sayan/Delphes-3.4.1/external/ExRootAnalysis/ExRootUtilities.h"


#include "/home/sayan/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 *branchJet = treeReader->UseBranch("Jet");

  //Book histograms
  TH1 *histMassT = new TH1F("Mass_{T}", "transverse Mass", 350, 0.0, 500.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);
 
      double mj2[2],ejt[2],x[2],y[2];

      //check if there is at least one jet
      if(branchJet->GetEntries() > 0){
         for(int i = 0; i <2 ; ++i)
	{
	 Jet *jet = (Jet*)branchJet->At(i);
         x[i] = jet->PT;
         y[i] = jet->Phi;
         TLorentzVector lorentz_jet = jet->P4();
         mj2[i] = ((lorentz_jet.E()*lorentz_jet.E()) - ((lorentz_jet.Px()*lorentz_jet.Px()) + (lorentz_jet.Py()*lorentz_jet.Py()) + (lorentz_jet.Pz()*lorentz_jet.Pz())));
         ejt[i] = sqrt(mj2[i] + x[i]*x[i]) ;
        }

        float fn2 = mj2[0] + mj2[1] + (2*ejt[0]*ejt[1]) - (2*x[0]*x[1]*cos(y[0]-y[1]));
        float 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();
}