//  
//  tempExtract.C
//  macro to extract the temperature raw data, plot it and store it in 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 <fstream>
#include <vector>
#include <map>
#include <stdlib.h>

#include <time.h>


//  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 tempExtract(int argc, char* argv[]);

double getMean(std::vector<double> v,int index, int len=20);
double getStdev(std::vector<double> v,int index, int len=20);


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


//  Executable method
int tempExtract(int argc, char* argv[]){

  std::vector<std::vector<std::string> > a_tt;
  std::vector<std::vector<std::string> > a_it;
  
  const char* homeVar = std::getenv ("CCEHOME");
  
  double uT2010 = 1262304000.0;
  double end2012 = uT2010+3*365*86400.0+86400.0*50.0;
  double endJuly2015 = 1438003626.0;
  double dec2015 = 1448928000.0;
  
  TString dn_lumi(Form("%s/data/temperature",homeVar));
  TString fn_tt(Form("%s/tttemp.txt",dn_lumi.Data()));
  TString fn_it(Form("%s/ittemp.txt",dn_lumi.Data()));
  
  
  std::vector<std::string> v_vname_tt;
  v_vname_tt.push_back("TT_AB5");
  v_vname_tt.push_back("TT_AB6");
  v_vname_tt.push_back("TT_AT5");
  v_vname_tt.push_back("TT_AT6");
  v_vname_tt.push_back("TT_CB5");
  v_vname_tt.push_back("TT_CB6");
  v_vname_tt.push_back("TT_CT5");
  v_vname_tt.push_back("TT_CT6");
  
  std::vector<std::string> v_vname_it;
  v_vname_it.push_back("IT_A1");
  v_vname_it.push_back("IT_B1");
  v_vname_it.push_back("IT_C1");
  v_vname_it.push_back("IT_T1");
  v_vname_it.push_back("IT_A2");
  v_vname_it.push_back("IT_B2");
  v_vname_it.push_back("IT_C2");
  v_vname_it.push_back("IT_T2");
  v_vname_it.push_back("IT_A3");
  v_vname_it.push_back("IT_B3");
  v_vname_it.push_back("IT_C3");
  v_vname_it.push_back("IT_T3");
  
  int nHeadLine = 2;
  int nOffRow   = 1;
  
  std::string dateFormat = "%d/%m/%Y %H:%M:%S";
 
  std::ifstream f_temp_tt(fn_tt.Data());
  std::ifstream f_temp_it(fn_it.Data());
  
  TFile* f_output = new TFile(Form("%s/temperature.root",dn_lumi.Data()),"RECREATE");
  
  // Writing each cell to array for TT temperatures
  while (f_temp_tt) {
    std::string line;
    if (!getline(f_temp_tt,line)) {
      break;
    }
    std::istringstream s_line(line);
    std::vector<std::string> v_line;
    while (s_line) {
      std::string bit;
      if (!getline(s_line,bit,',')) {
        break;
      }
      v_line.push_back(bit);
    }
    a_tt.push_back(v_line);
  }
  
  // Writing each cell to array for IT temperatures
  while (f_temp_it) {
    std::string line;
    if (!getline(f_temp_it,line)) {
      break;
    }
    std::istringstream s_line(line);
    std::vector<std::string> v_line;
    while (s_line) {
      std::string bit;
      if (!getline(s_line,bit,',')) {
        break;
      }
      v_line.push_back(bit);
    }
    a_it.push_back(v_line);
  }
  
  TMultiGraph* mg_TT = new TMultiGraph("mg_TT","mg_TT");
  TMultiGraph* mg_T1 = new TMultiGraph("mg_T1","mg_T1");
  TMultiGraph* mg_T2 = new TMultiGraph("mg_T2","mg_T2");
  TMultiGraph* mg_T3 = new TMultiGraph("mg_T3","mg_T3");
  
  TLegend* leg_TT = new TLegend(0.75,0.75,0.9,0.9);
  TLegend* leg_T1 = new TLegend(0.75,0.6,0.9,0.9);
  TLegend* leg_T2 = new TLegend(0.75,0.6,0.9,0.9);
  TLegend* leg_T3 = new TLegend(0.75,0.6,0.9,0.9);
  
  leg_TT->SetFillColor(0);
  leg_T1->SetFillColor(0);
  leg_T2->SetFillColor(0);
  leg_T3->SetFillColor(0);
  
  // Vectors
  
  std::vector<TVectorD> v_vr_plot_TT;
  std::vector<TVectorD> v_vr_plot_T1;
  std::vector<TVectorD> v_vr_plot_T2;
  std::vector<TVectorD> v_vr_plot_T3;
  
  // Reading temperatures from TT
  for (int i=0; i<v_vname_tt.size(); i++) {
    tm date;
    std::vector<double> v_time;
    std::vector<double> v_temp;
    int indTemp = 2*i+1+nOffRow;
    int indTime = 2*i+nOffRow;
    for (int j=nHeadLine; j<a_tt.size(); j++) {
      if (a_tt[j].size()<=indTemp) {
        Error("tempExtract","Too large index in row %d!   Skip it!",j);
        continue;
      }
      if (a_tt[j][indTemp].compare("")==0) {
        Warning("tempExtract","No temperature in row %d!    Skip it!",j);
        continue;
      }
      if (a_tt[j][indTime].compare("")==0) {
        Warning("tempExtract","No time in row %d!           Skip it!",j);
        continue;
      }
      if (TMath::Abs(atof(a_tt[j][indTemp].c_str())-5.0)>20.0) {
        Warning("tempExtract","Temperature off in row %d!   Skip it!",j);
        continue;
      }
      if (strptime((a_tt[j])[indTime].c_str(), &dateFormat[0], &date) == NULL) {
        Error("tempExtract","Time not readable in row %d! Skip it!");
      }
      v_temp.push_back(atof(a_tt[j][indTemp].c_str()));
      time_t time_val = mktime(&date);
      v_time.push_back((double)time_val);
      Info("tempExtract","%4.2f / %4.2f",atof(a_tt[j][indTemp].c_str()),(double)time_val);
    }
    
    std::vector<double>::iterator it=v_time.begin();
    for (std::vector<double>::iterator j=v_temp.begin(); j!=v_temp.end(); it++, j++) {
      if (getStdev(v_temp,(int)(j-v_temp.begin()),10)<20.0 && TMath::Abs(*j-getMean(v_temp,(int)(j-v_temp.begin()),10))>2.5) {
        v_time.erase(it);
        v_temp.erase(j);
        it--;
        j--;
      }
    }
    
    TVectorD vr_time = getROOTVector(v_time);
    TVectorD vr_temp = getROOTVector(v_temp);
    
    vr_time.Write(Form("v_%s_time",v_vname_tt[i].c_str()));
    vr_temp.Write(Form("v_%s_temp",v_vname_tt[i].c_str()));
    
    if (v_vname_tt[i].compare("TT_AT6")==0) {
      v_vr_plot_TT.push_back(vr_time);
      v_vr_plot_TT.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kRed);
      leg_TT->AddEntry(g_temp,"TT A-Side","l");
      mg_TT->Add(g_temp,"l");
    }else if (v_vname_tt[i].compare("TT_CT6")==0) {
      v_vr_plot_TT.push_back(vr_time);
      v_vr_plot_TT.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kBlue);
      leg_TT->AddEntry(g_temp,"TT C-Side","l");
      mg_TT->Add(g_temp,"l");
    }
  }
  
  // Reading temperatures from IT
  for (int i=0; i<v_vname_it.size(); i++) {
    tm date;
    std::vector<double> v_time;
    std::vector<double> v_temp;
    int indTemp = 2*i+1+nOffRow;
    int indTime = 2*i+nOffRow;
    for (int j=nHeadLine; j<a_it.size(); j++) {
      if (a_it[j].size()<=indTemp) {
        Error("tempExtract","Too large index in row %d!   Skip it!",j);
        continue;
      }
      if (a_it[j][indTemp].compare("")==0) {
        Warning("tempExtract","No temperature in row %d!    Skip it!",j);
        continue;
      }
      if (a_it[j][indTime].compare("")==0) {
        Warning("tempExtract","No time in row %d!           Skip it!",j);
        continue;
      }
      if (TMath::Abs(atof(a_it[j][indTemp].c_str())-10.0)>25.0) {
        Warning("tempExtract","Temperature off in row %d!   Skip it!",j);
        continue;
      }
      if (strptime((a_it[j])[indTime].c_str(), &dateFormat[0], &date) == NULL) {
        Error("tempExtract","Time not readable in row %d! Skip it!");
      }
      v_temp.push_back(atof(a_it[j][indTemp].c_str()));
      time_t time_val = mktime(&date);
      v_time.push_back((double)time_val);
      Info("tempExtract","%4.2f / %4.2f",atof(a_it[j][indTemp].c_str()),(double)time_val);
    }
    std::vector<double>::iterator it=v_time.begin();
    for (std::vector<double>::iterator j=v_temp.begin(); j!=v_temp.end(); it++, j++) {
      if (getStdev(v_temp,(int)(j-v_temp.begin()),10)<20.0 && TMath::Abs(*j-getMean(v_temp,(int)(j-v_temp.begin()),10))>2.5) {
        v_time.erase(it);
        v_temp.erase(j);
        it--;
        j--;
      }
    }
    
    TVectorD vr_time = getROOTVector(v_time);
    TVectorD vr_temp = getROOTVector(v_temp);
    
    vr_time.Write(Form("v_%s_time",v_vname_it[i].c_str()));
    vr_temp.Write(Form("v_%s_temp",v_vname_it[i].c_str()));
    
    
    if (v_vname_it[i].compare("IT_A1")==0) {
      v_vr_plot_T1.push_back(vr_time);
      v_vr_plot_T1.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kRed);
      leg_T1->AddEntry(g_temp,"T1 Box A","l");
      mg_T1->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_B1")==0) {
      v_vr_plot_T1.push_back(vr_time);
      v_vr_plot_T1.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kBlue);
      leg_T1->AddEntry(g_temp,"T1 Box B","l");
      mg_T1->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_C1")==0) {
      v_vr_plot_T2.push_back(vr_time);
      v_vr_plot_T2.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kGreen+3);
      leg_T1->AddEntry(g_temp,"T1 Box C","l");
      mg_T1->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_T1")==0) {
      v_vr_plot_T1.push_back(vr_time);
      v_vr_plot_T1.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kBlack);
      leg_T1->AddEntry(g_temp,"T1 Box T","l");
      mg_T1->Add(g_temp,"l");
    }
    
    
    if (v_vname_it[i].compare("IT_A2")==0) {
      v_vr_plot_T2.push_back(vr_time);
      v_vr_plot_T2.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kRed);
      leg_T2->AddEntry(g_temp,"T2 Box A","l");
      mg_T2->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_B2")==0) {
      v_vr_plot_T2.push_back(vr_time);
      v_vr_plot_T2.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kBlue);
      leg_T2->AddEntry(g_temp,"T2 Box B","l");
      mg_T2->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_C2")==0) {
      v_vr_plot_T2.push_back(vr_time);
      v_vr_plot_T2.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kGreen+3);
      leg_T2->AddEntry(g_temp,"T2 Box C","l");
      mg_T2->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_T2")==0) {
      v_vr_plot_T2.push_back(vr_time);
      v_vr_plot_T2.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kBlack);
      leg_T2->AddEntry(g_temp,"T2 Box T","l");
      mg_T2->Add(g_temp,"l");
    }
    
    
    
    if (v_vname_it[i].compare("IT_A3")==0) {
      v_vr_plot_T3.push_back(vr_time);
      v_vr_plot_T3.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kRed);
      leg_T3->AddEntry(g_temp,"T3 Box A","l");
      mg_T3->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_B3")==0) {
      v_vr_plot_T3.push_back(vr_time);
      v_vr_plot_T3.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kBlue);
      leg_T3->AddEntry(g_temp,"T3 Box B","l");
      mg_T3->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_C3")==0) {
      v_vr_plot_T3.push_back(vr_time);
      v_vr_plot_T3.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kGreen+3);
      leg_T3->AddEntry(g_temp,"T3 Box C","l");
      mg_T3->Add(g_temp,"l");
    }else if (v_vname_it[i].compare("IT_T3")==0) {
      v_vr_plot_T3.push_back(vr_time);
      v_vr_plot_T3.push_back(vr_temp);
      TGraph* g_temp = new TGraph(vr_time,vr_temp);
      g_temp->SetLineColor(kBlack);
      leg_T3->AddEntry(g_temp,"T3 Box T","l");
      mg_T3->Add(g_temp,"l");
    }
  }
  
  TCanvas* c_TT = new TCanvas("c_TT","c_TT",1500,600);
  
  mg_TT->Draw("A");
  mg_TT->GetHistogram()->SetMinimum(-10.0);
  mg_TT->GetHistogram()->SetMaximum( 30.0);
  mg_TT->GetXaxis()->SetTitle("date");
  mg_TT->GetXaxis()->SetTimeDisplay(1);
  mg_TT->GetXaxis()->SetTimeOffset(0.0);
  mg_TT->GetXaxis()->SetTimeFormat("%b/%Y");
  mg_TT->GetXaxis()->SetLimits(uT2010, dec2015);
  
  mg_TT->GetYaxis()->SetTitle("Temperature [#circ C]");
  mg_TT->GetYaxis()->SetTitleOffset(1.0);
  
  leg_TT->Draw();
  
  mg_TT->GetHistogram()->Draw("axissame");

  
  TCanvas* c_T1 = new TCanvas("c_T1","c_T1",1500,600);
    
  mg_T1->Draw("A");
  mg_T1->GetHistogram()->SetMinimum(  0.0);
  mg_T1->GetHistogram()->SetMaximum( 40.0);
  mg_T1->GetXaxis()->SetTitle("date");
  mg_T1->GetXaxis()->SetTimeDisplay(1);
  mg_T1->GetXaxis()->SetTimeOffset(0.0);
  mg_T1->GetXaxis()->SetTimeFormat("%b/%Y");
  mg_T1->GetXaxis()->SetLimits(uT2010, dec2015);
  
  mg_T1->GetYaxis()->SetTitle("Temperature [#circ C]");
  mg_T1->GetYaxis()->SetTitleOffset(1.0);
  
  leg_T1->Draw();
  
  mg_T1->GetHistogram()->Draw("axissame");
  
  TCanvas* c_T2 = new TCanvas("c_T2","c_T2",1500,600);
  
  mg_T2->Draw("A");
  mg_T2->GetHistogram()->SetMinimum(  0.0);
  mg_T2->GetHistogram()->SetMaximum( 40.0);
  mg_T2->GetXaxis()->SetTitle("date");
  mg_T2->GetXaxis()->SetTimeDisplay(1);
  mg_T2->GetXaxis()->SetTimeOffset(0.0);
  mg_T2->GetXaxis()->SetTimeFormat("%b/%Y");
  mg_T2->GetXaxis()->SetLimits(uT2010, dec2015);
  
  mg_T2->GetYaxis()->SetTitle("Temperature [#circ C]");
  mg_T2->GetYaxis()->SetTitleOffset(1.0);
  
  leg_T2->Draw();
  
  mg_T2->GetHistogram()->Draw("axissame");
  
  TCanvas* c_T3 = new TCanvas("c_T3","c_T3",1500,600);
  
  mg_T3->Draw("A");
  mg_T3->GetHistogram()->SetMinimum(  0.0);
  mg_T3->GetHistogram()->SetMaximum( 40.0);
  mg_T3->GetXaxis()->SetTitle("date");
  mg_T3->GetXaxis()->SetTimeDisplay(1);
  mg_T3->GetXaxis()->SetTimeOffset(0.0);
  mg_T3->GetXaxis()->SetTimeFormat("%b/%Y");
  mg_T3->GetXaxis()->SetLimits(uT2010, dec2015);
  
  mg_T3->GetYaxis()->SetTitle("Temperature [#circ C]");
  mg_T3->GetYaxis()->SetTitleOffset(1.0);
  
  leg_T3->Draw();
  
  mg_T3->GetHistogram()->Draw("axissame");

  
  f_output->Close();
  return 0;                             
                               
}

double getMean(std::vector<double> v,int index, int len){
  int start = index-len>0?index-len:0;
  int stop  = index+len>v.size()?v.size():index+len;
  
  int n = 0;
  double mean = 0;
  
  for (int i=start; i<stop; i++) {
    mean += v[i];
    n++;
  }
  
  return n>0?mean/n:0.0;
  
}

double getStdev(std::vector<double> v,int index, int len){
  int start = index-len>0?index-len:0;
  int stop  = index+len>v.size()?v.size():index+len;
  
  int n = 0;
  double mean = getMean(v,index,len);
  double x2 = 0.0;
  
  for (int i=start; i<stop; i++) {
    x2 += v[i]*v[i];
    n++;
  }
  
  return n>0?TMath::Sqrt(x2/n-mean*mean):0.0;
  
}