JTimeConverter.cc

Go to the documentation of this file.

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

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Hit.hh"
#include "km3net-dataformat/offline/io_online.hh"
#include "km3net-dataformat/tools/time_converter.hh"
#include "JDAQ/JDAQEventIO.hh"

#include "JAAnet/JAAnetToolkit.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"

#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

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


namespace {

  /**
   * Count number of coincident hits.
   *
   * \param  data1       first  data set
   * \param  data2       second data set
   * \param  Tmax_ns     time window [ns]
   * \return             number of hits
   */
  inline int get_count(const std::set<double>& data1,
                       const std::set<double>& data2,
                       const double Tmax_ns)
  {
    int number_of_hits = 0;

    for (std::set<double>::const_iterator p = data1.begin(), q = data2.begin(); ; ++p, ++q) {

      while (p != data1.end() && q != data2.end() && *p < *q - Tmax_ns) { ++p; }
      while (p != data1.end() && q != data2.end() && *q < *p - Tmax_ns) { ++q; }

      if (p == data1.end() || q == data2.end()) {
        break;
      }
      
      if (fabs(*p - *q) <= Tmax_ns) {
        ++number_of_hits;
      }
    }

    return number_of_hits;
  }
}


/**
 * \file
 *
 * Example program to test conversion between Monte Carlo and DAQ times.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JTriggeredFileScanner<> inputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         detectorFile;
  double         Tmax_ns;
  int            debug;

  try { 

    JParser<> zap("Example program to test conversion between Monte Carlo and DAQ times.");
    
    zap['f'] = make_field(inputFile);
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['T'] = make_field(Tmax_ns)             = 20.0;
    zap['d'] = make_field(debug)               = 3;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace KM3NETDAQ;

  cout.tie(&cerr);


  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  const JModuleRouter router(detector);


  //typedef JDAQSnapshotHit JHit_t;
  typedef JDAQTriggeredHit JHit_t;

  while (inputFile.hasNext()) {

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

    JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();

    const JDAQEvent* tev   = ps;
    const Evt*       event = ps;

    ASSERT(!tev->empty<JHit_t>());

    typedef map<int, set<double> > map_type;

    map_type mc;
    map_type data;
    
    for (vector<Hit>::const_iterator i = event->mc_hits.begin(); i != event->mc_hits.end(); ++i) {

      if (!is_noise(*i)) {

        DEBUG("Hit (MC):  " << setw(5) << i->pmt_id << ' ' << FIXED(6,1) << getTime(*i) << endl);

        mc[i->pmt_id].insert(getTime(*i));
      }
    }

    for (int test : {1, 2} ) {

      DEBUG("Test "<< test << endl);

      data.clear();

      time_converter converter;

      if        (test == 1) {

        converter = time_converter(*event, *tev);

        for (JDAQEvent::const_iterator<JHit_t> i = tev->begin<JHit_t>(); i != tev->end<JHit_t>(); ++i) {

          const JModule& module = router.getModule(i->getModuleID());
          const JPMT&    pmt    = module.getPMT(i->getPMT());

          data[pmt.getID()].insert(converter.getTime(getTime(*i, pmt)));
        }

      } else if (test == 2) {

        Evt evt = *event;

        read(evt, *tev);
        
        converter = time_converter(evt);

        for (vector<Hit>::const_iterator i = evt.hits.begin(); i != evt.hits.end(); ++i) {

          if (i->trig != 0) {

            const JModule& module = router.getModule(i->dom_id);
            const JPMT&    pmt    = module.getPMT(i->channel_id);

            data[pmt.getID()].insert(converter.getTime(getTime(i->tdc, pmt)));
            //data[pmt.getID()].insert(converter.getTime(getTime(i->tdc, pmt)));
          }
        }
      }

      for (map_type::const_iterator i = data.begin(); i != data.end(); ++i) {
        for (set<double>::const_iterator hit = i->second.begin(); hit != i->second.end(); ++hit) {
          DEBUG("Hit (DAQ): " << setw(5) << i->first << ' ' << FIXED(6,1) << *hit << endl);
        }
      }

      int number_of_hits = 0;

      for (map_type::const_iterator p = mc.begin(), q = data.begin(); ; ++p, ++q) {

        while (p != mc.end() && q != data.end() && p->first < q->first) { ++p; }
        while (p != mc.end() && q != data.end() && q->first < p->first) { ++q; }

        if (p == mc.end() || q == data.end()) {
          break;
        }
      
        if (p->first == q->first) {
          number_of_hits += get_count(p->second, q->second, Tmax_ns);
        }
      }

      NOTICE("Number of hits " << setw(5) << number_of_hits << '/' << setw(5) << tev->size<JHit_t>() << endl);
    
      ASSERT(number_of_hits != 0);
    }
  }

  ASSERT(inputFile.getCounter() != 0);
  
  return 0;
}