//  
//  systVoltage.C
//  <Enter discription here>
//
//  Created by Christian Elsasser on 31.05.13.
//  University of Zurich, elsasser@cern.ch
//  Copyright only for commercial use
//

//  General include
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <stdlib.h>
#include <limits>
#include <string>
#include <time.h>


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

#include "TExMap.h"


#include "basicROOTTools.h"


#include "lhcbstyle.C"


#include "deplTool.h"





int systVoltage(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 = systVoltage(argc,argv);
  
  
  
  Printf("End");
  if (!runBatch) {
    theApp->Run(!runBatch);
  }
  
  
  return exit;
}


//  Executable method
int systVoltage(int argc, char* argv[]){
  
  TVectorD v_res(5);
  
  gStyle->SetPadLeftMargin(0.20);
  gStyle->SetTitleOffset(1.5,"Y");
  gStyle->SetTitleOffset(1.0,"X");

  const char* diskVar = std::getenv ("DISK");
  const char* homeVar = std::getenv ("HOME");
  
  TString det("TT");
  TString lay("TTaU");
  
  
  // Reading command line arguments
  for (int i=1; i<argc; i++) {
    if (strcmp(argv[i],"-d")==0 && i+1<argc) { // Detector
      det = argv[++i];
    }else if (strcmp(argv[i],"-l")==0 && i+1<argc) { // Layer
      lay = argv[++i];
    }
  }
  
  
  
  TString dn_fill(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_fill(Form("%s/lumi.root",dn_fill.Data()));
  
  TString dn_deplV(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_deplV(Form("%s/vdepl.root",dn_deplV.Data()));

  TFile* f_fill   = TFile::Open(fn_fill.Data());
  TFile* f_deplV  = TFile::Open(fn_deplV.Data());
  
  TVectorD* vp_fill  = (TVectorD*)f_fill->Get("v_fill");
  TVectorD* vp_sector  = (TVectorD*)f_deplV->Get(Form("v_sector_%s",lay.Data()));
  
  if (!vp_sector || !vp_fill) {
    Error("systVoltage","Fill and/or sector vector are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_sector  = *vp_sector;
  TVectorD v_fill    = *vp_fill;
  
  TString dn_ratio(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db"
                        ,homeVar));
  TString fn_syst(Form("%s/systVolt.root",
                       dn_ratio.Data()));
  
  TFile* f_map = TFile::Open(fn_syst.Data());
  TExMap* m_lin = (TExMap*)f_map->Get(Form("m_lin_%s",lay.Data())); // Deviations from linear approach
  TExMap* m_sig = (TExMap*)f_map->Get(Form("m_sig_%s",lay.Data())); // Deviations from sigmoid approach

  TH1F* h_lin = new TH1F("h_lin","systematic histogram",
                         50,-50.0,50.0);
  h_lin->SetXTitle("#Delta#it{V} [V]");
  h_lin->SetYTitle(Form("Entries / (%4.2f V)",h_lin->GetBinWidth(1)));
  
  TH1F* h_sig = (TH1F*)h_lin->Clone("h_sig");
  
  for (int i=0; i<v_fill.GetNoElements(); i++) {
    for (int j=0; j<v_sector.GetNoElements(); j++) {
      int sector = (int)v_sector(j);
      int fill   = (int)v_fill(i);
      if (STTool::IsExcluded(sector)) {
        Warning("systVoltage","Sector excluded");
        continue;
      }
      double v_sig = *(reinterpret_cast<double*>(&(*m_sig)(fill*10000+sector)));
      if (TMath::Abs(v_sig)>1e-6) {
        h_sig->Fill(v_sig);
      }
      double v_lin = *(reinterpret_cast<double*>(&(*m_lin)(fill*10000+sector)));
      if (TMath::Abs(v_lin)>1e-6) {
        h_lin->Fill(v_lin);
      }
      
    }
  }
  
  Info("systVoltage","Systematic linear:  %8.4f V",h_lin->GetRMS());
  Info("systVoltage","Systematic sigmoid: %8.4f V",h_sig->GetRMS());
  
  TCanvas* c_lin = new TCanvas("c_lin","Canvas linear",800,600);
  c_lin->cd();
  h_lin->Draw("PE");
  
  TCanvas* c_sig = new TCanvas("c_sig","Canvas sigmoid",800,600);
  c_sig->cd();
  h_sig->Draw("PE");
  
  
  Printf("=========================================\n");

  return EXIT_SUCCESS;
}