#include <math.h>

Double_t PI=3.141592654;

Double_t pilarErf(Double_t *x,Double_t *par)
{
  Double_t result = par[0]*(1.+ROOT::Math::erf((x[0]-par[1])/par[2]))+par[3];
  return result;
}

Double_t pilarErfDouble(Double_t x,Double_t *par)
{
  Double_t result = par[0]*(1.+ROOT::Math::erf((x-par[1])/par[2]))+par[3];
  return result;
}

Double_t pilar(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t p3 = par[3];

  Double_t result = 0.;
  result = x[0]-p1;
  result = result/p2;
  result = atan(result);
  result = p0*result;
  result = result+p3;

  return result;
}

Double_t errorFitUp(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t x00 = par[3];
  Double_t x01 = par[4];
  Double_t x02 = par[5];
  Double_t x10 = par[6];
  Double_t x11 = par[7];
  Double_t x12 = par[8];
  Double_t x20 = par[9];
  Double_t x21 = par[10];
  Double_t x22 = par[11];

  Double_t zx0 = atan((x[0]-p1)/p2);
  Double_t zx1 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*(-1./p2); 
  Double_t zx2 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*((p1-x[0])/p2**2); 

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ 
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22; 
  return p0*atan((x[0]-p1)/p2) + sqrt(result);
}

Double_t errorFitUp4(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t p3 = par[3];
  Double_t x00 = par[4];
  Double_t x01 = par[5];
  Double_t x02 = par[6];
  Double_t x03 = par[7];
  Double_t x10 = par[8];
  Double_t x11 = par[9];
  Double_t x12 = par[10];
  Double_t x13 = par[11];
  Double_t x20 = par[12];
  Double_t x21 = par[13];
  Double_t x22 = par[14];
  Double_t x23 = par[15];
  Double_t x30 = par[16];
  Double_t x31 = par[17];
  Double_t x32 = par[18];
  Double_t x33 = par[19];

  Double_t zx0 = atan((x[0]-p1)/p2);
  Double_t zx1 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*(-1./p2); 
  Double_t zx2 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*((p1-x[0])/p2**2); 
  Double_t zx3 = 1.;

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ zx0*zx3*x03+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ zx1*zx3*x13+
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22+ zx2*zx3*x23+  
                    zx3*zx0*x30 + zx3*zx1*x31 + zx3*zx2*x32+ zx3*zx3*x33;
  
  return p0*atan((x[0]-p1)/p2)+p3 + sqrt(result);
}

Double_t errorFitDown(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t x00 = par[3];
  Double_t x01 = par[4];
  Double_t x02 = par[5];
  Double_t x10 = par[6];
  Double_t x11 = par[7];
  Double_t x12 = par[8];
  Double_t x20 = par[9];
  Double_t x21 = par[10];
  Double_t x22 = par[11];

  Double_t zx0 = atan((x[0]-p1)/p2);
  Double_t zx1 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*(-1./p2); 
  Double_t zx2 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*((p1-x[0])/p2**2); 

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ 
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22; 
  return p0*atan((x[0]-p1)/p2) - sqrt(result);
}

Double_t errorFitDown4(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t p3 = par[3];
  Double_t x00 = par[4];
  Double_t x01 = par[5];
  Double_t x02 = par[6];
  Double_t x03 = par[7];
  Double_t x10 = par[8];
  Double_t x11 = par[9];
  Double_t x12 = par[10];
  Double_t x13 = par[11];
  Double_t x20 = par[12];
  Double_t x21 = par[13];
  Double_t x22 = par[14];
  Double_t x23 = par[15];
  Double_t x30 = par[16];
  Double_t x31 = par[17];
  Double_t x32 = par[18];
  Double_t x33 = par[19];

  Double_t zx0 = atan((x[0]-p1)/p2);
  Double_t zx1 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*(-1./p2); 
  Double_t zx2 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*((p1-x[0])/p2**2); 
  Double_t zx3 = 1.;

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ zx0*zx3*x03+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ zx1*zx3*x13+
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22+ zx2*zx3*x23+  
                    zx3*zx0*x30 + zx3*zx1*x31 + zx3*zx2*x32+ zx3*zx3*x33;
  
  return p0*atan((x[0]-p1)/p2)+p3 - sqrt(result);
}

