MESetTrend.cc

Go to the documentation of this file.

#include "DQM/EcalCommon/interface/MESetTrend.h"
 
#include <ctime>
 
namespace ecaldqm {
 
  MESetTrend::MESetTrend(std::string const &_fullPath,
                         binning::ObjectType _otype,
                         binning::BinningType _btype,
                         MonitorElement::Kind _kind,
                         binning::AxisSpecs const *_xaxis /* = 0*/,
                         binning::AxisSpecs const *_yaxis /* = 0*/)
      : MESetEcal(_fullPath, _otype, _btype, _kind, 1, _xaxis, _yaxis),
        minutely_(false),
        shiftAxis_(false),
        currentBin_(-1) {
    switch (kind_) {
      case MonitorElement::Kind::TH1F:
      case MonitorElement::Kind::TH2F:
      case MonitorElement::Kind::TPROFILE:
      case MonitorElement::Kind::TPROFILE2D:
        break;
      default:
        throw_("Unsupported MonitorElement kind");
    }
  }
 
  MESetTrend::MESetTrend(MESetTrend const &_orig)
      : MESetEcal(_orig), minutely_(_orig.minutely_), shiftAxis_(_orig.shiftAxis_), currentBin_(_orig.currentBin_) {}
 
  MESet &MESetTrend::operator=(MESet const &_rhs) {
    MESetTrend const *pRhs(dynamic_cast<MESetTrend const *>(&_rhs));
    if (pRhs) {
      minutely_ = pRhs->minutely_;
      shiftAxis_ = pRhs->shiftAxis_;
      currentBin_ = pRhs->currentBin_;
    }
 
    return MESetEcal::operator=(_rhs);
  }
 
  MESet *MESetTrend::clone(std::string const &_path /* = ""*/) const {
    std::string path(path_);
    if (!_path.empty())
      path_ = _path;
    MESet *copy(new MESetTrend(*this));
    path_ = path;
    return copy;
  }
 
  void MESetTrend::book(DQMStore::IBooker &_ibooker) {
    binning::AxisSpecs xaxis;
    if (xaxis_)
      xaxis = *xaxis_;
    else {
      xaxis.nbins = 200;
      xaxis.low = 0.;
      xaxis.high = 2000.;
    }
 
    if (minutely_) {
      time_t localTime(time(nullptr));
      struct tm timeBuffer;
      gmtime_r(&localTime, &timeBuffer);  // gmtime() is not thread safe
      unsigned utcTime(mktime(&timeBuffer));
 
      xaxis.low = utcTime;
      if (xaxis_)
        xaxis.high = utcTime + xaxis_->high - xaxis_->low;
      else
        xaxis.high = xaxis.low + 200 * 60.;
    }
 
    binning::AxisSpecs const *xaxisTemp(xaxis_);
    xaxis_ = &xaxis;
 
    MESetEcal::book(_ibooker);
 
    xaxis_ = xaxisTemp;
 
    if (minutely_) {
      for (unsigned iME(0); iME < mes_.size(); ++iME)
        mes_[iME]->getTH1()->GetXaxis()->SetTimeDisplay(1);
      setAxisTitle("UTC");
    } else
      setAxisTitle("LumiSections");
  }
 
  void MESetTrend::fill(DetId const &_id, double _t, double _wy /* = 1.*/, double _w /* = 1.*/) {
    if (!active_)
      return;
 
    unsigned iME(binning::findPlotIndex(otype_, _id));
    checkME_(iME);
 
    if (shift_(unsigned(_t)))
      fill_(iME, _t + 0.5, _wy, _w);
  }
 
  void MESetTrend::fill(EcalElectronicsId const &_id, double _t, double _wy /* = 1.*/, double _w /* = 1.*/) {
    if (!active_)
      return;
 
    unsigned iME(binning::findPlotIndex(otype_, _id));
    checkME_(iME);
 
    if (shift_(unsigned(_t)))
      fill_(iME, _t + 0.5, _wy, _w);
  }
 
  void MESetTrend::fill(int _dcctccid, double _t, double _wy /* = 1.*/, double _w /* = 1.*/) {
    if (!active_)
      return;
 
    unsigned iME(binning::findPlotIndex(otype_, _dcctccid, btype_));
    checkME_(iME);
 
    if (shift_(unsigned(_t)))
      fill_(iME, _t + 0.5, _wy, _w);
  }
 
  void MESetTrend::fill(double _t, double _wy /* = 1.*/, double _w /* = 1.*/) {
    if (!active_)
      return;
    if (mes_.size() != 1)
      throw_("MESet type incompatible");
 
    if (shift_(unsigned(_t)))
      fill_(0, _t + 0.5, _wy, _w);
  }
 
