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

/*
SCurve contains some useful functions related to s-curves:
- fLeft describes the falling s-curve;
- fRight describes the rising s-curve;
- fSigma describes the sigma as a function of the z position;
- assignParInfoSCurves sets the starting values and the limits of the parameters of the two erf functions describing the s-curves according to the filename of the measurement;
- assignParInfoSigma sets the starting values and the limits of the parameters of the two sigma functions according to the filename of the measurement;
- assignParInfoEta sets the starting values and the limits of the parameters of the erf function describing the charge sharing coefficient according to the filename of the measurement;
- fitSCurve fits the s-curve and returns the result of the fit;
- fitSigma fits the sigma as a function of the z position and returns the result of the fit.
*/

// Header guard.
#ifndef __SCURVE_C_INCLUDED__
#define __SCURVE_C_INCLUDED__

#include "Lib.C"
#include "lhcbStyle.C"
#include "Style.C"
#include "Par.C"



// Declarations.

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

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

Double_t fSigma(Double_t *x, Double_t *par);

void assignParInfoSCurves(string filename, Float_t zvalue, Float_t *minRight, Float_t *maxRight, Float_t *minLeft, Float_t *maxLeft, Float_t parSCurves[], Int_t lengthParSCurves);
void assignParInfoSigma(string filename, Float_t parSigma[], Int_t lengthParSigma);
void assignParInfoEta(string filename, Float_t parEta[], Int_t lengthParEta);

void fitSCurve(string edge, TGraph *graph, TF1 *f, Float_t parSCurves[], Int_t lengthParSCurves, Float_t *a, Float_t *mu, Float_t *sigma, Float_t *b, Float_t *ua, Float_t *umu, Float_t *usigma, Float_t *ub);

void fitSigma(string edge, TGraph *graph, TF1 *f, Float_t parSigma[], Int_t lengthParSigma, Float_t *q0, Float_t *q1, Float_t *q2, Float_t *uq0, Float_t *uq1, Float_t *uq2);

// Definitions.

// fLeft has four parameters:
// - par[0] is the y-spread around the middle of the curve;
// - par[1] is the x-coordinate of the middle of the curve;
// - par[2] is the steepness of the curve;
// - par[3] is the y-coordinate of the middle of the curve.
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:
// - par[0] is the y-spread around the middle of the curve;
// - par[1] is the x-coordinate of the middle of the curve;
// - par[2] is the steepness of the curve;
// - par[3] is the y-coordinate of the middle of the curve.
Double_t fRight(Double_t *x, Double_t *par)
{
return par[0]*erf((x[0]-par[1])/(sqrt(2.)*par[2]))+par[3];
}

// fSigma is used to fit the sigma as a function of the z position and has three parameters:
// - par[0] is the angular coefficient of the two straight lines;
// - par[1] is the x-coordinate of the lowest point;
// - par[2] is the y-coordinate of the lowest point.
Double_t fSigma(Double_t *x, Double_t *par)
{
  return par[0]*abs(x[0]-par[1])+par[2];
}

