DoublePeakDetector.cc

Go to the documentation of this file.

#include "DoublePeakDetector.h"

// Offline header
#include <fwk/CentralConfig.h>

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

#include <evt/Event.h>

#include <sevt/EventTrigger.h>
#include <sevt/Header.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/PMT.h>
#include <sevt/PMTCalibData.h>
#include <sevt/PMTRecData.h>
#include <sevt/StationRecData.h>
#include <sevt/StationTriggerData.h>
#include <sevt/StationCalibData.h>

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


#include <iostream>
#include <sstream>
#include <algorithm>


#include <config.h>

#include <AugerEvent.h>
#include <IoSdData.h>

using namespace fwk;
using namespace evt;
using namespace sevt;
using namespace utl;
using namespace std;

using namespace DoublePeakDetectorNS;

const double kBinSize = 25;
const double kBinThr = 0.02;
const double kPeakThr = 3;
const int kPeakWind = 10;
const int kGapWind = 15;
const double kMinScoreDiff = 0.15;

int DoublePeakDetector::FindPeaks(const TraceD &trace, vector<signalBlock_t> &signalBlocks, const int nTraceBins)
{
  for (int i = 0; i < nTraceBins - kPeakWind; i++)
  {
    double binValue = trace[i];
    if (binValue < kBinThr) continue;

    double winSum = binValue;
    double winPeak = binValue;
    int nBinsOverThr = 0;
    int blockStart = i - 1;
    int blockEnd = i;
    int coreBlockEnd = i;
    bool blockFound = false;

    for (int j = i + 1; j < i + kPeakWind; j++) {
      binValue = trace[j];
      if (binValue < kBinThr) continue;
      winSum += binValue;
      nBinsOverThr++;

      if (winPeak < binValue) winPeak = binValue;

      if (winSum > kPeakThr) {
        coreBlockEnd = j;
        blockFound = true;
        break;
      }
    }
    if (nBinsOverThr == kPeakWind) {
      coreBlockEnd = i + kPeakWind - 1;
      blockFound = true;
    }

    if (!blockFound) continue;
    
    blockEnd = coreBlockEnd;
    signalBlock_t sigBlock;
    sigBlock.start = blockStart * kBinSize;
    // Add bins to the block
    int nEmptyBins = 0;
    for (int j = coreBlockEnd + 1; j < coreBlockEnd + nTraceBins; j++)
    {
      binValue = trace[j];

      if (binValue < kBinThr) {
        nEmptyBins++;
        if (nEmptyBins == kGapWind) {
          break;
        }
      }
      else {
        nEmptyBins = 0;
        blockEnd = j;
        nBinsOverThr++;
        winSum += binValue;
      }
    }

    sigBlock.end = blockEnd * kBinSize;

    sigBlock.score = winSum * nBinsOverThr;
    signalBlocks.push_back(sigBlock);

    i = blockEnd;
  }

  return signalBlocks.size();
}

bool DoublePeakDetector::ComputeIsMultiple(const vector<signalBlock_t> &signalBlocks)
{

  unsigned int nBlocks = signalBlocks.size();
  if (nBlocks < 2) return false;

  double maxScore = -1;
  for (unsigned int i = 0; i < nBlocks; i++) {
    double score = signalBlocks[i].score;
    if (maxScore < score) {
      maxScore = score;
    }
  }

  if (maxScore < 0) {
    std::ostringstream error;
    error << "\nError in trace!\n\n";
    ERROR(error);
    // throw -1000;
    //return eFailure;
  }

  unsigned int nSimilar = 0;
  for (unsigned int i = 0; i < nBlocks; i++) {
    double score = signalBlocks[i].score;
    if (score / maxScore > kMinScoreDiff) {
      nSimilar++;
      if (nSimilar == 2)
      {
        return true;
      }
    }
  }

  return false;
}

VModule::ResultFlag DoublePeakDetector::Init()
{
  Branch topB = CentralConfig::GetInstance()->GetTopBranch("DoublePeakDetector");
  // TFile outFile("traceHistos.root","RECREATE");
  // outFile.Close();
  return eSuccess;
}

VModule::ResultFlag DoublePeakDetector::Run(evt::Event &event)
{

  const auto& sDet = det::Detector::GetInstance().GetSDetector();
  sevt::SEvent &sEvent = event.GetSEvent();

  // TFile outFile("traceHistos.root","UPDATE");
  for (auto& sStation : sEvent.StationsRange()) {
    if (sStation.IsCandidate() && sDet.GetStation(sStation.GetId()).IsInGrid()) {
      const StationRecData &stRec = sStation.GetRecData();
      const TraceD &trace = sStation.GetVEMTrace();
      ostringstream os;
      os << "Event:" << sEvent.GetHeader().GetId() << " Station: " << sStation.GetId();
      
      vector<signalBlock_t> signalBlocks;
      double binning = trace.GetBinning();
      int start = stRec.GetSignalStartSlot();
      int end = stRec.GetSignalEndSlot() + 1;
      const int traceLength = end - start;
      const int nTraceBins = traceLength * binning;

      FindPeaks(trace, signalBlocks, nTraceBins);
      bool isMultiple = ComputeIsMultiple(signalBlocks);

      if (isMultiple)
      {
        std::ostringstream info;
        info << "Deleting: " << os.str() << "\n";
        INFO(info);
        StationTriggerData &trig = sStation.GetTriggerData();
        trig.SetAlgorithm(StationTriggerData::eRandom);
        sStation.SetRejected(StationConstants::eRandom);
        trig.SetErrorCode(5001);
      }


      // Yann: I don't know where to select that a station has at least 2 active PMTs
      // I put it here but I should ask.
      // This should be something done by default
      int pmtCounter = 0;
      for (auto& pmt : sStation.PMTsRange()) {
        if (pmt.HasCalibData() && pmt.GetCalibData().IsTubeOk()) {
          pmtCounter++;
        }
      }
      if (pmtCounter < 2) {
        std::ostringstream info;
        info << "Deleting: " << os.str() << "\n"
             << "Too few PMTs"<< "\n";
        INFO(info);
        StationTriggerData &trig = sStation.GetTriggerData();
        trig.SetAlgorithm(StationTriggerData::eRandom);
        sStation.SetRejected(StationConstants::eRandom);
        trig.SetErrorCode(5001);
      }
      else if (pmtCounter == 2 || pmtCounter == 3){}
      else {
        return eFailure;
      }
    }
  }

  // outFile.Close();

  return eSuccess;
}

VModule::ResultFlag DoublePeakDetector::Finish()
{
  return eSuccess;
}

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. DoublePeakDetector/DoublePeakDetector.o -k"
// End: