JTuneHV.cc File Reference

Example program to fit the gain to the high-voltage. More...

#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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Example program to fit the gain to the high-voltage.

Author

acreusot, mdejong and bjung

Definition in file JTuneHV.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

< Minimal HV step [V]

< Minimal valid gain fit-value

< Maximal valid gain fit-value

< Target nominal gain value

Definition at line 216 of file JTuneHV.cc.

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