// assignParInfoSCurves has the following parameters:
// - string filename, the filename of the measurement;
// - Float_t zvalue, the position along z;
// - Float_t *minRight, the lower limit of the fitting range for the fRight function;
// - Float_t *maxRight, the upper limit of the fitting range for the fRight function;
// - Float_t *minLeft, the lower limit of the fitting range for the fLeft function;
// - Float_t *maxLeft, the upper limit of the fitting range for the fLeft function;
// - Float_t parSCurves[], the array containing the parameters of the s-curves;
// - Float_t lengthSCurves, the length of the array containing the parameters of the s-curves.
void assignParInfoSCurves(string filename, Float_t zvalue, Float_t *minRight, Float_t *maxRight, Float_t *minLeft, Float_t *maxLeft, Float_t parSCurves[], Int_t lengthParSCurves)
{
  if (filename == "ProcessRawData-201501093-216") {
    // Right.
    if (zvalue < 120.)
    {
      *minRight = 70.;
      *maxRight = 180.;
    }
    else
    {
      *minRight = 60.;
      *maxRight = 140.;
    }
    
    // Left.  
    if (zvalue < 40.)
    {
      *minLeft = 180.;
      *maxLeft = 280.;
    }
    else if (zvalue < 120.)
    {
      *minLeft = 180.;
      *maxLeft = 270.;
    }
    else
    {
      *minLeft = 180.;
      *maxLeft = 260.;
    }

    // Right.
    if (zvalue < 120.)
      parSCurves[0] = 125.;
    else 
      parSCurves[0] = 140.;
    parSCurves[1] = 95.;
    parSCurves[2] = 10.;
    parSCurves[3] = 200.;

    parSCurves[4] = 90.;
    parSCurves[5] = 170.;
    parSCurves[6] = 70.;
    parSCurves[7] = 120.;
    parSCurves[8] = 3.;
    parSCurves[9] = 20.;
    parSCurves[10] = 150.;
    parSCurves[11] = 220.;

    // Left.
    if (zvalue < 120.)
      parSCurves[12] = 125.;
    else
      parSCurves[12] = 140.;
    parSCurves[13] = 230.;
    parSCurves[14] = 10.;
    parSCurves[15] = 200.;

    // fitLeft->SetParLimits(0,90.,170.);
    parSCurves[18] = 215.;
    parSCurves[19] = 250.;
    parSCurves[20] = 3.;
    parSCurves[21] = 20.;
    // fitLeft->SetParLimits(3,150.,220.);
  }



  // It does not look very nice.
  else if (filename == "ProcessRawData-20150111-216") {
    // Right.
    if (zvalue < 105.)
    {
      *minRight = 60.;
      *maxRight = 120.;
    }
    else
    {
      *minRight = 50.;
      *maxRight = 120.;
    }
  
    // Left.
      *minLeft = 180.;
      *maxLeft = 250.;

    // Right.
    parSCurves[0] = 140.;
    parSCurves[1] = 85.;
    parSCurves[2] = 6.;
    parSCurves[3] = 175.;

    parSCurves[4] = 110.;
    parSCurves[5] = 170.;
    parSCurves[6] = 70.;
    parSCurves[7] = 100.;
    parSCurves[8] = 3.;
    parSCurves[9] = 10.;
    parSCurves[10] = 150.;
    parSCurves[11] = 200.;

    // Left.
    parSCurves[12] = 140.;
    parSCurves[13] = 220.;
    parSCurves[14] = 6.;
    parSCurves[15] = 175.;

    // fitLeft->SetParLimits(0,90.,170.);
    parSCurves[18] = 210.;
    parSCurves[19] = 230.;
    parSCurves[20] = 3.;
    parSCurves[21] = 10.;
    // fitLeft->SetParLimits(3,150.,220.);
  }



  else if (filename == "ProcessRawData-201501124-216") {
    // Right.
    if (zvalue < 80.)
    {
      *minRight = -8440.;
      *maxRight = -8350.;
    }
    else if (zvalue < 180.)
    {
      *minRight = -8450.;
      *maxRight = -8370.;
    }
    else
    {
      *minRight = -8460.;
      *maxRight = -8370.;
    }
 
    // Left.
    if (zvalue == 100.)
    {
      *minLeft = -8320.;
      *maxLeft = -8240.;
    }
    else if (zvalue < 120.)
    {
      *minLeft = -8350.;
      *maxLeft = -8240.;
    }
    else if (zvalue < 180.)
    {
      *minLeft = -8350.;
      *maxLeft = -8250.;
    }
    else
    {
      *minLeft = -8360.;
      *maxLeft = -8260.;
    }

    // Right.
    if (zvalue < 100.)
      parSCurves[0] = 100.;
    else
      parSCurves[0] = 130.;
    if (zvalue < 100.)
      parSCurves[1] = -8400.;
    else 
      parSCurves[1] = -8420.;
    parSCurves[2] = 10.;
    if (zvalue < 100.)
      parSCurves[3] = 200.;
    else
      parSCurves[3] = 170.;

    parSCurves[4] = 80.;
    parSCurves[5] = 140.;
    parSCurves[6] = -8450.;
    parSCurves[7] = -8370.;
    parSCurves[8] = 3.;
    parSCurves[9] = 20.;
    parSCurves[10] = 150.;
    parSCurves[11] = 250.;

    // Left.
    if (zvalue < 100.)
      parSCurves[12] = 100.;
    else
      parSCurves[12] = 135.;
    parSCurves[13] = -8280.;
    parSCurves[14] = 10.;
    if (zvalue < 100.)
      parSCurves[15] = 200.;
    else
      parSCurves[15] = 180.;

    // fitLeft->SetParLimits(0,90.,170.);
    parSCurves[18] = -8300.;
    parSCurves[19] = -8250.;
    parSCurves[20] = 3.;
    parSCurves[21] = 20.;
    // fitLeft->SetParLimits(3,150.,220.);
  }



  else if (filename == "ProcessRawData-20150221-216") {
    // Right.
    if (zvalue < 180.)
    {
      *minRight = 60.;
      *maxRight = 160.;
    }
    else
    {
      *minRight = 50.;
      *maxRight = 160.;
    }
  
    // Left.
    if (zvalue < 180.)
    {
      *minLeft = 160.;
      *maxLeft = 260.;
    }
    else
    {
      *minLeft = 160.;
      *maxLeft = 250.;
    }

    // Right.
    parSCurves[0] = 15.;
    parSCurves[1] = 100.;
    parSCurves[2] = 6.;
    parSCurves[3] = 12.;

    // parSCurves[4] = 10.;
    // parSCurves[5] = 20.;
    parSCurves[6] = 70.;
    parSCurves[7] = 110.;
    parSCurves[8] = 3.;
    parSCurves[9] = 15.;
    // parSCurves[10] = 8.;
    // parSCurves[11] = 18.;

    // Left.
    parSCurves[12] = 17.;
    parSCurves[13] = 220.;
    parSCurves[14] = 6.;
    parSCurves[15] = 10.;

    // fitLeft->SetParLimits(0,90.,170.);
    parSCurves[18] = 200.;
    parSCurves[19] = 230.;
    parSCurves[20] = 3.;
    parSCurves[21] = 15.;
    // fitLeft->SetParLimits(3,150.,220.);
  }



  else if (filename == "ProcessRawData-201502235-216") {
    // Right.
    if (zvalue < 80.)
    {
      *minRight = 120.;
      *maxRight = 240.;
    }
    else
    {
      *minRight = 100.;
      *maxRight = 230.;
    }
  
    // Left.
    if (zvalue < 80.)
    {
      *minLeft = 240.;
      *maxLeft = 330.;
    }
    else
    {
      *minLeft = 230.;
      *maxLeft = 320.;
    }

    // Right.
    parSCurves[0] = 16.;
    parSCurves[1] = 160.;
    parSCurves[2] = 6.;
    if (zvalue < 60.)
      parSCurves[3] = 12.;
    else
      parSCurves[3] = 10.;

    // parSCurves[4] = 10.;
    // parSCurves[5] = 20.;
    parSCurves[6] = 120.;
    parSCurves[7] = 200.;
    parSCurves[8] = 3.;
    parSCurves[9] = 15.;
    // parSCurves[10] = 8.;
    // parSCurves[11] = 18.;

    // Left.
    parSCurves[12] = 16.;
    parSCurves[13] = 280.;
    parSCurves[14] = 6.;
    parSCurves[15] = 12.;

    // fitLeft->SetParLimits(0,90.,170.);
    parSCurves[18] = 260.;
    parSCurves[19] = 300.;
    parSCurves[20] = 3.;
    parSCurves[21] = 15.;
    // fitLeft->SetParLimits(3,150.,220.);
  }



  else if ((filename == "ProcessRawData-201502253-216") || (filename == "ProcessRawData-20150226-216") || (filename == "ProcessRawData-201502262-216")) {
    // Right.
      *minRight = 30.;
      *maxRight = 140.;
  
    // Left.
      if (zvalue < 150.) {
	*minLeft = 140.;
	*maxLeft = 235;
      }
      else if (zvalue < 230.) {
	*minLeft = 140.;
	*maxLeft = 240.;
      }
      else {
	*minLeft = 140.;
	*maxLeft = 230.;
      }

    // Right.
      if (zvalue < 130.)
	parSCurves[0] = 11.;
      else
	parSCurves[0] = 16.;
      if (zvalue < 110.)
	parSCurves[1] = 80.;
      else if (zvalue < 150.)
	parSCurves[1] = 70.;
      else
	parSCurves[1] = 60.;
    if (zvalue<130.)
      parSCurves[2] = 10.;
    else
      parSCurves[2] = 6.;
    parSCurves[3] = 10.;

    // parSCurves[4] = 10.;
    // parSCurves[5] = 20.;
    parSCurves[6] = 20.;
    parSCurves[7] = 100.;
    parSCurves[8] = 3.;
    parSCurves[9] = 20.;
    // parSCurves[10] = 8.;
    // parSCurves[11] = 18.;

    // Left.
    if (zvalue < 130.)
      parSCurves[12] = 11.;
    else
      parSCurves[12] = 16.;
    parSCurves[13] = 190.;
    if (zvalue < 130.)
      parSCurves[14] = 10.;
    else
      parSCurves[14] = 6.;
    parSCurves[15] = 10.;

    // fitLeft->SetParLimits(0,90.,170.);
    parSCurves[18] = 170.;
    parSCurves[19] = 210.;
    parSCurves[20] = 3.;
    parSCurves[21] = 20.;
    // fitLeft->SetParLimits(3,150.,220.);
  }

  return;
}

