CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
L1TOccupancyClient.cc
Go to the documentation of this file.
2 
10 #include <cstdio>
11 #include <sstream>
12 #include <cmath>
13 #include <vector>
14 #include <TMath.h>
15 #include <climits>
16 #include <TFile.h>
17 #include <TDirectory.h>
18 #include <TProfile.h>
19 
20 using namespace std;
21 using namespace edm;
22 
23 //____________________________________________________________________________
24 // Function: L1TOccupancyClient
25 // Description: This is the constructor, basic variable initialization
26 // Inputs:
27 // * const edm::ParameterSet& ps = Parameter for this analyzer
28 //____________________________________________________________________________
30  // Get parameters
31  parameters_ = ps;
32  verbose_ = ps.getParameter<bool>("verbose");
33  tests_ = ps.getParameter<std::vector<ParameterSet> >("testParams");
34 
35  if (verbose_) {
36  cout << "[L1TOccupancyClient:] Called constructor" << endl;
37  }
38 }
39 
40 //____________________________________________________________________________
41 // Function: ~L1TOccupancyClient
42 // Description: This is the destructor, basic variable deletion
43 //____________________________________________________________________________
45  if (verbose_) {
46  cout << "[L1TOccupancyClient:] Called destructor" << endl;
47  }
48 }
49 
50 //____________________________________________________________________________
51 // Function: beginRun
52 // Description: This is will be run at the begining of each run
53 // Inputs:
54 // * const Run& r = Run information
55 // * const EventSetup& context = Event Setup information
56 //____________________________________________________________________________
58  hservice_ = new L1TOccupancyClientHistogramService(parameters_, ibooker, verbose_);
59 
60  if (verbose_) {
61  cout << "[L1TOccupancyClient:] Called beginRun" << endl;
62 
63  // In verbose mode we will produce an extra output file with several tests
64  file_ = TFile::Open("DQM_L1TOccupancyClient_Snapshots_LS.root", "RECREATE");
65  }
66 
67  ibooker.setCurrentFolder("L1T/L1TOccupancy/");
68  //dbe_->setCurrentFolder("L1T/L1TOccupancy/Results");
69  //dbe_->setCurrentFolder("L1T/L1TOccupancy/BadCellValues");
70  //dbe_->setCurrentFolder("L1T/L1TOccupancy/Certification");
71 
72  // Loop over all tests in defined
73  for (vector<ParameterSet>::iterator it = tests_.begin(); it != tests_.end(); it++) {
74  // If the test algorithm is XYSymmetry we create the necessary histograms
75  if ((*it).getUntrackedParameter<string>("algoName", "XYSymmetry") == "XYSymmetry") {
76  // Getting Parameters for the test
77  string testName = (*it).getParameter<string>("testName");
78  ParameterSet algoParameters = (*it).getParameter<ParameterSet>("algoParams");
79  string histPath = algoParameters.getParameter<string>("histPath");
80 
81  if (verbose_) {
82  cout << "[L1TOccupancyClient:] Monitored histogram path: " << histPath << endl;
83 
84  // Creating verbose file directory structure
85  // test_name/test_name_Results,
86  // test_name/test_name_Histos
87  // TDirectory *td = file_->mkdir(testName.c_str() ,testName.c_str());
88  //FIXME: sub never used gcc361 warning
89  //TDirectory *sub = td ->mkdir((testName+"_Results").c_str(),string("_Results").c_str());
90 
91  //sub = td->mkdir((testName+"_Histos").c_str() ,(testName+"_Histos").c_str());
92  //sub = td->mkdir((testName+"_Histos_AllLS").c_str(),(testName+"_Histos_AllLS").c_str());
93  }
94 
95  // Load histograms in service instance
96  if (hservice_->loadHisto(igetter, testName, histPath)) {
97  // Mask channels specified in python file
98  hservice_->setMaskedBins(testName, algoParameters.getParameter<vector<ParameterSet> >("maskedAreas"));
99 
100  // Book MonitorElements
101  // * Test results
102  ibooker.setCurrentFolder("L1T/L1TOccupancy/Results");
103  string title = testName;
104  MonitorElement* m = ibooker.book2D(title.c_str(), hservice_->getDifferentialHistogram(testName));
105  m->setTitle(title);
106  m->Reset();
107  meResults[title] = m;
108 
109  // * Which cells are masked as bad
110  ibooker.setCurrentFolder("L1T/L1TOccupancy/HistogramDiff");
111  title = testName;
112  m = ibooker.book2D(title.c_str(), hservice_->getDifferentialHistogram(testName));
113  m->Reset();
114  m->setTitle(title);
115  meDifferential[title] = m;
116 
117  // * Fraction of bad cells
118  ibooker.setCurrentFolder("L1T/L1TOccupancy/Certification");
119  title = testName;
120  m = ibooker.book1D(title.c_str(), title.c_str(), 2500, -.5, 2500. - .5);
121  m->setTitle(title);
122  meCertification[title] = m;
123 
124  mValidTests.push_back(&(*it));
125  }
126  }
127  }
128 }
129 
130 //____________________________________________________________________________
131 // Function: endRun
132 // Description: This is will be run at the end of each run
133 // Inputs:
134 // * const Run& r = Run information
135 // * const EventSetup& context = Event Setup information
136 //____________________________________________________________________________
138  book(ibooker, igetter);
139 
140  if (verbose_) {
141  cout << "[L1TOccupancyClient:] Called endRun()" << endl;
142  }
143 
144  // Loop over every test in python
145  for (std::vector<ParameterSet*>::iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
146  ParameterSet& test = (**it);
147  string algo_name = test.getUntrackedParameter<string>("algoName", "XYSymmetry");
148  string test_name = test.getParameter<string>("testName");
149 
150  if (verbose_) {
151  cout << "[L1TOccupancyClient:] Starting calculations for: " << algo_name << " on: " << test_name << endl;
152  }
153 
154  if (algo_name == "XYSymmetry") {
155  ParameterSet ps = (**it).getParameter<ParameterSet>("algoParams");
156  string histPath = ps.getParameter<string>("histPath");
157 
158  vector<pair<int, double> > deadChannels;
159  vector<pair<int, double> > statDev;
160  bool enoughStats = false;
161 
162  // Make final block
163  hservice_->updateHistogramEndRun(test_name);
164 
165  // Perform the test
166  double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
167  stringstream str;
168  str << test_name << "_cumu_LS_EndRun";
169 
170  if (verbose_) {
171  TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
172 
173  cumulative_save->SetTitle(str.str().c_str());
174 
175  TDirectory* td = file_->GetDirectory(test_name.c_str());
176 
177  td->cd(string(test_name + "_Histos_AllLS").c_str());
178 
179  cumulative_save->Write();
180  }
181  // If we have enough statistics, we can write test result
182  if (enoughStats) {
183  // Make the result histogram
184  printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
185 
186  if (verbose_) {
187  TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
188  cumulative_save->SetTitle(str.str().c_str());
189  TDirectory* td = file_->GetDirectory(("DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
190  td->cd(string(test_name + "_Histos").c_str());
191  cumulative_save->Write();
192 
193  // save the result histo
194  TH2F* h2f = meResults[test_name]->getTH2F();
195  stringstream str2;
196  str2 << test_name << "_result_LS_EndRun";
197  TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
198 
199  td->cd(string(test_name + "_Results").c_str());
200  dead_save->SetTitle(str2.str().c_str());
201  dead_save->Write();
202  }
203 
204  // Updating test results
205  meDifferential[test_name]->Reset();
206  meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
207 
208  vector<int> lsCertification = hservice_->getLSCertification(test_name);
209 
210  // Fill fraction of dead channels
211  for (unsigned int i = 0; i < lsCertification.size(); i++) {
212  int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
213  meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
214  }
215 
216  // Reset differential histo
217  hservice_->resetHisto(test_name);
218 
219  if (verbose_) {
220  cout << "Now we have enough statstics for " << test_name << endl;
221  }
222 
223  } else {
224  if (verbose_) {
225  cout << "we don't have enough statstics for " << test_name << endl;
226  }
227 
228  // Getting LS which this test monitored
229  vector<int> lsCertification = hservice_->getLSCertification(test_name);
230 
231  // Fill fraction of dead channels
232  for (unsigned int i = 0; i < lsCertification.size(); i++) {
233  int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
234  meCertification[test_name]->getTH1()->SetBinContent(bin, -1);
235  }
236  }
237  } else {
238  if (verbose_) {
239  cout << "No valid algorithm" << std::endl;
240  }
241  }
242  }
243 
244  if (verbose_) {
245  file_->Close();
246  }
247 
248  delete hservice_;
249 }
250 
251 //____________________________________________________________________________
252 // Function: beginLuminosityBlock
253 // Description: This is will be run at the begining of each luminosity block
254 // Inputs:
255 // * const LuminosityBlock& lumiSeg = Luminosity Block information
256 // * const EventSetup& context = Event Setup information
257 //____________________________________________________________________________
258 
259 //____________________________________________________________________________
260 // Function: endLuminosityBlock
261 // Description: This is will be run at the end of each luminosity block
262 // Inputs:
263 // * const LuminosityBlock& lumiSeg = Luminosity Block information
264 // * const EventSetup& context = Event Setup information
265 //____________________________________________________________________________
267  DQMStore::IGetter& igetter,
268  const edm::LuminosityBlock& lumiSeg,
269  const edm::EventSetup& c) {
270  book(ibooker, igetter);
271 
272  int eventLS = lumiSeg.id().luminosityBlock();
273 
274  if (verbose_) {
275  cout << "[L1TOccupancyClient:] Called endLuminosityBlock()" << endl;
276  cout << "[L1TOccupancyClient:] Lumisection: " << eventLS << endl;
277  }
278 
279  // Loop over every test in python
280  for (std::vector<ParameterSet*>::const_iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
281  ParameterSet& test = (**it);
282  string algo_name = test.getUntrackedParameter<string>("algoName", "XYSymmetry");
283  string test_name = test.getParameter<string>("testName");
284 
285  if (verbose_) {
286  cout << "[L1TOccupancyClient:] Starting calculations for " << algo_name << " on:" << test_name << endl;
287  }
288 
289  if (algo_name == "XYSymmetry") {
290  ParameterSet ps = (**it).getParameter<ParameterSet>("algoParams");
291  string histPath = ps.getParameter<string>("histPath");
292 
293  vector<pair<int, double> > deadChannels;
294  vector<pair<int, double> > statDev;
295  bool enoughStats = false;
296 
297  // Update histo's data with data of this LS
298  hservice_->updateHistogramEndLS(igetter, test_name, histPath, eventLS);
299 
300  // Perform the test
301  double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
302  stringstream str;
303  str << test_name << "_cumu_LS_" << eventLS;
304 
305  if (verbose_) {
306  TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
307  cumulative_save->SetTitle(str.str().c_str());
308  TDirectory* td = file_->GetDirectory(test_name.c_str());
309  td->cd(string(test_name + "_Histos_AllLS").c_str());
310  cumulative_save->Write();
311  }
312 
313  // If we have enough statistics, we can write test result
314  if (enoughStats) {
315  // Make the result histogram
316  printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
317 
318  if (verbose_) {
319  TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
320  cumulative_save->SetTitle(str.str().c_str());
321  TDirectory* td = file_->GetDirectory(("DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
322  td->cd(string(test_name + "_Histos").c_str());
323  cumulative_save->Write();
324 
325  // save the result histo
326  TH2F* h2f = meResults[test_name]->getTH2F();
327  stringstream str2;
328  str2 << test_name << "_result_LS_" << eventLS;
329  TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
330 
331  td->cd(string(test_name + "_Results").c_str());
332  dead_save->SetTitle(str2.str().c_str());
333  dead_save->Write();
334  }
335 
336  // Updating test results
337  meDifferential[test_name]->Reset();
338  meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
339 
340  vector<int> lsCertification = hservice_->getLSCertification(test_name);
341 
342  // Fill fraction of dead channels
343  for (unsigned int i = 0; i < lsCertification.size(); i++) {
344  int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
345  meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
346  }
347 
348  // Reset differential histo
349  hservice_->resetHisto(test_name);
350 
351  if (verbose_) {
352  cout << "Now we have enough statstics for " << test_name << endl;
353  }
354 
355  } else {
356  if (verbose_) {
357  cout << "we don't have enough statstics for " << test_name << endl;
358  }
359  }
360  } else {
361  if (verbose_) {
362  cout << "No valid algorithm" << std::endl;
363  }
364  }
365  }
366 }
367 
368 //____________________________________________________________________________
369 // Function: analyze
370 // Description: This is will be run for every event
371 // Inputs:
372 // * const Event& e = Event information
373 // * const EventSetup& context = Event Setup information
374 //____________________________________________________________________________
375 //void L1TOccupancyClient::analyze(const Event& e, const EventSetup& context){}
376 
377 //____________________________________________________________________________
378 // Function: xySymmetry
379 // Description: This method preforms the XY Symmetry test
380 // Inputs:
381 // * ParameterSet ps = Parameters for the test
382 // * std::string test_name = Test name of the test to be executed
383 // * std::vector< pair<int,double> >& deadChannels = Vector of
384 // * std::vector< pair<int,double> >& statDev =
385 // * bool& enoughStats =
386 // Outputs:
387 // * double = fraction of bins that failed test, DeadChannels in vector, in: ParameterSet of test parameters
388 //____________________________________________________________________________
390  string iTestName,
391  vector<pair<int, double> >& deadChannels,
392  vector<pair<int, double> >& statDev,
393  bool& enoughStats) {
394  // Getting differential histogram for this this thes
395  TH2F* diffHist = hservice_->getDifferentialHistogram(iTestName);
396 
397  int pAxis = ps.getUntrackedParameter<int>("axis", 1);
398  int pAverageMode = ps.getUntrackedParameter<int>("averageMode", 2); // 1=arith. mean, 2=median
399  int nBinsX = diffHist->GetNbinsX(); // actual number of bins x
400  int nBinsY = diffHist->GetNbinsY(); // actual number of bins y
401 
402  // Axis==1 : Means symmetry axis is vertical
403  if (pAxis == 1) {
404  int maxBinStrip, centralBinStrip; // x-coordinate of strips
405 
406  maxBinStrip = nBinsX;
407 
408  // If takeCenter=true determine central bin of the pAxis
409  // If takeCenter=false determine the bin to use based user input
410  if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
411  centralBinStrip = nBinsX / 2 + 1;
412  } else {
413  double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
414  getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 1);
415  }
416 
417  // Assuming odd number of strips --> first comparison is middle strip to itself
418  int upBinStrip = centralBinStrip;
419  int lowBinStrip = centralBinStrip;
420 
421  // If even number decrease lowBinstrip by one
422  if (nBinsX % 2 == 0) {
423  lowBinStrip--;
424  }
425 
426  // Do we have enough statistics? Min(Max(strip_i,strip_j))>threshold
427  std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
428 
429  int nActualStrips = 0; //number of strips that are not fully masked
430  for (int i = 0, j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; i++, j++, k--) {
431  double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsY, j, pAverageMode);
432  double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsY, k, pAverageMode);
433 
434  // Protection for when both strips are masked
435  if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
436  maxAvgs[i] = TMath::Max(avg1, avg2);
437  nActualStrips++;
438  }
439  }
440 
441  vector<double> defaultMu0up;
442  defaultMu0up.push_back(13.7655);
443  defaultMu0up.push_back(184.742);
444  defaultMu0up.push_back(50735.3);
445  defaultMu0up.push_back(-97.6793);
446 
447  TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
448  vector<double> params = ps.getUntrackedParameter<vector<double> >("params_mu0_up", defaultMu0up);
449  for (unsigned int i = 0; i < params.size(); i++) {
450  tf.SetParameter(i, params[i]);
451  }
452  int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
453 
454  vector<double> defaultMu0low;
455  defaultMu0low.push_back(2.19664);
456  defaultMu0low.push_back(1.94546);
457  defaultMu0low.push_back(-99.3263);
458  defaultMu0low.push_back(19.388);
459 
460  params = ps.getUntrackedParameter<vector<double> >("params_mu0_low", defaultMu0low);
461  for (unsigned int i = 0; i < params.size(); i++) {
462  tf.SetParameter(i, params[i]);
463  }
464  int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
465 
466  if (verbose_) {
467  cout << "nbins: " << hservice_->getNBinsHistogram(iTestName) << endl;
468  cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
469  }
470 
471  enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
472  if (verbose_) {
473  cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
474  << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
475  << ", threshold: " << TMath::Max(statsup, statslow) << endl;
476  }
477 
478  //if enough statistics
479  //make the test
480  if (enoughStats) {
481  for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
482  double avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, upBinStrip, pAverageMode);
483  compareWithStrip(
484  diffHist, iTestName, lowBinStrip, nBinsY, pAxis, avg, ps, deadChannels); //compare with lower side
485 
486  avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, lowBinStrip, pAverageMode);
487  compareWithStrip(
488  diffHist, iTestName, upBinStrip, nBinsY, pAxis, avg, ps, deadChannels); //compare with upper side
489  }
490  }
491  }
492 
493  // pAxis==2 : Means symetry pAxis is horizontal
494  else if (pAxis == 2) {
495  int maxBinStrip, centralBinStrip; //x-coordinate of strips
496 
497  maxBinStrip = nBinsY;
498 
499  // Determine center of diagram: either with set pAxis or middle of diagram
500  if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
501  centralBinStrip = nBinsY / 2 + 1;
502  } else {
503  double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
504  getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 2);
505  }
506 
507  //assuming odd number of strips --> first comparison is middle strip to itself
508  int lowBinStrip = centralBinStrip, upBinStrip = centralBinStrip;
509 
510  //if even number
511  if (nBinsX % 2 == 0) {
512  //decrease lowBinstrip by one
513  lowBinStrip--;
514  }
515 
516  //do we have enough statistics? Min(Max(strip_i,strip_j))>threshold
517  std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
518  int nActualStrips = 0;
519  for (int i = 0, j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; i++, j++, k--) {
520  double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsX, j, pAverageMode);
521  double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsX, k, pAverageMode);
522  if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
523  maxAvgs[i] = TMath::Max(avg1, avg2);
524  nActualStrips++;
525  }
526  }
527 
528  vector<double> defaultMu0up;
529  defaultMu0up.push_back(13.7655);
530  defaultMu0up.push_back(184.742);
531  defaultMu0up.push_back(50735.3);
532  defaultMu0up.push_back(-97.6793);
533 
534  vector<double> params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_up", defaultMu0up);
535  TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
536  for (unsigned int i = 0; i < params.size(); i++) {
537  tf.SetParameter(i, params[i]);
538  }
539  int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
540 
541  vector<double> defaultMu0low;
542  defaultMu0low.push_back(2.19664);
543  defaultMu0low.push_back(1.94546);
544  defaultMu0low.push_back(-99.3263);
545  defaultMu0low.push_back(19.388);
546 
547  params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_low", defaultMu0low);
548  for (unsigned int i = 0; i < params.size(); i++) {
549  tf.SetParameter(i, params[i]);
550  }
551  int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
552  if (verbose_) {
553  cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
554  }
555  enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
556  if (verbose_) {
557  cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
558  << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
559  << ", threshold: " << TMath::Max(statsup, statslow) << endl;
560  }
561 
562  //if we have enough statistics
563  //make the test
564  if (enoughStats) {
565  for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
566  double avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, upBinStrip, pAverageMode);
567  compareWithStrip(
568  diffHist, iTestName, lowBinStrip, nBinsX, pAxis, avg, ps, deadChannels); //compare with lower side
569 
570  avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, lowBinStrip, pAverageMode);
571  compareWithStrip(
572  diffHist, iTestName, upBinStrip, nBinsX, pAxis, avg, ps, deadChannels); //compare with upper side
573  }
574  }
575  } else {
576  if (verbose_) {
577  cout << "Invalid axis" << endl;
578  }
579  }
580 
581  return (deadChannels.size() - hservice_->getNBinsMasked(iTestName)) * 1.0 / hservice_->getNBinsHistogram(iTestName);
582 }
583 
584 //____________________________________________________________________________
585 // Function: getAvrg
586 // Description: Calculate strip average with method iAvgMode, where strip is
587 // prependicular to iAxis at bin iBinStrip of histogram iHist
588 // Inputs:
589 // * TH2F* iHist = Histogram to be tested
590 // * string iTestName = Name of the test
591 // * int iAxis = Axis prependicular to plot symmetry
592 // * int iNBins = Number of bins in the strip
593 // * int iBinStrip = Bin corresponding to the strip in iAxis
594 // * int iAvgMode = Type of average mode 1) Average 2) Median
595 // Outputs:
596 // * double = Average of input strip
597 //____________________________________________________________________________
598 double L1TOccupancyClient::getAvrg(TH2F* iHist, string iTestName, int iAxis, int iNBins, int iBinStrip, int iAvgMode) {
599  double avg = 0.0;
600  TH1D* proj = nullptr;
601  TH2F* histo = new TH2F(*iHist);
602 
603  std::vector<double> values;
604  int marked;
605 
606  if (iAxis == 1) {
607  switch (iAvgMode) {
608  // arithmetic average
609  case 1:
610  marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
611  proj = histo->ProjectionX();
612  avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
613  break;
614 
615  // median
616  case 2:
617  marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
618  proj = histo->ProjectionY("_py", iBinStrip, iBinStrip);
619  for (int i = 0; i < iNBins; i++) {
620  values.push_back(proj->GetBinContent(i + 1));
621  }
622  avg = TMath::Median(iNBins, &values[0]);
623  break;
624  default:
625  if (verbose_) {
626  cout << "Invalid averaging mode!" << endl;
627  }
628  break;
629  }
630  } else if (iAxis == 2) {
631  switch (iAvgMode) {
632  // arithmetic average
633  case 1:
634  marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
635  proj = histo->ProjectionY();
636  avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
637  break;
638  // median
639  case 2:
640  marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
641  proj = histo->ProjectionX("_px", iBinStrip, iBinStrip);
642  for (int i = 0; i < iNBins; i++) {
643  values.push_back(proj->GetBinContent(i + 1));
644  }
645 
646  avg = TMath::Median(iNBins, &values[0]);
647  break;
648  default:
649  if (verbose_) {
650  cout << "invalid averaging mode!" << endl;
651  }
652  break;
653  }
654  } else {
655  if (verbose_) {
656  cout << "invalid axis" << endl;
657  }
658  }
659  delete proj;
660  delete histo;
661  return avg;
662 }
663 
664 //____________________________________________________________________________
665 // Function: printDeadChannels
666 // Description:
667 // Inputs:
668 // * vector< pair<int,double> > iDeadChannels = List of bin that are masked of failed tthe test
669 // * TH2F* oHistDeadChannels = Histogram where test results should be printed
670 // * vector< pair<int,double> > statDev = ???
671 // * string iTestName = Name of the test
672 //____________________________________________________________________________
673 void L1TOccupancyClient::printDeadChannels(const vector<pair<int, double> >& iDeadChannels,
674  TH2F* oHistDeadChannels,
675  const vector<std::pair<int, double> >& statDev,
676  string iTestName) {
677  // Reset the dead channels histogram
678  oHistDeadChannels->Reset();
679  if (verbose_) {
680  cout << "suspect or masked channels of " << iTestName << ": ";
681  }
682 
683  int x, y, z;
684  float chi2 = 0.0;
685 
686  // put all bad (value=1) and masked (value=-1) cells in histo
687  for (std::vector<pair<int, double> >::const_iterator it = iDeadChannels.begin(); it != iDeadChannels.end(); it++) {
688  int bin = (*it).first;
689  oHistDeadChannels->GetBinXYZ(bin, x, y, z);
690 
691  if (hservice_->isMasked(iTestName, x, y)) {
692  oHistDeadChannels->SetBinContent(bin, -1);
693  if (verbose_) {
694  printf("(%4i,%4i) Masked\n", x, y);
695  }
696  } else {
697  oHistDeadChannels->SetBinContent(bin, 1);
698  if (verbose_) {
699  printf("(%4i,%4i) Failed test\n", x, y);
700  }
701  }
702  }
703 
704  // FIXME: Is this needed?
705  for (std::vector<pair<int, double> >::const_iterator it = statDev.begin(); it != statDev.end(); it++) {
706  double dev = (*it).second;
707  chi2 += dev;
708  }
709  //put total chi2 in float
710 
711  if (verbose_) {
712  cout << "total number of suspect channels: " << (iDeadChannels.size() - (hservice_->getNBinsMasked(iTestName)))
713  << endl;
714  }
715 }
716 
717 //____________________________________________________________________________
718 // Function: compareWithStrip
719 // Description: Evaluates statistical compatibility of a strip (cell by cell) against a given average
720 // Inputs:
721 // * TH2F* iHist = Histogram to be tested
722 // * string iTestName = Which test to apply
723 // * int iBinStrip = Bin Coordinate (in bin units) of the stripo
724 // * int iNBins = Number of Bins in the strip
725 // * int iAxis = Which Axis is prependicular to the plot symmetry.
726 // * double iAvg = Average of the strip
727 // * ParameterSet iPS = Parameters for the test
728 // * vector<pair<int,double> >& oChannels = Output of bin that are masked or failed the test
729 // Outputs:
730 // * int = Number of dead channels
731 //____________________________________________________________________________
733  string iTestName,
734  int iBinStrip,
735  int iNBins,
736  int iAxis,
737  double iAvg,
738  const ParameterSet& iPS,
739  vector<pair<int, double> >& oChannels) {
740  int dead = 0;
741 
742  //
743  if (iAxis == 1) {
744  // Get and set parameters for working curves
745  TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
746  TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
747  fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
748  fmuup->SetParameter(1, iAvg);
749  fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
750  fmulow->SetParameter(1, iAvg);
751 
752  TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
753 
754  // Evaluate sigma up
755  vector<double> defaultChi2up;
756  defaultChi2up.push_back(5.45058e-05);
757  defaultChi2up.push_back(0.268756);
758  defaultChi2up.push_back(-11.7515);
759 
760  vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
761  for (unsigned int i = 0; i < params.size(); i++) {
762  fchi->SetParameter(i, params[i]);
763  }
764  double sigma_up = fchi->Eval(iAvg);
765 
766  // Evaluate sigma low
767  vector<double> defaultChi2low;
768  defaultChi2low.push_back(4.11095e-05);
769  defaultChi2low.push_back(0.577451);
770  defaultChi2low.push_back(-10.378);
771 
772  params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
773  for (unsigned int i = 0; i < params.size(); i++) {
774  fchi->SetParameter(i, params[i]);
775  }
776  double sigma_low = fchi->Eval(iAvg);
777 
778  if (verbose_) {
779  cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
780  }
781 
782  for (int i = 1; i <= iNBins; i++) {
783  if (verbose_) {
784  cout << " " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(iBinStrip, i))
785  << " low: " << fmulow->Eval(iHist->GetBinContent(iBinStrip, i)) << endl;
786  }
787 
788  // Evaluate chi2 for cells
789  double muup = fmuup->Eval(iHist->GetBinContent(iBinStrip, i));
790  double mulow = fmulow->Eval(iHist->GetBinContent(iBinStrip, i));
791 
792  // If channel is masked -> set it to value -1
793  if (hservice_->isMasked(iTestName, iBinStrip, i)) {
794  oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
795  }
796  //else perform test
797  else if (muup > sigma_up || mulow > sigma_low ||
798  ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
799  (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
800  dead++;
801  oChannels.push_back(
802  pair<int, double>(iHist->GetBin(iBinStrip, i), abs(iHist->GetBinContent(iBinStrip, i) - iAvg) / iAvg));
803  }
804  }
805  }
806  //
807  else if (iAxis == 2) {
808  //get and set parameters for working curves
809  TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
810  TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
811  fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
812  fmuup->SetParameter(1, iAvg);
813  fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
814  fmulow->SetParameter(1, iAvg);
815 
816  TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
817 
818  // Evaluate sigma up
819  vector<double> defaultChi2up;
820  defaultChi2up.push_back(5.45058e-05);
821  defaultChi2up.push_back(0.268756);
822  defaultChi2up.push_back(-11.7515);
823 
824  vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
825  for (unsigned int i = 0; i < params.size(); i++) {
826  fchi->SetParameter(i, params[i]);
827  }
828  double sigma_up = fchi->Eval(iAvg);
829 
830  // Evaluate sigma low
831  vector<double> defaultChi2low;
832  defaultChi2low.push_back(4.11095e-05);
833  defaultChi2low.push_back(0.577451);
834  defaultChi2low.push_back(-10.378);
835 
836  params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
837  for (unsigned int i = 0; i < params.size(); i++) {
838  fchi->SetParameter(i, params[i]);
839  }
840  double sigma_low = fchi->Eval(iAvg);
841 
842  if (verbose_) {
843  cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
844  }
845 
846  for (int i = 1; i <= iNBins; i++) {
847  if (verbose_) {
848  cout << " " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(i, iBinStrip))
849  << " low: " << fmulow->Eval(iHist->GetBinContent(i, iBinStrip)) << endl;
850  }
851 
852  //evaluate chi2 for cells
853  double muup = fmuup->Eval(iHist->GetBinContent(i, iBinStrip));
854  double mulow = fmulow->Eval(iHist->GetBinContent(i, iBinStrip));
855 
856  //if channel is masked -> set it to value -1
857  if (hservice_->isMasked(iTestName, i, iBinStrip)) {
858  oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
859  }
860  //else perform test
861  else if (muup > sigma_up || mulow > sigma_low ||
862  ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
863  (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
864  dead++;
865  oChannels.push_back(
866  pair<int, double>(iHist->GetBin(i, iBinStrip), abs(iHist->GetBinContent(i, iBinStrip) - iAvg) / iAvg));
867  }
868  }
869  } else {
870  if (verbose_) {
871  cout << "invalid axis" << endl;
872  }
873  }
874 
875  return dead;
876 }
877 
878 //____________________________________________________________________________
879 // Function: getBinCoordinateOnAxisWithValue
880 // Description: Returns the bin global bin number with the iValue in the iAxis
881 // Inputs:
882 // * TH2F* iHist = Histogram to be tested
883 // * double iValue = Value to be evaluated in the histogram iHist
884 // * int& oBinCoordinate = (output) bin number (X or Y) for iValue
885 // * int iAxis = Axis to be used
886 //____________________________________________________________________________
887 void L1TOccupancyClient::getBinCoordinateOnAxisWithValue(TH2F* iHist, double iValue, int& oBinCoordinate, int iAxis) {
888  int nBinsX = iHist->GetNbinsX(); //actual number of bins x
889  int nBinsY = iHist->GetNbinsY(); //actual number of bins y
890 
891  if (iAxis == 1) {
892  int global = iHist->GetXaxis()->FindFixBin(iValue);
893 
894  // If parameter exceeds axis' value: set to maximum number of bins in x-axis
895  if (global > nBinsX * nBinsY) {
896  global = iHist->GetXaxis()->GetLast();
897  }
898 
899  // Get coordinates of bin
900  int y, z;
901  iHist->GetBinXYZ(global, oBinCoordinate, y, z);
902  } else if (iAxis == 2) {
903  int global = iHist->GetYaxis()->FindFixBin(iValue);
904 
905  // If parameter exceeds axis' value: set to maximum number of bins in x-axis
906  if (global > nBinsX * nBinsY) {
907  global = iHist->GetYaxis()->GetLast();
908  }
909 
910  // Get coordinates of bin
911  int x, z;
912  iHist->GetBinXYZ(global, x, oBinCoordinate, z);
913  }
914 }
915 
916 //define this as a plug-in
LuminosityBlockID id() const
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX)
Definition: DQMStore.cc:239
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
virtual void setTitle(const std::string &title)
set (ie. change) histogram/profile title
void getBinCoordinateOnAxisWithValue(TH2F *h2f, double content, int &coord, int axis)
void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:418
void printDeadChannels(const std::vector< std::pair< int, double > > &deadChannels, TH2F *h2f, const std::vector< std::pair< int, double > > &statDev, std::string test_name)
double getAvrg(TH2F *h2f, std::string test, int axis, int nBins, int binStrip, int avrgMode)
double xySymmetry(const edm::ParameterSet &ps, std::string test_name, std::vector< std::pair< int, double > > &deadChannels, std::vector< std::pair< int, double > > &statDev, bool &enoughStats)
void dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) override
virtual void Reset()
reset ME (ie. contents, errors, etc)
void book(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
const double infinity
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
T Max(T a, T b)
Definition: MathUtil.h:44
LuminosityBlockNumber_t luminosityBlock() const
void dqmEndLuminosityBlock(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter, const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &c) override
HLT enums.
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Definition: DQMStore.cc:266
#define str(s)
int compareWithStrip(TH2F *histo, std::string test, int binStrip, int nBins, int axis, double avg, const edm::ParameterSet &ps, std::vector< std::pair< int, double > > &deadChannels)
L1TOccupancyClient(const edm::ParameterSet &ps)
Constructor.
~L1TOccupancyClient() override
Destructor.