//------------------------------------------------------------------------//
//                    P0DLI TDC analysis program                          //
//                    Tomasz Wachala, Sep 2012                            //
//------------------------------------------------------------------------//

#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <sstream>
#include <vector>
#include <map>
#include <string>

#include "TROOT.h"
#include "TGraph.h"
#include "TStyle.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TFile.h"
#include "TLegend.h"
#include "TMultiGraph.h"

int s2i( std::string s )
{
        std::stringstream istr(s);
        int i;
        istr >> i;
        return i;
}

int main(int argc, char* argv[]){

  //This is the ROOT's time format that will be used in the plots
  /*This is the ROOT's time format that will be used in the plots
  for date :
      %a abbreviated weekday name
      %b abbreviated month name
      %d day of the month (01-31)
      %m month (01-12)
      %y year without century

    for time :
      %H hour (24-hour clock)
      %I hour (12-hour clock)
      %p local equivalent of AM or PM
      %M minute (00-59)
      %S seconds (00-61)
      %% %
  */
  //const string timeFormat("%d-%H%F1995-01-01 9:00:00");
  std::string timeFormatHead("%d-%H");
  std::string timeFormatApp("%F2000-01-01 00:00:00");
  std::string timeFormat = timeFormatHead+timeFormatApp;

  /*Plot options for graphs
    Graphs can be drawn with the following options:

    "L"	A simple polyline is drawn
    "F"	A fill area is drawn ('CF' draw a smoothed fill area)
    "C"	A smooth Curve is drawn
    "*"	A Star is plotted at each point
    "P"	The current marker is plotted at each point
    "B"	A Bar chart is drawn
    "1"	When a graph is drawn as a bar chart, this option makes the bars start from the bottom of the pad. By default they start at 0.
    "X+"	The X-axis is drawn on the top side of the plot.
    "Y+"	The Y-axis is drawn on the right side of the plot. 
  */
  const std::string plotFormat("LP");

  gROOT->SetStyle("Plain");
  
  gStyle->SetHistLineWidth(1);
  gStyle->SetPalette(1);
  gStyle->SetPaperSize(20,26);
 
  // Turn off some borders
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetDrawBorder(0);
  gStyle->SetCanvasBorderSize(0);
  gStyle->SetFrameBorderSize(0);
  gStyle->SetPadBorderSize(1);
  gStyle->SetTitleBorderSize(0);
  
  // Say it in black and white!
  gStyle->SetAxisColor(1, "xyz");
  gStyle->SetCanvasColor(0);
  gStyle->SetFrameFillColor(0);
  gStyle->SetFrameLineColor(1);
  gStyle->SetHistFillColor(0);
  gStyle->SetHistLineColor(1);
  //gStyle->SetPadColor(1);
  gStyle->SetPadColor(kWhite);
  gStyle->SetStatColor(0);
  gStyle->SetStatTextColor(1);
  gStyle->SetTitleColor(1);
  gStyle->SetTitleTextColor(1);
  gStyle->SetLabelColor(1,"xyz");
  gStyle->SetLabelOffset(0.02);
  // Show functions in red...
  gStyle->SetFuncColor(2);

  //-------------------------------------------------------------------------------------
  // Get command line arguments
  //-------------------------------------------------------------------------------------
  std::vector<std::string> files;
  std::string outfile = "p0d-li-tdc-analysis.root";
  int ampCut = 0;
  int c;
  while ( (c = getopt(argc, argv, "a:o:")) != -1 ) {
    switch (c) {
    case 'a':
      ampCut = s2i(optarg);
      break;
    case 'o':
      outfile = optarg;
      break;
    }
  }
  while (optind<argc) {
    files.push_back( argv[optind++] );
  }

  if( files.size() <1 ) {
    std::cout << "Not enough arguments! " << std::endl;
    std::cout << "Usage: CalibOffline_li-tdc-analysis.exe -o <output> -a <amplitude cut> file1.root file2.root ..." << std::endl;
    std::cout << "  <amplitude cut>: a cut on the P0DLI amplitude setting" << std::endl;
    
    return -1;
  }

  //Load trees from the input files
  TChain *LITree = new TChain("LI");
  for ( std::vector<std::string>::iterator it = files.begin(); it != files.end(); it++ ) {
    LITree->Add( (*it).c_str() );
  }
  int n_entries = LITree->GetEntries();
  //  cout << n_entries << endl;

  //LITree->Print();

  std::cout << ">>>> LI tree entries: " << n_entries << std::endl;

  if(n_entries == 0){
    std::cout << "ERROR: LI tree has NO entries! Bailing..." << std::endl;
    exit(1); 
  }

  float fHighADCMean, fHighADCRMS, fLowADCMean, fLowADCRMS, fTDCMean, fTDCRMS, fTDCMaxBin;
  int fRMM,fTFB,fTripT,fChannel,fP0dule,fBar,fLayer,fPulser,fAmplitude,fFlashes,fRunId;
  double fTime;
  
  LITree->SetBranchAddress("RunNumber",&fRunId);
  LITree->SetBranchAddress("HighADCMean",&fHighADCMean);
  LITree->SetBranchAddress("HighADCRMS",&fHighADCRMS);
  LITree->SetBranchAddress("LowADCMean",&fLowADCMean);
  LITree->SetBranchAddress("LowADCRMS",&fLowADCRMS);
  LITree->SetBranchAddress("TDCMean",&fTDCMean);
  LITree->SetBranchAddress("TDCRMS",&fTDCRMS);
  LITree->SetBranchAddress("TDCMaxBin",&fTDCMaxBin);
  LITree->SetBranchAddress("RMM",&fRMM);
  LITree->SetBranchAddress("TFB",&fTFB);
  LITree->SetBranchAddress("TripT",&fTripT);
  LITree->SetBranchAddress("Channel",&fChannel);
  LITree->SetBranchAddress("P0dule",&fP0dule);
  LITree->SetBranchAddress("Bar",&fBar);
  LITree->SetBranchAddress("Layer",&fLayer);
  LITree->SetBranchAddress("Pulser",&fPulser);
  LITree->SetBranchAddress("Amplitude",&fAmplitude);
  LITree->SetBranchAddress("Flashes",&fFlashes);
  LITree->SetBranchAddress("Time",&fTime);

  ///////Get the number of data points in time////////

  std::cout << ">>>> Getting the number of points... " << std::endl;

  std::map<int,int> Datapoint_Time;
  std::vector<double> Time_Data;
  int counter=0;

  //Get the number of data points in time
  for(int i=0;i<n_entries;i++){

    LITree->GetEntry(i);

    if(fP0dule==0 && fBar==0 && fLayer==0 && fAmplitude==ampCut){
      Datapoint_Time[static_cast<int>(fTime)]=counter;
      Time_Data.push_back(fTime);

      std::cout << std::setprecision(15) << "Time = " << fTime << "<----->" << "Point = " << Datapoint_Time[static_cast<int>(fTime)] << std::endl;

      counter++;
    }

  }

  std::cout << ">>>> Booking histograms..." << std::endl;
 
  const int set_size= counter;

  char temptitle[50];
  TH1I *RMMtdc[6][set_size];
  TH1I *TFBtdc[6][48][set_size];

  for(int rmm=0;rmm<6;rmm++){
    for(int tfb=0;tfb<48;tfb++){
      for(int datapoint=0;datapoint<set_size;datapoint++){
  if(tfb==0){
    sprintf(temptitle,"RMM-%d-Datapoint-%d",rmm,datapoint);
    RMMtdc[rmm][datapoint] = new TH1I(temptitle,temptitle,300,500,800);
  }
  sprintf(temptitle,"RMM-%d-TFB-%d-Datapoint-%d",rmm,tfb,datapoint);
  TFBtdc[rmm][tfb][datapoint] = new TH1I(temptitle,temptitle,300,500,800);
      }
    }
  }

  std::cout << ">>>> Looping over " << n_entries << " entries..." << std::endl;  

  for(int i=0;i<n_entries;i++){

    LITree->GetEntry(i);

    if((i) % (n_entries / 10) == 0) std::cout << ">> Entries processed: " << i << "(" << (100*i)/n_entries << " %) " << std::endl;

    if( fAmplitude==ampCut && fTDCMaxBin != 0){
      //      std::cout << fTDCMaxBin << std::endl;

      RMMtdc[fRMM][Datapoint_Time[static_cast<int>(fTime)]]->Fill(fTDCMaxBin);
      TFBtdc[fRMM][fTFB][Datapoint_Time[static_cast<int>(fTime)]]->Fill(fTDCMaxBin);
      }
  }

  std::cout << "<<<< Done! " << std::endl;

  std::cout << "<<<< Plotting..." << std::endl;

  //Graph FIll and declaration
  TLegend *RMMLeg = new TLegend(0.8,0.8,0.99,0.99);
  TLegend *TFBLeg[6];
  TGraph *RMMtdc_graph[6];
  TGraph *TFBtdc_graph[6][48];

  for(int rmm=0;rmm<6;rmm++){

    TFBLeg[rmm] = new TLegend(0.9,0.5,0.99,0.99);

    for(int tfb=0;tfb<48;tfb++){
      int counter2=0;
      int counter3=0;

      if(tfb==0){
  RMMtdc_graph[rmm] = new TGraph();
  //RMMtdc_graph[rmm]->SetMarkerStyle(8);
  RMMtdc_graph[rmm]->SetMarkerStyle(7);
  RMMtdc_graph[rmm]->SetMarkerColor(rmm*2+1);
  RMMtdc_graph[rmm]->SetLineColor(rmm*2+1);
  RMMtdc_graph[rmm]->SetMinimum(400);
  RMMtdc_graph[rmm]->SetMaximum(800);

  sprintf(temptitle,"RMM-%d",rmm);
  RMMLeg->AddEntry(RMMtdc_graph[rmm],temptitle,"P");

  for(unsigned int j=0;j<Time_Data.size();j++){
    if(RMMtdc[rmm][j]->GetMean()!=0){
      RMMtdc_graph[rmm]->SetPoint(j-counter2,Time_Data[j],RMMtdc[rmm][j]->GetMean());

      //std::cout << j-counter2 << " " << Time_Data[j] << " " << RMMtdc[rmm][j]->GetMean() << std::endl;
    }
    else{
      counter2++;
    }
  }
  sprintf(temptitle,"RMM%d_Amplitude%d",rmm,ampCut);
  RMMtdc_graph[rmm]->SetTitle(temptitle);
      }

      if(TFBtdc[rmm][tfb][0]->GetMean()==0) continue;

      TFBtdc_graph[rmm][tfb] = new TGraph();
      //TFBtdc_graph[rmm][tfb]->SetMarkerStyle(8);
      TFBtdc_graph[rmm][tfb]->SetMarkerStyle(7);
      TFBtdc_graph[rmm][tfb]->SetMarkerColor(tfb+50);
      TFBtdc_graph[rmm][tfb]->SetLineColor(kBlack);
      TFBtdc_graph[rmm][tfb]->SetMinimum(400);
      TFBtdc_graph[rmm][tfb]->SetMaximum(800);

      sprintf(temptitle,"TFB-%d",tfb);
      TFBLeg[rmm]->AddEntry(TFBtdc_graph[rmm][tfb],temptitle,"P");

      //int misscount =0;
      for(unsigned int k=0;k<Time_Data.size();k++){
  if(TFBtdc[rmm][tfb][k]->GetMean()!=0){
    TFBtdc_graph[rmm][tfb]->SetPoint(k-counter3,Time_Data[k],TFBtdc[rmm][tfb][k]->GetMean());
  }
  else{
    counter3++;
  }
      }      
    }
  }

  TMultiGraph *allRMMtdc = new TMultiGraph(Form("RMM_MeanTDC_Amplitude%d",ampCut),Form("RMM_MeanTDC_Amplitude%d",ampCut));
  TMultiGraph *allTFBtdc[6];

  for(int rmm=0;rmm<6;rmm++){
    sprintf(temptitle,"TFB_MeanTDC_RMM%d_Amplitude%d",rmm,ampCut);
    allTFBtdc[rmm] = new TMultiGraph(temptitle,temptitle);
    //allRMMtdc->Add(RMMtdc_graph[rmm],"LP");
    allRMMtdc->Add(RMMtdc_graph[rmm],plotFormat.c_str());
    for(int tfb=0;tfb<48;tfb++){
      if(TFBtdc[rmm][tfb][0]->GetMean()==0) continue;      
      allTFBtdc[rmm]->Add(TFBtdc_graph[rmm][tfb],"P");
    }
  }

  ///////////////////////////////////////////////////////////////
  // Plotting
  ///////////////////////////////////////////////////////////////

  TFile out( outfile.c_str(), "RECREATE");

  TCanvas *c0 = new TCanvas("AllRMM","AllRMM",10,10,1000,700);
  allRMMtdc->Draw("ALP");
  allRMMtdc->SetMinimum(400);
  allRMMtdc->SetMaximum(800);
  allRMMtdc->GetXaxis()->SetTitle(Form("Time (%s)",timeFormatHead.c_str()));
  allRMMtdc->GetYaxis()->SetTitle("Mean TDC");
  allRMMtdc->GetXaxis()->SetTimeDisplay(1);
  allRMMtdc->GetXaxis()->SetTimeFormat(timeFormat.c_str());

  RMMLeg->Draw("SAME");

  /*
  if(savePlots){
    c0->Print((outputDirName+"/MeanTDC_AllRMM.eps").c_str());
    c0->Print((outputDirName+"/MeanTDC_AllRMM.png").c_str());
    c0->Print((outputDirName+"/MeanTDC_AllRMM.C").c_str());
  }
  */
  c0->Write();

  TCanvas *c1 = new TCanvas("RMMTDC","RMMTDC",10,10,1000,700);
  c1->Divide(3,2);

  for(int i=0;i<6;i++){
    c1->cd(i+1);
    RMMtdc_graph[i]->Draw((std::string("A")+plotFormat).c_str());
    RMMtdc_graph[i]->SetMinimum(400);
    RMMtdc_graph[i]->SetMaximum(800);
    RMMtdc_graph[i]->GetXaxis()->SetTitle(Form("Time (%s)",timeFormatHead.c_str()));
    RMMtdc_graph[i]->GetYaxis()->SetTitle("Mean TDC");
    RMMtdc_graph[i]->GetXaxis()->SetTimeDisplay(1);
    RMMtdc_graph[i]->GetXaxis()->SetTimeFormat(timeFormat.c_str());
    RMMLeg->Draw("SAME");
  }

  /*
  if(savePlots){
    c1->Print((outputDirName+"/MeanTDC_RMMTDC.eps").c_str());
    c1->Print((outputDirName+"/MeanTDC_RMMTDC.png").c_str());
    c1->Print((outputDirName+"/MeanTDC_RMMTDC.C").c_str());
  }
  */
  c1->Write();

  TCanvas *c2 = new TCanvas("AllTFB","AllTFB",10,10,1000,700);

  c2->Divide(3,2);

  for(int i=0;i<6;i++){
    c2->cd(i+1);
    allTFBtdc[i]->Draw("ALP");
    allTFBtdc[i]->SetMinimum(400);
    allTFBtdc[i]->SetMaximum(800);
    allTFBtdc[i]->GetXaxis()->SetTitle(Form("Time (%s)",timeFormatHead.c_str()));
    allTFBtdc[i]->GetYaxis()->SetTitle("Mean TDC");
    allTFBtdc[i]->GetXaxis()->SetTimeDisplay(1);
    allTFBtdc[i]->GetXaxis()->SetTimeFormat(timeFormat.c_str());
    TFBLeg[i]->Draw("SAME");
  }

  /*
  if(savePlots){
    c2->Print((outputDirName+"/MeanTDC_AllTFB.eps").c_str());
    c2->Print((outputDirName+"/MeanTDC_AllTFB.png").c_str());
    c2->Print((outputDirName+"/MeanTDC_AllTFB.C").c_str());
  }
  */
  c2->Write(); 

  /*
  TCanvas *c3 = new TCanvas("TDCMaxBin","TDCMaxBin",10,10,1000,700);
  c3->Divide(2,3);
  c3->cd(1);
  LITree->Draw("TDCMaxBin>>TDCMaxBin0",Form("RMM==0 && Amplitude==%d",ampCut));
  TH1F *h1 = ((TH1F*)gDirectory->Get("TDCMaxBin0"));
  h1->GetXaxis()->SetRangeUser(400,800);
  h1->Draw();

  c3->cd(2);
  LITree->Draw("TDCMaxBin>>TDCMaxBin1",Form("RMM==1 && Amplitude==%d",ampCut));
  TH1F *h2 = ((TH1F*)gDirectory->Get("TDCMaxBin1"));
  h2->GetXaxis()->SetRangeUser(400,800);
  h2->Draw();

  c3->cd(3);
  LITree->Draw("TDCMaxBin>>TDCMaxBin2",Form("RMM==2 && Amplitude==%d",ampCut));
  TH1F *h3 = ((TH1F*)gDirectory->Get("TDCMaxBin2"));
  h3->GetXaxis()->SetRangeUser(400,800);
  h3->Draw();

  c3->cd(4);
  LITree->Draw("TDCMaxBin>>TDCMaxBin3",Form("RMM==3 && Amplitude==%d",ampCut));
  TH1F *h4 = ((TH1F*)gDirectory->Get("TDCMaxBin3"));
  h4->GetXaxis()->SetRangeUser(400,800);
  h4->Draw();

  c3->cd(5);
  LITree->Draw("TDCMaxBin>>TDCMaxBin4",Form("RMM==4 && Amplitude==%d",ampCut));
  TH1F *h5 = ((TH1F*)gDirectory->Get("TDCMaxBin4"));
  h5->GetXaxis()->SetRangeUser(400,800);
  h5->Draw();

  c3->cd(6);
  LITree->Draw("TDCMaxBin>>TDCMaxBin5",Form("RMM==5 && Amplitude==%d",ampCut));
  TH1F *h6 = ((TH1F*)gDirectory->Get("TDCMaxBin5"));
  h6->GetXaxis()->SetRangeUser(400,800);
  h6->Draw();

  if(savePlots){
    c3->Print((outputDirName+"/TDCMaxBin.eps").c_str());
    c3->Print((outputDirName+"/TDCMaxBin.png").c_str());
    c3->Print((outputDirName+"/TDCMaxBin.C").c_str());
  }
  */ 
}