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