//************************************************
// Author: Federica Lionetto
// Created on: 03/11/2014
//************************************************

/*
ProcessRawData reads a ROOT file (containing the raw data to process) and a text file (containing the pedestals) and creates a ROOT file with the following information:
- all branches in the Header tree;
- all branches in the EventInfo tree, one per event;
- string InputFilename, the complete path to the input file (ROOT file);
- string PedFilename, the complete path to the pedestal file (text file);
- std::vector<float> Ped, the pedestals calculated by ComputePedestals;
- std::vector<float> ADCProcessed, one per event (ADC after pedestal subtraction);
- distribution of the ADC counts of each Beetle (TH1F *hADCBeetle1 and TH1F *hADCBeetle2, before pedestal subtraction, and TH1F *hADCPedSubBeetle1 and TH1F *hADCPedSubBeetle2, after pedestal subtraction);
- distribution of the ADC counts of each Beetle channel (TH1F *hADC[N], before pedestal subtraction, and TH1F *hADCPedSub[N], after pedestal subtraction);
- ADC counts as a function of the Beetle channel (event-by-event, displayed only for one particular event, TH1F *hADCvsStripBeetle1 and TH1F *hADCvsStripBeetle2, before pedestal subtraction, and TH1F *hADCPedSubvsStripBeetle1 and TH1F *hADCPedSubvsStripBeetle2, after pedestal subtraction).

Compile with:

make

Run with:

./ProcessRawData [input ROOT file] [input text file (pedestals)] [additional folder]

where:
- [input ROOT file] is the complete path, including the folder and the filename, of the ROOT file one wants to process;
- [input text file (pedestals)] is the complete path, including the folder and the filename, of the text file one wants to get the pedestals from;
- [additional folder] is the optional additional folder where the output will be saved.
*/

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

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

