//************************************************
// Author: Federica Lionetto
// Created on: 07/31/2014
//************************************************

/*
ComputePedestals reads a ROOT file and creates a ROOT file with the following information:
- distribution of the ADC counts of each Beetle;
- distribution of the ADC counts of each Beetle channel (hADCx, with x = {0,...,255});
- mean of the above distribution as a function of the Beetle channel (the pedestal);
- RMS of the above distribution as a function of the Beetle channel (the noise);
- distribution of the pedestals of each Beetle;
- distribution of the noise of each Beetle.

It also creates a text file with five columns:
- Beetle channel, from 0 to 255;
- pedestal;
- pedestal uncertainty;
- noise;
- noise uncertainty.

Compile with:

make

Run with:

./ComputePedestals [input ROOT file]

where [input ROOT file] is the complete path, including the folder and the filename, of the ROOT file one wants to process.
For example:

./ComputePedestals /.../run_xxxxxx_xxxx_xxxx.root

A folder named AnalysisResults will be created in a fixed location and a ROOT file will be saved in there. A folder named Figures will be created inside the folder named AnalysisResults, with some monitoring plots.

If needed, an additional folder can be specified, that will be created inside the folder named AnalysisResults.
In this case, run with:

./ComputePedestals [input ROOT file] [additional folder]
*/

#include "../Tools/Lib.C"
#include "../Tools/lhcbStyle.C"
#include "../Tools/Style.C"

#include "../Tools/Par.C"

void ComputePedestals(char *filename, char *externalPath = 0);

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

  if (argc < 2)
  {
    cout << "**************************************************" << endl;

    cout << "Error! Input file missing..." << endl;
    cout << "Please use the following format:" << endl;
    cout << "./ComputePedestals [1] [2]" << endl;
    cout << "with:" << endl;
    cout << "[1] = Input ROOT file, complete path;" << endl;
    cout << "[2] = Additional folder, optional." << endl;
    cout << "Type ./ComputePedestals --info for more information." << endl;
    
    cout << "**************************************************" << endl;

    return 0;
  }
  else if (string(argv[1]) == "--info")
  {
    cout << "**************************************************" << endl;

    cout << "ComputePedestals reads a ROOT file and creates a ROOT file with the following information:" << endl;
    cout << "- distribution of the ADC counts of each Beetle;" << endl;
    cout << "- distribution of the ADC counts of each Beetle channel (hADCx, with x = {0,...,255});" << endl;
    cout << "- mean of the above distribution as a function of the Beetle channel (the pedestal);" << endl;
    cout << "- RMS of the above distribution as a function of the Beetle channel (the noise);" << endl;
    cout << "- distribution of the pedestals of each Beetle;" << endl;
    cout << "- distribution of the noise of each Beetle." << endl;
    cout << "It also creates a text file with five columns:" << endl;
    cout << "- Beetle channel, from 0 to 255;" << endl;
    cout << "- pedestal;" << endl;
    cout << "- pedestal uncertainty;" << endl;
    cout << "- noise;" << endl;
    cout << "- noise uncertainty." << endl;
 
    cout << "**************************************************" << endl;

    cout << "Compile with:" << endl;
    cout << "make" << endl;
    cout << "Run with:" << endl;
    cout << "./ComputePedestals [input ROOT file]" << endl;
    cout << "where [input ROOT file] is the complete path, including the folder and the filename, of the file one wants to process." << endl;
    cout << "For example:" << endl;
    cout << "./ComputePedestals /.../run_xxxxxx_xxxx_xxxx.root" << endl;
    
    cout << "**************************************************" << endl;

    cout << "A folder named AnalysisResults will be created in a fixed location and a ROOT file will be saved in there. A folder named Figures will be created inside the folder named AnalysisResults, with some monitoring plots." << endl;

    cout << "If needed, an additional folder can be specified, that will be created inside the folder named AnalysisResults." << endl;
    cout << "In this case, run with:" << endl;
    cout << "./ComputePedestals [input ROOT file] [additional folder]" << endl;

    cout << "**************************************************" << endl;

    return 0;
  }
  else
  {
    cout << "File to process: " << argv[1] << endl;
    if (argc == 2)
      ComputePedestals(argv[1]);
    else if (argc == 3)
      ComputePedestals(argv[1],argv[2]);
    else
    {
      cout << "Error! Too many arguments given..." << endl;
      
      return 0;
    }
    
    return 0;  
  }
}



