ShowerFrontFunction.cc

Go to the documentation of this file.

#include "ShowerFrontFunction.h"
#include "SdHorizontalReconstruction.h"
#include "Utilities.h"

#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Math.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/MuonArrivalTime.h>
#include <utl/AugerException.h>

#include <tls/TankResponseUtilities.h>

#include <sdet/STimeVariance.h>

#include <limits>

using namespace SdHorizontalReconstructionNS;
using namespace utl;
using namespace std;
using namespace tls;


ShowerFrontFunction::ShowerFrontFunction(const SdHorizontalReconstruction& config,
                                 const Point& core,
                                 const StationList& list,
                                 const ExternalGeometryData& gd)
  : fConfig(config), fCore(core), fList(list), fExt(gd)
{
}


double
ShowerFrontFunction::operator()(const std::vector<double>& pars)
  const
{
  for (size_t i = 0; i < 4; ++i)
    if (std::isnan(pars[i]))
      return numeric_limits<double>::max();

  const double theta = pars[0];
  const double phi = pars[1];
  const double r = pars[2]*km;
  const double ct = pars[3]*km;

  double result = 0;

  const Point origin(r, theta, phi, fConfig.fBaryCS, Point::kSpherical);

  for (StationList::const_iterator
         sIt = fList.begin(),
         sEnd = fList.end();
       sIt != sEnd; ++sIt)
  {
    if (sIt->fRejected) continue;

    double meanCT, sigmaCT;
    Predict(meanCT, sigmaCT, *sIt, origin, ct);
    result += 0.5 * Sqr((sIt->fCTime - meanCT)/sigmaCT);
  }

  {
    const double delta[2] = {
      theta - fExt.fPar[eThetaExt],
      phi - fExt.fPar[ePhiExt]
    };

    for (size_t i = 0; i < 2; ++i)
      for (size_t j = 0; j < 2; ++j)
        result += 0.5 * delta[i] * fExt.fInvCov[eThetaExt+i][eThetaExt+j] * delta[j];
  }

  return result;
}


void
ShowerFrontFunction::Predict(double& meanCT, double& sigmaCT,
                             const StationData& sd,
                             const Point& originPos,
                             const double originCT)
  const
{
  const Vector posToOrigin = originPos - sd.fPos;
  const double sTheta = Angle(posToOrigin, (sd.fPos - fConfig.fEarthCenter));

  switch (fConfig.fShowerFrontModel)
  {
    case eSphere:
    {
      const sdet::STimeVariance& timeVariance = sdet::STimeVariance::GetInstance();

      meanCT = originCT + posToOrigin.GetMag();
      sigmaCT = sqrt(timeVariance.GetTimeSigma2(sd.fSignal, sd.fT50, cos(sTheta))) * kSpeedOfLight;
    }
    break;
    case eMuonArrival:
    {
      // uses model of muon arrival times from Lorenzo Cazon
      // model predicts time relative to plane shower front

      const Vector axis = originPos-fCore;
      fMat.SetThetaAndDistance(Angle(axis, fCore-fConfig.fEarthCenter), axis.GetMag());

      static const double kVerticalTrackLength = TankMeanTrackLength(0);
      static const double kStationTimeJitter = 20*ns; // according to USC code

      const int nmuon = round(sd.fSignal / (TankMeanTrackLength(sTheta)/kVerticalTrackLength));
      double rho, delta;
      GetRhoAndDelta(rho, delta, sd.fPos, fCore, originPos);
      fMat.SetCoordinates(rho, delta);
      double tmean, tsigma;
      fMat.GetApproximateMeanAndStDev(tmean, tsigma, nmuon < 1 ? 1: nmuon);

      meanCT = originCT + (originPos-fCore).GetMag() - delta + tmean*kSpeedOfLight;
      sigmaCT = sqrt(Sqr(tsigma)+Sqr(kStationTimeJitter))*kSpeedOfLight;
    }
    break;
    default:
    throw utl::NonExistentComponentException("shower front model not implemented");
  }
}