//  
//  lumiExtract_pA.C
//  macro to plot the integrated luminosity as a function of time for the proton-lead runs
//
//  Created by Christian Elsasser on 31.05.13.
//  University of Zurich, elsasser@cern.ch
//  Copyright only for commercial use
//

//  General include
#include <iostream>
#include <vector>
#include <map>


//  ROOT include
#include "incBasicROOT.h"
#include "incDrawROOT.h"
#include "incGraphROOT.h"
#include "incIOROOT.h"
#include "incLinAlgROOT.h"



#include "basicROOTTools.h"
#include "basicRooFitTools.h"

#include "lhcbstyle.C"








int lumiExtract(int argc, char* argv[]);

int main(int argc, char* argv[]){
  
  TStyle* style = lhcbStyle();
  
  int argc_copy = argc;
  
  TApplication theApp("Analysis", &argc, argv);
  argc = theApp.Argc();
  argv = theApp.Argv();
  
  
  
  int exit = lumiExtract(argc,argv);
  
  Printf("End");
  theApp.Run(kTRUE); 
  
  return exit;
}


//  Executable method
int lumiExtract(int argc, char* argv[]){
  
  TMultiGraph *mg_lumi = new TMultiGraph("mg_lumi","mg_lumi");
  
  TString dn_lumi("$HOME/cmtuser/Vetra_v13r2/ST/STAging/data/lumi");
  
  std::vector<TString> v_fn_lumi;
  v_fn_lumi.push_back(TString(Form("%s/lumi2013.root",dn_lumi.Data())));
  
  double uT2010 = 1262304000.0;
  
  std::vector<double> v_timeCorrFac;
  v_timeCorrFac.push_back(1.0);
  
  std::vector<double> v_timeCorrScale;
  v_timeCorrScale.push_back(0.0);
  
  
  std::vector<int> v_color;
  v_color.push_back(kRed);
  
  std::vector<double> v_cms;
  v_cms.push_back(5.02);
  
  std::vector<double> v_lumiScale;
  v_lumiScale.push_back(0.001);
  
  
  
  TFile* f_lumiout = new TFile(Form("%s/lumi_pA.root",dn_lumi.Data()),
                                     "RECREATE");
  
  std::vector<float>v_time_tot;
  std::vector<float>v_peak_tot;
  std::vector<float>v_stop_tot;
  std::vector<float>v_beam_tot;
  
  TLegend *leg_lumi = new TLegend(0.23,0.7,0.42,0.9);
  leg_lumi->SetFillColor(0);
  
  for (int i=0; i<v_fn_lumi.size(); i++) {
    TFile* f_lumi = TFile::Open(v_fn_lumi[i].Data());
    if (!f_lumi) {
      return 1;
    }
    TVectorF v_time = *(TVectorF*)f_lumi->Get("TimeofFill Vector");
    TVectorF v_peak = *(TVectorF*)f_lumi->Get("PeakLumi Vector");
    TVectorF v_deli = *(TVectorF*)f_lumi->Get("IntLumiDelivered Vector");
    TVectorF v_reco = *(TVectorF*)f_lumi->Get("IntLumiRecorded Vector");
    
    v_time *= v_timeCorrFac[i];
    v_time += v_timeCorrScale[i];
    
    v_time.Print();
    v_peak.Print();
    v_deli.Print();
    v_reco.Print();
    
    TVectorF v_deli_diff(v_deli.GetNoElements());
    for (int j=0; j<v_deli_diff.GetNoElements(); j++) {
      if (j!=0) {
        v_deli_diff(j) = v_deli(j) - v_deli(j-1);
      }else{
        v_deli_diff(j) = v_deli(j);
      }
    }

    
    TVectorF v_beam(v_peak.GetNoElements());
    v_beam += v_cms[i];
    v_beam.Print();
    
    TVectorF v_stop(v_deli_diff);
    v_stop *= 1e+6;
    
    TVectorF v_deli_plot(v_deli);
    v_deli_plot *= v_lumiScale[i];
    
    TVectorF v_reco_plot(v_reco);
    v_reco_plot *= v_lumiScale[i];
    
    for (int j=0; j<v_stop.GetNoElements(); j++) {
      v_stop(j) *= 1.0/v_peak(j);
      if(v_stop(j)<0.0){
        v_stop(j) = 0.0;
      }else if(TMath::IsNaN(v_stop(j)) || v_stop(j)==TMath::Infinity()){
        v_stop(j) = 0.0;
      }
    }
    
    int index_Ap = 0;
    
    
    
    v_stop += v_time;
    
    TVectorF v_time_plot(v_time);
    v_time_plot -= i*365*86400.0;
    
    for (int j=0; j<v_time_plot.GetNoElements(); j++) {
      if (v_time_plot(j)>(31+1+365*3)*86400.+uT2010) {
        index_Ap = j;
        break;
      }
    }
    
    TVectorF v_time_plot_pA = v_time_plot.GetSub(0,index_Ap);
    TVectorF v_time_plot_Ap = v_time_plot.GetSub(index_Ap,v_time_plot.GetNoElements()-1);
    
    TVectorF v_deli_plot_pA = v_deli_plot.GetSub(0,index_Ap);
    TVectorF v_deli_plot_Ap = v_deli_plot.GetSub(index_Ap,v_time_plot.GetNoElements()-1);
    v_deli_plot_Ap -= v_deli_plot(index_Ap);
    
    TVectorF v_reco_plot_pA = v_reco_plot.GetSub(0,index_Ap);
    TVectorF v_reco_plot_Ap = v_reco_plot.GetSub(index_Ap,v_time_plot.GetNoElements()-1);
    v_reco_plot_Ap -= v_reco_plot(index_Ap);
    
    v_time_plot_pA.Print();
    v_deli_plot_pA.Print();
    v_reco_plot_pA.Print();
    
    v_time_plot_Ap.Print();
    v_deli_plot_Ap.Print();
    v_reco_plot_Ap.Print();
    
    TGraph* g_deli_pA = new TGraph(v_time_plot_pA,v_deli_plot_pA);
    g_deli_pA->SetLineColor(kBlue);
    g_deli_pA->SetLineStyle(7);
    mg_lumi->Add(g_deli_pA,"l");
    
    TGraph* g_reco_pA = new TGraph(v_time_plot_pA,v_reco_plot_pA);
    g_reco_pA->SetLineColor(kBlue);
    g_reco_pA->SetLineStyle(0);
    mg_lumi->Add(g_reco_pA,"l");
    
    TGraph* g_deli_Ap = new TGraph(v_time_plot_Ap,v_deli_plot_Ap);
    g_deli_Ap->SetLineColor(kRed);
    g_deli_Ap->SetLineStyle(7);
    mg_lumi->Add(g_deli_Ap,"l");
    
    TGraph* g_reco_Ap = new TGraph(v_time_plot_Ap,v_reco_plot_Ap);
    g_reco_Ap->SetLineColor(kRed);
    g_reco_Ap->SetLineStyle(0);
    mg_lumi->Add(g_reco_Ap,"l");
    
    
    leg_lumi->AddEntry(g_reco_pA,"#it{p}Pb forward","l");
    leg_lumi->AddEntry(g_reco_Ap,"#it{p}Pb  backward","l");
    
    
    v_time_tot = addROOTVector(v_time_tot,v_time);
    v_peak_tot = addROOTVector(v_peak_tot,v_peak);
    v_stop_tot = addROOTVector(v_stop_tot,v_stop);
    v_beam_tot = addROOTVector(v_beam_tot,v_beam);
    
  }
  
  TCanvas* c1 = new TCanvas("c1","Int luminosity",800,600);
  
  TGraph* gd_2 = new TGraph();
  gd_2->SetLineStyle(7);
  leg_lumi->AddEntry(gd_2,"Delivered","l");
  TGraph* gd_1 = new TGraph();
  leg_lumi->AddEntry(gd_1,"Recorded","l");
  
  
  
  mg_lumi->Draw("A");
  mg_lumi->GetXaxis()->SetTitle("date");
  mg_lumi->GetXaxis()->SetTimeDisplay(1);
  mg_lumi->GetXaxis()->SetTimeFormat("%d/%b");
  mg_lumi->GetXaxis()->SetLimits(uT2010+3*365*86400.0+15*86400.,uT2010+3*365*86400.0+45*86400.);
  
  mg_lumi->GetYaxis()->SetTitle("Luminosity [nb^{-1}]");
  mg_lumi->GetYaxis()->SetTitleOffset(1.5);
  mg_lumi->GetHistogram()->SetMaximum(2.0);
  
  leg_lumi->Draw();
  
  mg_lumi->GetHistogram()->Draw("axissame");
  
  TVectorF vr_time_tot = getROOTVector(v_time_tot);
  TVectorF vr_peak_tot = getROOTVector(v_peak_tot);
  TVectorF vr_stop_tot = getROOTVector(v_stop_tot);
  TVectorF vr_beam_tot = getROOTVector(v_beam_tot);
  
  TMultiGraph *mg_peak = new TMultiGraph("mg_peak","mg_peak");
  
  
  TVectorF v_time_tot_plot(4*vr_time_tot.GetNoElements());
  TVectorF v_peak_tot_plot(4*vr_peak_tot.GetNoElements());
  
  for (int i=0; i<vr_time_tot.GetNoElements(); i++) {
    v_time_tot_plot(4*i)   = vr_time_tot(i);
    v_time_tot_plot(4*i+1) = vr_time_tot(i);
    v_time_tot_plot(4*i+2) = vr_stop_tot(i);
    v_time_tot_plot(4*i+3) = vr_stop_tot(i);
    
    v_peak_tot_plot(4*i)   = 0.0;
    v_peak_tot_plot(4*i+1) = vr_peak_tot(i);
    v_peak_tot_plot(4*i+2) = vr_peak_tot(i);
    v_peak_tot_plot(4*i+3) = 0.0;
  }
  
  TGraph* g_peak = new TGraph(v_time_tot_plot,v_peak_tot_plot);
  g_peak->SetLineColor(kBlue);
  g_peak->SetMarkerColor(kBlue);
  mg_peak->Add(g_peak,"l");
  

  
  TCanvas* c2 = new TCanvas("c2","Peak luminosity",1200,600);
  c2->SetLogy();
  
  mg_peak->Draw("A");
  //mg_peak->GetHistogram()->SetMinimum(1e-0);
  //mg_peak->GetHistogram()->SetMaximum(1e+3);
  mg_peak->GetHistogram()->Draw("axissame");
  
  mg_peak->GetXaxis()->SetTitle("date");
  mg_peak->GetXaxis()->SetTimeDisplay(1);
  mg_peak->GetXaxis()->SetTimeOffset(0.0);
  mg_peak->GetXaxis()->SetTimeFormat("%b/%Y");
  mg_peak->GetXaxis()->SetLimits(uT2010,uT2010+3*365*86400.0+100.0);
  
  mg_peak->GetYaxis()->SetTitle("Peak Inst Luminosity [10^{30} cm^{-2}s^{-1}]");
  mg_peak->GetYaxis()->SetTitleOffset(1.0);
  
  f_lumiout->cd();
  
  vr_time_tot.Print();
  vr_time_tot.Write("v_fillTime_tot");
  vr_beam_tot.Print();
  vr_beam_tot.Write("v_beamEner_tot");
  vr_stop_tot.Print();
  vr_stop_tot.Write("v_stopTime_tot");
  vr_peak_tot.Print();
  vr_peak_tot.Write("v_peakLumi_tot");
  
  f_lumiout->Close();
  return 0;                             
                               
}



