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 "JDAQ/JDAQEventIO.hh"

#include "JAAnet/JAAnetToolkit.hh"

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

#include "JTrigger/JTimeConverter.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>());
    
    const JTimeConverter converter(*event, *tev);

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

      DEBUG("Hit (DAQ): " << setw(5) << pmt.getID() << ' ' << FIXED(6,1) << converter.getTime(*i, pmt) << endl);

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

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