JAcousticsMonitorTest.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <map>
#include <locale>

#include "TROOT.h"
#include "TFile.h"
#include "TObject.h"
#include "TKey.h"
#include "TString.h"
#include "TRegexp.h"
#include "TGraph.h"
#include "TF1.h"
#include "TH2D.h"

#include "JGizmo/JRootObjectID.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JTools/JRange.hh"
#include "JTools/JHashCollection.hh"
#include "JTools/JHashMap.hh"
#include "JLang/JLangToolkit.hh"
#include "JLang/JColorFacet.hh"
#include "JLang/JVectorize.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JColor.hh"

#include "JROOT/JManager.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

namespace {

  using JTOOLS::JRange;

  /**
   * Auxilliary data structure with test criteria.
   */
  struct JParameters_t {

    static const char SKIPLINE  =  '#';           //!< skip line character

    /**
     * Default constructor.
     */
    JParameters_t() :
      working(1),
      expected_rate(0),
      range(JRange<double>::DEFAULT_RANGE()),
      number_of_outliers(0)
    {}


    /**
     * Read parameters from input stream.
     *
     * \param  input           input stream
     * \param  object          parameters
     * \return                 input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JParameters_t& object)
    {
      using namespace std;
      return in >> object.working >> object.expected_rate >> object.range >> object.number_of_outliers;
    }


    /**
     * Write parameters to output stream.
     *
     * \param  output          output stream
     * \param  object          parameters
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JParameters_t& object)
    {
      using namespace std;
      using namespace JPP;

      return out << setw(6)     << object.working     << ' '
                 << FIXED(2,6)  << object.expected_rate    << ' '
                 << FIXED(2,6)  << object.range.getLowerLimit() << ' '
                 << FIXED(2,6)  << object.range.getUpperLimit() << ' '
                 << setw(4)     << object.number_of_outliers;
    }

    int working;
    double  expected_rate;
    JRange<double>  range;
    int  number_of_outliers;
  };

  /*                                                                                                                          
   * Gets list of keys in a ROOT TDirectory and stores it on a vector.                                                        
   *                                                                                                                          
   * \param  dir        The ROOT directory
   * \param  buffer     Vector to store keys 
   */
  void readDir(TDirectory*                         dir,
               std::vector<JGIZMO::JRootObjectID>& buffer) {

    TIter iter(dir->GetListOfKeys());

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

      if (key->IsFolder()){

        dir->cd(key->GetName());

        TDirectory *subdir = gDirectory;
        readDir(subdir, buffer);

        dir->cd();

      } else {

        JGIZMO::JRootObjectID objectID(MAKE_STRING(dir->GetPath() << key->GetName()));

        buffer.push_back(objectID);
      }
    }
  }

}


/**
 *
 * Program to test:
 * - If there is acoustic data.
 * - If the rate recorded by each receiver is within the expected range.
 * - If emitters and receivers are working and compares it to what is expected.
 *
 * Reports failure if: no acoustic data, too large number of outliers, emitter expected to work is not working,
 * hydrophone expected to work is not working.
 *
 */

