//**********************************************
Double_t cteq61(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //  {1.56811e+10,-1.94498e+01, 2.87536e+01,-2.13869e+01, 7.04488e+00,-8.73467e-01, 4.20494e+00};
    //  {1.48779e+10, -1.94296e+01,  2.87521e+01, -2.13882e+01,  7.04434e+00, -8.73215e-01, 4.21074e+00};
  {1.22024e+13, -1.74694e+01, 1.82338e+01, -1.22027e+01, 3.86174e+00, -4.74714e-01, 5.10307e+00};

  Double_t val;
  val = p[0]*TMath::Power(10,(p[1]*TMath::log10(*x)+
			      p[2]*TMath::Power(TMath::log10(*x),2)+p[3]*TMath::Power(TMath::log10(*x),3)+
			      p[4]*TMath::Power(TMath::log10(*x),4)+p[5]*TMath::Power(TMath::log10(*x),5)))*
    TMath::Power((1-*x/900),p[6]);
  return val;

}


//***********
Double_t physR2Sys(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.50092E+07, -4.3622,  5.5546, -5.0231,  1.6308, -0.1875,  6.2428};
    {0.80608E+07, -4.6120,  5.5936, -4.9981,  1.6339, -0.1923,  5.5157};


  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}



//***********
Double_t sysEScaleDiffP(Double_t *x, Double_t *par) {
  Double_t val;
  // 5% escale error
  //  val = 2.18793e-01+3.30016e-04*(*x)+9.59740e-07*(*x)*(*x);
  // 3% escale error
  val = 1.42430e-01+4.64936e-06*(*x)+1.01157e-06*(*x)*(*x);
  
  return val;
}



//***********
Double_t sysEScaleDiffN(Double_t *x, Double_t *par) {
  Double_t val;
  // 5% escale error
  //  val = -1.89764e-01-3.29895e-04*(*x)-3.34663e-07*(*x)*(*x); 
  // 3% escale error
  val = -1.26656e-01-6.47329e-05*(*x)-6.80649e-07*(*x)*(*x); 


  return val;
}




//***********
Double_t sysUEP(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.40736E+07, -4.9159,  5.6363, -4.3766,  1.1988, -0.1038,  9.4034};
    //    {0.10346E+07, -3.8331,  5.5104, -5.0512,  1.6250, -0.1851,  5.9659};
    {0.11713E+07, -3.8952,  5.5194, -5.0456,  1.6259, -0.1861,  5.8748};
  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t sysUEN(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.47090E+06, -2.3809,  4.4818, -5.0336,  1.7747, -0.2141,  7.1747};
    //    {0.30257E+08, -4.9418,  5.4885, -4.9359,  1.6506, -0.2002,  5.1818};
    {0.33532E+08, -4.9488,  5.4396, -4.9170,  1.6568, -0.2033,  4.9526};
  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t sysUEDiffN(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (sysUEN(x, par) - physR2Sys(x, par))/physR2Sys(x, par);
  val = (sysUEP(x, par) - sysUEN(x, par))/2/physR2Sys(x, par);
  return val;
}

//***********
Double_t sysUEDiffP(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (sysUEP(x, par) - physR2Sys(x, par))/physR2Sys(x, par);
  val = (sysUEN(x, par) - sysUEP(x, par))/2/physR2Sys(x, par);
  return val;
}



//***********
Double_t sysFragP(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.36132E+06, -2.4681,  4.5247, -5.0151,  1.7763, -0.2183,  6.7994};
    //    {0.99804E+07, -4.7123,  5.5894, -4.9887,  1.6374, -0.1927,  6.1616};
    {0.11444E+08, -4.7721,  5.5890, -4.9802,  1.6394, -0.1941,  6.0315};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t sysFragN(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.94674E+07, -4.9981,  5.5442, -4.4042,  1.2491, -0.1164,  8.4303};
    //    {0.62701E+07, -4.4762,  5.5681, -5.0094,  1.6340, -0.1902,  5.8795};
    {0.15829E+08, -4.0946,  3.4749, -2.8896,  0.7780, -0.0645,  7.3061};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t sysFragDiffN(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - sysFragN(x, par))/physR2Sys(x, par);
  val = (sysFragP(x, par)-sysFragN(x, par))/2/physR2Sys(x, par);
  return val;
}