Double_t errorFit(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t x00 = par[3];
  Double_t x01 = par[4];
  Double_t x02 = par[5];
  Double_t x10 = par[6];
  Double_t x11 = par[7];
  Double_t x12 = par[8];
  Double_t x20 = par[9];
  Double_t x21 = par[10];
  Double_t x22 = par[11];

  Double_t zx0 = atan((x[0]-p1)/p2);
  Double_t zx1 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*(-1./p2); 
  Double_t zx2 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*((p1-x[0])/p2**2); 

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ 
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22; 
  return sqrt(result);
}

Double_t errorFit4(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t p3 = par[3];
  Double_t x00 = par[4];
  Double_t x01 = par[5];
  Double_t x02 = par[6];
  Double_t x03 = par[7];
  Double_t x10 = par[8];
  Double_t x11 = par[9];
  Double_t x12 = par[10];
  Double_t x13 = par[11];
  Double_t x20 = par[12];
  Double_t x21 = par[13];
  Double_t x22 = par[14];
  Double_t x23 = par[15];
  Double_t x30 = par[16];
  Double_t x31 = par[17];
  Double_t x32 = par[18];
  Double_t x33 = par[19];

  Double_t zx0 = atan((x[0]-p1)/p2);
  Double_t zx1 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*(-1./p2); 
  Double_t zx2 = p0*(1./( 1.+( (x[0]-p1)/p2  )**2  ))*((p1-x[0])/p2**2); 
  Double_t zx3 = 1.;

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ zx0*zx3*x03+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ zx1*zx3*x13+
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22+ zx2*zx3*x23+  
                    zx3*zx0*x30 + zx3*zx1*x31 + zx3*zx2*x32+ zx3*zx3*x33;
  
  return sqrt(result);
}

Double_t errorFitUpErf(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t x00 = par[3];
  Double_t x01 = par[4];
  Double_t x02 = par[5];
  Double_t x10 = par[6];
  Double_t x11 = par[7];
  Double_t x12 = par[8];
  Double_t x20 = par[9];
  Double_t x21 = par[10];
  Double_t x22 = par[11];

  Double_t zx0 = 1.+ROOT::Math::erf((x[0]-p1)/p2); 
  Double_t zx1 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*(-1./p2);
  Double_t zx2 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*((p1-x[0])/p2**2);

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ 
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22; 
  return par[0]*(1.+ROOT::Math::erf((x[0]-par[1])/par[2])) + sqrt(result);
}

Double_t errorFitUpErf4(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t p3 = par[3];
  Double_t x00 = par[4];
  Double_t x01 = par[5];
  Double_t x02 = par[6];
  Double_t x03 = par[7];
  Double_t x10 = par[8];
  Double_t x11 = par[9];
  Double_t x12 = par[10];
  Double_t x13 = par[11];
  Double_t x20 = par[12];
  Double_t x21 = par[13];
  Double_t x22 = par[14];
  Double_t x23 = par[15];
  Double_t x30 = par[16];
  Double_t x31 = par[17];
  Double_t x32 = par[18];
  Double_t x33 = par[19];
 
  Double_t zx0 = 1.+ROOT::Math::erf((x[0]-p1)/p2); 
  Double_t zx1 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*(-1./p2);
  Double_t zx2 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*((p1-x[0])/p2**2);
  Double_t zx3 = 1.;

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ zx0*zx3*x03+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ zx1*zx3*x13+
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22+ zx2*zx3*x23+  
                    zx3*zx0*x30 + zx3*zx1*x31 + zx3*zx2*x32+ zx3*zx3*x33;
 
  return par[0]*(1.+ROOT::Math::erf((x[0]-par[1])/par[2]))+par[3] + sqrt(result);
}


Double_t errorFitDownErf(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t x00 = par[3];
  Double_t x01 = par[4];
  Double_t x02 = par[5];
  Double_t x10 = par[6];
  Double_t x11 = par[7];
  Double_t x12 = par[8];
  Double_t x20 = par[9];
  Double_t x21 = par[10];
  Double_t x22 = par[11];

  Double_t zx0 = 1.+ROOT::Math::erf((x[0]-p1)/p2); 
  Double_t zx1 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*(-1./p2);
  Double_t zx2 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*((p1-x[0])/p2**2);

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ 
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22; 
  return par[0]*(1.+ROOT::Math::erf((x[0]-par[1])/par[2])) - sqrt(result);
}

