JRemovePMT.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <set>

#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQPMTIdentifier.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JTrigger/JTriggerParameters.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"

#include "JLang/JTypeSelector.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JLang/JObjectDemultiplexer.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

namespace {

  using KM3NETDAQ::JDAQPMTIdentifier;

  /**
   * Remove PMT(s) from data.
   *
   * \param  data        data (I/O)
   * \param  PMT         PMT(s) to be removed
   */
  template<class T>
  inline void remove(std::vector<T>& data, const std::set<JDAQPMTIdentifier>& PMT)
  {
    using namespace std;

    vector<T> buffer;

    for (typename vector<T>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
      if (PMT.find(JDAQPMTIdentifier(hit->getModuleID(), hit->getPMT())) == PMT.end() &&
          PMT.find(JDAQPMTIdentifier(-1,                 hit->getPMT())) == PMT.end()) {
        buffer.push_back(*hit);
      }
    }

    data.swap(buffer);
  }
}


/**
 * \file
 *
 * Example program to remove PMT(s) from data (and set corresponding rate to 0).
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
  JLimit_t&                       numberOfEvents = inputFile.getLimit();
  JFileRecorder<JTriggerTypes_t>  outputFile;
  set<JDAQPMTIdentifier>          PMT;
  int                             debug;

  try { 

    JParser<> zap("Example program to remove PMT(s) from data (and set corresponding rate to 0).");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "abc.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['P'] = make_field(PMT)                 = JPARSER::initialised();
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  cout.tie(&cerr);

  outputFile.open();
  
  {
    JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile);

    // put base class JDAQTimeslice at end of type list

    typedef JAppend<JRemove<JDAQTimesliceTypes_t, JDAQTimeslice>::typelist, JDAQTimeslice>::typelist daq_timeslicetypes_t;

    JObjectDemultiplexer<JDAQTimeslice, daq_timeslicetypes_t> out(outputFile);

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

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

      JDAQTimeslice* timeslice = in.next();
    
      for (JDAQTimeslice::iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {

        bool rm = false;

        for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
          
          if (PMT.find(JDAQPMTIdentifier(frame->getModuleID(), pmt)) != PMT.end() ||
              PMT.find(JDAQPMTIdentifier(-1,                   pmt)) != PMT.end()) {
            rm = true;
          }
        }

        if (rm) {

          vector<JDAQHit> buffer[NUMBER_OF_PMTS];

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

          for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
          
            if (PMT.find(JDAQPMTIdentifier(frame->getModuleID(), pmt)) == PMT.end() &&
                PMT.find(JDAQPMTIdentifier(-1,                   pmt)) == PMT.end()) {
              copy(buffer[pmt].begin(), buffer[pmt].end(), back_inserter(data));
            }
          }
          
          sort(data.begin(), data.end());
        
          frame->clear();
          frame->add(data.size(), data.data());
        }
      }
      
      out.put(*timeslice);
    }
    STATUS(endl);
  }
  {
    JMultipleFileScanner<JDAQEvent> in(inputFile);

    while (in.hasNext()) {

      STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);

      JDAQEvent* event = in.next();

      remove(event->getHits<JDAQTriggeredHit>(), PMT);
      remove(event->getHits<JDAQSnapshotHit> (), PMT);

      outputFile.put(*event);
    }
    STATUS(endl);
  }
  {
    JMultipleFileScanner<JDAQSummaryslice> in(inputFile);

    while (in.hasNext()) {

      STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);

      JDAQSummaryslice* summaryslice = in.next();

      for (JDAQSummaryslice::iterator frame = summaryslice->begin(); frame != summaryslice->end(); ++frame) {

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

          if (PMT.find(JDAQPMTIdentifier(frame->getModuleID(), pmt)) != PMT.end() ||
              PMT.find(JDAQPMTIdentifier(-1,                   pmt)) != PMT.end()) {
            (*frame)[pmt].setValue(0.0);
          }
        }
      }

      outputFile.put(*summaryslice);
    }
    STATUS(endl);
  }

  JMultipleFileScanner<JTriggerParameters> io(inputFile);

  io >> outputFile;

  outputFile.close();
}