RdStationPolarizationRejector.h

Go to the documentation of this file.

#ifndef _RdStationPolarizationRejector_RdStationPolarizationRejector_h_
#define _RdStationPolarizationRejector_RdStationPolarizationRejector_h_

#include <fwk/VModule.h>

#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/CoordinateSystem.h>
#include <utl/RadioGeometryUtilities.h>

namespace evt {

  class ShowerRRecData;

}

namespace RdStationPolarizationRejector {

  class RdStationPolarizationRejector : public fwk::VModule {

  public:

    fwk::VModule::ResultFlag Init() override;
    fwk::VModule::ResultFlag Run(evt::Event& event) override;
    fwk::VModule::ResultFlag Finish() override;

    void SetMaxChargeExcess(double aMax) { fMaxChargeExcess = aMax; }
    void SetMinChargeExcess(double aMin) { fMinChargeExcess = aMin; }

    double GetProbability() { return fProbability; }
    void SetProbability(double prob);

    void GetCovarianceMatrix(const evt::Event& event, double covarianceMatrix[3], const evt::ShowerRRecData& rShower,
                             const utl::CoordinateSystemPtr& cs);
    TVector3 GetStationPositionvxBvxvxB(const utl::Point stationPosition, const utl::Point& corePosition,
                                        const utl::RadioGeometryUtilities&);
    void GetLocalPolMaxima(double& pol1, double& pol2, const double angleToMagneticField, double maxChargeExcess);

  private:

    bool fConstantChargeExcess = false;
    double fMaxChargeExcess = 0;    //maximal acceptable charge excess strength
    double fMinChargeExcess = 0;    //minimal acceptable charge excess strength
    double fProbability = 0;        //percentile of region around SD core position
    double fChiSquared = 0;         //set depending on fProbability
    double fMaxBeta = 0;            //maximum discrepancy between expected and measured polarization with which a station is accepted immediately
    double fMaxSigmaBeta = 0;
    bool fZenithCut = false;        //If set to true, polarization rejection will only be performed for reasonable zenith angles
    int fNoiseErrorModel = 0;

    //storing charge excess limits
    double fMaxChargeExcess30[20] = { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.45, 0.38, 0.34, 0.32, 0.29,
                                      0.21, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3 };
    double fMinChargeExcess30[20] = { 0.0, 0.1, 0.1, 0.08, 0.05, 0.02, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0,
                                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
    double fMaxChargeExcess40[20] = { 0.37, 0.37, 0.32, 0.32, 0.33, 0.33, 0.33, 0.35, 0.35, 0.31,
                                      0.31, 0.28, 0.25, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3 };
    double fMinChargeExcess40[20] = { 0.0, 0.09, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.05, 0.05, 0.05, 0.0,
                                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
    double fMaxChargeExcess50[20] = { 0.25, 0.25, 0.25, 0.25, 0.26, 0.3, 0.3, 0.35, 0.35, 0.33, 0.33,
                                      0.32, 0.32, 0.24, 0.24, 0.23, 0.22, 0.22, 0.22, 0.22 };
    double fMinChargeExcess50[20] = { 0.0, 0.0, 0.0, 0.0, 0.05, 0.05, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
                                      0.1, 0.05, 0.05, 0.0, 0.0, 0.0, 0.0, 0.0 };
    double fMaxChargeExcess60[20] = { 0.03, 0.03, 0.05, 0.08, 0.1, 0.12, 0.13, 0.15, 0.17, 0.18, 0.19,
                                      0.2, 0.2, 0.2, 0.21, 0.21, 0.22, 0.22, 0.21, 0.21 };
    double fMinChargeExcess60[20] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
    double fMaxZenith[12] = { 60. * utl::degree, 70. * utl::degree, 70. * utl::degree, 70. * utl::degree,
                              70. * utl::degree, 60. * utl::degree, 50. * utl::degree, 50. * utl::degree,
                              50. * utl::degree, 50. * utl::degree, 50. * utl::degree, 50. * utl::degree };

    void GetCoreErrorEllipsevxBvxvxB(TVector3& semiMajorAxisvxBvxvxB, TVector3& semiMinorAxisvxBvxvxB,
                                     double& rotationAngleInVxB, const utl::Point& semiMajorAxis,
                                     const utl::Point& semiMinorAxis, const utl::Point& corePosition,
                                     const utl::RadioGeometryUtilities& rdGeometryUtilities);

    void GetCoreErrorEllipse(utl::Point& semiMajorAxis, utl::Point& semiMinorAxis, const utl::Point& corePosition,
                             const double covarianceMatrix[3], const utl::CoordinateSystemPtr& cs);

    bool IsStationInErrorEllipse(const TVector3& stationPositionvxBvxvxB, const TVector3& semiMajorAxisVxB,
                                 const TVector3& semiMinorAxisVxB, const double rotationAngle);

    void GetErrorEllipseTangentPoints(TVector3& ellipseTangentPoint1, TVector3& ellipseTangentPoint2,
                                      const TVector3& stationPositionvxBvxvxB, const TVector3& semiMajorAxisVxB,
                                      const TVector3& semiMinorAxisVxB, const double rotationAngle);

    double GetTangentPol(const TVector3& stationPositionvxBvxvxB, const TVector3& tangentPoint,
                         const double chargeExcess, const double angleToMagneticField);

    TVector3 GetEfieldInShowerPlane(const utl::Vector& efield, const utl::CoordinateSystemPtr& cs,
                                    const utl::RadioGeometryUtilities& rdGeometryUtilities);

    double GetPolAngleUncertainty(const double signalToNoise, const TVector3 efield,
                                  const double noiseVxB, const double noiseVxVxB);

    double GetMaxChargeExcess(const double zenith, const double distToAxis, const double semiMajorAxisLength);
    double GetMinChargeExcess(const double zenith, const double distToAxis, const double semiMajorAxisLength);

    bool CheckZenith(const double zenith, const double azimuth);

  REGISTER_MODULE("RdStationPolarizationRejector", RdStationPolarizationRejector);
  };

}

#endif