EcalCommon/src/MESetTrend.cc

0001 #include "DQM/EcalCommon/interface/MESetTrend.h"
0002
0003 #include <ctime>
0004
0005 namespace ecaldqm {
0006
0007   MESetTrend::MESetTrend(std::string const &_fullPath,
0008                          binning::ObjectType _otype,
0009                          binning::BinningType _btype,
0010                          MonitorElement::Kind _kind,
0011                          binning::AxisSpecs const *_xaxis /* = 0*/,
0012                          binning::AxisSpecs const *_yaxis /* = 0*/)
0013       : MESetEcal(_fullPath, _otype, _btype, _kind, 1, _xaxis, _yaxis),
0014         minutely_(false),
0015         shiftAxis_(false),
0016         currentBin_(-1) {
0017     switch (kind_) {
0018       case MonitorElement::Kind::TH1F:
0019       case MonitorElement::Kind::TH2F:
0020       case MonitorElement::Kind::TPROFILE:
0021       case MonitorElement::Kind::TPROFILE2D:
0022         break;
0023       default:
0024         throw_("Unsupported MonitorElement kind");
0025     }
0026   }
0027
0028   MESetTrend::MESetTrend(MESetTrend const &_orig)
0029       : MESetEcal(_orig), minutely_(_orig.minutely_), shiftAxis_(_orig.shiftAxis_), currentBin_(_orig.currentBin_) {}
0030
0031   MESet &MESetTrend::operator=(MESet const &_rhs) {
0032     MESetTrend const *pRhs(dynamic_cast<MESetTrend const *>(&_rhs));
0033     if (pRhs) {
0034       minutely_ = pRhs->minutely_;
0035       shiftAxis_ = pRhs->shiftAxis_;
0036       currentBin_ = pRhs->currentBin_;
0037     }
0038
0039     return MESetEcal::operator=(_rhs);
0040   }
0041
0042   MESet *MESetTrend::clone(std::string const &_path /* = ""*/) const {
0043     std::string path(path_);
0044     if (!_path.empty())
0045       path_ = _path;
0046     MESet *copy(new MESetTrend(*this));
0047     path_ = path;
0048     return copy;
0049   }
0050
0051   void MESetTrend::book(DQMStore::IBooker &_ibooker, EcalElectronicsMapping const *electronicsMap) {
0052     binning::AxisSpecs xaxis;
0053     if (xaxis_)
0054       xaxis = *xaxis_;
0055     else {
0056       xaxis.nbins = 200;
0057       xaxis.low = 0.;
0058       xaxis.high = 2000.;
0059     }
0060
0061     if (minutely_) {
0062       time_t localTime(time(nullptr));
0063       struct tm timeBuffer;
0064       gmtime_r(&localTime, &timeBuffer);  // gmtime() is not thread safe
0065       unsigned utcTime(mktime(&timeBuffer));
0066
0067       xaxis.low = utcTime;
0068       if (xaxis_)
0069         xaxis.high = utcTime + xaxis_->high - xaxis_->low;
0070       else
0071         xaxis.high = xaxis.low + 200 * 60.;
0072     }
0073
0074     binning::AxisSpecs const *xaxisTemp(xaxis_);
0075     xaxis_ = &xaxis;
0076
0077     MESetEcal::book(_ibooker, electronicsMap);
0078
0079     xaxis_ = xaxisTemp;
0080
0081     if (minutely_) {
0082       for (unsigned iME(0); iME < mes_.size(); ++iME)
0083         mes_[iME]->getTH1()->GetXaxis()->SetTimeDisplay(1);
0084       setAxisTitle("UTC");
0085     } else
0086       setAxisTitle("LumiSections");
0087   }
0088
0089   void MESetTrend::fill(
0090       EcalDQMSetupObjects const edso, DetId const &_id, double _t, double _wy /* = 1.*/, double _w /* = 1.*/) {
0091     if (!active_)
0092       return;
0093
0094     unsigned iME(binning::findPlotIndex(edso.electronicsMap, otype_, _id));
0095     checkME_(iME);
0096
0097     if (shift_(unsigned(_t)))
0098       fill_(iME, _t + 0.5, _wy, _w);
0099   }
0100
0101   void MESetTrend::fill(EcalDQMSetupObjects const edso,
0102                         EcalElectronicsId const &_id,
0103                         double _t,
0104                         double _wy /* = 1.*/,
0105                         double _w /* = 1.*/) {
0106     if (!active_)
0107       return;
0108
0109     unsigned iME(binning::findPlotIndex(edso.electronicsMap, otype_, _id));
0110     checkME_(iME);
0111
0112     if (shift_(unsigned(_t)))
0113       fill_(iME, _t + 0.5, _wy, _w);
0114   }
0115
0116   void MESetTrend::fill(
0117       EcalDQMSetupObjects const edso, int _dcctccid, double _t, double _wy /* = 1.*/, double _w /* = 1.*/) {
0118     if (!active_)
0119       return;
0120
0121     unsigned iME(binning::findPlotIndex(edso.electronicsMap, otype_, _dcctccid, btype_));
0122     checkME_(iME);
0123
0124     if (shift_(unsigned(_t)))
0125       fill_(iME, _t + 0.5, _wy, _w);
0126   }
0127
0128   void MESetTrend::fill(EcalDQMSetupObjects const edso, double _t, double _wy /* = 1.*/, double _w /* = 1.*/) {
0129     if (!active_)
0130       return;
0131     if (mes_.size() != 1)
0132       throw_("MESet type incompatible");
0133
0134     if (shift_(unsigned(_t)))
0135       fill_(0, _t + 0.5, _wy, _w);
0136   }
0137
0138   int MESetTrend::findBin(EcalDQMSetupObjects const edso, DetId const &_id, double _t, double _y /* = 0.*/) const {
0139     if (!active_)
0140       return -1;
0141
0142     unsigned iME(binning::findPlotIndex(edso.electronicsMap, otype_, _id));
0143     checkME_(iME);
0144
0145     return mes_[iME]->getTH1()->FindBin(_t + 0.5, _y);
0146   }
0147
0148   int MESetTrend::findBin(EcalDQMSetupObjects const edso,
0149                           EcalElectronicsId const &_id,
0150                           double _t,
0151                           double _y /* = 0.*/) const {
0152     if (!active_)
0153       return -1;
0154
0155     unsigned iME(binning::findPlotIndex(edso.electronicsMap, otype_, _id));
0156     checkME_(iME);
0157
0158     return mes_[iME]->getTH1()->FindBin(_t + 0.5, _y);
0159   }
0160
0161   int MESetTrend::findBin(EcalDQMSetupObjects const edso, int _dcctccid, double _t, double _y /* = 0.*/) const {
0162     if (!active_)
0163       return -1;
0164
0165     unsigned iME(binning::findPlotIndex(edso.electronicsMap, otype_, _dcctccid, btype_));
0166     checkME_(iME);
0167
0168     return mes_[iME]->getTH1()->FindBin(_t + 0.5, _y);
0169   }
0170
0171   int MESetTrend::findBin(EcalDQMSetupObjects const edso, double _t, double _y /* = 0.*/) const {
0172     if (!active_)
0173       return -1;
0174     if (mes_.size() != 1)
0175       throw_("MESet type incompatible");
0176
0177     return mes_[0]->getTH1()->FindBin(_t + 0.5, _y);
0178   }
0179
0180   void MESetTrend::setCumulative() {
0181     if (kind_ == MonitorElement::Kind::TPROFILE || kind_ == MonitorElement::Kind::TPROFILE2D)
0182       throw_("Cumulative flag set for a profile plot");
0183     currentBin_ = 1;
0184   }
0185
0186   bool MESetTrend::shift_(unsigned _t) {
0187     TAxis *tAxis(mes_[0]->getTH1()->GetXaxis());
0188     int nbinsX(tAxis->GetNbins());
0189     unsigned tLow(tAxis->GetBinLowEdge(1));
0190     unsigned tHigh(tAxis->GetBinUpEdge(nbinsX));
0191
0192     if (!shiftAxis_)
0193       return _t >= tLow && _t < tHigh;
0194
0195     int dBin(0);
0196     int unitsPerBin((tHigh - tLow) / nbinsX);
0197
0198     if (_t >= tLow && _t < tHigh) {
0199       if (currentBin_ > 0) {
0200         int thisBin(tAxis->FindBin(_t + 0.5));
0201         if (thisBin < currentBin_)
0202           return false;
0203         else if (thisBin > currentBin_) {
0204           for (unsigned iME(0); iME < mes_.size(); iME++) {
0205             MonitorElement *me(mes_[iME]);
0206             int nbinsY(me->getTH1()->GetNbinsY());
0207             for (int iy(1); iy <= nbinsY; ++iy) {
0208               int orig(me->getTH1()->GetBin(currentBin_, iy));
0209               double currentContent(me->getBinContent(orig));
0210               double currentError(me->getBinError(orig));
0211               for (int ix(currentBin_); ix <= thisBin; ++ix) {
0212                 int dest(me->getTH1()->GetBin(ix, iy));
0213                 me->setBinContent(dest, currentContent);
0214                 me->setBinError(dest, currentError);
0215               }
0216             }
0217           }
0218         }
0219       }
0220
0221       return true;
0222     } else if (_t >= tHigh) {
0223       dBin = (_t - tHigh) / unitsPerBin + 1;
0224       if (currentBin_ > 0)
0225         currentBin_ = nbinsX;
0226     } else if (_t < tLow) {
0227       if (currentBin_ > 0)
0228         return false;  // no going back in time in case of cumulative history
0229
0230       int maxBin(0);
0231
0232       for (unsigned iME(0); iME < mes_.size(); iME++) {
0233         MonitorElement *me(mes_[iME]);
0234
0235         bool filled(false);
0236         int iMax(nbinsX + 1);
0237         while (--iMax > 0 && !filled) {
0238           switch (kind_) {
0239             case MonitorElement::Kind::TH1F:
0240               if (me->getBinContent(iMax) != 0)
0241                 filled = true;
0242               break;
0243             case MonitorElement::Kind::TPROFILE:
0244               if (me->getBinEntries(iMax) != 0)
0245                 filled = true;
0246               break;
0247             case MonitorElement::Kind::TH2F:
0248               for (int iy(1); iy <= me->getNbinsY(); iy++) {
0249                 if (me->getBinContent(me->getTH1()->GetBin(iMax, iy)) != 0) {
0250                   filled = true;
0251                   break;
0252                 }
0253               }
0254               break;
0255             case MonitorElement::Kind::TPROFILE2D:
0256               for (int iy(1); iy <= me->getNbinsY(); iy++) {
0257                 if (me->getBinEntries(me->getTH1()->GetBin(iMax, iy)) != 0) {
0258                   filled = true;
0259                   break;
0260                 }
0261               }
0262               break;
0263             default:
0264               break;
0265           }
0266         }
0267
0268         if (iMax > maxBin)
0269           maxBin = iMax;
0270       }
0271
0272       if (_t < tLow - (nbinsX - maxBin) * unitsPerBin)
0273         return false;
0274
0275       dBin = (_t - tLow) / unitsPerBin - 1;
0276     }
0277
0278     int start(dBin > 0 ? dBin + 1 : nbinsX + dBin);
0279     int end(dBin > 0 ? nbinsX + 1 : 0);
0280     int step(dBin > 0 ? 1 : -1);
0281
0282     tLow += dBin * unitsPerBin;
0283     tHigh += dBin * unitsPerBin;
0284
0285     for (unsigned iME(0); iME < mes_.size(); iME++) {
0286       MonitorElement *me(mes_[iME]);
0287
0288       me->getTH1()->GetXaxis()->SetLimits(tLow, tHigh);
0289
0290       if ((end - start) / step < 0) {
0291         me->Reset();
0292         continue;
0293       }
0294
0295       me->setEntries(0.);
0296       double entries(0.);
0297
0298       switch (kind_) {
0299         case MonitorElement::Kind::TH1F: {
0300           int ix(start);
0301           for (; ix != end; ix += step) {
0302             double binContent(me->getBinContent(ix));
0303             entries += binContent;
0304             me->setBinContent(ix - dBin, binContent);
0305             me->setBinError(ix - dBin, me->getBinError(ix));
0306           }
0307           ix = end - dBin - 1 * step;
0308           double lastContent(currentBin_ > 0 ? me->getBinContent(ix) : 0.);
0309           double lastError(currentBin_ > 0 ? me->getBinContent(ix) : 0.);
0310           for (ix += step; ix != end; ix += step) {
0311             me->setBinContent(ix, lastContent);
0312             me->setBinError(ix, lastError);
0313           }
0314         } break;
0315         case MonitorElement::Kind::TPROFILE: {
0316           int ix(start);
0317           for (; ix != end; ix += step) {
0318             double binEntries(me->getBinEntries(ix));
0319             double binContent(me->getBinContent(ix));
0320             entries += binEntries;
0321             me->setBinEntries(ix - dBin, binEntries);
0322             me->setBinContent(ix - dBin, binContent * binEntries);
0323             if (binEntries > 0) {
0324               double rms(me->getBinError(ix) * std::sqrt(binEntries));
0325               double sumw2((rms * rms + binContent * binContent) * binEntries);
0326               me->setBinError(ix - dBin, std::sqrt(sumw2));
0327             } else
0328               me->setBinError(ix - dBin, 0.);
0329           }
0330           ix = end - dBin;
0331           for (; ix != end; ix += step) {
0332             me->setBinEntries(ix, 0.);
0333             me->setBinContent(ix, 0.);
0334             me->setBinError(ix, 0.);
0335           }
0336         } break;
0337         case MonitorElement::Kind::TH2F: {
0338           int ix(start);
0339           int nbinsY(me->getNbinsY());
0340           for (; ix != end; ix += step) {
0341             for (int iy(1); iy <= nbinsY; iy++) {
0342               int orig(me->getTH1()->GetBin(ix, iy));
0343               int dest(me->getTH1()->GetBin(ix - dBin, iy));
0344               double binContent(me->getBinContent(orig));
0345               entries += binContent;
0346               me->setBinContent(dest, binContent);
0347               me->setBinError(dest, me->getBinError(orig));
0348
0349               me->setBinContent(orig, 0.);
0350               me->setBinError(orig, 0.);
0351             }
0352           }
0353           ix = end - dBin - 1 * step;
0354           std::vector<double> lastContent;
0355           std::vector<double> lastError;
0356           for (int iy(1); iy <= nbinsY; iy++) {
0357             lastContent.push_back(currentBin_ > 0 ? me->getBinContent(ix, iy) : 0.);
0358             lastError.push_back(currentBin_ > 0 ? me->getBinError(ix, iy) : 0.);
0359           }
0360           for (ix += step; ix != end; ix += step) {
0361             for (int iy(1); iy <= nbinsY; iy++) {
0362               int bin(me->getTH1()->GetBin(ix, iy));
0363               me->setBinContent(bin, lastContent[iy - 1]);
0364               me->setBinError(bin, lastError[iy - 1]);
0365             }
0366           }
0367         } break;
0368         case MonitorElement::Kind::TPROFILE2D: {
0369           int ix(start);
0370           int nbinsY(me->getNbinsY());
0371           for (; ix != end; ix += step) {
0372             for (int iy(1); iy <= nbinsY; iy++) {
0373               int orig(me->getTH1()->GetBin(ix, iy));
0374               int dest(me->getTH1()->GetBin(ix - dBin, iy));
0375               double binEntries(me->getBinEntries(orig));
0376               double binContent(me->getBinContent(orig));
0377               entries += binEntries;
0378               me->setBinEntries(dest, binEntries);
0379               me->setBinContent(dest, binContent * binEntries);
0380               if (binEntries > 0) {
0381                 double rms(me->getBinError(orig) * std::sqrt(binEntries));
0382                 double sumw2((rms * rms + binContent * binContent) * binEntries);
0383                 me->setBinError(dest, std::sqrt(sumw2));
0384               } else
0385                 me->setBinError(dest, 0.);
0386             }
0387           }
0388           ix = end - dBin;
0389           for (; ix != end; ix += step) {
0390             for (int iy(1); iy <= nbinsY; iy++) {
0391               int bin(me->getTH1()->GetBin(ix, iy));
0392               me->setBinEntries(bin, 0.);
0393               me->setBinContent(bin, 0.);
0394               me->setBinError(bin, 0.);
0395             }
0396           }
0397         } break;
0398         default:
0399           break;
0400       }
0401
0402       me->setEntries(entries);
0403     }
0404
0405     return true;
0406   }
0407
0408 }  // namespace ecaldqm