Deprecated/RdStationSignalReconstructorLegacy/RdStationSignalReconstructor.h

Go to the documentation of this file.

#ifndef _RdStationSignalReconstructor_h_
#define _RdStationSignalReconstructor_h_

#include <fwk/VModule.h>

#include <revt/Station.h>
#include <revt/Channel.h>

#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>

#include <string>
#include <vector>

#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnUserParameterState.h"
#include <Minuit2/FCNBase.h>


namespace evt {
  class Event;
}

namespace revt {
  class Channel;
}

namespace RdStationSignalReconstructor {

  class RdStationSignalReconstructor : public fwk::VModule {

  public:
    RdStationSignalReconstructor() {}

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

  private:
    unsigned int fVectorialComponent = 0;
    bool fUseEnvelope = false;
    double fSignalIntegrationWindowSize = 0;
    
    std::string fNoiseSubtraction = "Undefined";
    std::string fSignalDef = "Undefined";
    std::string fNoiseDef = "Undefined";
    std::string fSignalToNoiseDef = "Undefined";
    std::string fPolarizationDef = "Undefined";
    
    bool fAddNoiseToUncertainty = false;

    double fMinSignal = 0;
    double fMinSignalToNoise = 0;
    double fRelativeAmplitudeUncertainty = 0;
    double fDefaultAntennaUncertainty = 0;
    double fButterflyAntennaUncertainty = 0;
    double fBlackSpiderAntennaUncertainty = 0;

    bool fSelfTriggeredSNcut = true;
    bool fAdjustSignalAmplitude = false;

    unsigned int fTotalNumberOfEvents = 0;
    unsigned int fNumberOfEvents1SignalStations = 0;
    unsigned int fNumberOfEvents2SignalStations = 0;
    unsigned int fNumberOfEvents3SignalStations = 0;
    unsigned int fNumberOfEvents4SignalStations = 0;
    unsigned int fNumberOfEvents5SignalStations = 0;
    unsigned int fNumberOfEvents5plusSignalStations = 0;

    int fInfoLevel = 0;
    bool fFitAnalyticSignal = 0;

    // calculates Mean and RMS in selected noise window
    void Noisefinder(const revt::ChannelTimeSeries& channeltrace,
                     double& RMSNoise, double& MeanNoise, double NoiseWindowStart, double NoiseWindowStop) const;

    // calculates the amplitude and the position of the maximum in the signal search window
    void Pulsefinder(const revt::ChannelTimeSeries& channeltrace,
                     double& PeakAmplitude, double& PeakTime, double& PeakTimeError,
                     double SignalSearchWindowStart, double SignalSearchWindowStop, unsigned int& sample) const;

    // calculates the FWHM and the integral of the signal in the FWHM
    void PulseFWHMIntegrator(const revt::ChannelTimeSeries& channeltrace,
                             unsigned int sample, const double PeakAmplitude, double& SignalFWHM, double& IntegratedSignal,
                             double& SignalWindowStart, double& SignalWindowStop) const;

    // calculates the integrated signal in a fixed window
    void PulseFixedWindowIntegrator(const revt::ChannelTimeSeries& channeltrace,
                                    unsigned int sample, double IntegrationTime, double& IntegratedSignal,
                                    double& SignalWindowStart, double& SignalWindowStop, const bool usePower) const;

    void PulseSlidingWindowIntegrator(const revt::ChannelTimeSeries& channeltrace,
                                      double IntegrationTime, double& IntegratedSignal,
                                      double& SignalWindowStart, double& SignalWindowStop,
                                      double SignalSearchWindowStart, double SignalSearchWindowStop) const;

    /* the formula is documented in GAP 2013-059
     * the input is the SNR defined as the maximum of the Hilbert envelope squared divided by the RMS of the Hilbert envelope squared
     */
    double GetSignalUncertainty(const double signalAmplitude, const double SNR) const;

    /* the formula is documented in PhD thesis of C. Glaser (RWTH Aachen University 2016)
     * the input is the SNR defined as the maximum of the Hilbert envelope squared divided by the RMS of the Hilbert envelope squared
     */
    double GetSignalTimeUncertainty(const double SNR) const;

    double GetAdjustedSignalAmplitude(const double signalAmpitude, const double SNR) const;

    ROOT::Minuit2::MnUserParameterState FitAnalyticSignal(utl::Trace<double>& trace,
                                                          const double noiseRMS,
                                                          revt::Station& /*station*/,
                                                          const int /*polarization*/,
                                                          const bool fixPhaseSlope,
                                                          const double phaseSlope);

    REGISTER_MODULE("RdStationSignalReconstructor", RdStationSignalReconstructor);

  };


  class SignalObjectiveFunction : public ROOT::Minuit2::FCNBase {

  public:
    SignalObjectiveFunction(utl::Trace<double> trace, const double noiseRMS) :
      fTrace(trace),
      fSigma(noiseRMS)
    { }

    double operator()(const std::vector<double>& pars) const;

    double Up() const { return fTheErrorDef; }

    void SetErrorDef(const double def) { fTheErrorDef = def; }

  private:
    utl::Trace<double> fTrace;
    double fSigma;
    double fTheErrorDef = 1;
 
  };

}


#endif