JHVTuning.cc

Go to the documentation of this file.

#include <string>
#include <vector>
#include <ctime>
#include <sys/time.h>
#include <unistd.h>

#include "JDAQ/JDAQTimesliceIO.hh"

#include "JTrigger/JTriggerParameters.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "Jeep/JTimer.hh"
#include "JLang/JString.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JROOT/JROOTClassSelector.hh"

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

#include "JDetector/JModuleRouter.hh"

#include "TH1D.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TAxis.h"
#include "TCanvas.h"
#include <TStyle.h>
#include <TLegend.h>
#include <TLegendEntry.h>

#include "dbclient/KM3NeTDBClient.h"
#include "JDB/JDB.hh"
#include "JDB/JSelector.hh"
#include "JDB/JDatabaseObjectIterator.hh"
#include "JDB/JDBincludes.hh"
#include "JDB/JDBToolkit.hh"


/**
 * \file
 * Program for production of the xml file with nominal high voltage
 * values. The calculation is based on the adjustment of the ToT
 * duration to a value of 26.4ns by scanning the high voltage. Inputs
 * are the root datafile list for each of the high voltage values
 * (option -f), the detector file (option -a), the username and password 
 * (option -u) for access to the database.
 * the detector file can be retrieved in the database:
 * https://km3netdbweb.in2p3.fr/detx/D_ARCA003?
 * 
 * \author acreusot
 */

// retrieve information in the database
using namespace std;
using namespace JPP;
using namespace KM3NETDAQ;


