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;
 
 
  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)));
          }
        }
      }
 
      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;
}