Double_t errorFitDownErf4(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t p3 = par[3];
  Double_t x00 = par[4];
  Double_t x01 = par[5];
  Double_t x02 = par[6];
  Double_t x03 = par[7];
  Double_t x10 = par[8];
  Double_t x11 = par[9];
  Double_t x12 = par[10];
  Double_t x13 = par[11];
  Double_t x20 = par[12];
  Double_t x21 = par[13];
  Double_t x22 = par[14];
  Double_t x23 = par[15];
  Double_t x30 = par[16];
  Double_t x31 = par[17];
  Double_t x32 = par[18];
  Double_t x33 = par[19];
 
  Double_t zx0 = 1.+ROOT::Math::erf((x[0]-p1)/p2); 
  Double_t zx1 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*(-1./p2);
  Double_t zx2 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*((p1-x[0])/p2**2);
  Double_t zx3 = 1.;

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ zx0*zx3*x03+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ zx1*zx3*x13+
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22+ zx2*zx3*x23+  
                    zx3*zx0*x30 + zx3*zx1*x31 + zx3*zx2*x32+ zx3*zx3*x33;
 
  return par[0]*(1.+ROOT::Math::erf((x[0]-par[1])/par[2]))+par[3] - sqrt(result);
}

Double_t errorFitErf(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t x00 = par[3];
  Double_t x01 = par[4];
  Double_t x02 = par[5];
  Double_t x10 = par[6];
  Double_t x11 = par[7];
  Double_t x12 = par[8];
  Double_t x20 = par[9];
  Double_t x21 = par[10];
  Double_t x22 = par[11];

  Double_t zx0 = 1.+ROOT::Math::erf((x[0]-p1)/p2); 
  Double_t zx1 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*(-1./p2);
  Double_t zx2 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*((p1-x[0])/p2**2);

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ 
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22; 
  return sqrt(result);
}

Double_t errorFitErfDouble(Double_t x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t x00 = par[3];
  Double_t x01 = par[4];
  Double_t x02 = par[5];
  Double_t x10 = par[6];
  Double_t x11 = par[7];
  Double_t x12 = par[8];
  Double_t x20 = par[9];
  Double_t x21 = par[10];
  Double_t x22 = par[11];

  Double_t zx0 = 1.+ROOT::Math::erf((x-p1)/p2); 
  Double_t zx1 = p0*(2./sqrt(PI))*exp((-1.)*( (x-p1) /p2 )**2)*(-1./p2);
  Double_t zx2 = p0*(2./sqrt(PI))*exp((-1.)*( (x-p1) /p2 )**2)*((p1-x)/p2**2);

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ 
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22; 
  return sqrt(result);
}



Double_t errorFitErf4(Double_t *x,Double_t *par)
{
  Double_t p0 = par[0];
  Double_t p1 = par[1];
  Double_t p2 = par[2];
  Double_t p3 = par[3];
  Double_t x00 = par[4];
  Double_t x01 = par[5];
  Double_t x02 = par[6];
  Double_t x03 = par[7];
  Double_t x10 = par[8];
  Double_t x11 = par[9];
  Double_t x12 = par[10];
  Double_t x13 = par[11];
  Double_t x20 = par[12];
  Double_t x21 = par[13];
  Double_t x22 = par[14];
  Double_t x23 = par[15];
  Double_t x30 = par[16];
  Double_t x31 = par[17];
  Double_t x32 = par[18];
  Double_t x33 = par[19];
 
  Double_t zx0 = 1.+ROOT::Math::erf((x[0]-p1)/p2); 
  Double_t zx1 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*(-1./p2);
  Double_t zx2 = p0*(2./sqrt(PI))*exp((-1.)*( (x[0]-p1) /p2 )**2)*((p1-x[0])/p2**2);
  Double_t zx3 = 1.;

  Double_t result = zx0*zx0*x00 + zx0*zx1*x01 + zx0*zx2*x02+ zx0*zx3*x03+ 
                    zx1*zx0*x10 + zx1*zx1*x11 + zx1*zx2*x12+ zx1*zx3*x13+
                    zx2*zx0*x20 + zx2*zx1*x21 + zx2*zx2*x22+ zx2*zx3*x23+  
                    zx3*zx0*x30 + zx3*zx1*x31 + zx3*zx2*x32+ zx3*zx3*x33;
 
  return sqrt(result);
}

Double_t sfFinal(Double_t *x, Double_t *par)
{
  Double_t par_data[3]= {par[0], par[1], par[2]};
  Double_t par_mc[3]= {par[12], par[13], par[14]};

  return pilar(x,par_data) / pilar(x,par_mc);
}