inline
void AccessDB(vector<vector<vector<string> > >& vvvupi, vector<vector<vector<double> > >& vvvHv, string& cookie, string& usr, string& pwd, string& query, const int detid, const string runId, const set<int>& lineIds, string& location) {
  try {
    JDB::reset(usr, pwd, cookie);
    usr = JDB::get()->User();
    if (runId == "-2") 
      cout << "opening connection as user " << usr << " for detector information" << endl;
    else if (runId == "-1") 
      cout << "opening connection as user " << usr << " for default HV values" << endl;
    else 
      cout << "opening connection as user " << usr << " for run " << runId << endl;
    const string detoid = getDetector(detid);
    JSelector selection;
    selection.add("detoid", detoid);
    if (query == "detectors") {
      //selection.add("detid", detid);
      JDatabaseObjectIterator<JDetectors> in(query, selection);
      while (in.hasNext()) {
          JDetectors* curdet = in.next();
          string curoid = curdet->OID;
          if (curoid == detoid) {
            location = curdet->LOCATIONID;
          }
      }
    } else if (query == "pmt_hv_run_settings") {
      selection.add("run", runId);
      JDatabaseObjectIterator<JPMTHVRunSettings> in(query, selection);
      while (in.hasNext()) {
        JPMTHVRunSettings* rs = in.next();
        int lineId = 0;
        for (set<int>::iterator it = lineIds.begin(); it != lineIds.end(); ++it) {
          if (*it == rs->DUID) break;
          ++lineId;
        }
        vvvHv[lineId][rs->FLOORID - 1][rs->PMTINTID] = rs->HV_VALUE;
      }
    } else if (query == "pmt_hv_settings") {
      JDatabaseObjectIterator<JPMTHVSettings> in(query, selection);
      while (in.hasNext()) {
        JPMTHVSettings* rs = in.next();
        int lineId = 0;
        for (set<int>::iterator it = lineIds.begin(); it != lineIds.end(); ++it) {
          if (*it == rs->DUID) break;
          ++lineId;
        }
        vvvHv[lineId][rs->FLOORID - 1][rs->PMTINTID] = rs->HV_VALUE;
        vvvupi[lineId][rs->FLOORID - 1][rs->PMTINTID] = rs->PMTUPI;
      }
    }
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
}


// solve the 2nd order polynom. If delta is negatif use a 1st order polynom.
inline
double SolveAnalytic(const double deg0, const double deg1, const double deg2, const double tot) {
  const double delta = deg1*deg1 - 4*deg2*(deg0 - tot);
  if ((delta < 0) || (deg2 == 0)) return deg1? (tot - deg0)/deg1 : 0;
  const double x1 = (-deg1 + sqrt(delta))/(2*deg2);
  const double x2 = (-deg1 - sqrt(delta))/(2*deg2);
  if (deg2 > 0) return (x1 < x2)? x1 : x2;
  else return (x1 < x2)? x2 : x1;
  //return (fabs(x1) < fabs(x2))? x1 : x2;
}

// solve the equation, whatever the function
inline
pair<double, double> Solve(vector<double>& vParameter, const string& curFunction, const double tot) {
  pair<double, double> curSol;
  if (curFunction == "pol1") {
    curSol.first = SolveAnalytic(vParameter[1], vParameter[2], 0, tot);
    curSol.second = vParameter[2]*curSol.first + vParameter[1];
  } else if (curFunction == "pol2") {
    curSol.first = SolveAnalytic(vParameter[1], vParameter[2], vParameter[3], tot);
    curSol.second = vParameter[3]*pow(curSol.first, 2) + vParameter[2]*curSol.first + vParameter[1];
  } else if (curFunction == "[0] + sqrt([1] - [2]*x)") {
    curSol.first = (vParameter[2] - pow(tot - vParameter[1], 2))/vParameter[3];
    curSol.second = vParameter[1] + sqrt(vParameter[2] - vParameter[3]*curSol.first);
  } else {
    curSol.first = -1;
    curSol.second = -1;
  }
  return curSol;
}

// optimize reduced chi2 and fill the ToT vs HV graph 
inline
void FillGraph(TGraphErrors*& graph, TH1D*& hChi2, TH1D*& histo, const double hv, int& localCounter) {
  const int firstBin = 15;
  const int lastBin = histo->GetXaxis()->GetLast() - 60;
  histo->GetXaxis()->SetRange(firstBin, lastBin);
  if (!histo->GetEntries()) return;
  if (!histo->Integral()) return;
  const int maxBin = histo->GetMaximumBin();
  double bestRedChi2 = 1e6;
  double bestMean = 0;
  for (int binInf = maxBin - 6; binInf < maxBin - 2; ++binInf) {
    for (int binSup = maxBin + 1; binSup <= maxBin + 6; ++binSup) {
      if (binSup - binInf < 4) continue;
      TF1* gausFit = new TF1("gausFit", "gaus", binInf, binSup);
      const double curInteg = histo->Integral(binInf, binSup);
      if (curInteg == 0) continue;
      histo->Fit(gausFit, "RQ");
      if (!gausFit->GetNDF()) {
        gausFit->Delete();
        continue;
      }
      const double fitMean = gausFit->GetParameter(1);
      const double fitSigma = gausFit->GetParameter(2);
      const double reducedChi2 = gausFit->GetChisquare()/gausFit->GetNDF();
      if ((std::fabs(reducedChi2 - 1) < std::fabs(bestRedChi2 - 1)) && (reducedChi2 > 0.5) && (fitSigma < 10)) {
        bestRedChi2 = reducedChi2;
        bestMean = fitMean;
      }
      gausFit->Delete();
    }
  }
  hChi2->Fill(bestRedChi2);
  histo->GetXaxis()->SetRange(0, lastBin);
  if (bestMean > 0) {
    graph->SetPoint(localCounter, hv, bestMean);
    ++localCounter;
  }
}

// find the HV range used to perform the scan
inline 
void FindHVRange(vector<vector<vector<pair<double, double> > > >& vvvHvRange, vector<vector<vector<vector<double> > > >& vvvRunHv, vector<vector<vector<double> > >& vvvHv) {
  for (unsigned int line = 0; line < vvvHvRange.size(); ++line) {
    for (unsigned int dom = 0; dom < vvvHvRange[line].size(); ++dom) {
      for (unsigned int pmt = 0; pmt < vvvHvRange[line][dom].size(); ++pmt) {
        double minHvDiff = 1e8;
        double maxHvDiff = 0;
        for (unsigned int run = 0; run < vvvRunHv.size(); ++run) {
          if ((vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt]) < minHvDiff)
            minHvDiff = vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt];
          if (((vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt]) > maxHvDiff) &&(vvvRunHv[run][line][dom][pmt]))
            maxHvDiff = vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt];
        }
        vvvHvRange[line][dom][pmt] = std::make_pair(minHvDiff, maxHvDiff);
      }
    }
  }
}