// assignParInfoSigma has the following parameters:
// - string filename, the filename of the measurement;
// - Float_t parSigma[], the array containing the parameters of the sigma function;
// - Float_t lengthParSigma, the length of the array containing the parameters of the sigma function.
void assignParInfoSigma(string filename, Float_t parSigma[], Int_t lengthParSigma)
{
  if ((filename == "ProcessRawData-201501093-216") || (filename == "ProcessRawData-20150111-216") || (filename == "ProcessRawData-201501124-216") || (filename == "ProcessRawData-20150221-216")) {
    // Sigma left.
    parSigma[0] = 5.;
    parSigma[1] = 600.;
    parSigma[2] = 0.025;
    
    parSigma[3] = 0.;
    parSigma[4] = 10.;
    
    parSigma[5] = 500.;
    parSigma[6] = 700.;
    
    // parSigma[7] = ?.;
    // parSigma[8] = ?.;

    // Sigma right.
    parSigma[9] = 5.;
    parSigma[10] = 600.;
    parSigma[11] = 0.025;
    
    parSigma[12] = 0.;
    parSigma[13] = 10.;
    
    parSigma[14] = 500.;
    parSigma[15] = 700.;
    
    // parSigma[16] = ?.;
    // parSigma[17] = ?.;
  }

  else if (filename == "ProcessRawData-201502235-216") {
    // Sigma left.
    parSigma[0] = 5.;
    parSigma[1] = 100.;
    parSigma[2] = 6.;
    
    parSigma[3] = 0.;
    parSigma[4] = 10.;
    
    parSigma[5] = 0.;
    parSigma[6] = 200.;
    
    // parSigma[7] = ?.;
    // parSigma[8] = ?.;

    // Sigma right.
    parSigma[9] = 5.;
    parSigma[10] = 100.;
    parSigma[11] = 6.;
    
    parSigma[12] = 0.;
    parSigma[13] = 10.;
    
    parSigma[14] = 0.;
    parSigma[15] = 200.;
    
    // parSigma[16] = ?.;
    // parSigma[17] = ?.;
  }

  else if ((filename == "ProcessRawData-201502253-216") || (filename == "ProcessRawData-20150226-216")) {
    // Sigma left.
    parSigma[0] = 0.01;
    parSigma[1] = 850.;
    parSigma[2] = 5.;
    
    parSigma[3] = 0.001;
    parSigma[4] = 1.;
    
    parSigma[5] = 800.;
    parSigma[6] = 900.;
    
    // parSigma[7] = ?.;
    // parSigma[8] = ?.;

    // Sigma right.
    parSigma[9] = 0.01;
    parSigma[10] = 850.;
    parSigma[11] = 5.;
    
    parSigma[12] = 0.001;
    parSigma[13] = 1.;
    
    parSigma[14] = 800.;
    parSigma[15] = 900.;
    
    // parSigma[16] = ?.;
    // parSigma[17] = ?.;
  }

  return;
}

