//************************************************
// Author: Federica Lionetto
// Created on: 19/12/2014
//************************************************

/*
Focusing reads all the ROOT files of a given focusing data taking and calculates the z coordinate of the focus.
For each z position, it fits the s-curve with two erf functions, one for the left side and one for the right side, and returns the mu and sigma, that will be used to determine the focus.

The focusing data taking is identified by the following information:
- <sensor>, that is, the type of sensor (Hans410, ...);
- <filename>, that is, the filename excluding the position value and the run type.
- <firstx>, that is, the first position along x;
- <lastx>, that is, the last position along x;
- <stepx>, that is, the step between two subsequent data acquisitions in x (1 step = 5 microns);
- <firstz>, that is, the first position along z;
- <lastz>, that is, the last position along z;
- <stepz>, that is, the step between two subsequent data acquisitions in z (1 step = 5 microns).

Compile with:

make

Run with:

./Focusing [sensor] [filename] [firstx] [lastx] [stepx] [firstz] [lastz] [stepz] [additional folder]

For example:

./Focusing Hans410 ProcessRawData-20141219 0 50 2 -400 -200 20 

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.
*/

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

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

void Focusing(char *sensor, char *filename, const Float_t firstx, const Float_t lastx, const Float_t stepx, const Float_t firstz, const Float_t lastz, const Float_t stepz, char *externalPath=0);

Double_t fLeft(Double_t *x, Double_t *par);
Double_t fRight(Double_t *x, Double_t *par);

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 < 9)
  {
    cout << "**************************************************" << endl;

    cout << "Error! Arguments missing..." << endl;
    cout << "Please use the following format:" << endl;
    cout << "./Focusing [1] [2] [3] [4] [5] [6] [7] [8] [9]" << endl;
    cout << "with:" << endl;
    cout << "[1] = Type of sensor (Hans410, ...);" << endl;
    cout << "[2] = Filename, excluding the position value and the run type;" << endl;
    cout << "[3] = First position along x;" << endl;
    cout << "[4] = Last position along x" << endl;
    cout << "[5] = Step between two subsequent data acquisitions in x (1 step = 5 microns);" << endl;
    cout << "[6] = First position along z;" << endl;
    cout << "[7] = Last position along z" << endl;
    cout << "[8] = Step between two subsequent data acquisitions in z (1 step = 5 microns);" << endl;
    cout << "[9] = Additional folder, optional." << endl;
    cout << "Type ./Focusing --info for more information." << endl;
    
    cout << "**************************************************" << endl;

    return 0;
  }
  else
  {
    cout << "Type of sensor: " << argv[1] << endl;
    cout << "Filename: " << argv[2] << endl;
    cout << "First position along x: " << argv[3] << endl;
    cout << "Last position along x: " << argv[4] << endl;
    cout << "Step between two subsequent data acquisitions in x: " << argv[5] << endl;
    cout << "First position along z: " << argv[6] << endl;
    cout << "Last position along z: " << argv[7] << endl;
    cout << "Step between two subsequent data acquisitions in z: " << argv[8] << endl;
    if (argc == 9)
      Focusing(argv[1],argv[2],atoi(argv[3]),atoi(argv[4]),atoi(argv[5]),atoi(argv[6]),atoi(argv[7]),atoi(argv[8]));
    else if (argc == 10)
      Focusing(argv[1],argv[2],atoi(argv[3]),atoi(argv[4]),atoi(argv[5]),atoi(argv[6]),atoi(argv[7]),atoi(argv[8]),argv[9]);
    else
    {
      cout << "Error! Too many arguments given..." << endl;
      
      return 0;
    }
    
    return 0;  
  }
}

