#include "header.h"
#include "TClonesArray.h"
// Header file for the classes stored in the TTree if any.
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <vector>
//#define newreader_cxx
//#include "newreader_subh.h"
#include <math.h>
#include <iostream>
#include <fstream>

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

//------------------------------------------------------------------------------
//using namespace std;

void subh_test(){


ofstream myfile1;
myfile1.open("subh_test.txt");
//Load shared library
gSystem->Load("~/RESEARCH/EVENT_GENERATORS/MADGRAPH5/MG5_aMC_v2_6_3_2/ExRootAnalysis/libExRootAnalysis.so");
gSystem->Load("libPhysics");
gSystem->Load("libMathCore");
//Define some global constant
double Pi=acos(-1.0);
double MwCntrl=80.385;
double MtCntrl=172.5;

//MULTIFILE loop starts here.
for (int jj = 1; jj <=1; jj++){
//Initialization of various parameters
int Ntt=0,Ntt1=0;
double px=0.0, py=0.0, pz=0.0, AE=0.0;
string coupling;

//Create chain of root trees
TChain chain("LHEF");
TString  rootchain;

if(jj==1)       coupling="autl";
if(jj==2)       coupling="autr";

                 // rootchain.Form("~/EE2TJ_LHE_FILES_WO_CUT/%s_lhe.root",coupling.str());
                  rootchain.Form("~/EE2TJ_LHE_FILES_WO_CUT/autl_lhe.root");
chain.Add(rootchain);

// Create object of class ExRootTreeReader
ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
Long64_t numberOfEntries = treeReader->GetEntries();
// Get pointers to branches used in this analysis
TClonesArray *branchParticle = treeReader->UseBranch("Particle");

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

// array for storing indices of different particles
int ine_pos[1];	        ine_pos[0] = 0;
int inebar_pos[1];	inebar_pos[0] = 0;
int l11p_pos[1];	l11p_pos[0] = 0;
int l12p_pos[1];	l12p_pos[0] = 0;
int l21_pos[1];	        l21_pos[0] = 0;
int l22p_pos[1];	l22p_pos[0] = 0;

if(branchParticle->GetEntries() > 0) {  // If branch particle is not empty
// Loop over all particle in a given event
for(Int_t i = 0; i < branchParticle->GetEntries(); ++i) {
// enlisting the index of generated particles
// i'th particle of one event
TRootLHEFParticle *p = (TRootLHEFParticle*) branchParticle->At(i);

//initial electron
if (p->PID == 11 ) 
	ine_pos[0]=i;
//initial positron
if (p->PID == -11  )
	inebar_pos[0]=i;
}

//Finalstate particles
for(Int_t ii = 2; ii < branchParticle->GetEntries(); ++ii) {
// enlisting the index of generated particles
// i'th particle of one event
	TRootLHEFParticle *p = (TRootLHEFParticle*) branchParticle->At(ii);
//bar{u},bar{c}
if ( p->PID == -2 || p->PID == -4)
	l21_pos[0]=ii;
//b
if ( p->PID ==5)
	l22p_pos[0]=ii;
//mu+
if (p->PID ==-13 )
	l11p_pos[0]=ii;
//vm
if (p->PID ==14 )
	l12p_pos[0]=ii;

//we are looking the distribution for l11 for the pol of V1

}

//Defining TLorentz vector
TRootLHEFParticle *ine = (TRootLHEFParticle*) branchParticle->At(ine_pos[0]);
TRootLHEFParticle *inebar = (TRootLHEFParticle*) branchParticle->At(inebar_pos[0]);

TRootLHEFParticle *l21 = (TRootLHEFParticle*) branchParticle->At(l21_pos[0]);
TRootLHEFParticle *l22p = (TRootLHEFParticle*) branchParticle->At(l22p_pos[0]);
TRootLHEFParticle *l11p = (TRootLHEFParticle*) branchParticle->At(l11p_pos[0]);
TRootLHEFParticle *l12p = (TRootLHEFParticle*) branchParticle->At(l12p_pos[0]);

TLorentzVector inev, inebarv, jet0, omupv0, obv0, oposv0, omunv0, testv, recotv,triletonv,ometvR1,ometvR2,qnup,qnum;

inev.SetPxPyPzE(ine->Px,ine->Py,ine->Pz,ine->E);
inebarv.SetPxPyPzE(inebar->Px,inebar->Py,inebar->Pz,inebar->E);

jet0.SetPxPyPzE(l21->Px,l21->Py,l21->Pz,l21->E);
obv0.SetPxPyPzE(l22p->Px,l22p->Py,l22p->Pz,l22p->E);
omupv0.SetPxPyPzE(l11p->Px,l11p->Py,l11p->Pz,l11p->E);
//l12pv.SetPxPyPzE(l12p->Px,l12p->Py,l12p->Pz,l12p->E);

/***********************************************************************************************/

   } // Loop over all events of one file ends.

} // If branch particle is not empty
    }  // MULTIFILE loop ends here.
myfile1.close();

}


#ifndef __CINT__
int main(){
   pxpypzAE_lheroot_EE2TJ();
    return 0;
};
#endif