UnfoldUtilities.icc

Go to the documentation of this file.

/*  Utilities to  unfold the vertical spectrum such as done for ICRC 2015  */

/* Author: I. Valiño, based on functions provided by A. Schulz */


TVectorD GetSDCalPars(bool infill = false)
{
  //Calibration constants
  TVectorD pCal(2);

  if (infill) {
    if (HeraldAnalysis) {
      pCal[0] = 1.62710e+16;
      pCal[1] = 0.9665210;
    } else {
      pCal[0] = 1.40738e+16;
      pCal[1] = 0.9999375;
    }
  } else {
    if (HeraldAnalysis) {
      pCal[0] = 2.23584e+17;
      pCal[1] = 0.984078;
    } else {
      pCal[0] = 1.78368e+17; 
      pCal[1] = 1.042043;
    }
  }

  return pCal;
}


double InverseSdCalibration(const double lgE, TVectorD pCal)
{
  return (lgE - log10(pCal[0])) / pCal[1]; /*  returns lgS38 or lgS35 */
}


double ScaledErrorFunction(const double x, const std::vector<double>& pars)
{
  return 0.5 * (1 + TMath::Erf((x - pars[0]) / pars[1]));
}


double SD1500TriggerEfficiency(const double lgE, const double zenith)
{

  const TVectorD pCal = GetSDCalPars(false);

  const double pA[4] = {1.18586522, -2.5940793, 3.25692833, -1.19604131};
  const double cos2th = pow(TMath::Cos(zenith), 2.);

  const double A = pA[0] + cos2th * pA[1] + cos2th * cos2th * pA[2] + cos2th * cos2th * cos2th * pA[3];
  const double B = 0.36868343;

  const double effpars[] = {A, B};
  const std::vector<double> parvec(effpars, effpars + 2);

  const double lgS38 = InverseSdCalibration(lgE, pCal);
  const double feff = ScaledErrorFunction(lgS38, parvec);

  return feff;
}


double SD750TriggerEfficiency(const double lgE, const double zenith)
{
  // From simulated Proton showers

  TVectorD pCal = GetSDCalPars(true);

  const double pA[4] = {2.38516877, -4.85552037, 4.10415064, -0.97515248};
  const double cos2th = pow(TMath::Cos(zenith), 2.);

  const double A = pA[0] + cos2th * pA[1] + cos2th * cos2th * pA[2] + cos2th * cos2th * cos2th * pA[3];
  const double B = 0.24852449;

  const double effpars[] = {A, B};
  const std::vector<double> parvec(effpars, effpars + 2);

  const double lgS35 = InverseSdCalibration(lgE, pCal);
  const double feff = ScaledErrorFunction(lgS35, parvec);

  return feff;
}


double EfficientZenithDistrib(double* x, double* par)
{

  const double zenith = x[0];
  const double lgE[1] = {par[0]};
  const bool infill = par[1];

  double eff;

  if (infill)
    eff = SD750TriggerEfficiency(lgE[0], x[0]);
  else
    eff = SD1500TriggerEfficiency(lgE[0], x[0]);

  return eff * TMath::Cos(zenith) * TMath::Sin(zenith);

}


double SDTriggerEfficiency(const double lgE, const bool infill)
{
  const double thmax = (infill ? 55. : 60.) * TMath::DegToRad();
  const double thmin = 0. * TMath::DegToRad();

  TF1 feff("feff", EfficientZenithDistrib, thmin, thmax, 2);
  feff.FixParameter(0, lgE);
  feff.FixParameter(1, infill);

  ROOT::Math::IntegratorOneDimOptions::SetDefaultAbsTolerance(1.E-6);
  ROOT::Math::IntegratorOneDimOptions::SetDefaultRelTolerance(1.E-6);

  ROOT::Math::WrappedTF1 weff(feff);
  ROOT::Math::Integrator ig(ROOT::Math::IntegrationOneDim::kADAPTIVE);
  ig.SetFunction(weff);

  double value_integral = ig.Integral(thmin, thmax);
  double fnorm = 0.5 * (pow(TMath::Sin(thmax), 2) - pow(TMath::Sin(thmin), 2));

  return value_integral / fnorm;
}


