//  
//  drawScheme.C
//  Draws color map to show the difference for the different read-out sectors with respect to the initial depletion voltage
//
//  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 "incDrawROOT.h"
#include "incGraphROOT.h"
#include "incHistROOT.h"
#include "incIOROOT.h"
#include "incRooFit.h"


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


#include "lhcbstyle.C"


#include "deplTool.h"
#include "hamburgTool.h"





int drawScheme(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();
    }
  }
  
  
  const Int_t NRGBs = 3;
  const Int_t NCont = 999;
  Double_t stops[NRGBs] = { 0.00, 0.50, 1.00 };
  Double_t red[NRGBs]   = { 0.00, 0.04, 0.80 };
  Double_t green[NRGBs] = { 0.00, 0.00, 0.00 };
  Double_t blue[NRGBs]  = { 0.60, 0.03, 0.00 };
  
  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
  style->SetNumberContours(NCont);

  
  style->SetPadLeftMargin(0.12);
  style->SetPadRightMargin(0.22);
  style->SetPadBottomMargin(0.15);
  style->SetPadTopMargin(0.05);
  
  int exit = drawScheme(argc,argv);
  
  
  
  Printf("End");
  if (!runBatch) {
    theApp->Run(!runBatch);
  }
  
  
  return exit;
}

// -d detector (TT/IT)
// -f fill number
// -ns Do not save the value in the root file


