extractEvent.cc

Go to the documentation of this file.

/*
  extract events from IoAuger file via
   - SD ID
   - Auger ID
   - Events in ADST file (using Auger IDs)

  Authors:
    Michael.Unger@kit.edu
    Fabian.Schuessler@ik.fzk.de
    Ralf.Ulrich@kit.edu
    Hans.Dembinski@kit.edu
    darko.veberic@kit.edu
*/

#include <RecEvent.h>
#include <RecEventFile.h>

#include <IoAuger.h>
#include <IoSdData.h>
#include <AugerEvent.h>

#include <boost/algorithm/string/trim.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/program_options.hpp>
#include <boost/filesystem/operations.hpp>

#include <iostream>
#include <sstream>
#include <string>
#include <fstream>
#include <algorithm>

using namespace std;


typedef unsigned long ulong;


namespace Global {
  bool gDebug = false;
}


set<ulong>
ParseIds(const string& file)
{
  ifstream ifs(file);
  if (!ifs.is_open()) {
    cerr << "cannot open file '" << file << "'!" << endl;
    exit(EXIT_FAILURE);
  }
  string line;
  ostringstream oss;
  int lineno = 0;
  set<ulong> ids;
  while (getline(ifs, line)) {
    ++lineno;
    if (Global::gDebug)
      cerr << "line " << lineno << ": \"" << line << "\"\n";
    try {
      boost::algorithm::trim(line);
      const auto id = boost::lexical_cast<ulong>(line);
      ids.insert(id);
      if (Global::gDebug)
        cerr << "got " << id << '\n';
    } catch (const boost::bad_lexical_cast& e) {
      cerr << "line " << lineno << " of the file '" << file << "' could not be converted to an ID" << endl;
      throw e;
    }
  }
  return ids;
}


void
DumpIds(const set<ulong>& ids)
{
  for (const auto& id : ids)
    cerr << id << '\n';
}


void
DumpKeys(const vector<string>& files)
{
  for (const auto& file : files) {
    const unique_ptr<TFile> f{TFile::Open(file.c_str(), "READ")};
    if (!f) {
      cerr << "cannot read file '" << file << "'" << endl;
      continue;
    }
    const auto* const keys = f->GetListOfKeys();
    const auto nKeys = keys->GetEntries();
    if (!nKeys) {
      cerr << "warning: file '" << file << "' has no keys!" << endl;
      continue;
    }
    for (int i = 0; i < nKeys; ++i) {
      const auto key = (TKey* const)keys->At(i);
      if (!key) {
        if (Global::gDebug)
          cerr << "key " << i << " is a nullptr!";
        continue;
      }
      const char* const name = key->GetName();
      if (!name || !('0' <= *name && *name <= '9')) {
        if (Global::gDebug)
          cerr << "could not convert key \"" << name << '"' << endl;
        continue;
      }
      cout << stoul(name) << endl;
    }
  }
}


