RdPolarizationReconstructor.h

Go to the documentation of this file.

#ifndef _RdPolarizationReconstructor_h_
#define _RdPolarizationReconstructor_h_

#include <tls/MuonProfile.h>
#include <tls/MuonProfileUtilities.h>
#include <fwk/VModule.h>
#include <revt/Station.h>
#include <revt/StationRecData.h>
#include <revt/Channel.h>
#include <utl/ParameterStorage.h>
#include <rdet/RDetector.h>
#include <utl/AxialVector.h>
#include <utl/GeometryUtilities.h>

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>

#include <string>
#include <map>


namespace evt {
  class Event;
}

namespace revt {
  class Channel;
}

namespace RdPolarizationReconstructor {

  const Complex I(0.0, 1.0);


  // Functionality for slicing the vectors is needed therefore we choose a boost
  // vector and a boost vector_range as a proxy in stead of an offline Trace
  typedef boost::numeric::ublas::vector<Complex> BoostVecC;
  typedef boost::numeric::ublas::vector_range<BoostVecC> BoostVecRangeC;
  typedef boost::numeric::ublas::vector<double> BoostVecD;
  typedef utl::TraceV3C::SizeType SizeT; // Shorthand


  // These maps are for storing intermediate values of the analytic error estimation method
  // The first map is used to store the covariances when the Analytic error calculation
  // method is used. It stores key value pairs with the value covariance (Complex) between
  // channel1 (char), channel2 (char) as the key. It is used in the
  // StokesParameters::GetSampleCovariance() method.
  // The second map is also only used for the analytic method. It stores the estimated
  // uncertainties of the single channels as a variance.
  typedef std::map< std::pair< char, char >, BoostVecC > SampleCovariancesT;
  typedef std::map< char, double > OneChannelPulseVariancesT;


  // The objects of this class contain the raw polarization data that are obtained from the background
  // and the raw polarization data that are obtained from the signal
  class RawStokesParameters {

  public:
    RawStokesParameters() { }

    template<class T>
    RawStokesParameters(const T& slice_x, const T& slice_y) :
      fSliceX(slice_x),
      fSliceY(slice_y),
      fSliceA((slice_x + slice_y)/sqrt(2)),
      fSliceB((slice_x - slice_y)/sqrt(2)),
      fSliceL((slice_x + I*slice_y)/sqrt(2)),
      fSliceR((slice_x - I*slice_y)/sqrt(2)),
      fAverageAbs2X(pow(norm_2(fSliceX), 2)/fSliceX.size()),
      fAverageAbs2Y(pow(norm_2(fSliceY), 2)/fSliceY.size()),
      fAverageAbs2A(pow(norm_2(fSliceA), 2)/fSliceA.size()),
      fAverageAbs2B(pow(norm_2(fSliceB), 2)/fSliceB.size()),
      fAverageAbs2L(pow(norm_2(fSliceL), 2)/fSliceL.size()),
      fAverageAbs2R(pow(norm_2(fSliceR), 2)/fSliceR.size()) //,
      //fRawStokesI(fAverageAbs2X + fAverageAbs2Y), // unused. LN.
      //fRawStokesQ(fAverageAbs2X - fAverageAbs2Y), // unused. LN.
      //fRawStokesU(fAverageAbs2A - fAverageAbs2B), // unused. LN.
      //fRawStokesV(fAverageAbs2L - fAverageAbs2R) // unused. LN.
    { }

    double GetRawStokesI() const { return fAverageAbs2X + fAverageAbs2Y; }
    double GetRawStokesQ() const { return fAverageAbs2X - fAverageAbs2Y; }
    double GetRawStokesU() const { return fAverageAbs2A - fAverageAbs2B; }
    double GetRawStokesV() const { return fAverageAbs2L - fAverageAbs2R; }
    SizeT GetSize() const { return fSliceX.size(); }
    BoostVecRangeC GetSliceX(SizeT pos, SizeT size); // this method should in principle be declared as const but boost doesn't like that
    BoostVecRangeC GetSliceY(SizeT pos, SizeT size); // idem
    Complex GetIndex(char channel, int index) const;
  private:
    BoostVecC fSliceX;
    BoostVecC fSliceY;
    BoostVecC fSliceA;
    BoostVecC fSliceB;
    BoostVecC fSliceL;
    BoostVecC fSliceR;
    double fAverageAbs2X;
    double fAverageAbs2Y;
    double fAverageAbs2A;
    double fAverageAbs2B;
    double fAverageAbs2L;
    double fAverageAbs2R;
    //double fRawStokesI; // unused. LN.
    //double fRawStokesQ; // unused. LN.
    //double fRawStokesU; // unused. LN.
    //double fRawStokesV; // unused. LN.
  };