TVectorD EnergyBias(TVectorD lgEs, TVectorD pCal, const bool infill)
{

  //const double pars[3] = {0.12287307, -0.10904358, 0.02411047};

  int ncol = lgEs.GetNoElements();
  TVectorD biasEnergy(ncol);

  for (int i = 0; i < ncol; ++i) {
    if (infill) {
      biasEnergy[i] = 1.; //No bias for infill
    } else {
      // double lgS38 = InverseSdCalibration(lgEs[i], pCal);
      // double biasS38 = 1 + pars[0] + pars[1] * lgS38 + pars[2] * lgS38 * lgS38;
      // biasEnergy[i] = pow(biasS38, pCal[1]);
      biasEnergy[i] = 1.; //No bias for SD
    }
  }

  return biasEnergy;
}



double SimulationSD1500EnergyResolution(const double lgE, const double zenith)
{
  /* obtained from 50-50 p/Fe simulations using QGSJetII-03 */

  double p0 =  1.090032150314515363e-01;
  double p1 =  4.353891452953017049e-01;

  double Esd = pow(10., lgE);

  double relErrs = p0 + p1 / sqrt(Esd / 1e+17) ;

  return relErrs;
}

double SimulationSD750EnergyResolution(const double lgE, const double zenith)
{
  //new infill resolution model from Daniela
  //determined together with fd energy bias correction and new invisible energy model
  //QGSJET 2.4 simulations, quality + FoV cuts, resolution determined for p and Fe

  const double pars_p[2] = {1.078638367648934449e-01, 1.595403269466228735e-01};
  const double pars_fe[2] = {4.836779040905066218e-02, 1.697359371453481258e-01};
  const double pars[2] = {(pars_p[0] + pars_fe[0]) / 2., (pars_p[1] + pars_fe[1]) / 2.};

  return pars[0] + pars[1] / sqrt(pow(10, lgE) / 1e17);
}

double DataSD1500EnergyResolution(const double lgE, const double zenith)
{
  /* obtained from data itself */
  TVectorD pCal = GetSDCalPars(false);
  const double lgS38 = InverseSdCalibration(lgE, pCal);

  const double pars[9] = {7.964648011936218460e-04,
                          4.118659698340852438e-01,
                          -1.517612358021901198e+00,
                          9.102368706313619384e+00,
                          -9.938644262370022631e+00,
                          4.120050818744602217e-01,
                          1.244124736151181560e-03,
                          -6.535475475395227107e+00,
                          1.039836139398396853e+01
                         };
  const double lgsectheta = TMath::Log10(1 / cos(zenith));
  double res = pars[0] + (pars[1] + pars[2] * lgsectheta + pars[3] * lgsectheta * lgsectheta + pars[4] * lgsectheta * lgsectheta * lgsectheta) * pow(10., -lgS38 * 0.5)
               + (pars[5] + pars[6] * lgsectheta + pars[7] * lgsectheta * lgsectheta + pars[8] * lgsectheta * lgsectheta * lgsectheta) * pow(10., -lgS38);

  return res * pCal[1];
}


double DataSD750EnergyResolution(const double lgE, const double zenith)
{
  TVectorD pCal = GetSDCalPars(true);
  const double lgS35 = InverseSdCalibration(lgE, pCal);

  const double pars[6] = {1.713893029281758462e+00, -1.115258876619560580e+00, 1.557604707934408061e+01,
                          6.861112716697782554e-01, -8.015754473911778089e-01, 8.852556156454199909e-01
                         };
  const double lgsectheta = TMath::Log10(1 / cos(zenith));

  return (pars[0] + pars[1] * lgsectheta + pars[2] * lgsectheta * lgsectheta)
         * (pow(10., (-pars[3] * lgS35)) + pars[4] * pow(10., (-pars[5] * lgS35)));
}


