JEvtMerge.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <deque>
#include <queue>
#include <algorithm>
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Evt.hh"
 
#include "JLang/JException.hh"
 
#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  /**
   * Auxiliary data structure for priority queue / heap.
   */
  struct entry_type {
    /**
     * Prioritise two entries.
     *
     * \param  first           first  entry
     * \param  second          second entry
     * \return                 true if first entry higher priority than second; else false
     */
    friend inline bool operator<(const entry_type& first, const entry_type& second)
    {
      return first.evt.mc_event_time.AsDouble() > second.evt.mc_event_time.AsDouble();
    }
 
    size_t index;
    Evt    evt;
  };
}
 
 
/**
 * \file
 * Auxiliary program for time sorting of Monte-Carlo events.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef typename JTYPELIST<Head, MultiHead, Evt, JMeta>::typelist  typelist;
  
  vector<string>          inputFile;
  JFileRecorder<typelist> outputFile;
  int                     debug;
 
  try {
    
    JParser<> zap("Auxiliary program for time sorting of Monte-Carlo events.");
    
    zap['f'] = make_field(inputFile,  "Input file (produced with JEvtSplit)");
    zap['o'] = make_field(outputFile);
    zap['d'] = make_field(debug)               = 2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  outputFile.open();
    
  outputFile.put(JMeta(argc, argv));
 
  const size_t N = inputFile.size();
 
  typedef JTreeScanner      <Evt>  JScanner_t;
  //typedef JSingleFileScanner<Evt>  JScanner_t;
 
  counter_type number_of_events = JTreeScanner<Evt>(inputFile).getEntries();
 
  vector<JScanner_t> in(N);
 
  for (size_t i = 0; i != N; ++i) {
 
    in[i].configure(inputFile[i]);
 
    STATUS("Processing " << inputFile[i] << endl);
  }
 
#ifdef PRIORITY_QUEUE
 
  priority_queue<entry_type, deque<entry_type> > buffer;    // multiplex
 
  for (size_t i = 0; i != N; ++i) {
 
    if (in[i].hasNext()) {
 
      buffer.push({i, *in[i].next()});
 
      --number_of_events;
    }
  }
 
  while (number_of_events != 0) {                           // process
 
    STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
 
    const entry_type& top = buffer.top();
 
    outputFile.put(top.evt);
 
    if (in[top.index].hasNext()) {
 
      buffer.push({top.index, *in[top.index].next()});
 
      --number_of_events;
    }
      
    buffer.pop();
  }
 
  while (!buffer.empty()) {                                 // flush
 
    const entry_type& top = buffer.top();
 
    outputFile.put(top.evt);
      
    buffer.pop();
  }
 
#else
 
  vector<entry_type> buffer;                                // multiplex
 
  for (size_t i = 0; i != N; ++i) {
 
    if (in[i].hasNext()) {
 
      buffer.push_back({i, *in[i].next()});
 
      --number_of_events;
    }
  }
 
  make_heap(buffer.begin(), buffer.end());
 
  while (number_of_events != 0) {                           // process
 
    STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
 
    const size_t index = buffer[0].index;
 
    outputFile.put(buffer[0].evt);
      
    pop_heap(buffer.begin(), buffer.end());
 
    buffer.pop_back();
 
    if (in[index].hasNext()) {
 
      buffer.push_back({index, *in[index].next()});
 
      push_heap(buffer.begin(), buffer.end());      
 
      --number_of_events;
    }
  }
 
  sort_heap(buffer.begin(), buffer.end());                  // flush
 
  for (vector<entry_type>::const_reverse_iterator top = buffer.rbegin(); top != buffer.rend(); ++top) {
    outputFile.put(top->evt);
  }
#endif
  STATUS(endl);
 
  int counter = 0;
 
  Head      h0;
  MultiHead m0;
 
  for (size_t i = 0; i != N; ++i) {
 
    try {
 
      const Head      h1 = getHeader     (inputFile[i]);
      const MultiHead m1 = getMultiHeader(inputFile[i]);
 
      if (counter == 0) {
 
        outputFile.put(h1);
        outputFile.put(m1);
 
      } else {
 
        if (h0 != h1) { FATAL("Inconsistent header "       << endl << h1 << endl << "!=" << endl << h0 << endl); }
        if (m0 != m1) { FATAL("Inconsistent multi header " << endl << m1 << endl << "!=" << endl << m0 << endl); }
      }
 
      counter += 1;
 
      h0 = h1;
      m0 = m1;
    }
    catch (const exception&) {}
  }
 
  if (counter == 0) {
    FATAL("Missing header." << endl);
  }
 
  for (size_t i = 0; i != N; ++i) {
 
    JSingleFileScanner<JMeta> io(inputFile[i]);
 
    io >> outputFile;
  }
 
  outputFile.close();
 
  return 0;
}