int
main(int argc, char* argv[])
{
  namespace bpo = boost::program_options;
  bpo::options_description optionSwitches("Options");
  optionSwitches.add_options()
    ("help,h",
     "Output this help.")
    ("debug",
     "Output some debug information.")
    ("ids,i",
     bpo::value<string>(),
     "File containing event IDs to select.")
    ("dump,d",
     "Dump all IDs in input IoAuger ROOT files.")
    ("native,n",
     "Search for events through native IoAuger interface (instead using faster ROOT keys).")
    ("output,o",
     bpo::value<string>(),
     "Output IoAuger ROOT file containing all selected events.")
    ;
  bpo::options_description optionHidden;
  optionHidden.add_options()
    ("input",
     bpo::value<vector<string>>()->required(),
     "Input ROOT file(s).");
  bpo::positional_options_description positional;
  positional.add("input", -1);
  bpo::options_description options;
  options.add(optionSwitches).add(optionHidden);
  bpo::variables_map vm;
  try {
    bpo::store(bpo::command_line_parser(argc, argv).options(options).positional(positional).run(), vm);
    if (vm.count("help")) {
      cerr << optionSwitches << endl;
      return EXIT_SUCCESS;
    }
  } catch (bpo::error& err) {
    cerr << "Command line error : " << err.what() << '\n'
         << optionSwitches << endl;
    return EXIT_FAILURE;
  }

  Global::gDebug = vm.count("debug");

  const bool native = vm.count("native");

  if (!vm.count("input")) {
    cerr << "no input files specified!" << endl;
    return EXIT_FAILURE;
  }

  const auto inputFiles = vm["input"].as<vector<string>>();

  if (vm.count("dump")) {
    DumpKeys(inputFiles);
    return EXIT_SUCCESS;
  }

  if (!vm.count("output")) {
    cerr << "no output file specified!" << endl;
    return EXIT_FAILURE;
  }

  const auto outputFile = vm["output"].as<string>();

  if (boost::filesystem::exists(outputFile)) {
    cerr << "not overwriting exisiting output file '" << outputFile << '\'' << endl;
    return EXIT_FAILURE;
  }

  if (!vm.count("ids")) {
    cerr << "no ID file specified!" << endl;
    return EXIT_FAILURE;
  }

  auto ids = ParseIds(vm["ids"].as<string>());

  const size_t nSearched = ids.size();

  if (!nSearched) {
    cerr << "no events to extract; IDs file is empty!" << endl;
    return EXIT_FAILURE;
  }

  if (Global::gDebug) {
    cerr << "IDs to find:\n";
    DumpIds(ids);
  }

  cout << "Events to extract: " << nSearched << "\n"
          "Using " << (native ? "native " : "key ") << "interface" << endl;

  AugerEvent event;

  unique_ptr<TFile> outputRootFile;
  unique_ptr<AugerFile> outputAugerFile;
  size_t nFound = 0;

  if (native) {

    // slow method: read each file and look for ID inside AugerEvent

    for (const auto& inFile : inputFiles) {
      cout << "Opening " << inFile << " ... [" << nFound << '/' << nSearched << ']' << endl;
      AugerFile auFile(inFile.c_str());
      while (auFile.ReadNext(event) != AugerFile::eEOF && nFound != nSearched) {
        auto it = ids.find(event.Id());
        if (it == ids.end() && event.HasSd())
          it = ids.find(event.Sd().Id);
        if (it != ids.end()) {
          ++nFound;
          if (!outputAugerFile)
            outputAugerFile = make_unique<AugerFile>(outputFile.c_str(), AugerFile::eWrite);
          outputAugerFile->Write(event);
          if (Global::gDebug)
            cout << "found ID " << *it << endl;
          ids.erase(it);
          if (ids.empty())
            break;
        }
      }
    }

  } else {

    // fast method: ID is encoded in the key name, just browse keys

    for (const auto& inFile : inputFiles) {

      cout << "Opening " << inFile << " ... [" << nFound << '/' << nSearched << ']' << endl;

      const unique_ptr<TFile> f{TFile::Open(inFile.c_str(), "READ")};
      if (!f) {
        cerr << "Error: cannot read file '" << inFile << '\'' << endl;
        continue;
      }

      const auto* const keys = f->GetListOfKeys();
      const auto nKeys = keys->GetEntries();
      if (!nKeys) {
        cerr << "warning: file '" << inFile << "' has no keys!" << endl;
        continue;
      }
      if (Global::gDebug)
        cerr << "file '" << inFile << "' has " << nKeys << " keys" << endl;
      for (int i = 0; i < nKeys; ++i) {
        const auto key = (TKey* const)keys->At(i);
        if (!key) {
          cerr << "key " << i << " is a nullptr!";
          continue;
        }
        const char* const name = key->GetName();
        if (!name || !('0' <= *name && *name <= '9')) {
          cerr << "key " << i << " is malformed: \"" << (name ? name : "nullptr") << '\"' << endl;
          continue;
        }
        const ulong id = stoul(name);
        if (Global::gDebug)
          cout << "id from key " << i << " \"" << name << "\" is " << id << endl;
        const auto it = ids.find(id);
        if (it != ids.end()) {
          ++nFound;
          if (Global::gDebug)
            cout << "found ID " << id << endl;
          if (!outputRootFile)
            outputRootFile.reset(TFile::Open(outputFile.c_str(), "RECREATE"));
          auto* const obj = key->ReadObj();
          outputRootFile->cd();
          obj->Write(name);
          ids.erase(it);
          if (ids.empty())
            break;
        }
      }

    }

  }

  if (Global::gDebug) {
    cerr << "IDs not found:\n";
    DumpIds(ids);
  }

  cout << "Events extracted: [" << nFound << "/" << nSearched << "]" << endl;

  return EXIT_SUCCESS;
}