//************************************************
// Author: Federica Lionetto
// Created on: 24/02/2015
//************************************************

/*
CCE reads all the ROOT files of a given CCE data taking and calculates the CCE.
In a CCE data taking the position of the laser along z is fixed while the position of the laser along x changes.
In the first part, the s-curve corresponding to a scan across two strips is fitted in order to have access to the subset of the full x range where the laser is not hitting the strips. This feature is not used at the moment.
In the second part, the CCE is plot as a function of the x position.
With CCE I mean:
- the S/N of the strip on the left of the laser;
- the S/N of the strip on the right of the laser;
- the sum of the S/N of the strips on the left and on the right of the laser; 
- the sum of the S/N of the four strips around the laser.
In the third part, the charge sharing coefficient eta is calculated at each x position.
The charge sharing coefficient eta is defined as (S/N_i+1 - S/N_i)/(S/N_i+1 + S/N_i), where S/N is the signal over noise ratio on strip i (left) or i+1 (right). If the charge carriers are collected only by the strip on the left then eta = -1 while if the charge carriers are collected only by the strip on the right then eta = 1.
The charge sharing coefficient eta as a function of the x position is fitted by the same function that describes the s-curve with the rising edge.
The fit results are written to a tree called tFitResults.

The uncertainty on S/N is calculated as
sigma_(S/N) = (S/N)*(sigma_S/S+sigma_N/N)
by using the fact that the relative uncertainty on S/N is the sum of the relative uncertainties on S and N.
The S/N of different strips are considered independent variables, so the uncertainty of a function of the S/N of different strips is calculated accordingly.

The CCE data taking is identified by the following information:
- <sensor>, that is, the type of sensor (Hans410, ...);
- <measurement>, that is, the type of measurement (CCEBiasVoltageScan, CCEAttenuationScan, CCEDelayScan, ...), defining the folder where to look for the input files;
- <filename>, that is, the filename excluding the position value and the run type.
- <z>, that is, the position along z;
- <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).

Compile with:

make

Run with:

./CCE [sensor] [measurement] [filename] [z] [firstx] [lastx] [stepx] [additional folder]

For example:

./CCE Hans410 CCEBiasVoltageScan ProcessRawData-20141121 0 0 75 3

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 CCE(char *sensor, char *measurement, char *filename, const Float_t z, const Float_t firstx, const Float_t lastx, const Float_t stepx, 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 < 8)
  {
    cout << "**************************************************" << endl;

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

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

void CCE(char *sensor, char *measurement, char *filename, const Float_t z, const Float_t firstx, const Float_t lastx, const Float_t stepx, char *externalPath)
{
  cout << "**************************************************" << endl;
  cout << "Measuring CCE..." << 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)+"/"+string(measurement);
  string outputDirectory = "/disk/groups/hep/flionett/TestStand/AnalysisResults/"+string(sensor)+"/"+string(measurement);
  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());

  ostringstream convertx;
  ostringstream convertz;
  string tempx;
  string tempz;

  convertz.str("");
  convertz << z;
  tempz = convertz.str();  

  cout << "Figures stored in: " << outputDirectory+"/Figures/"+filename+"-"+tempz+"z" << 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+"-"+tempz+"z";
  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);

  // Beetle channels where to look for signal.
  // In case of a CCEBiasVoltageScan, CCEAttenuationScan, or CCEDelayScan, we want these Beetle channels to be always the same, as defined from the measurement in the default configuration. So, an optional argument can be passed to CCE.C, corresponding to the full path of the text file with the values of ch1, ch2, ch3, ch4. If this optional argument is not passed, then CCE.C determines ch1, ch2, ch3, and ch4 based on the information available in the input files. 
  string inputFindStrip;
  string outputFindStrip;

  int strip;
  string direction;

  int ch1;
  int ch2;
  int ch3;
  int ch4;   
  
  char char_input_ROOT[200];
  string input_ROOT;
  string output_ROOT;

  Float_t signal[4]; // ch1, ch2, ch3, ch4.
  Float_t usignal[4]; // ch1, ch2, ch3, ch4.
  Float_t ursignal[4]; // ch1, ch2, ch3, ch4.
  Float_t noise[N];
  Float_t unoise[N];
  Float_t urnoise[N];

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

  // CCE.
  Float_t signalOverNoise1Left[steps];
  Float_t signalOverNoise1Right[steps];
  Float_t signalOverNoise2[steps];
  Float_t signalOverNoise4[steps];

  // The uncertainty on S/N is calculated as
  // sigma_(S/N) = (S/N)*(sigma_S/S+sigma_N/N)
  // by using the fact that the relative uncertainty on S/N is the sum of the relative uncertainties on S and N.
  Float_t usignalOverNoise1Left[steps];
  Float_t usignalOverNoise1Right[steps];
  Float_t usignalOverNoise2[steps];
  Float_t usignalOverNoise4[steps];

  // Charge sharing coefficient eta.
  Float_t etaSignalOverNoise[steps];
  Float_t uetaSignalOverNoise[steps];

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

  Float_t aRight;
  Float_t bRight;
  Float_t muRight;
  Float_t sigmaRight;

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

  Float_t uaRight;
  Float_t ubRight;
  Float_t umuRight;
  Float_t usigmaRight;

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

  // Open input data ROOT file.
  // I did not find a better way to create a string containing the bias voltage.
  int stepFindStrip = (int)steps/3;
  int xFindStrip = firstx+stepFindStrip*stepx;

  sprintf(char_input_ROOT,"%s/%s-%dx-%dz-las.root",inputDirectory.c_str(),filename,(int)xFindStrip,(int)z);
  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;
  }

  // Determine the four adjacent strips that the laser hits.  
  convertx.str("");
  convertx << xFindStrip;
  tempx = convertx.str();  

  inputFindStrip = input_ROOT;
  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;
  
  // Beetle channels where to look for signal.
  ch1 = strip;
  ch2 = strip+NSkip;
  ch3 = strip-NSkip;
  if (direction == "left")
    ch4 = strip-2*NSkip;
  else if (direction == "right")
    ch4 = strip+2*NSkip;
  
  cout << "Beetle channels: " << endl;
  cout << ch1 << endl;
  cout << ch2 << endl;
  cout << ch3 << endl;
  cout << ch4 << endl;
 
  for (Int_t step=0;step<steps;step++) {
    x[step] = firstx+step*stepx;
    ux[step] = 0.;

    // Open input data ROOT files.
    // I did not find a better way to create a string containing the bias voltage.
    sprintf(char_input_ROOT,"%s/%s-%dx-%dz-las.root",inputDirectory.c_str(),filename,(int)x[step],(int)z);
    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 #" << step+1 << ": " << input_ROOT << endl;
    TFile *input = TFile::Open(TString(input_ROOT));

    // 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]),&(unoise[iChannel]));
	// cout << iChannel << ", " << noise[iChannel] << ", " << unoise[iChannel] << endl;
	urnoise[iChannel] = unoise[iChannel]/noise[iChannel];
      }
  
      // Close input pedestal&noise text file.
      fclose(inputText);
    }
    
    // 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));

    signal[0] = hist1->GetMean();
    signal[1] = hist2->GetMean();
    signal[2] = hist3->GetMean();
    signal[3] = hist4->GetMean();

    usignal[0] = hist1->GetMeanError();
    usignal[1] = hist2->GetMeanError();
    usignal[2] = hist3->GetMeanError();
    usignal[3] = hist4->GetMeanError();

    ursignal[0] = usignal[0]/signal[0];
    ursignal[1] = usignal[1]/signal[1];
    ursignal[2] = usignal[2]/signal[2];
    ursignal[3] = usignal[3]/signal[3];

    if (direction == "left") {
      signalOverNoise1Left[step] = signal[2]/noise[ch3];
      signalOverNoise1Right[step] = signal[0]/noise[ch1];
      // Uncertainties.
      usignalOverNoise1Left[step] = signalOverNoise1Left[step]*(ursignal[2]+urnoise[ch3]);
      usignalOverNoise1Right[step] = signalOverNoise1Right[step]*(ursignal[0]+urnoise[ch1]);
    }
    else if (direction == "right") {
      signalOverNoise1Left[step] = signal[0]/noise[ch1];
      signalOverNoise1Right[step] = signal[1]/noise[ch2];
      // Uncertainties.
      usignalOverNoise1Left[step] = signalOverNoise1Left[step]*(ursignal[0]+urnoise[ch1]);
      usignalOverNoise1Right[step] = signalOverNoise1Right[step]*(ursignal[1]+urnoise[ch2]);
    }

    signalOverNoise2[step] = signalOverNoise1Left[step]+signalOverNoise1Right[step];
    signalOverNoise4[step] = signal[0]/noise[ch1]+signal[1]/noise[ch2]+signal[2]/noise[ch3]+signal[3]/noise[ch4];
    // Uncertainty.
    usignalOverNoise2[step] = usignalOverNoise1Left[step]+usignalOverNoise1Right[step];
    usignalOverNoise4[step] = signal[0]/noise[ch1]*(ursignal[0]+urnoise[ch1])+signal[1]/noise[ch2]*(ursignal[1]+urnoise[ch2])+signal[2]/noise[ch3]*(ursignal[2]+urnoise[ch3])+signal[3]/noise[ch4]*(ursignal[3]+urnoise[ch4]);



    etaSignalOverNoise[step] = (signalOverNoise1Right[step]-signalOverNoise1Left[step])/(signalOverNoise1Right[step]+signalOverNoise1Left[step]);
    uetaSignalOverNoise[step] = (2.*sqrt(pow((signalOverNoise1Left[step]*usignalOverNoise1Right[step]),2)+pow((signalOverNoise1Right[step]*usignalOverNoise1Left[step]),2)))/pow((signalOverNoise1Right[step]+signalOverNoise1Left[step]),2);

    // Close input data ROOT files.
    input->Close();
  }

  output->cd();

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

  // Four adjacent strips.
  TGraphErrors *gcheckAlignmentSignalOverNoise4 = new TGraphErrors(steps,x,signalOverNoise4,ux,usignalOverNoise4);
  InitGraphErrors(gcheckAlignmentSignalOverNoise4,"Check alignment - 4 adjacent strips","x (#mum)","S/N");

  // Fit the two s-curves and set the range for the CCE measurement.
  // The range for the CCE measurement is defined as [muRight+(muLeft-muRight)*0.05,muLeft-(muLeft-muRight)*0.05].
  // This is not used at the moment!!!

  // 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,&minRight,&maxRight,&minLeft,&maxLeft,parSCurves,lengthParSCurves);

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

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

  TCanvas *ccheckAlignmentSignalOverNoise4 = new TCanvas(Form("ccheckAlignmentSignalOverNoise4-%s",filename),"",400,300);
  DrawGraphErrorsFunc2(ccheckAlignmentSignalOverNoise4,gcheckAlignmentSignalOverNoise4,fitSCurveRight,fitSCurveLeft,"AP",path_to_figures);

  Float_t zMicrons = z*5.; // One step corresponds to 5 microns.

  cout << "Right side, z = " << zMicrons << ", a = " << aRight << ", mu = " << muRight << ", sigma = " << sigmaRight << ", b = " << bRight << endl;
  
  cout << "Left side, z = " << zMicrons << ", a = " << aLeft << ", mu = " << muLeft << ", sigma = " << sigmaLeft << ", b = " << bLeft << endl;
  
  Float_t minPosCCE = 0.;
  Float_t maxPosCCE = 0.;
  if ((string(sensor) == "Hans410") || (string(sensor) == "Hans320")) {
    minPosCCE = muRight+(muLeft-muRight)*0.05;
    maxPosCCE = muLeft-(muLeft-muRight)*0.05; 
  }
  else if (string(sensor) == "ATLAS") {
    minPosCCE = muLeft+(muRight-muLeft)*0.05;
    maxPosCCE = muRight-(muRight-muLeft)*0.05; 
  }

  // Set the four Beetle channels that will form the cluster where to look for collected charge.
  // These are ch1, ch2, ch3, and ch4.

  // CCE.
  // Charge collected by the different strips.
  TCanvas *cCCESignalOverNoise = new TCanvas(Form("cCCESignalOverNoise-%s",filename),"",550,300);

  TMultiGraph *mgCCESignalOverNoise = new TMultiGraph();

  TGraphErrors *gCCESignalOverNoise1Left = new TGraphErrors(steps,x,signalOverNoise1Left,ux,usignalOverNoise1Left);
  TGraphErrors *gCCESignalOverNoise1Right = new TGraphErrors(steps,x,signalOverNoise1Right,ux,usignalOverNoise1Right);
  TGraphErrors *gCCESignalOverNoise2 = new TGraphErrors(steps,x,signalOverNoise2,ux,usignalOverNoise2);
  TGraphErrors *gCCESignalOverNoise4 = new TGraphErrors(steps,x,signalOverNoise4,ux,usignalOverNoise4);

  InitGraphErrors(gCCESignalOverNoise1Left,"CCE","x (#mum)","S/N");
  InitGraphErrors(gCCESignalOverNoise1Right,"CCE","x (#mum)","S/N");
  InitGraphErrors(gCCESignalOverNoise2,"CCE","x (#mum)","S/N");
  InitGraphErrors(gCCESignalOverNoise4,"CCE","x (#mum)","S/N");

  gCCESignalOverNoise1Left->SetMarkerColor(kMagenta);
  gCCESignalOverNoise1Right->SetMarkerColor(kOrange);
  gCCESignalOverNoise2->SetMarkerColor(kRed);
  gCCESignalOverNoise4->SetMarkerColor(kRed+2);

  gCCESignalOverNoise1Left->SetLineColor(kMagenta);
  gCCESignalOverNoise1Right->SetLineColor(kOrange);
  gCCESignalOverNoise2->SetLineColor(kRed);
  gCCESignalOverNoise4->SetLineColor(kRed+2);

  gCCESignalOverNoise1Left->SetLineWidth(2);
  gCCESignalOverNoise1Right->SetLineWidth(2);
  gCCESignalOverNoise2->SetLineWidth(2);
  gCCESignalOverNoise4->SetLineWidth(2);

  TLegend *legCCESignalOverNoise;
  if ((string(sensor) == "Hans410") || (string(sensor) == "Hans320"))
    legCCESignalOverNoise = CreateLegend4(gCCESignalOverNoise1Left,"left strip",gCCESignalOverNoise1Right,"right strip",gCCESignalOverNoise2,"2 strips",gCCESignalOverNoise4,"4 strips","lpw",0.79,0.62,0.92,0.92);
  else if (string(sensor) == "ATLAS")
    legCCESignalOverNoise = CreateLegend4(gCCESignalOverNoise1Left,"left strip",gCCESignalOverNoise1Right,"right strip",gCCESignalOverNoise2,"2 strips",gCCESignalOverNoise4,"4 strips","lpw",0.77,0.21,0.90,0.51);

  mgCCESignalOverNoise->Add(gCCESignalOverNoise1Left);
  mgCCESignalOverNoise->Add(gCCESignalOverNoise1Right);
  mgCCESignalOverNoise->Add(gCCESignalOverNoise2);
  mgCCESignalOverNoise->Add(gCCESignalOverNoise4);

  if (string(sensor) == "Hans410") {
    mgCCESignalOverNoise->SetMinimum(-1.);
    mgCCESignalOverNoise->SetMaximum(30.);
  }
  else if (string(sensor) == "Hans320") {
    mgCCESignalOverNoise->SetMinimum(-1.);
    mgCCESignalOverNoise->SetMaximum(15.);
  }
  else if (string(sensor) == "ATLAS") {
    mgCCESignalOverNoise->SetMinimum(-65.);
    mgCCESignalOverNoise->SetMaximum(5.);
  }

  DrawGraphCompare(cCCESignalOverNoise,mgCCESignalOverNoise,legCCESignalOverNoise,"CCE","x (#mum)","S/N","APCE",path_to_figures);

  // Charge sharing coefficient eta.
  // The charge sharing coefficient eta is defined as (S/N_i+1 - S/N_i)/(S/N_i+1 + S/N_i), where S/N is the signal over noise ratio on strip i (left) or i+1 (right). If the charge carriers are collected only by the strip on the left then eta = -1 while if the charge carriers are collected only by the strip on the right then eta = 1.
  TCanvas *cEtaSignalOverNoise = new TCanvas(Form("cEtaSignalOverNoise-%s",filename),"",400,300);
  TGraphErrors *gEtaSignalOverNoise = new TGraphErrors(steps,x,etaSignalOverNoise,ux,uetaSignalOverNoise);
  InitGraphErrors(gEtaSignalOverNoise,"Charge sharing coefficient","x (#mum)","#eta S/N");

  // Fit with the fRight function (which describes the rising s-curve).
  TF1 *fEtaSignalOverNoise = new TF1("fEtaSignalOverNoise",fRight,minPosCCE,maxPosCCE,4);

  // Array of the parameters necessary for fitting with fitRight.
    // The array contains par0, par1, par2, par3, par0low, par0high, par1low, par1high, par2low, par2high, par3low, par3high.
  const Int_t lengthParEtaSignalOverNoise = 12;
  Float_t parEtaSignalOverNoise[lengthParEtaSignalOverNoise];

  assignParInfoEta(string(filename),parEtaSignalOverNoise,lengthParEtaSignalOverNoise);

  // Parameters of the fitting function.
  Float_t aEtaSignalOverNoise;
  Float_t muEtaSignalOverNoise;
  Float_t sigmaEtaSignalOverNoise;
  Float_t bEtaSignalOverNoise;

  // Uncertainties on the parameters of the fitting function.
  // u = uncertainty.
  Float_t uaEtaSignalOverNoise;
  Float_t umuEtaSignalOverNoise;
  Float_t usigmaEtaSignalOverNoise;
  Float_t ubEtaSignalOverNoise;

  fitSCurve(r,gEtaSignalOverNoise,fEtaSignalOverNoise,parEtaSignalOverNoise,lengthParEtaSignalOverNoise,&aEtaSignalOverNoise,&muEtaSignalOverNoise,&sigmaEtaSignalOverNoise,&bEtaSignalOverNoise,&uaEtaSignalOverNoise,&umuEtaSignalOverNoise,&usigmaEtaSignalOverNoise,&ubEtaSignalOverNoise);

  // Show fit results.
  gStyle->SetOptFit(1);
  gStyle->SetStatX(0.9);
  gStyle->SetStatY(0.55);
  gStyle->SetStatW(0.16);
  gStyle->SetStatH(0.23);
  gStyle->SetFitFormat(".2g");

  gEtaSignalOverNoise->GetXaxis()->SetLimits(minPosCCE,maxPosCCE);

  DrawGraphErrorsFunc(cEtaSignalOverNoise,gEtaSignalOverNoise,fEtaSignalOverNoise,"AP",path_to_figures);

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

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

  tFitResults->Branch("aEtaSignalOverNoise",&aEtaSignalOverNoise);
  tFitResults->Branch("muEtaSignalOverNoise",&muEtaSignalOverNoise);
  tFitResults->Branch("sigmaEtaSignalOverNoise",&sigmaEtaSignalOverNoise);
  tFitResults->Branch("bEtaSignalOverNoise",&bEtaSignalOverNoise);

  tFitResults->Branch("uaEtaSignalOverNoise",&uaEtaSignalOverNoise);
  tFitResults->Branch("umuEtaSignalOverNoise",&umuEtaSignalOverNoise);
  tFitResults->Branch("usigmaEtaSignalOverNoise",&usigmaEtaSignalOverNoise);
  tFitResults->Branch("ubEtaSignalOverNoise",&ubEtaSignalOverNoise);

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

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

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

  return;
}