FitterUtil.cc

Go to the documentation of this file.

#include <FitterUtil.h>
#include <TVector3.h>
#include <cmath>
#include <iostream>

#include <tls/UnivTimeKG.h>

using namespace std;

namespace FitterUtil {

  double getDistToFirstInteraction(double hfi, double sHeight, double zenith, double dist, double psi)
  {
    //dist, psi given in SP coordinates. This affects the x coordinate: x_G = x_SP/cos(zenith)
    TVector3 sPosition(dist * cos(psi) / cos(zenith), dist * sin(psi), sHeight);
    TVector3 core(0., 0., sHeight) ;
    TVector3 axis(sin(zenith), 0., cos(zenith));
    double d_core_fi = (hfi - sHeight) / cos(zenith);
    TVector3 fiPosition = core + axis * d_core_fi;
    double d_station_fi = (sPosition - fiPosition).Mag();
    return d_station_fi;
  }

  // compute the time difference between the arrival of the plane and curved shower front (see GAP 2013-105, fig. 4.17)
  double getFirstParticleArrivalTime(double r, double dx0)
  {
    //computes the delay referred to the plane front time
    //dist is given in SP coordinates.
    double a = 1. / (2 * utl::kSpeedOfLight); // just using a sphere that travels with the speed of light
    return a * (r * r / dx0);
  }

  // mean reconstructed Xmax bias = f(theta, log10(E))
  // derived from continuous library of QGSJetII-03 (mean of p+Fe)
  // CORSIKA showers from
  // http://augerdb.lal.in2p3.fr:8080/augerdb/simdb/Library-en.do?libraryId=11963710
  // http://augerdb.lal.in2p3.fr:8080/augerdb/simdb/Library-en.do?libraryId=9136660
  // see GAP 2013-105, fig. 6.6
  double getMeanXmaxBias(double loge, double theta, bool saturated)
  {
    double bias = 0.;
    if (saturated)
      bias = -1.655868e+03 + 8.171344e+01 * loge + 3.401878e+01 * theta;
    else
      bias = -9.512000e+02 + 4.997809e+01 * loge + -4.301163e+01 * theta;

    return bias * gpercm2;
  }



  // parametrization of the mean Xmax from simulations
  // mean of proton and iron, used as starting value for the Xmax(SD) reconstruction
  double getMeanParametrizedXmax(double loge)
  {

    double xmax_proton, xmax_iron;
    {
      // proton QGSJet II-03
      const double p0 = 7.36508298967330802e+02;
      const double p1 = 4.89094109886878101e+01;
      const double p2 = -2.27145404084728142e+00;
      xmax_proton = (p0 + p1 * (loge - 18) + p2 * pow(loge - 18, 2)) * gpercm2;
    }
    {
      // iron QGSJet II-03
      const double p0 = 6.45205996945714332e+02;
      const double p1 = 5.76560096429448592e+01;
      const double p2 = -3.60856107551782701e+00;
      xmax_iron = (p0 + p1 * (loge - 18) + p2 * pow(loge - 18, 2)) * gpercm2;
    }

    return (xmax_proton + xmax_iron) / 2.;
  }

  // fit to data from long Xmax paper 2014 (phys rev D 90, 122005 (2014))
  double getMeanParametrizedXmaxData(double loge)
  {
    double meanxmax;

    const double logE0 = 18.27, Xmax0 = 746.8, ER0 = 86.4, ER1 = 26.4;

    if (loge < logE0)
      meanxmax = Xmax0 + (loge - logE0) * ER0;
    else
      meanxmax = Xmax0 + (loge - logE0) * ER1;

    return meanxmax;
  }

  // parameterized with golden hybrids from 2004 - 2014
  double getMeanParametrizedNmuData(double loge, double theta)
  {
    // const double pars[3] = {2.18275629, 0.07113437, 0.37852463};
    const double pars[3] = {2.12606249, 0.09410723, 0.29323277};

    return pars[0] + pars[1] * (loge - 19.) + pars[2] * (1. / cos(theta) - 2);
  }

  double getMeanParametrizedNmuData(double loge, double theta, double Xmax)
  {
    // const double pars[2] = {2.29841005e-01, 7.57069565e+02};  // uncs: 0.35998355, 3.05089875
    // const double pars[2] = {9.43590516e-02, 7.56011955e+02};  // uncs: 0.35998355, 3.05089875
    const double pars[2] = {41.10517485, 756.92655277};  // uncs: 0.35974395, 3.07881652

    const double mval = getMeanParametrizedNmuData(loge, theta);
    const double xmaxval = (1 + 0.5 / M_PI * atan((pars[0] * loge - Xmax) / 40.)) /
                           (1 + 0.5 / M_PI * atan((pars[0] * loge - pars[1]) / 40.));
    return mval * xmaxval;
  }

  // adapted from Maximo
  double GetMeanMuonOffsetDataCalibrated(double r, double lgE, double theta)
  {
    // double OffsetM_Mu[3] = {0.06579588, 0.04701472, 0.11184198};  // uncs: 0.0467337, 0.03512474, 0.06394015
    double OffsetM_Mu[3] = {0.04308893, 0.0284907, 0.08383378};  // uncs: 0.0467337, 0.03512474, 0.06394015

    double fOffsetM_Mu = OffsetM_Mu[0] + OffsetM_Mu[1] * (lgE - 19.) + OffsetM_Mu[2] * (1 / cos(theta) - 2.);
    const double SysPar1 = 2.65;

    double C = fOffsetM_Mu / (1 - exp(-SysPar1));
    double mshift = C * (1 - exp(-SysPar1 * r / 2e5));

    return mshift;
  }

