17 #include <TDirectory.h> 33 tests_ = ps.
getParameter<std::vector<ParameterSet> >(
"testParams");
36 cout <<
"[L1TOccupancyClient:] Called constructor" << endl;
46 cout <<
"[L1TOccupancyClient:] Called destructor" << endl;
61 cout <<
"[L1TOccupancyClient:] Called beginRun" << endl;
64 file_ = TFile::Open(
"DQM_L1TOccupancyClient_Snapshots_LS.root",
"RECREATE");
73 for (vector<ParameterSet>::iterator it = tests_.begin(); it != tests_.end(); it++) {
75 if ((*it).getUntrackedParameter<
string>(
"algoName",
"XYSymmetry") ==
"XYSymmetry") {
77 string testName = (*it).getParameter<
string>(
"testName");
79 string histPath = algoParameters.getParameter<
string>(
"histPath");
82 cout <<
"[L1TOccupancyClient:] Monitored histogram path: " << histPath << endl;
96 if (hservice_->loadHisto(igetter, testName, histPath)) {
98 hservice_->setMaskedBins(testName, algoParameters.getParameter<vector<ParameterSet> >(
"maskedAreas"));
112 m = ibooker.
book2D(title.c_str(), hservice_->getDifferentialHistogram(testName));
115 meDifferential[
title] =
m;
120 m = ibooker.
book1D(title.c_str(), title.c_str(), 2500, -.5, 2500. - .5);
122 meCertification[
title] =
m;
124 mValidTests.push_back(&(*it));
138 book(ibooker, igetter);
141 cout <<
"[L1TOccupancyClient:] Called endRun()" << endl;
145 for (std::vector<ParameterSet*>::iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
148 string test_name = test.
getParameter<
string>(
"testName");
151 cout <<
"[L1TOccupancyClient:] Starting calculations for: " << algo_name <<
" on: " << test_name << endl;
154 if (algo_name ==
"XYSymmetry") {
158 vector<pair<int, double> > deadChannels;
159 vector<pair<int, double> > statDev;
160 bool enoughStats =
false;
163 hservice_->updateHistogramEndRun(test_name);
166 double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
168 str << test_name <<
"_cumu_LS_EndRun";
171 TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
173 cumulative_save->SetTitle(str.str().c_str());
175 TDirectory* td = file_->GetDirectory(test_name.c_str());
177 td->cd(
string(test_name +
"_Histos_AllLS").c_str());
179 cumulative_save->Write();
184 printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
187 TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
188 cumulative_save->SetTitle(str.str().c_str());
189 TDirectory* td = file_->GetDirectory((
"DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
190 td->cd(
string(test_name +
"_Histos").c_str());
191 cumulative_save->Write();
194 TH2F* h2f = meResults[test_name]->getTH2F();
196 str2 << test_name <<
"_result_LS_EndRun";
197 TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
199 td->cd(
string(test_name +
"_Results").c_str());
200 dead_save->SetTitle(str2.str().c_str());
205 meDifferential[test_name]->Reset();
206 meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
208 vector<int> lsCertification = hservice_->getLSCertification(test_name);
211 for (
unsigned int i = 0;
i < lsCertification.size();
i++) {
212 int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[
i]);
213 meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
217 hservice_->resetHisto(test_name);
220 cout <<
"Now we have enough statstics for " << test_name << endl;
225 cout <<
"we don't have enough statstics for " << test_name << endl;
229 vector<int> lsCertification = hservice_->getLSCertification(test_name);
232 for (
unsigned int i = 0;
i < lsCertification.size();
i++) {
233 int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[
i]);
234 meCertification[test_name]->getTH1()->SetBinContent(bin, -1);
239 cout <<
"No valid algorithm" << std::endl;
270 book(ibooker, igetter);
275 cout <<
"[L1TOccupancyClient:] Called endLuminosityBlock()" << endl;
276 cout <<
"[L1TOccupancyClient:] Lumisection: " << eventLS << endl;
280 for (std::vector<ParameterSet*>::const_iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
283 string test_name = test.
getParameter<
string>(
"testName");
286 cout <<
"[L1TOccupancyClient:] Starting calculations for " << algo_name <<
" on:" << test_name << endl;
289 if (algo_name ==
"XYSymmetry") {
293 vector<pair<int, double> > deadChannels;
294 vector<pair<int, double> > statDev;
295 bool enoughStats =
false;
298 hservice_->updateHistogramEndLS(igetter, test_name, histPath, eventLS);
301 double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
303 str << test_name <<
"_cumu_LS_" << eventLS;
306 TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
307 cumulative_save->SetTitle(str.str().c_str());
308 TDirectory* td = file_->GetDirectory(test_name.c_str());
309 td->cd(
string(test_name +
"_Histos_AllLS").c_str());
310 cumulative_save->Write();
316 printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
319 TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
320 cumulative_save->SetTitle(str.str().c_str());
321 TDirectory* td = file_->GetDirectory((
"DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
322 td->cd(
string(test_name +
"_Histos").c_str());
323 cumulative_save->Write();
326 TH2F* h2f = meResults[test_name]->getTH2F();
328 str2 << test_name <<
"_result_LS_" << eventLS;
329 TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
331 td->cd(
string(test_name +
"_Results").c_str());
332 dead_save->SetTitle(str2.str().c_str());
337 meDifferential[test_name]->Reset();
338 meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
340 vector<int> lsCertification = hservice_->getLSCertification(test_name);
343 for (
unsigned int i = 0;
i < lsCertification.size();
i++) {
344 int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[
i]);
345 meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
349 hservice_->resetHisto(test_name);
352 cout <<
"Now we have enough statstics for " << test_name << endl;
357 cout <<
"we don't have enough statstics for " << test_name << endl;
362 cout <<
"No valid algorithm" << std::endl;
391 vector<pair<int, double> >& deadChannels,
392 vector<pair<int, double> >& statDev,
395 TH2F* diffHist = hservice_->getDifferentialHistogram(iTestName);
399 int nBinsX = diffHist->GetNbinsX();
400 int nBinsY = diffHist->GetNbinsY();
404 int maxBinStrip, centralBinStrip;
406 maxBinStrip = nBinsX;
411 centralBinStrip = nBinsX / 2 + 1;
413 double pAxisSymmetryValue = ps.
getParameter<
double>(
"axisSymmetryValue");
414 getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 1);
418 int upBinStrip = centralBinStrip;
419 int lowBinStrip = centralBinStrip;
422 if (nBinsX % 2 == 0) {
427 std::unique_ptr<double[]> maxAvgs(
new double[maxBinStrip - upBinStrip + 1]);
429 int nActualStrips = 0;
430 for (
int i = 0,
j = upBinStrip,
k = lowBinStrip;
j <= maxBinStrip;
i++,
j++,
k--) {
431 double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsY, j, pAverageMode);
432 double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsY,
k, pAverageMode);
435 if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName,
k, pAxis)) {
441 vector<double> defaultMu0up;
442 defaultMu0up.push_back(13.7655);
443 defaultMu0up.push_back(184.742);
444 defaultMu0up.push_back(50735.3);
445 defaultMu0up.push_back(-97.6793);
447 TF1 tf(
"myFunc",
"[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
449 for (
unsigned int i = 0;
i < params.size();
i++) {
450 tf.SetParameter(
i, params[
i]);
452 int statsup = (
int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
454 vector<double> defaultMu0low;
455 defaultMu0low.push_back(2.19664);
456 defaultMu0low.push_back(1.94546);
457 defaultMu0low.push_back(-99.3263);
458 defaultMu0low.push_back(19.388);
461 for (
unsigned int i = 0;
i < params.size();
i++) {
462 tf.SetParameter(
i, params[
i]);
464 int statslow = (
int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
467 cout <<
"nbins: " << hservice_->getNBinsHistogram(iTestName) << endl;
468 cout <<
"statsup= " << statsup <<
", statslow= " << statslow << endl;
471 enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) >
TMath::Max(statsup, statslow);
473 cout <<
"stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
474 <<
", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
475 <<
", threshold: " <<
TMath::Max(statsup, statslow) << endl;
481 for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
482 double avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, upBinStrip, pAverageMode);
484 diffHist, iTestName, lowBinStrip, nBinsY, pAxis, avg, ps, deadChannels);
486 avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, lowBinStrip, pAverageMode);
488 diffHist, iTestName, upBinStrip, nBinsY, pAxis, avg, ps, deadChannels);
494 else if (pAxis == 2) {
495 int maxBinStrip, centralBinStrip;
497 maxBinStrip = nBinsY;
501 centralBinStrip = nBinsY / 2 + 1;
503 double pAxisSymmetryValue = ps.
getParameter<
double>(
"axisSymmetryValue");
504 getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 2);
508 int lowBinStrip = centralBinStrip, upBinStrip = centralBinStrip;
511 if (nBinsX % 2 == 0) {
517 std::unique_ptr<double[]> maxAvgs(
new double[maxBinStrip - upBinStrip + 1]);
518 int nActualStrips = 0;
519 for (
int i = 0,
j = upBinStrip,
k = lowBinStrip;
j <= maxBinStrip;
i++,
j++,
k--) {
520 double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsX, j, pAverageMode);
521 double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsX,
k, pAverageMode);
522 if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName,
k, pAxis)) {
528 vector<double> defaultMu0up;
529 defaultMu0up.push_back(13.7655);
530 defaultMu0up.push_back(184.742);
531 defaultMu0up.push_back(50735.3);
532 defaultMu0up.push_back(-97.6793);
535 TF1 tf(
"myFunc",
"[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
536 for (
unsigned int i = 0;
i < params.size();
i++) {
537 tf.SetParameter(
i, params[
i]);
539 int statsup = (
int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
541 vector<double> defaultMu0low;
542 defaultMu0low.push_back(2.19664);
543 defaultMu0low.push_back(1.94546);
544 defaultMu0low.push_back(-99.3263);
545 defaultMu0low.push_back(19.388);
548 for (
unsigned int i = 0;
i < params.size();
i++) {
549 tf.SetParameter(
i, params[
i]);
551 int statslow = (
int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
553 cout <<
"statsup= " << statsup <<
", statslow= " << statslow << endl;
555 enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) >
TMath::Max(statsup, statslow);
557 cout <<
"stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
558 <<
", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
559 <<
", threshold: " <<
TMath::Max(statsup, statslow) << endl;
565 for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
566 double avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, upBinStrip, pAverageMode);
568 diffHist, iTestName, lowBinStrip, nBinsX, pAxis, avg, ps, deadChannels);
570 avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, lowBinStrip, pAverageMode);
572 diffHist, iTestName, upBinStrip, nBinsX, pAxis, avg, ps, deadChannels);
577 cout <<
"Invalid axis" << endl;
581 return (deadChannels.size() - hservice_->getNBinsMasked(iTestName)) * 1.0 / hservice_->getNBinsHistogram(iTestName);
600 TH1D*
proj =
nullptr;
601 TH2F*
histo =
new TH2F(*iHist);
603 std::vector<double>
values;
610 marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
611 proj = histo->ProjectionX();
612 avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
617 marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
618 proj = histo->ProjectionY(
"_py", iBinStrip, iBinStrip);
619 for (
int i = 0;
i < iNBins;
i++) {
620 values.push_back(proj->GetBinContent(
i + 1));
622 avg = TMath::Median(iNBins, &values[0]);
626 cout <<
"Invalid averaging mode!" << endl;
630 }
else if (iAxis == 2) {
634 marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
635 proj = histo->ProjectionY();
636 avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
640 marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
641 proj = histo->ProjectionX(
"_px", iBinStrip, iBinStrip);
642 for (
int i = 0;
i < iNBins;
i++) {
643 values.push_back(proj->GetBinContent(
i + 1));
646 avg = TMath::Median(iNBins, &values[0]);
650 cout <<
"invalid averaging mode!" << endl;
656 cout <<
"invalid axis" << endl;
674 TH2F* oHistDeadChannels,
675 const vector<std::pair<int, double> >& statDev,
678 oHistDeadChannels->Reset();
680 cout <<
"suspect or masked channels of " << iTestName <<
": ";
687 for (std::vector<pair<int, double> >::const_iterator it = iDeadChannels.begin(); it != iDeadChannels.end(); it++) {
688 int bin = (*it).first;
689 oHistDeadChannels->GetBinXYZ(bin, x, y, z);
691 if (hservice_->isMasked(iTestName, x, y)) {
692 oHistDeadChannels->SetBinContent(bin, -1);
694 printf(
"(%4i,%4i) Masked\n", x, y);
697 oHistDeadChannels->SetBinContent(bin, 1);
699 printf(
"(%4i,%4i) Failed test\n", x, y);
705 for (std::vector<pair<int, double> >::const_iterator it = statDev.begin(); it != statDev.end(); it++) {
706 double dev = (*it).second;
712 cout <<
"total number of suspect channels: " << (iDeadChannels.size() - (hservice_->getNBinsMasked(iTestName)))
739 vector<pair<int, double> >& oChannels) {
745 TF1* fmuup =
new TF1(
"fmuup",
"TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
746 TF1* fmulow =
new TF1(
"fmulow",
"TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
748 fmuup->SetParameter(1, iAvg);
750 fmulow->SetParameter(1, iAvg);
752 TF1* fchi =
new TF1(
"fchi",
"[0]*x**2+[1]*x+[2]", 0., 1500.);
755 vector<double> defaultChi2up;
756 defaultChi2up.push_back(5.45058
e-05);
757 defaultChi2up.push_back(0.268756);
758 defaultChi2up.push_back(-11.7515);
761 for (
unsigned int i = 0;
i < params.size();
i++) {
762 fchi->SetParameter(
i, params[
i]);
764 double sigma_up = fchi->Eval(iAvg);
767 vector<double> defaultChi2low;
768 defaultChi2low.push_back(4.11095
e-05);
769 defaultChi2low.push_back(0.577451);
770 defaultChi2low.push_back(-10.378);
773 for (
unsigned int i = 0;
i < params.size();
i++) {
774 fchi->SetParameter(
i, params[
i]);
776 double sigma_low = fchi->Eval(iAvg);
779 cout <<
"binstrip= " << iBinStrip <<
", sigmaup= " << sigma_up <<
", sigmalow= " << sigma_low << endl;
782 for (
int i = 1;
i <= iNBins;
i++) {
784 cout <<
" " <<
i <<
" binContent: up:" << fmuup->Eval(iHist->GetBinContent(iBinStrip,
i))
785 <<
" low: " << fmulow->Eval(iHist->GetBinContent(iBinStrip,
i)) << endl;
789 double muup = fmuup->Eval(iHist->GetBinContent(iBinStrip,
i));
790 double mulow = fmulow->Eval(iHist->GetBinContent(iBinStrip,
i));
793 if (hservice_->isMasked(iTestName, iBinStrip,
i)) {
794 oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip,
i), -1.0));
797 else if (muup > sigma_up || mulow > sigma_low ||
802 pair<int, double>(iHist->GetBin(iBinStrip,
i),
abs(iHist->GetBinContent(iBinStrip,
i) - iAvg) / iAvg));
807 else if (iAxis == 2) {
809 TF1* fmuup =
new TF1(
"fmuup",
"TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
810 TF1* fmulow =
new TF1(
"fmulow",
"TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
812 fmuup->SetParameter(1, iAvg);
814 fmulow->SetParameter(1, iAvg);
816 TF1* fchi =
new TF1(
"fchi",
"[0]*x**2+[1]*x+[2]", 0., 1500.);
819 vector<double> defaultChi2up;
820 defaultChi2up.push_back(5.45058
e-05);
821 defaultChi2up.push_back(0.268756);
822 defaultChi2up.push_back(-11.7515);
825 for (
unsigned int i = 0;
i < params.size();
i++) {
826 fchi->SetParameter(
i, params[
i]);
828 double sigma_up = fchi->Eval(iAvg);
831 vector<double> defaultChi2low;
832 defaultChi2low.push_back(4.11095
e-05);
833 defaultChi2low.push_back(0.577451);
834 defaultChi2low.push_back(-10.378);
837 for (
unsigned int i = 0;
i < params.size();
i++) {
838 fchi->SetParameter(
i, params[
i]);
840 double sigma_low = fchi->Eval(iAvg);
843 cout <<
"binstrip= " << iBinStrip <<
", sigmaup= " << sigma_up <<
", sigmalow= " << sigma_low << endl;
846 for (
int i = 1;
i <= iNBins;
i++) {
848 cout <<
" " <<
i <<
" binContent: up:" << fmuup->Eval(iHist->GetBinContent(
i, iBinStrip))
849 <<
" low: " << fmulow->Eval(iHist->GetBinContent(
i, iBinStrip)) << endl;
853 double muup = fmuup->Eval(iHist->GetBinContent(
i, iBinStrip));
854 double mulow = fmulow->Eval(iHist->GetBinContent(
i, iBinStrip));
857 if (hservice_->isMasked(iTestName,
i, iBinStrip)) {
858 oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip,
i), -1.0));
861 else if (muup > sigma_up || mulow > sigma_low ||
866 pair<int, double>(iHist->GetBin(
i, iBinStrip),
abs(iHist->GetBinContent(
i, iBinStrip) - iAvg) / iAvg));
871 cout <<
"invalid axis" << endl;
888 int nBinsX = iHist->GetNbinsX();
889 int nBinsY = iHist->GetNbinsY();
892 int global = iHist->GetXaxis()->FindFixBin(iValue);
895 if (global > nBinsX * nBinsY) {
896 global = iHist->GetXaxis()->GetLast();
901 iHist->GetBinXYZ(global, oBinCoordinate, y, z);
902 }
else if (iAxis == 2) {
903 int global = iHist->GetYaxis()->FindFixBin(iValue);
906 if (global > nBinsX * nBinsY) {
907 global = iHist->GetYaxis()->GetLast();
912 iHist->GetBinXYZ(global, x, oBinCoordinate, z);
LuminosityBlockID id() const
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX)
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
virtual void setTitle(const std::string &title)
set (ie. change) histogram/profile title
void getBinCoordinateOnAxisWithValue(TH2F *h2f, double content, int &coord, int axis)
void setCurrentFolder(std::string const &fullpath)
void printDeadChannels(const std::vector< std::pair< int, double > > &deadChannels, TH2F *h2f, const std::vector< std::pair< int, double > > &statDev, std::string test_name)
double getAvrg(TH2F *h2f, std::string test, int axis, int nBins, int binStrip, int avrgMode)
double xySymmetry(const edm::ParameterSet &ps, std::string test_name, std::vector< std::pair< int, double > > &deadChannels, std::vector< std::pair< int, double > > &statDev, bool &enoughStats)
void dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) override
virtual void Reset()
reset ME (ie. contents, errors, etc)
void book(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter)
#define DEFINE_FWK_MODULE(type)
Abs< T >::type abs(const T &t)
LuminosityBlockNumber_t luminosityBlock() const
void dqmEndLuminosityBlock(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter, const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &c) override
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
int compareWithStrip(TH2F *histo, std::string test, int binStrip, int nBins, int axis, double avg, const edm::ParameterSet &ps, std::vector< std::pair< int, double > > &deadChannels)
L1TOccupancyClient(const edm::ParameterSet &ps)
Constructor.
~L1TOccupancyClient() override
Destructor.