///:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
///::::::::::::::::::::::::::::::::::	    FAST MATRIX ELEMENT TO VBF TOPOLOGY		::::::::::::::::::::::::::::::::::
///::::::::::::::::::::::::::::::::::	    Code Author: Miquéias M. de Almeida		::::::::::::::::::::::::::::::::::
///:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
///::															::
///:: This tool takes Data, MC Signal and Background tridimensional matrix containing leptons pT, eta and phi, and      ::
///:: compute the Data-MC event distance in a phase space defined by those leptons properties. In the 4e and 4mu final  ::
///:: state there's two choices, so the tool has two ways to make the comparison between the Data-MC leptons (minimum   ::
///:: distance between leptons and the media of the possible combinations). The tool use two branches to identify those ::
///:: special final states, one for indicate the final state and one for indicate what particle is being handle.	::
///:: The final states must be identify in the following way: 0 = 4e, 1 = 4mu and 2 = 2e2mu				::
///:: 															::
///:: To run the analysis make in Root (compilation mode is the fastest):						::
///:: root -l -b -q FastME.C+												::
///:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::


///To print informations about DR functions running
#define debug		false

#include <iostream>
#include <string>
#include <ctime>
#include <exception>
#include <vector>
#include <TString.h>
#include <TFile.h>
#include <TTree.h>
#include <TMath.h>
#include <TH1.h>
#include <TCanvas.h>

///Headers to TProcPool
#include <TTreeReader.h>
#include <TTreeReaderArray.h>
#include <TProcPool.h>
#include <PoolUtils.h>
#include <TSystem.h>



#define pi		3.14159265358979312
#define pedestal 	-99					///Reset Value to Variables


///Scale Factors to normalize Deltas
#define scale_dPt	50.
#define scale_dEta	5.
#define scale_dPhi	pi

///Final state
#define nObjs		6

using namespace std;



///========== Compute Distance Between Events - Minimum Distance Method ===============
///Takes the combination that gives minimum distance when >1 equals final states
Double_t ComputeDR_MinDist(TTreeReaderArray<Int_t> DataId, TTreeReaderArray<Double_t> DataPt, TTreeReaderArray<Double_t> DataPtRes,
                           TTreeReaderArray<Double_t> DataEta, TTreeReaderArray<Double_t> DataEtaRes, TTreeReaderArray<Int_t> McId,
			   TTreeReaderArray<Double_t> McPt, TTreeReaderArray<Double_t> McPtRes, TTreeReaderArray<Double_t> McEta, 
			   TTreeReaderArray<Double_t> McEtaRes){
  
  Double_t particles_distance, min_particles_distance;
  Double_t dPt=0, dEta=0;
  Double_t sum_dPt2=0, sum_dEta2=0, event_distance=-1;
  
  int min_imc, vmin_imc[6] = {-1,-1,-1,-1,-1,-1};
  for(int idt=0; idt<nObjs; idt++){
    min_imc = -1, particles_distance = -1; min_particles_distance = 1.E15;
    for(int imc=0; imc<nObjs; imc++){
      ///Avoid different Data-MC particles comparison
      if(DataId[idt] != McId[imc]) continue;
      if( debug ) cout<<"DataPos: "<<idt<<"  ID: "<<DataId[idt]<<"  MCPos: "<<imc<<"   ID: "<<McId[imc]<<" >>> ";

      dPt  = (DataPt[idt]-McPt[imc])/(scale_dPt*DataPtRes[idt]);
      dEta = (DataEta[idt]-McEta[imc])/(scale_dEta*DataEtaRes[idt]);
	
      particles_distance = sqrt(dPt*dPt + dEta*dEta);
      if( debug ) cout<<"particles_distance = "<<particles_distance<<endl;
      if(particles_distance < min_particles_distance && imc != vmin_imc[0] && imc != vmin_imc[1]
         && imc != vmin_imc[2] && imc != vmin_imc[3] && imc != vmin_imc[4]){
	min_imc = imc;
	min_particles_distance = particles_distance;
      }
    }
      
    ///Monitor of chosen MCs to avoid object recounting
    vmin_imc[idt] = min_imc;
    if( debug ) cout<<"Chosen->>  MCPos: "<<min_imc<<"   ID: "<<McId[min_imc]<<endl;
      
    dPt  = (DataPt[idt]-McPt[min_imc])/(scale_dPt*DataPtRes[idt]);
    dEta = (DataEta[idt]-McEta[min_imc])/(scale_dEta*DataEtaRes[idt]);
    sum_dPt2  += dPt*dPt;
    sum_dEta2 += dEta*dEta;
  }
    
  event_distance = sqrt(sum_dPt2 + sum_dEta2);
  
  if(event_distance == -1) throw exception();
  else return event_distance;
}
///================================================================================================
///================================================================================================


