JTuneHV.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>

#include "JSystem/JTime.hh"

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

#include "JSupport/JMeta.hh"

#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/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorCalibration.hh"

#include "JCalibrate/JFitToT.hh"

#include "JTools/JMap.hh"
#include "JTools/JMultiMap.hh"

#include "JROOT/JManager.hh"

#include "JSon/JSon.hh"

#include "TFile.h"
#include "TGraphErrors.h"


namespace {

  using namespace std;
  using namespace JPP;


  /**
   * Auxiliary data structure for interpolating (HV, gain)-data
   */
  struct interpolator {

    typedef map<double, pair<double, double>>    map_t;

    /**
     * Constructor
     *
     * \param  minHVdiff           minimal difference between two HV vlaues
     * \param  minGain             minimal valid gain value
     * \param  maxGain             maximal valid gain value
     * \param  targetGain          target gain value for interpolation
     */
    interpolator(const double minHVdiff, const double minGain, const double maxGain, const double targetGain) :
      dHV(minHVdiff), gainMin(minGain), gainMax(maxGain), gainTarget(targetGain) {}

    
    /*
     * Checks whether the HV-values of two map entries are different
     *
     * \param entry1              first  map entry
     * \param entry2              second map entry
     * \return                    true if absolute HV difference greater than dHV; else false
     */
    const bool checkHV(map_t::const_reference entry1, map_t::const_reference entry2) const {

      return (fabs(entry1.first - entry2.first) > dHV);
    }


    /*
     * Checks whether the gain of a map entries lies within the specified range
     *
     * \param entry               map entry
     * \return                    true if gain within gainRange; else false
     */
    const bool checkGain(map_t::const_reference entry) const {

      return (entry.second.first > gainMin && entry.second.first < gainMax);
    }


    /*
     * Checks whether the gain of two map entries are strictly increasing as function of the absolute HV
     *
     * \param entry1              first  map entry
     * \param entry2              second map entry
     * \return                    true if absolute HV difference greater than dHV; else false
     */
    const bool areIncreasing(map_t::const_reference entry1, map_t::const_reference entry2) const { 
      
      return ((entry1.first < entry2.first && entry1.second.first > entry2.second.first) ||
              (entry1.first > entry2.first && entry1.second.first < entry2.second.first));
    }


    /*
     * Checks whether two map entries are valid for interpolation (i.e. different HV and gain within specified range)
     *
     * \param entry1              first  map entry
     * \param entry2              second map entry
     * \return                    true if absolute HV difference greater than dHV; else false
     */
    const bool areValid(map_t::const_reference entry1,
                        map_t::const_reference entry2) const { 
      
      return (checkHV      (entry1, entry2) && checkGain(entry1) &&
              areIncreasing(entry1, entry2) && checkGain(entry2));
    }


    /*
     * Get interpolation points
     *
     * \param  data               HV-gain data
     * \return                    pair of interpolation points
     */
    const pair<map_t::const_reference,
               map_t::const_reference> getPoints(const map_t& data) const {

      map_t::const_iterator point0 = data.cbegin();
      map_t::const_iterator point1 = next(data.cbegin());
      
      for (map_t::const_iterator it = next(data.cbegin()); it != data.cend(); ++it) {
        
        if        (fabs(it->second.first - gainTarget) < fabs(point0->second.first - gainTarget)) {
          
          point1 = (areValid(*point0, *it) ? point0 : point1);
          point0 = it;
          
        } else if ((fabs(it->second.first - gainTarget) < fabs(point1->second.first - gainTarget) ||
                    !areValid(*point0, *point1)) && areValid(*point0, *it)) { 
          
          point1 = it;
        }
      }

      return pair<map_t::const_reference,
                  map_t::const_reference>(*point0, *point1);
    }
    

     /*
     * Inter-/Extrapolates the nominal high-voltage, given two map entries
     *
     * \param entry1              first  map entry
     * \param entry2              second map entry
     * \return                    inter-/extrapolated nominal high-voltage and error
     */
    const pair<double, double> getNominalHV(map_t::const_reference entry1,
                                            map_t::const_reference entry2) const {

      pair<double, double> result{0.0, 0.0};
      
      if (areValid(entry1, entry2)) {
        
        const double logHV0    = log(fabs(entry1.first));
        const double logHV1    = log(fabs(entry2.first));

        const double logG0     = log(entry1.second.first);
        const double logG1     = log(entry2.second.first);
        const double elogG0    = entry1.second.second / entry1.second.first;
        const double elogG1    = entry2.second.second / entry2.second.first;
      
        const double dlogG0    = log(gainTarget) - logG0;
        const double dlogG1    = log(gainTarget) - logG1;
      
        const double slope     = (logG1 - logG0) / (logHV1 - logHV0);
        const double logHVnom  = dlogG0 / slope + logHV0;
        const double elogHVnom = sqrt(dlogG1 * dlogG1 * elogG0 * elogG0 +
                                      dlogG0 * dlogG0 * elogG1 * elogG1) / fabs(slope * (logG1 - logG0));

        result.first  = -exp(logHVnom);
        result.second =  exp(logHVnom) * elogHVnom;
      }
      
      return result;
    }


