#include <iostream>
#include <string>
#include <cmath>

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

#include "common.h"

using namespace std;

const Float_t PI = 3.141592;

void muon_manipulation( const char* rdo_file_name );
Float_t delta_R( Float_t eta1, Float_t eta2, Float_t phi1, Float_t phi2 );

void muon_manipulation( const char* rdo_file_name ) {

  // Set up the file name:
  string file_name = rdo_file_name;
  unsigned int pos = file_name.rfind( '/' );
  file_name.erase( 0, pos + 1 );

  // Open the RDO file:
  TFile *frdo = new TFile( rdo_file_name, "READ" );

  // Get a hold of the tree(s):
  TTree *tt_truth  = ( TTree* )frdo->Get( "TRUTH/200" );

  //
  // Truth information
  //
  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 ];

  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 );

  //
  // Making the Canvas
  //
  TCanvas muon_event_window( "muon_event_window", "8500 events analysis", 800, 1100 );
  muon_event_window.cd();
  TH1I Muon_Event_hist( "Muon_Event_hist", "", 10, 0.0, 10.0 );
  Muon_Event_hist.AppendPad();

  TCanvas muon_deltaR_window( "muon_deltaR_window", "8500 events analysis", 800, 1100 );
  muon_deltaR_window.cd();
  TH1I deltaR_hist( "deltaR_hist", "", 200, 0.0, 1.0 );
  deltaR_hist.AppendPad();

  TCanvas muon_Pt_window( "muon_Pt_window", "8500 events analysis", 800, 1100 );
  muon_Pt_window.cd();
  TH1I Muon_Pt_hist( "Muon_Pt_hist", "", 50, 0.0, 50000.0 );
  Muon_Pt_hist.AppendPad();

  TCanvas fastest_muons_window( "fastest_muons_window", "8500 events analysis", 800, 1100 );
  fastest_muons_window.cd();
  TH1I fastest_muon_hist( "fastest_muon_hist", "", 50, 0.0, 50000.0 );
  fastest_muon_hist.SetLineColor( kRed );
  fastest_muon_hist.AppendPad();
  TH1I faster_muon_hist( "faster_muon_hist", "", 50, 0.0, 50000.0 );
  faster_muon_hist.SetLineColor( kBlue );
  faster_muon_hist.AppendPad();

  TCanvas muon_eta_window( "muon_eta_window", "8500 events analysis", 800, 1100 );
  muon_eta_window.cd();
  TH1I Muon_Eta_hist( "Muon_Eta_hist", "", 200, -5.0, 5.0 );
  Muon_Eta_hist.AppendPad();

  TCanvas muon_phi_window( "muon_phi_window", "8500 events analysis", 800, 1100 );
  muon_phi_window.cd();
  TH1I Muon_Phi_hist( "Muon_Phi_hist", "", 200, -3.2, 3.2 );
  Muon_Phi_hist.AppendPad();


  //
  // Asigning the Truth Muons information
  //
  unsigned int gen_events = 0;
  unsigned int muon_event = 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;

  for( Int_t i = 0; i < (Int_t)tt_truth->GetEntries(); ++i ) {

    tt_truth->GetEntry( i );

    for( Int_t j = 0; j < Spcl_Num; ++j ) {


      //Checking Which is Muon
      if( ( Spcl_ID[ j ] == 13 ) || ( Spcl_ID[ j ] == -13 ) ) {

	++gen_events;
	++muon_event;
	Muon_Eta_hist.Fill( Spcl_Eta[ j ] );
	Muon_Phi_hist.Fill( Spcl_Phi[ j ] );
	Muon_Pt_hist.Fill( Spcl_Pt[ j ] );

	// Checking which are the Muons with biggest Pt
	if ( Spcl_Pt[ j ] > P1 ) {

	  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 ];

	}

      }

    }

    Muon_Event_hist.Fill( muon_event );
    muon_event = 0;

    //
    // Checks if there were 2 muons
    //
    if ( P2!=0 ) {

      const Float_t dist = delta_R( Eta1, Eta2, Phi1, Phi2 );
      deltaR_hist.Fill( dist );
      fastest_muon_hist.Fill( P1 );
      faster_muon_hist.Fill( P2 );

    }

    P1 = 0;
    P2 = 0;
    Eta1 = 0;
    Eta2 = 0;
    Phi1 = 0;
    Phi2 = 0;

  }

  faster_muon_hist.Draw();
  fastest_muon_hist.Draw("same");
  fastest_muons_window.SaveAs("2_Fast_Muons.eps");
  muon_event_window.SaveAs("Muon_per_Event.eps");
  muon_deltaR_window.SaveAs("Muon_deltaR.eps");
  muon_Pt_window.SaveAs("Muon_Pt.eps");
  muon_eta_window.SaveAs("Muon_Eta.eps");
  muon_phi_window.SaveAs("Muon_Phi.eps");

  return;

}


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

}