Double_t sfFinal4(Double_t *x, Double_t *par)
{
  Double_t par_data[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc[4]= {par[20], par[21], par[22], par[23]};

  return pilar(x,par_data) / pilar(x,par_mc);
}

Double_t sfFinalUp(Double_t *x, Double_t *par)
{
  Double_t par12_data[12]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11]};
  Double_t par12_mc[12]= {par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19], par[20], par[21], par[22], par[23]};

  Double_t par_data[3]= {par[0], par[1], par[2]};
  Double_t par_mc[3]= {par[12], par[13], par[14]};

  Double_t scaleFactor = pilar(x,par_data) / pilar(x,par_mc);
 
  Double_t error_data = errorFit(x,par12_data); 
  Double_t error_mc = errorFit(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilar(x,par_data))**2 + (error_mc/pilar(x,par_mc))**2 );

  return scaleFactor + error;
}
Double_t sfFinalUp4(Double_t *x, Double_t *par)
{
  Double_t par12_data[20]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11],
                            par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19]};
  Double_t par12_mc[20]= {par[20], par[21], par[22], par[23], par[24], par[25], par[26], par[27], par[28], par[29], par[30], par[31],
                            par[32], par[33], par[34], par[35], par[36], par[37], par[38], par[39]};

  Double_t par_data[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc[4]= {par[20], par[21], par[22], par[23]};

  Double_t scaleFactor = pilar(x,par_data) / pilar(x,par_mc);
 
  Double_t error_data = errorFit4(x,par12_data); 
  Double_t error_mc = errorFit4(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilar(x,par_data))**2 + (error_mc/pilar(x,par_mc))**2 );

  return scaleFactor + error;
}

Double_t sfFinalDown(Double_t *x, Double_t *par)
{
  Double_t par12_data[12]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11]};
  Double_t par12_mc[12]= {par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19], par[20], par[21], par[22], par[23]};

  Double_t par_data[3]= {par[0], par[1], par[2]};
  Double_t par_mc[3]= {par[12], par[13], par[14]};
  Double_t scaleFactor = pilar(x,par_data) / pilar(x,par_mc);
 
  Double_t error_data = errorFit(x,par12_data); 
  Double_t error_mc = errorFit(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilar(x,par_data))**2 + (error_mc/pilar(x,par_mc))**2 );

  return scaleFactor - error;
}

Double_t sfFinalDown4(Double_t *x, Double_t *par)
{
  Double_t par12_data[20]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11],
                            par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19]};
  Double_t par12_mc[20]= {par[20], par[21], par[22], par[23], par[24], par[25], par[26], par[27], par[28], par[29], par[30], par[31],
                            par[32], par[33], par[34], par[35], par[36], par[37], par[38], par[39]};

  Double_t par_data[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc[4]= {par[20], par[21], par[22], par[23]};

  Double_t scaleFactor = pilar(x,par_data) / pilar(x,par_mc);
 
  Double_t error_data = errorFit4(x,par12_data); 
  Double_t error_mc = errorFit4(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilar(x,par_data))**2 + (error_mc/pilar(x,par_mc))**2 );

  return scaleFactor - error;
}


Double_t sfFinalErf(Double_t *x, Double_t *par)
{
  Double_t par_data[3]= {par[0], par[1], par[2]};
  Double_t par_mc[3]= {par[12], par[13], par[14]};

  return pilarErf(x,par_data) / pilarErf(x,par_mc);
}

Double_t sfFinalErf4(Double_t *x, Double_t *par)
{
  Double_t par_data[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc[4]= {par[20], par[21], par[22], par[23]};

  return pilarErf(x,par_data) / pilarErf(x,par_mc);
}

Double_t sfFinalErf4Simplified(Double_t *x, Double_t *par)
{
  Double_t par_data[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc[4]= {par[4], par[5], par[6], par[7]};

  return pilarErf(x,par_data) / pilarErf(x,par_mc);
}

Double_t sfFinalUpErf(Double_t *x, Double_t *par)
{
  Double_t par12_data[12]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11]};
  Double_t par12_mc[12]= {par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19], par[20], par[21], par[22], par[23]};

  Double_t par_data[3]= {par[0], par[1], par[2]};
  Double_t par_mc[3]= {par[12], par[13], par[14]};
  Double_t scaleFactor = pilarErf(x,par_data) / pilarErf(x,par_mc);
 
  Double_t error_data = errorFitErf(x,par12_data); 
  Double_t error_mc = errorFitErf(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilarErf(x,par_data))**2 + (error_mc/pilarErf(x,par_mc))**2 );

  return scaleFactor + error;
}