  int MESetTrend::findBin(DetId const &_id, double _t, double _y /* = 0.*/) const {
    if (!active_)
      return -1;
 
    unsigned iME(binning::findPlotIndex(otype_, _id));
    checkME_(iME);
 
    return mes_[iME]->getTH1()->FindBin(_t + 0.5, _y);
  }
 
  int MESetTrend::findBin(EcalElectronicsId const &_id, double _t, double _y /* = 0.*/) const {
    if (!active_)
      return -1;
 
    unsigned iME(binning::findPlotIndex(otype_, _id));
    checkME_(iME);
 
    return mes_[iME]->getTH1()->FindBin(_t + 0.5, _y);
  }
 
  int MESetTrend::findBin(int _dcctccid, double _t, double _y /* = 0.*/) const {
    if (!active_)
      return -1;
 
    unsigned iME(binning::findPlotIndex(otype_, _dcctccid, btype_));
    checkME_(iME);
 
    return mes_[iME]->getTH1()->FindBin(_t + 0.5, _y);
  }
 
  int MESetTrend::findBin(double _t, double _y /* = 0.*/) const {
    if (!active_)
      return -1;
    if (mes_.size() != 1)
      throw_("MESet type incompatible");
 
    return mes_[0]->getTH1()->FindBin(_t + 0.5, _y);
  }
 
  void MESetTrend::setCumulative() {
    if (kind_ == MonitorElement::Kind::TPROFILE || kind_ == MonitorElement::Kind::TPROFILE2D)
      throw_("Cumulative flag set for a profile plot");
    currentBin_ = 1;
  }
 