  // The objects of this class contain the corrected polarization data for the signal
  class StokesParameters {

  public:
    StokesParameters(const BoostVecRangeC& signal_x, const BoostVecRangeC& signal_y,
                     const BoostVecRangeC& noise_x, const BoostVecRangeC& noise_y, bool subtract_noise);
    double GetStokesI() const { return fStokesI; }
    double GetStokesQ() const { return fStokesQ; }
    double GetStokesU() const { return fStokesU; }
    double GetStokesV() const { return fStokesV; }
    double GetStokesRatioQdivI() const { return fStokesRatioQdivI; }
    double GetStokesRatioUdivI() const { return fStokesRatioUdivI; }
    double GetStokesRatioVdivI() const { return fStokesRatioVdivI; }
    double GetPolarizationAngle() const { return 0.5*atan(fStokesU/fStokesQ); }
    double GetNoiseAdditionCovariance(const revt::StationRRecDataQuantities, revt::StationRRecDataQuantities);
    double GetAnalyticMethodVariance(const revt::StationRRecDataQuantities);
  private:
    // For private and internal internal use only (no pun intended), when noise levels have already been calculated
    StokesParameters(const BoostVecC& signal_x, const BoostVecC& signal_y, RawStokesParameters& noise, double multiply_noise);
    void NoiseAdditionCovarianceHelper(BoostVecD& quanities, const revt::StationRRecDataQuantities quantity_enum,
                                       const StokesParameters& sliding_stokes, SizeT index);
    void NoiseAdditionCovariancePolarizationAngleHelper(double& angle, double origval) const;

    // Sample Covariances are not to be confused with the (Co)variances between the stokes parameters
    Complex AnalyticMethodSampleCovarianceHelper(char channel1, char channel2, SizeT index);
    double AnalyticMethodOneChannelVarianceHelper(char channel);

    RawStokesParameters fSignal;
    RawStokesParameters fNoise;

    // These maps are for storing intermediate values of the analytic error estimation method
    SampleCovariancesT sampleCovariances;
    OneChannelPulseVariancesT oneChannelPulseVariances;

    double fStokesI;
    double fStokesQ;
    double fStokesU;
    double fStokesV;
    double fStokesRatioQdivI;
    double fStokesRatioUdivI;
    double fStokesRatioVdivI;
  };


  class RdPolarizationReconstructor : public fwk::VModule {

  public:
    RdPolarizationReconstructor();
    ~RdPolarizationReconstructor() override { }

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

  private:
    enum InfoLevel {
      eNone = 0,
      eFinal = 1,
      eIntermediate = 2,
      eObscure = 3,
      eMinuit = 4
    };

    int fInfoLevel;
    int fProjectOntoGroundPlane;
    bool fQuadraticNoiseSubtraction;
    //bool fCorrectForOpeningAngle; // unused. LN.
    bool fAssumeGroundPlaneIsHorizontalAndFlat;
    std::string fGeometrySource;
    std::string fErrorCalculationType;

    double fSignalWindowStartOverride;
    double fSignalWindowStopOverride;
    double fNoiseWindowStartOverride;
    double fNoiseWindowStopOverride;

    void GetRotatedPolarizationTraces(const evt::Event& event, const revt::Station& station,
                                      const rdet::RDetector& rDetector, const utl::TraceV3C& stationTimeSeries,
                                      BoostVecC& traceX, BoostVecC& traceY, double& observerAngle,
                                      double& alpha) const;
    void CalculateAndFillStokesParameters(revt::Station& station, BoostVecC& traceX,
                                          BoostVecC& traceY, double observerAngle,
                                          double alpha);
    void CheckSignalAndNoiseRanges(SizeT signal_start_bin, SizeT signal_stop_bin,
                                   SizeT noise_start_bin, SizeT noise_stop_bin,
                                   SizeT size) const;
    void CalculateAndFillChargeExcessFraction(revt::Station& station, double& alpha,
                                              double& chxfraction);
    utl::Vector GetArrivalDirection(const evt::Event& event) const;
    utl::CoordinateSystemPtr GetCS(const evt::Event& event) const;
    void GetCoreXY(const evt::Event& event, const utl::CoordinateSystemPtr& CS, double& CoreX, double& CoreY) const;
    utl::Point GetCore(const evt::Event& event) const;
    utl::Vector GetPerpVector(const utl::Point& point, const utl::Line& line) const;

    REGISTER_MODULE("RdPolarizationReconstructor", RdPolarizationReconstructor);

  };

}


#endif