Jaws.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <bitset>

#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQPMTIdentifier.hh"

#include "JDAQ/JDAQTimesliceIO.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JTools/JRange.hh"
#include "JSystem/JKeypress.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {

  using KM3NETDAQ::JDAQHit;

  /**
   * Auxiliary class to search for periodic hits.
   */
  struct JPeriod {
    /**
     * Default constructor.
     */
    JPeriod() :
      period(-1),
      margin(-1)
    {}      


    /**
     * Constructor.
     *
     * \param  period           period [ns]
     * \param  margin           margin [ns]
     */
    JPeriod(int period,
            int margin) :
      period(period),
      margin(margin)
    {}

    /**
     * Check compliance with periodicity.
     *
     * \param  first            first  hit
     * \param  second           second hit
     * \return                  true if periodic; else false
     */
    bool operator()(const JDAQHit& first, const JDAQHit& second) const
    {
      const int t1 = (second.getT() - first.getT()) % period;

      return (t1 <= margin || t1 >= period - margin);
    }

    /**
     * Read periodicity from input stream.
     *
     * \param  in               input stream
     * \param  object           periodicity
     * \return                  input stream
     */
    friend std::istream& operator>>(std::istream& in, JPeriod& object)
    {
      return in >> object.period >> object.margin;
    }

    /**
     * Write periodicity to output stream.
     *
     * \param  out              output stream
     * \param  object           periodicity
     * \return                  output stream
     */
    friend std::ostream& operator<<(std::ostream& out, const JPeriod& object)
    {
      return out << object.period << ' ' << object.margin;
    }

    int period;
    int margin;
  };
}


/**
 * \file
 * Example program to search for spurious hits.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JRange<int>  JRange_t; 

  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  JROOTClassSelector     selector;
  JRange_t               ToT_ns;
  JPeriod                is_periodic;
  int                    numberOfHits;
  int                    numberOfColumns;
  set<JDAQPMTIdentifier> veto;
  int                    debug;

  try { 

    JParser<> zap("Example program to search for spurious hits.");
    
    zap['f'] = make_field(inputFile,         "input file).");
    zap['n'] = make_field(numberOfEvents)    = JLimit::max(); 
    zap['C'] = make_field(selector,          "timeslice selector, e.g. JDAQTimeslice.")   = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['T'] = make_field(ToT_ns,            "time-over-threshold range [ns]");
    zap['p'] = make_field(is_periodic,       "periodicity");
    zap['M'] = make_field(numberOfHits,      "minimum number of hits for periodicity");
    zap['N'] = make_field(numberOfColumns,   "number of columns for printing")            = 5;
    zap['V'] = make_field(veto,              "veto PMT from analysis")                    = JPARSER::initialised();
    zap['d'] = make_field(debug,             "debug.")                                    = 3;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);

  for (counter_type counter = 0; in.hasNext() && counter != inputFile.getLimit(); ++counter) {

    STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);

    const JDAQTimeslice* timeslice = in.next();

    for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {

      // organise data per PMT

      typedef vector<JDAQHit>        buffer1D_type;
      typedef vector<buffer1D_type>  buffer2D_type;

      buffer2D_type buffer(numeric_limits<JDAQHit::JPMT_t>::max(), buffer1D_type());

      for (JDAQFrame::const_iterator hit = frame->begin(); hit != frame->end(); ++hit) {
        buffer[hit->getPMT()].push_back(*hit);
      }

      // detection of perodic signal

      map<JDAQPMTIdentifier, vector<JDAQHit> > top;

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {

        if (!buffer[pmt].empty()) {

          const JDAQPMTIdentifier id(frame->getModuleID(), pmt);

          for (buffer1D_type::iterator i = buffer[pmt].begin(); i != buffer[pmt].end(); ) {

            if (ToT_ns(i->getToT())) {

              buffer1D_type::iterator p = i;

              for (++p; p != buffer[pmt].end() && ToT_ns(p->getToT()) && is_periodic(*i, *p); ++p) {}

              if (distance(i, p) >= numberOfHits) {
                top[id].push_back(*i);
              }

              i = p;

            } else {

              ++i;
            }
          }
        }
      }


      if (debug >= debug_t) {

        for (map<JDAQPMTIdentifier, vector<JDAQHit> >::const_iterator i = top.begin(); i != top.end(); ++i) {

          if (veto.count(i->first) == 0) {

            const int pmt = i->first.getPMTAddress();

            cout << timeslice->getRunNumber()         << '/' << FILL(6,'0') << timeslice->getFrameIndex() << FILL() << ' '
                 << setw(9) << i->first.getModuleID() << '.' << FILL(2,'0') << i->first.getPMTAddress()   << FILL() << ' '
                 << setw(1) << frame->testHighRateVeto(pmt) 
                 << setw(1) << frame->testFIFOStatus  (pmt)  << endl;

            int count = 0;

            string eol = "";

            for (JDAQFrame::const_iterator hit = frame->begin(); hit != frame->end(); ++hit) {

              const bool has_pmt  =  (hit->getPMT() == i->first.getPMTAddress());
              const bool has_hit  =  (find(i->second.begin(), i->second.end(), *hit) != i->second.end());

              if (has_pmt) {

                if (has_hit) {
                  eol = "*";
                }
                
                cout << setw(10) << hit->getT() << ' ' << setw(3) << (int) hit->getToT() << ' ' << RESET;

                if ((++count % numberOfColumns) == 0) {

                  cout << eol << endl;

                  eol = "";
                }
              }
            }

            cout << endl << "press any key to continue> " << flush;
          
            JKeypress(false).get();

            cout << endl;
          }
        }
      }
    }
  }
  STATUS(endl);
}