  bool MESetTrend::shift_(unsigned _t) {
    TAxis *tAxis(mes_[0]->getTH1()->GetXaxis());
    int nbinsX(tAxis->GetNbins());
    unsigned tLow(tAxis->GetBinLowEdge(1));
    unsigned tHigh(tAxis->GetBinUpEdge(nbinsX));
 
    if (!shiftAxis_)
      return _t >= tLow && _t < tHigh;
 
    int dBin(0);
    int unitsPerBin((tHigh - tLow) / nbinsX);
 
    if (_t >= tLow && _t < tHigh) {
      if (currentBin_ > 0) {
        int thisBin(tAxis->FindBin(_t + 0.5));
        if (thisBin < currentBin_)
          return false;
        else if (thisBin > currentBin_) {
          for (unsigned iME(0); iME < mes_.size(); iME++) {
            MonitorElement *me(mes_[iME]);
            int nbinsY(me->getTH1()->GetNbinsY());
            for (int iy(1); iy <= nbinsY; ++iy) {
              int orig(me->getTH1()->GetBin(currentBin_, iy));
              double currentContent(me->getBinContent(orig));
              double currentError(me->getBinError(orig));
              for (int ix(currentBin_); ix <= thisBin; ++ix) {
                int dest(me->getTH1()->GetBin(ix, iy));
                me->setBinContent(dest, currentContent);
                me->setBinError(dest, currentError);
              }
            }
          }
        }
      }
 
      return true;
    } else if (_t >= tHigh) {
      dBin = (_t - tHigh) / unitsPerBin + 1;
      if (currentBin_ > 0)
        currentBin_ = nbinsX;
    } else if (_t < tLow) {
      if (currentBin_ > 0)
        return false;  // no going back in time in case of cumulative history
 
      int maxBin(0);
 
      for (unsigned iME(0); iME < mes_.size(); iME++) {
        MonitorElement *me(mes_[iME]);
 
        bool filled(false);
        int iMax(nbinsX + 1);
        while (--iMax > 0 && !filled) {
          switch (kind_) {
            case MonitorElement::Kind::TH1F:
              if (me->getBinContent(iMax) != 0)
                filled = true;
              break;
            case MonitorElement::Kind::TPROFILE:
              if (me->getBinEntries(iMax) != 0)
                filled = true;
              break;
            case MonitorElement::Kind::TH2F:
              for (int iy(1); iy <= me->getNbinsY(); iy++) {
                if (me->getBinContent(me->getTH1()->GetBin(iMax, iy)) != 0) {
                  filled = true;
                  break;
                }
              }
              break;
            case MonitorElement::Kind::TPROFILE2D:
              for (int iy(1); iy <= me->getNbinsY(); iy++) {
                if (me->getBinEntries(me->getTH1()->GetBin(iMax, iy)) != 0) {
                  filled = true;
                  break;
                }
              }
              break;
            default:
              break;
          }
        }
 
        if (iMax > maxBin)
          maxBin = iMax;
      }
 
      if (_t < tLow - (nbinsX - maxBin) * unitsPerBin)
        return false;
 
      dBin = (_t - tLow) / unitsPerBin - 1;
    }
 
    int start(dBin > 0 ? dBin + 1 : nbinsX + dBin);
    int end(dBin > 0 ? nbinsX + 1 : 0);
    int step(dBin > 0 ? 1 : -1);
 
    tLow += dBin * unitsPerBin;
    tHigh += dBin * unitsPerBin;
 
    for (unsigned iME(0); iME < mes_.size(); iME++) {
      MonitorElement *me(mes_[iME]);
 
      me->getTH1()->GetXaxis()->SetLimits(tLow, tHigh);
 
      if ((end - start) / step < 0) {
        me->Reset();
        continue;
      }
 
      me->setEntries(0.);
      double entries(0.);
 
      switch (kind_) {
        case MonitorElement::Kind::TH1F: {
          int ix(start);
          for (; ix != end; ix += step) {
            double binContent(me->getBinContent(ix));
            entries += binContent;
            me->setBinContent(ix - dBin, binContent);
            me->setBinError(ix - dBin, me->getBinError(ix));
          }
          ix = end - dBin - 1 * step;
          double lastContent(currentBin_ > 0 ? me->getBinContent(ix) : 0.);
          double lastError(currentBin_ > 0 ? me->getBinContent(ix) : 0.);
          for (ix += step; ix != end; ix += step) {
            me->setBinContent(ix, lastContent);
            me->setBinError(ix, lastError);
          }
        } break;
        case MonitorElement::Kind::TPROFILE: {
          int ix(start);
          for (; ix != end; ix += step) {
            double binEntries(me->getBinEntries(ix));
            double binContent(me->getBinContent(ix));
            entries += binEntries;
            me->setBinEntries(ix - dBin, binEntries);
            me->setBinContent(ix - dBin, binContent * binEntries);
            if (binEntries > 0) {
              double rms(me->getBinError(ix) * std::sqrt(binEntries));
              double sumw2((rms * rms + binContent * binContent) * binEntries);
              me->setBinError(ix - dBin, std::sqrt(sumw2));
            } else
              me->setBinError(ix - dBin, 0.);
          }
          ix = end - dBin;
          for (; ix != end; ix += step) {
            me->setBinEntries(ix, 0.);
            me->setBinContent(ix, 0.);
            me->setBinError(ix, 0.);
          }
        } break;
        case MonitorElement::Kind::TH2F: {
          int ix(start);
          int nbinsY(me->getNbinsY());
          for (; ix != end; ix += step) {
            for (int iy(1); iy <= nbinsY; iy++) {
              int orig(me->getTH1()->GetBin(ix, iy));
              int dest(me->getTH1()->GetBin(ix - dBin, iy));
              double binContent(me->getBinContent(orig));
              entries += binContent;
              me->setBinContent(dest, binContent);
              me->setBinError(dest, me->getBinError(orig));
 
              me->setBinContent(orig, 0.);
              me->setBinError(orig, 0.);
            }
          }
          ix = end - dBin - 1 * step;
          std::vector<double> lastContent;
          std::vector<double> lastError;
          for (int iy(1); iy <= nbinsY; iy++) {
            lastContent.push_back(currentBin_ > 0 ? me->getBinContent(ix, iy) : 0.);
            lastError.push_back(currentBin_ > 0 ? me->getBinError(ix, iy) : 0.);
          }
          for (ix += step; ix != end; ix += step) {
            for (int iy(1); iy <= nbinsY; iy++) {
              int bin(me->getTH1()->GetBin(ix, iy));
              me->setBinContent(bin, lastContent[iy - 1]);
              me->setBinError(bin, lastError[iy - 1]);
            }
          }
        } break;
        case MonitorElement::Kind::TPROFILE2D: {
          int ix(start);
          int nbinsY(me->getNbinsY());
          for (; ix != end; ix += step) {
            for (int iy(1); iy <= nbinsY; iy++) {
              int orig(me->getTH1()->GetBin(ix, iy));
              int dest(me->getTH1()->GetBin(ix - dBin, iy));
              double binEntries(me->getBinEntries(orig));
              double binContent(me->getBinContent(orig));
              entries += binEntries;
              me->setBinEntries(dest, binEntries);
              me->setBinContent(dest, binContent * binEntries);
              if (binEntries > 0) {
                double rms(me->getBinError(orig) * std::sqrt(binEntries));
                double sumw2((rms * rms + binContent * binContent) * binEntries);
                me->setBinError(dest, std::sqrt(sumw2));
              } else
                me->setBinError(dest, 0.);
            }
          }
          ix = end - dBin;
          for (; ix != end; ix += step) {
            for (int iy(1); iy <= nbinsY; iy++) {
              int bin(me->getTH1()->GetBin(ix, iy));
              me->setBinEntries(bin, 0.);
              me->setBinContent(bin, 0.);
              me->setBinError(bin, 0.);
            }
          }
        } break;
        default:
          break;
      }
 
      me->setEntries(entries);
    }
 
    return true;
  }
 
}  // namespace ecaldqm