//  
//  lumiExtract.C
//  macro to extract the luminosity raw data, plot it and store it into TVectorD
//
//  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 "../../include/incBasicROOT.h"
#include "../../include/incDrawROOT.h"
#include "../../include/incGraphROOT.h"
#include "../../include/incIOROOT.h"
#include "../../include/incLinAlgROOT.h"



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

#include "../../include/lhcbstyle.C"








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

int main(int argc, char* argv[]){
  
  TStyle* style = lhcbStyle();
  
  int argc_copy = argc;
  
  TApplication* theApp = new TApplication("Analysis", &argc, argv);
  argc = theApp->Argc();
  argv = theApp->Argv();
  
  bool runBatch = false;
  
  for (int i = 1; i<argc; i++) {
    if (strcmp(argv[i],"-b")==0) {
      runBatch = true;
      gROOT->SetBatch();
    }
  }  
  
  int exit = lumiExtract(argc,argv);
  
  Printf("End");
  if (!runBatch) {
    theApp->Run(!runBatch);
  }
  
  return exit;
}


//  Executable method
int lumiExtract(int argc, char* argv[]){
  
  TMultiGraph *mg_lumi = new TMultiGraph("mg_lumi","mg_lumi");

  const char* homeVar = std::getenv ("CCEHOME");

  TString dn_lumi(Form("%s/data/lumi",homeVar));
  
  std::vector<TString> v_fn_lumi;
  v_fn_lumi.push_back(TString(Form("%s/lumi2010.root",dn_lumi.Data())));
  v_fn_lumi.push_back(TString(Form("%s/lumi2011.root",dn_lumi.Data())));
  v_fn_lumi.push_back(TString(Form("%s/lumi2012.root",dn_lumi.Data())));
  v_fn_lumi.push_back(TString(Form("%s/lumi2015.root",dn_lumi.Data())));
  v_fn_lumi.push_back(TString(Form("%s/lumi2016.root",dn_lumi.Data())));
  v_fn_lumi.push_back(TString(Form("%s/lumi2017.root",dn_lumi.Data())));
  
  double uT2010 = 1262304000.0;
  
  std::vector<double> v_timeCorrFac;
  v_timeCorrFac.push_back(86400.0);
  v_timeCorrFac.push_back(1.0);
  v_timeCorrFac.push_back(1.0);
  v_timeCorrFac.push_back(1.0);
  v_timeCorrFac.push_back(1.0);
  v_timeCorrFac.push_back(1.0);
  
  std::vector<double> v_timeCorrScale;
  v_timeCorrScale.push_back(uT2010);
  v_timeCorrScale.push_back(0.0);
  v_timeCorrScale.push_back(0.0);
  v_timeCorrScale.push_back(0.0);
  v_timeCorrScale.push_back(0.0);
  v_timeCorrScale.push_back(0.0);
  
  std::vector<double> v_lumiCorrFac;
  v_lumiCorrFac.push_back(1.0);
  v_lumiCorrFac.push_back(1.0);
  v_lumiCorrFac.push_back(1.0);
  v_lumiCorrFac.push_back(1.0e-9);
  v_lumiCorrFac.push_back(1.0e-9);
  //v_lumiCorrFac.push_back(1.0e-3);
  v_lumiCorrFac.push_back(1.0e-9);
  
  std::vector<int> v_color;
  v_color.push_back(kGreen+3);
  v_color.push_back(kBlue);
  v_color.push_back(kRed);
  v_color.push_back(kOrange);
  v_color.push_back(kCyan);
  v_color.push_back(kBlack);
  
  std::vector<double> v_cms;
  v_cms.push_back(7.0);
  v_cms.push_back(7.0);
  v_cms.push_back(8.0);
  v_cms.push_back(13.0);
  v_cms.push_back(13.0);
  v_cms.push_back(13.0);
  
  std::vector<double> v_lumiScale;
  v_lumiScale.push_back(10.0);
  v_lumiScale.push_back(1.0);
  v_lumiScale.push_back(1.0);
  //v_lumiScale.push_back(10.0);
  v_lumiScale.push_back(1.0);
  v_lumiScale.push_back(1.0);
  v_lumiScale.push_back(1.0);
  
  
  TFile* f_lumiout = new TFile(Form("%s/lumi.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_deli *= v_lumiCorrFac[i];
    v_reco *= v_lumiCorrFac[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;
      }
    }
    
    v_stop += v_time;
    
    TVectorF v_time_plot(v_time);
    
    if (i>2) {
      v_time_plot -= 365*86400.0; // subtraction for 2013
      v_time_plot -= 365*86400.0; // subtraction for 2014
    }
    v_time_plot -= i*365*86400.0; // conversion sec -> year
    
    TGraph* g_deli = new TGraph(v_time_plot,v_deli_plot);
    g_deli->SetLineColor(v_color[i]);
    g_deli->SetLineStyle(7);
    mg_lumi->Add(g_deli,"l");
    
    TGraph* g_reco = new TGraph(v_time_plot,v_reco_plot);
    g_reco->SetLineColor(v_color[i]);
    g_reco->SetLineStyle(0);
    mg_lumi->Add(g_reco,"l");
    
    if (i>2 && v_lumiScale[i]==1.0) {
      leg_lumi->AddEntry(g_reco,Form("%d",2012+i),"l");
    }else if (i>2) {
      leg_lumi->AddEntry(g_reco,Form("%d (%dx)",2012+i,(int)v_lumiScale[i]),"l");
    }else if (v_lumiScale[i]==1.0) {
      leg_lumi->AddEntry(g_reco,Form("%d",2010+i),"l");
    }else{
      leg_lumi->AddEntry(g_reco,Form("%d (%dx)",2010+i,(int)v_lumiScale[i]),"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("%b");
  mg_lumi->GetXaxis()->SetLimits(uT2010,uT2010+365*86400.0);
  
  mg_lumi->GetYaxis()->SetTitle("Luminosity [pb^{-1}]");
  mg_lumi->GetYaxis()->SetTitleOffset(1.5);
  mg_lumi->GetHistogram()->SetMaximum(2500.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(kRed);
  g_peak->SetMarkerColor(kRed);
  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+7*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");

  Info("lumiExtract", "Data written to %s/lumi.root", dn_lumi.Data());
  
  f_lumiout->Close();
  return EXIT_SUCCESS;                             
                               
}