    /*
     * Inter-/Extrapolates the nominal high-voltage, given HV-gain data
     *
     * \param data                HV-gain data
     * \return                    inter-/extrapolated nominal high-voltage and error
     */
    const pair<double, double> getNominalHV(const map_t& data) const {

      pair<map_t::const_reference,
           map_t::const_reference> intpPoints = getPoints(data);

      return getNominalHV(intpPoints.first, intpPoints.second);
    }
     
    
    const double dHV;         //!< minimal HV step value [V]
    const double gainMin;     //!< Minimal valid gain value
    const double gainMax;     //!< Maximal valid gain value
    const double gainTarget;  //!< target gain value for interpolation
  };
}



/**
 * \file
 *
 * Example program to fit the gain to the high-voltage
 * \author acreusot, mdejong and bjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef typename JMAPLIST<JMap, JMap, JMap, JMap, JMap>::maplist           mapList_t;

  typedef map<double, pair<double, double>>                                   runMap_t;
  typedef JMultiMap<int, runMap_t, mapList_t>                                 pmtMap_t;

  string           usr;
  string           pwd;
  string           cookie;
  
  map<int, string> inputFiles;
  string           detectorFile;
  string           outputFile;
  string           graphFile;
  string           testType;

  double           HVstepMin   = 2 * 3.14;                //!< Minimal HV step [V]
  double           gainMin     = FITTOT_GAIN_MIN + 0.1;   //!< Minimal valid gain fit-value
  double           gainMax     = FITTOT_GAIN_MAX - 0.1;   //!< Maximal valid gain fit-value
  double           gainTarget  = NOMINAL_GAIN;            //!< Target nominal gain value

  int              debug;
 

  try {

    JProperties properties;

    properties.insert(gmake_property(HVstepMin));
    properties.insert(gmake_property(gainMin));
    properties.insert(gmake_property(gainMax));
    properties.insert(gmake_property(gainTarget));
    
    JParser<> zap("Example program to find optimal HV-settings.");

    zap['u'] = make_field(usr)                                                        = "";
    zap['!'] = make_field(pwd)                                                        = "";
    zap['C'] = make_field(cookie)                                                     = "";
    zap['f'] = make_field(inputFiles,   "map of run numbers to file strings (i.e. to the output of JFitToT)");
    zap['a'] = make_field(detectorFile, "detector file");
    zap['T'] = make_field(testType,     "test type")                                  = "HV-TUNING-SEA-GAIN-v1";
    zap['O'] = make_field(outputFile,   "HV tuning calibration file (json format)");
    zap['o'] = make_field(graphFile,    "Graphical output           (root format)")   = "";
    zap['@'] = make_field(properties)                                                 = JPARSER::initialised();
    zap['d'] = make_field(debug,        "debug")                                      = 1;

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

  JDateAndTime           timer;
  
  JDetector              detector;
  JPersons               person;
  
  vector<string>         runNumbers;
  vector<int>            missingDOMs;
  pmtMap_t               data;

  JManager<string, TGraphErrors> manager(new TGraphErrors(), MAKE_STRING("%.HVxG"));  
    
  try {
    
    NOTICE("Loading database tables..." << endl);

    // Load detector file
    load(detectorFile, detector);

    JDB::reset(usr, pwd, cookie);
    
    // Extract user info
    {
      JSelector  selector = getSelector<JPersons>(JDB::get()->User());
      ResultSet& rs       = getResultSet(getTable<JPersons>(), selector);
      rs >> person;
      rs.Close();
    }
    
    // Extract gain and high-voltage data
    for (map<int, string>::const_iterator fileIt = inputFiles.cbegin(); fileIt != inputFiles.cend(); ++fileIt) {

      const int runNr = fileIt->first;
      
      NOTICE("Extracting Gain-/HV-data for run " << runNr << endl);
      runNumbers.push_back(to_string(runNr));

      {
        TFile fitData(fileIt->second.c_str(), "READ");
        
        JSelector  selector = getSelector<JPMTHVRunSettings>(getDetector(detector.getID()), runNr);
        ResultSet& rs       = getResultSet(getTable<JPMTHVRunSettings>(), selector);
        
        for (JPMTHVRunSettings table; rs >> table; ) {
          
          const int domID = detector.getModule(JLocation(table.DUID, table.FLOORID)).getID();
          
          const TH1* h1 = (TH1*) fitData.Get(MAKE_CSTRING(right << domID << FITTOT_SEPARATOR << FILL(2,'0') <<
                                                          table.PMTINTID << FITTOT_SEPARATOR << FITTOT_SUFFIX));
          const TF1* f1 = (h1 != NULL ? h1->GetFunction(FITTOT_FNAME.c_str()) : NULL);

          double gain      = (fabs(gainMin - gainTarget) > fabs(gainMax - gainTarget) ? gainMin : gainMax);
          double gainError = 0.0;
          
          if (f1 != NULL) {

            const int Ngain = f1->GetParNumber("gain");

            gain      = f1->GetParameter(Ngain);
            gainError = f1->GetParError (Ngain);
            
          } else if (h1 == NULL && find(missingDOMs.cbegin(), missingDOMs.cend(), domID) == missingDOMs.cend()) {
            
            missingDOMs.push_back(domID);
          }
          
          data[domID][table.PMTSERIAL][table.DUID][table.FLOORID][table.PMTINTID].insert(make_pair(table.HV_VALUE, make_pair(gain, gainError)));

          const string name  = MAKE_STRING(domID << '.' << FILL(2,'0') << table.PMTINTID);
          const int    index = distance(inputFiles.cbegin(), fileIt);
          
          manager[name]->SetPoint     (index, log10(fabs(table.HV_VALUE)), log10(gain));
          manager[name]->SetPointError(index, 0.0, gainError / gain / log(10));
        }
        
        rs.Close();
      }
    }
    
  } catch (const exception& error) {
    
    FATAL(error.what() << endl);
  }


  // Find optimal HV corresponding to a nominal gain of 1.0
  // via linear interpolation in log-log scale
  JHVCalibration HVcalibrations;
  
  NOTICE("Searching optimal high-voltage settings..." << endl <<
         LEFT(30) << "UPI" << CENTER(35) << "(identifer / location)" << RIGHT(35) << "Inter-/Extrapolated high-voltage [-]" << endl);

  for (pmtMap_t::super_const_iterator pmt = data.super_begin(); pmt != data.super_end(); ++pmt) {

    // Retrieve upi, location and serialID
    const int         domID   = pmt->first;
    const JUPI_t&     pmtUPI  = getUPI(PBS::PMT, pmt->second->first);
    const JLocation_t location(pmt->second->second->first,
                               pmt->second->second->second->first,
                               pmt->second->second->second->second->first);

    const runMap_t&   data    = pmt.getValue();
    
    NOTICE(LEFT(30) << pmtUPI << "(a.k.a. " << domID << '.' << FILL(2,'0') << location.position << " / " << location << "):");    
    
    // Check number of entries
    if (data.size() < 2 || find(missingDOMs.cbegin(), missingDOMs.cend(), domID) != missingDOMs.cend()) {

      NOTICE (setprecision(2) << right << FIXED(20,2) << 0.0 << " +/- " << FIXED(4,2) << 0.0 << endl);
      WARNING((data.size() < 2 ? "Insufficient data" : "Missing (dead?) DOM") << endl);

      HVcalibrations.push_back(JHVCalibration_t(pmtUPI, Failed_t, 0.0, 0.0, runNumbers));
      continue;
    }
    
    // Find interpolation points
    const interpolator interpolator(HVstepMin, gainMin, gainMax, gainTarget);
    
    const pair<runMap_t::const_reference,
               runMap_t::const_reference> points = interpolator.getPoints(data);
    const pair<double, double>            result = interpolator.getNominalHV(points.first, points.second);
    
    const bool valid = (interpolator.areValid(points.first, points.second) && result.first < 0.0 &&
                        ((log(fabs(result.first))        - log(fabs(points.first.first))) /
                         (log(fabs(points.second.first)) - log(fabs(points.first.first))) < 2.0));
                        
    NOTICE(setprecision(2) << right << FIXED(20,2) << result.first  << " +/- " << FIXED(4,2) << result.second << endl);

    if (result.first < data.cbegin()  -> first ||
        result.first > data.crbegin() -> first) {

      WARNING("Nominal gain " << gainTarget << " lies outside of HV-scan range; " <<
              "Extrapolation " << (valid ? "successful." : "failed; Setting zero.") << endl);
      
    } else if (!valid) {

      WARNING("Could not find optimal high-voltage; Setting zero." << endl);
    }
           
    const JHVCalibration_t HVcal(pmtUPI,
                                 valid ? OK_t : Failed_t,
                                 valid ? result.first : 0.0,
                                 valid ? NOMINAL_GAIN : 0.0,
                                 runNumbers);
    
    HVcalibrations.push_back(HVcal); 
  }


  if (! graphFile.empty()) {

    TFile out(graphFile.c_str(), "RECREATE");

    putObject(&out, JMeta(argc, argv));

    for (JManager<string, TGraphErrors>::const_iterator it = manager.cbegin(); it != manager.cend(); ++it) {
      it->second->SetMarkerStyle(kFullDotSmall);
    }
    
    out << manager;

    out.Close();
  }

  // Write json output
  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(HVcalibrations);

  ofstream ofs(outputFile.c_str());

  ofs << setw(2) << setprecision(8);
  ofs << js;

  ofs.close();
}