JAcousticsDisable.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <set>
#include <array>
#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 "JSystem/JGlob.hh"
#include "JSystem/JStat.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;
  typedef array<double, 2>                   array_type;
  
  vector<string>          inputFile;
  string                  outputFile;
  int                     Nmin;
  double                  Tmax_us;
  array_type              Q = { 0.1, 0.9 };
  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['o'] = make_field(outputFile)                                        = "disable.root";
    zap['N'] = make_field(Nmin,            "minimum number of entries")      = 1;
    zap['T'] = make_field(Tmax_us,         "maximal time [us]");
    zap['Q'] = make_field(Q,               "quantiles");
    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);
  }


  inputFile = getFilenames(inputFile);

  disable_container disable;

  if (disableFile != "" && getFileStatus(disableFile.c_str())) {
    disable.load(disableFile.c_str());
  }

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

  TFile out(outputFile.c_str(), "recreate");

  TH1D ha("ha", NULL, 1000, 0.0, 1.0e4);
  TH1D hb("hb", NULL, 1000, 0.0, 1.0e3);

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

    NOTICE("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;

            const int N = h1->GetEntries();

            if (N > 0 && N < Nmin) {

              NOTICE("disable: " << setw(2) << id.tx << ' ' << setw(8) << id.rx << ' ' << setw(6) << N << " < " << setw(6) << Nmin << endl);

              disable.insert(id);
            }

            ha.Fill(N);
            
            if (N > 0) {

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

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

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

              hb.Fill(T_us);

              DEBUG("transmission: " << setw(2) << id.tx << ' ' << setw(8) << id.rx << ' ' << setw(6) << N << ' ' << FIXED(6,1) << T_us << " [us]" << endl);
              
              if (T_us > Tmax_us) {

                NOTICE("disable: " << setw(2) << id.tx << ' ' << setw(8) << id.rx << ' ' << FIXED(6,1) << T_us << " > " << FIXED(6,1) << Tmax_us << " [us]" << endl);

                disable.insert(id);
              }
            }

          } else {

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

  if (disableFile != "") {
    disable.store(disableFile.c_str());
  }

  out.Write();
  out.Close();
}