JTuneHV.cc

Go to the documentation of this file.

#include <time.h>
#include <string>
#include <iostream>
#include <iomanip>

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

#include "JSupport/JMeta.hh"

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

#include "dbclient/KM3NeTDBClient.h"

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

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

#include "JCalibrate/JFitToT.hh"

#include "JGizmo/JManager.hh"

#include "JSon/JSon.hh"

#include "JROOT/JRootToolkit.hh"

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


/**
 * \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 JMAPLIST<JMap,JMap,JMap,JMap>::maplist                        JMapList_t;
  typedef JMultiMap<int, JUPI, JMapList_t>                             JMultiMap_t;

  typedef typename JMultiMap_t::super_const_iterator                    JSuperIt_t;

  typedef pair<double,double>                                            element_t;
  typedef map<JUPI, map<double, element_t>>                              datamap_t;


  string           usr;
  string           pwd;
  string           cookie;
  map<int, string> inputFiles;
  int              detid;
  string           outputFile;
  string           graphFile;
  string           testType;
  string           testProduct;
  double           HVerr;
  string           option;
  int              debug;
 

  try {

    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 and associated ToT-fit root files (JFitToT output)");
    zap['D'] = make_field(detid,       "detector ID");
    zap['T'] = make_field(testType,    "test type")                                 = "HV-TUNING-SEA-GAIN-v1";
    zap['P'] = make_field(testProduct, "test product")                              = "";  // For testing the DB-uploading
    zap['e'] = make_field(HVerr,       "Error on high-voltage measurements [V]")    = 5.0; // [V]
    zap['O'] = make_field(outputFile,  "HV tuning calibration file (json format)");
    zap['o'] = make_field(graphFile,   "Graphical output           (root format)")  = "";
    zap['d'] = make_field(debug,       "debug")                                     = 1;

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

  // Start time
  time_t startTime = time(0);

  vector<string> runNumbers;
  
  JPersons     person;
  JMultiMap_t  locMap;
  datamap_t   dataMap;
  
  try {
    
    // Extract user info
    JDB::reset(usr, pwd, cookie);
    {
      JSelector  selector = getSelector<JPersons>(JDB::get()->User());
      ResultSet& rs       = getResultSet(getTable<JPersons>(), selector);
      rs >> person;
      rs.Close();
    }

    
    // Extract detector integration info
    {
      JSelector  selector = getSelector<JDetectorIntegration>(getDetector(detid));
      ResultSet& rs       = getResultSet(getTable<JDetectorIntegration>(), selector); 
      
      for (JDetectorIntegration detint; rs >> detint; ) {
        
        if (detint.PMTID >= 0) {
          
          const JUPI upi = detint.PMTUPI;
          
          locMap[detint.DOMID][detint.DUID][detint.FLOORID][detint.PMTID] = upi;
        }
      }
      
      rs.Close();
    }
    

    // Extract gain and high-voltage data
    for (map<int,string>::const_iterator fileIt = inputFiles.begin(); fileIt != inputFiles.end(); ++fileIt) {

      const int runNr = fileIt->first;

      NOTICE("Extracting Gain-/HV-data for run " << runNr << endl);
      runNumbers.push_back(to_string(runNr));

      // Retrieve high-voltage settings
      map<JUPI, double> HVmap;

      {
        JSelector  selector = getSelector<JPMTHVRunSettings>(getDetector(detid), fileIt->first);
        ResultSet& rs       = getResultSet(getTable<JPMTHVRunSettings>(), selector);

        
        for (JPMTHVRunSettings parameters; rs >> parameters; ) {
          
          for (JSuperIt_t locIt = locMap.super_begin(); locIt != locMap.super_end(); ++locIt) {

            if (parameters.DUID     == locIt->second->first &&
                parameters.FLOORID  == locIt->second->second->first &&
                parameters.PMTINTID == locIt->second->second->second->first) {
              
              const JUPI upi = locIt.getValue();
              
              HVmap[upi] = log(fabs(parameters.HV_VALUE));
            }
          }
        }
        
        rs.Close();
      }
      
      
      // Retrieve fitted gains
      TFile in(fileIt->second.c_str(), "READ");
      TIter next(in.GetListOfKeys());

      while(TKey* key = (TKey*)next()) {

        const string name(key->GetName());
        if (!gROOT->GetClass(key->GetClassName())->InheritsFrom("TH1") ||
            name.find("1ToT") == string::npos) { continue; }

        const int pos   = name.find(".");
        const int modID = stoi(name.substr(0, pos));
        const int pmtID = stoi(name.substr(pos+1, 2));

        for (JSuperIt_t locIt = locMap.super_begin(); locIt != locMap.super_end(); ++locIt) {
          
          if (modID == locIt->first &&
              pmtID == locIt->second->second->second->first) {
              
            const JUPI   upi      = locIt.getValue();
            const double logAbsHV = HVmap[upi];
        
            const TH1*   h1       = (TH1*) key->ReadObj();
            const TF1*   f1       = h1->GetFunction(FITTOT_FNAME.c_str());

            double gain, gainErr;
            if (f1 != NULL) {

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

              gain      = f1->GetParameter(Ngain);
              gainErr   = f1->GetParError (Ngain);

            } else {

              const double prevG = (dataMap[upi].size() > 0 ?
                                    (--(dataMap[upi].lower_bound(logAbsHV)))->second.first : log(FITTOT_GAIN_MIN));

              gain    = (prevG > log(NOMINAL_GAIN) ? FITTOT_GAIN_MAX : FITTOT_GAIN_MIN);
              gainErr = 0.0;
            }

            dataMap[upi][logAbsHV] = element_t{log(gain), gainErr / gain};
          } 
        }
      }
    }
  }
  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;

  JManager<string, TGraphErrors> manager(new TGraphErrors(), '%', ios::fmtflags(), MAKE_STRING("%_HVxG"));
  
  NOTICE("Searching optimal high-voltage settings...\n\t" <<
         "High-Voltage [V] --- Gain [-]" << endl);

  for (datamap_t::const_iterator pmt  = (testProduct.empty() ? dataMap.cbegin() :    dataMap.find(testProduct)) ;
                                 pmt != (testProduct.empty() ? dataMap.cend()   : ++(dataMap.find(testProduct)));
                               ++pmt) {

    const JUPI                   upi  = pmt->first;
    const map<double, element_t> data = pmt->second;

    string location, serialID;
    for (JSuperIt_t locIt = locMap.super_begin(); locIt != locMap.super_end(); ++locIt) {

      if (locIt.getValue().getVariant() == upi.getVariant()) {

        char buffer[100];

        sprintf(buffer, "%02d.%02d.%02d",
                locIt->second->first,
                locIt->second->second->first,
                locIt->second->second->second->first);
        location = string(buffer);

        sprintf(buffer, "%08d.%02d",
                locIt->first,
                locIt->second->second->second->first);
        serialID = string(buffer);
      }
    }

    NOTICE(upi << " (a.k.a. " << serialID << " / " << location << "):\n\t");

    if (data.size() < 2) {

      WARNING("Insufficient data; skip!" << endl);
      continue;
    }

    manager[serialID]->Set(data.size());
    manager[serialID]->SetMarkerStyle(kFullDotSmall);

    const double logNominalGain = log(NOMINAL_GAIN);
    
    const double D1    = fabs(logNominalGain - log(FITTOT_GAIN_MAX));
    const double D2    = fabs(logNominalGain - log(FITTOT_GAIN_MIN));
    double minDistance = (D2 > D1 ? D2 : D1);
    
    map<double, element_t>::const_iterator k0 = data.cbegin();
    map<double, element_t>::const_iterator k1 = data.cbegin();
    
    for (map<double, element_t>::const_iterator it = data.cbegin(); it != data.cend(); ++it) {

      // Plot point
      const int dIt = distance(data.cbegin(), it);
      manager[serialID]->SetPoint(dIt, it->first / log(10), it->second.first / log(10));
      manager[serialID]->SetPointError(dIt, fabs(HVerr / exp(it->first) / log(10)), fabs(it->second.second / log(10)));
      
      // Find the point closest to the nominal gain-value (k0)
      // and a compatible auxiliary point (k1) for inter-/extrapolation
      const double distance = fabs(logNominalGain - it->second.first);
      
      if (distance < minDistance) {
        
        k0 = it;
        
        k1 = ( it->second.first < logNominalGain ?
              (it != --data.cend() ? next(it,1) : prev(it,1)) : // Closest point below nominal gain
              (it != data.cbegin() ? prev(it,1) : next(it,1)) ); // Closest point above nominal gain
        
        minDistance = distance;
      }
    }
    
    // Determine optimal high-voltage via linear inter-/extrapolation
    bool valid = false;
    pair<double, double> result{-exp(k0->first), exp(k0->second.first)};
    pair<double, double> errors{          HVerr, exp(k0->second.second)};
    
    if ((k1->first < k0->first && k1->second.first < k0->second.first) ||
        (k1->first > k0->first && k1->second.first > k0->second.first)) {

      // Logarithmic supply-voltage, logarithmic gain and logarithmic errors
      // for point closest to the nominal gain (0) and for the auxiliary point (1)
      const double logHV0       = k0->first;
      const double logHV1       = k1->first;

      const double logHV0err    = HVerr / exp(logHV0);
      const double logHV1err    = HVerr / exp(logHV1);
      
      const double logGain0     = k0->second.first;
      const double logGain1     = k1->second.first;
      
      const double logGain0err  = k0->second.second;
      const double logGain1err  = k1->second.second;

      // Interpolate the nominal logarithmic supply-voltage and compute error estimate
      const double slope        = (logGain1 - logGain0) / (logHV1 - logHV0); 
      const double logNominalHV = (logNominalGain - logGain0) / slope + logHV0;

      const double err1 = logHV0err * logHV0err;
      const double err2 = logGain0err * logGain0err / slope / slope;
      const double err3 = ( pow(logNominalGain - logGain0 / slope, 2.0) *
                            ((logGain0err*logGain0err + logGain1err*logGain1err) / pow(logGain1 - logGain0, 2.0) +
                             (logHV0err*logHV0err + logHV1err*logHV1err) / pow(logHV1 - logHV0, 2.0)) );
      
      const double logNominalHVerr = sqrt(err1 + err2 + err3);
      
      if (fabs(logNominalHV - logHV0) / fabs(logHV1 - logHV0) < 2.0) {
        
        valid = true;

        result.first  = -exp(logNominalHV);
        errors.first  =  exp(logNominalHV) * logNominalHVerr;
        result.second =  NOMINAL_GAIN;
        errors.second =  0.0;
      }

      if ((k0 == data.cbegin() && logNominalGain < logGain0) ||
          (k0 == --data.cend() && logNominalGain > logGain0)) {
        
        WARNING("Nominal gain " << NOMINAL_GAIN << " lies outside of HV-scan range; " <<
                "Extrapolation " << (valid ? "successful!" : "failed; Setting closest values!") << "\n\t");
      }

    } else {
    
      WARNING("Could not find optimal high-voltage; Setting closest values!\n\t");
    }

    NOTICE(setw(8) << setprecision(2) << fixed << left <<
           result.first  << setw(0) << " +/- " << setw(4) << errors.first << setw(10) << right << 
           result.second << setw(0) << " +/- " << setw(3) << errors.second << endl);
    
    const JHVCalibration_t HVcal(upi,
                                 valid ? OK_t : Failed_t,
                                 result.first,
                                 result.second,
                                 runNumbers);
    
    HVcalibrations.push_back(HVcal);
  }


  if (! graphFile.empty()) {

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

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

    out << manager;

    out.Close();
  }
  
  time_t stopTime = time(0);

  // Write json output
  json js;

  js[User_t]                  = person.LOGIN;
  js[Location_t]              = person.LOCATIONID;
  js[StartTime_t]             = get_date_string(startTime);
  js[EndTime_t]               = get_date_string(stopTime);
  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();
}