int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string inputFile;
  string parametersFile;
  string facet;
  string outputFile1;
  string outputFile2;
  string detectorFile;
  int    debug;
  int    number_of_failures = 0;
  int    run;
 
  try {

    JParser<> zap("Auxiliary program to apply test criteria to 2D histograms monitoring acoustic rate per emitter.");
    
    zap['f'] = make_field(inputFile,          "output root file from JAcousticsMonitor_short");
    zap['P'] = make_field(parametersFile,     "ASCII formatted input file with test criteria (acoustic_monitor_00000XXX.txt)");
    zap['F'] = make_field(facet,              "Color facet")                                       = get_keys(color_facets);
    zap['d'] = make_field(debug)              = 1;
    zap['w'] = make_field(outputFile1,         "output summary file");
    zap['t'] = make_field(outputFile2,         "output root file");
    zap['a'] = make_field(detectorFile);
    zap['r'] = make_field(run,                "run number");
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  ofstream out(outputFile1.c_str());
  out.imbue(locale(out.getloc(), color_facets[facet]->clone()));
  out << "ACOUSTIC MONITORING \nRun: " << run << endl;
  out << "\n(Note: red highlights are the reason for the warning)" << endl;

  // read parameters file

  typedef map<string, JParameters_t> map_type;

  map_type zmap;

  ifstream in(parametersFile.c_str());

  if (in) {

    string        key;
    JParameters_t parameters;

    for (string buffer; getline(in, buffer); ) {

      if (!buffer.empty() && buffer[0] != JParameters_t::SKIPLINE) {

        istringstream is(buffer);

        if (is >> key >> parameters) { zmap[key] = parameters; }
      }
    }

    in.close();

  } else {
    FATAL("Error opening file: " << parametersFile << endl);
  }

  // create hist for test output

  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  JHashMap<int, JLocation>   receivers;

  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {

    receivers[i->getID()] = i->getLocation();
  }

  const JHashCollection<int> string1(make_array(detector.begin(), detector.end(), &JModule::getString));
  const JRange         <int> floor1 (make_array(detector.begin(), detector.end(), &JModule::getFloor));

  JManager<int, TH2D> H3(new TH2D("H[%].rate-test",        NULL,
                                  string1.size(),             - 0.5, string1.size() - 0.5,
                                  floor1.getUpperLimit() + 1, - 0.5, floor1.getUpperLimit() + 0.5));

  for (Int_t i = 1; i <= H3->GetXaxis()->GetNbins(); ++i) {
    H3->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string1.at(i-1)));
  }

  for (Int_t i = 1; i <= H3->GetYaxis()->GetNbins(); ++i) {
    H3->GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i-1));
  }

  // read input file

  TFile* f = TFile::Open(inputFile.c_str());  

  vector<JRootObjectID> objectIDs;

  readDir(f,objectIDs);

  // check if the monitor.root file is empty  
  int no_data = 0;
  if(objectIDs.empty()) { 
    ++number_of_failures;
    no_data += 1;
    out << RED << "No acoustic data." << endl; 
    out << RESET;    
  }

  // loop over expected emitters

  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {

    if (i->first[0] == '0') { continue; } // skip receivers parameters file

    out << "\nEmitter: " << i->first << endl;

    const TRegexp regexp(i->first);

    TH2D* h3 = H3[stoi(i->first)];

    // if no data at all, all receivers out of range
    if (no_data > 0) {
      for (Int_t ix = 1; ix <= h3->GetXaxis()->GetNbins(); ++ix) {
        for (Int_t iy = 1; iy <= h3->GetYaxis()->GetNbins(); ++iy) {
          h3->Fill(ix-1, iy-1, 1.0);
        }
      }
    }

    // check if data from emitter is in input file    
    for (vector<JRootObjectID>::const_iterator objectID = objectIDs.cbegin() ; objectID != objectIDs.cend() ; ++objectID) {
    
      const TString& objectName = objectID->getFullObjectName();
   
      // data from emitter is in input file
      if (objectName.Index(regexp) != -1) {

        // check if expected status emitter changed
        if (!i->second.working) {
          out << (i->second.working != 0 ? RED : GREEN);
          out << "Emitter started working." << endl;
          out << RESET;   
          out << "(Acoustic rates not tested)" << endl;
        }

        // if emitter expected to work, test acoustic rate
        if (i->second.working) {
          TObject* p = (TObject*)f->Get(objectName);
          TH2*     h2 = NULL;

          int number_of_bins   = 0;
          int number_of_outliers = 0;

          double min_rate = i->second.expected_rate * i->second.range.getLowerLimit();
          double max_rate = i->second.expected_rate * i->second.range.getUpperLimit();
          int outliers = i->second.number_of_outliers;
        
          if (h2 == NULL && dynamic_cast<TH2*>(p) != NULL) { h2 = dynamic_cast<TH2*>(p); };

          if (h2 != NULL) {

            for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix) {
              for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {

                const Double_t z = h2->GetBinContent(ix, iy);
                string du = to_string(h2->GetXaxis()->GetBinLabel(ix));
                string floor = to_string(h2->GetYaxis()->GetBinLabel(iy));
                
                du.insert(du.begin(), 4 - du.length(), '0');
                floor.insert(floor.begin(), 3 - floor.length(), '0');
                const string id = du + '.' + floor;

                int working = 1;
                if (zmap.find(id)!=zmap.end()) { working = zmap.find(id)->second.working; }

                ++number_of_bins;

                // acoustic rate out of range
                if (z < min_rate || z > max_rate) {

                  h3->Fill(ix-1, iy-1, 1.0);

                  ++number_of_outliers;
                  if (working) { // receiver expected to work

                    if (floor == "000" && z == 0) { // if hydrophone doesn't work
                      out << "DU " << du << ", floor " << floor << " : Acoustic rate out of range -> " << z << " Hz." << RED << " Hydrophone stopped working." << endl;
                      out << RESET;
                      ++number_of_failures;
                    } else {
                      out << "DU " << du << ", floor " << floor << " : Acoustic rate out of range -> " << z << " Hz." << endl;
                    }
                  }
                } else if (!working) { // receiver expected to not work                                               
                  out << "DU " << du << ", floor " << floor << " : Working again -> " << z << " Hz." << endl;
                }
              }
            }
          
          } else { FATAL("Object at " << objectName << " is not TH2." << endl); }

          out << (number_of_outliers > outliers ? RED : GREEN) << "Number of outliers = " << number_of_outliers << "/" << number_of_bins << endl;
          out << RESET;

          if (number_of_outliers > outliers) {
          
            ++number_of_failures;
            out << (number_of_outliers > outliers ? RED : GREEN) << "Test failed." << endl;
            out << RESET;
          } else {
          
            out << (outliers >= number_of_outliers ? GREEN : RED) << "Test passed." << endl;
            out << RESET;
          }
        }
        break;
        
  
      // data from emitter is not in input file
      } else if (objectID + 1 == objectIDs.cend()) {

        // for all receivers rate out of range                                                                                        
        for (Int_t ix = 1; ix <= h3->GetXaxis()->GetNbins(); ++ix) {
          for (Int_t iy = 1; iy <= h3->GetYaxis()->GetNbins(); ++iy) {
            h3->Fill(ix-1, iy-1, 1.0);
          }
        }

        // check if expected status emitter changed
        if (i->second.working) {
          ++number_of_failures;
          out << (i->second.working != 0 ? RED : GREEN) << "Emitter stopped working." << endl;
          out << RESET;
        } else {
          out << "Emitter is expected to not work." << endl;
        }
      }
    }
  }
  out.close();

  TFile out2(outputFile2.c_str(), "recreate");

  out2 << H3;

  out2.Write();
  out2.Close();

  // if number_of_failures > 0 program output needs to be reported
  if (number_of_failures > 0) {
    number_of_failures = 1;
    cout << number_of_failures << endl;
  }
  return 0;
}