#define trial_cxx
#include "trial.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TTree.h>
#include "TROOT.h"
#include "TRint.h"
#include <TH1F.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <vector>
#include <string>
#include "TAxis.h"
#include <TChain.h>
#include <TChainElement.h>
#include "TMatrixDSymEigen.h"
#include "TMatrixDSym.h"
#include "TVectorD.h"
#include <TLorentzVector.h>


#ifdef __MAKECINT__
#pragma link C++ class std::vector < std::vector<float> >+;
#pragma link C++ class std::vector < std::vector<string> >+;
#pragma link C++ class std::vector < std::vector<int> >+;
#endif

using namespace ROOT;
using namespace std;

//struct isHarder {
//bool operator()(TLorentzVector a, TLorentzVector b){ return a.Pt() > b.Pt(); };
//};

void trial::Loop()
{


    struct isHarder {
  bool operator()(TLorentzVector a, TLorentzVector b){ return a.Pt() > b.Pt(); };
  };


  //////////////Histogram declaration/////////////////////////////////

  TFile*f2  = new TFile("tauVar_H+_130.root","RECREATE");

  TH1F *cosTheta = new TH1F ("cosTheta","cosTheta",100,0,1);
  TH1F *mT = new TH1F ("mT","mT",100,0,200);
  TH1F *TaupT = new TH1F ("TaupT","TaupT",100,0,200);
  TH1F *BpT = new TH1F ("BpT","BpT",100,0,500);
  TH1F *PipT = new TH1F ("pipT","pipT",100,0,50);
  TH1F *TaupTbyBpT = new TH1F ("TaupTbyBpT","TaupTbyBpT",100,0,6);
  TH1F *PipTbyTaupT = new TH1F ("PipTbyTaupT","PipTbyTaupT",100,0,2);

  TH1F * W_had =new TH1F ("Wmass","Wmass",100,0,120);
  TH1F * Top_had= new TH1F ("Topmass","Topmass",100,0,200);




  fChain->SetBranchStatus("*",0);  // disable all branches                                                                                               

  fChain->SetBranchStatus("jet_n",1);
  fChain->SetBranchStatus("jet_pt", 1);
  fChain->SetBranchStatus("jet_eta", 1);
  fChain->SetBranchStatus("jet_phi", 1);
  fChain->SetBranchStatus("jet_E",1);

  fChain->SetBranchStatus("jet_flavor_weight_JetFitterCOMBNN",1);

  fChain->SetBranchStatus("vxp_nTracks",1);
  /*
  float pionptx=0;
  float pionpty=0;
  float pionptz=0;
  float pionE=0;

  float Bptx=0;
  float Bpty=0;
  float Bptz=0;
  float BE=0;

  float mTop=171200;
  float mW=80398;

  float CosTheta=0;
  */

  //  vector<int> GoodPionIndex;
  vector<int> GoodBIndex;

   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
 
     
	
   //============================================Getting good jets and B-jets===============

      
     vector<TLorentzVector>  JetVec;
     double combNN;
     vector<TLorentzVector>  BJetVec;
     GoodBJetIndex.clear();
     
     for( int k = 0 ; k < jet_n ; ++k) {
       
       if ((jet_pt->at(k))>25000 && fabs(jet_eta->at(k))<2.5){
	 TLorentzVector j; 
	 
	 float_t p =jet_pt->at(k);
	 cout<<p<<endl;

	 float_t e =jet_eta->at(k);       
	 cout<<e<<endl;

	 float_t ph=jet_phi->at(k);
	 cout<<ph<<endl;

 	 float_t et=jet_E->at(k);
	 cout<<et<<endl;

	 j.SetPtEtaPhiE(p,e,ph,et);
	 JetVec.push_back(j);
	 
	 GoodJetIndex.push_back(k);
	 combNN=jet_flavor_weight_JetFitterCOMBNN->at(k);
	 
      
	 if (combNN>0.35){
	   TLorentzVector b;
	   b.SetPtEtaPhiE(jet_pt->at(k),jet_eta->at(k),jet_phi->at(k),jet_E->at(k));
	   GoodBJetIndex.push_back(k);
	   BJetVec.push_back(b);
	   GoodBJetIndex.push_back(k);
	 }
	 

     }
	 
}         

              

     TLorentzVector  J1 = JetVec.at(0);

     TLorentzVector J2 = JetVec.at(1);
   

      
     
 
  }      
      f2->cd();
      f2->Write();
      f2->Close();
      

}