IsLightning.cc

Go to the documentation of this file.

#include <sevt/Station.h>
#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <sevt/IsLightning.h>

#include <sstream>

using namespace std;
using namespace utl;


namespace sevt {

  namespace {
    const int kLightningThreshold = 1000;
    const int kLightningHysteresis = 4;
  }


  bool
  IsLightning(const Station& station)
  {
    if (!station.HasCalibData())
      return false;

    const auto& dStation = det::Detector::GetInstance().GetSDetector().GetStation(station);

    struct PMTInfo {
      const TraceI* trace;
      const TraceD* baseline;
    };

    PMTInfo activePMTs[dStation.GetNPMTs()];

    // find active PMTs and their baseline's
    int nActivePMTs = 0;

    for (const auto& pmt : station.PMTsRange()) {
      if (pmt.HasFADCTrace() && pmt.HasRecData()) {
        const auto& pmtRec = pmt.GetRecData();
        if (pmtRec.HasFADCBaseline(sdet::PMTConstants::eHighGain)) {
          activePMTs[nActivePMTs].trace = &(pmt.GetFADCTrace(sdet::PMTConstants::eHighGain));
          activePMTs[nActivePMTs].baseline = &(pmtRec.GetFADCBaseline(sdet::PMTConstants::eHighGain));
          ++nActivePMTs;
        }
      }
    }

    if (nActivePMTs == 3) {

      int count[3] = { 0 };
      bool cross[3] = { true, true, true };

      const int tracesize = dStation.GetFADCTraceLength();

      for (int pmt = 0; pmt < 3; ++pmt) {

        const PMTInfo& pmti = activePMTs[pmt];
        const TraceI& trace = *pmti.trace;
        const TraceD& baseline = *pmti.baseline;

        for (int i = 0; i < tracesize; ++i) {
          if (trace[i] > kLightningThreshold)
            return false;
          const double signal = trace[i] - baseline[i];
          if (cross[pmt] && signal < -kLightningHysteresis) {
            ++count[pmt];
            cross[pmt] = false;
          }
          if (!cross[pmt] && signal > kLightningHysteresis)
            cross[pmt] = true;
        }
      }

      // PMT 1 has smaller signal
      if (count[0] > 2 && count[1] > 3 && count[2] > 3) {
        ostringstream info;
        info << '(' << count[0] << ", " << count[1] << ", " << count[2]
             << ") oscillations in station " << station.GetId();
        INFO(info);
        return true;
      }
    }

    return false;
  }

}