/**
This Macro
1. Requires 2 extra jets in a given region of pT and eta

Need to specify
0. See Declare constants
1. Pass the name of the rootpla you want to process without the suffix .root (cut_data_OHLT_VBF_2jets.cc'('\"nameRootpla\"')')
*/
/////
//   Prepare Root 
/////
#include <fstream>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <algorithm>
#include "TTree.h"
#include "TTreePlayer.h"
#include "TMath.h"
#include <utility>
#include "TLorentzVector.h"
using namespace std;
using namespace TMath;
/////
//   Declare constants
/////
///Path for rootpla 
const string path_MC = "/home/francesco_romeo_86/ZPrime_TauTau/Trigger/Macro/Efficiency/BugRoot/";
//Requires 2 extra jets in a given region of pT and eta
const float jet_acc = 3; //Acceptance in jet eta
const float Jet_Th  = 30; //Threshold on jet pt 
//Should remain the same (It is needed for the macro)
const string suffisso   = "_cut_VBF_2jets";  
const string dotroot    = ".root";
/////
//   Main function
/////
void cut_data_OHLT_VBF_2jets(string rootpla_mc) {
 string rootpla_MC =  rootpla_mc;
 /////
 //1. Filter events where at least 2 hadronic taus in a given region of pT and eta are generated 
 /////
 //Call old File 
 string oldfil  = rootpla_MC+dotroot;
 TFile *oldfile = new TFile(oldfil.c_str());
 //Define new File
 string newfil  = rootpla_MC+suffisso+dotroot;
 TFile *newfile = new TFile(newfil.c_str(),"recreate");
 //Call old tree
 TTree *oldtree = (TTree*)oldfile->Get("HltTree");
 Int_t nentries = (Int_t)oldtree->GetEntries();
 cout<<nentries<<endl;
 //Variable for the cut
 const Int_t kMaxTrack = 5000;
 Int_t NTauPart;
 oldtree->SetBranchAddress("NTauPart", &NTauPart); 
 Float_t MCTauPt[kMaxTrack]; 
 oldtree->SetBranchAddress("MCTauPt", &MCTauPt);
 Float_t MCTauEta[kMaxTrack];
 oldtree->SetBranchAddress("MCTauEta", &MCTauEta); 
 Float_t MCTauPhi[kMaxTrack];
 oldtree->SetBranchAddress("MCTauPhi", &MCTauPhi);
 //For extra jet
 Int_t NrecoJetGen;
 oldtree->SetBranchAddress("NrecoJetGen", &NrecoJetGen);
 Float_t recoJetGenPt[kMaxTrack];
 oldtree->SetBranchAddress("recoJetGenPt", &recoJetGenPt);
 Float_t recoJetGenEta[kMaxTrack];
 oldtree->SetBranchAddress("recoJetGenEta", &recoJetGenEta);
 Float_t recoJetGenPhi[kMaxTrack];
 oldtree->SetBranchAddress("recoJetGenPhi", &recoJetGenPhi);
 Float_t recoJetGenE[kMaxTrack]; 
 oldtree->SetBranchAddress("recoJetGenE", &recoJetGenE);

 //Copy interested part of the tree from all the events
 TLorentzVector good_Jet0 = TLorentzVector(0,0,0,0);
 TLorentzVector good_Jet1 = TLorentzVector(0,0,0,0);
 TTree *newtree = oldtree->CloneTree(0);
 for (Int_t i=0; i<nentries; i++){
  oldtree->GetEntry(i);
  //Require extra jets
  bool goodJet = false;
  //Loop on gen Jet to find non tau Jet
  float dEta = 0; float dPhi = 0; float dR = 0;
  const int NAllJets = NrecoJetGen;
  int* goodJet_arr = new int[NAllJets];//(The size is for the best case when all recoJet are good)
  int goodJet_num = 0;
  for(int g=0; g<NrecoJetGen; g++){
   bool goodJet_num_bool = false;
   //Jet must be in a given region of pT and eta
   if(recoJetGenPt[g]>=Jet_Th && fabs(recoJetGenEta[g])<=jet_acc){
    //Verify the current Jet does not match any gen taus
    for(int h=0; h<NTauPart; h++){
     dEta = recoJetGenEta[g]-MCTauEta[h];
     dPhi = recoJetGenPhi[g]-MCTauPhi[h];
     dR = sqrt(pow(dEta,2)+pow(dPhi,2));
     if(dR>0.3){goodJet_num_bool = true;}
     else{goodJet_num_bool = false; break;} 
    }//End verify the current Jet does not match any gen taus
   }//End jet must be in a given region of pT and eta  
   if(goodJet_num_bool){goodJet_arr[goodJet_num]=g; goodJet_num++;}
  }//End loop on gen Jet to find non tau Jet
  //Now the combinatoric to find good non tau Jet
  if(goodJet_num>=2){
   for(int g=0; g<(goodJet_num-1); g++){  
    for(int h=g+1; h<goodJet_num; h++){
     if(fabs(recoJetGenEta[goodJet_arr[g]]-recoJetGenEta[goodJet_arr[h]])>4){
      good_Jet0.SetPtEtaPhiE(recoJetGenPt[goodJet_arr[g]],recoJetGenEta[goodJet_arr[g]],recoJetGenPhi[goodJet_arr[g]],recoJetGenE[goodJet_arr[g]]);
      good_Jet1.SetPtEtaPhiE(recoJetGenPt[goodJet_arr[h]],recoJetGenEta[goodJet_arr[h]],recoJetGenPhi[goodJet_arr[h]],recoJetGenE[goodJet_arr[h]]);
      double invMass = pow(good_Jet0.M(),2)+pow(good_Jet1.M(),2)+(2*good_Jet0.P()*good_Jet1.P()*(good_Jet0.Vect().Unit()*good_Jet1.Vect().Unit()))-2*good_Jet0.E()*good_Jet1.E();
      cout<<sqrt(-invMass)<<endl;
      if(sqrt(-invMass)>400) {goodJet = true; break;}
     }
    }
    if(goodJet) break;
   }   
  }//End combinatoric to find good non tau Jet
  //Fill new tree
  if(goodJet) newtree->Fill();
 }//End loop on all events
 cout<<newtree->GetEntries()<<endl;
 newtree->AutoSave();
 newfile->Write();
 newfile->Close();
 delete oldfile;
 delete newfile;
}