//  Executable method
int drawScheme(int argc, char* argv[]){
  
  
  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");
  int fill  = 1944;
  bool notSave = false;
  
  
  for (int i=1; i<argc; i++) {
    if (strcmp(argv[i],"-d")==0 && i+1<argc) {
      det = argv[++i];
    }else if (strcmp(argv[i],"-f")==0 && i+1<argc) {
      fill = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-ns")==0) {
      notSave = true;
    }
  }
  
  
  if (det.EqualTo("IT")) {
    lay = "T3X2";
  }else{
    lay = "TTaU";
  }
  
  
  // Directory and file name for the output graph
  TString dn_graph(Form("%s/data/ST/graphics/%s/%s",
                        diskVar,det.Data(),lay.Data()));
  TString fn_graph(Form("%s/map%d.pdf",
                        dn_graph.Data(),fill));
  
  // Directory and file name for the input data
  TString dn_data(Form("%s/data/ST/Aging",diskVar));
  TString fn_data(Form("%s/CCEScan.root",dn_data.Data()));
  
  TFile* f_data = TFile::Open(fn_data.Data());
  
  
  // Directory and file names for the measured Vdepl values, the fluence values,
  // the Vdepl values after production and the luminosity (i.e. development of the flux over time)
  
  
  TString dn_hit(Form("%s/data/ST/Aging",diskVar));
  TString fn_hit(Form("%s/%sHits.root",dn_hit.Data(),det.Data()));
                 
  TString tn_hit(Form("STHitIntegrator/%s",lay.Data()));
  
  TFile* f_hit = TFile::Open(fn_hit.Data());
  
  TTree* t_hit = (TTree*)f_hit->Get(tn_hit.Data());
                 
  TString dn_sector(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_sector(Form("%s/fluence.root",dn_sector.Data()));
  
  TFile* f_sector = TFile::Open(fn_sector.Data());
  
  TString dn_fill(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_fill(Form("%s/lumi.root",dn_fill.Data()));
  
  TFile* f_fill   = TFile::Open(fn_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_deplV  = TFile::Open(fn_deplV.Data());
  
  TVectorD* vp_secVdepl  = (TVectorD*)f_deplV->Get(Form("v_sector_%s",lay.Data()));
  TVectorD* vp_vdepl     = (TVectorD*)f_deplV->Get(Form("v_vdepl_%s",lay.Data()));
  
  if (!vp_secVdepl || !vp_vdepl) {
    Error("deplVphiN","Depletion Voltage vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_secVdepl  = *vp_secVdepl;
  TVectorD v_vdepl     = *vp_vdepl;
  
  
  // Defines the considered radii ranges
  const Int_t nRTT = 3;
  Double_t rTT[nRTT] = {-1.0,75.0,45.0};
  

  TVectorD* vp_fill  = (TVectorD*)f_fill->Get("v_fill");
  if (!vp_fill) {
    Error("drawScheme","Fill vector is not available");
    return EXIT_FAILURE;
  }
  TVectorD v_fill    = *vp_fill;
  
  
  TVectorD* vp_sector  = (TVectorD*)f_sector->Get(Form("v_sector_%s",lay.Data()));
  if (!vp_sector) {
    Error("drawScheme","Sector vector is not available");
    return EXIT_FAILURE;
  }
  TVectorD v_sector  = *vp_sector;
  
  
  Int_t effnR = nRTT;
  
  // No radii considered in case of IT
  if (det.EqualTo("IT")) {
    effnR = 1;
  }
  
  if (!t_hit) {
    Error("drawScheme","Hit tree %s is not available",tn_hit.Data());
    return EXIT_FAILURE;
  }
  
  // Define the number of scaned points
  int nEntries = t_hit->GetEntries();
  nEntries = 10000000;
  double x,y;
  UInt_t sec;
  // Extraction of the hit coordinates and the corresponding sector
  t_hit->SetBranchAddress("x",  &x  );
  t_hit->SetBranchAddress("y",  &y  );
  t_hit->SetBranchAddress("Sec",&sec);
  
  

  TH2F* h_map = new TH2F("h_map","Map histo",
                         200,-20,20,
                         200,-20,20);
  h_map->SetXTitle("#it{x} [cm]");
  h_map->SetYTitle("#it{y} [cm]");
  h_map->SetZTitle("#Delta#it{V}_{depl} [V]");
  for (int xbin=1; xbin<=h_map->GetNbinsX(); xbin++) {
    for (int ybin=1; ybin<=h_map->GetNbinsX(); ybin++) {
      h_map->SetBinContent(xbin,ybin,-1.0);
    }
  }
  
  
  h_map->SetMinimum(-0.01);
  h_map->SetMaximum(50.0);
  
  
  int fill0 = 1944;
  
  Info("drawScheme","Processing %d entries ...",nEntries);
  
  for (int i=0; i<nEntries; i++) {
    t_hit->GetEntry(i);
    TString vn_data(Form("%s/%s/%d/%d/v_volt_val%d",
                         det.Data(),lay.Data(),sec,fill,val));
    TString vn_data0(Form("%s/%s/%d/%d/v_volt_val%d",
                         det.Data(),lay.Data(),sec,fill0,val));
    // x and y are the coordinates of the hit
    double r2 = x*x+y*y;
    TString vn_rad = "";
    int bin = h_map->FindBin(x,y);
    if (h_map->GetBinContent(bin)>0.0) {
      continue;
    }
    
    for (int k=1; k<effnR; k++) {
      if (r2<rTT[k]*rTT[k]/100.0) {
        vn_rad = Form("_r%d",(int)rTT[k]);
      }
    }
    vn_data  += vn_rad;
    vn_data0 += vn_rad;
    
    int ind          = getClosestIndex(v_secVdepl,(double)sec);
    double vdeplProd = v_vdepl(ind);
    
    
    TVectorD* vp_data  = (TVectorD*)f_data->Get(vn_data.Data() );
    TVectorD* vp_data0 = (TVectorD*)f_data->Get(vn_data0.Data());
    
    if (!vp_data || !vp_data0){
      Warning("drawScheme","No vector found for %s",vn_data.Data());
      continue;
    }
    // Assume that the depletion voltage is decreasing with respect to
    // its initial value
    double deltaV = TMath::Max(vdeplProd-(*vp_data)(0),0.0);
    
    // Set the corresponding hit in the color map to the proper value
    h_map->SetBinContent(bin,deltaV);
    if (i%10000==0) {
      Info("drawScheme","%d entries processed",i);
    }
  }
  
  TCanvas* c_map = new TCanvas("c_map","Map canvas",1000,800);
  
  h_map->Draw("colz");
  
  if (!notSave) {
    // Save picture
    gSystem->Exec(Form("mkdir -p %s",dn_graph.Data()));
    gSystem->Exec(Form("rm %s &> /dev/null",fn_graph.Data()));
    c_map->SaveAs(fn_graph.Data());
  }
  
  Printf("=========================================\n");

  return EXIT_SUCCESS;
}