JDST.cc

Go to the documentation of this file.

#include <cstring>
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
 
#include "TError.h"
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TTree.h"
#include "TChain.h"
 
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/root.hh"
 
#include "JSystem/JGlob.hh"
#include "JSupport/JLimit.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  /**
   * Auxiliary data structure to tag event.
   */
  struct JTag {
    /**
     * Defaul constructor.
     */
    JTag() :
      run(-1),
      frame_index(-1),
      counter(-1)
    {}
 
    /**
     * Constructor.
     *
     * \param  run           run
     * \param  frame_index   frame index
     * \param  counter       trigger counter
     */
    JTag(const int run,
         const int frame_index,
         const int counter) :
      run(run),
      frame_index(frame_index),
      counter(counter)
    {}
 
    /**
     * Read tag from input stream.
     *
     * \param  in            input stream
     * \param  tag           tag
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JTag& tag)
    {
      return in >> tag.run
                >> tag.frame_index
                >> tag.counter;
    }
 
    /**
     * Write tag to output stream.
     *
     * \param  out           output stream
     * \param  tag           tag
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JTag& tag)
    {
      using namespace std;
 
      return out << setw(8) << tag.run         << ' '
                 << setw(6) << tag.frame_index << ' '
                 << setw(6) << tag.counter;
    }
 
    /**
     * Less-than operator for tag.
     *
     * \param  first         first  tag
     * \param  second        second tag
     * \return               true of first tag less than second; else false
     */
    friend inline bool operator<(const JTag& first, const JTag& second)
    {
      if (first.run == second.run) {
 
        if (first.frame_index == second.frame_index)
          return first.counter     < second.counter;
        else
          return first.frame_index < second.frame_index;
 
      } else {
 
        return first.run < second.run;
      }
    }
 
    int run;               //!< run number
    int frame_index;       //!< frame index
    int counter;           //!< trigger counter
  };
}
 
 
/**
 * \file
 * Example program to select events from DST files.
 * 
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  vector<string> inputFile;
  JLimit_t       numberOfEvents;
  string         selectFile;
  string         outputFile;
  int            debug;
 
  try { 
 
    JParser<> zap("Example program to select events from DST files.");
    
    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)       = JLimit::max();
    zap['<'] = make_field(selectFile,
                          "<run>  <frame index>  <trigger counter>"\
                          " (one header line)") = "";
    zap['o'] = make_field(outputFile);
    zap['d'] = make_field(debug)                = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  gErrorIgnoreLevel = kFatal;
 
  inputFile = getFilenames(inputFile);
 
  set<JTag> selection;
 
  if (selectFile != "") {
 
    ifstream in(selectFile.c_str());
 
    in.ignore(numeric_limits<streamsize>::max(), '\n');
    
    for (JTag tag; in >> tag; ) {
      selection.insert(tag);
    }
 
    in.close();
  }
 
  TFile* out = TFile::Open(outputFile.c_str(), "CREATE");
  
  if (out == NULL || ! out->IsOpen()) {
    FATAL("File " << outputFile << " not opened." << endl);
  }
 
  vector<TChain*> input;
  vector<TTree*>  output;
 
  for (vector<string>::const_iterator file_name = inputFile.begin(); file_name != inputFile.end(); ++file_name) {
    
    TFile* in = TFile::Open(file_name->c_str(), "EXISTS");
 
    if (in != NULL && in->IsOpen()) {
    
      TIter iter(in->GetListOfKeys(), kIterBackward);
 
      for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
 
        TKey* p = dynamic_cast<TKey*>(in->GetListOfKeys()->Before(key));
 
        if (p == NULL || strcmp(key->GetName(), p->GetName()) != 0) {                     // select last key
 
          TTree* t1 = dynamic_cast<TTree*>(key->ReadObj());
 
          if (t1 != NULL) {
 
            TChain* q = NULL;
            
            for (vector<TChain*>::iterator i = input.begin(); i != input.end(); ++i) {
 
              if (strcmp((*i)->GetName(), t1->GetName()) == 0) {
 
                q = *i;
 
                break;
              }
            }
 
            if (q == NULL) {
              input.push_back(q = new TChain(t1->GetName()));
            }
 
            q->Add(file_name->c_str());
          }
        }
      }
    }
  }
 
  // master
 
  TChain* tm = NULL;
 
  for (vector<TChain*>::iterator i = input.begin(); i != input.end(); ++i) {
 
    if (strcmp((*i)->GetName(), TTREE_OFFLINE_EVENT) == 0) {
 
      tm = *i;
 
      break;
    }
  }
 
  if (tm == NULL) {
    FATAL("Missing TTree " << TTREE_OFFLINE_EVENT << endl);
  }
    
  Evt* evt = NULL;
 
  tm->GetEntries();                                                  // load TTree
  tm->SetBranchAddress(TBRANCH_OFFLINE_EVENT, &evt);                 // 
      
  // select friend TTree's based on number of entries
      
  for (vector<TChain*>::iterator i = input.begin(); i != input.end(); ) {
    if (tm->GetEntries() == (*i)->GetEntries())
      ++i;
    else
      i = input.erase(i);
  }
      
  out->cd();
 
  for (vector<TChain*>::iterator i = input.begin(); i != input.end(); ++i) {
    output.push_back((*i)->GetTree()->CloneTree(0));
  }
 
  for (Long64_t i = numberOfEvents.getLowerLimit(); i < tm->GetEntries() && i < numberOfEvents.getUpperLimit(); ++i) {
 
    STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);
 
    tm->GetEntry(i);
 
    const JTag tag(evt->run_id, evt->frame_index, evt->trigger_counter);
 
    DEBUG("tag " << tag << ' ' << (selection.empty() || selection.count(tag)) << endl);
    
    if (selection.empty() || selection.count(tag)) {
        
      for (size_t m = 0; m != input.size(); ++m) {
        input [m]->GetEntry(i);
        output[m]->Fill();
      }
    }
  }
  STATUS(endl);
 
  out->Write();
  out->Close();
}