SiStripCommissioningSources/src/PedestalsTask.cc

0001 #include "DQM/SiStripCommissioningSources/interface/PedestalsTask.h"
0002 #include "DataFormats/SiStripCommon/interface/SiStripConstants.h"
0003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
0004 #include "DQMServices/Core/interface/DQMStore.h"
0005 #include "DQM/SiStripCommon/interface/ExtractTObject.h"
0006 #include "DQM/SiStripCommon/interface/UpdateTProfile.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include <algorithm>
0009 #include <cmath>
0010
0011 using namespace sistrip;
0012
0013 // -----------------------------------------------------------------------------
0014 //
0015 PedestalsTask::PedestalsTask(DQMStore* dqm, const FedChannelConnection& conn)
0016     : CommissioningTask(dqm, conn, "PedestalsTask"), peds_(), cm_() {
0017   LogTrace(mlDqmSource_) << "[PedestalsTask::" << __func__ << "]"
0018                          << " Constructing object...";
0019 }
0020
0021 // -----------------------------------------------------------------------------
0022 //
0023 PedestalsTask::~PedestalsTask() {
0024   LogTrace(mlDqmSource_) << "[PedestalsTask::" << __func__ << "]"
0025                          << " Destructing object...";
0026 }
0027
0028 // -----------------------------------------------------------------------------
0029 //
0030 void PedestalsTask::book() {
0031   LogTrace(mlDqmSource_) << "[PedestalsTask::" << __func__ << "]";
0032
0033   uint16_t nbins;
0034   std::string title;
0035   std::string extra_info;
0036   peds_.resize(2);
0037   nbins = 256;
0038
0039   // Pedestals histogram
0040   extra_info = sistrip::extrainfo::pedestals_;
0041   peds_[0].isProfile_ = true;
0042
0043   title = SiStripHistoTitle(sistrip::EXPERT_HISTO,
0044                             sistrip::PEDESTALS,
0045                             sistrip::FED_KEY,
0046                             fedKey(),
0047                             sistrip::LLD_CHAN,
0048                             connection().lldChannel(),
0049                             extra_info)
0050               .title();
0051
0052   peds_[0].histo(dqm()->bookProfile(title, title, nbins, -0.5, nbins * 1. - 0.5, 1025, 0., 1025.));
0053
0054   peds_[0].vNumOfEntries_.resize(nbins, 0);
0055   peds_[0].vSumOfContents_.resize(nbins, 0);
0056   peds_[0].vSumOfSquares_.resize(nbins, 0);
0057
0058   // Noise histogram
0059   extra_info = sistrip::extrainfo::noise_;
0060   peds_[1].isProfile_ = true;
0061
0062   title = SiStripHistoTitle(sistrip::EXPERT_HISTO,
0063                             sistrip::PEDESTALS,
0064                             sistrip::FED_KEY,
0065                             fedKey(),
0066                             sistrip::LLD_CHAN,
0067                             connection().lldChannel(),
0068                             extra_info)
0069               .title();
0070
0071   peds_[1].histo(dqm()->bookProfile(title, title, nbins, -0.5, nbins * 1. - 0.5, 1025, 0., 1025.));
0072
0073   peds_[1].vNumOfEntries_.resize(nbins, 0);
0074   peds_[1].vSumOfContents_.resize(nbins, 0);
0075   peds_[1].vSumOfSquares_.resize(nbins, 0);
0076
0077   // Common mode histograms
0078   cm_.resize(2);
0079   nbins = 1024;
0080   for (uint16_t iapv = 0; iapv < 2; iapv++) {
0081     title = SiStripHistoTitle(sistrip::EXPERT_HISTO,
0082                               sistrip::PEDESTALS,
0083                               sistrip::FED_KEY,
0084                               fedKey(),
0085                               sistrip::APV,
0086                               connection().i2cAddr(iapv),
0087                               sistrip::extrainfo::commonMode_)
0088                 .title();
0089
0090     cm_[iapv].histo(dqm()->book1D(title, title, nbins, -0.5, nbins * 1. - 0.5));
0091     cm_[iapv].isProfile_ = false;
0092
0093     cm_[iapv].vNumOfEntries_.resize(nbins, 0);
0094     cm_[iapv].vNumOfEntries_.resize(nbins, 0);
0095   }
0096 }
0097
0098 // -----------------------------------------------------------------------------
0099 //
0100 void PedestalsTask::fill(const SiStripEventSummary& summary, const edm::DetSet<SiStripRawDigi>& digis) {
0101   if (digis.data.size() != peds_[0].vNumOfEntries_.size()) {
0102     edm::LogWarning(mlDqmSource_) << "[PedestalsTask::" << __func__ << "]"
0103                                   << " Unexpected number of digis: " << digis.data.size();
0104     return;
0105   }
0106
0107   // Check number of digis
0108   uint16_t nbins = peds_[0].vNumOfEntries_.size();
0109   if (digis.data.size() < nbins) {
0110     nbins = digis.data.size();
0111   }
0112
0113   //@@ Inefficient!!!
0114   uint16_t napvs = nbins / 128;
0115   std::vector<uint32_t> cm;
0116   cm.resize(napvs, 0);
0117
0118   // Calc common mode for both APVs
0119   std::vector<uint16_t> adc;
0120   for (uint16_t iapv = 0; iapv < napvs; iapv++) {
0121     adc.clear();
0122     adc.reserve(128);
0123     for (uint16_t ibin = 0; ibin < 128; ibin++) {
0124       if ((iapv * 128) + ibin < nbins) {
0125         adc.push_back(digis.data[(iapv * 128) + ibin].adc());  //@@ VIRGIN RAW DATA (MUX, APV READOUT)
0126       }
0127     }
0128     sort(adc.begin(), adc.end());
0129     uint16_t index = adc.size() % 2 ? adc.size() / 2 : adc.size() / 2 - 1;
0130     if (!adc.empty()) {
0131       cm[iapv] = static_cast<uint32_t>(adc[index]);
0132     }
0133   }
0134
0135   for (uint16_t ibin = 0; ibin < nbins; ibin++) {
0136     updateHistoSet(peds_[0], ibin, digis.data[ibin].adc());  // peds and raw noise
0137     float diff = static_cast<float>(digis.data[ibin].adc()) - static_cast<float>(cm[ibin / 128]);
0138     updateHistoSet(peds_[1], ibin, diff);  // residuals and real noise
0139   }
0140
0141   if (cm.size() < cm_.size()) {
0142     edm::LogWarning(mlDqmSource_) << "[PedestalsTask::" << __func__ << "]"
0143                                   << " Fewer CM values than expected: " << cm.size();
0144   }
0145
0146   updateHistoSet(cm_[0], cm[0]);
0147   updateHistoSet(cm_[1], cm[1]);
0148 }
0149
0150 // -----------------------------------------------------------------------------
0151 //
0152 void PedestalsTask::update() {
0153   // Pedestals
0154   updateHistoSet(peds_[0]);
0155
0156   // Noise (cannot use HistoSet directly, as want to plot noise as "contents", not "error")
0157   TProfile* histo = ExtractTObject<TProfile>().extract(peds_[1].histo());
0158   for (uint16_t ii = 0; ii < peds_[1].vNumOfEntries_.size(); ++ii) {
0159     float mean = 0.;
0160     float spread = 0.;
0161     float entries = peds_[1].vNumOfEntries_[ii];
0162     if (entries > 0.) {
0163       mean = peds_[1].vSumOfContents_[ii] / entries;
0164       spread = sqrt(peds_[1].vSumOfSquares_[ii] / entries - mean * mean);
0165     }
0166
0167     float noise = spread;
0168     float error = 0;  // sqrt(entries) / entries;
0169
0170     UpdateTProfile::setBinContent(histo, ii + 1, entries, noise, error);
0171   }
0172
0173   // Common mode
0174   updateHistoSet(cm_[0]);
0175   updateHistoSet(cm_[1]);
0176 }