STimeVariance.cc

Go to the documentation of this file.

#include <utl/Math.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/ErrorLogger.h>
#include <sdet/STimeVariance.h>

using namespace utl;
using namespace sdet;


namespace sdet {

  inline
  double
  RelativeTrackLength(const double cosTheta)
  {
    const double kRadius = 1.8;
    const double kHeight = 1.2;
    const double sinTheta = sqrt(1 - Sqr(cosTheta));
    return 1 / (cosTheta + 2*kHeight*sinTheta/(kPi*kRadius));
  }


  inline
  double
  MeanT50(const double cosTheta, const double distance)
  {
    // this form can produce fake minima in geometry fit
    //return 0.8*nanosecond * (0.53*cosTheta - 0.11) * distance/meter;

    // alternative proposed by esteban but with better asymptotic behaviour
    constexpr double delta = 20;
    constexpr double p1 = delta * (-0.11);
    constexpr double p2 = delta * 0.53;
    const double p3 = p2 * cosTheta;
    const double p4 = std::exp(p1 + p3);
    constexpr double p5 = std::exp(0.1*delta);
    const double fz = std::log(p4 + p5) / delta;
    return (0.8*nanosecond/meter) * fz * distance;
  }


  inline
  double
  LimitT50(const double t50, const double cosTheta, const double distance)
  {
    // correction a la Xavier -- March 2009
    const double limit = 0.5 * MeanT50(cosTheta, distance);
    return (t50 < limit) ? limit : t50;
  }


  inline
  double
  VarianceForm2006(const double n, const double t50, const double a2, const double b2)
  {
    // C. Bonifazi, A. Letessier-Selvon: GAP2006-16
    return a2 * Sqr(2*t50/n) * (n - 1) / (n + 1) + b2;
  }


  inline
  double
  VarianceForm2007(const double n, const double t50, const double tmin, const double a2, const double b2)
  {
    // M. Horvat, D. Veberic: GAP2007-057
    return a2 * Sqr((t50 + tmin)/(n + 1)) * n / (n + 2) + b2;
  }


  inline
  double
  VarianceForm2012(const double n, const double t50, const double a2, const double b2)
  {
    // D. Harari, D. Melo, S. Mollerach, E. Roulet, and G. van Aar: GAP2012-145
    return a2 * Sqr(2*t50/(n - 1)) * n / (n + 2) + b2;
  }


  double
  STimeVariance::GetTimeSigma2(const double signal, const double t50, const double cosTheta, const double distance)
    const
  {
    const double n = signal / RelativeTrackLength(cosTheta);

    switch (fModel) {
    case eICRC2005:
      return (600*nanosecond2 + 1.2*Sqr(t50/signal)) * (0.4 + 1.2*cosTheta);
    case eGAP2006_016:
      return VarianceForm2006(n, t50, 1, 147*nanosecond2);
    case eNIMA:
      return VarianceForm2006(n, t50, 0.99, 147*nanosecond2);
    case eCDASv4r4:
      return VarianceForm2007(n, t50, 15*nanosecond, 5, 134*nanosecond2);
    case eGAP2007_057:
      return VarianceForm2007((n < 2) ? 2 : n, t50, 10*nanosecond, 2.4, 134*nanosecond2);
    case eCDASv4r6:
      return VarianceForm2006((n < 2) ? 2 : n, t50, 0.36, 212*nanosecond2);
    case eCDASv4r8:
      {
        const double dt = LimitT50(t50, cosTheta, distance);
        // March 2009
        const double a2 = 0.67612 + cosTheta*(0.16106 + cosTheta*(-0.47641));
        const double b2 = (128 + cosTheta*(-184 + cosTheta*413)) * nanosecond2;
        return VarianceForm2006((n < 2) ? 2 : n, dt, a2, b2);
      }
    case eCDASv5r0:
      {
        // GAP2018-048
        const double dt = LimitT50(t50, cosTheta, distance);
        const double a2 = 0.64871 + cosTheta*(0.22365 + cosTheta*(-0.49971));
        const double b2 = (141.24  + cosTheta*(-208.90 + cosTheta*412.19)) * nanosecond2;
        return VarianceForm2006((n < 2) ? 2 : n, dt, a2, b2);
     }
    case eGAP2012_145:
      {
        const double meanT50 = MeanT50(cosTheta, distance);
        const double dt = (n < 4 && t50 < meanT50) ? meanT50 : t50;
        // a and b are the same as in eCDASv5r0:
        const double a2 = 0.64871 + cosTheta*(0.22365 + cosTheta*(-0.49971));
        const double b2 = (141.24  + cosTheta*(-208.90 + cosTheta*412.19)) * nanosecond2;
        return VarianceForm2012((n < 2) ? 2 : n, dt, a2, b2);
      }
    case eNone:
    default:
      ERROR("You should specify a TimeVariance model!");
      return 0;
    }
  }

}