EcalCalibAlgos/src/ECALpedestalPCLHarvester.cc

0001
0002 // system include files
0003 #include <memory>
0004
0005 // user include files
0006 #include "Calibration/EcalCalibAlgos/interface/ECALpedestalPCLHarvester.h"
0007 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "CondFormats/EcalObjects/interface/EcalChannelStatusCode.h"
0010 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0011 #include <string>
0012
0013 ECALpedestalPCLHarvester::ECALpedestalPCLHarvester(const edm::ParameterSet& ps)
0014     : minEntries_(ps.getParameter<int>("MinEntries")),
0015       checkAnomalies_(ps.getParameter<bool>("checkAnomalies")),
0016       nSigma_(ps.getParameter<double>("nSigma")),
0017       thresholdAnomalies_(ps.getParameter<double>("thresholdAnomalies")),
0018       dqmDir_(ps.getParameter<std::string>("dqmDir")),
0019       labelG6G1_(ps.getParameter<std::string>("labelG6G1")),
0020       threshDiffEB_(ps.getParameter<double>("threshDiffEB")),
0021       threshDiffEE_(ps.getParameter<double>("threshDiffEE")),
0022       threshChannelsAnalyzed_(ps.getParameter<double>("threshChannelsAnalyzed")),
0023       channelsStatusToken_(esConsumes<edm::Transition::EndRun>()),
0024       pedestalsToken_(esConsumes<edm::Transition::EndRun>()),
0025       g6g1PedestalsToken_(esConsumes<edm::Transition::EndRun>(edm::ESInputTag("", labelG6G1_))),
0026       currentPedestals_(nullptr),
0027       g6g1Pedestals_(nullptr),
0028       channelStatus_(nullptr) {
0029   chStatusToExclude_ = StringToEnumValue<EcalChannelStatusCode::Code>(
0030       ps.getParameter<std::vector<std::string> >("ChannelStatusToExclude"));
0031 }
0032
0033 void ECALpedestalPCLHarvester::dqmEndJob(DQMStore::IBooker& ibooker_, DQMStore::IGetter& igetter_) {
0034   // calculate pedestals and fill db record
0035   EcalPedestals pedestals;
0036   float kBarrelSize = 61200;
0037   float kEndcapSize = 2 * 7324;
0038   float skipped_channels_EB = 0;
0039   float skipped_channels_EE = 0;
0040
0041   for (uint16_t i = 0; i < EBDetId::kSizeForDenseIndexing; ++i) {
0042     std::string hname = dqmDir_ + "/EB/" + std::to_string(int(i / 100)) + "/eb_" + std::to_string(i);
0043     MonitorElement* ch = igetter_.get(hname);
0044     if (ch == nullptr) {
0045       edm::LogWarning("MissingMonitorElement") << "failed to find MonitorElement " << hname;
0046       entriesEB_[i] = 0;
0047       continue;
0048     }
0049     double mean = ch->getMean();
0050     double rms = ch->getRMS();
0051     entriesEB_[i] = ch->getEntries();
0052
0053     DetId id = EBDetId::detIdFromDenseIndex(i);
0054     EcalPedestal ped;
0055     EcalPedestal oldped = *currentPedestals_->find(id.rawId());
0056     EcalPedestal g6g1ped = *g6g1Pedestals_->find(id.rawId());
0057
0058     ped.mean_x12 = mean;
0059     ped.rms_x12 = rms;
0060
0061     float diff = std::abs(mean - oldped.mean_x12);
0062
0063     // if bad channel or low stat skip or the difference is too large wrt to previous record
0064     if (ch->getEntries() < minEntries_ || !checkStatusCode(id) || diff > threshDiffEB_) {
0065       ped.mean_x12 = oldped.mean_x12;
0066       ped.rms_x12 = oldped.rms_x12;
0067
0068       skipped_channels_EB++;
0069     }
0070
0071     // copy g6 and g1 from the corressponding record
0072     ped.mean_x6 = g6g1ped.mean_x6;
0073     ped.rms_x6 = g6g1ped.rms_x6;
0074     ped.mean_x1 = g6g1ped.mean_x1;
0075     ped.rms_x1 = g6g1ped.rms_x1;
0076
0077     pedestals.setValue(id.rawId(), ped);
0078   }
0079
0080   for (uint16_t i = 0; i < EEDetId::kSizeForDenseIndexing; ++i) {
0081     std::string hname = dqmDir_ + "/EE/" + std::to_string(int(i / 100)) + "/ee_" + std::to_string(i);
0082
0083     MonitorElement* ch = igetter_.get(hname);
0084     if (ch == nullptr) {
0085       edm::LogWarning("MissingMonitorElement") << "failed to find MonitorElement " << hname;
0086       entriesEE_[i] = 0;
0087       continue;
0088     }
0089     double mean = ch->getMean();
0090     double rms = ch->getRMS();
0091     entriesEE_[i] = ch->getEntries();
0092
0093     DetId id = EEDetId::detIdFromDenseIndex(i);
0094     EcalPedestal ped;
0095     EcalPedestal oldped = *currentPedestals_->find(id.rawId());
0096     EcalPedestal g6g1ped = *g6g1Pedestals_->find(id.rawId());
0097
0098     ped.mean_x12 = mean;
0099     ped.rms_x12 = rms;
0100
0101     float diff = std::abs(mean - oldped.mean_x12);
0102
0103     // if bad channel or low stat skip or the difference is too large wrt to previous record
0104     if (ch->getEntries() < minEntries_ || !checkStatusCode(id) || diff > threshDiffEE_) {
0105       ped.mean_x12 = oldped.mean_x12;
0106       ped.rms_x12 = oldped.rms_x12;
0107
0108       skipped_channels_EE++;
0109     }
0110
0111     // copy g6 and g1 pedestals from corresponding record
0112     ped.mean_x6 = g6g1ped.mean_x6;
0113     ped.rms_x6 = g6g1ped.rms_x6;
0114     ped.mean_x1 = g6g1ped.mean_x1;
0115     ped.rms_x1 = g6g1ped.rms_x1;
0116
0117     pedestals.setValue(id.rawId(), ped);
0118   }
0119
0120   bool enough_stat = false;
0121   if ((kBarrelSize - skipped_channels_EB) / kBarrelSize > threshChannelsAnalyzed_ &&
0122       (kEndcapSize - skipped_channels_EE) / kEndcapSize > threshChannelsAnalyzed_) {
0123     enough_stat = true;
0124   }
0125
0126   dqmPlots(pedestals, ibooker_);
0127
0128   // check if there are large variations wrt exisiting pedstals
0129
0130   if (checkAnomalies_) {
0131     if (checkVariation(*currentPedestals_, pedestals)) {
0132       edm::LogError("Large Variations found wrt to old pedestals, no file created");
0133       return;
0134     }
0135   }
0136
0137   if (!enough_stat)
0138     return;
0139
0140   // write out pedestal record
0141   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0142
0143   if (!poolDbService.isAvailable()) {
0144     throw std::runtime_error("PoolDBService required.");
0145   }
0146
0147   poolDbService->writeOneIOV(pedestals, poolDbService->currentTime(), "EcalPedestalsRcd");
0148 }
0149
0150 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0151 void ECALpedestalPCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0152   edm::ParameterSetDescription desc;
0153   desc.setUnknown();
0154   descriptions.addDefault(desc);
0155 }
0156
0157 void ECALpedestalPCLHarvester::endRun(edm::Run const& run, edm::EventSetup const& isetup) {
0158   channelStatus_ = &isetup.getData(channelsStatusToken_);
0159   currentPedestals_ = &isetup.getData(pedestalsToken_);
0160   g6g1Pedestals_ = &isetup.getData(g6g1PedestalsToken_);
0161 }
0162
0163 bool ECALpedestalPCLHarvester::checkStatusCode(const DetId& id) {
0164   EcalChannelStatusMap::const_iterator dbstatusPtr;
0165   dbstatusPtr = channelStatus_->getMap().find(id.rawId());
0166   EcalChannelStatusCode::Code dbstatus = dbstatusPtr->getStatusCode();
0167
0168   std::vector<int>::const_iterator res = std::find(chStatusToExclude_.begin(), chStatusToExclude_.end(), dbstatus);
0169   if (res != chStatusToExclude_.end())
0170     return false;
0171
0172   return true;
0173 }
0174
0175 bool ECALpedestalPCLHarvester::isGood(const DetId& id) {
0176   EcalChannelStatusMap::const_iterator dbstatusPtr;
0177   dbstatusPtr = channelStatus_->getMap().find(id.rawId());
0178   if (dbstatusPtr == channelStatus_->getMap().end())
0179     edm::LogError("Invalid DetId supplied");
0180   EcalChannelStatusCode::Code dbstatus = dbstatusPtr->getStatusCode();
0181   if (dbstatus == 0)
0182     return true;
0183   return false;
0184 }
0185
0186 bool ECALpedestalPCLHarvester::checkVariation(const EcalPedestalsMap& oldPedestals,
0187                                               const EcalPedestalsMap& newPedestals) {
0188   uint32_t nAnomaliesEB = 0;
0189   uint32_t nAnomaliesEE = 0;
0190
0191   for (uint16_t i = 0; i < EBDetId::kSizeForDenseIndexing; ++i) {
0192     DetId id = EBDetId::detIdFromDenseIndex(i);
0193     const EcalPedestal& newped = *newPedestals.find(id.rawId());
0194     const EcalPedestal& oldped = *oldPedestals.find(id.rawId());
0195
0196     if (std::abs(newped.mean_x12 - oldped.mean_x12) > nSigma_ * oldped.rms_x12)
0197       nAnomaliesEB++;
0198   }
0199
0200   for (uint16_t i = 0; i < EEDetId::kSizeForDenseIndexing; ++i) {
0201     DetId id = EEDetId::detIdFromDenseIndex(i);
0202     const EcalPedestal& newped = *newPedestals.find(id.rawId());
0203     const EcalPedestal& oldped = *oldPedestals.find(id.rawId());
0204
0205     if (std::abs(newped.mean_x12 - oldped.mean_x12) > nSigma_ * oldped.rms_x12)
0206       nAnomaliesEE++;
0207   }
0208
0209   if (nAnomaliesEB > thresholdAnomalies_ * EBDetId::kSizeForDenseIndexing ||
0210       nAnomaliesEE > thresholdAnomalies_ * EEDetId::kSizeForDenseIndexing)
0211     return true;
0212
0213   return false;
0214 }
0215
0216 void ECALpedestalPCLHarvester::dqmPlots(const EcalPedestals& newpeds, DQMStore::IBooker& ibooker) {
0217   ibooker.cd();
0218   ibooker.setCurrentFolder(dqmDir_ + "/Summary");
0219
0220   MonitorElement* pmeb = ibooker.book2D("meaneb", "Pedestal Means EB", 360, 1., 361., 171, -85., 86.);
0221   MonitorElement* preb = ibooker.book2D("rmseb", "Pedestal RMS EB ", 360, 1., 361., 171, -85., 86.);
0222
0223   MonitorElement* pmeep = ibooker.book2D("meaneep", "Pedestal Means EEP", 100, 1, 101, 100, 1, 101);
0224   MonitorElement* preep = ibooker.book2D("rmseep", "Pedestal RMS EEP", 100, 1, 101, 100, 1, 101);
0225
0226   MonitorElement* pmeem = ibooker.book2D("meaneem", "Pedestal Means EEM", 100, 1, 101, 100, 1, 101);
0227   MonitorElement* preem = ibooker.book2D("rmseem", "Pedestal RMS EEM", 100, 1, 101, 100, 1, 101);
0228
0229   MonitorElement* pmebd = ibooker.book2D("meanebdiff", "Abs Rel Pedestal Means Diff EB", 360, 1., 361., 171, -85., 86.);
0230   MonitorElement* prebd = ibooker.book2D("rmsebdiff", "Abs Rel Pedestal RMS Diff E ", 360, 1., 361., 171, -85., 86.);
0231
0232   MonitorElement* pmeepd = ibooker.book2D("meaneepdiff", "Abs Rel Pedestal Means Diff EEP", 100, 1, 101, 100, 1, 101);
0233   MonitorElement* preepd = ibooker.book2D("rmseepdiff", "Abs Rel Pedestal RMS Diff EEP", 100, 1, 101, 100, 1, 101);
0234
0235   MonitorElement* pmeemd = ibooker.book2D("meaneemdiff", "Abs Rel Pedestal Means Diff EEM", 100, 1, 101, 100, 1, 101);
0236   MonitorElement* preemd = ibooker.book2D("rmseemdiff", "Abs RelPedestal RMS Diff EEM", 100, 1, 101, 100, 1, 101);
0237
0238   MonitorElement* poeb = ibooker.book2D("occeb", "Occupancy EB", 360, 1., 361., 171, -85., 86.);
0239   MonitorElement* poeep = ibooker.book2D("occeep", "Occupancy EEP", 100, 1, 101, 100, 1, 101);
0240   MonitorElement* poeem = ibooker.book2D("occeem", "Occupancy EEM", 100, 1, 101, 100, 1, 101);
0241
0242   MonitorElement* hdiffeb = ibooker.book1D("diffeb", "Pedestal Differences EB", 100, -2.5, 2.5);
0243   MonitorElement* hdiffee = ibooker.book1D("diffee", "Pedestal Differences EE", 100, -2.5, 2.5);
0244
0245   for (int hash = 0; hash < EBDetId::kSizeForDenseIndexing; ++hash) {
0246     EBDetId di = EBDetId::detIdFromDenseIndex(hash);
0247     float mean = newpeds[di].mean_x12;
0248     float rms = newpeds[di].rms_x12;
0249
0250     float cmean = (*currentPedestals_)[di].mean_x12;
0251     float crms = (*currentPedestals_)[di].rms_x12;
0252
0253     if (!isGood(di))
0254       continue;  // only good channels are plotted
0255
0256     pmeb->Fill(di.iphi(), di.ieta(), mean);
0257     preb->Fill(di.iphi(), di.ieta(), rms);
0258     if (cmean)
0259       pmebd->Fill(di.iphi(), di.ieta(), std::abs(mean - cmean) / cmean);
0260     if (crms)
0261       prebd->Fill(di.iphi(), di.ieta(), std::abs(rms - crms) / crms);
0262     poeb->Fill(di.iphi(), di.ieta(), entriesEB_[hash]);
0263     hdiffeb->Fill(mean - cmean);
0264   }
0265
0266   for (int hash = 0; hash < EEDetId::kSizeForDenseIndexing; ++hash) {
0267     EEDetId di = EEDetId::detIdFromDenseIndex(hash);
0268     float mean = newpeds[di].mean_x12;
0269     float rms = newpeds[di].rms_x12;
0270     float cmean = (*currentPedestals_)[di].mean_x12;
0271     float crms = (*currentPedestals_)[di].rms_x12;
0272
0273     if (!isGood(di))
0274       continue;  // only good channels are plotted
0275
0276     if (di.zside() > 0) {
0277       pmeep->Fill(di.ix(), di.iy(), mean);
0278       preep->Fill(di.ix(), di.iy(), rms);
0279       poeep->Fill(di.ix(), di.iy(), entriesEE_[hash]);
0280       if (cmean)
0281         pmeepd->Fill(di.ix(), di.iy(), std::abs(mean - cmean) / cmean);
0282       if (crms)
0283         preepd->Fill(di.ix(), di.iy(), std::abs(rms - crms) / crms);
0284     } else {
0285       pmeem->Fill(di.ix(), di.iy(), mean);
0286       preem->Fill(di.ix(), di.iy(), rms);
0287       if (cmean)
0288         pmeemd->Fill(di.ix(), di.iy(), std::abs(mean - cmean) / cmean);
0289       poeem->Fill(di.ix(), di.iy(), entriesEE_[hash]);
0290       if (crms)
0291         preemd->Fill(di.ix(), di.iy(), std::abs(rms - crms) / crms);
0292     }
0293     hdiffee->Fill(mean - cmean);
0294   }
0295 }