JTuneTH.cc

Go to the documentation of this file.

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

#include "dbclient/KM3NeTDBClient.h"

#include "JDB/JDB.hh"
#include "JDB/JDBToolkit.hh"
#include "JDB/JDBincludes.hh"
#include "JDB/JProductRouter.hh"
#include "JDB/JSelector.hh"
#include "JDB/JSelectorSupportkit.hh"

#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorCalibration.hh"

#include "JSupport/JMeta.hh"

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

#include "JCalibrate/JParameters.hh"

#include "JSon/JSon.hh"

#include "JROOT/JManager.hh"

/**
 * \file
 * Auxiliary program to compute PMT thresholds from scan.
 * \author acreusot
 */


int main(const int argc, const char * const argv[])
{
  using namespace std;
  using namespace JPP;

  string usr;
  string pwd;
  string cookie;
  string detectorFile;
  string inputFile;
  string outputFile;
  double stfCut;
  string testType;
  int debug;

  try {

    JParser<> zap("Auxiliary program to compute the PMT thresholds according to the small ToT fraction (i.e. the noise to signal ratio).");

    zap['u'] = make_field(usr)        = "";
    zap['!'] = make_field(pwd)        = "";
    zap['C'] = make_field(cookie)     = "";
    zap['a'] = make_field(detectorFile, "detector file");
    zap['f'] = make_field(inputFile, "output file of JIntegrateToT");
    zap['o'] = make_field(outputFile, "output file for json");
    zap['t'] = make_field(stfCut, "small ToT fraction cut-off") = 0.5;
    zap['T'] = make_field(testType, "test type") = "TH-TUNING-SEA-v1";
    zap['d'] = make_field(debug)      = 1;
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  JDateAndTime timer;
  // retrieve the detector 
  JDetector detector;
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
  const JModuleRouter moduleRouter(detector);
  const int detId = detector.getID();
 
  JPersons person;
  vector<string> runNumbers;
  vector<JParameters> vParameters;

  // deal with the database
  try {
    NOTICE("connecting to the database..." << endl);
    JPMTParametersMap parameters;
    JDB::reset(usr, pwd, cookie);
    {
      JSelector selector = getSelector<JPersons>(JDB::get()->User());
      ResultSet& rs = getResultSet(getTable<JPersons>(), selector);
      rs >> person;
      rs.Close();
    }

    // retrieve the detector integration parameters from the DB
    vector<JDetectorIntegration> detectorInt;
    try {
      JDetectorIntegration parameters;
      ResultSet& rs = getResultSet(getTable<JDetectorIntegration>(), getSelector<JDetectorIntegration>(detId));
      while (rs >> parameters) {
        detectorInt.push_back(parameters);
      }
      rs.Close();
    }
    catch(const exception& error) {
      FATAL(error.what() << endl);
    }

    // read the input file and file up the JParameters vector
    ifstream istr(inputFile.c_str());
    int prevRunId = -1;
    for (JParameters curPar; istr >> curPar;) {
      if (curPar.runId != prevRunId) {
        NOTICE("Extracting thresholds for run " << curPar.runId << endl);
        runNumbers.push_back(to_string(curPar.runId));
      }
      // insert this JParameters in the vector of JParameters
      vParameters.push_back(curPar);
      prevRunId = curPar.runId;
    }
    istr.close();

    for (unsigned int r = 0; r < runNumbers.size(); ++r) {
      // retrieve the OID from the database with the detId and the run number
      string rs_oid;
      try {
        JRuns parameters;
        ResultSet& rs = getResultSet(getTable<JRuns>(), getSelector<JRuns>(detId, stoi(runNumbers[r])));
        if (rs >> parameters) {
          rs_oid = parameters.RUNSETUPID;
        }
        rs.Close();
        DEBUG("Run setup " << rs_oid << endl);
      }
      catch(const exception& error) {
        FATAL(error.what() << endl);
      }

      // retrieve OPTICS parameters (for whatever it means) from the database 
      JAllParams upars;
      try {
        JAllParams parameters;
        ResultSet& rs = getResultSet(getTable<JAllParams>(), getSelector<JAllParams>("OPTICS"));
        while (rs >> parameters) {
          if (parameters.NAME.find("PMT_THRESHOLD") != string::npos) {
            upars = parameters;
          }
        }
        rs.Close();
      }
      catch(const exception& error) {
        FATAL(error.what() << endl);
      }

      // retrieve the parameters from the database with the OID and detId 
      // and with a cut in the OPTICS
      map<string, int> umap;
      try {
        JRunsetupParams parameters;
        ResultSet& rs = getResultSet(getTable<JRunsetupParams>(), getSelector<JRunsetupParams>(getDetector(detId), rs_oid));
        while (rs >> parameters) {
          if (parameters.PAR_OID == upars.OID && parameters.ISINPUT == 'Y') {
            if (parameters.VALUE != "") {
              umap[parameters.UPIFILTER] = to_value<int>(parameters.VALUE);
            }
          }
        }
        rs.Close();
      }
      catch(const exception& error) {
        FATAL(error.what() << endl);
      }

      // match the threshold with the detector integration parameters 
      for (vector<JDetectorIntegration>::const_iterator i = detectorInt.begin(); i != detectorInt.end(); ++i) {
        if (i->PMTID != -1) {
          const JUPI_t upi = i->PMTUPI;
          map<string, int>::const_iterator p = umap.find(upi.toString());
          if (p == umap.end()) {
            ostringstream os;
            os << upi.getPBS() << JUPI_t::SEPARATOR
               << upi.getVariant() << JUPI_t::SEPARATOR
               << JUPI_t::DOT;
            p = umap.find(os.str());
          }
          if (p != umap.end()) {
            const JPMTIdentifier id(i->DOMID, i->CABLEPOS);
            parameters[id].threshold = p->second;
          } else {
            ERROR("Missing threshold data PMT " << upi << endl);
          }
        }
      }

      // fill the JParameters vector with threshold information
      for (unsigned int pmt = 0; pmt < vParameters.size(); ++pmt) {
        const int domId = vParameters[pmt].domId;
        const int pmtId = vParameters[pmt].pmtId;
        const JPMTIdentifier id(domId, pmtId);
        const double threshold = parameters[id].threshold;
        if (vParameters[pmt].runId == stoi(runNumbers[r])) {
          vParameters[pmt].threshold = threshold;
        }
      }
    }

    // read the JParameters vector to find optimal threshold
    JTHCalibration THcalibrations;
    int totalPmts = 0;
    int validPmts = 0;
    bool newPmt = true;
    int previousSerial = -1;
    double optth = -1;
    double defaultth = -1;
    bool curbc = true;
    bool calibratedPmt = false;
    sort(vParameters.begin(), vParameters.end());
    for (unsigned int pmt = 0; pmt < vParameters.size(); ++pmt) {
      const int domId = vParameters[pmt].domId;
      const int pmtId = vParameters[pmt].pmtId;
      const int pmtSerial = moduleRouter.getModule(domId).getPMT(pmtId).getID();
      if (pmtSerial != previousSerial) {
        newPmt = true;
        ++totalPmts;
      }
      previousSerial = pmtSerial;
      if ((calibratedPmt == true) && (newPmt == false)) continue;
      const int duId = moduleRouter.getModule(domId).getString();
      const int floorId = moduleRouter.getModule(domId).getFloor();
      const JLocation_t location(duId, floorId, pmtId);
      if (newPmt == true) {
        defaultth = vParameters[pmt].threshold;
        curbc = true;
        calibratedPmt = false;
      }
      const double curth = vParameters[pmt].threshold;
      const double noise = vParameters[pmt].noise;
      const double signal = vParameters[pmt].signal;
      const double ratio = signal ? noise/signal : 0;
      if ((ratio > stfCut) || (ratio == 0)) {
        curbc = true;
      } else {
        ++validPmts;
        curbc = false;
        optth = curth;
        calibratedPmt = true; 
      }
      if (curbc == true) {
        if (vParameters[pmt].runId == stoi(runNumbers.back())) {
          NOTICE("Bad channel for pmt " << pmtId << "(upi " 
                 << pmtSerial << "), floor " << floorId << "(dom " << domId << "), du " 
                 << duId << endl);
        }
      } else {
        if (optth > defaultth) {
          NOTICE("Non default thresholds " << optth << " for pmt " << pmtId << "(upi " 
                 << pmtSerial << "), floor " << floorId << "(dom " << domId << "), du " 
                 << duId << " with a STF of " << ratio << endl);
        }
      }
      if (calibratedPmt == true) { 
        const JUPI_t& pmtUPI = getUPI(PBS::PMT, pmtSerial);
        const JTHCalibration_t THcal(pmtUPI, OK_t, optth, runNumbers);
        THcalibrations.push_back(THcal);
      }
      newPmt = false;
    }
    json js;
    js[User_t] = person.LOGIN;
    js[Location_t] = person.LOCATIONID;
    js[StartTime_t] = timer.toString();
    js[EndTime_t] = timer().toString();
    js[Test_t + Type_t] = testType;
    js[Tests_t] = json(THcalibrations);
    ofstream ofs(outputFile.c_str());
    ofs << setw(2) << setprecision(8);
    ofs << js;
    ofs.close();
    NOTICE(validPmts << " calibrated PMTs among " << totalPmts << endl);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
}