JDataQuality.cc

Go to the documentation of this file.

#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "TString.h"
#include "TNtuple.h"
#include "TRegexp.h"

#include "JDB/JDB.hh"
#include "JDB/JSelector.hh"
#include "JDB/JSelectorSupportkit.hh"
#include "JDB/JRunQuality.hh"
#include "JDB/JRunsetups.hh"
#include "JDB/JDBToolkit.hh"
#include "JDB/JGITTags.hh"

#include "JLang/JLangToolkit.hh"
#include "JTools/JRange.hh"
#include "JROOT/JRootClass.hh"
#include "JROOT/JRootPrinter.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JRootStreamObjectOutput.hh"
#include "JROOT/JRootDictionary.hh"
#include "JROOT/JASCIIFileStreamer.hh"
#include "JROOT/JManager.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JFilenameSupportkit.hh"

#include "Jeep/JeepToolkit.hh"
#include "Jeep/JComment.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

namespace {

  /**
   * Auxiliary data structure for run quality evaluation.
   *
   * A condition is composed of a ROOT TFormula compatible expression and a validity range for the result.
   */
  struct JCondition {

    static const char SEPARATOR = ';';

    /**
     * Read selection from input stream.
     *
     * \param  in           input stream
     * \param  object       selection
     * \return              input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JCondition& object)
    {
      std::getline(in, object.formula, SEPARATOR);

      object.range  = JTOOLS::JRange<double>();
      object.weight = 1.0;

      if ((in >> object.range) && ! (in >> object.weight)) {
        in.clear();
      }

      return in;
    }


    /**
     * Write selection to output stream.
     *
     * \param  out          output stream
     * \param  object       selection
     * \return              output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JCondition& object)
    {
      out << object.formula << SEPARATOR << object.range << ' ' << object.weight << std::endl;

      return out;
    }

    std::string            formula;
    JTOOLS::JRange<double> range;
    double                 weight;
  };
}


/**
 * \file
 *
 * Main program to evaluate quality data from data base.
 *
 * The evaluation of the data quality is based on the following procedure.\n
 * For each data taking run, the script JQAQC.sh is executed which produces a set of values. 
 * The values are subsequently uploaded into the database in table "runsummarynumbers". 
 * This table is a general purpose table for any value related to some data taking run.
 * With this application,
 * the QA/QC data are consistently downloaded from the database;
 * evaluated according some formula and
 * tested against a validity range.
 * The underlying data structure is JDATABASE::JRunQuality.
 *
 * The formula and validity range are referred to as a condition and can be specified with option -Q.
 * Each condition has an optional weight (default is 1).
 * The normalised total weight is used to define an overall quality, ranging from 0 = worst to 1 = best.
 * There is also the possibility to specify a veto (option -V),
 * which is similar to a condition but only the number of vetoes is counted.
 * The overall quality and the number of vetoes can be used to accept or reject a data taking run.
 *
 * By default, a list of the run numbers, the overall weight and the number of vetoes is printed.
 * If a ROOT file is specified for the output (i.e.\ filename extension is JSUPPORT::ROOT_FILE_FORMAT),
 * a file is written which contains a histogram for each condition and each veto; some summary histograms and a ROOT n-tuple.
 * If an ASCII file is specified for the output (i.e.\ filename extension is JSUPPORT::ASCII_FILE_FORMAT),
 * a text file is written which can subsequently be used as input file.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JRange<int>           JRange_t; 
  typedef vector<JCondition>    JCondition_t;

  JServer         server;
  string          usr;
  string          pwd;
  string          cookie;
  string          inputFile;
  string          outputFile;
  string          detid;
  JRange_t        runs;
  vector<string>  source;
  JCondition_t    condition;
  JCondition_t    veto;
  TRegexp         regexp(".");
  JSelector       selection;
  int             debug;

  try {

    JParser<> zap("Main program to evaluate quality data from data base."\
                  "\nThe condition and veto can be any TFormula compatible expression involving QA/QC parameters (see e.g. JQAQC.sh -h).");
    
    zap['s'] = make_field(server)        = getServernames();
    zap['u'] = make_field(usr)           = "";
    zap['!'] = make_field(pwd)           = "";
    zap['C'] = make_field(cookie)        = "";
    zap['f'] = make_field(inputFile,  "Optional input file instead of database.")       = "";
    zap['o'] = make_field(outputFile, "ROOT file with histograms and n-tuple or ASCII file with QA/QC data.") = "quality.root";
    zap['D'] = make_field(detid)         = "";
    zap['R'] = make_field(runs,       "Run range")                                      = JRange_t(1, JRange_t::getMaximum());
    zap['S'] = make_field(source,     "GIT versions")                                   = getGITTags(TRegexp("v[0-9]*\\.[0-9]*\\.[0-9]*$"), JGITTags_t::key_type("2019-04-12"));
    zap['Q'] = make_field(condition,  "User defined conditions")                        = JPARSER::initialised();
    zap['V'] = make_field(veto,       "User defined vetos.")                            = JPARSER::initialised();
    zap['r'] = make_field(regexp,     "TRegexp for selection of run setup names.")      = JPARSER::initialised();
    zap['@'] = make_field(selection)     = JPARSER::initialised();
    zap['d'] = make_field(debug,      "Debug level")                                    = 1;

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


  double W = 0.0;
  
  for (JCondition_t::const_iterator i = condition.begin(); i != condition.end(); ++i) {
    W += i->weight;
  }

  if (getFilenameExtension(outputFile) == ROOT_FILE_FORMAT) {
    if (W <= 0.0) {
      FATAL("Invalid total weight: " << W << endl);
    }
  }


  JRunsetups       setups;                  // run -> value
  set<JRunQuality> buffer;                  // storage


  if (inputFile == "") {                    // read data from database

    ASSERT(detid != "");

    try {

      JDB::reset(usr, pwd, cookie);
      
      const int ID = getDetector<int>   (detid);
      detid        = getDetector<string>(detid);

      // run setups

      selection += getSelector<JRuns>(ID);

      NOTICE("Extracting run information from database... " << flush);

      ResultSet& rs  = getResultSet(getTable<JRuns>(), selection);

      for (JRuns parameters; rs >> parameters; ) {
        if (TString(parameters.RUNSETUPNAME.c_str()).Contains(regexp)) {
          setups.put(parameters);
        }
      }

      rs.Close();

      NOTICE("OK" << endl);

      // run summary data

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

        typedef map<string, string>    data_type;
        typedef map<int, data_type>    map_type;

        map_type zmap;

        JSelector selector = getSelector<JRunSummaryNumbers>(detid, runs.getLowerLimit(), runs.getUpperLimit());

        selector.add(&JRunSummaryNumbers::SOURCE_NAME, *git);

        try {

          NOTICE("Extracting run summmary information with source " << *git << " from database... " << flush);

          ResultSet& rs  = getResultSet(getTable<JRunSummaryNumbers>(), selector);

          for (JRunSummaryNumbers parameters; rs >> parameters; ) {
            if (setups.has(parameters.RUN)) {
              zmap[parameters.RUN].insert(make_pair(parameters.PARAMETER_NAME, parameters.DATA_VALUE));
            }
          }

          rs.Close();

          NOTICE("OK" << endl);
        }
        catch(const exception& error) { NOTICE(endl); }

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

          JRunQuality quality;
          
          quality.GIT      = *git;
          quality.detector = ID;
          quality.run      = run->first;
          quality.name     = setups[run->first].name;
          quality.value    = setups[run->first].value;

          for (data_type::const_iterator i = run->second.begin(); i != run->second.end(); ++i) {
            quality.put(i->first, i->second);
          }

          quality.name     = replace(quality.name, ' ', '_');

          buffer.insert(quality);           // only insert if absent
        }
      }
    }
    catch(const exception& error) {
      FATAL(error.what() << endl);
    }

  } else {                                  // read data from file

    JASCIIFileReader<JRunQuality> in(inputFile.c_str(), JDBDictionary::getInstance());

    JComment comment;

    in >> comment;

    for (JRunQuality quality; in >> quality; ) {
      buffer.insert(quality);
    }

    in.close();
  }


  if (buffer.empty()) {
    FATAL("No data." << endl);
  }


  if (getFilenameExtension(outputFile) == ROOT_FILE_FORMAT) {

    runs = JRange_t(buffer.begin()->run, buffer.rbegin()->run);

    TH1D h0("h0", NULL, 1000, 0.0, 1.01);
    TH1D h1("h1", NULL, condition.size(), -0.5, condition.size() - 0.5);

    for (size_t i = 0; i != condition.size(); ++i) {
      h1.GetXaxis()->SetBinLabel(i+1, condition[i].formula.c_str());
    }

    h0.SetMinimum(0.0);
    h0.SetMaximum(1.1 * (runs.getUpperLimit() - runs.getLowerLimit()));

    h1.SetMinimum(0.0);
    h1.SetMaximum(1.1 * (runs.getUpperLimit() - runs.getLowerLimit()));

    const Double_t xmin = runs.getLowerLimit() - 0.5;
    const Double_t xmax = runs.getUpperLimit() + 0.5;

    JManager<TString, TH1D> H1(new TH1D("%", NULL, runs.getLength() + 1, xmin, xmax));

    // show validity range

    for (JCondition_t::const_iterator i = condition.begin(); i != condition.end();++i) {
      H1[i->formula]->GetListOfFunctions()->Add(new TF1(MAKE_CSTRING("C " << i->formula << ":upper"), MAKE_CSTRING(i->range.getUpperLimit()), xmin, xmax));
      H1[i->formula]->GetListOfFunctions()->Add(new TF1(MAKE_CSTRING("C " << i->formula << ":lower"), MAKE_CSTRING(i->range.getLowerLimit()), xmin, xmax));
    }

    for (JCondition_t::const_iterator i = veto.begin(); i != veto.end();++i) {
      H1[i->formula]->GetListOfFunctions()->Add(new TF1(MAKE_CSTRING("V " << i->formula << ":upper"), MAKE_CSTRING(i->range.getUpperLimit()), xmin, xmax));
      H1[i->formula]->GetListOfFunctions()->Add(new TF1(MAKE_CSTRING("V " << i->formula << ":lower"), MAKE_CSTRING(i->range.getLowerLimit()), xmin, xmax));
    }

    ostringstream os;                       // N-tuple format

    os << "run";

    for (size_t i = 0; i != condition.size(); ++i) {
      os << ":" << (char) ('a' + i);
    }

    os << ":Q:V:R";

    TNtuple n1("n1", "quality", os.str().c_str());


    // process data

    for (set<JRunQuality>::const_iterator quality = buffer.begin(); quality != buffer.end(); ++quality) {

      vector<Float_t> tuple(1, (Float_t) quality->run);

      double w = 0.0;

      for (size_t i = 0; i != condition.size(); ++i) {

        const JCondition& ps = condition[i];
        const double      y  = getResult(ps.formula, *quality);

        DEBUG(ps.formula << ' ' << y << endl);

        tuple.push_back((Float_t) y);

        TH1D* p = H1[ps.formula];

        p->SetBinContent(p->FindBin((double) quality->run), y);
        p->SetBinError  (p->FindBin((double) quality->run), numeric_limits<double>::epsilon());

        if (ps.range(y)) {
          w += ps.weight;
        }

        h1.AddBinContent(i + 1, ps.range(y) ? 1.0 : 0.0);
      }

      const double Q = w/W;

      tuple.push_back((Float_t) Q);

      int V = 0;

      for (size_t i = 0; i != veto.size(); ++i) {

        const JCondition& ps = veto[i];
        const double      y  = getResult(ps.formula, *quality);
        
        DEBUG(ps.formula << ' ' << y << endl);

        TH1D* p = NULL;

        p = H1[ps.formula];

        p->SetBinContent(p->FindBin((double) quality->run), y);
        p->SetBinError  (p->FindBin((double) quality->run), numeric_limits<double>::epsilon());

        p = H1[MAKE_STRING("VETO[" << ps.formula << "]")];

        p->SetBinContent(p->FindBin((double) quality->run), ps.range(y) ? 0.0 : 1.0);
        p->SetBinError  (p->FindBin((double) quality->run), numeric_limits<double>::epsilon());

        if (!ps.range(y)) {
          ++V;
        }
      }

      tuple.push_back((Float_t) V);
      tuple.push_back((Float_t) setups.get(quality->run));

      h0.Fill(Q);
      n1.Fill(tuple.data());

      cout << setw(8) << quality->run << ' ' << FIXED(5,3) << Q << ' ' << setw(2) << V << endl;
    }

    double w = 0.0;

    for (Int_t i = 0; i <= h0.GetXaxis()->GetNbins(); ++i) {
      h0.SetBinContent(i, (w += h0.GetBinContent(i)));
    }

    // store results

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

    out << h0 << h1 << n1 << H1;
    
    out.Write();
    out.Close();
  }


  if (getFilenameExtension(outputFile) == ASCII_FILE_FORMAT) {

    // store data

    JASCIIFileWriter<JRunQuality> out(outputFile.c_str(), JDBDictionary::getInstance());

    out.setf(ios::fixed);

    JComment comment;

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

    out << comment;

    for (set<JRunQuality>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
      out << *i << endl;
    }

    out.close();
  }
}