void ComputePedestals(char *filename, char *externalPath) 
{
  cout << "**************************************************" << endl;
  cout << "Computing pedestals..." << endl;
  cout << "**************************************************" << endl;
 
  // Do not comment this line.
  gROOT->ProcessLine("#include <vector>");

  // Do some fanciness to get the directory right.
  string analysis = "/disk/groups/hep/flionett/TestStand/AnalysisResults";
  if (externalPath!=0)
    analysis = string(analysis+"/"+externalPath);

  string filename_as_string = string(filename);
  // Check that the filename provided corresponds to a ROOT file.
  int found = filename_as_string.find(".root");
  if (found==string::npos)  {
    cout << "Error! The filename provided is not associated to a ROOT file." << endl;
    return;
  }

  string filename_as_string_no_path = filename_as_string.substr(filename_as_string.find_last_of('/')+1);
  string filename_as_string_no_path_no_extension = filename_as_string_no_path.substr(0,filename_as_string_no_path.length()-5);

  // Create a folder named AnalysisResults. Do not worry, nothing bad is going to happen if the folder already exists.
  string path_to_make = "mkdir -p "+analysis;
  system(path_to_make.c_str());
  cout << "Setting output to: " << path_to_make << endl;

  // Create a folder named Figures inside the folder named AnalysisResults.
  string path_to_make_figures = "mkdir "+analysis+"/Figures";
  system(path_to_make_figures.c_str());
  path_to_make_figures = "mkdir "+analysis+"/Figures"+"/ComputePedestals";
  system(path_to_make_figures.c_str());
  path_to_make_figures = "mkdir "+analysis+"/Figures"+"/ComputePedestals"+"/"+filename_as_string_no_path_no_extension;
  system(path_to_make_figures.c_str());
  cout << "Setting figures to: " << path_to_make_figures << endl;  

  TString path_to_figures = (string(path_to_make_figures)).substr((string(path_to_make_figures)).find_last_of(' ')+1);

  // Open input data ROOT file.
  TFile *input = TFile::Open(TString(filename));

  int NEvents;
  int RunType;

  std::vector<unsigned short> *ADC = 0;

  TBranch *b_NEvents = 0;
  TBranch *b_RunType = 0;
  
  TBranch *b_ADC = 0;

  // Get info from trees.
  TTree *Header = (TTree *)input->Get("Header");
  int EventsHeader = Header->GetEntries();
  if (EventsHeader != 1)
  {
    cout << "Error! The ROOT file has been corrupted..." << endl;

    return ;
  }
  else 
  {
    Header->SetBranchAddress("NEvents",&NEvents,&b_NEvents);
    Header->SetBranchAddress("RunType",&RunType,&b_RunType);

    Header->GetEntry();

    cout << "Number of events: " << NEvents << endl;
    cout << "RunType: " << RunType << endl;
    if (RunType != 5)
      cout << "Warning! Calculating pedestals on a non-pedestal run..." << endl;
  }
    
  // Get info from trees.
  TTree *EventInfo = (TTree *)input->Get("EventInfo");
  int EventsEventInfo = EventInfo->GetEntries();
  if (EventsEventInfo != NEvents)
  {
    cout << "Error! The ROOT file has been corrupted..." << endl;

    return;
  }
  else 
  {
    EventInfo->SetBranchAddress("ADC",&ADC,&b_ADC);
  }

  TFile *output = TFile::Open(analysis+"/ComputePedestals-"+TString(filename_as_string_no_path),"RECREATE");

  // Open output pedestal text file.
  FILE *output_text = fopen(analysis+"/ComputePedestals-"+TString(filename_as_string_no_path_no_extension)+".dat","w");

  TH1F *hADCBeetle1 = new TH1F("hADCBeetle1","",200,400,600);
  TCanvas *cADCBeetle1 = new TCanvas("cADCBeetle1","",400,300);
  InitHist(hADCBeetle1,"Raw ADC counts on Beetle 1","ADC counts","");

  TH1F *hADCBeetle2 = new TH1F("hADCBeetle2","",200,400,600);
  TCanvas *cADCBeetle2 = new TCanvas("cADCBeetle2","",400,300);
  InitHist(hADCBeetle2,"Raw ADC counts on Beetle 2","ADC counts","");
  
  // Fill histogram.
  for (int iEvent=0;iEvent<NEvents;++iEvent) {
    EventInfo->GetEntry(iEvent);
    for (int ch=0; ch<NBeetle; ++ch) 
    {
      hADCBeetle1->Fill((float)ADC->at(ch));
    }    
    for (int ch=NBeetle; ch<N; ++ch)
    {
      hADCBeetle2->Fill((float)ADC->at(ch));
    } 
  }

  // Change the margins.
  cADCBeetle1->SetTopMargin(0.1);
  cADCBeetle2->SetTopMargin(0.1);

  DrawHist(cADCBeetle1,hADCBeetle1,"",path_to_figures);
  DrawHist(cADCBeetle2,hADCBeetle2,"",path_to_figures);

  // Remove potential outliers. Use the mean and the RMS of the hADCTemp histograms and exclude everything that is outside <factor> times the RMS.
  float mean[N];
  float rms[N];
  
  float factor = 3.; // Entries further than factor*rms from the mean are excluded from the pedestal calculation.

  // Fill histograms, necessary to remove potential outliers.
  TH1F *hADCTemp[N];
  for (int ch=0;ch<N;++ch) {
    hADCTemp[ch] = new TH1F(Form("hADCTemp%d",ch),"",1024,0,1024);   
    InitHist(hADCTemp[ch],Form("Raw ADC counts on Beetle channel %d - temporary",ch),"ADC counts","");

    // Fill histograms.
    for (int iEvent=0;iEvent<NEvents;++iEvent) {
      EventInfo->GetEntry(iEvent);
      hADCTemp[ch]->Fill((float)ADC->at(ch));
      if (ADC->at(ch) > 1024.)
        cout << "Error! ADC counts of Beetle channel " << ch << ": " << ADC->at(ch) << endl;
    }

    TF1 *fGaussian = new TF1("fGaussian","gaus",400,600);
    hADCTemp[ch]->Fit("fGaussian","R");
    mean[ch] = fGaussian->GetParameter(1);
    rms[ch] = fGaussian->GetParameter(2);
  }

  // Fill histograms.
  TH1F *hADC[N];
  TCanvas *cADC[N];

  float pedestal[N];
  float upedestal[N]; // Uncertainty.
  float noise[N];
  float unoise[N]; // Uncertainty.

  for (int ch=0;ch<N;++ch) {
    hADC[ch] = new TH1F(Form("hADC%d",ch),"",1024,0,1024);
    cADC[ch] = new TCanvas(Form("cADC%d",ch),"",400,300);
    
    InitHist(hADC[ch],Form("Raw ADC counts on Beetle channel %d",ch),"ADC counts","");

    // Fill histograms.
    for (int iEvent=0;iEvent<NEvents;++iEvent) {
      EventInfo->GetEntry(iEvent);
      if ((abs((float)ADC->at(ch)-mean[ch]))<(factor*rms[ch]))
        hADC[ch]->Fill((float)ADC->at(ch));
      if (ADC->at(ch) > 1024.)
        cout << "Error! ADC counts of Beetle channel " << ch << ": " << ADC->at(ch) << endl;
    }
    
    // Get mean and RMS (that is, pedestal and noise).
    pedestal[ch] = hADC[ch]->GetMean();
    upedestal[ch] = hADC[ch]->GetMeanError();
    noise[ch] = hADC[ch]->GetRMS();
    unoise[ch] = hADC[ch]->GetRMSError();
  }  

  // Mean and RMS as a function of the strip number.
  float strip[N];
  for (int ch=0;ch<N;++ch)
    strip[ch] = (float)ch;

  TH1F *hpedestalvsstrip = new TH1F("hpedestalvsstrip","",256,0,256);
  TH1F *hnoisevsstrip = new TH1F("hnoisevsstrip","",256,0,256);
  
  InitHist(hpedestalvsstrip,"Pedestal as a function of Beetle channel","Beetle channel","ADC counts");
  InitHist(hnoisevsstrip,"Noise as a function of Beetle channel","Beetle channel","ADC counts");

  // Fill histograms.
  for (int ch=0;ch<N;++ch)
  {
    hpedestalvsstrip->Fill((float)ch,pedestal[ch]);
    hnoisevsstrip->Fill((float)ch,noise[ch]);
  }

  // Define range.
  hpedestalvsstrip->SetMinimum(400.);
  hpedestalvsstrip->SetMaximum(600.);
  hnoisevsstrip->SetMinimum(0.);
  hnoisevsstrip->SetMaximum(30.);

  TCanvas *cpedestalvsstrip = new TCanvas("cpedestalvsstrip","",400,300);
  TCanvas *cnoisevsstrip = new TCanvas("cnoisevsstrip","",400,300);

  // Style.
  hpedestalvsstrip->SetFillColor(kAzure);
  hnoisevsstrip->SetFillColor(kRed);
 
  // Stats.
  hpedestalvsstrip->SetStats(0);
  hnoisevsstrip->SetStats(0);

  DrawHist(cpedestalvsstrip,hpedestalvsstrip,"",path_to_figures);
  DrawHist(cnoisevsstrip,hnoisevsstrip,"",path_to_figures);
  
  // Distributions of mean and RMS.
  TH1F *hpedestalBeetle1 = new TH1F("hpedestalBeetle1","",100,400,600);
  TH1F *hnoiseBeetle1 = new TH1F("hnoiseBeetle1","",30,0,30);

  TH1F *hpedestalBeetle2 = new TH1F("hpedestalBeetle2","",100,400,600);
  TH1F *hnoiseBeetle2 = new TH1F("hnoiseBeetle2","",30,0,30);

  InitHist(hpedestalBeetle1,"Pedestals","ADC counts","");
  InitHist(hnoiseBeetle1,"Noise","ADC counts","");

  InitHist(hpedestalBeetle2,"Pedestals","ADC counts","");
  InitHist(hnoiseBeetle2,"Noise","ADC counts","");

  for (int ch=0;ch<NBeetle;++ch) {
    hpedestalBeetle1->Fill(pedestal[ch]);
    hnoiseBeetle1->Fill(noise[ch]);
  }

  for (int ch=NBeetle;ch<N;++ch) {
    hpedestalBeetle2->Fill(pedestal[ch]);
    hnoiseBeetle2->Fill(noise[ch]);
  }

  TCanvas *cpedestalBeetle1 = new TCanvas("cpedestalBeetle1","",400,300);
  TCanvas *cnoiseBeetle1 = new TCanvas("cnoiseBeetle1","",400,300);

  TCanvas *cpedestalBeetle2 = new TCanvas("cpedestalBeetle2","",400,300);
  TCanvas *cnoiseBeetle2 = new TCanvas("cnoiseBeetle2","",400,300);

  DrawHist(cpedestalBeetle1,hpedestalBeetle1,"",path_to_figures);
  DrawHist(cnoiseBeetle1,hnoiseBeetle1,"",path_to_figures);
  
  DrawHist(cpedestalBeetle2,hpedestalBeetle2,"",path_to_figures);
  DrawHist(cnoiseBeetle2,hnoiseBeetle2,"",path_to_figures);

  // Write output pedestal text file.
  for (int ch=0;ch<N;++ch) 
  {
    fprintf(output_text,"%d %.2f %.2f %.2f %.2f\n",ch,pedestal[ch],upedestal[ch],noise[ch],unoise[ch]);
  }
  
  // Write output ROOT file.
  output->Write();
  
  // Close input and output ROOT files.
  input->Close();      
  output->Close();

  // Close output text file.
  fclose(output_text);
  
  return;
}