JAcousticsDisable.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <set>
#include <limits>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TKey.h"
#include "TString.h"
#include "TRegexp.h"

#include "JAcoustics/JTransmission_t.hh"

#include "JSupport/JMeta.hh"

#include "Jeep/JContainer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

/**
 * \author mdejong
 * Auxiliary program to set disable status of transmission based on time-of-arrival histograms.
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JContainer< set<JTransmission_t> > disable_container;

  vector<string>          inputFile;
  int                     Nmin;
  double                  Tmax_us;
  string                  disableFile;
  int                     debug;

  try { 

    JParser<> zap("Auxiliary program to set disable status of transmission based on time-of-arrival histograms.");
    
    zap['f'] = make_field(inputFile,       "input file (output from JCanberra).");
    zap['N'] = make_field(Nmin,            "minimum number of entries")      = 1;
    zap['T'] = make_field(Tmax_us,         "maximal time [us]")              = 60.0;
    zap['!'] = make_field(disableFile,     "disable transmission file")      = "";
    zap['d'] = make_field(debug,           "debug.")                         = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  disable_container disable;

  if (disableFile != "") {
    
    ifstream in(disableFile.c_str());

    in >> disable;

    in.close();
  }

  disable.comment.add(JMeta(argc, argv));


  const vector<double> Q = { 0.333, 0.666 };

  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    DEBUG("Processing " << *i << endl) ;

    TFile in(i->c_str(), "read");
    
    TIter iter(in.GetListOfKeys());

    for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {

      if (TString(key->GetName()).EndsWith(".toa")) {

        TH1* h1 = dynamic_cast<TH1*>(key->ReadObj());

        if (h1 != NULL) {

          TString buffer(h1->GetName());

          const char* regexp = "[^0-9][0-9]* [0-9]*[^0-9]";

          buffer = buffer(TRegexp(regexp));
          buffer = buffer(1, buffer.Length() - 2);

          if (buffer.Length() > 0 && buffer.IsDigit()) {

            JTransmission_t id;

            istringstream(buffer.Data()) >> id;

            int    N    = h1->GetEntries();
            double T_us = numeric_limits<double>::max();

            if (N > 0) {

              vector<double> R(Q.size());

              h1->GetQuantiles(Q.size(), R.data(), Q.data());

              T_us = (*R.rbegin() - *R.begin()) * 1.0e6;
            }

            if (N < Nmin || T_us > Tmax_us) {

              NOTICE("disable: " << id << ' ' 
                     << setw(6)    << N     << (N    < Nmin    ? "*" : "") << "  " 
                     << FIXED(5,1) << T_us  << (T_us > Tmax_us ? "*" : "") << " [us]" << endl);

              disable.insert(id);
            }

          } else {

            ERROR("Histogram name " << h1->GetName() << " not compatible with regular expression " << regexp << "." << endl);
          }
        }
      }
    }
    
    in.Close();
  }

  if (disableFile != "") {
    
    ofstream out(disableFile.c_str());

    out << disable;

    out.close();
  }
}