JDataQuality.cc

Go to the documentation of this file.

#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#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/JRootPrinter.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JRootStreamObjectOutput.hh"
#include "JROOT/JRootDictionary.hh"
#include "JROOT/JManager.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 exoression and a range.
   */
  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
 *
 * Example program to plot quality data from data base.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

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

  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("Example program to plot 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['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, "Output file containing histograms and n-tuple.") = "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.*\\..*\\..*"), JGITTags_t::key_type("2019-04-12"));
    zap['Q'] = make_field(condition,  "User defined conditions");
    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 (W <= 0.0) {
    FATAL("Invalid total weight: " << W << endl);
  }

  JRunsetups       setup;

  set<JRunQuality> buffer;

  if (inputFile == "") {

    ASSERT(detid != "");

    int ID = -1;        // integer representation of detector identifier

    try {

      JDB::reset(usr, pwd, cookie);

      if (is_integer(detid))
        ID = to_value<int>(detid);
      else
        ID = getDetector(detid);

      if (is_integer(detid)) {
        detid = getDetector(to_value<int>(detid));
      }

      selection += getSelector<JRuns>(ID);

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

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

      rs.Close();

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

        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, *i);

        try {

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

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

          rs.Close();
        }
        catch(const exception& error) {}

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

          JRunQuality quality;

          quality.GIT      = *i;
          quality.detector = ID;
          quality.run      = run->first;

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

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

  } else {

    ifstream in(inputFile.c_str());

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

    in.close();
  }

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

  if (debug >= debug_t) {

    JRootWriter writer(cout, JEquationParameters(), JDBDictionary::getInstance());

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


  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);

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

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

  JManager<TString, TH1D> zmap(new TH1D("%", NULL, 
                                        runs.getLength() + 1,
                                        runs.getLowerLimit() - 0.5,
                                        runs.getUpperLimit() + 0.5));

  // show range

  for (JCondition_t::const_iterator i = condition.begin(); i != condition.end();++i) {

    const string buffer = i->formula;

    zmap[buffer]->GetListOfFunctions()->Add(new TF1(MAKE_CSTRING(buffer << ":upper"), MAKE_CSTRING(i->range.getUpperLimit()), runs.getLowerLimit(), runs.getUpperLimit()));
    zmap[buffer]->GetListOfFunctions()->Add(new TF1(MAKE_CSTRING(buffer << ":lower"), MAKE_CSTRING(i->range.getLowerLimit()), runs.getLowerLimit(), runs.getUpperLimit()));
  }

  for (JCondition_t::const_iterator i = veto.begin(); i != veto.end();++i) {

    const string buffer = MAKE_STRING("VETO[" << i->formula << "]");

    zmap[buffer]->GetListOfFunctions()->Add(new TF1(MAKE_CSTRING(buffer << ":upper"), MAKE_CSTRING(i->range.getUpperLimit()), runs.getLowerLimit(), runs.getUpperLimit()));
    zmap[buffer]->GetListOfFunctions()->Add(new TF1(MAKE_CSTRING(buffer << ":lower"), MAKE_CSTRING(i->range.getLowerLimit()), runs.getLowerLimit(), runs.getUpperLimit()));
  }

  ostringstream os;

  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());


  resetObject(&h0);
  resetObject(&h1);
  resetObject(&n1);
  resetObject(&zmap);

  // process data

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

    vector<Float_t> tuple;

    tuple.push_back((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 = zmap[ps.formula];

      p->SetBinContent(p->FindBin((double) quality->run), y);
      p->SetBinError  (p->FindBin((double) quality->run), 0.0);

      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);

      TH1D* p = zmap[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), 0.0);

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

    tuple.push_back((Float_t) V);
    tuple.push_back((Float_t) setup.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)));
  }
  
  // output

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

  out << h0 << h1 << n1 << zmap;

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