EnergyCalculation.cc

Go to the documentation of this file.

#include "EnergyCalculation.h"
#include <det/Detector.h>
#include <utl/Math.h>
#include <utl/PhysicalConstants.h>
#include <evt/Event.h>
#include <sevt/Meteo.h>
#include <sevt/SEvent.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CentralConfig.h>
#include <cmath>

using namespace evt;
using namespace sevt;
using namespace utl;
using namespace fwk;
using namespace std;


namespace EnergyCalculationPG {

  void
  EnergyCalculation::GetCorrectedS1000Geomagnetic(const double theta,
                                                  const double phi,
                                                  const double s1000,
                                                  const double s1000Err,
                                                  double& sGeo1000,
                                                  double& sGeo1000Err)
  {
    const double sinTheta = sin(theta);
    const double cosTheta = cos(theta);
    // From Moritz Munchmeyer's Thesis 2012, Eq. (4.3.6):
    const double thetaB = 54.8 * degree;
    const double phiB = 87.4 * degree;
    const double axisBMag =
      sinTheta * cos(phi) * sin(thetaB) * cos(phiB) +
      sinTheta * sin(phi) * sin(thetaB) * sin(phiB) +
      cosTheta * cos(thetaB);
    const double correctionFactor =
      1 + 4.2e-3 * pow(cosTheta, -2.8) * (1 - Sqr(axisBMag));

    sGeo1000 = s1000 / correctionFactor;
    sGeo1000Err = s1000Err / correctionFactor;
  }


  VModule::ResultFlag
  EnergyCalculation::Run(Event& event)
  {
    if (!event.HasSEvent() || !event.HasRecShower())
      return eContinueLoop;
    auto& recShower = event.GetRecShower();
    if (!recShower.HasSRecShower())
      return eContinueLoop;

    // get the needed variables from the event calculate the energy
    // and the energy error set the energy to the RecShower

    auto& sRecShower = recShower.GetSRecShower();

    const double showerSize = sRecShower.GetShowerSize();
    const double showerSizeErr = sRecShower.GetShowerSizeError();

    const Point core = sRecShower.GetCorePosition();
    const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(core);
    const auto& axis = sRecShower.GetAxis();
    const double theta = axis.GetTheta(coreCS);
    const double phi = axis.GetPhi(coreCS);

    double s1000Geo = 0;
    double s1000ErrGeo = 0;
    GetCorrectedS1000Geomagnetic(theta, phi, showerSize, showerSizeErr, s1000Geo, s1000ErrGeo);
    // only set corrected variable, don't overwrite s1000!
    sRecShower.SetSizeGeomagneticCorr(s1000Geo, s1000ErrGeo);

    return eSuccess;
  }

}