Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JHVTuning.cc File Reference

Program for production of the xml file with nominal high voltage values. More...

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

Go to the source code of this file.

Functions

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)
 
double SolveAnalytic (const double deg0, const double deg1, const double deg2, const double tot)
 
pair< double, double > Solve (vector< double > &vParameter, const string &curFunction, const double tot)
 
void FillGraph (TGraphErrors *&graph, TH1D *&hChi2, TH1D *&histo, const double hv, int &localCounter)
 
void FindHVRange (vector< vector< vector< pair< double, double > > > > &vvvHvRange, vector< vector< vector< vector< double > > > > &vvvRunHv, vector< vector< vector< double > > > &vvvHv)
 
void FillParameterFitFunction (vector< double > &vParameter, TGraphErrors *&localGraph, const string &curFunction, const string &curOptions, const pair< double, double > range)
 
int main (int argc, char **argv)
 

Detailed Description

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

Definition in file JHVTuning.cc.

Function Documentation

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 
)
inline

Definition at line 64 of file JHVTuning.cc.

64  {
65  try {
66  JDB::reset(usr, pwd, cookie);
67  usr = JDB::get()->User();
68  if (runId == "-2")
69  cout << "opening connection as user " << usr << " for detector information" << endl;
70  else if (runId == "-1")
71  cout << "opening connection as user " << usr << " for default HV values" << endl;
72  else
73  cout << "opening connection as user " << usr << " for run " << runId << endl;
74  const string detoid = getDetector(detid);
75  JSelector selection;
76  selection.add("detoid", detoid);
77  if (query == "detectors") {
78  //selection.add("detid", detid);
79  JDatabaseObjectIterator<JDetectors> in(query, selection);
80  while (in.hasNext()) {
81  JDetectors* curdet = in.next();
82  string curoid = curdet->OID;
83  if (curoid == detoid) {
84  location = curdet->LOCATIONID;
85  }
86  }
87  } else if (query == "pmt_hv_run_settings") {
88  selection.add("run", runId);
90  while (in.hasNext()) {
91  JPMTHVRunSettings* rs = in.next();
92  int lineId = 0;
93  for (set<int>::iterator it = lineIds.begin(); it != lineIds.end(); ++it) {
94  if (*it == rs->DUID) break;
95  ++lineId;
96  }
97  vvvHv[lineId][rs->FLOORID - 1][rs->PMTINTID] = rs->HV_VALUE;
98  }
99  } else if (query == "pmt_hv_settings") {
101  while (in.hasNext()) {
102  JPMTHVSettings* rs = in.next();
103  int lineId = 0;
104  for (set<int>::iterator it = lineIds.begin(); it != lineIds.end(); ++it) {
105  if (*it == rs->DUID) break;
106  ++lineId;
107  }
108  vvvHv[lineId][rs->FLOORID - 1][rs->PMTINTID] = rs->HV_VALUE;
109  vvvupi[lineId][rs->FLOORID - 1][rs->PMTINTID] = rs->PMTUPI;
110  }
111  }
112  }
113  catch(const exception& error) {
114  FATAL(error.what() << endl);
115  }
116 }
Object iteration from database.
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
Definition: JSirene.sh:45
T get(const JHead &header)
Get object from header.
Auxiliary class for specifying selection of database data.
std::string LOCATIONID
Definition: JDetectors.hh:25
JSelector & add(const JSelector &selection)
Add selection.
#define FATAL(A)
Definition: JMessage.hh:67
static const JDetectorsHelper & getDetector
Function object for mapping serial number and object identifier of detectors.
Definition: JDBToolkit.hh:131
void reset(T &value)
Reset value.
double SolveAnalytic ( const double  deg0,
const double  deg1,
const double  deg2,
const double  tot 
)
inline

Definition at line 121 of file JHVTuning.cc.

121  {
122  const double delta = deg1*deg1 - 4*deg2*(deg0 - tot);
123  if ((delta < 0) || (deg2 == 0)) return deg1? (tot - deg0)/deg1 : 0;
124  const double x1 = (-deg1 + sqrt(delta))/(2*deg2);
125  const double x2 = (-deg1 - sqrt(delta))/(2*deg2);
126  if (deg2 > 0) return (x1 < x2)? x1 : x2;
127  else return (x1 < x2)? x2 : x1;
128  //return (fabs(x1) < fabs(x2))? x1 : x2;
129 }
pair<double, double> Solve ( vector< double > &  vParameter,
const string &  curFunction,
const double  tot 
)
inline

Definition at line 133 of file JHVTuning.cc.

133  {
134  pair<double, double> curSol;
135  if (curFunction == "pol1") {
136  curSol.first = SolveAnalytic(vParameter[1], vParameter[2], 0, tot);
137  curSol.second = vParameter[2]*curSol.first + vParameter[1];
138  } else if (curFunction == "pol2") {
139  curSol.first = SolveAnalytic(vParameter[1], vParameter[2], vParameter[3], tot);
140  curSol.second = vParameter[3]*pow(curSol.first, 2) + vParameter[2]*curSol.first + vParameter[1];
141  } else if (curFunction == "[0] + sqrt([1] - [2]*x)") {
142  curSol.first = (vParameter[2] - pow(tot - vParameter[1], 2))/vParameter[3];
143  curSol.second = vParameter[1] + sqrt(vParameter[2] - vParameter[3]*curSol.first);
144  } else {
145  curSol.first = -1;
146  curSol.second = -1;
147  }
148  return curSol;
149 }
double SolveAnalytic(const double deg0, const double deg1, const double deg2, const double tot)
Definition: JHVTuning.cc:121
void FillGraph ( TGraphErrors *&  graph,
TH1D *&  hChi2,
TH1D *&  histo,
const double  hv,
int &  localCounter 
)
inline

