#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 "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/internal/BasicRandom.hh"
//#include "fastjet/tools/MassDropTagger.hh"
//#include "fastjet/tools/Filter.hh"
//#include "fastjet/contrib/Njettiness.hh"
//#include "fastjet/contrib/Nsubjettiness.hh"
//#include "fastjet/contrib/NjettinessPlugin.hh"
//#include  "fastjet/contrib/SoftDrop.hh"


#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;
using namespace fastjet;
//using namespace fastjet::contrib;


int main(int argc, char *argv[]){
    char commandcheck[50], filecheck[50];
    
    TString in_root, out_dat, n_events;
    int events_to_analyze(10);
    
    
    if ( argc < 3 ){
		cout<<"usage: "<< argv[0] <<" <Input Delphes ROOT File> <Output File>  \n" ;
		exit(1);
	}
	else{
        
        
        in_root = (TString) argv[1] ;
        out_dat = (TString) argv[2] ;
        
        
        }
    
    ofstream myfile;
    myfile.open (out_dat);
    

  // Jet Parameters
    
    double R= 0.8;
    RecombinationScheme recomb_scheme=E_scheme;
    //RecombinationScheme recomb_scheme= WTA_pt_scheme ;
    Strategy strategy = Best;
    JetDefinition jet_def(antikt_algorithm, R, recomb_scheme, strategy);
    double ptmin_jet(10.0);


  // Initialize run counters
    
    int event_counter=0;



  // Reading in Delphess tree from input file
    gSystem->Load("/home/sayan/Delphes-3.4.1/libDelphes");
    TChain chain("Delphes");
    chain.Add(in_root);


  // Create class ExRootTreeReader Object
    ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
    Long64_t numberOfEntries = treeReader->GetEntries();
    
    if(argc < 4) events_to_analyze = numberOfEntries;
    
    cout << "Total number of events in file :" << numberOfEntries << endl;
    cout << "Events to be analyzed : " << events_to_analyze << endl;

  // Get pointers to branches used in this analysis

    TClonesArray *branchTower = treeReader->UseBranch("Tower");
    TClonesArray *branchTrack = treeReader->UseBranch("Track");
    TClonesArray *branchMuon = treeReader->UseBranch("Muon");
    TClonesArray *branchMET = treeReader->UseBranch("MissingET");
    
    TObject *object;
    
    myfile<<"index, jL_e, jL_px, jL_py, jL_pz, jL_pt, jL_m"<<endl;

  //event loop
   for(Long64_t entry = 0; entry < events_to_analyze; ++entry){
        
        if (!(entry % 1000)) cout << "At event "<<entry<<endl;
        event_counter++;
        if (event_counter > events_to_analyze) break;
        
       //if (entry == -1 || entry == -5933 || entry == -6705 || entry == -8668 || entry == -9146 || entry == -9492 || entry == -6925) continue;
        
        
        
        
        // Load selected branches with data from specified event
        treeReader->ReadEntry(entry);
        
        // Jets from Towers
        vector<PseudoJet > towerjets;

        //tower loop
        for(Int_t i = 0; i < branchTower->GetEntriesFast(); ++i){
            Tower *delphes_tower = (Tower*) branchTower->At(i);
            
            TLorentzVector lorentz_jet = delphes_tower->P4();
            PseudoJet towerjet(lorentz_jet.Px(), lorentz_jet.Py(), lorentz_jet.Pz(),lorentz_jet.E());
            
            towerjets.push_back(towerjet);
            
        }
        
        towerjets = sorted_by_pt(towerjets);
 
       
        // Form jets and define filters
        vector<PseudoJet > input_particles = towerjets;
        ClusterSequence clust_seq(input_particles, jet_def);
        vector<PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(ptmin_jet));


    /*  // print the jets
  cout <<   "        pt y phi" << entry << endl;
  for (unsigned i = 0; i < jets.size(); i++) {
    cout << "jet " << i << ": "<< jets[i].pt() << " " 
                   << jets[i].rap() << " " << jets[i].phi() << endl;
    vector<PseudoJet> constituents = jets[i].constituents();
    for (unsigned j = 0; j < constituents.size(); j++) {
      cout << "    constituent " << j << "'s pt: " << constituents[j].pt()
           << endl;
    }
  } */


        // Write to file
        myfile <<entry<<", "<< jets[0].e()<< ", "<< jets[0].px()<<", "<< jets[0].py()<< ", "<<jets[0].pz()<< ", "<<jets[0].pt()<< ", "<<jets[0].m()<< ", ";

      }

  myfile.close();
    return 0;
}