// assignParInfoEta has the following parameters:
// - string filename, the filename of the measurement;
// - Float_t parEta[], the array containing the parameters of the erf function;
// - Float_t lengthParEta, the length of the array containing the parameters of the erf function.
void assignParInfoEta(string filename, Float_t parEta[], Int_t lengthParEta) {
 if ((filename == "ProcessRawData-201502262-216")) {
    parEta[0] = 1.5;
    parEta[1] = 140.;
    parEta[2] = 5.;
    parEta[3] = 0.;

    /*
    parEta[4] = ;
    parEta[5] = ;

    parEta[6] = ;
    parEta[7] = ;
    
    parEta[8] = ;
    parEta[9] = ;
   
    parEta[10] = ;
    parEta[11] = ;
    */
  }

  return;
}

void fitSCurve(string edge, TGraph *graph, TF1 *f, Float_t parSCurves[], Int_t lengthParSCurves, Float_t *a, Float_t *mu, Float_t *sigma, Float_t *b, Float_t *ua, Float_t *umu, Float_t *usigma, Float_t *ub) {
  // Set parameter names.
  f->SetParName(0,"a");
  f->SetParName(1,"#mu");
  f->SetParName(2,"#sigma");
  f->SetParName(3,"b");

  // Set parameters and parameter limits.
  if (edge == "right") {
    f->SetParameter(0,parSCurves[0]);
    f->SetParameter(1,parSCurves[1]);
    f->SetParameter(2,parSCurves[2]);
    f->SetParameter(3,parSCurves[3]);
    // f->SetParLimits(0,parSCurves[4],parSCurves[5]);
    f->SetParLimits(1,parSCurves[6],parSCurves[7]);
    f->SetParLimits(2,parSCurves[8],parSCurves[9]);
    // f->SetParLimits(3,parSCurves[10],parSCurves[11]);
  }
  else if ((edge == "left") && (lengthParSCurves>12)) {  
    f->SetParameter(0,parSCurves[12]);
    f->SetParameter(1,parSCurves[13]);
    f->SetParameter(2,parSCurves[14]);
    f->SetParameter(3,parSCurves[15]);
    // f->SetParLimits(0,parSCurves[16],parSCurves[17]);
    f->SetParLimits(1,parSCurves[18],parSCurves[19]);
    f->SetParLimits(2,parSCurves[20],parSCurves[21]);
    // f->SetParLimits(3,parSCurves[22],parSCurves[23]);
  }

  // Fit.
  graph->Fit(f,"BR");
  
  // Get parameters and parameter uncertainties.
  *a = f->GetParameter(0);
  *mu = f->GetParameter(1);
  *sigma = f->GetParameter(2);
  *b = f->GetParameter(3);
  
  *ua = f->GetParError(0);
  *umu = f->GetParError(1);
  *usigma = f->GetParError(2);
  *ub = f->GetParError(3);
  
  cout << edge << " side, a = " << *a << ", mu = " << *mu << ", sigma = " << *sigma << ", b = " << *b << endl;
  
  return;
}