//***********
Double_t sysFragDiffP(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - sysFragP(x, par))/physR2Sys(x, par);
  val = (sysFragN(x, par)-sysFragP(x, par))/2/physR2Sys(x, par);
  return val;
}



//**************
Double_t syselphP(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.11070E+08, -4.8078,  5.0650, -4.0845,  1.1628, -0.1072,  8.9551};
    //    {0.60349E+07, -4.4597,  5.5676, -5.0108,  1.6332, -0.1905,  5.6571};
    {0.61652E+07, -4.4681,  5.5674, -5.0104,  1.6334, -0.1904,  5.7447};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t syselphN(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.14922E+08, -4.5909,  4.2562, -3.3789,  0.9083, -0.0730,  9.4803};
    //    {0.61669E+07, -4.4798,  5.5667, -5.0101,  1.6337, -0.1903,  5.9263};
    {0.87191E+07, -4.6552,  5.5960, -4.9936,  1.6348, -0.1930,  5.5499};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t syselphDiffN(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - syselphN(x, par))/physR2Sys(x, par);
  val = (syselphN(x, par) - syselphP(x, par))/2/physR2Sys(x, par);
  return val;
}

//***********
Double_t syselphDiffP(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - syselphP(x, par))/physR2Sys(x, par);
  val = (syselphP(x, par) - syselphN(x, par))/2/physR2Sys(x, par);
  return val;
}

//**************
Double_t sysloptP(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.23461E+06, -2.2652,  4.4956, -5.0331,  1.7737, -0.2149,  7.1520};
    //    {0.10199E+08, -4.8187,  5.6928, -4.9927,  1.6143, -0.1861,  6.4173};
    {0.97543E+07, -4.7096,  5.5861, -4.9859,  1.6378, -0.1931,  6.0023};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t sysloptN(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.19367E+06, -2.1507,  4.4985, -5.0392,  1.7697, -0.2167,  5.7238};
    //    {0.55736E+07, -4.4132,  5.5681, -5.0135,  1.6315, -0.1913,  5.0920};
    {0.56948E+07, -4.4217,  5.5679, -5.0132,  1.6317, -0.1912,  5.1780};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t sysloptDiffN(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - sysloptN(x, par))/physR2Sys(x, par);
  val = (sysloptP(x, par) - sysloptN(x, par))/2/physR2Sys(x, par);
  return val;
}

//***********
Double_t sysloptDiffP(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - sysloptP(x, par))/physR2Sys(x, par);
  val = (sysloptN(x, par) - sysloptP(x, par))/2/physR2Sys(x, par);
  return val;
}


//**************
Double_t syshiptP(Double_t *x, Double_t *par) {
  Double_t p[] =
    //    {0.19652E+06, -2.1737,  4.4959, -5.0392,  1.7699, -0.2166,  6.0041};
    //    {0.55854E+07, -4.4333,  5.5660, -5.0137,  1.6315, -0.1913,  5.3019};
    {0.24919E+08, -3.7844,  2.3407, -1.9354,  0.4451, -0.0217,  7.8388};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t syshiptN(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.21768E+06, -2.2124,  4.4980, -5.0352,  1.7722, -0.2156,  6.4929};
    //    {0.62701E+07, -4.4762,  5.5681, -5.0094,  1.6340, -0.1902,  5.8795};
    {0.80306E+07, -4.6029,  5.5917, -4.9988,  1.6343, -0.1918,  5.6722};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t syshiptDiffN(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - syshiptN(x, par))/physR2Sys(x, par);
  val = (syshiptP(x, par) - syshiptN(x, par))/2/physR2Sys(x, par);
  return val;
}

