Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JHVTuning.cc
Go to the documentation of this file.
1 #include <string>
2 #include <vector>
3 #include <ctime>
4 #include <sys/time.h>
5 #include <unistd.h>
6 
7 #include "JDAQ/JDAQTimeslice.hh"
8 
9 #include "JSupport/JMultipleFileScanner.hh"
10 #include "JSupport/JSupport.hh"
11 #include "Jeep/JParser.hh"
12 #include "Jeep/JMessage.hh"
13 #include "Jeep/JTimer.hh"
14 #include "JLang/JString.hh"
15 #include "JLang/JObjectMultiplexer.hh"
16 #include "JROOT/JROOTClassSelector.hh"
17 
18 #include "JDetector/JDetector.hh"
19 #include "JDetector/JDetectorToolkit.hh"
20 
21 #include "JDetector/JModuleRouter.hh"
22 
23 #include "TH1D.h"
24 #include "TGraphErrors.h"
25 #include "TF1.h"
26 #include "TAxis.h"
27 #include "TCanvas.h"
28 #include <TStyle.h>
29 #include <TLegend.h>
30 #include <TLegendEntry.h>
31 
32 #include "KM3NeTDBClient.h"
33 
34 
35 /**
36  * \file
37  * Program for production of the xml file with nominal high voltage
38  * values. The calculation is based on the adjustment of the ToT
39  * duration to a value of 26.4ns by scanning the high voltage. Inputs
40  * are the root datafile list for each of the high voltage values
41  * (option -f), the detector file (option -a), the username and password
42  * (option -u) for access to the database.
43  * the detector file can be retrieved in the database:
44  * https://km3netdbweb.in2p3.fr/detx/D_ARCA003?
45  *
46  * \author acreusot
47  */
48 
49 // retrieve information in the database
50 using namespace KM3NeT;
51 using namespace KM3NeT::DB;
52 inline
53 //void AccessDB(std::vector<std::vector<std::vector<int> > >& vvvupi, std::vector<std::vector<std::vector<int> > >& vvvupiExt, std::vector<std::vector<std::vector<double> > >& vvvHv, std::string& username, std::string& password, std::string& streamName, std::string& detectorId, const int runId, const std::string& hvDefault) {
54 void AccessDB(std::vector<std::vector<std::vector<int> > >& vvvupi, std::vector<std::vector<std::vector<int> > >& vvvupiExt, std::vector<std::vector<std::vector<double> > >& vvvHv, std::string& cookie, std::string& username, std::string& streamName, std::string& detectorId, const int runId, const std::string& hvDefault, std::vector<int>& lineIds) {
55  bool validRun = false;
56  int tryCounter = 0;
57  while ((validRun == false) && (tryCounter < 10)) {
58  try {
59  std::vector<Selector> vSel;
60  vSel.clear();
61  if (streamName == "vendorhv") {
62  Selector detIdSel("detid", detectorId.c_str());
63  vSel.push_back(detIdSel);
64  }
65  if (streamName == "pmt_hv_settings") {
66  Selector detIdSel("detid", detectorId.c_str());
67  vSel.push_back(detIdSel);
68  }
69  if (streamName == "pmt_hv_run_settings") {
70  std::stringstream runIdSStr;
71  runIdSStr << runId;
72  std::string runIdStr = runIdSStr.str();
73  Selector runIdSel("run", runIdStr.c_str());
74  vSel.push_back(runIdSel);
75  }
76  //auto client = Client::Create(Server::Default, username.c_str(), password.c_str());
77  auto client = Client::Create(Server::Default, cookie.c_str());
78  username = client->User();
79  if (runId == -1) std::cout << "opening connection as user " << username << " for default HV values" << std::endl;
80  else std::cout << "opening connection as user " << username << " for run " << runId << std::endl;
81  auto &rs = client->StreamDS(streamName.c_str(), vSel);
82  while (rs.Next()) {
83  int lineId = -1;
84  int domId = -1;
85  int pmtId = -1;
86  int upi = -1;
87  int upiExt = -1;
88  int hv = 0;
89  bool test = true;
90  bool validHv = false;
91  for (unsigned int i = 0; i < rs.FieldCount(); ++i) {
92  std::istringstream valueSStr(rs.GetString(i));
93  int value;
94  valueSStr >> value;
95  std::string fieldName = rs.FieldName(i);
96  if (fieldName == "DUID") {
97  for (unsigned int i = 0; i < lineIds.size(); ++i) {
98  if (lineIds[i] == value) {
99  lineId = i;
100  }
101  }
102  } else if (fieldName == "FLOORID") {
103  domId = value - 1;
104  } else if (fieldName == "CABLEPOS") {
105  pmtId = value;
106  } else if (fieldName == "PMTSERIAL") {
107  upi = value;
108  } else if (fieldName == "PMTUPI") {
109  std::string fullUpi = valueSStr.str().substr(22, 5);
110  std::istringstream fullUpiSStr(fullUpi);
111  fullUpiSStr >> upi;
112  std::string fullUpiExt = valueSStr.str().substr(20, 2);
113  std::istringstream fullUpiExtSStr(fullUpiExt);
114  fullUpiExtSStr >> upiExt;
115  } else if (fieldName == "CG_NAME") {
116  if (valueSStr.str() == hvDefault) validHv = true;
117  } else if (fieldName == "HV_VALUE") {
118  hv = value;
119  if (streamName == "pmt_hv_run_settings") validHv = true;
120  } else if (fieldName == "PMT_SUPPLY_VOLTAGE") {
121  validHv = true;
122  hv = value;
123  } else if (fieldName == "DETID") {
124  if (valueSStr.str() != detectorId) test = false;
125  } else if (fieldName == "PMTEXT") {
126  upiExt = value;
127  }
128  }
129  //std::cout << lineId << ' ' << domId << ' ' << pmtId << ' ' << upi << ' ' << hv << std::endl;
130  if (test == true) {
131  if (validHv == true) vvvHv[lineId][domId][pmtId] = hv;
132  vvvupi[lineId][domId][pmtId] = upi;
133  if (streamName == "pmt_hv_settings") vvvupiExt[lineId][domId][pmtId] = upiExt;
134  }
135  }
136  rs.Close();
137  client->Close();
138  validRun = true;
139  }
140  catch (DBException x) {
141  std::cout << x.what() << std::endl;
142  ++tryCounter;
143  }
144  }
145 }
146 
147 // solve the 2nd order polynom. If delta is negatif use a 1st order polynom.
148 inline
149 double SolveAnalytic(const double deg0, const double deg1, const double deg2, const double tot) {
150  const double delta = deg1*deg1 - 4*deg2*(deg0 - tot);
151  if ((delta < 0) || (deg2 == 0)) return deg1? (tot - deg0)/deg1 : 0;
152  const double x1 = (-deg1 + sqrt(delta))/(2*deg2);
153  const double x2 = (-deg1 - sqrt(delta))/(2*deg2);
154  if (deg2 > 0) return (x1 < x2)? x1 : x2;
155  else return (x1 < x2)? x2 : x1;
156  //return (fabs(x1) < fabs(x2))? x1 : x2;
157 }
158 
159 // solve the equation, whatever the function
160 inline
161 std::pair<double, double> Solve(std::vector<double>& vParameter, const std::string& curFunction, const double tot) {
162  std::pair<double, double> curSol;
163  if (curFunction == "pol1") {
164  curSol.first = SolveAnalytic(vParameter[1], vParameter[2], 0, tot);
165  curSol.second = vParameter[2]*curSol.first + vParameter[1];
166  } else if (curFunction == "pol2") {
167  curSol.first = SolveAnalytic(vParameter[1], vParameter[2], vParameter[3], tot);
168  curSol.second = vParameter[3]*pow(curSol.first, 2) + vParameter[2]*curSol.first + vParameter[1];
169  } else if (curFunction == "[0] + sqrt([1] - [2]*x)") {
170  curSol.first = (vParameter[2] - pow(tot - vParameter[1], 2))/vParameter[3];
171  curSol.second = vParameter[1] + sqrt(vParameter[2] - vParameter[3]*curSol.first);
172  } else {
173  curSol.first = -1;
174  curSol.second = -1;
175  }
176  return curSol;
177 }
178 
179 // optimize reduced chi2 and fill the ToT vs HV graph
180 inline
181 void FillGraph(TGraphErrors*& graph, TH1D*& hChi2, TH1D*& histo, const double hv, int& localCounter) {
182  const int firstBin = 15;
183  const int lastBin = histo->GetXaxis()->GetLast() - 60;
184  histo->GetXaxis()->SetRange(firstBin, lastBin);
185  if (!histo->GetEntries()) return;
186  if (!histo->Integral()) return;
187  const int maxBin = histo->GetMaximumBin();
188  double bestRedChi2 = 1e6;
189  double bestMean = 0;
190  for (int binInf = maxBin - 6; binInf < maxBin; ++binInf) {
191  for (int binSup = maxBin + 1; binSup <= maxBin + 6; ++binSup) {
192  if (binSup - binInf < 4) continue;
193  TF1* gausFit = new TF1("gausFit", "gaus", binInf, binSup);
194  const double curInteg = histo->Integral(binInf, binSup);
195  if (curInteg == 0) continue;
196  histo->Fit(gausFit, "RQ");
197  if (!gausFit->GetNDF()) {
198  gausFit->Delete();
199  continue;
200  }
201  const double fitMean = gausFit->GetParameter(1);
202  const double fitSigma = gausFit->GetParameter(2);
203  const double reducedChi2 = gausFit->GetChisquare()/gausFit->GetNDF();
204  if ((std::fabs(reducedChi2 - 1) < std::fabs(bestRedChi2 - 1)) && (reducedChi2 > 0.5) && (fitSigma < 10)) {
205  bestRedChi2 = reducedChi2;
206  bestMean = fitMean;
207  }
208  gausFit->Delete();
209  }
210  }
211  hChi2->Fill(bestRedChi2);
212  histo->GetXaxis()->SetRange(0, lastBin);
213  if (bestMean > 0) {
214  graph->SetPoint(localCounter, hv, bestMean);
215  ++localCounter;
216  }
217 }
218 
219 // find the HV range used to perform the scan
220 inline
221 void FindHVRange(std::vector<std::vector<std::vector<std::pair<double, double> > > >& vvvHvRange, std::vector<std::vector<std::vector<std::vector<double> > > >& vvvRunHv, std::vector<std::vector<std::vector<double> > >& vvvHv) {
222  for (unsigned int line = 0; line < vvvHvRange.size(); ++line) {
223  for (unsigned int dom = 0; dom < vvvHvRange[line].size(); ++dom) {
224  for (unsigned int pmt = 0; pmt < vvvHvRange[line][dom].size(); ++pmt) {
225  double minHvDiff = 1e8;
226  double maxHvDiff = 0;
227  for (unsigned int run = 0; run < vvvRunHv.size(); ++run) {
228  if ((vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt]) < minHvDiff)
229  minHvDiff = vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt];
230  if (((vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt]) > maxHvDiff) &&(vvvRunHv[run][line][dom][pmt]))
231  maxHvDiff = vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt];
232  }
233  vvvHvRange[line][dom][pmt] = std::make_pair(minHvDiff, maxHvDiff);
234  }
235  }
236  }
237 }
238 
239 // fill the parameter list from the function fit
240 inline
241 void FillParameterFitFunction(std::vector<double>& vParameter, TGraphErrors*& localGraph, const std::string& curFunction, const std::string& curOptions, const std::pair<double, double> range) {
242  TF1* fLocal = new TF1(curFunction.c_str(), curFunction.c_str(), range.first, range.second);
243  if (curFunction == "[0] + sqrt([1] - [2]*x)")
244  fLocal->SetParameter(1, range.second + 25.0);
245  localGraph->Fit(fLocal, curOptions.c_str());
246  const int ndf = fLocal->GetNDF();
247  const double chi2 = fLocal->GetChisquare();
248  double reducedChi2 = 0;
249  if (ndf) {
250  reducedChi2 = chi2/ndf;
251  }
252  vParameter.push_back(reducedChi2);
253  for (int p = 0; p < fLocal->GetNpar(); ++p) {
254  vParameter.push_back(fLocal->GetParameter(p));
255  }
256 }
257 
258 int main(int argc, char **argv)
259 {
260  using namespace std;
261  using namespace JPP;
262  using namespace KM3NETDAQ;
263 
264  gStyle->SetLegendBorderSize(0);
265  struct timeval startTimeTv;
266  gettimeofday(&startTimeTv, NULL);
267  time_t startTime = (time_t)startTimeTv.tv_sec;
268  tm* startGmTime = gmtime(&startTime);
269  const string startUTCTime = asctime(startGmTime);
270 
271  JMultipleFileScanner<JDAQTimesliceTypes_t> vInputFiles;
272  int debug;
273  std::string detectorFile;
274  std::string outXMLFile;
275  double tot;
276  //std::string usernameAndPassword;
277  std::string mycookies;
278  std::string testType;
279  std::string location;
280  double redChi2Max;
281  double slopeMax;
282  std::string hvDefault;
283  try {
284  JParser<> zap;
285  zap['f'] = make_field(vInputFiles);
286  zap['a'] = make_field(detectorFile);
287  zap['o'] = make_field(outXMLFile);
288  zap['d'] = make_field(debug) = 1;
289  zap['t'] = make_field(tot) = 26.4;
290  //zap['u'] = make_field(usernameAndPassword);
291  zap['u'] = make_field(mycookies);
292  zap['T'] = make_field(testType) = "HV-TUNING-V2";
293  zap['w'] = make_field(location) = "A00070004";
294  zap['C'] = make_field(redChi2Max) = 10;
295  zap['P'] = make_field(slopeMax) = 1;
296  zap['D'] = make_field(hvDefault) = "VENDOR";
297  zap(argc, argv);
298  }
299  catch(const exception& error) {
300  FATAL(error.what() << endl);
301  }
302 
303  JDetector detector;
304  try {
305  load(detectorFile, detector);
306  }
307  catch(const JException& error) {
308  FATAL(error);
309  }
310  const JModuleRouter moduleRouter(detector);
311 
312  vector<int> lineIds;
313  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
314  const int curId = module->getString();
315  if (lineIds.size() == 0) {
316  lineIds.push_back(curId);
317  } else {
318  bool newLine = true;
319  for (unsigned int l = 0; l < lineIds.size(); ++l) {
320  if (curId == lineIds[l]) {
321  newLine = false;
322  }
323  }
324  if (newLine == true) {
325  lineIds.push_back(curId);
326  }
327  }
328  }
329  const unsigned int lineNb = lineIds.size();
330 
331  const int AllDomNb = detector.size();
332  const int runNb = vInputFiles.size();
333  const int pmtNb = NUMBER_OF_PMTS;
334  const int domNb = AllDomNb/lineNb;
335  const double diffHVMax = 250;
336 
337  //stringstream ssUser;
338  //ssUser.str(usernameAndPassword);
339  //string sTemp;
340  //vector<string> username;
341  //while (getline(ssUser, sTemp, ':'))
342  // username.push_back(sTemp);
343 
344  string cookiesContent;
345  ifstream incookies(mycookies.c_str(), ofstream::in);
346  getline(incookies, cookiesContent);
347  string username;
348 
349  // stream the vendor hvs from the db
350  vector<vector<vector<vector<double> > > > vvvRunHv;
351  vvvRunHv.resize(runNb, vector<vector<vector<double> > > (lineNb, vector<vector<double> > (domNb, vector<double>(pmtNb))));
352  vector<vector<vector<double> > > vvvHv;
353  vvvHv.resize(lineNb, vector<vector<double> >(domNb, vector<double>(pmtNb)));
354  vector<vector<vector<int> > > vvvupi;
355  vvvupi.resize(lineNb, vector<vector<int> >(domNb, vector<int>(pmtNb)));
356  vector<vector<vector<int> > > vvvupiExt;
357  vvvupiExt.resize(lineNb, vector<vector<int> >(domNb, vector<int>(pmtNb)));
358  stringstream detectorIdSStr;
359  int detectorId = detector.getID();
360  detectorIdSStr << detectorId;
361  string detectorIdStr = detectorIdSStr.str();
362  //string streamName = "vendorhv";
363  string streamName = "pmt_hv_settings";
364  //AccessDB(vvvupi, vvvupiExt, vvvHv, username[0], username[1], streamName, detectorIdStr, -1, hvDefault, lineIds);
365  AccessDB(vvvupi, vvvupiExt, vvvHv, cookiesContent, username, streamName, detectorIdStr, -1, hvDefault, lineIds);
366 
367  // prepare the ToT histo
368  vector<vector<vector<vector<TH1D*> > > > vhToT;
369  vhToT.resize(runNb, vector<vector<vector<TH1D*> > >(lineNb, vector<vector<TH1D*> >(domNb, vector<TH1D*>(pmtNb))));
370  vector<TH1D*> vhToTReducedChi2;
371  vhToTReducedChi2.resize(runNb);
372  for (int run = 0; run < runNb; ++run) {
373  ostringstream plotname;
374  plotname << "hToTReducedChi2-run" << run + 1;
375  vhToTReducedChi2[run] = new TH1D(plotname.str().c_str(), plotname.str().c_str(), 101, -0.05, 10.05);
376  for (unsigned int line = 0; line < lineNb; ++line) {
377  for (int dom = 0; dom < domNb; ++dom) {
378  for (int pmt = 0; pmt < pmtNb; ++pmt) {
379  plotname.str("");
380  plotname.clear();
381  plotname << "hToT-run" << run + 1 << "-line" << line + 1 << "-dom" << dom + 1 << "-pmt" << pmt;
382  vhToT[run][line][dom][pmt] = new TH1D(plotname.str().c_str(), plotname.str().c_str(), 301, -0.5, 300.5);
383  }
384  }
385  }
386  }
387 
388  // read the data and fill the ToT histo
389  vector<int> runList;
390  int runCounter = 0;
391  int oldRun = -1;
392  JROOTClassSelector selector("JDAQTimesliceL1");
393  JLANG::JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(vInputFiles, selector);
394  counter_type counter = 0;
395  for ( ; in.hasNext() && counter != vInputFiles.getLimit(); ++counter) {
396  JDAQTimeslice* timeslice = in.next();
397  const int newRun = timeslice->getRunNumber();
398  if (newRun != oldRun) {
399  streamName = "pmt_hv_run_settings";
400  //AccessDB(vvvupi, vvvupiExt, vvvRunHv[runCounter], username[0], username[1], streamName, detectorIdStr, newRun, hvDefault, lineIds);
401  AccessDB(vvvupi, vvvupiExt, vvvRunHv[runCounter], cookiesContent, username, streamName, detectorIdStr, newRun, hvDefault, lineIds);
402  ++runCounter;
403  runList.push_back(newRun);
404  }
405  oldRun = newRun;
406  for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
407  const int domId = super_frame->getModuleID();
408  const JModule module = moduleRouter.getModule(domId);
409  const int lineIndex = module.getString() - 1;
410  const int domIndex = module.getFloor() - 1;
411  for(JDAQSuperFrame::const_iterator hit = super_frame->begin(); hit!=super_frame->end(); ++hit){
412  const double hitToT = hit->getToT();
413  const int pmtIndex = hit->getPMT();
414  const int pmtId = module.getPMT(pmtIndex).getID();
415  if (pmtId != vvvupi[lineIndex][domIndex][pmtIndex]) {
416  DEBUG("mismatch between pmt id of detector file " << pmtId << " and the xml file "
417  << vvvupi[lineIndex][domIndex][pmtIndex] << endl);
418  continue;
419  }
420  vhToT[runCounter - 1][lineIndex][domIndex][pmtIndex]->Fill(hitToT);
421  }
422  }
423  }
424 
425  vector<vector<vector<TGraphErrors*> > > vgToTvsHV;
426  vgToTvsHV.resize(lineNb, vector<vector<TGraphErrors*> >(domNb, vector<TGraphErrors*>(pmtNb)));
427  for (unsigned int line = 0; line < lineNb; ++line) {
428  for (int dom = 0; dom < domNb; ++dom) {
429  int fullDom = 0;
430  for (int pmt = 0; pmt < pmtNb; ++pmt) {
431  vgToTvsHV[line][dom][pmt] = new TGraphErrors();
432  int localCounter = 0;
433  for (int run = 0; run < runNb; ++run) {
434  if (vvvRunHv[run][line][dom][pmt] == 0) continue;
435  FillGraph(vgToTvsHV[line][dom][pmt], vhToTReducedChi2[run], vhToT[run][line][dom][pmt], vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt], localCounter);
436  }
437  if ((localCounter < runNb) && (localCounter > 0))
438  DEBUG("only " << localCounter << " runs out of " << runNb << " have ToT data for line " << line + 1 << " - dom " << dom + 1 << " - pmt " << pmt << endl);
439  if (localCounter == 0) ++fullDom;
440  }
441  if (fullDom > 0)
442  DEBUG(fullDom << " pmt do not have data in dom " << dom + 1 << " of line " << line + 1 << endl);
443  }
444  }
445 
446  // create a xml file with the final high voltage values.
447  fstream outCalibXml(outXMLFile.c_str(), ofstream::out);
448  struct timeval endTimeTv;
449  gettimeofday(&endTimeTv, NULL);
450  time_t endTime = (time_t)endTimeTv.tv_sec;
451  tm* endGmTime = gmtime(&endTime);
452  const string endUTCTime = asctime(endGmTime);
453 
454  ostringstream endTimeStr;
455  endTimeStr << 1900 + endGmTime->tm_year << "-" << std::setw(2) << std::setfill('0')
456  << endGmTime->tm_mon + 1 << "-" << std::setw(2) << std::setfill('0')
457  << endGmTime->tm_mday << "T" << std::setw(2) << std::setfill('0')
458  << endGmTime->tm_hour << ":" << std::setw(2) << std::setfill('0')
459  << endGmTime->tm_min << ":" << std::setw(2) << std::setfill('0')
460  << endGmTime->tm_sec << "." << std::setw(6) << std::setfill('0')
461  << endTimeTv.tv_usec << "+01:00";
462  ostringstream startTimeStr;
463  startTimeStr << 1900 + startGmTime->tm_year << "-" << std::setw(2) << std::setfill('0')
464  << startGmTime->tm_mon + 1 << "-" << std::setw(2) << std::setfill('0')
465  << startGmTime->tm_mday << "T" << std::setw(2) << std::setfill('0')
466  << startGmTime->tm_hour << ":" << std::setw(2) << std::setfill('0')
467  << startGmTime->tm_min << ":" << std::setw(2) << std::setfill('0')
468  << startGmTime->tm_sec << "." << std::setw(6) << std::setfill('0')
469  << startTimeTv.tv_usec << "+01:00";
470  outCalibXml << "<ProductTestSession>\r\n";
471  //outCalibXml << "\t<User>" << username[0] << "</User>\r\n";
472  outCalibXml << "\t<User>" << username << "</User>\r\n";
473  outCalibXml << "\t<Location>" << location << "</Location>\r\n";
474  outCalibXml << "\t<StartTime>" << startTimeStr.str() << "</StartTime>\r\n";
475  outCalibXml << "\t<EndTime>" << endTimeStr.str() << "</EndTime>\r\n";
476  outCalibXml << "\t<TestType>" << testType << "</TestType>\r\n";
477  outCalibXml << "\t<Tests>\r\n";
478  TCanvas* cToTvsHV = new TCanvas("cToTvsHV", "cToTvsHV", 800, 600);
479  TH1D* hToTvsHVReducedChi2 = new TH1D("hToTvsHVReducedChi2", "hToTvsHVReducedChi2", 101, -0.05, 10.05);
480  TH1D* hToTExpected = new TH1D("hToTExpected", "hToTExpected", 100001, -0.0005, 100.0005);
481  bool firstEvent = true;
482  vector<vector<vector<double> > > vvvHvOpt;
483  vvvHvOpt.resize(lineNb, vector<vector<double> >(domNb, vector<double>(pmtNb)));
484  vector<vector<vector<pair<double, double> > > > vvvHvRange;
485  vvvHvRange.resize(lineNb, vector<vector<pair<double, double> > >(domNb, vector<pair<double, double> >(pmtNb)));
486  FindHVRange(vvvHvRange, vvvRunHv, vvvHv);
487  int counterLinear = 0;
488  int counter2Order = 0;
489  int counter2OrderPos = 0;
490  int counter2OrderNeg = 0;
491  int counterSqrt = 0;
492  double aveLinearSlope = 0;
493  for (unsigned int line = 0; line < lineNb; ++line) {
494  for (int dom = 0; dom < domNb; ++dom) {
495  for (int pmt = 0; pmt < pmtNb; ++pmt) {
496  stringstream filename;
497  filename.str("");
498  filename.clear();
499  filename << "ToTvsHV-line" << line + 1 << "-dom" << dom + 1 << "-pmt" << pmt;
500  vgToTvsHV[line][dom][pmt]->SetName(filename.str().c_str());
501  const int pointNb = vgToTvsHV[line][dom][pmt]->GetN();
502  //if ((line == 0) && (dom == 17) && (pmt == 24)) firstEvent = true;
503  if (pointNb > 0) {
504  vgToTvsHV[line][dom][pmt]->SetLineColor(line*18 + dom*31 + pmt);
505  vgToTvsHV[line][dom][pmt]->SetMarkerColor(line*18 + dom*31 + pmt);
506  vgToTvsHV[line][dom][pmt]->SetMarkerStyle(2);
507  vgToTvsHV[line][dom][pmt]->SetMarkerSize(1);
508  if (firstEvent == true) {
509  firstEvent = false;
510  vgToTvsHV[line][dom][pmt]->SetTitle("");
511  vgToTvsHV[line][dom][pmt]->GetXaxis()->SetTitle("HV_{offset} [V]");
512  vgToTvsHV[line][dom][pmt]->GetXaxis()->CenterTitle();
513  vgToTvsHV[line][dom][pmt]->GetXaxis()->SetTitleOffset(1.2);
514  vgToTvsHV[line][dom][pmt]->GetXaxis()->SetRangeUser(-200, +200);
515  vgToTvsHV[line][dom][pmt]->GetYaxis()->SetTitle("ToT values [ns]");
516  vgToTvsHV[line][dom][pmt]->GetYaxis()->CenterTitle();
517  vgToTvsHV[line][dom][pmt]->GetYaxis()->SetTitleOffset(1.2);
518  vgToTvsHV[line][dom][pmt]->GetYaxis()->SetRangeUser(0, 80);
519  vgToTvsHV[line][dom][pmt]->Draw("APL");
520  } else {
521  vgToTvsHV[line][dom][pmt]->GetXaxis()->SetRangeUser(-200, +200);
522  vgToTvsHV[line][dom][pmt]->Draw("samePL");
523  }
524  }
525  pair<double, double> hvAndToTOpt = make_pair(0, 0);
526  string functionFormula = "";
527  if (pointNb > 2) {
528  vector<double> vParameterPol2;
529  const string pol2Name = "pol2";
530  FillParameterFitFunction(vParameterPol2, vgToTvsHV[line][dom][pmt], pol2Name, "0RQ", vvvHvRange[line][dom][pmt]);
531  vector<double> vParameterPol1;
532  const string pol1Name = "pol1";
533  FillParameterFitFunction(vParameterPol1, vgToTvsHV[line][dom][pmt], pol1Name, "0RQ", vvvHvRange[line][dom][pmt]);
534  vector<double> vParameterSqrt;
535  const string sqrtName = "[0] + sqrt([1] - [2]*x)";
536  FillParameterFitFunction(vParameterSqrt, vgToTvsHV[line][dom][pmt], sqrtName, "0RQ", vvvHvRange[line][dom][pmt]);
537  vector<double> vParameter;
538  const double testChi2 = min({fabs(vParameterPol1[0] - 1), fabs(vParameterPol2[0] - 1), fabs(vParameterSqrt[0] - 1)});
539  if (testChi2 == fabs(vParameterPol1[0] - 1)) {
540  vParameter = vParameterPol1;
541  functionFormula = pol1Name;
542  aveLinearSlope += vParameterPol1[2];
543  ++counterLinear;
544  } else if (testChi2 == fabs(vParameterPol2[0] - 1)) {
545  vParameter = vParameterPol2;
546  functionFormula = pol2Name;
547  ++counter2Order;
548  if (vParameterPol2[3] > 0) ++counter2OrderPos;
549  else ++counter2OrderNeg;
550  } else if (testChi2 == fabs(vParameterSqrt[0] - 1)) {
551  vParameter = vParameterSqrt;
552  functionFormula = sqrtName;
553  ++counterSqrt;
554  } else {
555  DEBUG("problem in the function definition, check the code at about line 500");
556  }
557  vgToTvsHV[line][dom][pmt]->Fit(functionFormula.c_str(), "RQ");
558  hToTvsHVReducedChi2->Fill(vParameter[0]);
559  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);
560  hvAndToTOpt = Solve(vParameter, functionFormula, tot);
561  if ((hvAndToTOpt.first == -1) && (hvAndToTOpt.second == -1))
562  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);
563  if ((fabs(vParameter[2]) > slopeMax) && ((functionFormula == "pol1") || (functionFormula == "pol2")))
564  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);
565  hToTExpected->Fill(hvAndToTOpt.second);
566  }
567  if (fabs(hvAndToTOpt.first) > diffHVMax) hvAndToTOpt.first = 0;
568  vvvHvOpt[line][dom][pmt] = hvAndToTOpt.first;
569  outCalibXml << "\t\t<ProductTest>\r\n";
570  outCalibXml << "\t\t\t<UPI>3.4.2.3/HAMA-R12199/" << vvvupiExt[line][dom][pmt] << "." << vvvupi[line][dom][pmt] << "</UPI>\r\n";
571  outCalibXml << "\t\t\t<TestResult>OK</TestResult>\r\n";
572  outCalibXml << "\t\t\t<TestParameters>\r\n";
573  outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
574  outCalibXml << "\t\t\t\t\t<Name>PMT_Supply_Voltage</Name>\r\n";
575  outCalibXml << "\t\t\t\t\t<Values>\r\n";
576  outCalibXml << "\t\t\t\t\t\t<string>" << hvAndToTOpt.first + vvvHv[line][dom][pmt] << "</string>\r\n";
577  outCalibXml << "\t\t\t\t\t</Values>\r\n";
578  outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
579  outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
580  outCalibXml << "\t\t\t\t\t<Name>RUN_NUMBER</Name>\r\n";
581  outCalibXml << "\t\t\t\t\t<Values>\r\n";
582  for (int run = 0; run < runNb; ++run) {
583  outCalibXml << "\t\t\t\t\t\t<string>" << runList[run] << "</string>\r\n";
584  }
585  outCalibXml << "\t\t\t\t\t</Values>\r\n";
586  outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
587  outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
588  outCalibXml << "\t\t\t\t\t<Name>PMT_Time_over_Threshold</Name>\r\n";
589  outCalibXml << "\t\t\t\t\t<Values>\r\n";
590  outCalibXml << "\t\t\t\t\t\t<string>" << tot*1e-9 << "</string>\r\n";
591  outCalibXml << "\t\t\t\t\t</Values>\r\n";
592  outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
593  outCalibXml << "\t\t\t</TestParameters>\r\n";
594  outCalibXml << "\t\t</ProductTest>\r\n";
595  }
596  }
597  }
598  cout << "linear:" << counterLinear << " - 2nd order (pos, neg):" << counter2Order
599  << "(" << counter2OrderPos << "," << counter2OrderNeg << ")"
600  << " - sqrt:" << counterSqrt << " - linear slope average: "
601  << aveLinearSlope/counterLinear << endl;
602  outCalibXml << "\t</Tests>\r\n";
603  outCalibXml << "</ProductTestSession>\r\n";
604  outCalibXml.close();
605  cToTvsHV->SaveAs("figs/cToTvsHV.png");
606  cToTvsHV->SaveAs("figs/cToTvsHV.root");
607 
608  ostringstream osstxt;
609  osstxt.clear();
610  osstxt.str("");
611  osstxt << "HVTuning-det" << detectorId << ".txt";
612 
613  ofstream outHv(osstxt.str().c_str(), ofstream::out);
614  for (unsigned int line = 0; line < lineNb; ++line) {
615  for (int dom = 0; dom < domNb; ++dom) {
616  for (int pmt = 0; pmt < pmtNb; ++pmt) {
617  outHv << line + 1 << ' ' << dom + 1 << ' ' << pmt << ' ' << vvvupiExt[line][dom][pmt] << ' ' << vvvupi[line][dom][pmt] << ' ' << vvvHv[line][dom][pmt] << ' ' << vvvHvOpt[line][dom][pmt] + vvvHv[line][dom][pmt] << endl;
618  }
619  }
620  }
621  outHv.close();
622 
623  TCanvas* cToTs[lineNb][domNb];
624  for (unsigned int line = 0; line < lineNb; ++line) {
625  for (int dom = 0; dom < domNb; ++dom) {
626  ostringstream iss;
627  iss.clear();
628  iss.str("");
629  iss << "cToTs-line" << line + 1 << "-dom" << dom + 1;
630  cToTs[line][dom] = new TCanvas(iss.str().c_str(), iss.str().c_str(), 800, 600);
631  cToTs[line][dom]->Divide(8,4);
632  for (int pmt = 0; pmt < pmtNb; ++pmt) {
633  cToTs[line][dom]->cd(pmt + 1);
634  double localMax = 0;
635  for (int run = 0; run < runNb; ++run) {
636  vhToT[run][line][dom][pmt]->SetLineColor(run + 1);
637  if (vhToT[run][line][dom][pmt]->GetMaximum() > localMax) localMax = vhToT[run][line][dom][pmt]->GetMaximum();
638  if (run == 0) {
639  vhToT[run][line][dom][pmt]->SetTitle("");
640  vhToT[run][line][dom][pmt]->SetStats(0);
641  vhToT[run][line][dom][pmt]->GetXaxis()->SetTitle("ToT values [ns]");
642  vhToT[run][line][dom][pmt]->GetXaxis()->CenterTitle();
643  vhToT[run][line][dom][pmt]->GetXaxis()->SetTitleOffset(1.2);
644  vhToT[run][line][dom][pmt]->GetXaxis()->SetRange(0, 100);
645  vhToT[run][line][dom][pmt]->Draw();
646  } else {
647  vhToT[run][line][dom][pmt]->Draw("same");
648  }
649  gPad->SetLogy();
650  }
651  vhToT[0][line][dom][pmt]->GetYaxis()->SetRangeUser(0.1, 2*localMax);
652  }
653  iss.clear();
654  iss.str("");
655  iss << "figs/cToTs-line" << line + 1 << "-dom" << dom + 1 << ".png";
656  cToTs[line][dom]->SaveAs(iss.str().c_str());
657  iss.clear();
658  iss.str("");
659  iss << "figs/cToTs-line" << line + 1 << "-dom" << dom + 1 << ".root";
660  cToTs[line][dom]->SaveAs(iss.str().c_str());
661  }
662  }
663 
664  TCanvas* cToTReducedChi2 = new TCanvas("cToTReducedChi2", "cToTReducedChi2", 800, 600);
665  for (int run = 0; run < runNb; ++run) {
666  if (run == 0) {
667  vhToTReducedChi2[run]->SetTitle("");
668  vhToTReducedChi2[run]->SetStats(0);
669  vhToTReducedChi2[run]->GetXaxis()->SetTitle("reduced #chi^{2}");
670  vhToTReducedChi2[run]->GetXaxis()->CenterTitle();
671  vhToTReducedChi2[run]->GetXaxis()->SetTitleOffset(1.2);
672  vhToTReducedChi2[run]->Draw();
673  } else {
674  vhToTReducedChi2[run]->Draw("same");
675  }
676  vhToTReducedChi2[run]->SetLineColor(run + 1);
677  }
678  cToTReducedChi2->SetLogy();
679  cToTReducedChi2->SaveAs("figs/cToT-ReducedChi2.png");
680  cToTReducedChi2->SaveAs("figs/cToT-ReducedChi2.root");
681 
682  TCanvas* cToTvsHVReducedChi2 = new TCanvas("cToTvsHVReducedChi2", "cToTvsHVReducedChi2", 800, 600);
683  hToTvsHVReducedChi2->SetTitle("");
684  hToTvsHVReducedChi2->SetStats(0);
685  hToTvsHVReducedChi2->GetXaxis()->SetTitle("reduced #chi^{2}");
686  hToTvsHVReducedChi2->GetXaxis()->CenterTitle();
687  hToTvsHVReducedChi2->GetXaxis()->SetTitleOffset(1.2);
688  hToTvsHVReducedChi2->Draw();
689  cToTvsHVReducedChi2->SetLogy();
690  cToTvsHVReducedChi2->SaveAs("figs/cToTvsHV-ReducedChi2.png");
691  cToTvsHVReducedChi2->SaveAs("figs/cToTvsHV-ReducedChi2.root");
692 
693  TCanvas* cToTExpected = new TCanvas("cToTExpected", "cToTExpected", 800, 600);
694  hToTExpected->SetTitle("");
695  hToTExpected->SetStats(0);
696  hToTExpected->GetXaxis()->SetTitle("Expected ToT value [ns]");
697  hToTExpected->GetXaxis()->CenterTitle();
698  hToTExpected->GetXaxis()->SetTitleOffset(1.2);
699  hToTExpected->Draw();
700  //cToTExpected->SetLogy();
701  cToTExpected->SaveAs("figs/cToTExpected.png");
702  cToTExpected->SaveAs("figs/cToTExpected.root");
703 
704  const int lineId1 = 0;
705  const int domId1 = 2;
706  const int pmtId1 = 5;
707  const int lineId2 = 0;
708  const int domId2 = 5;
709  const int pmtId2 = 7;
710  ostringstream issTemp;
711  issTemp << "cToTsdom" << domId1 + 1 << "pmt" << pmtId1;
712  TCanvas* cToTsdomxpmtx = new TCanvas(issTemp.str().c_str(), issTemp.str().c_str(), 800, 600);
713  bool pmtTest = true;
714  TLegend* leg1 = new TLegend(0.70,0.55,0.90,0.90);
715  leg1->SetFillStyle(0);
716  double localMaxTemp = 0;
717  if ((lineId1 > lineNb - 1) || (domId1 > domNb - 1) || (pmtId1 > pmtNb - 1))
718  pmtTest = false;
719  if (pmtTest == false) {
720  DEBUG("wrong PMT1 identification for sample plot");
721  } else {
722  ostringstream plotname0;
723  for (int run = 0; run < runNb; ++run) {
724  vhToT[run][lineId1][domId1][pmtId1]->SetLineColor(run + 1);
725  if (vhToT[run][lineId1][domId1][pmtId1]->GetMaximum() > localMaxTemp) localMaxTemp = vhToT[run][lineId1][domId1][pmtId1]->GetMaximum();
726  if (run == 0) {
727  vhToT[run][lineId1][domId1][pmtId1]->SetTitle("");
728  vhToT[run][lineId1][domId1][pmtId1]->SetStats(0);
729  vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetTitle("ToT values [ns]");
730  vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->CenterTitle();
731  vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetTitleOffset(1.2);
732  vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetRange(0, 100);
733  vhToT[run][lineId1][domId1][pmtId1]->Draw();
734  } else {
735  vhToT[run][lineId1][domId1][pmtId1]->Draw("same");
736  }
737  plotname0.str("");
738  plotname0.clear();
739  plotname0 << "HV: " << setprecision(2) << fixed
740  << vvvRunHv[run][lineId1][domId1][pmtId1] << "V";
741  leg1->AddEntry(vhToT[run][lineId1][domId1][pmtId1], plotname0.str().c_str(), "l");
742  gPad->SetLogy();
743  }
744  vhToT[0][lineId1][domId1][pmtId1]->GetYaxis()->SetRangeUser(0.1, 2*localMaxTemp);
745  leg1->Draw();
746  issTemp.clear();
747  issTemp.str("");
748  issTemp << "figs/cToTs-line" << lineId1 + 1 << "-dom" << domId1 + 1 << "-pmt" << pmtId1 << ".png";
749  cToTsdomxpmtx->SaveAs(issTemp.str().c_str());
750  issTemp.clear();
751  issTemp.str("");
752  issTemp << "figs/cToTs-line" << lineId1 + 1 << "-dom" << domId1 + 1 << "-pmt" << pmtId1 << ".root";
753  cToTsdomxpmtx->SaveAs(issTemp.str().c_str());
754  }
755 
756  issTemp.clear();
757  issTemp.str("");
758  issTemp << "cToTsdom" << domId2 + 1 << "pmt" << pmtId2;
759  TCanvas* cToTsdomxpmty = new TCanvas(issTemp.str().c_str(), issTemp.str().c_str(), 800, 600);
760  pmtTest = true;
761  TLegend* leg2 = new TLegend(0.70,0.55,0.90,0.90);
762  if ((lineId2 > lineNb - 1) || (domId2 > domNb - 1) || (pmtId2 > pmtNb - 1))
763  pmtTest = false;
764  if (pmtTest == false) {
765  DEBUG("wrong PMT2 identification for sample plot");
766  } else {
767  ostringstream plotname0;
768  for (int run = 0; run < runNb; ++run) {
769  vhToT[run][lineId2][domId2][pmtId2]->SetLineColor(run + 1);
770  if (vhToT[run][lineId2][domId2][pmtId2]->GetMaximum() > localMaxTemp) localMaxTemp = vhToT[run][lineId2][domId2][pmtId2]->GetMaximum();
771  if (run == 0) {
772  vhToT[run][lineId2][domId2][pmtId2]->SetTitle("");
773  vhToT[run][lineId2][domId2][pmtId2]->SetStats(0);
774  vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetTitle("ToT values [ns]");
775  vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->CenterTitle();
776  vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetTitleOffset(1.2);
777  vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetRange(0, 100);
778  vhToT[run][lineId2][domId2][pmtId2]->Draw();
779  } else {
780  vhToT[run][lineId2][domId2][pmtId2]->Draw("same");
781  }
782  plotname0.str("");
783  plotname0.clear();
784  plotname0 << "HV: " << setprecision(2) << fixed
785  << vvvRunHv[run][lineId2][domId2][pmtId2] << "V";
786  leg2->AddEntry(vhToT[run][lineId2][domId2][pmtId2], plotname0.str().c_str(), "l");
787  gPad->SetLogy();
788  }
789  vhToT[0][lineId2][domId2][pmtId2]->GetYaxis()->SetRangeUser(0.1, 2*localMaxTemp);
790  leg2->Draw();
791  issTemp.clear();
792  issTemp.str("");
793  issTemp << "figs/cToTs-line" << lineId2 + 1 << "-dom" << domId2 + 1 << "-pmt" << pmtId2 << ".png";
794  cToTsdomxpmty->SaveAs(issTemp.str().c_str());
795  issTemp.clear();
796  issTemp.str("");
797  issTemp << "figs/cToTs-line" << lineId2 + 1 << "-dom" << domId2 + 1 << "-pmt" << pmtId2 << ".root";
798  cToTsdomxpmty->SaveAs(issTemp.str().c_str());
799  }
800 
801  DEBUG("analysis duration: " << difftime(endTime, startTime) << endl);
802 }
AccessDB
void AccessDB(std::vector< std::vector< std::vector< int > > > &vvvupi, std::vector< std::vector< std::vector< int > > > &vvvupiExt, std::vector< std::vector< std::vector< double > > > &vvvHv, std::string &cookie, std::string &username, std::string &streamName, std::string &detectorId, const int runId, const std::string &hvDefault, std::vector< int > &lineIds)
Definition: JHVTuning.cc:54
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1410
FillParameterFitFunction
void FillParameterFitFunction(std::vector< double > &vParameter, TGraphErrors *&localGraph, const std::string &curFunction, const std::string &curOptions, const std::pair< double, double > range)
Definition: JHVTuning.cc:241
FindHVRange
void FindHVRange(std::vector< std::vector< std::vector< std::pair< double, double > > > > &vvvHvRange, std::vector< std::vector< std::vector< std::vector< double > > > > &vvvRunHv, std::vector< std::vector< std::vector< double > > > &vvvHv)
Definition: JHVTuning.cc:221
JROOT::counter_type
Long64_t counter_type
Type definition for counter.
Definition: JCounter.hh:24
JLANG::JObjectMultiplexer
Auxiliary class for multiplexing object iterators.
Definition: JObjectIterator.hh:35
SolveAnalytic
double SolveAnalytic(const double deg0, const double deg1, const double deg2, const double tot)
Definition: JHVTuning.cc:149
KM3NETDAQ::JDAQChronometer::getRunNumber
int getRunNumber() const
Get run number.
Definition: JDAQChronometer.hh:121
JDetector.hh
Data structure for detector geometry and calibration.
KM3NETDAQ::JDAQHit
Hit data structure.
Definition: JDAQHit.hh:36
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
KM3NETDAQ::JDAQTimeslice
Data time slice.
Definition: JDAQTimeslice.hh:36
JLANG::getline
std::istream & getline(std::istream &in, JString &object)
Read string from input stream until end of line.
Definition: JString.hh:468
JLANG::JObjectMultiplexer::next
virtual const pointer_type & next()
Get next element.
Definition: JObjectMultiplexer.hh:76
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:444
debug
int debug
debug level
Definition: JSirene.cc:59
JMessage.hh
General purpose messaging.
FATAL
#define FATAL(A)
Definition: JMessage.hh:65
JMultipleFileScanner.hh
Scanning of objects from multiple files according a format that follows from the extension of each fi...
JModuleRouter.hh
Direct access to module in detector data structure.
Solve
std::pair< double, double > Solve(std::vector< double > &vParameter, const std::string &curFunction, const double tot)
Definition: JHVTuning.cc:161
FillGraph
void FillGraph(TGraphErrors *&graph, TH1D *&hChi2, TH1D *&histo, const double hv, int &localCounter)
Definition: JHVTuning.cc:181
JParser.hh
Utility class to parse command line options.
JSupport.hh
ROOT TTree parameter settings.
JLANG::JObjectMultiplexer::hasNext
virtual bool hasNext()
Check availability of next element.
Definition: JObjectMultiplexer.hh:62
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
main
int main(int argc, char *argv[])
Definition: Main.cpp:15
Generated by   doxygen 1.8.5