Double_t sfFinalUpErf4(Double_t *x, Double_t *par)
{
  Double_t par12_data[20]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11],
                            par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19]};
  Double_t par12_mc[20]= {par[20], par[21], par[22], par[23], par[24], par[25], par[26], par[27], par[28], par[29], par[30], par[31],
                            par[32], par[33], par[34], par[35], par[36], par[37], par[38], par[39]};

  Double_t par_data[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc[4]= {par[20], par[21], par[22], par[23]};


  Double_t scaleFactor = pilarErf(x,par_data) / pilarErf(x,par_mc);
 
  Double_t error_data = errorFitErf4(x,par12_data); 
  Double_t error_mc = errorFitErf4(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilarErf(x,par_data))**2 + (error_mc/pilarErf(x,par_mc))**2 );

  return scaleFactor + error;
}

Double_t sfFinalDownErf(Double_t *x, Double_t *par)
{
  Double_t par12_data[12]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11]};
  Double_t par12_mc[12]= {par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19], par[20], par[21], par[22], par[23]};

  Double_t par_data[3]= {par[0], par[1], par[2]};
  Double_t par_mc[3]= {par[12], par[13], par[14]};
  Double_t scaleFactor = pilarErf(x,par_data) / pilarErf(x,par_mc);
 
  Double_t error_data = errorFitErf(x,par12_data); 
  Double_t error_mc = errorFitErf(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilarErf(x,par_data))**2 + (error_mc/pilarErf(x,par_mc))**2 );

  return scaleFactor - error;
}

Double_t sfFinalDownErf4(Double_t *x, Double_t *par)
{
  Double_t par12_data[20]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11],
                            par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19]};
  Double_t par12_mc[20]= {par[20], par[21], par[22], par[23], par[24], par[25], par[26], par[27], par[28], par[29], par[30], par[31],
                            par[32], par[33], par[34], par[35], par[36], par[37], par[38], par[39]};

  Double_t par_data[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc[4]= {par[20], par[21], par[22], par[23]};


  Double_t scaleFactor = pilarErf(x,par_data) / pilarErf(x,par_mc);
 
  Double_t error_data = errorFitErf4(x,par12_data); 
  Double_t error_mc = errorFitErf4(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilarErf(x,par_data))**2 + (error_mc/pilarErf(x,par_mc))**2 );

  return scaleFactor - error;
}

Double_t sfFinalRatio4(Double_t *x, Double_t *par)
{
  Double_t par_data1[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc1[4]= {par[4], par[5], par[6], par[7]};
  Double_t par_data2[4]= {par[8], par[9], par[10], par[11]};
  Double_t par_mc2[4]= {par[12], par[13], par[14], par[15]};


  Double_t scaleFactor1 = pilarErf(x,par_data1) / pilarErf(x,par_mc1);
  Double_t scaleFactor2 = pilarErf(x,par_data2) / pilarErf(x,par_mc2);
 
  return scaleFactor1/scaleFactor2;

}
Double_t efficFinalRatio4(Double_t *x, Double_t *par)
{
  Double_t par_data1[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_data2[4]= {par[4], par[5], par[6], par[7]};

  Double_t result = pilarErf(x,par_data1) / pilarErf(x,par_data2);
 
  return result;

}

Double_t sfFinalErrorErf4(Double_t *x, Double_t *par)
{
  Double_t par12_data[20]= {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], par[9], par[10], par[11],
                            par[12], par[13], par[14], par[15], par[16], par[17], par[18], par[19]};
  Double_t par12_mc[20]= {par[20], par[21], par[22], par[23], par[24], par[25], par[26], par[27], par[28], par[29], par[30], par[31],
                            par[32], par[33], par[34], par[35], par[36], par[37], par[38], par[39]};

  Double_t par_data[4]= {par[0], par[1], par[2], par[3]};
  Double_t par_mc[4]= {par[20], par[21], par[22], par[23]};


  Double_t scaleFactor = pilarErf(x,par_data) / pilarErf(x,par_mc);
 
  Double_t error_data = errorFitErf4(x,par12_data); 
  Double_t error_mc = errorFitErf4(x,par12_mc); 

  Double_t error = scaleFactor*sqrt( (error_data/pilarErf(x,par_data))**2 + (error_mc/pilarErf(x,par_mc))**2 );

  return error;
}