//#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 ) {
void muon_manipulation( ) {

  // 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" );
  TFile *frdo = new TFile("jpsi_test.root", "READ" );

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

  //
  // Truth information
  //
  Int_t   Spcl_Num;
  const Int_t MAX_Particle = 1000;
  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=new TCanvas( "muon_event_window", "8500 events analysis", 800, 1100 );
  muon_event_window->cd();
  TH1I *Muon_Event_hist=new TH1I( "Muon_Event_hist", "", 10, 0.0, 10.0 );
  Muon_Event_hist->Draw();

  TCanvas *muon_deltaR_window=new TCanvas( "muon_deltaR_window", "8500 events analysis", 800, 1100 );
  muon_deltaR_window->cd();
  TH1I *deltaR_hist=new TH1I( "deltaR_hist", "", 200, 0.0, 1.0 );
  deltaR_hist->Draw();

  TCanvas *muon_Pt_window=new TCanvas( "muon_Pt_window", "8500 events analysis", 800, 1100 );
  muon_Pt_window->cd();
  TH1I *Muon_Pt_hist=new TH1I( "Muon_Pt_hist", "", 50, 0.0, 50000.0 );
  Muon_Pt_hist->Draw();

  TCanvas *fastest_muons_window=new TCanvas( "fastest_muons_window", "8500 events analysis", 800, 1100 );
  fastest_muons_window->cd();
  TH1I *fastest_muon_hist=new TH1I( "fastest_muon_hist", "", 50, 0.0, 50000.0 );
  fastest_muon_hist->SetLineColor( kRed );
  TH1I *faster_muon_hist=new TH1I( "faster_muon_hist", "", 50, 0.0, 50000.0 );
  faster_muon_hist->SetLineColor( kBlue );
  faster_muon_hist->Draw();
  fastest_muon_hist->Draw("same");

  TCanvas *muon_eta_window=new TCanvas( "muon_eta_window", "8500 events analysis", 800, 1100 );
  muon_eta_window->cd();
  TH1I *Muon_Eta_hist=new TH1I( "Muon_Eta_hist", "", 200, -5.0, 5.0 );
  Muon_Eta_hist->Draw();

  TCanvas *muon_phi_window=new TCanvas( "muon_phi_window", "8500 events analysis", 800, 1100 );
  muon_phi_window->cd();
  TH1I *Muon_Phi_hist=new TH1I( "Muon_Phi_hist", "", 200, -3.2, 3.2 );
  Muon_Phi_hist->Draw();


  //
  // 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 = TMath::Abs( phi1 - phi2 );
  if ( d_phi > TMath::Pi() )
    d_phi = 2 * TMath::Pi() - d_phi;
  Float_t d_eta = TMath::Abs( eta1 - eta2 );

  //Float_t distance = sqrt( pow( d_phi, 2 ) + pow( d_eta, 2 ) );
  Float_t distance = TMath::Sqrt(d_phi*d_phi + d_eta*d_eta );

  return distance;

}