double ResolutionZenithDistrib(double* x, double* par)
{

  const double zenith = x[0];
  const double lgE = par[0];
  const bool data = bool(par[1]);
  const bool infill = bool(par[2]);

  double res;
  if (data) {
    if (infill)
      res = DataSD750EnergyResolution(lgE, x[0]);
    else
      res = DataSD1500EnergyResolution(lgE, x[0]);
  } else {
    if (infill)
      res = SimulationSD750EnergyResolution(lgE, x[0]);
    else
      res = SimulationSD1500EnergyResolution(lgE, x[0]);
  }
  return res * TMath::Cos(zenith) * TMath::Sin(zenith);
}


double IntegrateEnergyResolution(const double lgE, const bool data, const bool infill)
{
  const double thmax = (infill ? 55. : 60.) * TMath::DegToRad();
  const double thmin = 0. * TMath::DegToRad();

  TF1 fres("fres", ResolutionZenithDistrib, thmin, thmax, 3);
  fres.FixParameter(0, lgE);
  fres.FixParameter(1, 1. * data);
  fres.FixParameter(2, 1. * infill);

  ROOT::Math::IntegratorOneDimOptions::SetDefaultAbsTolerance(1.E-6);
  ROOT::Math::IntegratorOneDimOptions::SetDefaultRelTolerance(1.E-6);

  ROOT::Math::WrappedTF1 wres(fres);
  ROOT::Math::Integrator ig(ROOT::Math::IntegrationOneDim::kADAPTIVE);
  ig.SetFunction(wres);

  double value_integral = ig.Integral(thmin, thmax);
  double fnorm = 0.5 * (pow(TMath::Sin(thmax), 2) - pow(TMath::Sin(thmin), 2));

  return value_integral / fnorm;
}


double CustomShowerToShowerFluctuation(const double lgE)
{
  const double xmin = 16.5;
  const double xmax = 20.5;
  const double z = (lgE - xmin) / (xmax - xmin);

  const double pars[2] = {0.2, 0.0783};
  double shsh = (1 - z) * pars[0] + z * pars[1];
  return shsh;
}


// TODO
double CustomInfillShowerToShowerFluctuation(const double lgE)
{
  double shsh = 0.13;
  return shsh;
}


TVectorD EnergyResolution(TVectorD lgEs, const bool data, const bool infill)
{
  int ncol = lgEs.GetNoElements();
  TVectorD EnergyRes(ncol);

  for (int i = 0; i < ncol; ++i) {
    double lgE = lgEs[i];
    double Esd = pow(10., lgE);

    double relErrs1 = IntegrateEnergyResolution(lgE, data, infill);
    double relErrs2 = 0;

    if (data) {
      if (infill)
        relErrs2 = CustomInfillShowerToShowerFluctuation(lgE);
      else
        relErrs2 = CustomShowerToShowerFluctuation(lgE);
    }

    double relErrs = TMath::Sqrt(pow(relErrs1, 2)  + pow(relErrs2, 2));

    EnergyRes[i] = relErrs * Esd;
  }

  return EnergyRes;
}


TMatrixD kResolutionMatrix(TVectorD pCal, TVectorD xos, double* resolutionOffset, const bool useResFromData, const bool infill)
{

  int N = xos.GetNoElements();

  TVectorD ses = EnergyResolution(xos, useResFromData, infill);

  TVectorD Ebias = EnergyBias(xos, pCal, infill);

  TMatrixD tmp_kmatrix(N, N);
  for (int i = 0; i < N; ++i) {

    double ey = pow(10, xos[i]);
    for (int j = 0; j < N; ++j) {

      double ex = pow(10, xos[j]);
      tmp_kmatrix[i][j] = TMath::Gaus(ey, ex * Ebias[j], ses[j], true) * ey * log(10.) * SDTriggerEfficiency(xos[j], infill);
    }
  }
  return tmp_kmatrix;
}