#include <iostream>
#include <fstream>
#include <cmath>

#include "../../lut_creators/Calculators/UnifiedMultiplicityCalculator.h"

#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TText.h>
#include <TH1.h>

#include "../../lut_creators/Common/common.h"

using namespace std;

const Float_t PI = 3.141592;

int main();
Float_t delta_R( Float_t eta1, Float_t eta2, Float_t phi1, Float_t phi2 );

int main() {

  ifstream lut_file( "./luts/all_inclusive_lut.data" );
  LUT::UnifiedMultiplicityCalculator calc( lut_file );

  const Int_t MAX_DataWord = 1024;
  const Int_t MAX_Particle = 50;
  const char* file_name = "/home/adevir/double_muon/jpsi_test.root";

  TFile *frdo = new TFile( file_name, "READ" );

  // TTree(s):
  TTree *tt_roi;
  TTree *tt_muctpi;
  TTree *tt_truth;

  tt_truth  = ( TTree* )frdo->Get( "TRUTH/200" );
  tt_roi    = ( TTree* )frdo->Get( "MUONROI/200" );
  tt_muctpi = ( TTree* )frdo->Get( "MUCTPI/200" );

  // Truth Data
  Int_t   Spcl_Num;
  Int_t   Spcl_ID [ MAX_Particle ];
  Float_t Spcl_Phi[ MAX_Particle ];
  Float_t Spcl_Eta[ MAX_Particle ];
  Float_t Spcl_Pt [ MAX_Particle ];
  // Muctpi RDO Data
  Int_t canMulti;
  Int_t nDataWord;
  Int_t dataWord[ MAX_DataWord ];
  Int_t nCandidates;
  //  Double_t etaCand[ MAX_RoI ];
  //  Double_t phiCand[ MAX_RoI ];

  // Assigning the Data
  tt_truth->SetBranchAddress( "Spcl_Num", &Spcl_Num );
  tt_truth->SetBranchAddress( "Spcl_ID", Spcl_ID );
  tt_truth->SetBranchAddress( "Spcl_Phi", Spcl_Phi );
  tt_truth->SetBranchAddress( "Spcl_Eta", Spcl_Eta );
  tt_truth->SetBranchAddress( "Spcl_Pt", Spcl_Pt );
  tt_muctpi->SetBranchAddress( "canMulti", &canMulti );
  tt_muctpi->SetBranchAddress( "nDataWord", &nDataWord );
  tt_muctpi->SetBranchAddress( "dataWord", dataWord );
  tt_roi->SetBranchAddress( "nCandidates", &nCandidates );
  //  tt_roi->SetBranchAddress( "etaCand", etaCand );
  //  tt_roi->SetBranchAddress( "phiCand", phiCand );

  // Making the Canvases
  TCanvas ref_deltaR_window( "ref_deltaR_window", "8500 events analysis", 800, 1100 );
  ref_deltaR_window.cd();
  TH1I Ref_deltaR_dist( "Ref_deltaR_dist", "", 100, 0.0, 1.0 );
  Ref_deltaR_dist.AppendPad();

  TCanvas one_deltaR_window( "one_deltaR_window", "8500 events analysis", 800, 1100 );
  one_deltaR_window.cd();
  TH1I One_deltaR_dist( "One_deltaR_dist", "", 100, 0.0, 1.0 );
  One_deltaR_dist.AppendPad();

  TCanvas two_deltaR_window( "two_deltaR_window", "8500 events analysis", 800, 1100 );
  two_deltaR_window.cd();
  TH1I Two_deltaR_dist( "Two_deltaR_dist", "", 100, 0.0, 1.0 );
  Two_deltaR_dist.AppendPad();

  TCanvas three_deltaR_window( "three_deltaR_window", "8500 events analysis", 800, 1100 );
  three_deltaR_window.cd();
  TH1I Three_deltaR_dist( "Three_deltaR_dist", "", 100, 0.0, 1.0 );
  Three_deltaR_dist.AppendPad();

  // Variables
  Int_t muon_number = 0;
  Int_t total_number = 0;
  Int_t one_muon = 0;
  Int_t double_muon = 0;
  Int_t triple_muon = 0;
  Int_t quatro_muon = 0;
  Int_t five_muon = 0;
  Float_t procent = 0;
  Float_t P1 = 0;
  Float_t P2 = 0;
  Float_t Eta1 = 0;
  Float_t Eta2 = 0;
  Float_t Phi1 = 0;
  Float_t Phi2 = 0;
  Float_t dist = 0;

  for( Int_t i = 0; i < (Int_t)tt_muctpi->GetEntries(); ++i ) {
    tt_muctpi->GetEntry( i );
    tt_roi->GetEntry( i );
    tt_truth->GetEntry( i );
    vector< unsigned int > dataWords;

    for( Int_t j = 0; j < nDataWord; ++j ) dataWords.push_back( dataWord[ j ] );
    muon_number = calc.multiplicity( dataWords );
    total_number += muon_number;
    if ( muon_number >= 1) ++one_muon;
    if ( muon_number >= 2) ++double_muon;
    if ( muon_number >= 3) ++triple_muon;
    if ( muon_number >= 4) ++quatro_muon;
    if ( muon_number >= 5) ++five_muon;

    for( Int_t j = 0; j < Spcl_Num; ++j ) {
      if( ( Spcl_ID[ j ] == 13 ) || ( Spcl_ID[ j ] == -13 ) ) {   //Checking Which is Muon
	if ( Spcl_Pt[ j ] > P1 ) { // Checking which are the Muons with Highest Pt
	  P1 = Spcl_Pt[ j ];
	  Eta1 = Spcl_Eta[ j ];
	  Phi1 = Spcl_Phi[ j ];
	}
	else if ( Spcl_Pt[ j ] > P2 ){
	  P2 = Spcl_Pt[ j ];
	  Eta2 = Spcl_Eta[ j ];
	  Phi2 = Spcl_Phi[ j ];
	}
      }
    }
    if ( P2 != 0 ) {
      dist = delta_R( Eta1, Eta2, Phi1, Phi2 );
      Ref_deltaR_dist.Fill( dist );
      if ( muon_number == 1 ) One_deltaR_dist.Fill( dist );
      else if ( muon_number == 2 ) Two_deltaR_dist.Fill( dist );
      else if ( muon_number == 3 ) Three_deltaR_dist.Fill(dist );
    }

    dataWords.clear();
    P1 = 0;
    P2 = 0;
    muon_number = 0;

  }
  One_deltaR_dist.Divide( Ref_deltaR_dist );
  ref_deltaR_window.SaveAs("./hist/ref_dR.eps");
  one_deltaR_window.SaveAs("./hist/1_dR.eps");
  two_deltaR_window.SaveAs("./hist/2_dR.eps");
  three_deltaR_window.SaveAs("./hist/3_dR.eps");

  Float_t average = ((Float_t)total_number)/((Int_t)tt_muctpi->GetEntries());
  cout << "There were " << (Int_t)tt_muctpi->GetEntries() <<" events in which there were " << total_number << " muons in total."<<endl;
  cout << "The average muons' number per event was "<< average << " muons in an event."<< endl;
  procent = ( (Float_t)one_muon )/( (Int_t)tt_muctpi->GetEntries() )*100;
  cout << "There were " << one_muon << " events in which there were at least 1 muon. They make " << procent << "% from all the " << (Int_t)tt_muctpi->GetEntries() << " events." << endl;
  procent = ( (Float_t)double_muon )/( (Int_t)tt_muctpi->GetEntries() )*100;
  cout << "There were " << double_muon << " events in which there were at least 2 muons. They make " << procent << "% from all the " << (Int_t)tt_muctpi->GetEntries() << " events." << endl;
  procent = ( (Float_t)triple_muon )/( (Int_t)tt_muctpi->GetEntries() )*100;
  cout << "There were " << triple_muon << " events in which there were at least 3 muons. They make " << procent << "% from all the " << (Int_t)tt_muctpi->GetEntries() << " events." << endl;
  procent = ( (Float_t)quatro_muon )/( (Int_t)tt_muctpi->GetEntries() )*100;
  cout << "There were " << quatro_muon << " events in which there were at least 4 muons. They make " << procent << "% from all the " << (Int_t)tt_muctpi->GetEntries() << " events." << endl;
  procent = ( (Float_t)five_muon )/( (Int_t)tt_muctpi->GetEntries() )*100;
  cout << "There were " << five_muon << " events in which there were at least 5 muons. They make " << procent << "% from all the " << (Int_t)tt_muctpi->GetEntries() << " events." << endl;

  return 0;

}

//
// Calculate the distance between two objects
//
Float_t delta_R( Float_t eta1, Float_t eta2, Float_t phi1, Float_t phi2 ) {

  Float_t d_phi = fabs( phi1 - phi2 );
  if ( d_phi > PI )
    d_phi = 2 * PI - d_phi;
  Float_t d_eta = fabs( eta1 - eta2 );

  Float_t distance = sqrt( pow( d_phi, 2 ) + pow( d_eta, 2 ) );

  return distance;

}