// fill the parameter list from the function fit
inline
void FillParameterFitFunction(vector<double>& vParameter, TGraphErrors*& localGraph, const string& curFunction, const string& curOptions, const pair<double, double> range) {
  TF1* fLocal = new TF1(curFunction.c_str(), curFunction.c_str(), range.first, range.second);
  if (curFunction == "[0] + sqrt([1] - [2]*x)")
    fLocal->SetParameter(1, range.second + 25.0);
  localGraph->Fit(fLocal, curOptions.c_str());
  const int ndf = fLocal->GetNDF();
  const double chi2 = fLocal->GetChisquare();
  double reducedChi2 = 0;
  if (ndf) {
    reducedChi2 = chi2/ndf;
  }
  vParameter.push_back(reducedChi2);
  for (int p = 0; p < fLocal->GetNpar(); ++p) {
    vParameter.push_back(fLocal->GetParameter(p));
  }
}

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

  gStyle->SetLegendBorderSize(0);
  struct timeval startTimeTv;
  gettimeofday(&startTimeTv, NULL);
  time_t startTime = (time_t)startTimeTv.tv_sec;
  tm* startGmTime = gmtime(&startTime);
  const string startUTCTime = asctime(startGmTime);

  JMultipleFileScanner<JDAQTimesliceTypes_t> vInputFiles; 
  int debug;
  string detectorFile;
  string outXMLFile;
  double tot;
  string usr;
  string pwd;
  string mycookies;
  string testType;
  double redChi2Max;
  double slopeMax;
  try {
    JParser<> zap;
    zap['f'] = make_field(vInputFiles);
    zap['a'] = make_field(detectorFile);  
    zap['o'] = make_field(outXMLFile);  
    zap['d'] = make_field(debug) = 1;
    zap['t'] = make_field(tot) = 26.4;
    zap['u'] = make_field(usr) = "";
    zap['!'] = make_field(pwd) = "";
    zap['C'] = make_field(mycookies) = "";
    zap['T'] = make_field(testType) = "HV-TUNING-V2";
    zap['c'] = make_field(redChi2Max) = 10;
    zap['P'] = make_field(slopeMax) = 1;
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  JTriggerParameters parameters;
  try {
    parameters = getTriggerParameters(vInputFiles);
  }
  catch(const JException& error) {
    FATAL("No trigger parameters from input:" << error.what() << endl);
  }
  cout.tie(&cerr);

  JDetector detector;
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  const JModuleRouter moduleRouter(detector);
  const int lineNb = getNumberOfStrings(detector);
  const set<int> lineIds = getStringIDs(detector);
  const int runNb = vInputFiles.size();
  const int pmtNb = NUMBER_OF_PMTS;
  const int domNb = getNumberOfFloors(detector);
  const double diffHVMax = 250;

  string cookiesContent;
  ifstream incookies(mycookies.c_str(), ofstream::in);
  getline(incookies, cookiesContent);

  // stream the vendor hvs from the db
  vector<vector<vector<double> > > vvvHv;
  vvvHv.resize(lineNb, vector<vector<double> > (domNb, vector<double>(pmtNb)));
  vector<vector<vector<vector<double> > > > vvvRunHv;
  vvvRunHv.resize(runNb, vector<vector<vector<double> > > (lineNb, vector<vector<double> > (domNb, vector<double>(pmtNb))));
  vector<vector<vector<string> > > vvvupi;
  vvvupi.resize(lineNb, vector<vector<string> >(domNb, vector<string>(pmtNb)));
  int detid = detector.getID();
  string streamName = "detectors";
  string location = "";
  AccessDB(vvvupi, vvvHv, cookiesContent, usr, pwd, streamName, detid, "-2", lineIds, location);
  streamName = "pmt_hv_settings";
  AccessDB(vvvupi, vvvHv, cookiesContent, usr, pwd, streamName, detid, "-1", lineIds, location);

  // prepare the ToT histo
  vector<vector<vector<vector<TH1D*> > > > vhToT;
  vhToT.resize(runNb, vector<vector<vector<TH1D*> > >(lineNb, vector<vector<TH1D*> >(domNb, vector<TH1D*>(pmtNb))));
  vector<TH1D*> vhToTReducedChi2;
  vhToTReducedChi2.resize(runNb);
  for (int run = 0; run < runNb; ++run) {
    ostringstream plotname;
    plotname << "hToTReducedChi2-run" << run + 1;
    vhToTReducedChi2[run] = new TH1D(plotname.str().c_str(), plotname.str().c_str(), 101, -0.05, 10.05);
    for (int line = 0; line < lineNb; ++line) {
      for (int dom = 0; dom < domNb; ++dom) {
        for (int pmt = 0; pmt < pmtNb; ++pmt) {
          plotname.str("");
          plotname.clear();
          plotname << "hToT-run" << run + 1 << "-line" << line + 1 << "-dom" << dom + 1 << "-pmt" << pmt;
          vhToT[run][line][dom][pmt] = new TH1D(plotname.str().c_str(), plotname.str().c_str(), 301, -0.5, 300.5);
        }
      }
    }
  }

  // read the data and fill the ToT histo
  vector<int> runList;
  int runCounter = 0;
  int oldRun = -1;
  JROOTClassSelector selector("JDAQTimesliceL1");
  JLANG::JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(vInputFiles, selector);
  counter_type counter = 0;
  for ( ; in.hasNext() && counter != vInputFiles.getLimit(); ++counter) {
    JDAQTimeslice* timeslice = in.next();
    const int newRun = timeslice->getRunNumber();
    if (newRun != oldRun) {
      streamName = "pmt_hv_run_settings";
      stringstream runIdSStr;
      runIdSStr << newRun;
      const string newRunStr = runIdSStr.str();
      AccessDB(vvvupi, vvvRunHv[runCounter], cookiesContent, usr, pwd, streamName, detid, newRunStr, lineIds, location);
      ++runCounter;
      runList.push_back(newRun);
    }
    oldRun = newRun;
    for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
      const int domId = super_frame->getModuleID();
      const JModule module = moduleRouter.getModule(domId);
      const int lineIndex = module.getString() - 1; 
      const int domIndex = module.getFloor() - 1;
      for(JDAQSuperFrame::const_iterator hit = super_frame->begin(); hit!=super_frame->end(); ++hit){
        const double hitToT = hit->getToT();
        const int pmtIndex = hit->getPMT();
        const int pmtId = module.getPMT(pmtIndex).getID();
        const string pmtIdDBStr = vvvupi[lineIndex][domIndex][pmtIndex].substr(22);
        const int pmtIdDB = atoi(pmtIdDBStr.c_str());
        if (pmtId != pmtIdDB) {
          DEBUG("mismatch between pmt ids of detector file " << pmtId << " and database "
                << vvvupi[lineIndex][domIndex][pmtIndex] << endl);
          continue;
        }
        vhToT[runCounter - 1][lineIndex][domIndex][pmtIndex]->Fill(hitToT);
      }
    }     
  }

  vector<vector<vector<TGraphErrors*> > > vgToTvsHV;
  vgToTvsHV.resize(lineNb, vector<vector<TGraphErrors*> >(domNb, vector<TGraphErrors*>(pmtNb)));
  for (int line = 0; line < lineNb; ++line) {
    for (int dom = 0; dom < domNb; ++dom) {
      int fullDom = 0;
      for (int pmt = 0; pmt < pmtNb; ++pmt) {
        vgToTvsHV[line][dom][pmt] = new TGraphErrors();
        int localCounter = 0;
        for (int run = 0; run < runNb; ++run) {
          if (vvvRunHv[run][line][dom][pmt] == 0) continue;
          FillGraph(vgToTvsHV[line][dom][pmt], vhToTReducedChi2[run], vhToT[run][line][dom][pmt], vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt], localCounter);
        }  
        if ((localCounter < runNb) && (localCounter > 0))
          DEBUG("only " << localCounter << " runs out of " << runNb << " have ToT data for line " << line + 1 << " - dom " << dom + 1 << " - pmt " << pmt << endl);
        if (localCounter == 0) ++fullDom;
      }
      if (fullDom > 0)
        DEBUG(fullDom << " pmt do not have data in dom " << dom + 1 << " of line " << line + 1 << endl);
    }
  }

  // create a xml file with the final high voltage values.
  fstream outCalibXml(outXMLFile.c_str(), ofstream::out);
  struct timeval endTimeTv;
  gettimeofday(&endTimeTv, NULL);
  time_t endTime = (time_t)endTimeTv.tv_sec;
  tm* endGmTime = gmtime(&endTime);
  const string endUTCTime = asctime(endGmTime);

  ostringstream endTimeStr;
  endTimeStr << 1900 + endGmTime->tm_year << "-" << std::setw(2) << std::setfill('0') 
             << endGmTime->tm_mon + 1 << "-" << std::setw(2) << std::setfill('0')  
             << endGmTime->tm_mday << "T" << std::setw(2) << std::setfill('0') 
             << endGmTime->tm_hour << ":" << std::setw(2) << std::setfill('0') 
             << endGmTime->tm_min << ":" << std::setw(2) << std::setfill('0') 
             << endGmTime->tm_sec << "." << std::setw(6) << std::setfill('0') 
             << endTimeTv.tv_usec << "+01:00";
  ostringstream startTimeStr; 
  startTimeStr << 1900 + startGmTime->tm_year << "-" << std::setw(2) << std::setfill('0')
               << startGmTime->tm_mon + 1 << "-" << std::setw(2) << std::setfill('0')
               << startGmTime->tm_mday << "T" << std::setw(2) << std::setfill('0')
               << startGmTime->tm_hour << ":" << std::setw(2) << std::setfill('0')
               << startGmTime->tm_min << ":" << std::setw(2) << std::setfill('0')
               << startGmTime->tm_sec << "." << std::setw(6) << std::setfill('0')
               << startTimeTv.tv_usec << "+01:00";
  outCalibXml << "<ProductTestSession>\r\n";
  outCalibXml << "\t<User>" << usr << "</User>\r\n";
  outCalibXml << "\t<Location>" << location << "</Location>\r\n";
  outCalibXml << "\t<StartTime>" << startTimeStr.str() << "</StartTime>\r\n";
  outCalibXml << "\t<EndTime>" << endTimeStr.str() << "</EndTime>\r\n";
  outCalibXml << "\t<TestType>" << testType << "</TestType>\r\n";
  outCalibXml << "\t<Tests>\r\n";
  TCanvas* cToTvsHV = new TCanvas("cToTvsHV", "cToTvsHV", 800, 600);
  TH1D* hToTvsHVReducedChi2 = new TH1D("hToTvsHVReducedChi2", "hToTvsHVReducedChi2", 101, -0.05, 10.05);
  TH1D* hToTExpected = new TH1D("hToTExpected", "hToTExpected", 100001, -0.0005, 100.0005);
  bool firstEvent = true;
  vector<vector<vector<double> > > vvvHvOpt;
  vvvHvOpt.resize(lineNb, vector<vector<double> >(domNb, vector<double>(pmtNb)));
  vector<vector<vector<pair<double, double> > > > vvvHvRange;
  vvvHvRange.resize(lineNb, vector<vector<pair<double, double> > >(domNb, vector<pair<double, double> >(pmtNb)));
  FindHVRange(vvvHvRange, vvvRunHv, vvvHv);
  int counterLinear = 0;
  int counter2Order = 0;
  int counter2OrderPos = 0;
  int counter2OrderNeg = 0;
  int counterSqrt = 0;
  double aveLinearSlope = 0;
  for (int line = 0; line < lineNb; ++line) {
    for (int dom = 0; dom < domNb; ++dom) {
      for (int pmt = 0; pmt < pmtNb; ++pmt) {
        stringstream filename;
        filename.str("");
        filename.clear();
        filename << "ToTvsHV-line" << line + 1 << "-dom" << dom + 1 << "-pmt" << pmt;
        vgToTvsHV[line][dom][pmt]->SetName(filename.str().c_str());
        const int pointNb = vgToTvsHV[line][dom][pmt]->GetN();
        //if ((line == 0) && (dom == 17) && (pmt == 24)) firstEvent = true;
        if (pointNb > 0) {
          vgToTvsHV[line][dom][pmt]->SetLineColor(line*18 + dom*31 + pmt);
          vgToTvsHV[line][dom][pmt]->SetMarkerColor(line*18 + dom*31 + pmt);
          vgToTvsHV[line][dom][pmt]->SetMarkerStyle(2);
          vgToTvsHV[line][dom][pmt]->SetMarkerSize(1);
          if (firstEvent == true) {
            firstEvent = false;
            vgToTvsHV[line][dom][pmt]->SetTitle("");
            vgToTvsHV[line][dom][pmt]->GetXaxis()->SetTitle("HV_{offset}  [V]");
            vgToTvsHV[line][dom][pmt]->GetXaxis()->CenterTitle();
            vgToTvsHV[line][dom][pmt]->GetXaxis()->SetTitleOffset(1.2);
            vgToTvsHV[line][dom][pmt]->GetXaxis()->SetRangeUser(-200, +200);
            vgToTvsHV[line][dom][pmt]->GetYaxis()->SetTitle("ToT values  [ns]");
            vgToTvsHV[line][dom][pmt]->GetYaxis()->CenterTitle();
            vgToTvsHV[line][dom][pmt]->GetYaxis()->SetTitleOffset(1.2);
            vgToTvsHV[line][dom][pmt]->GetYaxis()->SetRangeUser(0, 80);
            vgToTvsHV[line][dom][pmt]->Draw("APL");
          } else {
            vgToTvsHV[line][dom][pmt]->GetXaxis()->SetRangeUser(-200, +200);
            vgToTvsHV[line][dom][pmt]->Draw("samePL");
          }
        }
        pair<double, double> hvAndToTOpt = make_pair(0, 0);
        string functionFormula = "";
        if (pointNb > 2) {
          vector<double> vParameterPol2;
          const string pol2Name = "pol2";
          FillParameterFitFunction(vParameterPol2, vgToTvsHV[line][dom][pmt], pol2Name, "0RQ", vvvHvRange[line][dom][pmt]); 
          vector<double> vParameterPol1;
          const string pol1Name = "pol1";
          FillParameterFitFunction(vParameterPol1, vgToTvsHV[line][dom][pmt], pol1Name, "0RQ", vvvHvRange[line][dom][pmt]);
          vector<double> vParameterSqrt;
          const string sqrtName = "[0] + sqrt([1] - [2]*x)";
          FillParameterFitFunction(vParameterSqrt, vgToTvsHV[line][dom][pmt], sqrtName, "0RQ", vvvHvRange[line][dom][pmt]);
          vector<double> vParameter;
          const double testChi2 = min({fabs(vParameterPol1[0] - 1), fabs(vParameterPol2[0] - 1), fabs(vParameterSqrt[0] - 1)});
          if (testChi2 == fabs(vParameterPol1[0] - 1)) {
            vParameter = vParameterPol1;
            functionFormula = pol1Name;
            aveLinearSlope += vParameterPol1[2];
            ++counterLinear;
          } else if (testChi2 == fabs(vParameterPol2[0] - 1)) {
            vParameter = vParameterPol2;
            functionFormula = pol2Name;
            ++counter2Order;
            if (vParameterPol2[3] > 0) ++counter2OrderPos;
            else ++counter2OrderNeg;
          } else if (testChi2 == fabs(vParameterSqrt[0] - 1)) {
            vParameter = vParameterSqrt;
            functionFormula = sqrtName;
            ++counterSqrt;
          } else {
            DEBUG("problem in the function definition, check the code at about line 500");
          }
          vgToTvsHV[line][dom][pmt]->Fit(functionFormula.c_str(), "RQ");
          hToTvsHVReducedChi2->Fill(vParameter[0]);
          if (vParameter[0] > redChi2Max) DEBUG("ToTvsHV fit may be bad for pmt " << pmt << " of dom " << dom + 1 << " of line " << line + 1 << ", the reduced chi square is " << vParameter[0] << endl);
          hvAndToTOpt = Solve(vParameter, functionFormula, tot);
          if ((hvAndToTOpt.first == -1) && (hvAndToTOpt.second == -1))
            DEBUG("function used for fitting " << functionFormula << " is not handled by the code. Change the code or use a predefined function such as pol1 or pol2" << endl);
          if ((fabs(vParameter[2]) > slopeMax) && ((functionFormula == "pol1") || (functionFormula == "pol2")))
            DEBUG("ToTvsHV fit may be bad for pmt " << pmt << " of dom " << dom + 1 << " of line " << line + 1 << ", the slope at 0 is " << vParameter[2] << endl);
          hToTExpected->Fill(hvAndToTOpt.second);
        }
        if (fabs(hvAndToTOpt.first) > diffHVMax) hvAndToTOpt.first = 0;
        vvvHvOpt[line][dom][pmt] = hvAndToTOpt.first;
        outCalibXml << "\t\t<ProductTest>\r\n";
        outCalibXml << "\t\t\t<UPI>" << vvvupi[line][dom][pmt] << "</UPI>\r\n";
        outCalibXml << "\t\t\t<TestResult>OK</TestResult>\r\n";
        outCalibXml << "\t\t\t<TestParameters>\r\n";
        outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
        outCalibXml << "\t\t\t\t\t<Name>PMT_Supply_Voltage</Name>\r\n";
        outCalibXml << "\t\t\t\t\t<Values>\r\n";
        outCalibXml << "\t\t\t\t\t\t<string>" << hvAndToTOpt.first + vvvHv[line][dom][pmt] << "</string>\r\n";
        outCalibXml << "\t\t\t\t\t</Values>\r\n";
        outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
        outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
        outCalibXml << "\t\t\t\t\t<Name>RUN_NUMBER</Name>\r\n";
        outCalibXml << "\t\t\t\t\t<Values>\r\n";
        for (int run = 0; run < runNb; ++run) {
          outCalibXml << "\t\t\t\t\t\t<string>" << runList[run] << "</string>\r\n";
        }
        outCalibXml << "\t\t\t\t\t</Values>\r\n";
        outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
        outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
        outCalibXml << "\t\t\t\t\t<Name>PMT_Time_over_Threshold</Name>\r\n";
        outCalibXml << "\t\t\t\t\t<Values>\r\n";
        outCalibXml << "\t\t\t\t\t\t<string>" << tot*1e-9 << "</string>\r\n";
        outCalibXml << "\t\t\t\t\t</Values>\r\n";
        outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
        outCalibXml << "\t\t\t</TestParameters>\r\n";
        outCalibXml << "\t\t</ProductTest>\r\n";
      }
    }
  }
  cout << "linear:" << counterLinear << " - 2nd order (pos, neg):" << counter2Order 
       << "(" << counter2OrderPos << "," << counter2OrderNeg << ")"
       << " - sqrt:" << counterSqrt << " - linear slope average: " 
       << aveLinearSlope/counterLinear << endl;
  outCalibXml << "\t</Tests>\r\n";
  outCalibXml << "</ProductTestSession>\r\n";
  outCalibXml.close();
  cToTvsHV->SaveAs("figs/cToTvsHV.png");
  cToTvsHV->SaveAs("figs/cToTvsHV.root");

  ostringstream osstxt;
  osstxt.clear();
  osstxt.str("");
  osstxt << "results/HVTuning-det" << detid << ".txt";
  
  ofstream outHv(osstxt.str().c_str(), ofstream::out);
  for (int line = 0; line < lineNb; ++line) {
    for (int dom = 0; dom < domNb; ++dom) {
      for (int pmt = 0; pmt < pmtNb; ++pmt) {
        outHv << line + 1 << ' ' << dom + 1 << ' ' << pmt << ' ' << vvvupi[line][dom][pmt] << ' ' << vvvHv[line][dom][pmt] << ' ' << vvvHvOpt[line][dom][pmt] + vvvHv[line][dom][pmt] << endl;
      }
    }
  }
  outHv.close();

  TCanvas* cToTs[lineNb][domNb];
  for (int line = 0; line < lineNb; ++line) {
    for (int dom = 0; dom < domNb; ++dom) {
      ostringstream iss;
      iss.clear();
      iss.str("");
      iss << "cToTs-line" << line + 1 << "-dom" << dom + 1; 
      cToTs[line][dom] = new TCanvas(iss.str().c_str(), iss.str().c_str(), 800, 600);
      cToTs[line][dom]->Divide(8,4);
      for (int pmt = 0; pmt < pmtNb; ++pmt) {
        cToTs[line][dom]->cd(pmt + 1);
        double localMax = 0;
        for (int run = 0; run < runNb; ++run) { 
          vhToT[run][line][dom][pmt]->SetLineColor(run + 1);
          if (vhToT[run][line][dom][pmt]->GetMaximum() > localMax) localMax = vhToT[run][line][dom][pmt]->GetMaximum();
          if (run == 0) {
            vhToT[run][line][dom][pmt]->SetTitle("");
            vhToT[run][line][dom][pmt]->SetStats(0);
            vhToT[run][line][dom][pmt]->GetXaxis()->SetTitle("ToT values  [ns]");
            vhToT[run][line][dom][pmt]->GetXaxis()->CenterTitle();
            vhToT[run][line][dom][pmt]->GetXaxis()->SetTitleOffset(1.2);
            vhToT[run][line][dom][pmt]->GetXaxis()->SetRange(0, 100);
            vhToT[run][line][dom][pmt]->Draw();
          } else {
            vhToT[run][line][dom][pmt]->Draw("same");      
          }
          gPad->SetLogy();
        }
        vhToT[0][line][dom][pmt]->GetYaxis()->SetRangeUser(0.1, 2*localMax);
      }
      iss.clear();
      iss.str("");
      iss << "figs/cToTs-line" << line + 1 << "-dom" << dom + 1 << ".png"; 
      cToTs[line][dom]->SaveAs(iss.str().c_str());
      iss.clear();
      iss.str("");
      iss << "figs/cToTs-line" << line + 1 << "-dom" << dom + 1 << ".root";  
      cToTs[line][dom]->SaveAs(iss.str().c_str());
    }
  }

  TCanvas* cToTReducedChi2 = new TCanvas("cToTReducedChi2", "cToTReducedChi2", 800, 600);
  for (int run = 0; run < runNb; ++run) {
    if (run == 0) { 
      vhToTReducedChi2[run]->SetTitle("");
      vhToTReducedChi2[run]->SetStats(0);
      vhToTReducedChi2[run]->GetXaxis()->SetTitle("reduced #chi^{2}");
      vhToTReducedChi2[run]->GetXaxis()->CenterTitle();
      vhToTReducedChi2[run]->GetXaxis()->SetTitleOffset(1.2);
      vhToTReducedChi2[run]->Draw();
    } else {
      vhToTReducedChi2[run]->Draw("same");
    }
    vhToTReducedChi2[run]->SetLineColor(run + 1);
  }
  cToTReducedChi2->SetLogy();
  cToTReducedChi2->SaveAs("figs/cToT-ReducedChi2.png");
  cToTReducedChi2->SaveAs("figs/cToT-ReducedChi2.root");

  TCanvas* cToTvsHVReducedChi2 = new TCanvas("cToTvsHVReducedChi2", "cToTvsHVReducedChi2", 800, 600);
  hToTvsHVReducedChi2->SetTitle("");
  hToTvsHVReducedChi2->SetStats(0);
  hToTvsHVReducedChi2->GetXaxis()->SetTitle("reduced #chi^{2}");
  hToTvsHVReducedChi2->GetXaxis()->CenterTitle();
  hToTvsHVReducedChi2->GetXaxis()->SetTitleOffset(1.2);
  hToTvsHVReducedChi2->Draw();
  cToTvsHVReducedChi2->SetLogy();
  cToTvsHVReducedChi2->SaveAs("figs/cToTvsHV-ReducedChi2.png");
  cToTvsHVReducedChi2->SaveAs("figs/cToTvsHV-ReducedChi2.root");

  TCanvas* cToTExpected = new TCanvas("cToTExpected", "cToTExpected", 800, 600);
  hToTExpected->SetTitle("");
  hToTExpected->SetStats(0);
  hToTExpected->GetXaxis()->SetTitle("Expected ToT value  [ns]");
  hToTExpected->GetXaxis()->CenterTitle();
  hToTExpected->GetXaxis()->SetTitleOffset(1.2);
  hToTExpected->Draw();
  //cToTExpected->SetLogy();
  cToTExpected->SaveAs("figs/cToTExpected.png");
  cToTExpected->SaveAs("figs/cToTExpected.root");
  
  const int lineId1 = 0;
  const int domId1 = 2;
  const int pmtId1 = 5;
  const int lineId2 = 0;
  const int domId2 = 5;
  const int pmtId2 = 7;
  ostringstream issTemp;
  issTemp << "cToTsdom" << domId1 + 1 << "pmt" << pmtId1;
  TCanvas* cToTsdomxpmtx = new TCanvas(issTemp.str().c_str(), issTemp.str().c_str(), 800, 600);
  bool pmtTest = true;
  TLegend* leg1 = new TLegend(0.70,0.55,0.90,0.90);
  leg1->SetFillStyle(0);
  double localMaxTemp = 0;
  if ((lineId1 > lineNb - 1) || (domId1 > domNb - 1) || (pmtId1 > pmtNb - 1))
    pmtTest = false;
  if (pmtTest == false) {
    DEBUG("wrong PMT1 identification for sample plot");
  } else {
    ostringstream plotname0;
    for (int run = 0; run < runNb; ++run) {
      vhToT[run][lineId1][domId1][pmtId1]->SetLineColor(run + 1);
      if (vhToT[run][lineId1][domId1][pmtId1]->GetMaximum() > localMaxTemp) localMaxTemp = vhToT[run][lineId1][domId1][pmtId1]->GetMaximum();
      if (run == 0) {
        vhToT[run][lineId1][domId1][pmtId1]->SetTitle("");
        vhToT[run][lineId1][domId1][pmtId1]->SetStats(0);
        vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetTitle("ToT values  [ns]");
        vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->CenterTitle();
        vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetTitleOffset(1.2);
        vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetRange(0, 100);
        vhToT[run][lineId1][domId1][pmtId1]->Draw();
      } else {
        vhToT[run][lineId1][domId1][pmtId1]->Draw("same");
      }
      plotname0.str("");
      plotname0.clear();
      plotname0 << "HV: " << setprecision(2) << fixed
                << vvvRunHv[run][lineId1][domId1][pmtId1] << "V";
      leg1->AddEntry(vhToT[run][lineId1][domId1][pmtId1], plotname0.str().c_str(), "l");
      gPad->SetLogy();
    }
    vhToT[0][lineId1][domId1][pmtId1]->GetYaxis()->SetRangeUser(0.1, 2*localMaxTemp);
    leg1->Draw();
    issTemp.clear();
    issTemp.str("");
    issTemp << "figs/cToTs-line" << lineId1 + 1 << "-dom" << domId1 + 1 << "-pmt" << pmtId1 << ".png";
    cToTsdomxpmtx->SaveAs(issTemp.str().c_str());
    issTemp.clear();
    issTemp.str("");
    issTemp << "figs/cToTs-line" << lineId1 + 1 << "-dom" << domId1 + 1 << "-pmt" << pmtId1 << ".root";
    cToTsdomxpmtx->SaveAs(issTemp.str().c_str());
  }

  issTemp.clear();
  issTemp.str("");
  issTemp << "cToTsdom" << domId2 + 1 << "pmt" << pmtId2;
  TCanvas* cToTsdomxpmty = new TCanvas(issTemp.str().c_str(), issTemp.str().c_str(), 800, 600);
  pmtTest = true;
  TLegend* leg2 = new TLegend(0.70,0.55,0.90,0.90);
  if ((lineId2 > lineNb - 1) || (domId2 > domNb - 1) || (pmtId2 > pmtNb - 1))
    pmtTest = false;
  if (pmtTest == false) {
    DEBUG("wrong PMT2 identification for sample plot");
  } else {
    ostringstream plotname0;
    for (int run = 0; run < runNb; ++run) {
      vhToT[run][lineId2][domId2][pmtId2]->SetLineColor(run + 1);
      if (vhToT[run][lineId2][domId2][pmtId2]->GetMaximum() > localMaxTemp) localMaxTemp = vhToT[run][lineId2][domId2][pmtId2]->GetMaximum();
      if (run == 0) {
        vhToT[run][lineId2][domId2][pmtId2]->SetTitle("");
        vhToT[run][lineId2][domId2][pmtId2]->SetStats(0);
        vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetTitle("ToT values  [ns]");
        vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->CenterTitle();
        vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetTitleOffset(1.2);
        vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetRange(0, 100);
        vhToT[run][lineId2][domId2][pmtId2]->Draw();
      } else {
        vhToT[run][lineId2][domId2][pmtId2]->Draw("same");
      }
      plotname0.str("");
      plotname0.clear();
      plotname0 << "HV: " << setprecision(2) << fixed
                << vvvRunHv[run][lineId2][domId2][pmtId2] << "V";
      leg2->AddEntry(vhToT[run][lineId2][domId2][pmtId2], plotname0.str().c_str(), "l");
      gPad->SetLogy();
    }
    vhToT[0][lineId2][domId2][pmtId2]->GetYaxis()->SetRangeUser(0.1, 2*localMaxTemp);
    leg2->Draw();
    issTemp.clear();
    issTemp.str("");
    issTemp << "figs/cToTs-line" << lineId2 + 1 << "-dom" << domId2 + 1 << "-pmt" << pmtId2 << ".png";
    cToTsdomxpmty->SaveAs(issTemp.str().c_str());
    issTemp.clear();
    issTemp.str("");
    issTemp << "figs/cToTs-line" << lineId2 + 1 << "-dom" << domId2 + 1 << "-pmt" << pmtId2 << ".root";
    cToTsdomxpmty->SaveAs(issTemp.str().c_str());
  }

  DEBUG("analysis duration: " << difftime(endTime, startTime) << endl);
}