///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>

#include "TStopwatch.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

//#define nData     100
#define nData     1000
#define nMonteCarlo		10000

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;
}
///================================================================================================



void FastME_ProcPool2n(UInt_t ncores = 2){
   
  ///--------------     Files     -------------------------------------------------------------
  std::vector<std::string> MCs;
#if 0
  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");
#else
  std::string Data_Path("data/NewVBF_QED4_QCD0_SIG_FastME_1.root");  
  MCs.push_back("data/NewVBF_QED4_QCD0_SIG_FastME_2.root");
  MCs.push_back("data/NewVBF_QED4_QCD0_BKG_FastME_2.root");
#endif
  ///------------------------------------------------------------------------------------------
  
  ///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");

    //cout<<"\ntreadEntries: "<<tread.GetEntries(true)<<"\t\trefReaderEntries: "<<refReader.GetEntries(true)<<endl;
    for(int dt=0; dt<nData; dt++){
      if(dt!= 0 && dt%1000 == 0) cout<<":: Remaining Data: "<<nData-dt<<endl;
      refReader.SetEntry(dt); ///Move on loop on Data
      
      for(int mc=0; mc<nMonteCarlo; mc++){
      	tread.SetEntry(mc); ///Move on loop on MC
        ///=================================================================================
        ///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 ) continue;

        //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);
      ///=================================================================================
		
          //cout<<"DataEvent: "<<dt<<"\t\tMcEvent: "<<mc<<"\t\tEventsDistance: "<<event_distance<<endl;
	  //cout<<"(#e,u - dt dt/mc mc): "<<dt_el<<" "<<dt_mu<<"\t"<<mc_el<<" "<<mc_mu<<endl;
	      h->Fill(event_distance);
      }      
    }
    
    delete fData;
    return h;
  };
  
  ///For timming the process
  TStopwatch t;
  t.Start();
  ///--------------------------

  TH1D *hf = 0;
  if (ncores > 0) {
     cout<<"::: Starting Analysis with TProcPool :::"<<endl;
     // Use procpool
     TProcPool workers(ncores);
     hf = (TH1D *) workers.ProcTree(MCs, workItem);
  } else {
     cout<<"::: Starting Analysis in sequential mode :::"<<endl;
     // Sequential
     TH1D *hfx = 0;
     TList *hl = new TList;
     ///This is for MC
     for (int i = 0; i < 2; i++) {
        TFile *fMC = TFile::Open(MCs[0].c_str());
        TTreeReader MCReader("VBF",fMC);
        hfx = (TH1D *) workItem(MCReader);
        if (!hf) {
          hf = hfx;
        } else {
          hl->Add(hfx);
        }
        delete fMC;
      }
      hf->Merge(hl); 
  }
  
  //hf->Print();
  hf->DrawClone();
  
  ///Stoping timming
  t.Stop();
  t.Print();
  
}
///====================================================================================================