Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-11 04:32:56 

0001 #include "DQMServices/Core/interface/QTest.h"
0002 #include "DQMServices/Core/interface/DQMStore.h"
0003 #include "Math/ProbFuncMathCore.h"
0004 #include "TMath.h"
0005 #include <cmath>
0006 #include <iostream>
0007 #include <sstream>
0008 
0009 using namespace std;
0010 
0011 const float QCriterion::ERROR_PROB_THRESHOLD = 0.50;
0012 const float QCriterion::WARNING_PROB_THRESHOLD = 0.90;
0013 
0014 // initialize values
0015 void QCriterion::init() {
0016   errorProb_ = ERROR_PROB_THRESHOLD;
0017   warningProb_ = WARNING_PROB_THRESHOLD;
0018   setAlgoName("NO_ALGORITHM");
0019   status_ = dqm::qstatus::DID_NOT_RUN;
0020   message_ = "NO_MESSAGE";
0021   verbose_ = 0;  // 0 = silent, 1 = algorithmic failures, 2 = info
0022 }
0023 
0024 float QCriterion::runTest(const MonitorElement* /* me */) { return 0.; }
0025 //===================================================//
0026 //================ QUALITY TESTS ====================//
0027 //==================================================//
0028 
0029 //----------------------------------------------------//
0030 //--------------- ContentsXRange ---------------------//
0031 //----------------------------------------------------//
0032 float ContentsXRange::runTest(const MonitorElement* me) {
0033   badChannels_.clear();
0034 
0035   if (!me)
0036     return -1;
0037   if (!me->getRootObject())
0038     return -1;
0039   TH1* h = nullptr;
0040 
0041   if (verbose_ > 1)
0042     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
0043   // -- TH1F
0044   if (me->kind() == MonitorElement::Kind::TH1F) {
0045     h = me->getTH1F();
0046   }
0047   // -- TH1S
0048   else if (me->kind() == MonitorElement::Kind::TH1S) {
0049     h = me->getTH1S();
0050   }
0051   // -- TH1D
0052   else if (me->kind() == MonitorElement::Kind::TH1D) {
0053     h = me->getTH1D();
0054   } else {
0055     if (verbose_ > 0)
0056       std::cout << "QTest:ContentsXRange"
0057                 << " ME " << me->getFullname() << " does not contain TH1F/TH1S/TH1D, exiting\n";
0058     return -1;
0059   }
0060 
0061   if (!rangeInitialized_) {
0062     if (h->GetXaxis())
0063       setAllowedXRange(h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
0064     else
0065       return -1;
0066   }
0067   int ncx = h->GetXaxis()->GetNbins();
0068   // use underflow bin
0069   int first = 0;  // 1
0070   // use overflow bin
0071   int last = ncx + 1;  // ncx
0072   // all entries
0073   double sum = 0;
0074   // entries outside X-range
0075   double fail = 0;
0076   int bin;
0077   for (bin = first; bin <= last; ++bin) {
0078     double contents = h->GetBinContent(bin);
0079     double x = h->GetBinCenter(bin);
0080     sum += contents;
0081     if (x < xmin_ || x > xmax_)
0082       fail += contents;
0083   }
0084 
0085   if (sum == 0)
0086     return 1;
0087   // return fraction of entries within allowed X-range
0088   return (sum - fail) / sum;
0089 }
0090 
0091 //-----------------------------------------------------//
0092 //--------------- ContentsYRange ---------------------//
0093 //----------------------------------------------------//
0094 float ContentsYRange::runTest(const MonitorElement* me) {
0095   badChannels_.clear();
0096 
0097   if (!me)
0098     return -1;
0099   if (!me->getRootObject())
0100     return -1;
0101   TH1* h = nullptr;
0102 
0103   if (verbose_ > 1)
0104     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
0105 
0106   if (me->kind() == MonitorElement::Kind::TH1F) {
0107     h = me->getTH1F();  //access Test histo
0108   } else if (me->kind() == MonitorElement::Kind::TH1S) {
0109     h = me->getTH1S();  //access Test histo
0110   } else if (me->kind() == MonitorElement::Kind::TH1D) {
0111     h = me->getTH1D();  //access Test histo
0112   } else {
0113     if (verbose_ > 0)
0114       std::cout << "QTest:ContentsYRange"
0115                 << " ME " << me->getFullname() << " does not contain TH1F/TH1S/TH1D, exiting\n";
0116     return -1;
0117   }
0118 
0119   if (!rangeInitialized_ || !h->GetXaxis())
0120     return 1;  // all bins are accepted if no initialization
0121   int ncx = h->GetXaxis()->GetNbins();
0122   // do NOT use underflow bin
0123   int first = 1;
0124   // do NOT use overflow bin
0125   int last = ncx;
0126   // bins outside Y-range
0127   int fail = 0;
0128   int bin;
0129 
0130   if (useEmptyBins_)  ///Standard test !
0131   {
0132     for (bin = first; bin <= last; ++bin) {
0133       double contents = h->GetBinContent(bin);
0134       bool failure = false;
0135       failure = (contents < ymin_ || contents > ymax_);  // allowed y-range: [ymin_, ymax_]
0136       if (failure) {
0137         DQMChannel chan(bin, 0, 0, contents, h->GetBinError(bin));
0138         badChannels_.push_back(chan);
0139         ++fail;
0140       }
0141     }
0142     // return fraction of bins that passed test
0143     return 1. * (ncx - fail) / ncx;
0144   } else  ///AS quality test !!!
0145   {
0146     for (bin = first; bin <= last; ++bin) {
0147       double contents = h->GetBinContent(bin);
0148       bool failure = false;
0149       if (contents)
0150         failure = (contents < ymin_ || contents > ymax_);  // allowed y-range: [ymin_, ymax_]
0151       if (failure)
0152         ++fail;
0153     }
0154     // return fraction of bins that passed test
0155     return 1. * (ncx - fail) / ncx;
0156   }  ///end of AS quality tests
0157 }
0158 
0159 //-----------------------------------------------------//
0160 //------------------ DeadChannel ---------------------//
0161 //----------------------------------------------------//
0162 float DeadChannel::runTest(const MonitorElement* me) {
0163   badChannels_.clear();
0164   if (!me)
0165     return -1;
0166   if (!me->getRootObject())
0167     return -1;
0168   TH1* h1 = nullptr;
0169   TH2* h2 = nullptr;  //initialize histogram pointers
0170 
0171   if (verbose_ > 1)
0172     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
0173   //TH1F
0174   if (me->kind() == MonitorElement::Kind::TH1F) {
0175     h1 = me->getTH1F();  //access Test histo
0176   }
0177   //TH1S
0178   else if (me->kind() == MonitorElement::Kind::TH1S) {
0179     h1 = me->getTH1S();  //access Test histo
0180   }
0181   //TH1D
0182   else if (me->kind() == MonitorElement::Kind::TH1D) {
0183     h1 = me->getTH1D();  //access Test histo
0184   }
0185   //-- TH2F
0186   else if (me->kind() == MonitorElement::Kind::TH2F) {
0187     h2 = me->getTH2F();  // access Test histo
0188   }
0189   //-- TH2S
0190   else if (me->kind() == MonitorElement::Kind::TH2S) {
0191     h2 = me->getTH2S();  // access Test histo
0192   }
0193   //-- TH2D
0194   else if (me->kind() == MonitorElement::Kind::TH2D) {
0195     h2 = me->getTH2D();  // access Test histo
0196   } else {
0197     if (verbose_ > 0)
0198       std::cout << "QTest:DeadChannel"
0199                 << " ME " << me->getFullname() << " does not contain TH1F/TH1S/TH1D/TH2F/TH2S/TH2D, exiting\n";
0200     return -1;
0201   }
0202 
0203   int fail = 0;  // number of failed channels
0204 
0205   //--------- do the quality test for 1D histo ---------------//
0206   if (h1 != nullptr) {
0207     if (!rangeInitialized_ || !h1->GetXaxis())
0208       return 1;  // all bins are accepted if no initialization
0209     int ncx = h1->GetXaxis()->GetNbins();
0210     int first = 1;
0211     int last = ncx;
0212     int bin;
0213 
0214     /// loop over all channels
0215     for (bin = first; bin <= last; ++bin) {
0216       double contents = h1->GetBinContent(bin);
0217       bool failure = false;
0218       failure = contents <= ymin_;  // dead channel: equal to or less than ymin_
0219       if (failure) {
0220         DQMChannel chan(bin, 0, 0, contents, h1->GetBinError(bin));
0221         badChannels_.push_back(chan);
0222         ++fail;
0223       }
0224     }
0225     //return fraction of alive channels
0226     return 1. * (ncx - fail) / ncx;
0227   }
0228   //----------------------------------------------------------//
0229 
0230   //--------- do the quality test for 2D -------------------//
0231   else if (h2 != nullptr) {
0232     int ncx = h2->GetXaxis()->GetNbins();  // get X bins
0233     int ncy = h2->GetYaxis()->GetNbins();  // get Y bins
0234 
0235     /// loop over all bins
0236     for (int cx = 1; cx <= ncx; ++cx) {
0237       for (int cy = 1; cy <= ncy; ++cy) {
0238         double contents = h2->GetBinContent(h2->GetBin(cx, cy));
0239         bool failure = false;
0240         failure = contents <= ymin_;  // dead channel: equal to or less than ymin_
0241         if (failure) {
0242           DQMChannel chan(cx, cy, 0, contents, h2->GetBinError(h2->GetBin(cx, cy)));
0243           badChannels_.push_back(chan);
0244           ++fail;
0245         }
0246       }
0247     }
0248     //return fraction of alive channels
0249     return 1. * (ncx * ncy - fail) / (ncx * ncy);
0250   } else {
0251     if (verbose_ > 0)
0252       std::cout << "QTest:DeadChannel"
0253                 << " TH1/TH2F are NULL, exiting\n";
0254     return -1;
0255   }
0256 }
0257 
0258 //-----------------------------------------------------//
0259 //----------------  NoisyChannel ---------------------//
0260 //----------------------------------------------------//
0261 // run the test (result: fraction of channels not appearing noisy or "hot")
0262 // [0, 1] or <0 for failure
0263 float NoisyChannel::runTest(const MonitorElement* me) {
0264   badChannels_.clear();
0265   if (!me)
0266     return -1;
0267   if (!me->getRootObject())
0268     return -1;
0269   TH1* h = nullptr;   //initialize histogram pointer
0270   TH2* h2 = nullptr;  //initialize histogram pointer
0271 
0272   if (verbose_ > 1)
0273     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
0274 
0275   int nbins = 0;
0276   int nbinsX = 0, nbinsY = 0;
0277   //-- TH1F
0278   if (me->kind() == MonitorElement::Kind::TH1F) {
0279     nbins = me->getTH1F()->GetXaxis()->GetNbins();
0280     h = me->getTH1F();  // access Test histo
0281   }
0282   //-- TH1S
0283   else if (me->kind() == MonitorElement::Kind::TH1S) {
0284     nbins = me->getTH1S()->GetXaxis()->GetNbins();
0285     h = me->getTH1S();  // access Test histo
0286   }
0287   //-- TH1D
0288   else if (me->kind() == MonitorElement::Kind::TH1D) {
0289     nbins = me->getTH1D()->GetXaxis()->GetNbins();
0290     h = me->getTH1D();  // access Test histo
0291   }
0292   //-- TH2
0293   else if (me->kind() == MonitorElement::Kind::TH2F) {
0294     nbinsX = me->getTH2F()->GetXaxis()->GetNbins();
0295     nbinsY = me->getTH2F()->GetYaxis()->GetNbins();
0296     h2 = me->getTH2F();  // access Test histo
0297   }
0298   //-- TH2
0299   else if (me->kind() == MonitorElement::Kind::TH2S) {
0300     nbinsX = me->getTH2S()->GetXaxis()->GetNbins();
0301     nbinsY = me->getTH2S()->GetYaxis()->GetNbins();
0302     h2 = me->getTH2S();  // access Test histo
0303   }
0304   //-- TH2
0305   else if (me->kind() == MonitorElement::Kind::TH2D) {
0306     nbinsX = me->getTH2F()->GetXaxis()->GetNbins();
0307     nbinsY = me->getTH2F()->GetYaxis()->GetNbins();
0308     h2 = me->getTH2D();  // access Test histo
0309   } else {
0310     if (verbose_ > 0)
0311       std::cout << "QTest:NoisyChannel"
0312                 << " ME " << me->getFullname() << " does not contain TH1F/TH1S/TH1D or TH2F/TH2S/TH2D, exiting\n";
0313     return -1;
0314   }
0315 
0316   //--  QUALITY TEST itself
0317 
0318   // do NOT use underflow bin
0319   int first = 1;
0320   // do NOT use overflow bin
0321   int last = nbins;
0322   int lastX = nbinsX, lastY = nbinsY;
0323   // bins outside Y-range
0324   int fail = 0;
0325   int bin;
0326   int binX, binY;
0327   if (h != nullptr) {
0328     if (!rangeInitialized_ || !h->GetXaxis()) {
0329       return 1;  // all channels are accepted if tolerance has not been set
0330     }
0331     for (bin = first; bin <= last; ++bin) {
0332       double contents = h->GetBinContent(bin);
0333       double average = getAverage(bin, h);
0334       bool failure = false;
0335       if (average != 0)
0336         failure = (((contents - average) / std::abs(average)) > tolerance_);
0337 
0338       if (failure) {
0339         ++fail;
0340         DQMChannel chan(bin, 0, 0, contents, h->GetBinError(bin));
0341         badChannels_.push_back(chan);
0342       }
0343     }
0344 
0345     // return fraction of bins that passed test
0346     return 1. * (nbins - fail) / nbins;
0347   } else if (h2 != nullptr) {
0348     for (binY = first; binY <= lastY; ++binY) {
0349       for (binX = first; binX <= lastX; ++binX) {
0350         double contents = h2->GetBinContent(binX, binY);
0351         double average = getAverage2D(binX, binY, h2);
0352         bool failure = false;
0353         if (average != 0)
0354           failure = (((contents - average) / std::abs(average)) > tolerance_);
0355         if (failure) {
0356           ++fail;
0357           DQMChannel chan(binX, 0, 0, contents, h2->GetBinError(binX));
0358           badChannels_.push_back(chan);
0359         }
0360       }  //end x loop
0361     }  //end y loop
0362     // return fraction of bins that passed test
0363     return 1. * ((nbinsX * nbinsY) - fail) / (nbinsX * nbinsY);
0364   }  //end nullptr conditional
0365   else {
0366     if (verbose_ > 0)
0367       std::cout << "QTest:NoisyChannel"
0368                 << " TH1/TH2F are NULL, exiting\n";
0369     return -1;
0370   }
0371 }
0372 
0373 // get average for bin under consideration
0374 // (see description of method setNumNeighbors)
0375 double NoisyChannel::getAverage(int bin, const TH1* h) const {
0376   int first = 1;                        // Do NOT use underflow bin
0377   int ncx = h->GetXaxis()->GetNbins();  // Do NOT use overflow bin
0378   double sum = 0;
0379   int bin_start, bin_end;
0380   int add_right = 0;
0381   int add_left = 0;
0382 
0383   bin_start = bin - numNeighbors_;  // First bin in integral
0384   bin_end = bin + numNeighbors_;    // Last bin in integral
0385 
0386   if (bin_start < first) {          // If neighbors take you outside of histogram range shift integral right
0387     add_right = first - bin_start;  // How much to shift remembering we are not using underflow
0388     bin_start = first;              // Remember to reset the starting bin
0389     bin_end += add_right;
0390     if (bin_end > ncx)
0391       bin_end = ncx;  // If the test would be larger than histogram just sum histogram without overflow
0392   }
0393 
0394   if (bin_end > ncx) {  // Now we make sure doesn't run off right edge of histogram
0395     add_left = bin_end - ncx;
0396     bin_end = ncx;
0397     bin_start -= add_left;
0398     if (bin_start < first)
0399       bin_start = first;  // If the test would be larger than histogram just sum histogram without underflow
0400   }
0401 
0402   sum += h->Integral(bin_start, bin_end);
0403   sum -= h->GetBinContent(bin);
0404 
0405   int dimension = 2 * numNeighbors_ + 1;
0406   if (dimension > h->GetNbinsX())
0407     dimension = h->GetNbinsX();
0408 
0409   return sum / (dimension - 1);
0410 }
0411 
0412 double NoisyChannel::getAverage2D(int binX, int binY, const TH2* h2) const {
0413   /// Do NOT use underflow or overflow bins
0414   int firstX = 1;
0415   int firstY = 1;
0416   double sum = 0;
0417   int ncx = h2->GetXaxis()->GetNbins();
0418   int ncy = h2->GetYaxis()->GetNbins();
0419 
0420   int neighborsX, neighborsY;  // Convert unsigned input to int so we can use comparators
0421   neighborsX = numNeighbors_;
0422   neighborsY = numNeighbors_;
0423   int bin_startX, bin_endX;
0424   int add_rightX = 0;  // Start shifts at 0
0425   int add_leftX = 0;
0426   int bin_startY, bin_endY;
0427   int add_topY = 0;
0428   int add_downY = 0;
0429 
0430   bin_startX = binX - neighborsX;  // First bin in X
0431   bin_endX = binX + neighborsX;    // Last bin in X
0432 
0433   if (bin_startX < firstX) {           // If neighbors take you outside of histogram range shift integral right
0434     add_rightX = firstX - bin_startX;  // How much to shift remembering we are no using underflow
0435     bin_startX = firstX;               // Remember to reset the starting bin
0436     bin_endX += add_rightX;
0437     if (bin_endX > ncx)
0438       bin_endX = ncx;
0439   }
0440 
0441   if (bin_endX > ncx) {  // Now we make sure doesn't run off right edge of histogram
0442     add_leftX = bin_endX - ncx;
0443     bin_endX = ncx;
0444     bin_startX -= add_leftX;
0445     if (bin_startX < firstX)
0446       bin_startX = firstX;  // If the test would be larger than histogram just sum histogram without underflow
0447   }
0448 
0449   bin_startY = binY - neighborsY;  // First bin in Y
0450   bin_endY = binY + neighborsY;    // Last bin in Y
0451 
0452   if (bin_startY < firstY) {         // If neighbors take you outside of histogram range shift integral up
0453     add_topY = firstY - bin_startY;  // How much to shift remembering we are no using underflow
0454     bin_startY = firstY;             // Remember to reset the starting bin
0455     bin_endY += add_topY;
0456     if (bin_endY > ncy)
0457       bin_endY = ncy;
0458   }
0459 
0460   if (bin_endY > ncy) {  // Now we make sure doesn't run off top edge of histogram
0461     add_downY = bin_endY - ncy;
0462     bin_endY = ncy;
0463     bin_startY -= add_downY;
0464     if (bin_startY < firstY)
0465       bin_startY = firstY;  // If the test would be larger than histogram just sum histogram without underflow
0466   }
0467 
0468   sum += h2->Integral(bin_startX, bin_endX, bin_startY, bin_endY);
0469   sum -= h2->GetBinContent(binX, binY);
0470 
0471   int dimensionX = 2 * neighborsX + 1;
0472   int dimensionY = 2 * neighborsY + 1;
0473 
0474   if (dimensionX > h2->GetNbinsX())
0475     dimensionX = h2->GetNbinsX();
0476   if (dimensionY > h2->GetNbinsY())
0477     dimensionY = h2->GetNbinsY();
0478 
0479   return sum / (dimensionX * dimensionY - 1);  // Average is sum over the # of bins used
0480 
0481 }  // End getAverage2D
0482 
0483 //-----------------------------------------------------//
0484 //-----Content Sigma (Emma Yeager and Chad Freer)------//
0485 //----------------------------------------------------//
0486 // run the test (result: fraction of channels with sigma that is not noisy or hot)
0487 
0488 float ContentSigma::runTest(const MonitorElement* me) {
0489   badChannels_.clear();
0490   if (!me)
0491     return -1;
0492   if (!me->getRootObject())
0493     return -1;
0494   TH1* h = nullptr;  //initialize histogram pointer
0495 
0496   if (verbose_ > 1)
0497     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
0498 
0499   unsigned nbinsX;
0500   unsigned nbinsY;
0501 
0502   //-- TH1F
0503   if (me->kind() == MonitorElement::Kind::TH1F) {
0504     nbinsX = me->getTH1F()->GetXaxis()->GetNbins();
0505     nbinsY = me->getTH1F()->GetYaxis()->GetNbins();
0506     h = me->getTH1F();  // access Test histo
0507   }
0508   //-- TH1S
0509   else if (me->kind() == MonitorElement::Kind::TH1S) {
0510     nbinsX = me->getTH1S()->GetXaxis()->GetNbins();
0511     nbinsY = me->getTH1S()->GetYaxis()->GetNbins();
0512     h = me->getTH1S();  // access Test histo
0513   }
0514   //-- TH1D
0515   else if (me->kind() == MonitorElement::Kind::TH1D) {
0516     nbinsX = me->getTH1D()->GetXaxis()->GetNbins();
0517     nbinsY = me->getTH1D()->GetYaxis()->GetNbins();
0518     h = me->getTH1D();  // access Test histo
0519   }
0520   //-- TH2
0521   else if (me->kind() == MonitorElement::Kind::TH2F) {
0522     nbinsX = me->getTH2F()->GetXaxis()->GetNbins();
0523     nbinsY = me->getTH2F()->GetYaxis()->GetNbins();
0524     h = me->getTH2F();  // access Test histo
0525   }
0526   //-- TH2
0527   else if (me->kind() == MonitorElement::Kind::TH2S) {
0528     nbinsX = me->getTH2S()->GetXaxis()->GetNbins();
0529     nbinsY = me->getTH2S()->GetYaxis()->GetNbins();
0530     h = me->getTH2S();  // access Test histo
0531   }
0532   //-- TH2
0533   else if (me->kind() == MonitorElement::Kind::TH2D) {
0534     nbinsX = me->getTH2D()->GetXaxis()->GetNbins();
0535     nbinsY = me->getTH2D()->GetYaxis()->GetNbins();
0536     h = me->getTH2D();  // access Test histo
0537   } else {
0538     if (verbose_ > 0)
0539       std::cout << "QTest:ContentSigma"
0540                 << " ME " << me->getFullname() << " does not contain TH1F/TH1S/TH1D or TH2F/TH2S/TH2D, exiting\n";
0541     return -1;
0542   }
0543 
0544   //--  QUALITY TEST itself
0545 
0546   if (!rangeInitialized_ || !h->GetXaxis())
0547     return 1;  // all channels are accepted if tolerance has not been set
0548 
0549   int fail = 0;       // initialize bin failure count
0550   unsigned xMin = 1;  //initialize minimums and maximums with expected values
0551   unsigned yMin = 1;
0552   unsigned xMax = nbinsX;
0553   unsigned yMax = nbinsY;
0554   unsigned XBlocks = numXblocks_;  //Initialize xml inputs blocks and neighbors
0555   unsigned YBlocks = numYblocks_;
0556   unsigned neighborsX = numNeighborsX_;
0557   unsigned neighborsY = numNeighborsY_;
0558   unsigned Xbinnum = 1;
0559   unsigned Ybinnum = 1;
0560   unsigned XWidth = 0;
0561   unsigned YWidth = 0;
0562 
0563   if (neighborsX == 999) {
0564     neighborsX = 0;
0565   }
0566   if (neighborsY == 999) {
0567     neighborsY = 0;
0568   }
0569 
0570   //give users option for automatic mininum and maximum selection by inputting 0 to any of the parameters
0571   // check that user's parameters are completely in agreement with histogram
0572   // for instance, if inputted xMax is out of range xMin will automatically be ignored
0573   if (xMin_ != 0 && xMax_ != 0) {
0574     if ((xMax_ <= nbinsX) && (xMin_ <= xMax_)) {  // rescale area of histogram being analyzed
0575       nbinsX = xMax_ - xMin_ + 1;
0576       xMax = xMax_;  // do NOT use overflow bin
0577       xMin = xMin_;  // do NOT use underflow bin
0578     }
0579   }
0580   //give users option for automatic mininum and maximum selection by inputting 0 to any of the parameters
0581   if (yMin_ != 0 && yMax_ != 0) {
0582     if ((yMax_ <= nbinsY) && (yMin_ <= yMax_)) {
0583       nbinsY = yMax_ - yMin_ + 1;
0584       yMax = yMax_;
0585       yMin = yMin_;
0586     }
0587   }
0588 
0589   if (neighborsX * 2 >= nbinsX) {  //make sure neighbor check does not overlap with bin under consideration
0590     if (nbinsX % 2 == 0) {
0591       neighborsX = nbinsX / 2 - 1;  //set neighbors for no overlap
0592     } else {
0593       neighborsX = (nbinsX - 1) / 2;
0594     }
0595   }
0596 
0597   if (neighborsY * 2 >= nbinsY) {
0598     if (nbinsY % 2 == 0) {
0599       neighborsY = nbinsY / 2 - 1;
0600     } else {
0601       neighborsY = (nbinsY - 1) / 2;
0602     }
0603   }
0604 
0605   if (XBlocks == 999) {  //Setting 999 prevents blocks and does quality tests by bins only
0606     XBlocks = nbinsX;
0607   }
0608   if (YBlocks == 999) {
0609     YBlocks = nbinsY;
0610   }
0611 
0612   Xbinnum = nbinsX / XBlocks;
0613   Ybinnum = nbinsY / YBlocks;
0614   for (unsigned groupx = 0; groupx < XBlocks; ++groupx) {  //Test over all the blocks
0615     for (unsigned groupy = 0; groupy < YBlocks; ++groupy) {
0616       double blocksum = 0;
0617       for (unsigned binx = 0; binx < Xbinnum; ++binx) {  //Sum the contents of the block in question
0618         for (unsigned biny = 0; biny < Ybinnum; ++biny) {
0619           if (groupx * Xbinnum + xMin + binx <= xMax && groupy * Ybinnum + yMin + biny <= yMax) {
0620             blocksum += abs(h->GetBinContent(groupx * Xbinnum + xMin + binx, groupy * Ybinnum + yMin + biny));
0621           }
0622         }
0623       }
0624 
0625       double sum = getNeighborSum(groupx, groupy, XBlocks, YBlocks, neighborsX, neighborsY, h);
0626       sum -= blocksum;  //remove center block to test
0627 
0628       if (neighborsX > groupx) {  //Find correct average at the edges
0629         XWidth = neighborsX + groupx + 1;
0630       } else if (neighborsX > (XBlocks - (groupx + 1))) {
0631         XWidth = (XBlocks - groupx) + neighborsX;
0632       } else {
0633         XWidth = 2 * neighborsX + 1;
0634       }
0635       if (neighborsY > groupy) {
0636         YWidth = neighborsY + groupy + 1;
0637       } else if (neighborsY > (YBlocks - (groupy + 1))) {
0638         YWidth = (YBlocks - groupy) + neighborsY;
0639       } else {
0640         YWidth = 2 * neighborsY + 1;
0641       }
0642 
0643       double average = sum / (XWidth * YWidth - 1);
0644       double sigma = getNeighborSigma(average, groupx, groupy, XBlocks, YBlocks, neighborsX, neighborsY, h);
0645       //get rid of block being tested just like we did with the average
0646       sigma -= (average - blocksum) * (average - blocksum);
0647       double sigma_2 = sqrt(sigma) / sqrt(XWidth * YWidth - 2);  //N-1 where N=XWidth*YWidth - 1
0648       double sigma_real = sigma_2 / (XWidth * YWidth - 1);
0649       //double avg_uncrt = average*sqrt(sum)/sum;//Obsolete now(Chad Freer)
0650       double avg_uncrt = 3 * sigma_real;
0651 
0652       double probNoisy = ROOT::Math::poisson_cdf_c(blocksum - 1, average + avg_uncrt);
0653       double probDead = ROOT::Math::poisson_cdf(blocksum, average - avg_uncrt);
0654       double thresholdNoisy = ROOT::Math::normal_cdf_c(toleranceNoisy_);
0655       double thresholdDead = ROOT::Math::normal_cdf(-1 * toleranceDead_);
0656 
0657       int failureNoisy = 0;
0658       int failureDead = 0;
0659 
0660       if (average != 0) {
0661         if (noisy_ && dead_) {
0662           if (blocksum > average) {
0663             failureNoisy = probNoisy < thresholdNoisy;
0664           } else {
0665             failureDead = probDead < thresholdDead;
0666           }
0667         } else if (noisy_) {
0668           if (blocksum > average) {
0669             failureNoisy = probNoisy < thresholdNoisy;
0670           }
0671         } else if (dead_) {
0672           if (blocksum < average) {
0673             failureDead = probDead < thresholdDead;
0674           }
0675         } else {
0676           std::cout << "No test type selected!\n";
0677         }
0678         //Following lines useful for debugging using verbose (Chad Freer)
0679         //string histName = h->GetName();
0680         //if (histName == "emtfTrackBX") {
0681         //   std::printf("Chad says: %i XBlocks, %i XBlocks, %f Blocksum, %f Average", XBlocks,YBlocks,blocksum,average);}
0682       }
0683 
0684       if (failureNoisy || failureDead) {
0685         ++fail;
0686         //DQMChannel chan(groupx*Xbinnum+xMin+binx, 0, 0, blocksum, h->GetBinError(groupx*Xbinnum+xMin+binx));
0687         //badChannels_.push_back(chan);
0688       }
0689     }
0690   }
0691   return 1. * ((XBlocks * YBlocks) - fail) / (XBlocks * YBlocks);
0692 }
0693 
0694 //Gets the sum of the bins surrounding the block to be tested (Chad Freer)
0695 double ContentSigma::getNeighborSum(unsigned groupx,
0696                                     unsigned groupy,
0697                                     unsigned Xblocks,
0698                                     unsigned Yblocks,
0699                                     unsigned neighborsX,
0700                                     unsigned neighborsY,
0701                                     const TH1* h) const {
0702   double sum = 0;
0703   unsigned nbinsX = h->GetXaxis()->GetNbins();
0704   unsigned nbinsY = h->GetYaxis()->GetNbins();
0705   unsigned xMin = 1;
0706   unsigned yMin = 1;
0707   unsigned xMax = nbinsX;
0708   unsigned yMax = nbinsY;
0709   unsigned Xbinnum = 1;
0710   unsigned Ybinnum = 1;
0711 
0712   //give users option for automatic mininum and maximum selection by inputting 0 to any of the parameters
0713   // check that user's parameters are completely in agreement with histogram
0714   // for instance, if inputted xMax is out of range xMin will automatically be ignored
0715   if (xMin_ != 0 && xMax_ != 0) {
0716     if ((xMax_ <= nbinsX) && (xMin_ <= xMax_)) {
0717       nbinsX = xMax_ - xMin_ + 1;
0718       xMax = xMax_;  // do NOT use overflow bin
0719       xMin = xMin_;  // do NOT use underflow bin
0720     }
0721   }
0722   if (yMin_ != 0 && yMax_ != 0) {
0723     if ((yMax_ <= nbinsY) && (yMin_ <= yMax_)) {
0724       nbinsY = yMax_ - yMin_ + 1;
0725       yMax = yMax_;
0726       yMin = yMin_;
0727     }
0728   }
0729 
0730   if (Xblocks == 999) {  //Check to see if blocks should be ignored
0731     Xblocks = nbinsX;
0732   }
0733   if (Yblocks == 999) {
0734     Yblocks = nbinsY;
0735   }
0736 
0737   Xbinnum = nbinsX / Xblocks;
0738   Ybinnum = nbinsY / Yblocks;
0739 
0740   unsigned xLow, xHi, yLow, yHi;  //Define the neighbor blocks edges to be summed
0741   if (groupx > neighborsX) {
0742     xLow = (groupx - neighborsX) * Xbinnum + xMin;
0743     if (xLow < xMin) {
0744       xLow = xMin;  //If the neigbor block would go outside the histogram edge, set it the edge
0745     }
0746   } else {
0747     xLow = xMin;
0748   }
0749   xHi = (groupx + 1 + neighborsX) * Xbinnum + xMin;
0750   if (xHi > xMax) {
0751     xHi = xMax;
0752   }
0753   if (groupy > neighborsY) {
0754     yLow = (groupy - neighborsY) * Ybinnum + yMin;
0755     if (yLow < yMin) {
0756       yLow = yMin;
0757     }
0758   } else {
0759     yLow = yMin;
0760   }
0761   yHi = (groupy + 1 + neighborsY) * Ybinnum + yMin;
0762   if (yHi > yMax) {
0763     yHi = yMax;
0764   }
0765 
0766   for (unsigned xbin = xLow; xbin <= xHi; ++xbin) {  //now sum over all the bins
0767     for (unsigned ybin = yLow; ybin <= yHi; ++ybin) {
0768       sum += h->GetBinContent(xbin, ybin);
0769     }
0770   }
0771   return sum;
0772 }
0773 
0774 //Similar to algorithm  above but returns a version of standard deviation. Additional operations to return real standard deviation used above (Chad Freer)
0775 double ContentSigma::getNeighborSigma(double average,
0776                                       unsigned groupx,
0777                                       unsigned groupy,
0778                                       unsigned Xblocks,
0779                                       unsigned Yblocks,
0780                                       unsigned neighborsX,
0781                                       unsigned neighborsY,
0782                                       const TH1* h) const {
0783   double sigma = 0;
0784   unsigned nbinsX = h->GetXaxis()->GetNbins();
0785   unsigned nbinsY = h->GetYaxis()->GetNbins();
0786   unsigned xMin = 1;
0787   unsigned yMin = 1;
0788   unsigned xMax = nbinsX;
0789   unsigned yMax = nbinsY;
0790   unsigned Xbinnum = 1;
0791   unsigned Ybinnum = 1;
0792   double block_sum;
0793 
0794   if (xMin_ != 0 && xMax_ != 0) {
0795     if ((xMax_ <= nbinsX) && (xMin_ <= xMax_)) {
0796       nbinsX = xMax_ - xMin_ + 1;
0797       xMax = xMax_;
0798       xMin = xMin_;
0799     }
0800   }
0801   if (yMin_ != 0 && yMax_ != 0) {
0802     if ((yMax_ <= nbinsY) && (yMin_ <= yMax_)) {
0803       nbinsY = yMax_ - yMin_ + 1;
0804       yMax = yMax_;
0805       yMin = yMin_;
0806     }
0807   }
0808   if (Xblocks == 999) {
0809     Xblocks = nbinsX;
0810   }
0811   if (Yblocks == 999) {
0812     Yblocks = nbinsY;
0813   }
0814 
0815   Xbinnum = nbinsX / Xblocks;
0816   Ybinnum = nbinsY / Yblocks;
0817 
0818   unsigned xLow, xHi, yLow, yHi;
0819   for (unsigned x_block_count = 0; x_block_count <= 2 * neighborsX; ++x_block_count) {
0820     for (unsigned y_block_count = 0; y_block_count <= 2 * neighborsY; ++y_block_count) {
0821       //Sum over blocks. Need to find standard deviation of average of blocksums. Set up low and hi values similar to sum but for blocks now.
0822       if (groupx + x_block_count > neighborsX) {
0823         xLow = (groupx + x_block_count - neighborsX) * Xbinnum + xMin;
0824         if (xLow < xMin) {
0825           xLow = xMin;
0826         }
0827       } else {
0828         xLow = xMin;
0829       }
0830       xHi = xLow + Xbinnum;
0831       if (xHi > xMax) {
0832         xHi = xMax;
0833       }
0834       if (groupy + y_block_count > neighborsY) {
0835         yLow = (groupy + y_block_count - neighborsY) * Ybinnum + yMin;
0836         if (yLow < yMin) {
0837           yLow = yMin;
0838         }
0839       } else {
0840         yLow = yMin;
0841       }
0842       yHi = yLow + Ybinnum;
0843       if (yHi > yMax) {
0844         yHi = yMax;
0845       }
0846       block_sum = 0;
0847       for (unsigned x_block_bin = xLow; x_block_bin <= xHi; ++x_block_bin) {
0848         for (unsigned y_block_bin = yLow; y_block_bin <= yHi; ++y_block_bin) {
0849           block_sum += h->GetBinContent(x_block_bin, y_block_bin);
0850         }
0851       }
0852       sigma += (average - block_sum) * (average - block_sum);  //will sqrt and divide by sqrt(N-1) outside of function
0853     }
0854   }
0855   return sigma;
0856 }
0857 
0858 //-----------------------------------------------------------//
0859 //----------------  ContentsWithinExpected ---------------------//
0860 //-----------------------------------------------------------//
0861 // run the test (result: fraction of channels that passed test);
0862 // [0, 1] or <0 for failure
0863 float ContentsWithinExpected::runTest(const MonitorElement* me) {
0864   badChannels_.clear();
0865   if (!me)
0866     return -1;
0867   if (!me->getRootObject())
0868     return -1;
0869   TH1* h = nullptr;  //initialize histogram pointer
0870 
0871   if (verbose_ > 1)
0872     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
0873 
0874   int ncx;
0875   int ncy;
0876 
0877   if (useEmptyBins_) {
0878     //-- TH2
0879     if (me->kind() == MonitorElement::Kind::TH2F) {
0880       ncx = me->getTH2F()->GetXaxis()->GetNbins();
0881       ncy = me->getTH2F()->GetYaxis()->GetNbins();
0882       h = me->getTH2F();  // access Test histo
0883     }
0884     //-- TH2S
0885     else if (me->kind() == MonitorElement::Kind::TH2S) {
0886       ncx = me->getTH2S()->GetXaxis()->GetNbins();
0887       ncy = me->getTH2S()->GetYaxis()->GetNbins();
0888       h = me->getTH2S();  // access Test histo
0889     }
0890     //-- TH2D
0891     else if (me->kind() == MonitorElement::Kind::TH2D) {
0892       ncx = me->getTH2D()->GetXaxis()->GetNbins();
0893       ncy = me->getTH2D()->GetYaxis()->GetNbins();
0894       h = me->getTH2D();  // access Test histo
0895     }
0896     //-- TProfile
0897     else if (me->kind() == MonitorElement::Kind::TPROFILE) {
0898       ncx = me->getTProfile()->GetXaxis()->GetNbins();
0899       ncy = 1;
0900       h = me->getTProfile();  // access Test histo
0901     }
0902     //-- TProfile2D
0903     else if (me->kind() == MonitorElement::Kind::TPROFILE2D) {
0904       ncx = me->getTProfile2D()->GetXaxis()->GetNbins();
0905       ncy = me->getTProfile2D()->GetYaxis()->GetNbins();
0906       h = me->getTProfile2D();  // access Test histo
0907     } else {
0908       if (verbose_ > 0)
0909         std::cout << "QTest:ContentsWithinExpected"
0910                   << " ME does not contain TH2F/TH2S/TH2D/TPROFILE/TPROFILE2D, exiting\n";
0911       return -1;
0912     }
0913 
0914     int nsum = 0;
0915     double sum = 0.0;
0916     double average = 0.0;
0917 
0918     if (checkMeanTolerance_) {  // calculate average value of all bin contents
0919 
0920       for (int cx = 1; cx <= ncx; ++cx) {
0921         for (int cy = 1; cy <= ncy; ++cy) {
0922           if (me->kind() == MonitorElement::Kind::TH2F) {
0923             sum += h->GetBinContent(h->GetBin(cx, cy));
0924             ++nsum;
0925           } else if (me->kind() == MonitorElement::Kind::TH2S) {
0926             sum += h->GetBinContent(h->GetBin(cx, cy));
0927             ++nsum;
0928           } else if (me->kind() == MonitorElement::Kind::TH2D) {
0929             sum += h->GetBinContent(h->GetBin(cx, cy));
0930             ++nsum;
0931           } else if (me->kind() == MonitorElement::Kind::TPROFILE) {
0932             if (me->getTProfile()->GetBinEntries(h->GetBin(cx)) >= minEntries_ / (ncx)) {
0933               sum += h->GetBinContent(h->GetBin(cx));
0934               ++nsum;
0935             }
0936           } else if (me->kind() == MonitorElement::Kind::TPROFILE2D) {
0937             if (me->getTProfile2D()->GetBinEntries(h->GetBin(cx, cy)) >= minEntries_ / (ncx * ncy)) {
0938               sum += h->GetBinContent(h->GetBin(cx, cy));
0939               ++nsum;
0940             }
0941           }
0942         }
0943       }
0944 
0945       if (nsum > 0)
0946         average = sum / nsum;
0947 
0948     }  // calculate average value of all bin contents
0949 
0950     int fail = 0;
0951 
0952     for (int cx = 1; cx <= ncx; ++cx) {
0953       for (int cy = 1; cy <= ncy; ++cy) {
0954         bool failMean = false;
0955         bool failRMS = false;
0956         bool failMeanTolerance = false;
0957 
0958         if (me->kind() == MonitorElement::Kind::TPROFILE &&
0959             me->getTProfile()->GetBinEntries(h->GetBin(cx)) < minEntries_ / (ncx))
0960           continue;
0961 
0962         if (me->kind() == MonitorElement::Kind::TPROFILE2D &&
0963             me->getTProfile2D()->GetBinEntries(h->GetBin(cx, cy)) < minEntries_ / (ncx * ncy))
0964           continue;
0965 
0966         if (checkMean_) {
0967           double mean = h->GetBinContent(h->GetBin(cx, cy));
0968           failMean = (mean < minMean_ || mean > maxMean_);
0969         }
0970 
0971         if (checkRMS_) {
0972           double rms = h->GetBinError(h->GetBin(cx, cy));
0973           failRMS = (rms < minRMS_ || rms > maxRMS_);
0974         }
0975 
0976         if (checkMeanTolerance_) {
0977           double mean = h->GetBinContent(h->GetBin(cx, cy));
0978           failMeanTolerance = (std::abs(mean - average) > toleranceMean_ * std::abs(average));
0979         }
0980 
0981         if (failMean || failRMS || failMeanTolerance) {
0982           if (me->kind() == MonitorElement::Kind::TH2F) {
0983             DQMChannel chan(cx, cy, 0, h->GetBinContent(h->GetBin(cx, cy)), h->GetBinError(h->GetBin(cx, cy)));
0984             badChannels_.push_back(chan);
0985           } else if (me->kind() == MonitorElement::Kind::TH2S) {
0986             DQMChannel chan(cx, cy, 0, h->GetBinContent(h->GetBin(cx, cy)), h->GetBinError(h->GetBin(cx, cy)));
0987             badChannels_.push_back(chan);
0988           } else if (me->kind() == MonitorElement::Kind::TH2D) {
0989             DQMChannel chan(cx, cy, 0, h->GetBinContent(h->GetBin(cx, cy)), h->GetBinError(h->GetBin(cx, cy)));
0990             badChannels_.push_back(chan);
0991           } else if (me->kind() == MonitorElement::Kind::TPROFILE) {
0992             DQMChannel chan(
0993                 cx, cy, int(me->getTProfile()->GetBinEntries(h->GetBin(cx))), 0, h->GetBinError(h->GetBin(cx)));
0994             badChannels_.push_back(chan);
0995           } else if (me->kind() == MonitorElement::Kind::TPROFILE2D) {
0996             DQMChannel chan(cx,
0997                             cy,
0998                             int(me->getTProfile2D()->GetBinEntries(h->GetBin(cx, cy))),
0999                             h->GetBinContent(h->GetBin(cx, cy)),
1000                             h->GetBinError(h->GetBin(cx, cy)));
1001             badChannels_.push_back(chan);
1002           }
1003           ++fail;
1004         }
1005       }
1006     }
1007     return 1. * (ncx * ncy - fail) / (ncx * ncy);
1008   }  /// end of normal Test
1009 
1010   else  /// AS quality test !!!
1011   {
1012     if (me->kind() == MonitorElement::Kind::TH2F) {
1013       ncx = me->getTH2F()->GetXaxis()->GetNbins();
1014       ncy = me->getTH2F()->GetYaxis()->GetNbins();
1015       h = me->getTH2F();  // access Test histo
1016     } else if (me->kind() == MonitorElement::Kind::TH2S) {
1017       ncx = me->getTH2S()->GetXaxis()->GetNbins();
1018       ncy = me->getTH2S()->GetYaxis()->GetNbins();
1019       h = me->getTH2S();  // access Test histo
1020     } else if (me->kind() == MonitorElement::Kind::TH2D) {
1021       ncx = me->getTH2D()->GetXaxis()->GetNbins();
1022       ncy = me->getTH2D()->GetYaxis()->GetNbins();
1023       h = me->getTH2D();  // access Test histo
1024     } else {
1025       if (verbose_ > 0)
1026         std::cout << "QTest:ContentsWithinExpected AS"
1027                   << " ME does not contain TH2F/TH2S/TH2D, exiting\n";
1028       return -1;
1029     }
1030 
1031     // if (!rangeInitialized_) return 0; // all accepted if no initialization
1032     int fail = 0;
1033     for (int cx = 1; cx <= ncx; ++cx) {
1034       for (int cy = 1; cy <= ncy; ++cy) {
1035         bool failure = false;
1036         double Content = h->GetBinContent(h->GetBin(cx, cy));
1037         if (Content)
1038           failure = (Content < minMean_ || Content > maxMean_);
1039         if (failure)
1040           ++fail;
1041       }
1042     }
1043     return 1. * (ncx * ncy - fail) / (ncx * ncy);
1044   }  /// end of AS quality test
1045 }
1046 /// set expected value for mean
1047 void ContentsWithinExpected::setMeanRange(double xmin, double xmax) {
1048   if (xmax < xmin)
1049     if (verbose_ > 0)
1050       std::cout << "QTest:ContentsWitinExpected"
1051                 << " Illogical range: (" << xmin << ", " << xmax << "\n";
1052   minMean_ = xmin;
1053   maxMean_ = xmax;
1054   checkMean_ = true;
1055 }
1056 
1057 /// set expected value for mean
1058 void ContentsWithinExpected::setRMSRange(double xmin, double xmax) {
1059   if (xmax < xmin)
1060     if (verbose_ > 0)
1061       std::cout << "QTest:ContentsWitinExpected"
1062                 << " Illogical range: (" << xmin << ", " << xmax << "\n";
1063   minRMS_ = xmin;
1064   maxRMS_ = xmax;
1065   checkRMS_ = true;
1066 }
1067 
1068 //----------------------------------------------------------------//
1069 //--------------------  MeanWithinExpected  ---------------------//
1070 //---------------------------------------------------------------//
1071 // run the test;
1072 //   (a) if useRange is called: 1 if mean within allowed range, 0 otherwise
1073 //   (b) is useRMS or useSigma is called: result is the probability
1074 //   Prob(chi^2, ndof=1) that the mean of histogram will be deviated by more than
1075 //   +/- delta from <expected_mean>, where delta = mean - <expected_mean>, and
1076 //   chi^2 = (delta/sigma)^2. sigma is the RMS of the histogram ("useRMS") or
1077 //   <expected_sigma> ("useSigma")
1078 //   e.g. for delta = 1, Prob = 31.7%
1079 //  for delta = 2, Prob = 4.55%
1080 //   (returns result in [0, 1] or <0 for failure)
1081 float MeanWithinExpected::runTest(const MonitorElement* me) {
1082   if (!me)
1083     return -1;
1084   if (!me->getRootObject())
1085     return -1;
1086   TH1* h = nullptr;
1087 
1088   if (verbose_ > 1)
1089     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
1090 
1091   if (minEntries_ != 0 && me->getEntries() < minEntries_)
1092     return -1;
1093 
1094   if (me->kind() == MonitorElement::Kind::TH1F) {
1095     h = me->getTH1F();  //access Test histo
1096   } else if (me->kind() == MonitorElement::Kind::TH1S) {
1097     h = me->getTH1S();  //access Test histo
1098   } else if (me->kind() == MonitorElement::Kind::TH1D) {
1099     h = me->getTH1D();  //access Test histo
1100   } else {
1101     if (verbose_ > 0)
1102       std::cout << "QTest:MeanWithinExpected"
1103                 << " ME " << me->getFullname() << " does not contain TH1F/TH1S/TH1D, exiting\n";
1104     return -1;
1105   }
1106 
1107   if (useRange_) {
1108     double mean = h->GetMean();
1109     if (mean <= xmax_ && mean >= xmin_)
1110       return 1;
1111     else
1112       return 0;
1113   } else if (useSigma_) {
1114     if (sigma_ != 0.) {
1115       double chi = (h->GetMean() - expMean_) / sigma_;
1116       return TMath::Prob(chi * chi, 1);
1117     } else {
1118       if (verbose_ > 0)
1119         std::cout << "QTest:MeanWithinExpected"
1120                   << " Error, sigma_ is zero, exiting\n";
1121       return 0;
1122     }
1123   } else if (useRMS_) {
1124     if (h->GetRMS() != 0.) {
1125       double chi = (h->GetMean() - expMean_) / h->GetRMS();
1126       return TMath::Prob(chi * chi, 1);
1127     } else {
1128       if (verbose_ > 0)
1129         std::cout << "QTest:MeanWithinExpected"
1130                   << " Error, RMS is zero, exiting\n";
1131       return 0;
1132     }
1133   } else {
1134     if (verbose_ > 0)
1135       std::cout << "QTest:MeanWithinExpected"
1136                 << " Error, neither Range, nor Sigma, nor RMS, exiting\n";
1137     return -1;
1138   }
1139 }
1140 
1141 void MeanWithinExpected::useRange(double xmin, double xmax) {
1142   useRange_ = true;
1143   useSigma_ = useRMS_ = false;
1144   xmin_ = xmin;
1145   xmax_ = xmax;
1146   if (xmin_ > xmax_)
1147     if (verbose_ > 0)
1148       std::cout << "QTest:MeanWithinExpected"
1149                 << " Illogical range: (" << xmin_ << ", " << xmax_ << "\n";
1150 }
1151 void MeanWithinExpected::useSigma(double expectedSigma) {
1152   useSigma_ = true;
1153   useRMS_ = useRange_ = false;
1154   sigma_ = expectedSigma;
1155   if (sigma_ == 0)
1156     if (verbose_ > 0)
1157       std::cout << "QTest:MeanWithinExpected"
1158                 << " Expected sigma = " << sigma_ << "\n";
1159 }
1160 
1161 void MeanWithinExpected::useRMS() {
1162   useRMS_ = true;
1163   useSigma_ = useRange_ = false;
1164 }
1165 
1166 //----------------------------------------------------------------//
1167 //------------------------  CompareToMedian  ---------------------------//
1168 //----------------------------------------------------------------//
1169 /* 
1170 Test for TProfile2D
1171 For each x bin, the median value is calculated and then each value is compared with the median.
1172 This procedure is repeated for each x-bin of the 2D profile
1173 The parameters used for this comparison are:
1174 MinRel and MaxRel to identify outliers wrt the median value
1175 An absolute value (MinAbs, MaxAbs) on the median is used to identify a full region out of specification 
1176 */
1177 float CompareToMedian::runTest(const MonitorElement* me) {
1178   int32_t nbins = 0, failed = 0;
1179   badChannels_.clear();
1180 
1181   if (!me)
1182     return -1;
1183   if (!me->getRootObject())
1184     return -1;
1185   TH1* h = nullptr;
1186 
1187   if (verbose_ > 1) {
1188     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
1189     std::cout << "\tMin = " << _min << "; Max = " << _max << "\n";
1190     std::cout << "\tMinMedian = " << _minMed << "; MaxMedian = " << _maxMed << "\n";
1191     std::cout << "\tUseEmptyBins = " << (_emptyBins ? "Yes" : "No") << "\n";
1192   }
1193 
1194   if (me->kind() == MonitorElement::Kind::TPROFILE2D) {
1195     h = me->getTProfile2D();  // access Test histo
1196   } else {
1197     if (verbose_ > 0)
1198       std::cout << "QTest:ContentsWithinExpected"
1199                 << " ME does not contain TPROFILE2D, exiting\n";
1200     return -1;
1201   }
1202 
1203   nBinsX = h->GetNbinsX();
1204   nBinsY = h->GetNbinsY();
1205   int entries = 0;
1206   float median = 0.0;
1207 
1208   //Median calculated with partially sorted vector
1209   for (int binX = 1; binX <= nBinsX; binX++) {
1210     reset();
1211     // Fill vector
1212     for (int binY = 1; binY <= nBinsY; binY++) {
1213       int bin = h->GetBin(binX, binY);
1214       auto content = (double)h->GetBinContent(bin);
1215       if (content == 0 && !_emptyBins)
1216         continue;
1217       binValues.push_back(content);
1218       entries = me->getTProfile2D()->GetBinEntries(bin);
1219     }
1220     if (binValues.empty())
1221       continue;
1222     nbins += binValues.size();
1223 
1224     //calculate median
1225     if (!binValues.empty()) {
1226       int medPos = (int)binValues.size() / 2;
1227       nth_element(binValues.begin(), binValues.begin() + medPos, binValues.end());
1228       median = *(binValues.begin() + medPos);
1229     }
1230 
1231     // if median == 0, use the absolute cut
1232     if (median == 0) {
1233       if (verbose_ > 0) {
1234         std::cout << "QTest: Median is 0; the fixed cuts: [" << _minMed << "; " << _maxMed << "]  are used\n";
1235       }
1236       for (int binY = 1; binY <= nBinsY; binY++) {
1237         int bin = h->GetBin(binX, binY);
1238         auto content = (double)h->GetBinContent(bin);
1239         entries = me->getTProfile2D()->GetBinEntries(bin);
1240         if (entries == 0)
1241           continue;
1242         if (content > _maxMed || content < _minMed) {
1243           DQMChannel chan(binX, binY, 0, content, h->GetBinError(bin));
1244           badChannels_.push_back(chan);
1245           failed++;
1246         }
1247       }
1248       continue;
1249     }
1250 
1251     //Cut on stat: will mask rings with no enought of statistics
1252     if (median * entries < _statCut)
1253       continue;
1254 
1255     // If median is off the absolute cuts, declare everything bad (if bin has non zero entries)
1256     if (median > _maxMed || median < _minMed) {
1257       for (int binY = 1; binY <= nBinsY; binY++) {
1258         int bin = h->GetBin(binX, binY);
1259         auto content = (double)h->GetBinContent(bin);
1260         entries = me->getTProfile2D()->GetBinEntries(bin);
1261         if (entries == 0)
1262           continue;
1263         DQMChannel chan(binX, binY, 0, content / median, h->GetBinError(bin));
1264         badChannels_.push_back(chan);
1265         failed++;
1266       }
1267       continue;
1268     }
1269 
1270     // Test itself
1271     float minCut = median * _min;
1272     float maxCut = median * _max;
1273     for (int binY = 1; binY <= nBinsY; binY++) {
1274       int bin = h->GetBin(binX, binY);
1275       auto content = (double)h->GetBinContent(bin);
1276       entries = me->getTProfile2D()->GetBinEntries(bin);
1277       if (entries == 0)
1278         continue;
1279       if (content > maxCut || content < minCut) {
1280         DQMChannel chan(binX, binY, 0, content / median, h->GetBinError(bin));
1281         badChannels_.push_back(chan);
1282         failed++;
1283       }
1284     }
1285   }
1286 
1287   if (nbins == 0) {
1288     if (verbose_ > 0)
1289       std::cout << "QTest:CompareToMedian: Histogram is empty" << std::endl;
1290     return 1.;
1291   }
1292   return 1 - (float)failed / nbins;
1293 }
1294 //----------------------------------------------------------------//
1295 //------------------------  CompareLastFilledBin -----------------//
1296 //----------------------------------------------------------------//
1297 /* 
1298 Test for TH1F and TH2F
1299 For the last filled bin the value is compared with specified upper and lower limits. If 
1300 it is outside the limit the test failed test result is returned
1301 The parameters used for this comparison are:
1302 MinRel and MaxRel to check identify outliers wrt the median value
1303 */
1304 float CompareLastFilledBin::runTest(const MonitorElement* me) {
1305   if (!me)
1306     return -1;
1307   if (!me->getRootObject())
1308     return -1;
1309   TH1* h1 = nullptr;
1310   TH2* h2 = nullptr;
1311   if (verbose_ > 1) {
1312     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
1313     std::cout << "\tMin = " << _min << "; Max = " << _max << "\n";
1314   }
1315   if (me->kind() == MonitorElement::Kind::TH1F) {
1316     h1 = me->getTH1F();  // access Test histo
1317   } else if (me->kind() == MonitorElement::Kind::TH2F) {
1318     h2 = me->getTH2F();  // access Test histo
1319   } else {
1320     if (verbose_ > 0)
1321       std::cout << "QTest:ContentsWithinExpected"
1322                 << " ME does not contain TH1F or TH2F, exiting\n";
1323     return -1;
1324   }
1325   int lastBinX = 0;
1326   int lastBinY = 0;
1327   float lastBinVal;
1328 
1329   //--------- do the quality test for 1D histo ---------------//
1330   if (h1 != nullptr) {
1331     lastBinX = h1->FindLastBinAbove(_average, 1);
1332     lastBinVal = h1->GetBinContent(lastBinX);
1333     if (h1->GetEntries() == 0 || lastBinVal < 0)
1334       return 1;
1335   } else if (h2 != nullptr) {
1336     lastBinX = h2->FindLastBinAbove(_average, 1);
1337     lastBinY = h2->FindLastBinAbove(_average, 2);
1338     if (h2->GetEntries() == 0 || lastBinX < 0 || lastBinY < 0)
1339       return 1;
1340     lastBinVal = h2->GetBinContent(h2->GetBin(lastBinX, lastBinY));
1341   } else {
1342     if (verbose_ > 0)
1343       std::cout << "QTest:" << getAlgoName() << " Histogram does not exist" << std::endl;
1344     return 1;
1345   }
1346   if (verbose_ > 0)
1347     std::cout << "Min and Max values " << _min << " " << _max << " Av value " << _average << " lastBinX " << lastBinX
1348               << " lastBinY " << lastBinY << " lastBinVal " << lastBinVal << std::endl;
1349   if (lastBinVal > _min && lastBinVal <= _max)
1350     return 1;
1351   else
1352     return 0;
1353 }
1354 //----------------------------------------------------//
1355 //--------------- CheckVariance ---------------------//
1356 //----------------------------------------------------//
1357 float CheckVariance::runTest(const MonitorElement* me) {
1358   badChannels_.clear();
1359 
1360   if (!me)
1361     return -1;
1362   if (!me->getRootObject())
1363     return -1;
1364   TH1* h = nullptr;
1365 
1366   if (verbose_ > 1)
1367     std::cout << "QTest:" << getAlgoName() << "::runTest called on " << me->getFullname() << "\n";
1368   // -- TH1F
1369   if (me->kind() == MonitorElement::Kind::TH1F) {
1370     h = me->getTH1F();
1371   }
1372   // -- TH1D
1373   else if (me->kind() == MonitorElement::Kind::TH1D) {
1374     h = me->getTH1D();
1375   } else if (me->kind() == MonitorElement::Kind::TPROFILE) {
1376     h = me->getTProfile();  // access Test histo
1377   } else {
1378     if (verbose_ > 0)
1379       std::cout << "QTest:CheckVariance"
1380                 << " ME " << me->getFullname() << " does not contain TH1F/TH1D/TPROFILE, exiting\n";
1381     return -1;
1382   }
1383 
1384   int ncx = h->GetXaxis()->GetNbins();
1385 
1386   double sum = 0;
1387   double sum2 = 0;
1388   for (int bin = 1; bin <= ncx; ++bin) {
1389     double contents = h->GetBinContent(bin);
1390     sum += contents;
1391   }
1392   if (sum == 0)
1393     return -1;
1394   double avg = sum / ncx;
1395 
1396   for (int bin = 1; bin <= ncx; ++bin) {
1397     double contents = h->GetBinContent(bin);
1398     sum2 += (contents - avg) * (contents - avg);
1399   }
1400 
1401   double Variance = TMath::Sqrt(sum2 / ncx);
1402   return Variance;
1403 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!