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