File indexing completed on 2021-02-14 13:10:50
0001 #include "DQM/L1TMonitorClient/interface/L1TOccupancyClient.h"
0002
0003 #include "FWCore/ServiceRegistry/interface/Service.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include <cstdio>
0010 #include <sstream>
0011 #include <cmath>
0012 #include <vector>
0013 #include <TMath.h>
0014 #include <climits>
0015 #include <TFile.h>
0016 #include <TDirectory.h>
0017 #include <TProfile.h>
0018
0019 using namespace std;
0020 using namespace edm;
0021
0022
0023
0024
0025
0026
0027
0028 L1TOccupancyClient::L1TOccupancyClient(const edm::ParameterSet& ps) {
0029
0030 parameters_ = ps;
0031 verbose_ = ps.getParameter<bool>("verbose");
0032 tests_ = ps.getParameter<std::vector<ParameterSet> >("testParams");
0033
0034 if (verbose_) {
0035 cout << "[L1TOccupancyClient:] Called constructor" << endl;
0036 }
0037 }
0038
0039
0040
0041
0042
0043 L1TOccupancyClient::~L1TOccupancyClient() {
0044 if (verbose_) {
0045 cout << "[L1TOccupancyClient:] Called destructor" << endl;
0046 }
0047 }
0048
0049
0050
0051
0052
0053
0054
0055
0056 void L1TOccupancyClient::book(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter) {
0057 hservice_ = new L1TOccupancyClientHistogramService(parameters_, ibooker, verbose_);
0058
0059 if (verbose_) {
0060 cout << "[L1TOccupancyClient:] Called beginRun" << endl;
0061
0062
0063 file_ = TFile::Open("DQM_L1TOccupancyClient_Snapshots_LS.root", "RECREATE");
0064 }
0065
0066 ibooker.setCurrentFolder("L1T/L1TOccupancy/");
0067
0068
0069
0070
0071
0072 for (vector<ParameterSet>::iterator it = tests_.begin(); it != tests_.end(); it++) {
0073
0074 if ((*it).getUntrackedParameter<string>("algoName", "XYSymmetry") == "XYSymmetry") {
0075
0076 string testName = (*it).getParameter<string>("testName");
0077 ParameterSet algoParameters = (*it).getParameter<ParameterSet>("algoParams");
0078 string histPath = algoParameters.getParameter<string>("histPath");
0079
0080 if (verbose_) {
0081 cout << "[L1TOccupancyClient:] Monitored histogram path: " << histPath << endl;
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092 }
0093
0094
0095 if (hservice_->loadHisto(igetter, testName, histPath)) {
0096
0097 hservice_->setMaskedBins(testName, algoParameters.getParameter<vector<ParameterSet> >("maskedAreas"));
0098
0099
0100
0101 ibooker.setCurrentFolder("L1T/L1TOccupancy/Results");
0102 string title = testName;
0103 MonitorElement* m = ibooker.book2D(title.c_str(), hservice_->getDifferentialHistogram(testName));
0104 m->setTitle(title);
0105 m->Reset();
0106 meResults[title] = m;
0107
0108
0109 ibooker.setCurrentFolder("L1T/L1TOccupancy/HistogramDiff");
0110 title = testName;
0111 m = ibooker.book2D(title.c_str(), hservice_->getDifferentialHistogram(testName));
0112 m->Reset();
0113 m->setTitle(title);
0114 meDifferential[title] = m;
0115
0116
0117 ibooker.setCurrentFolder("L1T/L1TOccupancy/Certification");
0118 title = testName;
0119 m = ibooker.book1D(title.c_str(), title.c_str(), 2500, -.5, 2500. - .5);
0120 m->setTitle(title);
0121 meCertification[title] = m;
0122
0123 mValidTests.push_back(&(*it));
0124 }
0125 }
0126 }
0127 }
0128
0129
0130
0131
0132
0133
0134
0135
0136 void L1TOccupancyClient::dqmEndJob(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter) {
0137 book(ibooker, igetter);
0138
0139 if (verbose_) {
0140 cout << "[L1TOccupancyClient:] Called endRun()" << endl;
0141 }
0142
0143
0144 for (std::vector<ParameterSet*>::iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
0145 ParameterSet& test = (**it);
0146 string algo_name = test.getUntrackedParameter<string>("algoName", "XYSymmetry");
0147 string test_name = test.getParameter<string>("testName");
0148
0149 if (verbose_) {
0150 cout << "[L1TOccupancyClient:] Starting calculations for: " << algo_name << " on: " << test_name << endl;
0151 }
0152
0153 if (algo_name == "XYSymmetry") {
0154 ParameterSet ps = (**it).getParameter<ParameterSet>("algoParams");
0155 string histPath = ps.getParameter<string>("histPath");
0156
0157 vector<pair<int, double> > deadChannels;
0158 vector<pair<int, double> > statDev;
0159 bool enoughStats = false;
0160
0161
0162 hservice_->updateHistogramEndRun(test_name);
0163
0164
0165 double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
0166 stringstream str;
0167 str << test_name << "_cumu_LS_EndRun";
0168
0169 if (verbose_) {
0170 TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0171
0172 cumulative_save->SetTitle(str.str().c_str());
0173
0174 TDirectory* td = file_->GetDirectory(test_name.c_str());
0175
0176 td->cd(string(test_name + "_Histos_AllLS").c_str());
0177
0178 cumulative_save->Write();
0179 }
0180
0181 if (enoughStats) {
0182
0183 printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
0184
0185 if (verbose_) {
0186 TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0187 cumulative_save->SetTitle(str.str().c_str());
0188 TDirectory* td = file_->GetDirectory(("DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
0189 td->cd(string(test_name + "_Histos").c_str());
0190 cumulative_save->Write();
0191
0192
0193 TH2F* h2f = meResults[test_name]->getTH2F();
0194 stringstream str2;
0195 str2 << test_name << "_result_LS_EndRun";
0196 TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
0197
0198 td->cd(string(test_name + "_Results").c_str());
0199 dead_save->SetTitle(str2.str().c_str());
0200 dead_save->Write();
0201 }
0202
0203
0204 meDifferential[test_name]->Reset();
0205 meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
0206
0207 vector<int> lsCertification = hservice_->getLSCertification(test_name);
0208
0209
0210 for (unsigned int i = 0; i < lsCertification.size(); i++) {
0211 int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
0212 meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
0213 }
0214
0215
0216 hservice_->resetHisto(test_name);
0217
0218 if (verbose_) {
0219 cout << "Now we have enough statstics for " << test_name << endl;
0220 }
0221
0222 } else {
0223 if (verbose_) {
0224 cout << "we don't have enough statstics for " << test_name << endl;
0225 }
0226
0227
0228 vector<int> lsCertification = hservice_->getLSCertification(test_name);
0229
0230
0231 for (unsigned int i = 0; i < lsCertification.size(); i++) {
0232 int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
0233 meCertification[test_name]->getTH1()->SetBinContent(bin, -1);
0234 }
0235 }
0236 } else {
0237 if (verbose_) {
0238 cout << "No valid algorithm" << std::endl;
0239 }
0240 }
0241 }
0242
0243 if (verbose_) {
0244 file_->Close();
0245 }
0246
0247 delete hservice_;
0248 }
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265 void L1TOccupancyClient::dqmEndLuminosityBlock(DQMStore::IBooker& ibooker,
0266 DQMStore::IGetter& igetter,
0267 const edm::LuminosityBlock& lumiSeg,
0268 const edm::EventSetup& c) {
0269 book(ibooker, igetter);
0270
0271 int eventLS = lumiSeg.id().luminosityBlock();
0272
0273 if (verbose_) {
0274 cout << "[L1TOccupancyClient:] Called endLuminosityBlock()" << endl;
0275 cout << "[L1TOccupancyClient:] Lumisection: " << eventLS << endl;
0276 }
0277
0278
0279 for (std::vector<ParameterSet*>::const_iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
0280 ParameterSet& test = (**it);
0281 string algo_name = test.getUntrackedParameter<string>("algoName", "XYSymmetry");
0282 string test_name = test.getParameter<string>("testName");
0283
0284 if (verbose_) {
0285 cout << "[L1TOccupancyClient:] Starting calculations for " << algo_name << " on:" << test_name << endl;
0286 }
0287
0288 if (algo_name == "XYSymmetry") {
0289 ParameterSet ps = (**it).getParameter<ParameterSet>("algoParams");
0290 string histPath = ps.getParameter<string>("histPath");
0291
0292 vector<pair<int, double> > deadChannels;
0293 vector<pair<int, double> > statDev;
0294 bool enoughStats = false;
0295
0296
0297 hservice_->updateHistogramEndLS(igetter, test_name, histPath, eventLS);
0298
0299
0300 double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
0301 stringstream str;
0302 str << test_name << "_cumu_LS_" << eventLS;
0303
0304 if (verbose_) {
0305 TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0306 cumulative_save->SetTitle(str.str().c_str());
0307 TDirectory* td = file_->GetDirectory(test_name.c_str());
0308 td->cd(string(test_name + "_Histos_AllLS").c_str());
0309 cumulative_save->Write();
0310 }
0311
0312
0313 if (enoughStats) {
0314
0315 printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
0316
0317 if (verbose_) {
0318 TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0319 cumulative_save->SetTitle(str.str().c_str());
0320 TDirectory* td = file_->GetDirectory(("DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
0321 td->cd(string(test_name + "_Histos").c_str());
0322 cumulative_save->Write();
0323
0324
0325 TH2F* h2f = meResults[test_name]->getTH2F();
0326 stringstream str2;
0327 str2 << test_name << "_result_LS_" << eventLS;
0328 TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
0329
0330 td->cd(string(test_name + "_Results").c_str());
0331 dead_save->SetTitle(str2.str().c_str());
0332 dead_save->Write();
0333 }
0334
0335
0336 meDifferential[test_name]->Reset();
0337 meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
0338
0339 vector<int> lsCertification = hservice_->getLSCertification(test_name);
0340
0341
0342 for (unsigned int i = 0; i < lsCertification.size(); i++) {
0343 int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
0344 meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
0345 }
0346
0347
0348 hservice_->resetHisto(test_name);
0349
0350 if (verbose_) {
0351 cout << "Now we have enough statstics for " << test_name << endl;
0352 }
0353
0354 } else {
0355 if (verbose_) {
0356 cout << "we don't have enough statstics for " << test_name << endl;
0357 }
0358 }
0359 } else {
0360 if (verbose_) {
0361 cout << "No valid algorithm" << std::endl;
0362 }
0363 }
0364 }
0365 }
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388 double L1TOccupancyClient::xySymmetry(const ParameterSet& ps,
0389 string iTestName,
0390 vector<pair<int, double> >& deadChannels,
0391 vector<pair<int, double> >& statDev,
0392 bool& enoughStats) {
0393
0394 TH2F* diffHist = hservice_->getDifferentialHistogram(iTestName);
0395
0396 int pAxis = ps.getUntrackedParameter<int>("axis", 1);
0397 int pAverageMode = ps.getUntrackedParameter<int>("averageMode", 2);
0398 int nBinsX = diffHist->GetNbinsX();
0399 int nBinsY = diffHist->GetNbinsY();
0400
0401
0402 if (pAxis == 1) {
0403 int maxBinStrip, centralBinStrip;
0404
0405 maxBinStrip = nBinsX;
0406
0407
0408
0409 if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
0410 centralBinStrip = nBinsX / 2 + 1;
0411 } else {
0412 double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
0413 getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 1);
0414 }
0415
0416
0417 int upBinStrip = centralBinStrip;
0418 int lowBinStrip = centralBinStrip;
0419
0420
0421 if (nBinsX % 2 == 0) {
0422 lowBinStrip--;
0423 }
0424
0425
0426 std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
0427
0428 int nActualStrips = 0;
0429 for (int i = 0, j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; i++, j++, k--) {
0430 double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsY, j, pAverageMode);
0431 double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsY, k, pAverageMode);
0432
0433
0434 if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
0435 maxAvgs[i] = TMath::Max(avg1, avg2);
0436 nActualStrips++;
0437 }
0438 }
0439
0440 vector<double> defaultMu0up;
0441 defaultMu0up.push_back(13.7655);
0442 defaultMu0up.push_back(184.742);
0443 defaultMu0up.push_back(50735.3);
0444 defaultMu0up.push_back(-97.6793);
0445
0446 TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
0447 vector<double> params = ps.getUntrackedParameter<vector<double> >("params_mu0_up", defaultMu0up);
0448 for (unsigned int i = 0; i < params.size(); i++) {
0449 tf.SetParameter(i, params[i]);
0450 }
0451 int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0452
0453 vector<double> defaultMu0low;
0454 defaultMu0low.push_back(2.19664);
0455 defaultMu0low.push_back(1.94546);
0456 defaultMu0low.push_back(-99.3263);
0457 defaultMu0low.push_back(19.388);
0458
0459 params = ps.getUntrackedParameter<vector<double> >("params_mu0_low", defaultMu0low);
0460 for (unsigned int i = 0; i < params.size(); i++) {
0461 tf.SetParameter(i, params[i]);
0462 }
0463 int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0464
0465 if (verbose_) {
0466 cout << "nbins: " << hservice_->getNBinsHistogram(iTestName) << endl;
0467 cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
0468 }
0469
0470 enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
0471 if (verbose_) {
0472 cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
0473 << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
0474 << ", threshold: " << TMath::Max(statsup, statslow) << endl;
0475 }
0476
0477
0478
0479 if (enoughStats) {
0480 for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
0481 double avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, upBinStrip, pAverageMode);
0482 compareWithStrip(
0483 diffHist, iTestName, lowBinStrip, nBinsY, pAxis, avg, ps, deadChannels);
0484
0485 avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, lowBinStrip, pAverageMode);
0486 compareWithStrip(
0487 diffHist, iTestName, upBinStrip, nBinsY, pAxis, avg, ps, deadChannels);
0488 }
0489 }
0490 }
0491
0492
0493 else if (pAxis == 2) {
0494 int maxBinStrip, centralBinStrip;
0495
0496 maxBinStrip = nBinsY;
0497
0498
0499 if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
0500 centralBinStrip = nBinsY / 2 + 1;
0501 } else {
0502 double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
0503 getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 2);
0504 }
0505
0506
0507 int lowBinStrip = centralBinStrip, upBinStrip = centralBinStrip;
0508
0509
0510 if (nBinsX % 2 == 0) {
0511
0512 lowBinStrip--;
0513 }
0514
0515
0516 std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
0517 int nActualStrips = 0;
0518 for (int i = 0, j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; i++, j++, k--) {
0519 double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsX, j, pAverageMode);
0520 double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsX, k, pAverageMode);
0521 if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
0522 maxAvgs[i] = TMath::Max(avg1, avg2);
0523 nActualStrips++;
0524 }
0525 }
0526
0527 vector<double> defaultMu0up;
0528 defaultMu0up.push_back(13.7655);
0529 defaultMu0up.push_back(184.742);
0530 defaultMu0up.push_back(50735.3);
0531 defaultMu0up.push_back(-97.6793);
0532
0533 vector<double> params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_up", defaultMu0up);
0534 TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
0535 for (unsigned int i = 0; i < params.size(); i++) {
0536 tf.SetParameter(i, params[i]);
0537 }
0538 int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0539
0540 vector<double> defaultMu0low;
0541 defaultMu0low.push_back(2.19664);
0542 defaultMu0low.push_back(1.94546);
0543 defaultMu0low.push_back(-99.3263);
0544 defaultMu0low.push_back(19.388);
0545
0546 params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_low", defaultMu0low);
0547 for (unsigned int i = 0; i < params.size(); i++) {
0548 tf.SetParameter(i, params[i]);
0549 }
0550 int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0551 if (verbose_) {
0552 cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
0553 }
0554 enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
0555 if (verbose_) {
0556 cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
0557 << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
0558 << ", threshold: " << TMath::Max(statsup, statslow) << endl;
0559 }
0560
0561
0562
0563 if (enoughStats) {
0564 for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
0565 double avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, upBinStrip, pAverageMode);
0566 compareWithStrip(
0567 diffHist, iTestName, lowBinStrip, nBinsX, pAxis, avg, ps, deadChannels);
0568
0569 avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, lowBinStrip, pAverageMode);
0570 compareWithStrip(
0571 diffHist, iTestName, upBinStrip, nBinsX, pAxis, avg, ps, deadChannels);
0572 }
0573 }
0574 } else {
0575 if (verbose_) {
0576 cout << "Invalid axis" << endl;
0577 }
0578 }
0579
0580 return (deadChannels.size() - hservice_->getNBinsMasked(iTestName)) * 1.0 / hservice_->getNBinsHistogram(iTestName);
0581 }
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597 double L1TOccupancyClient::getAvrg(TH2F* iHist, string iTestName, int iAxis, int iNBins, int iBinStrip, int iAvgMode) {
0598 double avg = 0.0;
0599 TH1D* proj = nullptr;
0600 TH2F* histo = new TH2F(*iHist);
0601
0602 std::vector<double> values;
0603 int marked;
0604
0605 if (iAxis == 1) {
0606 switch (iAvgMode) {
0607
0608 case 1:
0609 marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0610 proj = histo->ProjectionX();
0611 avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
0612 break;
0613
0614
0615 case 2:
0616 marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0617 proj = histo->ProjectionY("_py", iBinStrip, iBinStrip);
0618 for (int i = 0; i < iNBins; i++) {
0619 values.push_back(proj->GetBinContent(i + 1));
0620 }
0621 avg = TMath::Median(iNBins, &values[0]);
0622 break;
0623 default:
0624 if (verbose_) {
0625 cout << "Invalid averaging mode!" << endl;
0626 }
0627 break;
0628 }
0629 } else if (iAxis == 2) {
0630 switch (iAvgMode) {
0631
0632 case 1:
0633 marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0634 proj = histo->ProjectionY();
0635 avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
0636 break;
0637
0638 case 2:
0639 marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0640 proj = histo->ProjectionX("_px", iBinStrip, iBinStrip);
0641 for (int i = 0; i < iNBins; i++) {
0642 values.push_back(proj->GetBinContent(i + 1));
0643 }
0644
0645 avg = TMath::Median(iNBins, &values[0]);
0646 break;
0647 default:
0648 if (verbose_) {
0649 cout << "invalid averaging mode!" << endl;
0650 }
0651 break;
0652 }
0653 } else {
0654 if (verbose_) {
0655 cout << "invalid axis" << endl;
0656 }
0657 }
0658 delete proj;
0659 delete histo;
0660 return avg;
0661 }
0662
0663
0664
0665
0666
0667
0668
0669
0670
0671
0672 void L1TOccupancyClient::printDeadChannels(const vector<pair<int, double> >& iDeadChannels,
0673 TH2F* oHistDeadChannels,
0674 const vector<std::pair<int, double> >& statDev,
0675 string iTestName) {
0676
0677 oHistDeadChannels->Reset();
0678 if (verbose_) {
0679 cout << "suspect or masked channels of " << iTestName << ": ";
0680 }
0681
0682 int x, y, z;
0683 float chi2 = 0.0;
0684
0685
0686 for (std::vector<pair<int, double> >::const_iterator it = iDeadChannels.begin(); it != iDeadChannels.end(); it++) {
0687 int bin = (*it).first;
0688 oHistDeadChannels->GetBinXYZ(bin, x, y, z);
0689
0690 if (hservice_->isMasked(iTestName, x, y)) {
0691 oHistDeadChannels->SetBinContent(bin, -1);
0692 if (verbose_) {
0693 printf("(%4i,%4i) Masked\n", x, y);
0694 }
0695 } else {
0696 oHistDeadChannels->SetBinContent(bin, 1);
0697 if (verbose_) {
0698 printf("(%4i,%4i) Failed test\n", x, y);
0699 }
0700 }
0701 }
0702
0703
0704 for (std::vector<pair<int, double> >::const_iterator it = statDev.begin(); it != statDev.end(); it++) {
0705 double dev = (*it).second;
0706 chi2 += dev;
0707 }
0708
0709
0710 if (verbose_) {
0711 cout << "total number of suspect channels: " << (iDeadChannels.size() - (hservice_->getNBinsMasked(iTestName)))
0712 << endl;
0713 }
0714 }
0715
0716
0717
0718
0719
0720
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731 int L1TOccupancyClient::compareWithStrip(TH2F* iHist,
0732 string iTestName,
0733 int iBinStrip,
0734 int iNBins,
0735 int iAxis,
0736 double iAvg,
0737 const ParameterSet& iPS,
0738 vector<pair<int, double> >& oChannels) {
0739 int dead = 0;
0740
0741
0742 if (iAxis == 1) {
0743
0744 TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0745 TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0746 fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
0747 fmuup->SetParameter(1, iAvg);
0748 fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
0749 fmulow->SetParameter(1, iAvg);
0750
0751 TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
0752
0753
0754 vector<double> defaultChi2up;
0755 defaultChi2up.push_back(5.45058e-05);
0756 defaultChi2up.push_back(0.268756);
0757 defaultChi2up.push_back(-11.7515);
0758
0759 vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
0760 for (unsigned int i = 0; i < params.size(); i++) {
0761 fchi->SetParameter(i, params[i]);
0762 }
0763 double sigma_up = fchi->Eval(iAvg);
0764
0765
0766 vector<double> defaultChi2low;
0767 defaultChi2low.push_back(4.11095e-05);
0768 defaultChi2low.push_back(0.577451);
0769 defaultChi2low.push_back(-10.378);
0770
0771 params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
0772 for (unsigned int i = 0; i < params.size(); i++) {
0773 fchi->SetParameter(i, params[i]);
0774 }
0775 double sigma_low = fchi->Eval(iAvg);
0776
0777 if (verbose_) {
0778 cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
0779 }
0780
0781 for (int i = 1; i <= iNBins; i++) {
0782 if (verbose_) {
0783 cout << " " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(iBinStrip, i))
0784 << " low: " << fmulow->Eval(iHist->GetBinContent(iBinStrip, i)) << endl;
0785 }
0786
0787
0788 double muup = fmuup->Eval(iHist->GetBinContent(iBinStrip, i));
0789 double mulow = fmulow->Eval(iHist->GetBinContent(iBinStrip, i));
0790
0791
0792 if (hservice_->isMasked(iTestName, iBinStrip, i)) {
0793 oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
0794 }
0795
0796 else if (muup > sigma_up || mulow > sigma_low ||
0797 ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
0798 (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
0799 dead++;
0800 oChannels.push_back(
0801 pair<int, double>(iHist->GetBin(iBinStrip, i), abs(iHist->GetBinContent(iBinStrip, i) - iAvg) / iAvg));
0802 }
0803 }
0804 }
0805
0806 else if (iAxis == 2) {
0807
0808 TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0809 TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0810 fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
0811 fmuup->SetParameter(1, iAvg);
0812 fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
0813 fmulow->SetParameter(1, iAvg);
0814
0815 TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
0816
0817
0818 vector<double> defaultChi2up;
0819 defaultChi2up.push_back(5.45058e-05);
0820 defaultChi2up.push_back(0.268756);
0821 defaultChi2up.push_back(-11.7515);
0822
0823 vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
0824 for (unsigned int i = 0; i < params.size(); i++) {
0825 fchi->SetParameter(i, params[i]);
0826 }
0827 double sigma_up = fchi->Eval(iAvg);
0828
0829
0830 vector<double> defaultChi2low;
0831 defaultChi2low.push_back(4.11095e-05);
0832 defaultChi2low.push_back(0.577451);
0833 defaultChi2low.push_back(-10.378);
0834
0835 params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
0836 for (unsigned int i = 0; i < params.size(); i++) {
0837 fchi->SetParameter(i, params[i]);
0838 }
0839 double sigma_low = fchi->Eval(iAvg);
0840
0841 if (verbose_) {
0842 cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
0843 }
0844
0845 for (int i = 1; i <= iNBins; i++) {
0846 if (verbose_) {
0847 cout << " " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(i, iBinStrip))
0848 << " low: " << fmulow->Eval(iHist->GetBinContent(i, iBinStrip)) << endl;
0849 }
0850
0851
0852 double muup = fmuup->Eval(iHist->GetBinContent(i, iBinStrip));
0853 double mulow = fmulow->Eval(iHist->GetBinContent(i, iBinStrip));
0854
0855
0856 if (hservice_->isMasked(iTestName, i, iBinStrip)) {
0857 oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
0858 }
0859
0860 else if (muup > sigma_up || mulow > sigma_low ||
0861 ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
0862 (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
0863 dead++;
0864 oChannels.push_back(
0865 pair<int, double>(iHist->GetBin(i, iBinStrip), abs(iHist->GetBinContent(i, iBinStrip) - iAvg) / iAvg));
0866 }
0867 }
0868 } else {
0869 if (verbose_) {
0870 cout << "invalid axis" << endl;
0871 }
0872 }
0873
0874 return dead;
0875 }
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886 void L1TOccupancyClient::getBinCoordinateOnAxisWithValue(TH2F* iHist, double iValue, int& oBinCoordinate, int iAxis) {
0887 int nBinsX = iHist->GetNbinsX();
0888 int nBinsY = iHist->GetNbinsY();
0889
0890 if (iAxis == 1) {
0891 int global = iHist->GetXaxis()->FindFixBin(iValue);
0892
0893
0894 if (global > nBinsX * nBinsY) {
0895 global = iHist->GetXaxis()->GetLast();
0896 }
0897
0898
0899 int y, z;
0900 iHist->GetBinXYZ(global, oBinCoordinate, y, z);
0901 } else if (iAxis == 2) {
0902 int global = iHist->GetYaxis()->FindFixBin(iValue);
0903
0904
0905 if (global > nBinsX * nBinsY) {
0906 global = iHist->GetYaxis()->GetLast();
0907 }
0908
0909
0910 int x, z;
0911 iHist->GetBinXYZ(global, x, oBinCoordinate, z);
0912 }
0913 }
0914
0915
0916 DEFINE_FWK_MODULE(L1TOccupancyClient);