RdSimulationRadioTrigger.cc

Go to the documentation of this file.

#include "RdSimulationRadioTrigger.h"

#include <fwk/CentralConfig.h>

#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/TraceAlgorithm.h>
#include <utl/String.h>
#include <utl/Math.h>

#include <evt/Event.h>

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

#include <rdet/RDetector.h>

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>
#include <sevt/StationTriggerData.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <sdet/StationTriggerUtil.h>


using namespace std;
using namespace revt;
using namespace fwk;
using namespace utl;
using namespace det;
using namespace sdet::Trigger;


namespace RdSimulationRadioTrigger {

  VModule::ResultFlag
  RdSimulationRadioTrigger::Init()
  {
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdSimulationRadioTrigger");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("ThresholdMethod").GetData(fThresholdMethod);
    topBranch.GetChild("RadioT2Threshold").GetData(fRadioT2Threshold);
    topBranch.GetChild("RejectRFIThreshold").GetData(fRejectRFIThreshold);
    topBranch.GetChild("RejectRFIMultiplicity").GetData(fRejectRFIMultiplicity);

    ostringstream info;
    info << "\n\tUsing "<< (fThresholdMethod ? "Geometric sum" : "Maximum channel") << " as RD trigger method \n"
         << "\tThreshold for a radio T2 trigger: "<< fRadioT2Threshold << " ADC counts";
    INFOFinal(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdSimulationRadioTrigger::Run(evt::Event& event)
  {
    // Check if there are events at all
    if(!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    ostringstream msg;

    vector<int> RDT2Stations;
    REvent& rEvent = event.GetREvent();
    const auto& eventTime = event.GetHeader().GetTime();

    if(!event.HasSEvent()) {
      event.MakeSEvent();
    }

    sevt::SEvent& sEvent = event.GetSEvent();

    const Detector& Det = Detector::GetInstance();
    const rdet::RDetector& radioDet = Det.GetRDetector();

    // loop through stations and for every station through every channel
    for (auto& rStation : rEvent.StationsRange()) {

      std::pair<ChannelTimeSeries, ChannelTimeSeries> channelAdcTrace;

      for (const auto& channel : rStation.ChannelsRange()) {

        const unsigned int channelID = channel.GetId();

        // Checking if the channel of the processed station contribute to the event.
        // If not, continue with the next channel.
        if (!rStation.HasChannel(channelID) || !channel.IsActive())
          continue;

        if (channelID == 0)
          GetTraceWithoutDCOffset(channelAdcTrace.first, channel);
        else
          GetTraceWithoutDCOffset(channelAdcTrace.second, channel);
      }

      if (channelAdcTrace.first.GetSize() != channelAdcTrace.second.GetSize() ) {
        msg.str("");
        msg << "Station " << rStation.GetId()
            << " is rejected as its channel traces do not have the same length!";
        WARNING(msg);

        rStation.SetExcludedReason(eHardwareMismatch);
        continue;
      }

      // binning should be the same for both channels, just take Ch1
      const int binning = rStation.GetChannel(1).GetChannelADCTimeSeries().GetBinning();

      double peakAmplitude = 0;
      double peakTime = 0;
      int sample = 0;

      std::vector<unsigned int> t2Bins;
      std::vector<unsigned int> rfiBins;
      bool signalT2 = false;

      // calculate geometric sum of both channels for threshold trigger
      std::vector<double> finalAdcTrace;

      const int maxCh0 = *max_element(channelAdcTrace.first.Begin(),
            channelAdcTrace.first.End(), AbsValue);
      const int maxCh1 = *max_element(channelAdcTrace.second.Begin(),
            channelAdcTrace.second.End(), AbsValue);

      if (fThresholdMethod == 0) {  // maximum of both channels
        if (maxCh0 > maxCh1)
          std::copy(channelAdcTrace.first.Begin(), channelAdcTrace.first.End(),
                    std::back_inserter(finalAdcTrace));
        else
          std::copy(channelAdcTrace.second.Begin(), channelAdcTrace.second.End(),
                    std::back_inserter(finalAdcTrace));
      } else if (fThresholdMethod == 1) {  // geometric sum
        for (unsigned int i = 0; i < channelAdcTrace.first.GetSize() ; ++i) {
          finalAdcTrace.push_back(sqrt(Sqr(channelAdcTrace.first[i]) + Sqr(channelAdcTrace.second[i])));
        }
      } else {
        ERROR("Threshold method not existent.");
        return eFailure;
      }

      for (unsigned int i = 0; i < finalAdcTrace.size(); ++i) {
        if (std::abs(finalAdcTrace[i]) > fRadioT2Threshold)
          t2Bins.push_back(i);
        if (std::abs(finalAdcTrace[i]) > fRejectRFIThreshold)
          rfiBins.push_back(i);
      }

      // if no trigger just move on with next station
      if (t2Bins.empty()) {
        msg.str("");
        msg << "Station " << rStation.GetId()
            << " has no T2 trigger (not above threshold).";
        INFODebug(msg);
        continue;
      }

      /* see if it is really a signal T2 or likely RFI.
         For that check that there are not too many bins over some threshold within 500ns
         after 1mus after the trigger */
      for (auto const& t2bin: t2Bins) {
        const int cnt = std::count_if(std::begin(rfiBins), std::end(rfiBins),
                        [&t2bin](unsigned int rfibin){return (rfibin >= t2bin+25 && rfibin <= t2bin+275);});
        if (cnt <= fRejectRFIMultiplicity) {
          peakAmplitude = finalAdcTrace[t2bin];
          peakTime = t2bin * binning / 1000;
          sample = t2bin;
          signalT2 = true;
          break;
        }
      }

      if (!signalT2) {
        msg.str("");
        msg << "Station " << rStation.GetId()
            << " is rejected as it is a RFI trace!";
        WARNING(msg);

        rStation.SetRejectedReason(eNoisy);
        continue;
      }

      // here comes the SEvent part. Pushing it to the SEvent and make WCD traces
      const int sStID = rStation.GetId() - radioDet.GetRdSdStationIdLink();
      RDT2Stations.push_back(sStID);

      if (!sEvent.HasStation(sStID))
        sEvent.MakeStation(sStID);

      sevt::Station& sStation = sEvent.GetStation(sStID);

      if (!sStation.HasSimData())
        sStation.MakeSimData();

      // okay now pushing some stuff into RD
      // this is ugly, I know, but until now we do not know where the time difference
      // between particles and radio comes from (electronics, physics or both) and how it behaves
      // in simulations vs data. When we investigated that, we should change this. But it works here
      rStation.GetTriggerData().SetSelfTriggerTime(
        (peakTime + rStation.GetRecData().GetParameter(eTraceStartTime)/1000 + 13.2) * 1000 * nanosecond);
      rStation.GetTriggerData().SetTriggerSource(StationTriggerData::eSelf);
      rStation.GetTriggerData().SetSelfTriggerType(StationTriggerData::eT2Threshold);
      rStation.GetTriggerData().SetSelfTriggerBin(sample);
      rStation.GetTriggerData().SetSelfTriggerAmp(peakAmplitude);

      vector<const TimeDistributionI*> traces;

      int startBin = 0;
      int stopBin = 0;
      PrepareTraces(sStation, startBin, stopBin, traces);

      if (traces.size() < 3)
        sStation.SetRejected(sevt::StationConstants::eNoCalibData);

      const auto& dStation = Detector::GetInstance().GetSDetector().GetStation(sStation);
      const int WCDTrigBin = rStation.GetTriggerData().GetSelfTriggerTime() / dStation.GetFADCBinSize();
      const int traceStart = std::max(WCDTrigBin - (int)dStation.GetFADCTraceLength(), startBin);
      Buffer(sStation, eventTime, sevt::StationTriggerData::eRDThreshold, sevt::StationTriggerData::ePLDRD,
             traceStart, WCDTrigBin);
    }

    if (!RDT2Stations.empty()) {
      ostringstream info;
      info << "\nAdding stations:" << '\n' << String::OfSortedIds(RDT2Stations) << " as RD T2s";
      INFO(info);
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdSimulationRadioTrigger::Finish()
  {
    ostringstream info;
    info << "\n\tThere were " << fNbOfSaturatedChannels << " saturated channels!";
    INFO(info);

    return eSuccess;
  }


  void
  RdSimulationRadioTrigger::GetTraceWithoutDCOffset(ChannelTimeSeries& currTrace, const Channel& channel)
  {
    const auto adcTrace = channel.GetChannelADCTimeSeries();
    const double mean = TraceAlgorithm::Mean(adcTrace, 0, adcTrace.GetSize() - 1);
    for (auto& adcSignal : adcTrace)
      currTrace.PushBack((double)adcSignal - mean);
  }


  bool
  RdSimulationRadioTrigger::AbsValue(double a , double b)
  {
    return abs(a) < abs(b);
  }

}