Definition at line 153 of file JHVTuning.cc.

153  {
154  const int firstBin = 15;
155  const int lastBin = histo->GetXaxis()->GetLast() - 60;
156  histo->GetXaxis()->SetRange(firstBin, lastBin);
157  if (!histo->GetEntries()) return;
158  if (!histo->Integral()) return;
159  const int maxBin = histo->GetMaximumBin();
160  double bestRedChi2 = 1e6;
161  double bestMean = 0;
162  for (int binInf = maxBin - 6; binInf < maxBin - 2; ++binInf) {
163  for (int binSup = maxBin + 1; binSup <= maxBin + 6; ++binSup) {
164  if (binSup - binInf < 4) continue;
165  TF1* gausFit = new TF1("gausFit", "gaus", binInf, binSup);
166  const double curInteg = histo->Integral(binInf, binSup);
167  if (curInteg == 0) continue;
168  histo->Fit(gausFit, "RQ");
169  if (!gausFit->GetNDF()) {
170  gausFit->Delete();
171  continue;
172  }
173  const double fitMean = gausFit->GetParameter(1);
174  const double fitSigma = gausFit->GetParameter(2);
175  const double reducedChi2 = gausFit->GetChisquare()/gausFit->GetNDF();
176  if ((std::fabs(reducedChi2 - 1) < std::fabs(bestRedChi2 - 1)) && (reducedChi2 > 0.5) && (fitSigma < 10)) {
177  bestRedChi2 = reducedChi2;
178  bestMean = fitMean;
179  }
180  gausFit->Delete();
181  }
182  }
183  hChi2->Fill(bestRedChi2);
184  histo->GetXaxis()->SetRange(0, lastBin);
185  if (bestMean > 0) {
186  graph->SetPoint(localCounter, hv, bestMean);
187  ++localCounter;
188  }
189 }
void FindHVRange ( vector< vector< vector< pair< double, double > > > > &  vvvHvRange,
vector< vector< vector< vector< double > > > > &  vvvRunHv,
vector< vector< vector< double > > > &  vvvHv 
)
inline

Definition at line 193 of file JHVTuning.cc.

193  {
194  for (unsigned int line = 0; line < vvvHvRange.size(); ++line) {
195  for (unsigned int dom = 0; dom < vvvHvRange[line].size(); ++dom) {
196  for (unsigned int pmt = 0; pmt < vvvHvRange[line][dom].size(); ++pmt) {
197  double minHvDiff = 1e8;
198  double maxHvDiff = 0;
199  for (unsigned int run = 0; run < vvvRunHv.size(); ++run) {
200  if ((vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt]) < minHvDiff)
201  minHvDiff = vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt];
202  if (((vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt]) > maxHvDiff) &&(vvvRunHv[run][line][dom][pmt]))
203  maxHvDiff = vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt];
204  }
205  vvvHvRange[line][dom][pmt] = std::make_pair(minHvDiff, maxHvDiff);
206  }
207  }
208  }
209 }
esac $JPP_DIR examples JDetector JTransitTime o $OUTPUT_FILE n N $NPE T $TTS_NS d $DEBUG for HISTOGRAM in tts tt2 pmt
Definition: JTransitTime.sh:36
void FillParameterFitFunction ( vector< double > &  vParameter,
TGraphErrors *&  localGraph,
const string &  curFunction,
const string &  curOptions,
const pair< double, double >  range 
)
inline

Definition at line 213 of file JHVTuning.cc.

213  {
214  TF1* fLocal = new TF1(curFunction.c_str(), curFunction.c_str(), range.first, range.second);
215  if (curFunction == "[0] + sqrt([1] - [2]*x)")
216  fLocal->SetParameter(1, range.second + 25.0);
217  localGraph->Fit(fLocal, curOptions.c_str());
218  const int ndf = fLocal->GetNDF();
219  const double chi2 = fLocal->GetChisquare();
220  double reducedChi2 = 0;
221  if (ndf) {
222  reducedChi2 = chi2/ndf;
223  }
224  vParameter.push_back(reducedChi2);
225  for (int p = 0; p < fLocal->GetNpar(); ++p) {
226  vParameter.push_back(fLocal->GetParameter(p));
227  }
228 }
int main ( int  argc,
char **  argv 
)

Definition at line 230 of file JHVTuning.cc.

