software/JCompass/JAHRSMonitor.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"

#include "JLang/JComparator.hh"
#include "JLang/JComparison.hh"

#include "JDB/JSupport.hh"
#include "JDB/JAHRS.hh"
#include "JDB/JAHRSCalibration_t.hh"
#include "JDB/JAHRSToolkit.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSystem/JStat.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JRootToolkit.hh"

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

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


/**
 * \file
 *
 * Program to monitor AHRS data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<JAHRS>  inputFile;
  JLimit_t&       numberOfEvents = inputFile.getLimit();
  string          outputFile;
  string          detectorFile;
  string          ahrsFile;
  double          precision;
  vector<double>  Q;
  int             qaqc;
  int             debug;

  try {

    JParser<> zap("Program to monitor AHRS data.");
    
    zap['f'] = make_field(inputFile,    "output of JConvertDB -q ahrs");
    zap['n'] = make_field(numberOfEvents)              = JLimit::max();
    zap['o'] = make_field(outputFile,   "");
    zap['a'] = make_field(detectorFile);
    zap['c'] = make_field(ahrsFile,     "output of JAHRSCalibration");
    zap['p'] = make_field(precision,    "precision")   = 1.0e-6;
    zap['q'] = make_field(Q,            "quantile");
    zap['Q'] = make_field(qaqc)                        = 0;
    zap['d'] = make_field(debug)                       = 2;

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


  JDetector detector;

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

  const JModuleRouter router(detector);

  const JAHRSCalibration_t calibration(ahrsFile.c_str());
  const JAHRSValidity      is_valid;

  const floor_range   range = getRangeOfFloors(detector);

  const JStringRouter string(detector);

  TH1D  h0("h0", NULL, 100, 0.0, 5.0);
  TH1D  h1("h1", NULL, 100, 0.0, 5.0);
  TH2D  h2("h2", NULL, 100, 0.0, 5.0, 100, 0.0, 5.0);
  TH1D  hn("hn", NULL, 100, 0.0, 1.0e4);
  
  TH2D* ha = new TH2D("ha", NULL, 
                      string.size(), -0.5, string.size() - 0.5,
                      range.getUpperLimit(), 1 - 0.5, range.getUpperLimit() + 0.5);

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

  TH2D* hb = (TH2D*) ha->Clone("hb");
  TH2D* hc = (TH2D*) ha->Clone("hc");
  TH2D* hd = (TH2D*) ha->Clone("hd");

  typedef map<int, int>  map_type;

  map_type data;

  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
    if (module->getFloor() != 0) {
      data[module->getID()] = 0;
    }
  }

  map<int, vector<JAHRS> > buffer;

  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    const JAHRS* parameters = inputFile.next();

    if (router.hasModule(parameters->DOMID)) {

      const JLocation& location = router.getModule(parameters->DOMID);

      if (location.getFloor() != 0) {

        const double A = sqrt(parameters->AHRS_A0 * parameters->AHRS_A0  +
                              parameters->AHRS_A1 * parameters->AHRS_A1  +
                              parameters->AHRS_A2 * parameters->AHRS_A2);
        const double H = sqrt(parameters->AHRS_H0 * parameters->AHRS_H0  +
                              parameters->AHRS_H1 * parameters->AHRS_H1  +
                              parameters->AHRS_H2 * parameters->AHRS_H2);
        
        h0.Fill(A);
        h1.Fill(H);
        h2.Fill(A,H);
        
        ha->Fill((double) string.getIndex(location.getString()), (double) location.getFloor());

        if (is_valid(*parameters)) {

          hb->Fill((double) string.getIndex(location.getString()), (double) location.getFloor());

          if (calibration.has(parameters->DOMID)) {
            hc->Fill((double) string.getIndex(location.getString()), (double) location.getFloor());
          }
        }

        if (is_valid(*parameters) && calibration.has(parameters->DOMID))
          buffer[parameters->DOMID].push_back(*parameters);
        else
          DEBUG("AHRS "
                << setw(16)   << parameters->UNIXTIME << ' '
                << setw(4)    << parameters->DUID     << ' '
                << setw(2)    << parameters->FLOORID  << ' '
                << FIXED(7,3) << parameters->AHRS_A0  << ' '
                << FIXED(7,3) << parameters->AHRS_A1  << ' '
                << FIXED(7,3) << parameters->AHRS_A2  << ' '
                << FIXED(7,3) << parameters->AHRS_H0  << ' '
                << FIXED(7,3) << parameters->AHRS_H1  << ' '
                << FIXED(7,3) << parameters->AHRS_H2  << ' '
                << is_valid(*parameters)
                << calibration.has(parameters->DOMID) << endl);
      }
    }
  }
  STATUS(endl);


  for (map<int, vector<JAHRS> >::iterator i = buffer.begin(); i != buffer.end(); ++i) {

    data[i->first] = 1;
    
    if (i->second.size() >= 2u) {

      sort(i->second.begin(), i->second.end(), make_comparator(&JAHRS::UNIXTIME));

      for (vector<JAHRS>::const_iterator q = i->second.begin(), p = q++; q != i->second.end(); ++p, ++q) {

        if ((fabs(p->AHRS_A0 - q->AHRS_A0) > precision ||
             fabs(p->AHRS_A1 - q->AHRS_A1) > precision ||
             fabs(p->AHRS_A2 - q->AHRS_A2) > precision)
            &&
            (fabs(p->AHRS_H0 - q->AHRS_H0) > precision ||
             fabs(p->AHRS_H1 - q->AHRS_H1) > precision ||
             fabs(p->AHRS_H2 - q->AHRS_H2) > precision)) {

          data[i->first] += 1;
        }
      }
    }
  }

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

    const JLocation& location = router.getModule(i->first);
    
    hn.Fill(i->second);

    hd->Fill((double) string.getIndex(location.getString()), (double) location.getFloor(), i->second);
  }

  if (debug >= status_t) {

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

      map_type::const_iterator p = data.find(module->getID());

      if (p != data.end()) {
        cout << getLabel(module->getLocation()) << ' ' << setw(8) << module->getID() << ' ' << setw(6) << p->second << endl;
      }
    }
  }
  
  if (outputFile != "") {

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

    out << h0  << h1  << h2 << hn;
    out << *ha << *hb << *hc << *hd;

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

  JQuantile q1("", true);

  for (map_type::const_iterator i = data.begin(); i != data.end(); ++i) {
    q1.put(i->second);
  }

  for (vector<double>::const_iterator i = Q.begin(); i != Q.end(); ++i) {
    QAQC(' ' << FIXED(8,3) << q1.getQuantile(*i));
  }
  QAQC(endl);

  return 0;
}