void Focusing(char *sensor, char *filename, const Float_t firstx, const Float_t lastx, const Float_t stepx, const Float_t firstz, const Float_t lastz, const Float_t stepz, char *externalPath)
{
  cout << "**************************************************" << endl;
  cout << "Focusing..." << endl;
  cout << "**************************************************" << endl;

  // Do not comment this line.
  gROOT->ProcessLine("#include <vector>");

  // Do some fanciness to get the directory right.
  string inputDirectory = "/disk/groups/hep/flionett/TestStand/AnalysisResults/"+string(sensor)+"/Alignment";  
  string outputDirectory = "/disk/groups/hep/flionett/TestStand/AnalysisResults/"+string(sensor)+"/Focusing";
  if (externalPath!=0)
    outputDirectory = string(outputDirectory+"/"+externalPath);
  cout << "The input directory is: " << inputDirectory << endl;
  cout << "The output directory is: " << outputDirectory << endl;

  // Create the outputDirectory directory if it does not exist.
  string path_to_make = "mkdir -p "+outputDirectory;
  system(path_to_make.c_str());

  // Necessary to access the input ROOT files.
  ostringstream convertx;
  ostringstream convertz;
  string tempx;
  string tempz;

  cout << "Figures stored in: " << outputDirectory+"/Figures/"+filename << endl;

  // Create a directory named Figures inside the directory named outputDirectory if it does not exist.
  string path_to_make_figures = "mkdir -p "+outputDirectory+"/Figures/"+filename;
  system(path_to_make_figures.c_str());

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

  // Necessary to find the strips hit by the laser.
  int strip;
  string direction;

  char char_input_ROOT[200];
  string input_ROOT;
  string output_ROOT;

  string inputFindStrip;
  string outputFindStrip;

  const Int_t stepsx = (Int_t)((lastx-firstx)/stepx+1);
  Float_t x[stepsx];
  Float_t mean4ADC[stepsx];

  const Int_t stepsz = (Int_t)((lastz-firstz)/stepz+1);
  Float_t z[stepsz];
  Float_t muLeft[stepsz];
  Float_t sigmaLeft[stepsz];
  Float_t muRight[stepsz];
  Float_t sigmaRight[stepsz];



  // Open output ROOT file.
  output_ROOT = outputDirectory+"/Focusing-"+filename+".root";
  TFile *output = TFile::Open(TString(output_ROOT),"RECREATE");

  TGraph *gcheckAlignment4[stepsz];
  TCanvas *ccheckAlignment4[stepsz];

  for (Int_t zcounter=0;zcounter<stepsz;zcounter++) {
    z[zcounter] = firstz+zcounter*stepz;



    convertz.str("");
    convertz << z[zcounter];
    tempz = convertz.str();

    

    for (Int_t xcounter=0;xcounter<stepsx;xcounter++) {
      x[xcounter] = firstx+xcounter*stepx;

      cout << "Position along z: " << z[zcounter] << endl;
      cout << "Position along x: " << x[xcounter] << endl;
      
      // Open input data ROOT files.
      sprintf(char_input_ROOT,"%s/%s-%dx-%dz-las.root",inputDirectory.c_str(),filename,(int)x[xcounter],(int)z[zcounter]);
      input_ROOT = string(char_input_ROOT);
      // Check that the filename provided corresponds to a ROOT file.
      int found = input_ROOT.find(".root");
      if (found==string::npos)  {
	cout << "Error! The filename provided is not associated to a ROOT file." << endl;
	return;
      }
      
      cout << "Open ROOT file #" << xcounter+1 << ": " << input_ROOT << endl;
      TFile *input = TFile::Open(TString(input_ROOT));
      
      inputFindStrip = input_ROOT;
      convertx.str("");
      convertx << x[xcounter];
      tempx = convertx.str();
      outputFindStrip = outputDirectory+"/FindStrip-"+filename+"-"+tempx+"x-"+tempz+"z-las.root";
      FindStrip(inputFindStrip,outputFindStrip,path_to_figures,&strip,&direction);  
      
      cout << "Histograms: " << endl;
      cout << Form("hADCPedSub%d",strip) << endl;
      cout << Form("hADCPedSub%d",strip+NSkip) << endl;
      cout << Form("hADCPedSub%d",strip-NSkip) << endl;
      cout << Form("hADCPedSub%d",strip-2*NSkip) << endl;
      cout << Form("hADCPedSub%d",strip+2*NSkip) << endl;
      
      // Four adjacent strips.
      TH1D *hist1 = (TH1D *)input->Get(Form("hADCPedSub%d",strip));
      TH1D *hist2 = (TH1D *)input->Get(Form("hADCPedSub%d",strip+NSkip));
      TH1D *hist3 = (TH1D *)input->Get(Form("hADCPedSub%d",strip-NSkip));
      TH1D *hist4;
      if (direction == "left")
	hist4 = (TH1D *)input->Get(Form("hADCPedSub%d",strip-2*NSkip));
      else if (direction == "right")
	hist4 = (TH1D *)input->Get(Form("hADCPedSub%d",strip+2*NSkip));
      mean4ADC[xcounter] = hist1->GetMean()+hist2->GetMean()+hist3->GetMean()+hist4->GetMean();
      
      // Close input data ROOT files.
      input->Close();
    }

    output->cd();
    
    for (Int_t xcounter=0;xcounter<stepsx;xcounter++)
      x[xcounter] = x[xcounter]*5.; // One step corresponds to 5 microns.

    // Four adjacent strips.
    gcheckAlignment4[zcounter] = new TGraph(stepsx,x,mean4ADC);
    InitGraph(gcheckAlignment4[zcounter],Form("Focusing - 4 adjacent strips - z = %d",(int)z[zcounter]),"x (#mum)","ADC counts");

    // Fit.

    // Right side.
    Float_t minRight = 40.;
    Float_t maxRight = 140.;
    TF1 *fitRight = new TF1("fitRight",fRight,minRight,maxRight,4);
    fitRight->SetParameter(0,100.);
    fitRight->SetParameter(1,80.);
    fitRight->SetParameter(2,10.);
    fitRight->SetParameter(3,200.);
    fitRight->SetParLimits(1,60.,100.);
    fitRight->SetParLimits(2,0.,20.);
    gcheckAlignment4[zcounter]->Fit(fitRight,"BR");

    muRight[zcounter] = fitRight->GetParameter(1);
    sigmaRight[zcounter] = fitRight->GetParameter(2);
    cout << "Right side, z = " << zcounter << ", mu = " << muRight[zcounter] << ", sigma = " << sigmaRight[zcounter] << endl;


    // Left side.
    Float_t minLeft = 140.;
    Float_t maxLeft = 240.;
    TF1 *fitLeft = new TF1("fitLeft",fLeft,minLeft,maxLeft,4);
    fitLeft->SetParameter(0,100.);
    fitLeft->SetParameter(1,210.);
    fitLeft->SetParameter(2,10.);
    fitLeft->SetParameter(3,200.);
    fitLeft->SetParLimits(1,200.,240.);
    fitLeft->SetParLimits(2,0.,20.);
    gcheckAlignment4[zcounter]->Fit(fitLeft,"BR");

    muLeft[zcounter] = fitLeft->GetParameter(1);
    sigmaLeft[zcounter] = fitLeft->GetParameter(2);
    cout << "Left side, z = " << zcounter << ", mu = " << muLeft[zcounter] << ", sigma = " << sigmaLeft[zcounter] << endl;



    ccheckAlignment4[zcounter] = new TCanvas(Form("ccheckAlignment4-%s-%dz",filename,(int)z[zcounter]),"",400,300);
    DrawGraphFunc2(ccheckAlignment4[zcounter],gcheckAlignment4[zcounter],fitRight,fitLeft,"AP",path_to_figures);



  }



  for (Int_t zcounter=0;zcounter<stepsz;zcounter++)
    z[zcounter] = z[zcounter]*5.; // One step corresponds to 5 microns.

  // Mu.
  TGraph *gmuLeft = new TGraph(stepsz,z,muLeft);
  InitGraph(gmuLeft,"Focusing - left side","z (#mum)","#mu_{left}");
  TCanvas *cmuLeft = new TCanvas("cmuLeft","",400,300);
  DrawGraph(cmuLeft,gmuLeft,"AP",path_to_figures);

  TGraph *gmuRight = new TGraph(stepsz,z,muRight);
  InitGraph(gmuRight,"Focusing - right side","z (#mum)","#mu_{right}");
  TCanvas *cmuRight = new TCanvas("cmuRight","",400,300);
  DrawGraph(cmuRight,gmuRight,"AP",path_to_figures);

  // Sigma.
  TGraph *gsigmaLeft = new TGraph(stepsz,z,sigmaLeft);
  InitGraph(gsigmaLeft,"Focusing - left side","z (#mum)","#sigma_{left}");
  TCanvas *csigmaLeft = new TCanvas("csigmaLeft","",400,300);
  DrawGraph(csigmaLeft,gsigmaLeft,"AP",path_to_figures);

  TGraph *gsigmaRight = new TGraph(stepsz,z,sigmaRight);
  InitGraph(gsigmaRight,"Focusing - right side","z (#mum)","#sigma_{right}");
  TCanvas *csigmaRight = new TCanvas("csigmaRight","",400,300);
  DrawGraph(csigmaRight,gsigmaRight,"AP",path_to_figures);



  // Write output ROOT file.
  output->Write();

  // Close output ROOT file.
  output->Close();

  return;
}

// fLeft has four parameters.
Double_t fLeft(Double_t *x, Double_t *par)
{
  return par[0]*(1.-erf((x[0]-par[1])/(sqrt(2.)*par[2])))+par[3];
}

// fRight has four parameters.
Double_t fRight(Double_t *x, Double_t *par)
{
return par[0]*erf((x[0]-par[1])/(sqrt(2.)*par[2]))+par[3];
}