EcalPreshowerMonitorModule/src/ESTrendTask.cc

0001 #include <iostream>
0002 #include <fstream>
0003 #include <vector>
0004 #include <cmath>
0005
0006 #include "FWCore/Framework/interface/Run.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0013 #include "DQMServices/Core/interface/DQMStore.h"
0014 #include "DQM/EcalPreshowerMonitorModule/interface/ESTrendTask.h"
0015 #include "DataFormats/EcalRawData/interface/ESDCCHeaderBlock.h"
0016 #include "DataFormats/EcalRawData/interface/ESKCHIPBlock.h"
0017
0018 using namespace cms;
0019 using namespace edm;
0020 using namespace std;
0021
0022 ESTrendTask::ESTrendTask(const ParameterSet& ps) {
0023   prefixME_ = ps.getUntrackedParameter<string>("prefixME", "EcalPreshower");
0024   rechittoken_ = consumes<ESRecHitCollection>(ps.getParameter<InputTag>("RecHitLabel"));
0025   dccCollections_ = consumes<ESRawDataCollection>(ps.getParameter<InputTag>("ESDCCCollections"));
0026
0027   for (int i = 0; i < 2; ++i)
0028     for (int j = 0; j < 2; ++j) {
0029       hESRecHitTrend_[i][j] = nullptr;
0030       hESRecHitTrendHr_[i][j] = nullptr;
0031     }
0032
0033   hESSLinkErrTrend_ = nullptr;
0034   hESFiberErrTrend_ = nullptr;
0035   hESSLinkErrTrendHr_ = nullptr;
0036   hESFiberErrTrendHr_ = nullptr;
0037
0038   start_time_ = 0;
0039   current_time_ = 0;
0040   last_time_ = 0;
0041
0042   ievt_ = 0;
0043 }
0044
0045 void ESTrendTask::dqmBeginRun(const Run& r, const EventSetup& c) {
0046   start_time_ = (r.beginTime()).unixTime();
0047
0048   //std::cout << "start time : " << start_time_ << std::endl;
0049 }
0050
0051 void ESTrendTask::bookHistograms(DQMStore::IBooker& iBooker, Run const&, EventSetup const&) {
0052   char histo[200];
0053
0054   iBooker.setCurrentFolder(prefixME_ + "/ESTrendTask");
0055
0056   for (int i = 0; i < 2; ++i)
0057     for (int j = 0; j < 2; ++j) {
0058       int iz = (i == 0) ? 1 : -1;
0059       sprintf(histo, "ES Trending RH Occ per 5 mins Z %d P %d", iz, j + 1);
0060       hESRecHitTrend_[i][j] = iBooker.bookProfile(histo, histo, 36, 0.0, 180.0, 100, 0.0, 1.0e6, "s");
0061       hESRecHitTrend_[i][j]->setAxisTitle("Elapse time (Minutes)", 1);
0062       hESRecHitTrend_[i][j]->setAxisTitle("ES RecHit Occupancy / 5 minutes", 2);
0063
0064       sprintf(histo, "ES Trending RH Occ per hour Z %d P %d", iz, j + 1);
0065       hESRecHitTrendHr_[i][j] = iBooker.bookProfile(histo, histo, 24, 0.0, 24.0, 100, 0.0, 1.0e6, "s");
0066       hESRecHitTrendHr_[i][j]->setAxisTitle("Elapse time (Hours)", 1);
0067       hESRecHitTrendHr_[i][j]->setAxisTitle("ES RecHit Occupancy / hour", 2);
0068     }
0069
0070   sprintf(histo, "ES Trending SLink CRC Error per 5 mins");
0071   hESSLinkErrTrend_ = iBooker.bookProfile(histo, histo, 36, 0.0, 180.0, 100, 0.0, 1.0e6, "s");
0072   hESSLinkErrTrend_->setAxisTitle("Elapse time (Minutes)", 1);
0073   hESSLinkErrTrend_->setAxisTitle("ES SLink CRC Err / 5 minutes", 2);
0074
0075   sprintf(histo, "ES Trending Fiber Error per 5 mins");
0076   hESFiberErrTrend_ = iBooker.bookProfile(histo, histo, 36, 0.0, 180.0, 100, 0.0, 1.0e6, "s");
0077   hESFiberErrTrend_->setAxisTitle("Elapse time (Minutes)", 1);
0078   hESFiberErrTrend_->setAxisTitle("ES Fiber Err / 5 minutes", 2);
0079
0080   sprintf(histo, "ES Trending SLink CRC Error per hour");
0081   hESSLinkErrTrendHr_ = iBooker.bookProfile(histo, histo, 24, 0.0, 24.0, 100, 0.0, 1.0e6, "s");
0082   hESSLinkErrTrendHr_->setAxisTitle("Elapse time (Hours)", 1);
0083   hESSLinkErrTrendHr_->setAxisTitle("ES SLink CRC Err / hour", 2);
0084
0085   sprintf(histo, "ES Trending Fiber Error per hour");
0086   hESFiberErrTrendHr_ = iBooker.bookProfile(histo, histo, 24, 0.0, 24.0, 100, 0.0, 1.0e6, "s");
0087   hESFiberErrTrendHr_->setAxisTitle("Elapse time (Hours)", 1);
0088   hESFiberErrTrendHr_->setAxisTitle("ES Fiber Err / hour", 2);
0089 }
0090
0091 void ESTrendTask::analyze(const Event& e, const EventSetup& c) {
0092   ievt_++;
0093
0094   // Collect time information
0095   updateTime(e);
0096
0097   long int diff_current_start = current_time_ - start_time_;
0098   long int diff_last_start = last_time_ - start_time_;
0099   //LogInfo("ESTrendTask") << "time difference is negative in " << ievt_ << " events\n"
0100   //<< "\tcurrent - start time = " << diff_current_start
0101   //<< ", \tlast - start time = " << diff_last_start << endl;
0102
0103   //  std::cout << "current_time : " << current_time_ << ", diff : " << diff_current_start << std::endl;
0104
0105   // Calculate time interval and bin width
0106   //  int minuteBinWidth = int(nBasicClusterMinutely_->getTProfile()->GetXaxis()->GetBinWidth(1));
0107   int minuteBinWidth = 5;
0108   long int minuteBinDiff = diff_current_start / 60 / minuteBinWidth - diff_last_start / 60 / minuteBinWidth;
0109   long int minuteDiff = (current_time_ - last_time_) / 60;
0110
0111   //  int hourBinWidth = int(nBasicClusterHourly_->getTProfile()->GetXaxis()->GetBinWidth(1));
0112   int hourBinWidth = 1;
0113   long int hourBinDiff = diff_current_start / 3600 / hourBinWidth - diff_last_start / 3600 / hourBinWidth;
0114   long int hourDiff = (current_time_ - last_time_) / 3600;
0115
0116   if (minuteDiff >= minuteBinWidth) {
0117     while (minuteDiff >= minuteBinWidth)
0118       minuteDiff -= minuteBinWidth;
0119   }
0120   if (hourDiff >= hourBinWidth) {
0121     while (hourDiff >= hourBinWidth)
0122       hourDiff -= hourBinWidth;
0123   }
0124
0125   // ES DCC
0126   Int_t slinkCRCErr = 0;
0127   Int_t fiberErr = 0;
0128   vector<int> fiberStatus;
0129   Handle<ESRawDataCollection> dccs;
0130   if (e.getByToken(dccCollections_, dccs)) {
0131     for (ESRawDataCollection::const_iterator dccItr = dccs->begin(); dccItr != dccs->end(); ++dccItr) {
0132       const ESDCCHeaderBlock& dcc = (*dccItr);
0133
0134       if (dcc.getDCCErrors() == 101)
0135         slinkCRCErr++;
0136
0137       fiberStatus = dcc.getFEChannelStatus();
0138       for (unsigned int i = 0; i < fiberStatus.size(); ++i) {
0139         if (fiberStatus[i] == 4 || fiberStatus[i] == 8 || fiberStatus[i] == 10 || fiberStatus[i] == 11 ||
0140             fiberStatus[i] == 12)
0141           fiberErr++;
0142       }
0143     }
0144   }
0145
0146   shift2Right(hESSLinkErrTrend_->getTProfile(), minuteBinDiff);
0147   hESSLinkErrTrend_->Fill(minuteDiff, slinkCRCErr);
0148
0149   shift2Right(hESFiberErrTrend_->getTProfile(), minuteBinDiff);
0150   hESFiberErrTrend_->Fill(minuteDiff, fiberErr);
0151
0152   shift2Right(hESSLinkErrTrendHr_->getTProfile(), hourBinDiff);
0153   hESSLinkErrTrendHr_->Fill(hourDiff, slinkCRCErr);
0154
0155   shift2Right(hESFiberErrTrendHr_->getTProfile(), hourBinDiff);
0156   hESFiberErrTrendHr_->Fill(hourDiff, fiberErr);
0157
0158   // ES RecHits
0159   int zside, plane;
0160   int nrh[2][2];
0161   for (int i = 0; i < 2; i++)
0162     for (int j = 0; j < 2; j++) {
0163       nrh[i][j] = 0;
0164     }
0165
0166   Handle<ESRecHitCollection> ESRecHit;
0167   if (e.getByToken(rechittoken_, ESRecHit)) {
0168     for (ESRecHitCollection::const_iterator hitItr = ESRecHit->begin(); hitItr != ESRecHit->end(); ++hitItr) {
0169       ESDetId id = ESDetId(hitItr->id());
0170
0171       zside = id.zside();
0172       plane = id.plane();
0173
0174       int i = (zside == 1) ? 0 : 1;
0175       int j = plane - 1;
0176
0177       nrh[i][j]++;
0178     }
0179   } else {
0180     LogWarning("ESTrendTask") << "RecHitCollection not available";
0181   }
0182
0183   for (int i = 0; i < 2; ++i)
0184     for (int j = 0; j < 2; ++j) {
0185       shift2Right(hESRecHitTrend_[i][j]->getTProfile(), minuteBinDiff);
0186       hESRecHitTrend_[i][j]->Fill(minuteDiff, nrh[i][j] / (1072 * 32.));
0187
0188       shift2Right(hESRecHitTrendHr_[i][j]->getTProfile(), hourBinDiff);
0189       hESRecHitTrendHr_[i][j]->Fill(hourDiff, nrh[i][j] / (1072 * 32.));
0190     }
0191 }
0192
0193 void ESTrendTask::updateTime(const edm::Event& e) {
0194   last_time_ = current_time_;
0195   current_time_ = e.time().unixTime();
0196 }
0197
0198 void ESTrendTask::shift2Right(TProfile* p, int bins) {
0199   if (bins <= 0)
0200     return;
0201
0202   if (!p->GetSumw2())
0203     p->Sumw2();
0204   int nBins = p->GetXaxis()->GetNbins();
0205
0206   // by shifting n bin to the right, the number of entries are
0207   // reduced by the number in n bins including the overflow bin.
0208   double nentries = p->GetEntries();
0209   for (int i = 0; i < bins; i++)
0210     nentries -= p->GetBinEntries(nBins + 1 - bins);
0211   p->SetEntries(nentries);
0212
0213   // the last bin goes to overflow
0214   // each bin moves to the right
0215
0216   TArrayD* sumw2 = p->GetSumw2();
0217
0218   for (int i = nBins + 1; i > bins; i--) {
0219     // GetBinContent return binContent/binEntries
0220     p->SetBinContent(i, p->GetBinContent(i - bins) * p->GetBinEntries(i - bins));
0221     p->SetBinEntries(i, p->GetBinEntries(i - bins));
0222     sumw2->SetAt(sumw2->GetAt(i - bins), i);
0223   }
0224 }
0225
0226 void ESTrendTask::shift2Left(TProfile* p, int bins) {
0227   if (bins <= 0)
0228     return;
0229
0230   if (!p->GetSumw2())
0231     p->Sumw2();
0232   int nBins = p->GetXaxis()->GetNbins();
0233
0234   // by shifting n bin to the left, the number of entries are
0235   // reduced by the number in n bins including the underflow bin.
0236   double nentries = p->GetEntries();
0237   for (int i = 0; i < bins; i++)
0238     nentries -= p->GetBinEntries(i);
0239   p->SetEntries(nentries);
0240
0241   // the first bin goes to underflow
0242   // each bin moves to the right
0243
0244   TArrayD* sumw2 = p->GetSumw2();
0245
0246   for (int i = 0; i <= nBins + 1 - bins; i++) {
0247     // GetBinContent return binContent/binEntries
0248     p->SetBinContent(i, p->GetBinContent(i + bins) * p->GetBinEntries(i + bins));
0249     p->SetBinEntries(i, p->GetBinEntries(i + bins));
0250     sumw2->SetAt(sumw2->GetAt(i + bins), i);
0251   }
0252 }
0253
0254 DEFINE_FWK_MODULE(ESTrendTask);