void fitSigma(string edge, TGraph *graph, TF1 *f, Float_t parSigma[], Int_t lengthParSigma, Float_t *q0, Float_t *q1, Float_t *q2, Float_t *uq0, Float_t *uq1, Float_t *uq2) {
  // Set parameter names.
  f->SetParName(0,"{q}_0");
  f->SetParName(1,"{q}_1");
  f->SetParName(2,"{q}_2");

  // Set parameters and parameter limits.
  if (edge == "left") {  
  f->SetParameter(0,parSigma[0]);
  f->SetParameter(1,parSigma[1]);
  f->SetParameter(2,parSigma[2]);
  f->SetParLimits(0,parSigma[3],parSigma[4]);
  f->SetParLimits(1,parSigma[5],parSigma[6]);
  // f->SetParLimits(2,parSigma[7],parSigma[8]);
  }
  else if (edge == "right") {
  f->SetParameter(0,parSigma[9]);
  f->SetParameter(1,parSigma[10]);
  f->SetParameter(2,parSigma[11]);
  f->SetParLimits(0,parSigma[12],parSigma[13]);
  f->SetParLimits(1,parSigma[14],parSigma[15]);
  // f->SetParLimits(2,parSigma[16],parSigma[17]);
  }

  // Fit.
  f->SetLineColor(kRed);
  graph->Fit(f,"BR");
  
  // Get parameters and parameter uncertainties.
  *q0 = f->GetParameter(0);
  *q1 = f->GetParameter(1);
  *q2 = f->GetParameter(2);

  *uq0 = f->GetParError(0);
  *uq1 = f->GetParError(1);
  *uq2 = f->GetParError(2);

  return;
}

#endif