JAcousticsMonitorRateTest.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 output_summary;
  string output_status;
  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['i'] = make_field(inputFile,          "output root file from JAcousticsMonitorRateSummary");
    zap['@'] = make_field(parametersFile,     "ASCII formatted input file with test criteria (acoustics_monitor)");
    zap['F'] = make_field(facet,              "Color facet")                              = get_keys(color_facets);
    zap['d'] = make_field(debug)              = 1;
    zap['s'] = make_field(output_summary,     "output summary file");
    zap['o'] = make_field(output_status,      "output information file working/not working");
    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_summary(output_summary.c_str());
  out_summary.imbue(locale(out_summary.getloc(), color_facets[facet]->clone()));
  out_summary << "ACOUSTIC MONITORING \nRun: " << run << endl;
  out_summary << "\n(Note: red highlights are the reason for the warning)" << endl;
  
  ofstream out_status;
  out_status.open(output_status.c_str());
  out_status << "INFORMATION ABOUT RECEIVER STATUS \nRun: " << run << 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_summary << RED << "No acoustic data." << endl; 
    out_summary << 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_summary << "\nEmitter: " << i->first << endl;
    out_status << "\nEmitter: " << i->first << endl;
    const TRegexp regexp("H\\[" + i->first + "\\]");
 
    TH2D* h3 = H3[stoi(i->first)];
 
    // if no data at all, all receivers out of range
    if (no_data > 0) {
      out_summary << "No data at all." << endl;
      out_status << "expected rate: " << i->second.expected_rate << endl;
      out_status << "min rate fraction: " << i->second.range.getLowerLimit() << endl; 
      out_status << "max rate fraction: " << i->second.range.getUpperLimit() << endl;
      out_status << "max allowed outliers: " << i->second.number_of_outliers << endl;
      out_status << "\nDU\tFloor\tOutOfRange\tSupposedToWork\tRate" << endl;
      for (Int_t ix = 1; ix <= h3->GetXaxis()->GetNbins(); ++ix) {
        for (Int_t iy = 1; iy <= h3->GetYaxis()->GetNbins(); ++iy) {
          string du = to_string(h3->GetXaxis()->GetBinLabel(ix));
          string floor = to_string(h3->GetYaxis()->GetBinLabel(iy));
          h3->Fill(ix-1, iy-1, 1.0);
          out_status << du << "\t" << floor << "\t" << 1 << "\t" << i->second.working << "\t" << 0 << endl; // for all receivers OutOfRange=1, Rate=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_summary << (i->second.working != 0 ? RED : GREEN);
          out_summary << "Emitter started working." << endl;
          out_summary << RESET;
          out_summary << "(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) {
            out_summary << "Testing acoustic rate." << endl;
            out_status << "expected rate: " << i->second.expected_rate << endl;
            out_status << "min rate fraction: " << i->second.range.getLowerLimit() << endl; 
            out_status << "max rate fraction: " << i->second.range.getUpperLimit() << endl;
            out_status << "max allowed outliers: " << i->second.number_of_outliers << endl;
            out_status << "\nDU\tFloor\tOutOfRange\tSupposedToWork\tRate" << endl;
 
            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;
                int status = 0; // assume receiver is working
                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);
                  status = 1; // set receiver to not working
                
                  ++number_of_outliers;
                  
                  if (working) { // receiver expected to work
                
                    if (floor == "000" && z == 0) { // if hydrophone doesn't work
                      out_summary << "DU " << du << ", floor " << floor << " : Acoustic rate out of range -> " << z << " Hz." << RED << " Hydrophone stopped working." << endl;
                      out_summary << RESET;
                      ++number_of_failures;
                    } else {
                      out_summary << "DU " << du << ", floor " << floor << " : Acoustic rate out of range -> " << z << " Hz." << endl;
                      }
                  }
                } else if (!working) { // receiver expected to not work                                               
                  out_summary << "DU " << du << ", floor " << floor << " : Working again -> " << z << " Hz." << endl;
                }
              out_status << du << "\t" << floor << "\t" << status << "\t" << working << "\t" << z <<  endl;
              }
            }
          } else { FATAL("Object at " << objectName << " is not TH2." << endl); }
            
          out_summary << (number_of_outliers > outliers ? RED : GREEN) << "Number of outliers = " << number_of_outliers << "/" << number_of_bins << endl;
          out_summary << RESET;
            
          if (number_of_outliers > outliers) {
          
            ++number_of_failures;
            out_summary << (number_of_outliers > outliers ? RED : GREEN) << "Test failed." << endl;
            out_summary << RESET;
          } else {
        
          out_summary << (outliers >= number_of_outliers ? GREEN : RED) << "Test passed." << endl;
          out_summary << RESET;
          }
        }
        break;
 
      // data from emitter is not in input file
      } else if (objectID + 1 == objectIDs.cend()) {
        out_summary << "Data from emitter not in input file." << endl;
        out_status << "expected rate: " << i->second.expected_rate << endl;
        out_status << "min rate fraction: " << i->second.range.getLowerLimit() << endl; 
        out_status << "max rate fraction: " << i->second.range.getUpperLimit() << endl;
        out_status << "max allowed outliers: " << i->second.number_of_outliers << endl;
        out_status << "\nDU\tFloor\tOutOfRange\tSupposedToWork\tRate" << endl;
        // 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) {
            string du = to_string(h3->GetXaxis()->GetBinLabel(ix));
            string floor = to_string(h3->GetYaxis()->GetBinLabel(iy));
            h3->Fill(ix-1, iy-1, 1.0);
            out_status << du << "\t" << floor << "\t" << 1 << "\t" << i->second.working << "\t" << 0 << endl; // for all receivers OutOfRange=1, Rate=0
          }
        }
        
        // check if expected status emitter changed
        if (i->second.working) {
          ++number_of_failures;
          out_summary << (i->second.working != 0 ? RED : GREEN) << "Emitter stopped working." << endl;
          out_summary << RESET;
        } else {
          out_summary << "Emitter is expected to not work." << endl;
        }
      }
    }
  }
  out_summary.close();
  out_status.close();
 
  return 0;
}