//***********
Double_t syshiptDiffP(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - syshiptP(x, par))/physR2Sys(x, par);
  val = (syshiptN(x, par) - syshiptP(x, par))/2/physR2Sys(x, par);
  return val;
}



//**************
Double_t sysResP(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.70104E+07, -4.9984,  5.7187, -4.5947,  1.3312, -0.1283,  8.9305};
    //    {0.31863E+07, -4.2462,  5.5565, -5.0271,  1.6291, -0.1879,  5.9871};
    {0.36615E+07, -4.3133,  5.5638, -5.0203,  1.6304, -0.1892,  5.8693};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t sysResN(Double_t *x, Double_t *par) {
  Double_t p[] = 
    //    {0.40779E+06, -2.4523,  4.5214, -5.0179,  1.7750, -0.2187,  6.0437};
    //    {0.39754E+08, -3.6642,  1.9028, -1.6664,  0.3873, -0.0175,  7.9398};
    {0.15407E+08, -4.8184,  5.5762, -4.9707,  1.6399, -0.1961,  5.1748};

  Double_t ecm = 1960.;
  Double_t val;
  Double_t L10Et,pol,log_term,xt,xt_term,pt_term;

  L10Et =TMath::log10(*x);

  pol = p[1]*L10Et
    + p[2]*L10Et**2
    + p[3]*L10Et**3
    + p[4]*L10Et**4
    + p[5]*L10Et**5;

  if (fabs(pol)>50.0) {
    if(pol<0) {
      pol=-50.0;
    } else {
      pol=+50.0;
    }
  }

  log_term = TMath::Power(10,pol);
  xt = 2.0*(*x)/ecm;
  xt_term = TMath::Power((1-xt),p[6]);
  pt_term = 1.0;
  val = p[0]*pt_term*xt_term*log_term;
  return val;

}


//***********
Double_t sysResDiffN(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - sysResN(x, par))/physR2Sys(x, par);
  val = (sysResP(x, par) - sysResN(x, par))/2/physR2Sys(x, par);
  return val;
}

//***********
Double_t sysResDiffP(Double_t *x, Double_t *par) {
  Double_t val;
  //  val = (physR2Sys(x, par) - sysResP(x, par))/physR2Sys(x, par);
  val = (sysResN(x, par) - sysResP(x, par))/2/physR2Sys(x, par);
  return val;
}

//***********
Double_t sysTotDiffN(Double_t *x, Double_t *par) {
  Double_t val;
  val = (TMath::Power(sysEScaleDiffN(x,par),2) +
	 TMath::Power(syshiptDiffN(x,par),2) +
	 TMath::Power(sysloptDiffN(x,par),2) +
	 TMath::Power(sysUEDiffN(x,par),2) +
	 TMath::Power(syselphDiffN(x,par),2) +
	 TMath::Power(sysFragDiffN(x,par),2) +
	 TMath::Power(sysResDiffN(x,par),2));
  
  val = -TMath::Power(val,0.5);
  //  val = sysEScaleDiffN(x,par);
  return val;
}

//***********
Double_t sysTotDiffP(Double_t *x, Double_t *par) {
  Double_t val;
  val = (TMath::Power(sysEScaleDiffP(x,par),2) +
	 TMath::Power(syshiptDiffP(x,par),2) +
	 TMath::Power(sysloptDiffP(x,par),2) +
	 TMath::Power(sysUEDiffP(x,par),2) +
	 TMath::Power(syselphDiffP(x,par),2) +
	 TMath::Power(sysFragDiffP(x,par),2) +
	 TMath::Power(sysResDiffP(x,par),2));

  //  val = sysEScaleDiffP(x,par);
  val = TMath::Power(val,0.5);

  return val;
}