MTimeVariance.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 <mdet/MTimeVariance.h>

using namespace utl;
using namespace mdet;


namespace mdet {

  inline
  double
  VarianceForm1(unsigned int n, double r, double a, double b, double c)
  {

    // approximating the median variance by the mean variance 
    // valid for an exponential distribution in the large sample limit
    const double var = (a*a + r*r * b*b)/n + c*c;

    return var;
  }


  // Variance model from Dec 28 2017
  // Parameters:
  //    n: number of muons
  //    r: core distance
  //    sigma0: width of the single muon time distribution at the reference distance
  //    alpha: slope of the muon width with the distance
  //    sigma2: all contributions independent of the number of muons (RMS GPS SD and MD sinchronisation)  

  inline
  double
  VarianceForm2(unsigned int n, double r, double sigma0, double alpha, double sigma2)
  {

    // Width of the time distribution of a single muon
    // Modelled as a broken function:
    //  constant if distance is less than 200 m
    //  linear with r if distance is greater than 200 m 

    double sigma1;
    //const double rcut = 200;   
    const double r0 = 450;  // reference distance

    if (r < 200) {
      sigma1 = 25; 
    } else { 
      sigma1 = sigma0 + alpha * (r - r0);
    }

    // Conversion factor between the variance of the mean to the variance of the median
    if (n <= 2) {
      alpha = 1;
    } else {
      alpha = double(n) / double(n+1) * kPi / 2.;
    }

    // Variance propagation
    const double var = alpha*(sigma1*sigma1)/n + sigma2*sigma2;

    return var;
  }

}


double
MTimeVariance::GetTimeSigma2(const unsigned int nmuons, const double distance)
  const
{
  switch (fModel) {
  case eVersion1: {
    const double a = 21;  // in ns
    const double b = 0.1204; // in ns/m 
    const double sigma2 = 17; // in ns
    const double var = VarianceForm1(nmuons, distance, a, b, sigma2);
    return var;
  }
  case eVersion2: {
    const double sigma0 = 58.1;  // in ns
    const double alpha = 0.130; // in ns/m 
    const double sigma2 = 17; // in ns
    const double var = VarianceForm2(nmuons, distance, sigma0, alpha, sigma2);
    return var;
  }
  case eNone:
  default:
    ERROR("You should specify a TimeVariance model!");
    return 0;
  }
}