//  
//  pulseFit.C
//  Obtaining the integral of the pulse shape (between the roots)
//  representing the amount of collected charge
//
//  Created by Christian Elsasser on 31.05.13.
//  University of Zurich, elsasser@cern.ch
//  Copyright only for commercial use
//
//  Compile with:
//  g++ -o pulseFit pulseFit.C `root-config --glibs --cflags` -lRooFit -lRooFitCore -lEG
//

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


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


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


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


#include "deplTool.h"


// -d detector (TT/IT)
// -l layer (TTaX,TTaU,TTbV,TTbX, T[1-3]{X1,U,V,X2})
// -s read-out sector number
// -f fill number
// -c Voltage step [0-10]
// -r consider only hit in a distance closer than r mm from the beam axis
// -v number of summed up strips
// -ns Do not save the value in the root file
// -ss Perform a temporary save


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


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

  const char* diskVar = std::getenv ("DISK");
  
  int fill      = 2083;
  TString det("TT");
  TString lay("TTaU");
  int sector    = 2634;
  int caliStep  = 0;
  int val       = 3;
  bool notSave  = false;
  bool temSave  = false;
  double radius = -10.0;
  
  for (int i=1; i<argc; i++) {
    if (strcmp(argv[i],"-d")==0 && i+1<argc) {
      det = argv[++i];
    }else if (strcmp(argv[i],"-l")==0 && i+1<argc) {
      lay = argv[++i];
    }else if (strcmp(argv[i],"-s")==0 && i+1<argc) {
      sector = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-f")==0 && i+1<argc) {
      fill = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-v")==0 && i+1<argc) {
      val = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-c")==0 && i+1<argc) {
      caliStep = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-r")==0 && i+1<argc) {
      radius = std::atof(argv[++i]);
    }else if (strcmp(argv[i],"-ns")==0) {
      notSave = true;
    }else if (strcmp(argv[i],"-ss")==0) {
      temSave = true;
    }
  }
  
  
  
  int voltStep  = caliStep;
  
  double volt   = 0.0;
  
  
  // Vectors to store the ADC and time values obtained in the landau fit
  std::vector<double> v_time_val;
  std::vector<double> v_time_err;
  std::vector<double> v_adc_val;
  std::vector<double> v_adc_err;
  
  
  
  double time   = STTool::GetTime(caliStep, fill);
  
  
  // Get the correct voltage value
  if (det.EqualTo("TT")) {
    volt = STTool::GetTTVoltage(caliStep*6);
  }else{
    volt = STTool::GetITVoltage(caliStep*6);
  }
  
  // Name of the vector where the result is stored
  TString vn_result(Form("v_pulse%d_val%d",caliStep,val));
  
  // Directory and file name for the plots
  TString dn_graph(Form("%s/data/ST/graphics/%s/%s/%d/%d",
                        diskVar,det.Data(),lay.Data(),sector,fill));
  TString fn_graph_p(Form("%s/pulse_%d_val%d",
                          dn_graph.Data(),caliStep,val));
  
  
  // Dealing with the radius
  if (radius<0.0) {
  }else{
    vn_result  +=Form("_r%d",(int)radius);
    fn_graph_p +=Form("_r%d",(int)radius);
  }
  
  fn_graph_p += ".pdf";
  
  
  Printf("=========================================");
  Info("pulseFit","Analyze: ");
  Printf("               Sector:    %d",sector);
  if (radius>0.0) {
    Printf("               Radius:    %4.2f mm",radius);
  }
  Printf("               Fill:      %d",fill);
  Printf("               Detector:  %s",det.Data());
  Printf("               Layer:     %s",lay.Data());
  Printf("               Voltage:   %4.2f V",volt);
  
  
  // 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()));
  TString fn_dummy(Form("%s/check",dn_data.Data()));
  
  if (temSave) {
    fn_data= TString(Form("%s/temp_%s_%s_%d_%d_r%d.root",
                          dn_data.Data(),det.Data(),lay.Data(),
                          sector,fill,(int)radius));
  }

  TFile* f_data       = TFile::Open(fn_data.Data());
  
  if (!f_data) {
    Error("pulseFit","Data file does not exist!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  // Loading the values from the Landau distribution fits and the corresponding time
  Info("pulseFit","Loading data:");
  TString dn_input(Form("%s/%s/%d/%d",
                        det.Data(),
                        lay.Data(),
                        sector,
                        fill));
  for (int i=6*caliStep; i<6*caliStep+6; i++) {
    TString vn_input(Form("%s/v_%d_val%d",dn_input.Data(),i,val));
    if (radius>0.0) {
      vn_input  +=Form("_r%d",(int)radius);
    }
    TVectorD* v_input = (TVectorD*)f_data->Get(vn_input.Data());
    if (!v_input) {
      Warning("pulseFit","Vector %s not available",vn_input.Data());
      continue;
    }
    Printf("        %8.2f ns: %8.4f+/-%8.4f",STTool::GetTime(i,fill),(*v_input)(0),(*v_input)(1));
    if (i==6*caliStep) {
      v_res(6) = (*v_input)(0);
      v_res(7) = (*v_input)(1);
    }
    
    // Excluding steps bad data quality
    if (i==0 && (fill==2797 || fill==3108)) {
      Info("pulseFit","Bad calibration step! Skip it!");
      continue;
    }
    
    v_adc_val.push_back((*v_input)(0));
    v_adc_err.push_back((*v_input)(1));
    v_time_val.push_back(STTool::GetTime(i,fill));
    v_time_err.push_back(0.0);
    
  }
  
  if (v_adc_val.size()<3) {
    Error("pulseFit","Too few (%d) data points",v_adc_val.size());
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  
  // Transform STL vectors into ROOT vectors
  TVectorD vr_adc_val  = getROOTVector(v_adc_val);
  TVectorD vr_adc_err  = getROOTVector(v_adc_err);
  TVectorD vr_time_val = getROOTVector(v_time_val);
  TVectorD vr_time_err = getROOTVector(v_time_err);
  
  TGraphErrors* ge_pulse = new TGraphErrors(vr_time_val,vr_adc_val,vr_time_err,vr_adc_err);
  
  // Fit the pulse shape
  TF1* fu_pulse = STTool::halfGauss;
  ge_pulse->Fit(fu_pulse,"M0Q","0");
  
  double chi2ndf = fu_pulse->GetChisquare()/fu_pulse->GetNDF();
  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse->GetChisquare(),fu_pulse->GetNDF());

  
  // Make PDG-like correction of the uncertainties
  //Commented out by Michele - 8/31 2017
  //if (chi2ndf>1.0) {
  //  Info("pulseFit","scaling");
  //  vr_adc_err*=TMath::Sqrt(chi2ndf);
  //}
  //Michele

  TMultiGraph* mg_pulse = new TMultiGraph("mg_pulse","Pulse multigraph");
  
  ge_pulse = new TGraphErrors(vr_time_val,vr_adc_val,vr_time_err,vr_adc_err);
  ge_pulse->Fit(fu_pulse,"M0Q","0");
  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse->GetChisquare(),fu_pulse->GetNDF());



  //Evaluation of the systematic error on the integral of the pulse shape using a polinomial
  //Define the interval were most likely the zeros should be
  //TO DO: Set min and max from the time vector
  double low_limit = -50;
  double high_limit = 50;
  //HERE
  //POLY3

  cout << "================ poly3 ====================" << endl<<endl;

  TF1 *poly3 = new TF1("poly3","[0] + [1]*x +[2]*x**2 + [3]*x**3 ",-40.,40.);
  poly3->SetLineColor(1);
  ge_pulse->Fit("poly3","M");
  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly3->GetChisquare(),poly3->GetNDF());


  std::vector<double> v_max_p3;
  std::vector<double> v_time_srtd;
  bool b_max_p3 = true;

  //Copy the vector of interest
  for (int i=0; i<v_time_val.size(); i++)
    v_time_srtd.push_back(v_time_val[i]);
  //sort it
  std::sort(v_time_srtd.begin(), v_time_srtd.end()); 
  double max_tmp = 0;
  //fill a vector with local maxima
  for (int i=0; i<v_time_srtd.size()-1; i++){

    //Get local Maximum
    max_tmp = poly3->GetMaximum(v_time_srtd[i], v_time_srtd[i+1]);

    //Checks on the maximum -> do not accept maxima at the boundary
    if( (max_tmp <=  poly3->Eval(v_time_srtd[i])) || (max_tmp <=  poly3->Eval(v_time_srtd[i+1]))   ){
      continue;
    }else{
    v_max_p3.push_back(max_tmp);
    }
  }

  if (v_max_p3.size() == 1){

    //Final check

    cout << "Found the maximum of interest!It is "<<v_max_p3[0] <<endl;
  }else{
    cout << "Found "<<v_max_p3.size()<<" maxima within intervals, not a good pulse shape!"<<endl;
    b_max_p3 = false;
  }
    



  /*
  //POLY4

  cout << "============ poly4 ==============" << endl<<endl;

  TF1 *poly4 = new TF1("poly4","[0] + [1]*x +[2]*x**2 + [3]*x**3 + [4]*x**4 ",-40.,40.);
  poly4->SetLineColor(4);
  ge_pulse->Fit("poly4","M0Q","0");
  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly4->GetChisquare(),poly4->GetNDF());

  double max_poly4 = poly4->GetMaximum(low_limit, high_limit);


  //POLY5
  cout << "================  poly5 =================" << endl<<endl;

  TF1 *poly5 = new TF1("poly5","[0] + [1]*x +[2]*x**2 + [3]*x**3 + [4]*x**4 + [5]*x**5",-40.,40.);
  poly5->SetLineColor(6);
  ge_pulse->Fit("poly5","M");
  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly5->GetChisquare(),poly5->GetNDF()); 

  double max_poly5 = poly5->GetMaximum(low_limit, high_limit);
  */
  //MIchele - End
  /*
  TF1* fu_pulse2 = STTool::halfGauss2;
  ge_pulse->Fit(fu_pulse2,"MQ","0");
  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse2->GetChisquare(),fu_pulse2->GetNDF());
  */

  
  TVirtualFitter *fit_volt = TVirtualFitter::GetFitter();
  Double_t* cov_volt       = fit_volt->GetCovarianceMatrix();
  
  double maxPos = 3-TMath::Sqrt(3);
  double maxFac = TMath::Exp(-maxPos)*(maxPos*maxPos/2.0-maxPos*maxPos*maxPos/6.0);
  
  double max_val = maxFac*fu_pulse->GetParameter(0);
  double max_err = maxFac*fu_pulse->GetParError(0);
  
  double tau_val = fu_pulse->GetParameter(2);
  double tau_err = fu_pulse->GetParError(2);
  
  double t0_val  = fu_pulse->GetParameter(1);
  double t0_err  = fu_pulse->GetParError(1);
  
  // Calculate integrale; factor 1/20 is included such that the maximal value of the pulse and the integral have
  // a similar numerical value
  double int_val = 9.0/2.0*TMath::Exp(-3.0)*fu_pulse->GetParameter(0)*tau_val/20.0;
  double int_err = int_val*TMath::Sqrt(TMath::Power(fu_pulse->GetParError(0)/fu_pulse->GetParameter(0),2)+
                                       TMath::Power(tau_err/tau_val,2)+
                                       2.0*cov_volt[2]/(tau_val*fu_pulse->GetParameter(0))
                                       );
                                       

  //Michele 3/9/2017
  //Estimate the  systematics: # if a polynomial with zeros is found take as sys the semidifference 
  //                             between its integral and the integral obtained from the fitted theoretical model
  //                           # if no polynomial behaved properly set the systematic error to -999
  //Michele - end
  double max_sys = -999;
  if (b_max_p3 == true){
    max_sys = TMath::Abs(max_val - v_max_p3[0])/2.;
  }


  
  mg_pulse->Add(ge_pulse,"PE");
  
  
  
  Info("pulseFit","pulse Height: %8.4f+/-%8.4f",max_val,max_err);
  Info("pulseFit","         tau: %8.4f+/-%8.4f ns",tau_val,tau_err);

  //Michele
  //TO DO: control
  if(b_max_p3 == true){
    Info("pulseFit","   charge eq: %8.4f +/- %8.4f(stat) +/- %8.4f (sys)",max_val,max_err, max_sys);
  }else{

    Info("pulseFit","   charge eq: %8.4f +/- %8.4f(stat) +/- No sys eval.", max_val,max_err);
  }
  //
  // Plot pulse
  TCanvas* c_pul = new TCanvas("c_pul","Pulse Canvas",800,600);
  
  mg_pulse->Draw("A");
  fu_pulse->Draw("Same");
  //fu_pulse2->Draw("Same");
  poly3->Draw("Same");
  //poly4->Draw("Same");
  //poly5->Draw("Same");
  mg_pulse->Draw("");
  mg_pulse->GetXaxis()->SetLimits(-30.0,30.0);
  mg_pulse->GetXaxis()->SetTitle("#delta#it{t} [ns]");
  mg_pulse->SetMaximum( 1.2*max_val);
  mg_pulse->SetMinimum(0.0);
  mg_pulse->GetYaxis()->SetTitle("MPV [ADC value]");
  TLine* line = new TLine(-30.0,0.0,30.0,0.0);
  //line->Draw();
  mg_pulse->GetHistogram()->Draw("axissame");
  
  // 0,1: Maximal signal height
  v_res(0) = max_val;
  v_res(1) = max_err;
  //v_res(1) = max_sys;

  // 2,3: tau (width of pulse)
  v_res(2) = tau_val;
  v_res(3) = tau_err;
  // 4,5: t0 (first root)b
  v_res(4) = t0_val;
  v_res(5) = t0_err;
  // 6,7: covariance (tau,A0) (t0,A0) (t0,tau)
  v_res(8) = cov_volt[1];
  v_res(9) = cov_volt[2];
  v_res(10)= cov_volt[5];
  v_res(11) = max_sys;
  
  
  // Write the data and the plots

  while (!gSystem->Exec(Form("ls %s &> /dev/null",fn_dummy.Data()))) {
    Info("landauFit","Waiting until file is closed");
    gSystem->Sleep(1000);
  }

  gSystem->Exec(Form("touch %s &> /dev/null",fn_dummy.Data()));
  
  f_data = new TFile(fn_data.Data(),"UPDATE");
  
  TString dn_result = TString::Format("%s/%s/%d/%d",
                                      det.Data(),
                                      lay.Data(),
                                      sector,
                                      fill);
  
  // Get right directory
  bool isThere = f_data->GetListOfKeys()->Contains(det.Data());
  if (!isThere) {
    f_data->mkdir(det.Data());
  }
  
  isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s",
                                                   det.Data(),
                                                   lay.Data()));
  if (!isThere) {
    f_data->mkdir(Form("%s/%s",
                       det.Data(),
                       lay.Data()));
  }
  
  isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s/%d",
                                                   det.Data(),
                                                   lay.Data(),
                                                   sector));
  if (!isThere) {
    f_data->mkdir(Form("%s/%s/%d",
                       det.Data(),
                       lay.Data(),
                       sector));
  }
  
  isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s/%d/%d",
                                                   det.Data(),
                                                   lay.Data(),
                                                   sector,
                                                   fill));
  if (!isThere) {
    f_data->mkdir(Form("%s/%s/%d/%d",
                       det.Data(),
                       lay.Data(),
                       sector,
                       fill));
  }
  
  
  f_data->cd(Form("%s/%s/%d/%d",
                  det.Data(),
                  lay.Data(),
                  sector,
                  fill));
  
  if (!notSave) {
    gDirectory->Delete(Form("%s;*",
                            vn_result.Data()));
    
   v_res.Write(Form("%s",
                    vn_result.Data()));
    Info("pulseFit","Result written to %s",fn_data.Data());
   
  }
  
  f_data->Close();
  gSystem->Exec(Form("rm -f %s &> /dev/null",fn_dummy.Data()));
  if (!notSave) {
    // Save picture
    gSystem->Exec(Form("mkdir -p %s",dn_graph.Data()));
    gSystem->Exec(Form("rm %s &> /dev/null",fn_graph_p.Data()));
    c_pul->SaveAs(fn_graph_p.Data());
  }
  
  Printf("=========================================\n");

  return EXIT_SUCCESS;
}