void ProcessRawData(char *filename, char *filenamePed, char *externalPath=0);

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

  if ((argc ==2) && (string(argv[1]) == "--info"))
  {
    cout << "**************************************************" << endl;

    cout << "Some comments." << endl;
    
    cout << "**************************************************" << endl;

    return 0;
  }
  else if (argc < 3)
  {
    cout << "**************************************************" << endl;
    
    cout << "Error! Input file or pedestal file missing..." << endl;
    cout << "Please use the following format:" << endl;
    cout << "./ProcessRawData [1] [2] [3]" << endl;
    cout << "with:" << endl;
    cout << "[1] = Input ROOT file, complete path;" << endl;
    cout << "[2] = Input text file (pedestals), complete path;" << endl;
    cout << "[3] = Additional folder, optional." << endl;
    cout << "Type ./ProcessRawData --info for more information." << endl;

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

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

void ProcessRawData(char *filename, char *filenamePed, char *externalPath)
{
  cout << "**************************************************" << endl;
  cout << "Processing raw data..." << endl;
  cout << "**************************************************" << endl;
 
  // Do not comment this line.
  gROOT->ProcessLine("#include <vector>");

  // Do some fanciness to get the directory right.
  // string analysis = "/home/hep/flionett/TestStand/AnalysisResults";
  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"+"/ProcessRawData";
  system(path_to_make_figures.c_str());
  path_to_make_figures = "mkdir "+analysis+"/Figures"+"/ProcessRawData"+"/"+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;

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

  TBranch *b_NEvents = 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->GetEntry();

    cout << "Number of events: " << NEvents << 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);
  }

  // Open input pedestal text file.
  FILE *inputText = fopen(TString(filenamePed),"r");

  // Read pedestals.
  float ped[N];
  
  for (int iChannel=0; iChannel<N; ++iChannel) {
    fscanf(inputText,"%*d%*c%f%*c%*f%*c%*f%*c%*f%*c",&(ped[iChannel]));
    // cout << iChannel << ", " << ped[iChannel] << endl;
  }
  
  // Close input pedestal text file.
  fclose(inputText);

  // Open output data ROOT file.
  TFile *output = TFile::Open(analysis+"/ProcessRawData-"+TString(filename_as_string_no_path),"RECREATE");

  // Copy the Header tree to the output data ROOT file and add some more information.
  string InputFilename;
  string PedFilename;
  std::vector<float> Ped;

  InputFilename = string(filename);
  PedFilename = string(filenamePed);

  for (int iChannel=0; iChannel<N; ++iChannel)
  {
    Ped.push_back(ped[iChannel]);
  }
  TTree *HeaderProcessed = Header->CloneTree(0);
  HeaderProcessed->Branch("InputFilename",&InputFilename);
  HeaderProcessed->Branch("PedFilename",&PedFilename);
  HeaderProcessed->Branch("Ped",&Ped);
  HeaderProcessed->Fill();

  // Show the distributions for one particular event.
  int sample = 970; // Modify it.
  
  // Copy the EventInfo tree to the output data ROOT file and add some more information.
  std::vector<float> ADCProcessed; // Raw ADC - pedestal.

  TTree *EventInfoProcessed = EventInfo->CloneTree(0);
  EventInfoProcessed->Branch("ADCProcessed",&ADCProcessed);
  
  // Distribution of the ADC counts of each Beetle (TH1F *hADCBeetle1 and TH1F *hADCBeetle2, before pedestal subtraction, and TH1F *hADCPedSubBeetle1 and TH1F *hADCPedSubBeetle2, after pedestal subtraction);

  TH1F *hADCBeetle1 = new TH1F("hADCBeetle1","",400,300,700); // Modify it.
  TCanvas *cADCBeetle1 = new TCanvas("cADCBeetle1","",400,300);
  InitHist(hADCBeetle1,"Raw ADC counts on Beetle 1","ADC counts","");

  TH1F *hADCPedSubBeetle1 = new TH1F("hADCPedSubBeetle1","",200,-100,100); // Modify it.
  TCanvas *cADCPedSubBeetle1 = new TCanvas("cADCPedSubBeetle1","",400,300);
  InitHist(hADCPedSubBeetle1,"Pedestal subtracted ADC counts on Beetle 1","ADC counts","");

  TH1F *hADCBeetle2 = new TH1F("hADCBeetle2","",400,300,700); // Modify it.
  TCanvas *cADCBeetle2 = new TCanvas("cADCBeetle2","",400,300);
  InitHist(hADCBeetle2,"Raw ADC counts on Beetle 2","ADC counts","");

  TH1F *hADCPedSubBeetle2 = new TH1F("hADCPedSubBeetle2","",200,-100,100); // Modify it.
  TCanvas *cADCPedSubBeetle2 = new TCanvas("cADCPedSubBeetle2","",400,300);
  InitHist(hADCPedSubBeetle2,"Pedestal subtracted ADC counts on Beetle 2","ADC counts","");

  // Distribution of the ADC counts of each Beetle channel (TH1F *hADC[N], before pedestal subtraction, and TH1F *hADCPedSub[N], after pedestal subtraction);

  TH1F *hADC[N];
  TH1F *hADCPedSub[N];

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

    hADCPedSub[iChannel] = new TH1F(Form("hADCPedSub%d",iChannel),"",1024,-512,512);
    InitHist(hADCPedSub[iChannel],Form("Pedestal subtracted ADC counts on Beetle channel %d",iChannel),"ADC counts","");
  }

  // ADC counts as a function of the Beetle channel (event-by-event, displayed only for one particular event, TH1F *hADCvsStripBeetle1 and TH1F *hADCvsStripBeetle2, before pedestal subtraction, and TH1F *hADCPedSubvsStripBeetle1 and TH1F *hADCPedSubvsStripBeetle2, after pedestal subtraction).

  TH1F *hADCvsStripBeetle1 = new TH1F("hADCvsStripBeetle1","",128,0,128);
  TCanvas *cADCvsStripBeetle1 = new TCanvas("cADCvsStripBeetle1","",400,300);

  TH1F *hADCvsStripBeetle2 = new TH1F("hADCvsStripBeetle2","",128,128,256);
  TCanvas *cADCvsStripBeetle2 = new TCanvas("cADCvsStripBeetle2","",400,300);
 
  InitHist(hADCvsStripBeetle1,Form("Event %d - Beetle 1, raw",sample),"Beetle channel","ADC counts");
 
  InitHist(hADCvsStripBeetle2,Form("Event %d - Beetle 2, raw",sample),"Beetle channel","ADC counts");

  TH1F *hADCPedSubvsStripBeetle1 = new TH1F("hADCPedSubvsStripBeetle1","",128,0,128);
  TCanvas *cADCPedSubvsStripBeetle1 = new TCanvas("cADCPedSubvsStripBeetle1","",400,300);

  TH1F *hADCPedSubvsStripBeetle2 = new TH1F("hADCPedSubvsStripBeetle2","",128,128,256);
  TCanvas *cADCPedSubvsStripBeetle2 = new TCanvas("cADCPedSubvsStripBeetle2","",400,300);
 
  InitHist(hADCPedSubvsStripBeetle1,Form("Event %d - Beetle 1, ped sub",sample),"Beetle channel","ADC counts");
 
  InitHist(hADCPedSubvsStripBeetle2,Form("Event %d - Beetle 2, ped sub",sample),"Beetle channel","ADC counts");

  // Style and range for Beetle 1 plots.
  hADCvsStripBeetle1->SetFillColor(kAzure);
  hADCPedSubvsStripBeetle1->SetFillColor(kAzure);
  hADCvsStripBeetle1->SetMinimum(0);
  hADCvsStripBeetle1->SetMaximum(1200);
  hADCPedSubvsStripBeetle1->SetMinimum(-150);
  hADCPedSubvsStripBeetle1->SetMaximum(900);

  // Style and range for Beetle 2 plots.
  hADCvsStripBeetle2->SetFillColor(kRed);
  hADCPedSubvsStripBeetle2->SetFillColor(kRed);
  hADCvsStripBeetle2->SetMinimum(0);
  hADCvsStripBeetle2->SetMaximum(1200);
  hADCPedSubvsStripBeetle2->SetMinimum(-150);
  hADCPedSubvsStripBeetle2->SetMaximum(900);

  // Fill histograms.
  for (int iEvent=0;iEvent<NEvents;++iEvent) {
    if (iEvent%1000==0)
      cout << "Processing Event " << iEvent << endl;

    // Initialize histograms.
    if (iEvent>0) 
    {
      ADCProcessed.clear();
    }    
    
    EventInfo->GetEntry(iEvent);
    
    // Fill histograms.
    for (int iChannel=0;iChannel<N;++iChannel) { 
      if (ADC->at(iChannel) > 1024.)
	cout << "Error! ADC counts of Beetle channel " << iChannel << ": " << ADC->at(iChannel) << endl; 
      
      ADCProcessed.push_back((float)ADC->at(iChannel)-ped[iChannel]);
      
      if ((iChannel>=0) && (iChannel<NBeetle))
      {
	hADCBeetle1->Fill((float)ADC->at(iChannel));
	hADCPedSubBeetle1->Fill((float)ADC->at(iChannel)-ped[iChannel]);
	if (iEvent == sample) {
	  hADCvsStripBeetle1->Fill(iChannel,(float)ADC->at(iChannel));
	  hADCPedSubvsStripBeetle1->Fill(iChannel,(float)ADC->at(iChannel)-ped[iChannel]);
	}
      }
      else if ((iChannel>=NBeetle) && (iChannel<N))
      {
	hADCBeetle2->Fill((float)ADC->at(iChannel));
	hADCPedSubBeetle2->Fill((float)ADC->at(iChannel)-ped[iChannel]);
	if (iEvent == sample) {
	  hADCvsStripBeetle2->Fill(iChannel,(float)ADC->at(iChannel));
	  hADCPedSubvsStripBeetle2->Fill(iChannel,(float)ADC->at(iChannel)-ped[iChannel]);
	}
      }
      
      hADC[iChannel]->Fill((float)ADC->at(iChannel));
      hADCPedSub[iChannel]->Fill((float)ADC->at(iChannel)-ped[iChannel]);
    }
    
    EventInfoProcessed->Fill();
  }   

  DrawHist(cADCBeetle1,hADCBeetle1,"",path_to_figures);
  DrawHist(cADCPedSubBeetle1,hADCPedSubBeetle1,"",path_to_figures);

  DrawHist(cADCBeetle2,hADCBeetle2,"",path_to_figures);
  DrawHist(cADCPedSubBeetle2,hADCPedSubBeetle2,"",path_to_figures);

  DrawHist(cADCvsStripBeetle1,hADCvsStripBeetle1,"",path_to_figures);      
  DrawHist(cADCPedSubvsStripBeetle1,hADCPedSubvsStripBeetle1,"",path_to_figures);      
      
  DrawHist(cADCvsStripBeetle2,hADCvsStripBeetle2,"",path_to_figures);      
  DrawHist(cADCPedSubvsStripBeetle2,hADCPedSubvsStripBeetle2,"",path_to_figures);      

  // Write output ROOT file.
  output->Write();
  
  // Close input and output ROOT files.
  input->Close();
  output->Close();
  
  return;
}