///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> McEta){
///================================================================================================
  ///Checks the final states
  int dt_el= 0, dt_mu= 0, mc_el= 0, mc_mu=0;
  for(int u=0; u<nObjs; u++){
    if(abs(DataId[u]) == 11) dt_el++;
    if(abs(DataId[u]) == 13) dt_mu++;
    if(abs(McId[u])   == 11) mc_el++;
    if(abs(McId[u])   == 13) mc_mu++;
  }
  
  ///To distiguish between right and wrong framework
  if( dt_el != mc_el || dt_mu != mc_mu ) return -1;
  
  ///If Data and MC final state are equals compute the distance between events
  if( debug ) cout<<"\n\n(#e,#u)>>>> Data ("<<dt_el<<", "<<dt_mu<<")\t\tMc ("<<mc_el<<", "<<mc_mu<<")"<<endl;
  bool mc_approved;
  Double_t particles_distance, min_particles_distance;
  Double_t dPt=0, dEta=0, dPhi=0;
  Double_t sum_dPt2=0, sum_dEta2=0, sum_dPhi2=0, event_distance=-1;
  
  int min_imc;
  std::vector<int> vmin_imc;
  vmin_imc.clear();
  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++){
      mc_approved = true;
      if(int(vmin_imc.size()) > 0)
	for(int g=0; g<int(vmin_imc.size()); g++)
	  if(imc == vmin_imc[g]) mc_approved = false;
	
      if(mc_approved == true){
	///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){
	  min_imc = imc;
	  min_particles_distance = particles_distance;
	}
      }
    }
      
    ///Monitor of chosen MCs to avoid object recounting
    vmin_imc.push_back(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);
///================================================================================================  
  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/StructToTProcPool/NewVBF_QED4_QCD0_SIG_FastME_1.root");  
  MCs.push_back("../VBF_Samples/StructToTProcPool/NewVBF_QED4_QCD0_SIG_FastME_2.root");
  MCs.push_back("../VBF_Samples/StructToTProcPool/NewVBF_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, 20);
    h->SetDirectory(0);

    ///This is for MC
    TTreeReaderArray<Int_t>    McId(tread, "Particle.id");
    TTreeReaderArray<Double_t> McPt(tread, "Particle.pt");
    //TTreeReaderArray<Double_t> McPtRes(tread, "Particle.pt_res");
    TTreeReaderArray<Double_t> McEta(tread, "Particle.eta");
    //TTreeReaderArray<Double_t> McEtaRes(tread, "Particle.eta_res");
    
    ///This is for Data
    TFile *fData = TFile::Open((TString)Data_Path);
    TTreeReader refReader("VBF",fData);
    TTreeReaderArray<Int_t>    DataId(refReader, "Particle.id");
    TTreeReaderArray<Double_t> DataPt(refReader, "Particle.pt");
    TTreeReaderArray<Double_t> DataPtRes(refReader, "Particle.pt_res");
    TTreeReaderArray<Double_t> DataEta(refReader, "Particle.eta");
    TTreeReaderArray<Double_t> DataEtaRes(refReader, "Particle.eta_res");

    Double_t event_distance = -1;
    int iData = 0, nMc = 0;
    cout<<"treadEntries: "<<tread.GetEntries(true)<<"\t\trefReaderEntries: "<<refReader.GetEntries(true)<<endl;
    while(tread.Next()){
      iData++;
      while(refReader.Next()){
	nMc++;
  ///=================================================================================
  ///Checks the final states
  int dt_el= 0, dt_mu= 0, mc_el= 0, mc_mu=0;
  for(int u=0; u<nObjs; u++){
    if(abs(DataId[u]) == 11) dt_el++;
    if(abs(DataId[u]) == 13) dt_mu++;
    if(abs(McId[u])   == 11) mc_el++;
    if(abs(McId[u])   == 13) mc_mu++;
  }

  if(iData < 3 && nMc < 2){
      cout<<"Data("<<iData<<") >> Pt: "<<DataPt[0]<<"\t\tEta: "<<DataEta[0]<<endl;
      cout<<"Mc("<<nMc<<") >> Pt: "<<McPt[0]<<"\t\tEta: "<<McEta[0]<<endl;
  }

  bool mc_approved;
  Double_t particles_distance, min_particles_distance;
  Double_t dPt=0, dEta=0, dPhi=0;
  Double_t sum_dPt2=0, sum_dEta2=0, sum_dPhi2=0;// event_distance=-1;
    
  int min_imc;
  std::vector<int> vmin_imc;
  vmin_imc.clear();
  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++){
      mc_approved = true;
      if(int(vmin_imc.size()) > 0)
	for(int g=0; g<int(vmin_imc.size()); g++)
	  if(imc == vmin_imc[g]) mc_approved = false;
	
      if(mc_approved == true){
	///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){
	  min_imc = imc;
	  min_particles_distance = particles_distance;
	}
      }
    }
      
    ///Monitor of chosen MCs to avoid object recounting
    vmin_imc.push_back(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 < min_dr_sig) min_dr_sig = event_distance;//h->Fill(event_distance);
        //h->Fill( ComputeDR_MinDist(DataId, DataPt, DataPtRes, DataEta, DataEtaRes, McId, McPt, McEta) );
	///To distiguish between right and wrong framework
	if( dt_el == mc_el || dt_mu == mc_mu )
	  h->Fill(event_distance);
      }
      //h->Fill(event_distance);
      //vminDist.push_back(min_dr_sig);
      
    }
    
    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();      
  
  ///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
}
///====================================================================================================