  // todo: add pars
  double
  GetMeanMuonOffsetMCCalibrated(double /*r*/, double /*lgE*/, double /*theta*/)
  {
    // update
    /*const double OffsetM_Mu[3] = { 0.11, 0., 0. };

    const double fOffsetM_Mu = OffsetM_Mu[0] + OffsetM_Mu[1] * (lgE - 19) + OffsetM_Mu[2] * (1 / cos(theta) - 2);
    const double SysPar1 = 2.65;

    const double C = fOffsetM_Mu / (1 - exp(-SysPar1));
    double mshift = C * (1 - exp(-SysPar1 * r / 2e5));*/

    return 0;
  }


  double GetMeanMuonOffset(double r, double offset0)
  {
    const double SysPar1 = 2.65;

    double C = offset0 / (1 - exp(-SysPar1));
    double mshift = C * (1 - exp(-SysPar1 * r / 2e5));

    return mshift;
  }

  // see GAP 2013-105, sec. 5.5
  double getX0(int version, double xmax, double loge)
  {
    // derivation of the X0-Xmax correlation function in x0_xmax_conversion.py
    // it is done separately for fixed energies.
    // The parameters are interpolated linearly in loge

    // The available configurations are:
    // 0: qgsjet, p + Fe
    // 1: epos, p + Fe
    // 2: qgsjet + epos, p + Fe
    // 3: qgsjet photons
    // 4: qgsjet, p
    // 5: epos, p

    const int nModels = 6;

    if (version < 0 || version > nModels - 1) {
      cout << "################\nbad version for x0-xmax correlation\n################\n";
      return 0;
    }

    const int nBins = 4;
    double loges[nBins] = { 1.860000e+01, 1.900000e+01, 1.950000e+01, 2.000000e+01 };

    const double pars0[nModels][nBins] = {
      { 6.456017e+02, 6.469736e+02, 6.489524e+02, 6.476612e+02 },
      { 6.453341e+02, 6.485834e+02, 6.501094e+02, 6.492983e+02 },
      { 6.445297e+02, 6.474811e+02, 6.491586e+02, 6.502373e+02 },
      { 6.432335e+02, 6.522376e+02, 6.632606e+02, 6.955606e+02 },
      { 6.474369e+02, 6.502122e+02, 6.488559e+02, 6.493190e+02 },
      { 6.542931e+02, 6.521569e+02, 6.479219e+02, 6.500125e+02 }
    };

    const double pars1[nModels][nBins] = {
      { 1.616533e-03, 1.444867e-03, 1.252005e-03, 1.582818e-03 },
      { 8.559603e-04, 1.039581e-03, 7.401688e-04, 8.037290e-04 },
      { 1.029642e-03, 1.117932e-03, 1.046979e-03, 9.735569e-04 },
      { 1.036600e-03, 9.843582e-04, 2.213315e-04, 2.842102e-04 },
      { 1.960041e-03, 1.456550e-03, 1.139915e-03, 1.162492e-03 },
      { 1.178549e-03, 8.926028e-04, 1.382968e-03, 1.341381e-03 }
    };

    const double pars2[nModels][nBins] = {
      { 4.319217e+00, 2.778060e+00, 1.106250e+00, 5.520477e+00 },
      { 3.101907e+00, 1.945493e+00, -3.072746e+00, -2.944935e+00 },
      { 3.807706e+00, 2.476609e+00, 9.049038e-01, 1.516479e+00 },
      { -1.954682e+01, 1.482052e+01, -4.562288e+01, 4.857960e+01 },
      { 5.610586e+00, 9.524124e-01, -2.232481e+00, -1.119290e+00 },
      { -5.986844e+00, -1.302564e+01, 9.704845e-01, 3.654311e+00 }
    };

    const double pars3[nModels][nBins] = {
      { 1.021706e-01, 5.903550e-02, 2.015838e-02, -1.283983e-01 },
      { 1.533218e-01, 6.349976e-02, 7.258369e-02, 3.835097e-03 },
      { 1.374037e-01, 7.072211e-02, 2.301778e-02, -2.622843e-02 },
      { -7.461706e-02, -2.455179e-01, 1.211286e-01, -1.878918e-01 },
      { 5.551634e-02, 7.942683e-02, 6.273052e-02, -1.132578e-02 },
      { 2.103373e-01, 1.940762e-01, -8.057426e-02, -1.497254e-01 }
    };

    xmax /= utl::gram / utl::cm2;
    double x0;

    double p0 = UnivTimeKG::int1(loges, pars0[version], loge, nBins);
    double p1 = UnivTimeKG::int1(loges, pars1[version], loge, nBins);
    double p2 = UnivTimeKG::int1(loges, pars2[version], loge, nBins);
    double p3 = UnivTimeKG::int1(loges, pars3[version], loge, nBins);

    double xmin = p0 - p3 / (2.*p1);
    double ymin = p2 + p3 * (xmin - p0) + p1 * pow(xmin - p0, 2);

    if (xmax < xmin)
      x0 = ymin;
    else
      x0 = p2 + p3 * (xmax - p0) + p1 * pow(xmax - p0, 2);

    if (std::isnan(x0) || x0 < 0) x0 = 0;
    return x0 * (utl::gram / utl::cm2);
  }
}