231 {
232 
233  gStyle->SetLegendBorderSize(0);
234  struct timeval startTimeTv;
235  gettimeofday(&startTimeTv, NULL);
236  time_t startTime = (time_t)startTimeTv.tv_sec;
237  tm* startGmTime = gmtime(&startTime);
238  const string startUTCTime = asctime(startGmTime);
239 
241  int debug;
242  string detectorFile;
243  string outXMLFile;
244  double tot;
245  string usr;
246  string pwd;
247  string mycookies;
248  string testType;
249  double redChi2Max;
250  double slopeMax;
251  try {
252  JParser<> zap;
253  zap['f'] = make_field(vInputFiles);
254  zap['a'] = make_field(detectorFile);
255  zap['o'] = make_field(outXMLFile);
256  zap['d'] = make_field(debug) = 1;
257  zap['t'] = make_field(tot) = 26.4;
258  zap['u'] = make_field(usr) = "";
259  zap['!'] = make_field(pwd) = "";
260  zap['C'] = make_field(mycookies) = "";
261  zap['T'] = make_field(testType) = "HV-TUNING-V2";
262  zap['c'] = make_field(redChi2Max) = 10;
263  zap['P'] = make_field(slopeMax) = 1;
264  zap(argc, argv);
265  }
266  catch(const exception& error) {
267  FATAL(error.what() << endl);
268  }
269 
271  try {
272  parameters = getTriggerParameters(vInputFiles);
273  }
274  catch(const JException& error) {
275  FATAL("No trigger parameters from input:" << error.what() << endl);
276  }
277  cout.tie(&cerr);
278 
280  try {
281  load(detectorFile, detector);
282  }
283  catch(const JException& error) {
284  FATAL(error);
285  }
286 
287  const JModuleRouter moduleRouter(detector);
288  const int lineNb = getNumberOfStrings(detector);
289  const set<int> lineIds = getStringIDs(detector);
290  const int runNb = vInputFiles.size();
291  const int pmtNb = NUMBER_OF_PMTS;
292  const int domNb = getNumberOfFloors(detector);
293  const double diffHVMax = 250;
294 
295  string cookiesContent;
296  ifstream incookies(mycookies.c_str(), ofstream::in);
297  getline(incookies, cookiesContent);
298 
299  // stream the vendor hvs from the db
301  vvvHv.resize(lineNb, vector<vector<double> > (domNb, vector<double>(pmtNb)));
303  vvvRunHv.resize(runNb, vector<vector<vector<double> > > (lineNb, vector<vector<double> > (domNb, vector<double>(pmtNb))));
304  vector<vector<vector<string> > > vvvupi;
305  vvvupi.resize(lineNb, vector<vector<string> >(domNb, vector<string>(pmtNb)));
306  int detid = detector.getID();
307  string streamName = "detectors";
308  string location = "";
309  AccessDB(vvvupi, vvvHv, cookiesContent, usr, pwd, streamName, detid, "-2", lineIds, location);
310  streamName = "pmt_hv_settings";
311  AccessDB(vvvupi, vvvHv, cookiesContent, usr, pwd, streamName, detid, "-1", lineIds, location);
312 
313  // prepare the ToT histo
315  vhToT.resize(runNb, vector<vector<vector<TH1D*> > >(lineNb, vector<vector<TH1D*> >(domNb, vector<TH1D*>(pmtNb))));
316  vector<TH1D*> vhToTReducedChi2;
317  vhToTReducedChi2.resize(runNb);
318  for (int run = 0; run < runNb; ++run) {
319  ostringstream plotname;
320  plotname << "hToTReducedChi2-run" << run + 1;
321  vhToTReducedChi2[run] = new TH1D(plotname.str().c_str(), plotname.str().c_str(), 101, -0.05, 10.05);
322  for (int line = 0; line < lineNb; ++line) {
323  for (int dom = 0; dom < domNb; ++dom) {
324  for (int pmt = 0; pmt < pmtNb; ++pmt) {
325  plotname.str("");
326  plotname.clear();
327  plotname << "hToT-run" << run + 1 << "-line" << line + 1 << "-dom" << dom + 1 << "-pmt" << pmt;
328  vhToT[run][line][dom][pmt] = new TH1D(plotname.str().c_str(), plotname.str().c_str(), 301, -0.5, 300.5);
329  }
330  }
331  }
332  }
333 
334  // read the data and fill the ToT histo
335  vector<int> runList;
336  int runCounter = 0;
337  int oldRun = -1;
338  JROOTClassSelector selector("JDAQTimesliceL1");
340  counter_type counter = 0;
341  for ( ; in.hasNext() && counter != vInputFiles.getLimit(); ++counter) {
342  JDAQTimeslice* timeslice = in.next();
343  const int newRun = timeslice->getRunNumber();
344  if (newRun != oldRun) {
345  streamName = "pmt_hv_run_settings";
346  stringstream runIdSStr;
347  runIdSStr << newRun;
348  const string newRunStr = runIdSStr.str();
349  AccessDB(vvvupi, vvvRunHv[runCounter], cookiesContent, usr, pwd, streamName, detid, newRunStr, lineIds, location);
350  ++runCounter;
351  runList.push_back(newRun);
352  }
353  oldRun = newRun;
354  for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
355  const int domId = super_frame->getModuleID();
356  const JModule module = moduleRouter.getModule(domId);
357  const int lineIndex = module.getString() - 1;
358  const int domIndex = module.getFloor() - 1;
359  for(JDAQSuperFrame::const_iterator hit = super_frame->begin(); hit!=super_frame->end(); ++hit){
360  const double hitToT = hit->getToT();
361  const int pmtIndex = hit->getPMT();
362  const int pmtId = module.getPMT(pmtIndex).getID();
363  const string pmtIdDBStr = vvvupi[lineIndex][domIndex][pmtIndex].substr(22);
364  const int pmtIdDB = atoi(pmtIdDBStr.c_str());
365  if (pmtId != pmtIdDB) {
366  DEBUG("mismatch between pmt ids of detector file " << pmtId << " and database "
367  << vvvupi[lineIndex][domIndex][pmtIndex] << endl);
368  continue;
369  }
370  vhToT[runCounter - 1][lineIndex][domIndex][pmtIndex]->Fill(hitToT);
371  }
372  }
373  }
374 
376  vgToTvsHV.resize(lineNb, vector<vector<TGraphErrors*> >(domNb, vector<TGraphErrors*>(pmtNb)));
377  for (int line = 0; line < lineNb; ++line) {
378  for (int dom = 0; dom < domNb; ++dom) {
379  int fullDom = 0;
380  for (int pmt = 0; pmt < pmtNb; ++pmt) {
381  vgToTvsHV[line][dom][pmt] = new TGraphErrors();
382  int localCounter = 0;
383  for (int run = 0; run < runNb; ++run) {
384  if (vvvRunHv[run][line][dom][pmt] == 0) continue;
385  FillGraph(vgToTvsHV[line][dom][pmt], vhToTReducedChi2[run], vhToT[run][line][dom][pmt], vvvRunHv[run][line][dom][pmt] - vvvHv[line][dom][pmt], localCounter);
386  }
387  if ((localCounter < runNb) && (localCounter > 0))
388  DEBUG("only " << localCounter << " runs out of " << runNb << " have ToT data for line " << line + 1 << " - dom " << dom + 1 << " - pmt " << pmt << endl);
389  if (localCounter == 0) ++fullDom;
390  }
391  if (fullDom > 0)
392  DEBUG(fullDom << " pmt do not have data in dom " << dom + 1 << " of line " << line + 1 << endl);
393  }
394  }
395 
396  // create a xml file with the final high voltage values.
397  fstream outCalibXml(outXMLFile.c_str(), ofstream::out);
398  struct timeval endTimeTv;
399  gettimeofday(&endTimeTv, NULL);
400  time_t endTime = (time_t)endTimeTv.tv_sec;
401  tm* endGmTime = gmtime(&endTime);
402  const string endUTCTime = asctime(endGmTime);
403 
404  ostringstream endTimeStr;
405  endTimeStr << 1900 + endGmTime->tm_year << "-" << std::setw(2) << std::setfill('0')
406  << endGmTime->tm_mon + 1 << "-" << std::setw(2) << std::setfill('0')
407  << endGmTime->tm_mday << "T" << std::setw(2) << std::setfill('0')
408  << endGmTime->tm_hour << ":" << std::setw(2) << std::setfill('0')
409  << endGmTime->tm_min << ":" << std::setw(2) << std::setfill('0')
410  << endGmTime->tm_sec << "." << std::setw(6) << std::setfill('0')
411  << endTimeTv.tv_usec << "+01:00";
412  ostringstream startTimeStr;
413  startTimeStr << 1900 + startGmTime->tm_year << "-" << std::setw(2) << std::setfill('0')
414  << startGmTime->tm_mon + 1 << "-" << std::setw(2) << std::setfill('0')
415  << startGmTime->tm_mday << "T" << std::setw(2) << std::setfill('0')
416  << startGmTime->tm_hour << ":" << std::setw(2) << std::setfill('0')
417  << startGmTime->tm_min << ":" << std::setw(2) << std::setfill('0')
418  << startGmTime->tm_sec << "." << std::setw(6) << std::setfill('0')
419  << startTimeTv.tv_usec << "+01:00";
420  outCalibXml << "<ProductTestSession>\r\n";
421  outCalibXml << "\t<User>" << usr << "</User>\r\n";
422  outCalibXml << "\t<Location>" << location << "</Location>\r\n";
423  outCalibXml << "\t<StartTime>" << startTimeStr.str() << "</StartTime>\r\n";
424  outCalibXml << "\t<EndTime>" << endTimeStr.str() << "</EndTime>\r\n";
425  outCalibXml << "\t<TestType>" << testType << "</TestType>\r\n";
426  outCalibXml << "\t<Tests>\r\n";
427  TCanvas* cToTvsHV = new TCanvas("cToTvsHV", "cToTvsHV", 800, 600);
428  TH1D* hToTvsHVReducedChi2 = new TH1D("hToTvsHVReducedChi2", "hToTvsHVReducedChi2", 101, -0.05, 10.05);
429  TH1D* hToTExpected = new TH1D("hToTExpected", "hToTExpected", 100001, -0.0005, 100.0005);
430  bool firstEvent = true;
431  vector<vector<vector<double> > > vvvHvOpt;
432  vvvHvOpt.resize(lineNb, vector<vector<double> >(domNb, vector<double>(pmtNb)));
434  vvvHvRange.resize(lineNb, vector<vector<pair<double, double> > >(domNb, vector<pair<double, double> >(pmtNb)));
435  FindHVRange(vvvHvRange, vvvRunHv, vvvHv);
436  int counterLinear = 0;
437  int counter2Order = 0;
438  int counter2OrderPos = 0;
439  int counter2OrderNeg = 0;
440  int counterSqrt = 0;
441  double aveLinearSlope = 0;
442  for (int line = 0; line < lineNb; ++line) {
443  for (int dom = 0; dom < domNb; ++dom) {
444  for (int pmt = 0; pmt < pmtNb; ++pmt) {
445  stringstream filename;
446  filename.str("");
447  filename.clear();
448  filename << "ToTvsHV-line" << line + 1 << "-dom" << dom + 1 << "-pmt" << pmt;
449  vgToTvsHV[line][dom][pmt]->SetName(filename.str().c_str());
450  const int pointNb = vgToTvsHV[line][dom][pmt]->GetN();
451  //if ((line == 0) && (dom == 17) && (pmt == 24)) firstEvent = true;
452  if (pointNb > 0) {
453  vgToTvsHV[line][dom][pmt]->SetLineColor(line*18 + dom*31 + pmt);
454  vgToTvsHV[line][dom][pmt]->SetMarkerColor(line*18 + dom*31 + pmt);
455  vgToTvsHV[line][dom][pmt]->SetMarkerStyle(2);
456  vgToTvsHV[line][dom][pmt]->SetMarkerSize(1);
457  if (firstEvent == true) {
458  firstEvent = false;
459  vgToTvsHV[line][dom][pmt]->SetTitle("");
460  vgToTvsHV[line][dom][pmt]->GetXaxis()->SetTitle("HV_{offset} [V]");
461  vgToTvsHV[line][dom][pmt]->GetXaxis()->CenterTitle();
462  vgToTvsHV[line][dom][pmt]->GetXaxis()->SetTitleOffset(1.2);
463  vgToTvsHV[line][dom][pmt]->GetXaxis()->SetRangeUser(-200, +200);
464  vgToTvsHV[line][dom][pmt]->GetYaxis()->SetTitle("ToT values [ns]");
465  vgToTvsHV[line][dom][pmt]->GetYaxis()->CenterTitle();
466  vgToTvsHV[line][dom][pmt]->GetYaxis()->SetTitleOffset(1.2);
467  vgToTvsHV[line][dom][pmt]->GetYaxis()->SetRangeUser(0, 80);
468  vgToTvsHV[line][dom][pmt]->Draw("APL");
469  } else {
470  vgToTvsHV[line][dom][pmt]->GetXaxis()->SetRangeUser(-200, +200);
471  vgToTvsHV[line][dom][pmt]->Draw("samePL");
472  }
473  }
474  pair<double, double> hvAndToTOpt = make_pair(0, 0);
475  string functionFormula = "";
476  if (pointNb > 2) {
477  vector<double> vParameterPol2;
478  const string pol2Name = "pol2";
479  FillParameterFitFunction(vParameterPol2, vgToTvsHV[line][dom][pmt], pol2Name, "0RQ", vvvHvRange[line][dom][pmt]);
480  vector<double> vParameterPol1;
481  const string pol1Name = "pol1";
482  FillParameterFitFunction(vParameterPol1, vgToTvsHV[line][dom][pmt], pol1Name, "0RQ", vvvHvRange[line][dom][pmt]);
483  vector<double> vParameterSqrt;
484  const string sqrtName = "[0] + sqrt([1] - [2]*x)";
485  FillParameterFitFunction(vParameterSqrt, vgToTvsHV[line][dom][pmt], sqrtName, "0RQ", vvvHvRange[line][dom][pmt]);
486  vector<double> vParameter;
487  const double testChi2 = min({fabs(vParameterPol1[0] - 1), fabs(vParameterPol2[0] - 1), fabs(vParameterSqrt[0] - 1)});
488  if (testChi2 == fabs(vParameterPol1[0] - 1)) {
489  vParameter = vParameterPol1;
490  functionFormula = pol1Name;
491  aveLinearSlope += vParameterPol1[2];
492  ++counterLinear;
493  } else if (testChi2 == fabs(vParameterPol2[0] - 1)) {
494  vParameter = vParameterPol2;
495  functionFormula = pol2Name;
496  ++counter2Order;
497  if (vParameterPol2[3] > 0) ++counter2OrderPos;
498  else ++counter2OrderNeg;
499  } else if (testChi2 == fabs(vParameterSqrt[0] - 1)) {
500  vParameter = vParameterSqrt;
501  functionFormula = sqrtName;
502  ++counterSqrt;
503  } else {
504  DEBUG("problem in the function definition, check the code at about line 500");
505  }
506  vgToTvsHV[line][dom][pmt]->Fit(functionFormula.c_str(), "RQ");
507  hToTvsHVReducedChi2->Fill(vParameter[0]);
508  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);
509  hvAndToTOpt = Solve(vParameter, functionFormula, tot);
510  if ((hvAndToTOpt.first == -1) && (hvAndToTOpt.second == -1))
511  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);
512  if ((fabs(vParameter[2]) > slopeMax) && ((functionFormula == "pol1") || (functionFormula == "pol2")))
513  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);
514  hToTExpected->Fill(hvAndToTOpt.second);
515  }
516  if (fabs(hvAndToTOpt.first) > diffHVMax) hvAndToTOpt.first = 0;
517  vvvHvOpt[line][dom][pmt] = hvAndToTOpt.first;
518  outCalibXml << "\t\t<ProductTest>\r\n";
519  outCalibXml << "\t\t\t<UPI>" << vvvupi[line][dom][pmt] << "</UPI>\r\n";
520  outCalibXml << "\t\t\t<TestResult>OK</TestResult>\r\n";
521  outCalibXml << "\t\t\t<TestParameters>\r\n";
522  outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
523  outCalibXml << "\t\t\t\t\t<Name>PMT_Supply_Voltage</Name>\r\n";
524  outCalibXml << "\t\t\t\t\t<Values>\r\n";
525  outCalibXml << "\t\t\t\t\t\t<string>" << hvAndToTOpt.first + vvvHv[line][dom][pmt] << "</string>\r\n";
526  outCalibXml << "\t\t\t\t\t</Values>\r\n";
527  outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
528  outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
529  outCalibXml << "\t\t\t\t\t<Name>RUN_NUMBER</Name>\r\n";
530  outCalibXml << "\t\t\t\t\t<Values>\r\n";
531  for (int run = 0; run < runNb; ++run) {
532  outCalibXml << "\t\t\t\t\t\t<string>" << runList[run] << "</string>\r\n";
533  }
534  outCalibXml << "\t\t\t\t\t</Values>\r\n";
535  outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
536  outCalibXml << "\t\t\t\t<ProductTestParameter>\r\n";
537  outCalibXml << "\t\t\t\t\t<Name>PMT_Time_over_Threshold</Name>\r\n";
538  outCalibXml << "\t\t\t\t\t<Values>\r\n";
539  outCalibXml << "\t\t\t\t\t\t<string>" << tot*1e-9 << "</string>\r\n";
540  outCalibXml << "\t\t\t\t\t</Values>\r\n";
541  outCalibXml << "\t\t\t\t</ProductTestParameter>\r\n";
542  outCalibXml << "\t\t\t</TestParameters>\r\n";
543  outCalibXml << "\t\t</ProductTest>\r\n";
544  }
545  }
546  }
547  cout << "linear:" << counterLinear << " - 2nd order (pos, neg):" << counter2Order
548  << "(" << counter2OrderPos << "," << counter2OrderNeg << ")"
549  << " - sqrt:" << counterSqrt << " - linear slope average: "
550  << aveLinearSlope/counterLinear << endl;
551  outCalibXml << "\t</Tests>\r\n";
552  outCalibXml << "</ProductTestSession>\r\n";
553  outCalibXml.close();
554  cToTvsHV->SaveAs("figs/cToTvsHV.png");
555  cToTvsHV->SaveAs("figs/cToTvsHV.root");
556 
557  ostringstream osstxt;
558  osstxt.clear();
559  osstxt.str("");
560  osstxt << "results/HVTuning-det" << detid << ".txt";
561 
562  ofstream outHv(osstxt.str().c_str(), ofstream::out);
563  for (int line = 0; line < lineNb; ++line) {
564  for (int dom = 0; dom < domNb; ++dom) {
565  for (int pmt = 0; pmt < pmtNb; ++pmt) {
566  outHv << line + 1 << ' ' << dom + 1 << ' ' << pmt << ' ' << vvvupi[line][dom][pmt] << ' ' << vvvHv[line][dom][pmt] << ' ' << vvvHvOpt[line][dom][pmt] + vvvHv[line][dom][pmt] << endl;
567  }
568  }
569  }
570  outHv.close();
571 
572  TCanvas* cToTs[lineNb][domNb];
573  for (int line = 0; line < lineNb; ++line) {
574  for (int dom = 0; dom < domNb; ++dom) {
575  ostringstream iss;
576  iss.clear();
577  iss.str("");
578  iss << "cToTs-line" << line + 1 << "-dom" << dom + 1;
579  cToTs[line][dom] = new TCanvas(iss.str().c_str(), iss.str().c_str(), 800, 600);
580  cToTs[line][dom]->Divide(8,4);
581  for (int pmt = 0; pmt < pmtNb; ++pmt) {
582  cToTs[line][dom]->cd(pmt + 1);
583  double localMax = 0;
584  for (int run = 0; run < runNb; ++run) {
585  vhToT[run][line][dom][pmt]->SetLineColor(run + 1);
586  if (vhToT[run][line][dom][pmt]->GetMaximum() > localMax) localMax = vhToT[run][line][dom][pmt]->GetMaximum();
587  if (run == 0) {
588  vhToT[run][line][dom][pmt]->SetTitle("");
589  vhToT[run][line][dom][pmt]->SetStats(0);
590  vhToT[run][line][dom][pmt]->GetXaxis()->SetTitle("ToT values [ns]");
591  vhToT[run][line][dom][pmt]->GetXaxis()->CenterTitle();
592  vhToT[run][line][dom][pmt]->GetXaxis()->SetTitleOffset(1.2);
593  vhToT[run][line][dom][pmt]->GetXaxis()->SetRange(0, 100);
594  vhToT[run][line][dom][pmt]->Draw();
595  } else {
596  vhToT[run][line][dom][pmt]->Draw("same");
597  }
598  gPad->SetLogy();
599  }
600  vhToT[0][line][dom][pmt]->GetYaxis()->SetRangeUser(0.1, 2*localMax);
601  }
602  iss.clear();
603  iss.str("");
604  iss << "figs/cToTs-line" << line + 1 << "-dom" << dom + 1 << ".png";
605  cToTs[line][dom]->SaveAs(iss.str().c_str());
606  iss.clear();
607  iss.str("");
608  iss << "figs/cToTs-line" << line + 1 << "-dom" << dom + 1 << ".root";
609  cToTs[line][dom]->SaveAs(iss.str().c_str());
610  }
611  }
612 
613  TCanvas* cToTReducedChi2 = new TCanvas("cToTReducedChi2", "cToTReducedChi2", 800, 600);
614  for (int run = 0; run < runNb; ++run) {
615  if (run == 0) {
616  vhToTReducedChi2[run]->SetTitle("");
617  vhToTReducedChi2[run]->SetStats(0);
618  vhToTReducedChi2[run]->GetXaxis()->SetTitle("reduced #chi^{2}");
619  vhToTReducedChi2[run]->GetXaxis()->CenterTitle();
620  vhToTReducedChi2[run]->GetXaxis()->SetTitleOffset(1.2);
621  vhToTReducedChi2[run]->Draw();
622  } else {
623  vhToTReducedChi2[run]->Draw("same");
624  }
625  vhToTReducedChi2[run]->SetLineColor(run + 1);
626  }
627  cToTReducedChi2->SetLogy();
628  cToTReducedChi2->SaveAs("figs/cToT-ReducedChi2.png");
629  cToTReducedChi2->SaveAs("figs/cToT-ReducedChi2.root");
630 
631  TCanvas* cToTvsHVReducedChi2 = new TCanvas("cToTvsHVReducedChi2", "cToTvsHVReducedChi2", 800, 600);
632  hToTvsHVReducedChi2->SetTitle("");
633  hToTvsHVReducedChi2->SetStats(0);
634  hToTvsHVReducedChi2->GetXaxis()->SetTitle("reduced #chi^{2}");
635  hToTvsHVReducedChi2->GetXaxis()->CenterTitle();
636  hToTvsHVReducedChi2->GetXaxis()->SetTitleOffset(1.2);
637  hToTvsHVReducedChi2->Draw();
638  cToTvsHVReducedChi2->SetLogy();
639  cToTvsHVReducedChi2->SaveAs("figs/cToTvsHV-ReducedChi2.png");
640  cToTvsHVReducedChi2->SaveAs("figs/cToTvsHV-ReducedChi2.root");
641 
642  TCanvas* cToTExpected = new TCanvas("cToTExpected", "cToTExpected", 800, 600);
643  hToTExpected->SetTitle("");
644  hToTExpected->SetStats(0);
645  hToTExpected->GetXaxis()->SetTitle("Expected ToT value [ns]");
646  hToTExpected->GetXaxis()->CenterTitle();
647  hToTExpected->GetXaxis()->SetTitleOffset(1.2);
648  hToTExpected->Draw();
649  //cToTExpected->SetLogy();
650  cToTExpected->SaveAs("figs/cToTExpected.png");
651  cToTExpected->SaveAs("figs/cToTExpected.root");
652 
653  const int lineId1 = 0;
654  const int domId1 = 2;
655  const int pmtId1 = 5;
656  const int lineId2 = 0;
657  const int domId2 = 5;
658  const int pmtId2 = 7;
659  ostringstream issTemp;
660  issTemp << "cToTsdom" << domId1 + 1 << "pmt" << pmtId1;
661  TCanvas* cToTsdomxpmtx = new TCanvas(issTemp.str().c_str(), issTemp.str().c_str(), 800, 600);
662  bool pmtTest = true;
663  TLegend* leg1 = new TLegend(0.70,0.55,0.90,0.90);
664  leg1->SetFillStyle(0);
665  double localMaxTemp = 0;
666  if ((lineId1 > lineNb - 1) || (domId1 > domNb - 1) || (pmtId1 > pmtNb - 1))
667  pmtTest = false;
668  if (pmtTest == false) {
669  DEBUG("wrong PMT1 identification for sample plot");
670  } else {
671  ostringstream plotname0;
672  for (int run = 0; run < runNb; ++run) {
673  vhToT[run][lineId1][domId1][pmtId1]->SetLineColor(run + 1);
674  if (vhToT[run][lineId1][domId1][pmtId1]->GetMaximum() > localMaxTemp) localMaxTemp = vhToT[run][lineId1][domId1][pmtId1]->GetMaximum();
675  if (run == 0) {
676  vhToT[run][lineId1][domId1][pmtId1]->SetTitle("");
677  vhToT[run][lineId1][domId1][pmtId1]->SetStats(0);
678  vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetTitle("ToT values [ns]");
679  vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->CenterTitle();
680  vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetTitleOffset(1.2);
681  vhToT[run][lineId1][domId1][pmtId1]->GetXaxis()->SetRange(0, 100);
682  vhToT[run][lineId1][domId1][pmtId1]->Draw();
683  } else {
684  vhToT[run][lineId1][domId1][pmtId1]->Draw("same");
685  }
686  plotname0.str("");
687  plotname0.clear();
688  plotname0 << "HV: " << setprecision(2) << fixed
689  << vvvRunHv[run][lineId1][domId1][pmtId1] << "V";
690  leg1->AddEntry(vhToT[run][lineId1][domId1][pmtId1], plotname0.str().c_str(), "l");
691  gPad->SetLogy();
692  }
693  vhToT[0][lineId1][domId1][pmtId1]->GetYaxis()->SetRangeUser(0.1, 2*localMaxTemp);
694  leg1->Draw();
695  issTemp.clear();
696  issTemp.str("");
697  issTemp << "figs/cToTs-line" << lineId1 + 1 << "-dom" << domId1 + 1 << "-pmt" << pmtId1 << ".png";
698  cToTsdomxpmtx->SaveAs(issTemp.str().c_str());
699  issTemp.clear();
700  issTemp.str("");
701  issTemp << "figs/cToTs-line" << lineId1 + 1 << "-dom" << domId1 + 1 << "-pmt" << pmtId1 << ".root";
702  cToTsdomxpmtx->SaveAs(issTemp.str().c_str());
703  }
704 
705  issTemp.clear();
706  issTemp.str("");
707  issTemp << "cToTsdom" << domId2 + 1 << "pmt" << pmtId2;
708  TCanvas* cToTsdomxpmty = new TCanvas(issTemp.str().c_str(), issTemp.str().c_str(), 800, 600);
709  pmtTest = true;
710  TLegend* leg2 = new TLegend(0.70,0.55,0.90,0.90);
711  if ((lineId2 > lineNb - 1) || (domId2 > domNb - 1) || (pmtId2 > pmtNb - 1))
712  pmtTest = false;
713  if (pmtTest == false) {
714  DEBUG("wrong PMT2 identification for sample plot");
715  } else {
716  ostringstream plotname0;
717  for (int run = 0; run < runNb; ++run) {
718  vhToT[run][lineId2][domId2][pmtId2]->SetLineColor(run + 1);
719  if (vhToT[run][lineId2][domId2][pmtId2]->GetMaximum() > localMaxTemp) localMaxTemp = vhToT[run][lineId2][domId2][pmtId2]->GetMaximum();
720  if (run == 0) {
721  vhToT[run][lineId2][domId2][pmtId2]->SetTitle("");
722  vhToT[run][lineId2][domId2][pmtId2]->SetStats(0);
723  vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetTitle("ToT values [ns]");
724  vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->CenterTitle();
725  vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetTitleOffset(1.2);
726  vhToT[run][lineId2][domId2][pmtId2]->GetXaxis()->SetRange(0, 100);
727  vhToT[run][lineId2][domId2][pmtId2]->Draw();
728  } else {
729  vhToT[run][lineId2][domId2][pmtId2]->Draw("same");
730  }
731  plotname0.str("");
732  plotname0.clear();
733  plotname0 << "HV: " << setprecision(2) << fixed
734  << vvvRunHv[run][lineId2][domId2][pmtId2] << "V";
735  leg2->AddEntry(vhToT[run][lineId2][domId2][pmtId2], plotname0.str().c_str(), "l");
736  gPad->SetLogy();
737  }
738  vhToT[0][lineId2][domId2][pmtId2]->GetYaxis()->SetRangeUser(0.1, 2*localMaxTemp);
739  leg2->Draw();
740  issTemp.clear();
741  issTemp.str("");
742  issTemp << "figs/cToTs-line" << lineId2 + 1 << "-dom" << domId2 + 1 << "-pmt" << pmtId2 << ".png";
743  cToTsdomxpmty->SaveAs(issTemp.str().c_str());
744  issTemp.clear();
745  issTemp.str("");
746  issTemp << "figs/cToTs-line" << lineId2 + 1 << "-dom" << domId2 + 1 << "-pmt" << pmtId2 << ".root";
747  cToTsdomxpmty->SaveAs(issTemp.str().c_str());
748  }
749 
750  DEBUG("analysis duration: " << difftime(endTime, startTime) << endl);
751 }
Utility class to parse command line options.
Definition: JParser.hh:1493
General exception.
Definition: JException.hh:23
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
Data structure for a composite optical module.
Definition: JModule.hh:50
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)
Definition: JHVTuning.cc:64
Detector data structure.
Definition: JDetector.hh:80
Auxiliary class to select ROOT class based on class name.
Router for direct addressing of module data in detector data structure.
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings IDs.
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
Definition: JSirene.sh:45
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
void FillParameterFitFunction(vector< double > &vParameter, TGraphErrors *&localGraph, const string &curFunction, const string &curOptions, const pair< double, double > range)
Definition: JHVTuning.cc:213
Long64_t counter_type
Type definition for counter.
Auxiliary class for multiplexing object iterators.
int getRunNumber() const
Get run number.
esac $JPP_DIR examples JDetector JTransitTime o $OUTPUT_FILE n N $NPE T $TTS_NS d $DEBUG for HISTOGRAM in tts tt2 pmt
Definition: JTransitTime.sh:36
Hit data structure.
Definition: JDAQHit.hh:34
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
int getID() const
Get identifier.
Definition: JObjectID.hh:55
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
Data time slice.
std::istream & getline(std::istream &in, JString &object)
Read string from input stream until end of line.
Definition: JString.hh:478
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:61
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:174
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
#define FATAL(A)
Definition: JMessage.hh:67
int getString() const
Get string number.
Definition: JLocation.hh:134
void FindHVRange(vector< vector< vector< pair< double, double > > > > &vvvHvRange, vector< vector< vector< vector< double > > > > &vvvRunHv, vector< vector< vector< double > > > &vvvHv)
Definition: JHVTuning.cc:193
void FillGraph(TGraphErrors *&graph, TH1D *&hChi2, TH1D *&histo, const double hv, int &localCounter)
Definition: JHVTuning.cc:153
General purpose class for object reading from a list of file names.
virtual const char * what() const
Get error message.
Definition: JException.hh:48
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
pair< double, double > Solve(vector< double > &vParameter, const string &curFunction, const double tot)
Definition: JHVTuning.cc:133