#include <TSystem.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TH2.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TLine.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TFitResult.h>
#include <vector>
#include <iostream>
#include <stdlib.h>

Double_t Pi = TMath::Pi();
Double_t ndof = 4; //(four independent sigmas)

///------------------------------ configs to run --------------------------------------
///Path to get/store input/output file
//TString path = "../NeutrinoPU140/";
//TString path = "../SingleMuonPU0/";
TString path = "../SingleMuonPU140/";
//TString path = "../SinglePionPU0/";
//TString path = "../SinglePionPU140/";
//TString path = "../SingleElectronPU0/";

///For single particle events (save time)
bool  mc_singleTrack = true;
bool  isSignal = true;

///Select logic (6 to 5/6, 8 to 6/6 and -1 to deactive the flag- takes all of them)
Int_t use_logic = -1;

///min/maxChi2 (-1 to deactive the flag)
Double_t minChi2 = -1;
Double_t maxChi2 =  3.;

///min/maxPt (-1 to deactive the flag)
Double_t minPt = 3;
Double_t maxPt = -1;

///Restrict to tt27
bool tt27_filter = true;

///tt27 Physic Definition 
Double_t phi_i = (3/8.)*Pi-Pi/16.;//Pi/4.;
Double_t phi_f = (3/8.)*Pi+Pi/16.;//Pi/2.;
Double_t eta_i = 0.;
Double_t eta_f = 0.733333;

///maxEntries (-1 to deactive the flag)
Int_t maxEntries = -1;
///------------------------------------------------------------------------------------


///Functions from fitting single pion (no PU) resolutions                                                      
///q/pT sigma function of q/pT                                                                                 
Double_t r_qbpT(Double_t qbpT){
  const Double_t Const   =  0.001648;
  const Double_t G_Const = -0.001465;
  const Double_t G_Media =  0.000920;
  const Double_t G_Sigma =  0.146079;
  Double_t r = Const + G_Const*TMath::Gaus(qbpT, G_Media, G_Sigma);

  return r;
}

///phi0 sigma funtion of q/pT                                                                                  
Double_t r_phi0(Double_t qbpT){
  const Double_t Const   =  0.000642;
  const Double_t G_Const = -0.000538;
  const Double_t G_Media = -0.001350;
  const Double_t G_Sigma =  0.151268;
  Double_t r = Const + G_Const*TMath::Gaus(qbpT, G_Media, G_Sigma);

  return r;
}

///z0 sigma function of q/pT                                                                                   
Double_t r_z0(Double_t qbpT){
  const Double_t Const   =  0.100000;
  const Double_t G_Const = -0.021369;
  const Double_t G_Media =  0.015085;
  const Double_t G_Sigma =  0.370718;
  Double_t r = Const + G_Const*TMath::Gaus(qbpT, G_Media, G_Sigma);

  return r;
}

///cotTheta sigma function of q/pT                                                                             
Double_t r_cotTheta(Double_t qbpT){
  const Double_t Const   =  0.002985;
  const Double_t G_Const = -0.000892;
  const Double_t G_Media =  0.003922;
  const Double_t G_Sigma =  0.249607;
  Double_t r = Const + G_Const*TMath::Gaus(qbpT, G_Media, G_Sigma);

  return r;
}
///--------------------------------------------------------------------------------------                       


