JTunePMTThreshold.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <vector>

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

#include "km3net-dataformat/online/JDAQHeader.hh"

#include "JDB/JDB.hh"
#include "JDB/JSelector.hh"
#include "JDB/JSelectorSupportkit.hh"
#include "JDB/JDBToolkit.hh"
#include "JDB/JPersons.hh"
#include "JDB/JPMTThreshold.hh"
#include "JDB/JDetectorIntegration.hh"
#include "JCalibrate/JCalibrateToT.hh"
#include "JDetector/JDetectorCalibration.hh"
#include "JTools/JRange.hh"
#include "JLang/JLangToolkit.hh"
#include "JLang/JComparator.hh"
#include "JROOT/JRootFileReader.hh"
#include "JSon/JSon.hh"

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

namespace {

  using KM3NETDAQ::JDAQHeader;
  using JLANG::JException;

  /**
   * Get DAQ header from file.
   *
   * \param  file_name      file name
   * \return                DAQ header
   */
  inline JDAQHeader getDAQHeader(const char* const file_name)
  {
    using namespace JPP;

    JRootFileReader<JDAQHeader> in(file_name);

    if (in.hasNext()) {

      const JDAQHeader* p = in.next();

      if (!in.hasNext()) {
        return *p;
      }
    }

    THROW(JException, "Error reading DAQ header from file " << file_name);
  }
}


/**
 * \file
 * Auxiliary program to compute the PMT thresholds according to the small time-over-threshold fraction.
 *
 * \author acreusot, mdejong
 */
int main(const int argc, const char * const argv[])
{
  using namespace std;
  using namespace JPP;

  typedef JRange<double> JRange_t;

  JServer         server;
  string          usr;
  string          pwd;
  string          cookie;
  vector<string>  inputFile;
  string          outputFile;
  JRange_t        range;
  double          fraction;
  string          testType;
  int             debug;

  try {

    JParser<> zap("Auxiliary program to compute the PMT thresholds according to the small time-over-threshold fraction.");

    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,         "list of file names (output of JCalibrateToT)");
    zap['o'] = make_field(outputFile,        "output file for JSon")                          = "";
    zap['x'] = make_field(range,             "integration range for noise.")                  = JRange_t(0.0, 8.0);
    zap['t'] = make_field(fraction,          "maximal fraction of signal allowed for noise.") = 0.5;
    zap['T'] = make_field(testType,          "test type")                                     = "TH-TUNING-SEA-v1";
    zap['d'] = make_field(debug,             "debug.")                                        = 2;

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


  if (inputFile.empty()) {
    FATAL("No input files.");
  }

  JDateAndTime timer;
  JPersons     person;
  
  try {

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

    ResultSet& rs = getResultSet(getTable<JPersons>(), getSelector<JPersons>(JDB::get()->User()));

    rs >> person;

    rs.Close();
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  const JDAQHeader header = getDAQHeader(inputFile[0].c_str());
  const int        ID     = header.getDetectorID();

  typedef map<int, JDetectorIntegration>  module_type;
  typedef map<int, module_type>           detector_type;

  detector_type detector;

  try {

    ResultSet& rs = getResultSet(getTable<JDetectorIntegration>(), getSelector<JDetectorIntegration>(ID));

    for (JDetectorIntegration parameters; rs >> parameters; ) {
      detector[parameters.DOMID][parameters.CABLEPOS] = parameters;
    }

    rs.Close();
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  struct parameters_type {
    int    threshold;
    double signal;
    double noise;
  };

  typedef vector<parameters_type>  data_type;
  typedef map<JUPI_t, data_type>   map_type;

  map_type     data;
  vector<int>  runs;

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

    const JDAQHeader header = getDAQHeader(i->c_str());

    ASSERT(header.getDetectorID() == ID, "file: " << *i);

    runs.push_back(header.getRunNumber());

    const JPMTThreshold getPMTThreshold(header.getDetectorID(), header.getRunNumber());

    TFile* in = TFile::Open(i->c_str(), "exist");

    if (in == NULL || !in->IsOpen()) {
      FATAL("File: " << *i << " not opened." << endl);
    }

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

      if (!module->second.empty()) {

        TH2D* h2s = (TH2D*) in->Get(MAKE_CSTRING(module->first << _2SToT));

        if (h2s != NULL) {

          for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {

            double noise  = 0.0;
            double signal = 0.0;

            for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {

              const double x = h2s->GetYaxis()->GetBinCenter(iy);
              const double y = h2s->GetBinContent(ix, iy);

              signal += y;

              if (range(x)) {
                noise += y;
              }
            }

            const JUPI_t upi = module->second.at(ix-1).PMTUPI;

            data[upi].push_back({ getPMTThreshold(upi).value, signal, noise});
          }

        } else {

          WARNING("No histogram for module " << module->first << "; skip." << endl);
        }
      }
    }

    in->Close();
  }


  JPMTThresholdCalibration calibration;

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

    sort(i->second.begin(), i->second.end(), make_comparator(&parameters_type::threshold));

    int threshold = i->second.begin()->threshold;

    for (data_type::const_iterator p = i->second.begin(); p != i->second.end(); ++p) {

      if (p->noise <= fraction * p->signal) {

        threshold = p->threshold;

        break;
      }
    }

    if (debug >= debug_t || threshold > i->second.begin()->threshold) {

      cout << "PMT " << left << setw(32) << i->first << " -> " << right << setw(3) << threshold << " (" << i->second.begin()->threshold << ")" << endl;

      for (data_type::const_iterator p = i->second.begin(); p != i->second.end(); ++p) {
        DEBUG(setw(3) << p->threshold << ' ' << FIXED(7,0) << p->noise << "/" << FIXED(7,0) << p->signal << ' ' << (p->noise <= fraction * p->signal) << endl);
      }
    }
     
    calibration.push_back(JPMTThresholdCalibration_t(i->first.toString(), OK_t, threshold, runs));
  }


  if (outputFile != "") {

    json js;

    js[User_t]           = person.LOGIN;
    js[Location_t]       = person.LOCATIONID;
    js[Start_t + Time_t] = timer.toString();
    js[End_t   + Time_t] = timer().toString();
    js[Test_t  + Type_t] = testType;
    js[Tests_t]          = json(calibration);

    ofstream out(outputFile.c_str());

    out << setw(2) << setprecision(8);
    out << js;
    
    out.close();
  }
}