///========= Compute Discriminant Value ==========================================================
///Based on Distance
Double_t PsbD(Double_t min_dr_sig, Double_t min_dr_bkg){
  Double_t DD = min_dr_bkg/(min_dr_sig + min_dr_bkg);
  return DD;
}
///================================================================================================



int FastME_ProcPool2(UInt_t ncores = 2){
   
  ///--------------     Files     -----------------------------------------
  std::vector<std::string> MCs;
  std::string Data_Path("../VBF_Samples/VBF_QED4_QCD0_BKG_FastME_1.root");  
  MCs.push_back("../VBF_Samples/VBF_QED4_QCD0_SIG_FastME_2.root");
  MCs.push_back("../VBF_Samples/VBF_QED4_QCD0_BKG_FastME_2.root");
  ///------------------------------------------------------------------------------------
  
  ///TProcPool declarations
  auto workItem = [Data_Path](TTreeReader &tread) -> TObject* {
    
    TH1D *h = new TH1D("Distances", "Distances between Data and MC events", 200, 0, 100);
    h->SetDirectory(0);
    
    TTreeReaderArray<Int_t>    DataId(tread, "Particle.id");
    TTreeReaderArray<Double_t> DataPt(tread, "Particle.pt");
    TTreeReaderArray<Double_t> DataPtRes(tread, "Particle.pt_res");
    TTreeReaderArray<Double_t> DataEta(tread, "Particle.eta");
    TTreeReaderArray<Double_t> DataEtaRes(tread, "Particle.eta_res");

    
    TFile *fData = TFile::Open((TString)Data_Path);
    TTreeReader refReader("VBF",fData);
    TTreeReaderArray<Int_t>    McId(refReader, "Particle.id");
    TTreeReaderArray<Double_t> McPt(refReader, "Particle.pt");
    TTreeReaderArray<Double_t> McPtRes(refReader, "Particle.pt_res");
    TTreeReaderArray<Double_t> McEta(refReader, "Particle.eta");
    TTreeReaderArray<Double_t> McEtaRes(refReader, "Particle.eta_res");

    
    while(tread.Next()){
      while(refReader.Next()){
	//refReader.SetEntriesRange(0,10);
	Double_t dr_test = -1.;
	dr_test = ComputeDR_MinDist(DataId, DataPt, DataPtRes, DataEta, DataEtaRes,
	                            McId, McPt, McPtRes, McEta, McEtaRes);
	h->Fill(dr_test);
      }
    }
    
    delete fData;
    return h;
  };
  
  ///For timming the process
  time_t start, stop, delaied;
  double seconds, elapsed_time;
  string unity = "s";
  time(&start);
  ///--------------------------


  ///The process
  TProcPool workers(ncores);
  auto hf = workers.ProcTree(MCs, workItem);
  
  hf->Print();
  hf->DrawClone();  
    
/*    
    ///Getting the discriminant value based on event distance
    if(min_dr_sig == 0 && min_dr_bkg == 0){ 
      min_dr_sig = 1.;
      min_dr_bkg = 1.;
    }
    psb_distance = PsbD(min_dr_sig, min_dr_bkg);
    if(psb_distance < 0 || psb_distance > 1){
      cout<<"[Warning] PSB_Distance: "<<psb_distance<<"    SigD: "<<min_dr_sig<<"  BkgD: "<<min_dr_bkg<<endl;
    }
    minDR_toSig  = min_dr_sig;
    minDR_toBkg  = min_dr_bkg;
    
  }
*/  
  
  ///Stoping timming
  time(&stop);
  cout<<":::::::::: Process Finished ::::::::::::"<<endl;    
  seconds = difftime(stop,start);
  if(seconds < 60) elapsed_time = seconds;
  if(seconds >= 60){
    elapsed_time = seconds/60.;
    unity = "min.";
  }
  if(seconds >= 3600){
    elapsed_time = seconds/3600.;
    unity = "h";
  }
  cout<<":: Time Consumed: "<<elapsed_time<<unity<<endl;
  cout<<"::--------------------------------------\n"<<endl;
  
  
  return 0; //if process well finished
}
///====================================================================================================