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