RdLDFChargeExcessCorrector.cc

Go to the documentation of this file.

#include "RdLDFChargeExcessCorrector.h"

#include <fwk/CentralConfig.h>
#include <fwk/MagneticFieldModel.h>

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

#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/BasicVector.h>
#include <utl/GeometryUtilities.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/AxialVector.h>

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

#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <det/VManager.h>

#include <TF1.h>

#include <math.h>

using namespace revt;
using namespace fwk;
using namespace det;

namespace RdLDFChargeExcessCorrector {

  RdLDFChargeExcessCorrector::RdLDFChargeExcessCorrector() :
      fInfoLevel(eNone),
      fMeanAsymmetry(0.),
      fMeanAsymmetryError(0.),
      fAsymmetryParametrization(""),
      fAsymmetryDefinition("Mean")
  {
  }

  RdLDFChargeExcessCorrector::~RdLDFChargeExcessCorrector()
  {
  }

  VModule::ResultFlag RdLDFChargeExcessCorrector::Init()
  {
    INFO("RdLDFChargeExcessCorrector::Init()");

    utl::Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdLDFChargeExcessCorrector");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("MeanAsymmetry").GetData(fMeanAsymmetry);
    topBranch.GetChild("MeanAsymmetryError").GetData(fMeanAsymmetryError);
    topBranch.GetChild("AsymmetryParametrization").GetData(fAsymmetryParametrization);
    topBranch.GetChild("AsymmetryDefinition").GetData(fAsymmetryDefinition);

    return eSuccess;

  }

  VModule::ResultFlag RdLDFChargeExcessCorrector::Run(evt::Event& event)
  {
    if(!event.HasREvent()) {
      WARNING("No radio event found!");
      return eSuccess;
    }

    REvent& rEvent = event.GetREvent();
    evt::ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

    if(!showerrrec.HasCorePosition() || !showerrrec.HasAxis()) {
        WARNING("No core or axis. Charge-excess correction is not performed.");
        return eSuccess;
    }


    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    const utl::Vector showerAxis = showerrrec.GetAxis();

    utl::CoordinateSystemPtr baryCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

    const utl::Point core = showerrrec.GetCorePosition();

    utl::Vector localMagneticField = fwk::MagneticFieldModel::instance().GetMagneticFieldVector(event.GetRecShower().GetRRecShower().GetCoordinateOrigin(), event.GetREvent().GetHeader().GetTime().GetGPSSecond());

    utl::RadioGeometryUtilities transform(showerAxis,baryCS,fwk::MagneticFieldModel::instance().GetMagneticFieldVector(event.GetRecShower().GetRRecShower().GetCoordinateOrigin(), event.GetREvent().GetHeader().GetTime().GetGPSSecond()));

    const double geomag_angle = utl::Angle(localMagneticField,showerAxis);

    for(REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin(); sIt != rEvent.CandidateStationsEnd(); ++sIt) {

      Station& station = *sIt;
      const StationRecData& sRec = station.GetRecData();

      if(!sRec.GetPulseFound())
        continue;

      const utl::Point sPos = rDetector.GetStation(station).GetPosition();

      const utl::Vector sVector = sPos - core;

      double magx(0), magy(0), magz(0);
      transform.GetVectorInShowerPlaneVxB(magx, magy, magz, utl::Point(sVector.GetX(baryCS), sVector.GetY(baryCS), sVector.GetZ(baryCS), baryCS));

      const double t_angle = atan2(magy,magx);

      double eps = 0.0;
      double eps_err = 0.0;

      if(fAsymmetryDefinition == "Mean") {
        eps = fMeanAsymmetry;
        eps_err = fMeanAsymmetryError;
      }

      if(fAsymmetryDefinition == "Polarization") {

        //  Warning / TODO: #warning: not yet fully implemented
        //eps = station.GetRecData().GetParameter(eChargeExcessAsymmetry);
        //eps_err = station.GetRecData().GetParameterError(eChargeExcessAsymmetry);
      }

      if(fAsymmetryDefinition == "Geometry") {
      }

      if(fAsymmetryDefinition == "Parametrization") {
        const double dist = utl::RadioGeometryUtilities::GetDistanceToAxis(showerAxis, core, sPos);
        TF1 *asymm_parametrization = new TF1("asymmetry_parametrization",fAsymmetryParametrization.c_str(),0,500);
        // Warning / TODO: #warning: harcoded to 500 m!
        eps = asymm_parametrization->Eval(dist);
        // Warning / TODO: #warning: for parametrization the uncertainty is set to zero!
        eps_err = 0.0;
        asymm_parametrization->Delete();
      }

      const double Amp = station.GetRecData().GetParameter(eSignal);
      const double Amp_err = station.GetRecData().GetParameterError(eSignal);

      double A0 = Amp / sqrt(pow(eps,2) + 2*eps*cos(t_angle)*sin(geomag_angle) + pow(sin(geomag_angle),2));
      double A0_err = sqrt(  pow(Amp_err / sqrt(pow(eps,2) + 2*eps*cos(t_angle)*sin(geomag_angle) + pow(sin(geomag_angle),2)),2)
                           + pow( (eps_err * Amp * (eps + eps*cos(t_angle)*sin(geomag_angle)) ) / ( pow(pow(eps,2) + 2*eps*cos(t_angle)*sin(geomag_angle)
                           + pow(sin(geomag_angle),2),1.5) ) ,2) );
      station.GetRecData().SetParameter(eSignal, A0);
      station.GetRecData().SetParameterError(eSignal, A0_err);
    }

    return eSuccess;
  }

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

}