ShowerFRecData.cc

Go to the documentation of this file.

#include <evt/ShowerFRecData.h>
#include <utl/ErrorLogger.h>
#include <utl/PhysicalConstants.h>

#include <utl/LameShadowPtr_imp.h>
LAMESHADOWPTR_INST(evt::ShowerFRecData);

#include <sstream>

using namespace utl;
using namespace evt;
using namespace std;


inline
double
QuadraticSum(const ShowerFRecData::UncertaintyMap& uncertainties)
{
  double variance = 0.;
  for (ShowerFRecData::UncertaintyMapConstIterator iE = uncertainties.begin();
       iE != uncertainties.end(); ++iE)
    variance += pow(iE->second, 2);

  return variance;
}

ShowerFRecData::ShowerFRecData() :
  fEmEnergy(0),
  fTotalEnergy(0),
  fEnergyCutoff(0),
  fCorrelationXY(0),
  fCoreTimeError(0),
  fCorrelationThetaPhi(0),
  fSDTimeResidual(0)
{
}

void
ShowerFRecData::MakeEnergyDeposit()
{
  if (fEnergyDeposit) {
    ERROR("ShowerFRecData::EnergyDeposit already exists - Not replacing");
  }
  else
    fEnergyDeposit = new TabulatedFunctionErrors;
}


void
ShowerFRecData::MakeLongitudinalProfile()
{
  if (fLongitudinalProfile) {
    ERROR("ShowerFRecData::LongitudinalProfile already exists - Not replacing");
  }
  else
    fLongitudinalProfile = new TabulatedFunctionErrors;
}


void
ShowerFRecData::MakeFluorescencePhotons()
{
  if (fFluorescencePhotons) {
    ERROR("ShowerFRecData::FluorescenecePhotons already exists - Not replacing");
  }
  else
    fFluorescencePhotons = new TabulatedFunctionErrors;
}


void
ShowerFRecData::MakeCherenkovPhotons()
{
  if (fCherenkovPhotons) {
    ERROR("ShowerFRecData::CherenkovPhotons already exists - Not replacing");
  }
  else
    fCherenkovPhotons = new TabulatedFunction;
}


void
ShowerFRecData::MakeGHParameters(const VGaisserHillasParameter& gh)
{
  if (fGHParameters) {
    ERROR("ShowerFRecData::GHParameters already exists - Not replacing");
  }
  else
    fGHParameters = gh.Clone();
}


double
ShowerFRecData::GetEmEnergyError(const EUncertaintyType type)
  const
{
  if (fEmEnergyError.find(type) !=  fEmEnergyError.end())
    return fEmEnergyError.find(type)->second;
  else {
    if (type == eTotal)
      return sqrt(QuadraticSum(fEmEnergyError));
    else {
      ostringstream err;
      err << " No uncertainty type " << type
          << " available ";
      ERROR(err.str());
      throw utl::MissingEventDataException(err.str());
    }
  }
}

bool
ShowerFRecData::HasEmEnergyError(const EUncertaintyType type)
  const
{
  if (type == eTotal)
    return !fEmEnergyError.empty();
  else
    return fEmEnergyError.find(type) !=  fEmEnergyError.end();
}

double
ShowerFRecData::GetTotalEnergyError(const EUncertaintyType type)
  const
{
  if (fTotalEnergyError.find(type) !=  fTotalEnergyError.end())
    return fTotalEnergyError.find(type)->second;
  else {
    // quadratically sum up errors for eTotal
    if (type == eTotal)
      return sqrt(QuadraticSum(fTotalEnergyError));
    else {
      ostringstream err;
      err << " No uncertainty type " << type
          << " available ";
      ERROR(err.str());
      throw utl::MissingEventDataException(err.str());
    }
  }
}


bool
ShowerFRecData::HasTotalEnergyError(const EUncertaintyType type)
  const
{
  if (type == eTotal)
    return !fTotalEnergyError.empty();
  else
    return fTotalEnergyError.find(type) !=  fTotalEnergyError.end();
}


double
ShowerFRecData::GetXmaxError(const EUncertaintyType type)
  const
{
  if (fXmaxError.find(type) !=  fXmaxError.end())
    return fXmaxError.find(type)->second;
  else {
    if (type == eTotal)
      return sqrt(QuadraticSum(fXmaxError));
    else {
      ostringstream err;
      err << " No uncertainty type " << type
          << " available ";
      ERROR(err.str());
      throw utl::MissingEventDataException(err.str());
    }
  }
}


bool
ShowerFRecData::HasXmaxError(const EUncertaintyType type)
  const
{
  if (type == eTotal)
    return !fXmaxError.empty();
  else
    return fXmaxError.find(type) !=  fXmaxError.end();
}


void
ShowerFRecData::SetCoreTime(const utl::TimeStamp& eyeTriggerTime,
                            const double rp, const double chi0,
                            const double t0)
{

  fCoreTime = eyeTriggerTime + TimeInterval(t0-rp/tan(chi0)/kSpeedOfLight);

}

utl::Point
ShowerFRecData::CalculatePointOnShower(const utl::TimeStamp& timeAtTelescope,
                                       const utl::Point& telescopePosition)
  const
{
  /*

    for distance from shower core to emission point (a)
    law of cosines: b^2 = a^2 + c^2 - 2 * a * c * cos(beta)
    c = distance(telescope, core)
    b - a = (telTime - coreTime)/c_light = delta

  */


  const Vector coreTelescopeVec  = telescopePosition - fCore;
  const double c = coreTelescopeVec.GetMag();
  const double cosBeta = CosAngle(fAxis, coreTelescopeVec);
  const double delta = (timeAtTelescope - fCoreTime).GetInterval() * kSpeedOfLight;
  const double distance =
    -(c * c - delta * delta)/(2 * (delta + c * cosBeta));
  const utl::Point pointOnShower = fCore - distance * fAxis;
  return pointOnShower;

}