JFeedback.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <set>
 
#include "TROOT.h"
#include "TFile.h"
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/weightlist.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
 
#include "flux/Flux.hh"
 
#include "JAAnet/JHead.hh"
#include "JLang/JManip.hh"
 
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  using JLANG::JObjectIterator;
 
  const Flux_Atmospheric  conventional;
 
 
  /**
   * Auxiliary data structure for run statistics.
   */
  struct stats_t {
 
    std::string name;
    double events     = 0;      //!< number of events
    double range_s    = 0.0;    //!< time difference between first and last event
    double livetime_s = 0.0;    //!< livetime of events 
    double daq_s      = 0.0;    //!< livetime of data taking run
    double gap_s      = 0.0;    //!< maximal time between two consecutive events
 
    friend inline std::ostream& operator<<(std::ostream& out, const stats_t& object)
    {
      using namespace std;
 
      out << object.name  << endl;
      out << "events:   " << FIXED(9,1) << object.events               << endl;
      out << "range:    " << FIXED(7,1) << object.range_s    << " [s]" << endl;
      out << "livetime: " << FIXED(7,1) << object.livetime_s << " [s]" << endl;
      out << "DAQ:      " << FIXED(7,1) << object.daq_s      << " [s]" << endl;
      out << "gap:      " << FIXED(7,1) << object.gap_s      << " [s]" << endl;
 
      return out;
    }
  };
 
 
  /**
   * Get run statistics for real data.
   *
   * \param  head                    header
   * \param  in                      event input
   * \return                         statistics
   */
  inline stats_t getStatsFromData(const JHead& head, JObjectIterator<Evt>& in)
  {
    using namespace std;
    using namespace JPP;
    using namespace KM3NETDAQ;
 
    stats_t object;
 
    set<int> buffer;
    double   events = 0;
    
    while (in.hasNext()) {
 
      const Evt* evt = in.next();
 
      buffer.insert(evt->frame_index);
 
      events += 1;
    }
 
    if (buffer.size() < 2) {
      FATAL("Insufficient events " << buffer.size() << endl);
    }
 
    int gap = 0;
    
    for (set<int>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {
      if (*q - *p > gap) {
        gap = *q - *p;
      }
    }
 
    object.name       = "data";
    object.events     = events;
    object.range_s    = (*buffer.rbegin() - *buffer.begin()) * getFrameTime() * 1.0e-9;
    object.livetime_s =   buffer.size()                      * getFrameTime() * 1.0e-9;
    object.daq_s      = head.DAQ.livetime_s;
    object.gap_s      = gap                                  * getFrameTime() * 1.0e-9;
 
    return  object;
  }
 
 
  /**
   * Get run statistics for MUPAGE.
   *
   * \param  head                    header
   * \param  in                      event input
   * \return                         statistics
   */
  inline stats_t getStatsFromMUPAGE(const JHead& head, JObjectIterator<Evt>& in)
  {
    using namespace std;
    using namespace JPP;
 
    stats_t object;
 
    set<double> buffer;
    double      events = 0;
    
    while (in.hasNext()) {
 
      const Evt* evt = in.next();
 
      buffer.insert(evt->mc_event_time.AsDouble());
 
      events += 1;
    }
 
    if (buffer.size() < 2) {
      FATAL("Insufficient events " << buffer.size() << endl);
    }
 
    double gap = 0;
    
    for (set<double>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {
      if (*q - *p > gap) {
        gap = *q - *p;
      }
    }
 
    object.name       = "MUPAGE";
    object.events     = events * head.DAQ.livetime_s / head.livetime.numberOfSeconds;
    object.range_s    = (*buffer.rbegin() - *buffer.begin());
    object.livetime_s = head.livetime.numberOfSeconds;
    object.daq_s      = head.DAQ.livetime_s;
    object.gap_s      = gap;
 
    return  object;
  }
 
 
  /**
   * Get run statistics for gSeaGen.
   *
   * \param  head                    header
   * \param  in                      event input
   * \return                         statistics
   */
  inline stats_t getStatsFromGSEAGEN(const JHead& head, JObjectIterator<Evt>& in)
  {
    using namespace std;
    using namespace JPP;
 
    stats_t object;
 
    set<double> buffer;
    double      events = 0;
    
    while (in.hasNext()) {
 
      const Evt* evt = in.next();
 
      buffer.insert(evt->mc_event_time.AsDouble());
 
      const Trk trk = get_neutrino(*evt);
 
      events += evt->w[WEIGHTLIST_DIFFERENTIAL_EVENT_RATE] * evt->w[WEIGHTLIST_NORMALISATION] * conventional.dNdEdOmega(trk.type, log10(trk.E), trk.dir.z);
    }
 
    if (buffer.size() < 2) {
      FATAL("Insufficient events " << buffer.size() << endl);
    }
 
    double gap = 0;
 
    for (set<double>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {
      if (*q - *p > gap) {
        gap = *q - *p;
      }
    }
 
    object.name       = "gSeaGen";
    object.events     = events * head.DAQ.livetime_s;
    object.range_s    = (*buffer.rbegin() - *buffer.begin());
    object.livetime_s = 0.0;
    object.daq_s      = head.DAQ.livetime_s;
    object.gap_s      = gap;
 
    return  object;
  }
 
 
  /**
   * Get run statistics for pure-noise.
   *
   * \param  head                    header
   * \param  in                      event input
   * \return                         statistics
   */
  inline stats_t getStatsFromK40(const JHead& head, JObjectIterator<Evt>& in)
  {
    using namespace std;
    using namespace JPP;
 
    stats_t object;
 
    set<int> buffer;
    double   events = 0;
    
    while (in.hasNext()) {
 
      const Evt* evt = in.next();
 
      buffer.insert(evt->frame_index);
 
      events += 1;
    }
 
    if (buffer.size() < 2) {
      FATAL("Insufficient events " << buffer.size() << endl);
    }
 
    int gap = 0;
    
    for (set<int>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {
      if (*q - *p > gap) {
        gap = *q - *p;
      }
    }
 
    object.name       = "noise";
    object.events     = events * head.DAQ.livetime_s / head.K40.livetime_s;
    object.range_s    = (*buffer.rbegin() - *buffer.begin()) * getFrameTime() * 1.0e-9;
    object.livetime_s = head.K40.livetime_s;
    object.daq_s      = head.DAQ.livetime_s;
    object.gap_s      = gap                                  * getFrameTime() * 1.0e-9;
 
    return  object;
  }
 
 
  /**
   * Get run statistics.
   *
   * \param  head                    header
   * \param  in                      event input
   * \return                         statistics
   */
  inline stats_t getStats(const JHead& head, JObjectIterator<Evt>& in)
  {
    using namespace std;
    using namespace JPP;
 
    const bool mupage  = is_mupage(head);
    const bool gseagen = !mupage && head.is_valid(&JHead::depth) && (head.calibration.buffer == calibration::statical());
    const bool k40     = head.is_valid(&JHead::K40);
 
    if      (mupage)
      return getStatsFromMUPAGE (head, in);
    else if (gseagen)
      return getStatsFromGSEAGEN(head, in);
    else if (k40)
      return getStatsFromK40    (head, in);
    else
      return getStatsFromData   (head, in);
  }
}
 
 
/**
 * \file
 *
 * Example program to analyse track fit results from Evt formatted data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  JSingleFileScanner<Evt> inputFile;
  JLimit_t&  numberOfEvents = inputFile.getLimit();
  int        debug;
 
  try {
 
    JParser<> zap("Example program to analyse track fit results from Evt formatted data.");
 
    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['d'] = make_field(debug)               =  2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  const JHead   head  = getHeader(inputFile);
  const stats_t stats = getStats(head, inputFile);
 
  DEBUG(inputFile << endl << head << endl << stats << endl);
 
  cout << setw(132) << left << inputFile << right << ' '
       << FIXED(9,1) << stats.events              << ' '
       << FIXED(7,1) << stats.range_s             << ' '
       << FIXED(7,1) << stats.livetime_s          << ' '
       << FIXED(7,1) << stats.daq_s               << ' '
       << FIXED(7,1) << stats.gap_s               << endl;
}