void getNormalizedDeltas(void){
  
  ///File to be analised
  TString file = path+"tracks_from_LTF_logic5b6.root";
  ///----------- list configs of run ----------------------------
  cout<<"Computing deltas in pT slices..."<<endl;
  cout<<"Using logic "<<use_logic<<endl;
  cout<<Form("[minPt,maxPt]: [%.1f, %.1f]",minPt,maxPt)<<endl;
  cout<<Form("[Phi_i,Phi_f]: [%.5f, %.5f]",phi_i,phi_f)<<endl;
  cout<<Form("[Eta_i,Eta_f]: [%.5f, %.5f]",eta_i,eta_f)<<endl;
  cout<<Form("[minChi2, maxChi2]: [%.2f, %.2f]",minChi2,maxChi2)<<endl;
  cout<<"Analysing "<<maxEntries<<" from file "<<file<<endl;
  cout<<endl;
  ///------------------------------------------------------------

  
  Int_t bin_ipt = 20;
  Double_t iptxi = -0.4;
  Double_t iptxf = 0.4;

  ///Stores the parameters delta
  TH2D *d_invPt = new TH2D("d_invPt","#deltaq/p_{T}",bin_ipt,iptxi,iptxf,100,-6,6);
  d_invPt->GetXaxis()->SetTitle("AM q/p_{T}");
  d_invPt->GetYaxis()->SetTitle("#deltaq/p_{T}");
  TH2D *d_invPt_z = new TH2D("d_invPt_z","#deltaq/p_{T}",bin_ipt,-0.05,0.05,100,-6,6);
  d_invPt_z->GetXaxis()->SetTitle("AM q/p_{T}");
  d_invPt_z->GetYaxis()->SetTitle("#deltaq/p_{T}");
  TH2D *d_phi0 = new TH2D("d_phi0","#delta#phi_{0}",bin_ipt,iptxi,iptxf,100,-6,6);
  d_phi0->GetXaxis()->SetTitle("AM q/p_{T}");
  d_phi0->GetYaxis()->SetTitle("#delta#phi_{0}");
  TH2D *d_z0 = new TH2D("d_z0","#deltaz_{0}",bin_ipt,iptxi,iptxf,100,-6,6);
  d_z0->GetXaxis()->SetTitle("AM q/p_{T}");
  d_z0->GetYaxis()->SetTitle("#deltaz_{0}");
  TH2D *d_cotTheta = new TH2D("d_cotTheta","#deltacot#theta ",bin_ipt,iptxi,iptxf,100,-6,6);
  d_cotTheta->GetXaxis()->SetTitle("AM q/p_{T}");
  d_cotTheta->GetYaxis()->SetTitle("#deltacot#theta");

  TH1D *TChi2 = new TH1D("TChi2","Match chi2",500,0,10);
  TH1D *TChi2_ndof = new TH1D("TChi2_ndof","Match chi2/ndof",500,0,10);

  
  vector<bool> *trkParts_signal=0;
  vector<float> *trkParts_pt=0, *trkParts_eta=0, *trkParts_phi=0, *trkParts_vz=0;
  vector<float> *AMTTTracks_pt=0, *AMTTTracks_invPt=0, *AMTTTracks_phi0=0, *AMTTTracks_z0=0, *AMTTTracks_cottheta=0, *AMTTTracks_chi2=0;
  vector<int> *trkParts_charge=0, *AMTTTracks_ndof=0;

  TFile *infile = TFile::Open(file);  
  TTree *intree = (TTree*)infile->Get("ntupler/tree");
  intree->SetBranchAddress("trkParts_charge",&trkParts_charge);
  intree->SetBranchAddress("trkParts_pt",&trkParts_pt);
  intree->SetBranchAddress("trkParts_eta",&trkParts_eta);
  intree->SetBranchAddress("trkParts_phi",&trkParts_phi);
  intree->SetBranchAddress("trkParts_vz",&trkParts_vz);
  intree->SetBranchAddress("trkParts_signal",&trkParts_signal);
  intree->SetBranchAddress("AMTTTracks_pt",&AMTTTracks_pt);
  intree->SetBranchAddress("AMTTTracks_invPt",&AMTTTracks_invPt);
  intree->SetBranchAddress("AMTTTracks_phi0",&AMTTTracks_phi0);
  intree->SetBranchAddress("AMTTTracks_z0",&AMTTTracks_z0);
  intree->SetBranchAddress("AMTTTracks_cottheta",&AMTTTracks_cottheta);
  intree->SetBranchAddress("AMTTTracks_chi2",&AMTTTracks_chi2);
  intree->SetBranchAddress("AMTTTracks_ndof",&AMTTTracks_ndof);
  Long64_t Nentries = intree->GetEntries();
  if(maxEntries != -1 && maxEntries < Nentries) Nentries = maxEntries;
  else if(maxEntries > Nentries) cout<<"maxEntries > Nentries! Using maximum entries: "<<Nentries<<endl;

    
  Double_t lgChi2 = 0;
  for(Long64_t i=0; i<Nentries; i++){
    intree->GetEntry(i);
    if(i%(Nentries/10) == 0) 
      cout<<"Event: "<<i<<endl;
    
    Int_t mc_tracks = trkParts_charge->size();
    Int_t am_tracks = AMTTTracks_invPt->size();

    //if(mc_tracks == 0 || am_tracks == 0) continue;
    //cout<<Form("\n\nEvent(%i): nMCtracks= %i\t\tnAMtracks= %i",i,mc_tracks,am_tracks)<<endl;

    ///Loop over all MC tracks inside Trigger Tower 27
    for(Int_t mc_itrack=0; mc_itrack<mc_tracks; mc_itrack++){
      bool is_Signal = (*trkParts_signal)[mc_itrack];
      if(mc_singleTrack == true && is_Signal != isSignal) continue;
      if(abs((*trkParts_charge)[mc_itrack]) == 0) continue;
      if( (minPt!=-1 && (*trkParts_pt)[mc_itrack]<minPt) || (maxPt!=-1 && (*trkParts_pt)[mc_itrack]>maxPt) )
	continue;
      ///Limit to tower 27
      if(tt27_filter == true){
	if( (*trkParts_phi)[mc_itrack]<phi_i || (*trkParts_phi)[mc_itrack]>phi_f ) continue;
	if( (*trkParts_eta)[mc_itrack]<eta_i || (*trkParts_eta)[mc_itrack]>eta_f ) continue;
      }

      //float trkParts_pti      = (*trkParts_pt)[mc_itrack];
      float trkParts_invPt    = (*trkParts_charge)[mc_itrack]/(*trkParts_pt)[mc_itrack];
      float trkParts_phi0     = (*trkParts_phi)[mc_itrack];
      float trkParts_z0	      = (*trkParts_vz)[mc_itrack];
      float trkParts_cottheta = TMath::SinH((*trkParts_eta)[mc_itrack]);
      if(fabs(trkParts_z0) > 5) continue;

      ///To check AM inefficiency
      if(mc_tracks != 0 && am_tracks == 0){
        TChi2->Fill(-1);
        TChi2_ndof->Fill(-1);
      }      
      
      ///Loop over all AM tracks inside Trigger Tower 27
      Double_t min_dinvPt = 999, min_dphi0 = 999, min_dz0 = 999, min_dcotTheta = 999, min_tChi2 = 1.e15;
      Double_t am_ipt_min_dinvPt = -1, am_ipt_min_dphi0 = -1, am_ipt_min_dz0 = -1, am_ipt_min_dcotTheta = -1;
      for(Int_t am_itrack=0; am_itrack<am_tracks; am_itrack++){
        double track_chi2 = (*AMTTTracks_chi2)[am_itrack]/float((*AMTTTracks_ndof)[am_itrack]);
        if(minChi2 != -1 && track_chi2 < minChi2) continue;
        if(maxChi2 != -1 && track_chi2 > maxChi2) continue;
	if( (minPt != -1 && (*AMTTTracks_pt)[am_itrack]<minPt) || (maxPt != -1 && (*AMTTTracks_pt)[am_itrack]>maxPt) ) continue;
	if( use_logic != -1 && (*AMTTTracks_ndof)[am_itrack] != use_logic ) continue;
	///Limit to tower 27
	if(tt27_filter == true){
	  if( (*AMTTTracks_phi0)[am_itrack]<phi_i || (*AMTTTracks_phi0)[am_itrack]>phi_f ) continue;
	  if( TMath::SinH((*AMTTTracks_cottheta)[am_itrack])<eta_i || TMath::SinH((*AMTTTracks_cottheta)[am_itrack])>eta_f ) continue;
	}

	//float AM_pt         = (*AMTTTracks_pt)[am_itrack];
	float AM_invPt      = (*AMTTTracks_invPt)[am_itrack];
	float AM_phi0	    = (*AMTTTracks_phi0)[am_itrack];
	float AM_z0	    = (*AMTTTracks_z0)[am_itrack];
	float AM_cottheta   = (*AMTTTracks_cottheta)[am_itrack];
	if(fabs(AM_z0) > 5) continue;

        ///Compute parameters delta
	//float dpt       = fabs(AM_pt)-fabs(trkParts_pti);
  	float dinvPt    = AM_invPt-trkParts_invPt;
        float dphi0     = AM_phi0-trkParts_phi0;
        float dz0       = AM_z0-trkParts_z0;	
        float dcotTheta = AM_cottheta-trkParts_cottheta;
		

	///Finds the minimum delta for each parameter
	//if( dpt/fabs(AM_pt)< min_dpt ){
	  //min_dpt = dpt/fabs(AM_pt);
	  //am_ipt_min_dpt = AM_invPt;
	//}
	if( dinvPt < min_dinvPt ){
	  min_dinvPt = dinvPt;
	  am_ipt_min_dinvPt = AM_invPt;
	}
	if( dphi0 < min_dphi0 ){
	  min_dphi0 = dphi0;
	  am_ipt_min_dphi0 = AM_invPt;
	}
	if( dz0 < min_dz0 ){
	  min_dz0 = dz0;
	  am_ipt_min_dz0 = AM_invPt;
	}
	if( dcotTheta < min_dcotTheta ){
	  min_dcotTheta = dcotTheta;
	  am_ipt_min_dcotTheta = AM_invPt;
	}

	///Pseudo chi2
	Double_t tChi2 = pow(dinvPt/r_qbpT(AM_invPt),2) + pow(dphi0/r_phi0(AM_invPt),2) + pow(dz0/r_z0(AM_invPt),2) + pow(dcotTheta/r_cotTheta(AM_invPt),2);
	if(tChi2 < min_tChi2) min_tChi2 = tChi2;
      }///End on AM tracks loop


      ///Normalized deltas
      if(min_tChi2 != 1.e15){
	Double_t norm_dinvPt    = min_dinvPt/r_qbpT(am_ipt_min_dinvPt);
	Double_t norm_dphi0     = min_dphi0/r_phi0(am_ipt_min_dphi0);
	Double_t norm_dz0       = min_dz0/r_z0(am_ipt_min_dz0);
	Double_t norm_dcotTheta = min_dcotTheta/r_cotTheta(am_ipt_min_dcotTheta);
	
	d_invPt->Fill( am_ipt_min_dinvPt, norm_dinvPt );
	d_invPt_z->Fill( am_ipt_min_dinvPt, norm_dinvPt );
	d_phi0->Fill( am_ipt_min_dphi0, norm_dphi0 );
	d_z0->Fill( am_ipt_min_dz0, norm_dz0 );
	d_cotTheta->Fill( am_ipt_min_dcotTheta, norm_dcotTheta );
            
	///Minimum match chi2
	TChi2->Fill(min_tChi2);
	TChi2_ndof->Fill(min_tChi2/ndof);
	if(min_tChi2 > lgChi2) lgChi2 = min_tChi2;
      }

      if(mc_singleTrack == true) break;///Takes only one track (I saw that signal points to more than one!)
    }///End on trkParts loop
     
 }///End on events loop
 cout<<"Largest Chi2= "<<lgChi2<<"\tLargest Chi2/ndof= "<<lgChi2/ndof<<endl;
 
 delete trkParts_signal;
 delete trkParts_pt;
 delete trkParts_eta;
 delete trkParts_phi;
 delete trkParts_vz;
 delete AMTTTracks_pt;
 delete AMTTTracks_invPt;
 delete AMTTTracks_phi0;
 delete AMTTTracks_z0;
 delete AMTTTracks_cottheta;
 delete AMTTTracks_chi2;
 delete trkParts_charge;
 delete AMTTTracks_ndof;
 
  ///Storing the things...
  TFile *f = new TFile(path+"histos_normalized_deltas_allogics_smallphi_smallz0.root","recreate");
  d_invPt->Write();
  d_invPt_z->Write();
  d_phi0->Write();
  d_z0->Write();
  d_cotTheta->Write();
  TChi2->Write();
  TChi2_ndof->Write();
    
  TString paths[5] = {"qbpT","qbpT_z","phi0","z0","cotTheta"};
  TH2D *Ohists[5] = {d_invPt, d_invPt_z, d_phi0, d_z0, d_cotTheta};
  for(Int_t ih=0; ih<5; ih++){
    f->mkdir(paths[ih]);
    f->cd(paths[ih]);        
    cout<<"Saving q/pT slices for "<<paths[ih]<<endl;
    for(Int_t rgs=0; rgs<Ohists[ih]->GetNbinsX(); rgs++){
      Double_t low = Ohists[ih]->GetBinLowEdge(rgs+1);
      Double_t upp = low + Ohists[ih]->GetBinWidth(rgs+1);
      TH1D *pT_slice = Ohists[ih]->ProjectionY(paths[ih]+Form("_ipT_%.3f_%.3f",low,upp),rgs+1,rgs+1);
      if(pT_slice->Integral() == 0) continue;
      ///Storing slice plots
      pT_slice->SetLineColor(9);
      pT_slice->Write();
    }
  }
  
  f->Close();


}///End of program