RdAirShowerReconstruction.cc

Go to the documentation of this file.

#include "RdAirShowerReconstruction.h"

#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerSRecData.h>

#include <revt/REvent.h>

#include <utl/config.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>

using namespace fwk;
using namespace evt;
using namespace utl;
using namespace revt;

namespace RdAirShowerReconstruction {

  RdAirShowerReconstruction::RdAirShowerReconstruction() :
      fInfoLevel(eNone),
      fEnergyFitId(0),
      fEnergyFitParameter(0),
      fEnergyKappa(0),
      fEnergyB(0),
      fXmaxFitId(0),
      fXmaxFitParameter(0),
      fXmaxGroundDepth(0),
      fXmaxBbar(0),
      fXmaxParamA(0),
      fXmaxParamB(0),
      fZenith(0)
  { }


  VModule::ResultFlag
  RdAirShowerReconstruction::Init()
  {
    INFO("RdAirShowerReconstruction::Init()");
    utl::Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdAirShowerReconstruction");

    // Read in the configurations of the xml file
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    // Energy reconstruction parameters
    topBranch.GetChild("energyFitId").GetData(fEnergyFitId);
    topBranch.GetChild("energyFitParameter").GetData(fEnergyFitParameter);
    topBranch.GetChild("energyKappa").GetData(fEnergyKappa);
    topBranch.GetChild("energyB").GetData(fEnergyB);

    // Xmax reconstruction parameters
    topBranch.GetChild("xmaxFitId").GetData(fXmaxFitId);
    topBranch.GetChild("xmaxFitParameter").GetData(fXmaxFitParameter);
    topBranch.GetChild("xmaxGroundDepth").GetData(fXmaxGroundDepth);
    topBranch.GetChild("xmaxBbar").GetData(fXmaxBbar);
    topBranch.GetChild("xmaxParamA").GetData(fXmaxParamA);
    topBranch.GetChild("xmaxParamB").GetData(fXmaxParamB);

    return eSuccess;
  }


  VModule::ResultFlag
  RdAirShowerReconstruction::Run(evt::Event& event)
  {
    // Checking the existance of the reconstructed shower
    if (!event.HasRecShower()) {
      WARNING("No RecShower found!");
      return eSuccess;
    }
    if (!event.GetRecShower().HasRRecShower()) {
      WARNING("No RRecShowerFound!");
      return eSuccess;
    }

    // Getting variables
    ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();
    fZenith = showerrrec.GetZenith();

    // Getting results of the RdLDFMultiFitter module
    const std::vector<TF1>& fitResult = showerrrec.GetRdLDF();

    // Checking the existance of functions with required indices
    // ... for energy
    if (fitResult.size() < fEnergyFitId + 1) {
      WARNING("Fit with EnergyFitId does not exist!");
      return eSuccess;
    }
    // ... for Xmax
    if (fitResult.size() < fXmaxFitId + 1) {
      WARNING("Fit with XmaxFitId does not exist!");
      return eSuccess;
    }

    // Checking the existance of functions with required indices
    // ... for energy
    if (fitResult[fEnergyFitId].GetNpar() < int(fEnergyFitParameter + 1)) {
      WARNING("EnergyFitParameter does not exist!");
      return eSuccess;
    }
    if (fitResult[fEnergyFitId].GetNpar() < int(fXmaxFitParameter + 1)) {
      WARNING("XmaxFitParamter does not exist!");
      return eSuccess;
    }

    // Assigning variables
    const double E0 = fitResult[fEnergyFitId].GetParameter(fEnergyFitParameter);
    const double eta = fitResult[fXmaxFitId].GetParameter(fXmaxFitParameter);
    const double E0_err = fitResult[fEnergyFitId].GetParError(fEnergyFitParameter);
    const double eta_err = fitResult[fXmaxFitId].GetParError(fXmaxFitParameter);

    // Calculating energy and Xmax
    double energy_err = 0;
    const double energy = CalculateEnergy(E0, E0_err, energy_err);
    double Xmax_err = 0;
    const double Xmax = CalculateXmax(eta, eta_err, Xmax_err);

    if (std::isnan(energy) || std::isnan(energy_err)) {
      WARNING("Can not reconstruct energy!");
      return eSuccess;
    }
    if (std::isnan(Xmax) || std::isnan(Xmax_err)) {
      WARNING("Can not reconstruct Xmax!");
      return eSuccess;
    }

    showerrrec.SetParameter(eReconstructedEnergyTunkaRex, energy);
    showerrrec.SetParameter(eReconstructedXmaxTunkaRex, Xmax);
    showerrrec.SetParameterError(eReconstructedEnergyTunkaRex, energy_err);
    showerrrec.SetParameterError(eReconstructedXmaxTunkaRex, Xmax_err);

    return eSuccess;
  }


  double
  RdAirShowerReconstruction::CalculateEnergy(const double E0, const double E0_err, double& energy_err)
  {
    energy_err = fabs(fEnergyKappa * fEnergyB * pow(E0, fEnergyB - 1) * E0_err);
    return fEnergyKappa * pow(E0, fEnergyB);
  }


  double
  RdAirShowerReconstruction::CalculateXmax(const double eta, const double eta_err, double& Xmax_err)
  {
    Xmax_err = fabs(fXmaxParamB / (fXmaxBbar - eta) * eta_err);
    const double distance = fXmaxParamA + fXmaxParamB*log(fXmaxBbar - eta);
    return fXmaxGroundDepth / cos(fZenith) - distance;
  }


  VModule::ResultFlag
  RdAirShowerReconstruction::Finish()
  {
    INFO("RdAirShowerReconstruction::Finish()");
    return eSuccess;
  }

}