//************************************************
// 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 calculate the z coordinate of the focus. The starting values and the limits of the parameters of the two erf functions are set according to the filename of the measurement, by using the "assignParInfoSCurves" procedure.
The mu and sigma as a function of the z position are fitted by the functions
mu = p0+p1*z and sigma = q0*|z-q1|+q2, respectively, and the values returned by the fit are printed out on the screen.
The starting values and the limits of the parameters of the two sigma functions are set according to the filename of the measurement, by using the "assignParInfoSigma" procedure. 
Knowing p1, the angle alpha = arctan(p1) between the laser and the normal to the surface of the sensor is also determined. Its uncertainty is calculated from its derivative, which is a good estimate for the moment, but not fully correct.
The fit results are written to a tree called tFitResults.

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/SCurve.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);

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

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

    cout << "Focusing reads all the ROOT files of a given focusing data taking and calculates the z coordinate of the focus." << endl;
    cout << "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 calculate the z coordinate of the focus. The starting values and the limits of the parameters of the two erf functions are set according to the filename of the measurement, by using the 'assignParInfoSCurves' procedure." << endl;
    cout << "The mu and sigma as a function of the z position are fitted by the functions mu = p0+p1*z and sigma = q0*|z-q1|+q2, respectively, and the values returned by the fit are printed out on the screen." << endl;
    cout << "The starting values and the limits of the parameters of the two sigma functions are set according to the filename of the measurement, by using the 'assignParInfoSigma' procedure." << endl; 
    cout << "Knowing p1, the angle alpha = arctan(p1) between the laser and the normal to the surface of the sensor is also determined. Its uncertainty is calculated from its derivative, which is a good estimate for the moment, but not fully correct." << endl;
    cout << "The fit results are written to a tree called tFitResults." << endl;

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

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

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

    cout << "Compile with:" << endl;

    cout << "make" << endl;

    cout << "Run with:" << endl;

    cout << "./Focusing [sensor] [filename] [firstx] [lastx] [stepx] [firstz] [lastz] [stepz] [additional folder]" << endl;

    cout << "For example:" << endl;

    cout << "./Focusing Hans410 ProcessRawData-20141219 0 50 2 -400 -200 20" << 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 << "**************************************************" << 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;

  Float_t noise[N];

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

  const Int_t stepsz = (Int_t)((lastz-firstz)/stepz+1);
  Float_t z[stepsz];
  Float_t uz[stepsz];

  for (int i=0;i<stepsz;i++)
    uz[i] = 0.;

  // Parameters of the fitting functions.
  Float_t aLeft[stepsz];
  Float_t bLeft[stepsz];
  Float_t muLeft[stepsz];
  Float_t sigmaLeft[stepsz];

  Float_t aRight[stepsz];
  Float_t bRight[stepsz];
  Float_t muRight[stepsz];
  Float_t sigmaRight[stepsz];

  // Uncertainties on the parameters of the fitting functions.
  // u = uncertainty.
  Float_t uaLeft[stepsz];
  Float_t ubLeft[stepsz];
  Float_t umuLeft[stepsz];
  Float_t usigmaLeft[stepsz];

  Float_t uaRight[stepsz];
  Float_t ubRight[stepsz];
  Float_t umuRight[stepsz];
  Float_t usigmaRight[stepsz];

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

  TGraph *gcheckAlignmentSignalOverNoise4[stepsz];
  TCanvas *ccheckAlignmentSignalOverNoise4[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;
      
      // Find the pedestal run from which pedestals have been calculated and get the noise of each Beetle channel.
      string *PedFilename = new string();
      TBranch *b_PedFilename = 0;
      
      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("PedFilename",&PedFilename,&b_PedFilename);
	
	Header->GetEntry();
	
	cout << "Pedestals calculated from: " << *PedFilename << endl; 
	
	// Get the noise of each Beetle channel.
	// Open input pedestal&noise text file.
	FILE *inputText = fopen(PedFilename->c_str(),"r");
	
	// Read noise.  
	for (int iChannel=0; iChannel<N; ++iChannel) {
	  fscanf(inputText,"%*d%*c%*f%*c%*f%*c%f%*c%*f%*c",&(noise[iChannel]));
	  // cout << iChannel << ", " << noise[iChannel] << endl;
	}
	
	// Close input pedestal&noise text file.
	fclose(inputText);
      }
      
      // Beetle channels where to look for signal.
      int ch1 = strip;
      int ch2 = strip+NSkip;
      int ch3 = strip-NSkip;
      int ch4;   
      if (direction == "left")
	ch4 = strip-2*NSkip;
      else if (direction == "right")
	ch4 = strip+2*NSkip;
      
      // Four adjacent strips.
      TH1D *hist1 = (TH1D *)input->Get(Form("hADCPedSub%d",ch1));
      TH1D *hist2 = (TH1D *)input->Get(Form("hADCPedSub%d",ch2));
      TH1D *hist3 = (TH1D *)input->Get(Form("hADCPedSub%d",ch3));
      TH1D *hist4 = (TH1D *)input->Get(Form("hADCPedSub%d",ch4));

      signalOverNoise4[xcounter] = (hist1->GetMean())/noise[ch1]+(hist2->GetMean())/noise[ch2]+(hist3->GetMean())/noise[ch3]+(hist4->GetMean())/noise[ch4];

      // 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.
    gcheckAlignmentSignalOverNoise4[zcounter] = new TGraph(stepsx,x,signalOverNoise4);
    InitGraph(gcheckAlignmentSignalOverNoise4[zcounter],Form("Focusing - 4 adjacent strips - z = %d",(int)z[zcounter]),"x (#mum)","S/N");

    // Fit.

    // Right side (increasing y).
    Float_t minRight;
    Float_t maxRight;

    // Left side (decreasing y).
    Float_t minLeft;
    Float_t maxLeft;

    // Array of the parameters necessary for fitting with fitRight and fitLeft.
    // The array contains par0, par1, par2, par3, par0low, par0high, par1low, par1high, par2low, par2high, par3low, par3high for fitRight and the same for fitLeft.
    const Int_t lengthParSCurves = 24;
    Float_t parSCurves[lengthParSCurves];

    assignParInfoSCurves(string(filename),z[zcounter],&minRight,&maxRight,&minLeft,&maxLeft,parSCurves,lengthParSCurves);

    // Right side (increasing y).
    TF1 *fitSCurveRight = new TF1("fitSCurveRight",fRight,minRight,maxRight,4);
    fitSCurve(r,gcheckAlignmentSignalOverNoise4[zcounter],fitSCurveRight,parSCurves,lengthParSCurves,&aRight[zcounter],&muRight[zcounter],&sigmaRight[zcounter],&bRight[zcounter],&uaRight[zcounter],&umuRight[zcounter],&usigmaRight[zcounter],&ubRight[zcounter]);

    // Left side (decreasing y).
    TF1 *fitSCurveLeft = new TF1("fitSCurveLeft",fLeft,minLeft,maxLeft,4);
    fitSCurve(l,gcheckAlignmentSignalOverNoise4[zcounter],fitSCurveLeft,parSCurves,lengthParSCurves,&aLeft[zcounter],&muLeft[zcounter],&sigmaLeft[zcounter],&bLeft[zcounter],&uaLeft[zcounter],&umuLeft[zcounter],&usigmaLeft[zcounter],&ubLeft[zcounter]);

    ccheckAlignmentSignalOverNoise4[zcounter] = new TCanvas(Form("ccheckAlignmentSignalOverNoise4-%s-%dz",filename,(int)z[zcounter]),"",400,300);
    DrawGraphFunc2(ccheckAlignmentSignalOverNoise4[zcounter],gcheckAlignmentSignalOverNoise4[zcounter],fitSCurveRight,fitSCurveLeft,"AP",path_to_figures);
  }



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

    cout << "Right side, z = " << zcounter << ", a = " << aRight[zcounter] << ", mu = " << muRight[zcounter] << ", sigma = " << sigmaRight[zcounter] << ", b = " << bRight[zcounter] << endl;

    cout << "Left side, z = " << zcounter << ", a = " << aLeft[zcounter] << ", mu = " << muLeft[zcounter] << ", sigma = " << sigmaLeft[zcounter] << ", b = " << bLeft[zcounter] << endl;
  }

  // Show fit results for mu and sigma.
  gStyle->SetOptFit(1);
  gStyle->SetFitFormat(".2g");
  
  // Mu - parameters.
  // y = p0+p1*x, with x = z position (microns) and y = mu (microns).
  // The angle <alpha> is the angle with respect to the normal, so it is zero if the laser is perpendicular to the surface of the sensor. 
  Float_t p0Left;
  Float_t up0Left;
  Float_t p1Left;
  Float_t up1Left;

  Float_t alphaLeft;
  Float_t ualphaLeft;

  Float_t p0Right;
  Float_t up0Right;
  Float_t p1Right;
  Float_t up1Right;

  Float_t alphaRight;
  Float_t ualphaRight;

  // Mu - range.
  Float_t minMu = z[0];
  Float_t maxMu = z[stepsz-1];

  // Mu - left edge.

  // Mu - style for left fit results.
  /*
    gStyle->SetStatX(0.65);
    gStyle->SetStatY(0.9);
  */
  gStyle->SetStatX(0.93);
  gStyle->SetStatY(0.9);

  TGraphErrors *gmuLeft = new TGraphErrors(stepsz,z,muLeft,uz,umuLeft);
  InitGraphErrors(gmuLeft,"Focusing - left side","z (#mum)","#mu_{left} (#mum)");
  TCanvas *cmuLeft = new TCanvas("cmuLeft","",400,300);
  gmuLeft->GetXaxis()->SetLimits(z[0]-100.,z[stepsz-1]+100.);

  TF1 *fitMuLeft = new TF1("fitMuLeft","pol1",minMu,maxMu);
  fitMu(l,gmuLeft,fitMuLeft,&p0Left,&p1Left,&up0Left,&up1Left);

  alphaLeft = atan(p1Left);
  ualphaLeft = up1Left/(1+p1Left*p1Left);

  DrawGraphErrors(cmuLeft,gmuLeft,"AP",path_to_figures);

  // Mu - right edge.

  // Mu - style for right fit results.
  gStyle->SetStatX(0.93);
  gStyle->SetStatY(0.9);

  TGraphErrors *gmuRight = new TGraphErrors(stepsz,z,muRight,uz,umuRight);
  InitGraphErrors(gmuRight,"Focusing - right side","z (#mum)","#mu_{right} (#mum)");
  TCanvas *cmuRight = new TCanvas("cmuRight","",400,300);
  gmuRight->GetXaxis()->SetLimits(z[0]-100.,z[stepsz-1]+100.);

  TF1 *fitMuRight = new TF1("fitMuRight","pol1",minMu,maxMu);
  fitMu(r,gmuRight,fitMuRight,&p0Right,&p1Right,&up0Right,&up1Right);

  alphaRight = atan(p1Right);
  ualphaRight = up1Right/(1+p1Right*p1Right);

  DrawGraphErrors(cmuRight,gmuRight,"AP",path_to_figures);

  // Sigma - style for fit results.
  gStyle->SetStatX(0.65);
  gStyle->SetStatY(0.9);

  // Sigma - parameters.
  // y = p0*|x-p1|+p2, with x = z position (microns) and y = sigma (microns).
  Float_t q0Left;
  Float_t uq0Left;
  Float_t q1Left;
  Float_t uq1Left;
  Float_t q2Left;
  Float_t uq2Left;

  Float_t q0Right;
  Float_t uq0Right;
  Float_t q1Right;
  Float_t uq1Right;
  Float_t q2Right;
  Float_t uq2Right;

  // Sigma - range.
  Float_t minSigma = z[0];
  Float_t maxSigma = z[stepsz-1];

  // Array of the parameters necessary for fitting with fitSigmaLeft and fitSigmaRight.
  // The array contains par0, par1, par2, par0low, par0high, par1low, par1high, par2low, and par2high for fitSigmaLeft and the same for fitSigmaRight.
  const Int_t lengthParSigma = 18;
  Float_t parSigma[lengthParSigma];
  
  assignParInfoSigma(string(filename),parSigma,lengthParSigma);

  // Sigma - left edge.
  TGraphErrors *gsigmaLeft = new TGraphErrors(stepsz,z,sigmaLeft,uz,usigmaLeft);
  InitGraphErrors(gsigmaLeft,"Focusing - left side","z (#mum)","#sigma_{left} (#mum)");
  gsigmaLeft->GetXaxis()->SetLimits(z[0]-100.,z[stepsz-1]+100.);

  TF1 *fitSigmaLeft = new TF1("fitSigmaLeft",fSigma,minSigma,maxSigma,3);
  fitSigma(l,gsigmaLeft,fitSigmaLeft,parSigma,lengthParSigma,&q0Left,&q1Left,&q2Left,&uq0Left,&uq1Left,&uq2Left);

  TCanvas *csigmaLeft = new TCanvas("csigmaLeft","",400,300);
  DrawGraphErrors(csigmaLeft,gsigmaLeft,"AP",path_to_figures);

  // Sigma - right edge.
  TGraphErrors *gsigmaRight = new TGraphErrors(stepsz,z,sigmaRight,uz,usigmaRight);
  InitGraphErrors(gsigmaRight,"Focusing - right side","z (#mum)","#sigma_{right} (#mum)");
  gsigmaRight->GetXaxis()->SetLimits(z[0]-100.,z[stepsz-1]+100.);

  TF1 *fitSigmaRight = new TF1("fitSigmaRight",fSigma,minSigma,maxSigma,3);
  fitSigma(r,gsigmaRight,fitSigmaRight,parSigma,lengthParSigma,&q0Right,&q1Right,&q2Right,&uq0Right,&uq1Right,&uq2Right);

  TCanvas *csigmaRight = new TCanvas("csigmaRight","",400,300);
  DrawGraphErrors(csigmaRight,gsigmaRight,"AP",path_to_figures);

  gStyle->SetOptFit(0);
  gStyle->SetFitFormat("5.4g");

  // Save parameters.
  TTree *tFitResults = new TTree("tFitResults","tFitResults");

  // Left edge.
  tFitResults->Branch("p0Left",&p0Left);
  tFitResults->Branch("p1Left",&p1Left);
  tFitResults->Branch("alphaLeft",&alphaLeft);

  tFitResults->Branch("up0Left",&up0Left);
  tFitResults->Branch("up1Left",&up1Left);
  tFitResults->Branch("ualphaLeft",&ualphaLeft);

  tFitResults->Branch("q0Left",&q0Left);
  tFitResults->Branch("q1Left",&q1Left);
  tFitResults->Branch("q2Left",&q2Left);

  tFitResults->Branch("uq0Left",&uq0Left);
  tFitResults->Branch("uq1Left",&uq1Left);
  tFitResults->Branch("uq2Left",&uq2Left);

  // Right edge.
  tFitResults->Branch("p0Right",&p0Right);
  tFitResults->Branch("p1Right",&p1Right);
  tFitResults->Branch("alphaRight",&alphaRight);

  tFitResults->Branch("up0Right",&up0Right);
  tFitResults->Branch("up1Right",&up1Right);
  tFitResults->Branch("ualphaRight",&ualphaRight);

  tFitResults->Branch("q0Right",&q0Right);
  tFitResults->Branch("q1Right",&q1Right);
  tFitResults->Branch("q2Right",&q2Right);

  tFitResults->Branch("uq0Right",&uq0Right);
  tFitResults->Branch("uq1Right",&uq1Right);
  tFitResults->Branch("uq2Right",&uq2Right);

  tFitResults->Fill();
  tFitResults->Print();

  // Print parameters.
  cout << "**********************************" << endl;

  cout << "Fit results - left edge" << endl;

  cout << "" << endl;

  cout << "mu = p0+p1*z, with:" << endl;
  cout << "- p0 = (" << p0Left << " p/m " << up0Left << ") microns" << endl;
  cout << "- p1 = " << p1Left << " p/m " << up1Left << endl;
  cout << "Angle with respect to the normal: (" << alphaLeft << " p/m " << ualphaLeft << ") rad" << endl; 

  cout << "" << endl;

  cout << "sigma = q0*|z-q1|+q2, with:" << endl;
  cout << "- q0 = " << q0Left << " p/m " << uq0Left << endl;
  cout << "- q1 = (" << q1Left << " p/m " << uq1Left << ") microns ---> position of the focus along z" << endl;
  cout << "- q2 = (" << q2Left << " p/m " << uq2Left << ") microns ---> minimum sigma" << endl;

  cout << "" << endl;

  cout << "Fit results - right edge" << endl;
  cout << "mu = p0+p1*z, with:" << endl;
  cout << "- p0 = (" << p0Right << " p/m " << up0Right << ") microns" << endl;
  cout << "- p1 = " << p1Right << " p/m " << up1Right << endl;
  cout << "Angle with respect to the normal: (" << alphaRight << " p/m " << ualphaRight << ") rad" << endl; 

  cout << "" << endl;

  cout << "sigma = q0*|z-q1|+q2, with:" << endl;
  cout << "- q0 = " << q0Right << " p/m " << uq0Right << endl;
  cout << "- q1 = (" << q1Right << " p/m " << uq1Right << ") microns ---> position of the focus along z" << endl;
  cout << "- q2 = (" << q2Right << " p/m " << uq2Right << ") microns ---> minimum sigma" << endl;

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



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

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

  return;
}