SiStripCommissioningAnalysis/src/PedestalsAlgorithm.cc

0001 #include "DQM/SiStripCommissioningAnalysis/interface/PedestalsAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/PedestalsAnalysis.h"
0003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
0004 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "TProfile.h"
0007 #include "TH1.h"
0008 #include <iostream>
0009 #include <iomanip>
0010 #include <cmath>
0011
0012 using namespace sistrip;
0013
0014 // ----------------------------------------------------------------------------
0015 //
0016 PedestalsAlgorithm::PedestalsAlgorithm(const edm::ParameterSet& pset, PedestalsAnalysis* const anal)
0017     : CommissioningAlgorithm(anal),
0018       hPeds_(nullptr, ""),
0019       hNoise_(nullptr, ""),
0020       deadStripMax_(pset.getParameter<double>("DeadStripMax")),
0021       noisyStripMin_(pset.getParameter<double>("NoisyStripMin")) {
0022   LogDebug(mlCommissioning_) << "[PedestalsAlgorithm::" << __func__ << "]"
0023                              << " Set maximum noise deviation for dead strip determination to: " << deadStripMax_;
0024   LogDebug(mlCommissioning_) << "[PedestalsAlgorithm::" << __func__ << "]"
0025                              << " Set minimal noise deviation for noisy strip determination to: " << noisyStripMin_;
0026 }
0027
0028 // ----------------------------------------------------------------------------
0029 //
0030 void PedestalsAlgorithm::extract(const std::vector<TH1*>& histos) {
0031   if (!anal()) {
0032     edm::LogWarning(mlCommissioning_) << "[PedestalsAlgorithm::" << __func__ << "]"
0033                                       << " NULL pointer to Analysis object!";
0034     return;
0035   }
0036
0037   // Check number of histograms
0038   if (histos.size() != 2) {
0039     anal()->addErrorCode(sistrip::numberOfHistos_);
0040   }
0041
0042   // Extract FED key from histo title
0043   if (!histos.empty()) {
0044     anal()->fedKey(extractFedKey(histos.front()));
0045   }
0046
0047   // Extract histograms
0048   std::vector<TH1*>::const_iterator ihis = histos.begin();
0049   for (; ihis != histos.end(); ihis++) {
0050     // Check for NULL pointer
0051     if (!(*ihis)) {
0052       continue;
0053     }
0054
0055     // Check run type
0056     SiStripHistoTitle title((*ihis)->GetName());
0057     if (title.runType() != sistrip::PEDESTALS) {
0058       anal()->addErrorCode(sistrip::unexpectedTask_);
0059       continue;
0060     }
0061
0062     // Extract peds and noise histos (check for legacy names first!)
0063     if (title.extraInfo().find(sistrip::extrainfo::pedsAndRawNoise_) != std::string::npos) {
0064       hPeds_.first = *ihis;
0065       hPeds_.second = (*ihis)->GetName();
0066       PedestalsAnalysis* a = dynamic_cast<PedestalsAnalysis*>(const_cast<CommissioningAnalysis*>(anal()));
0067       if (a) {
0068         a->legacy_ = true;
0069       }
0070     } else if (title.extraInfo().find(sistrip::extrainfo::pedsAndCmSubNoise_) != std::string::npos) {
0071       hNoise_.first = *ihis;
0072       hNoise_.second = (*ihis)->GetName();
0073       PedestalsAnalysis* a = dynamic_cast<PedestalsAnalysis*>(const_cast<CommissioningAnalysis*>(anal()));
0074       if (a) {
0075         a->legacy_ = true;
0076       }
0077     } else if (title.extraInfo().find(sistrip::extrainfo::pedestals_) != std::string::npos) {
0078       hPeds_.first = *ihis;
0079       hPeds_.second = (*ihis)->GetName();
0080     } else if (title.extraInfo().find(sistrip::extrainfo::noise_) != std::string::npos) {
0081       hNoise_.first = *ihis;
0082       hNoise_.second = (*ihis)->GetName();
0083     } else if (title.extraInfo().find(sistrip::extrainfo::commonMode_) != std::string::npos) {
0084       //@@ something here for CM plots?
0085     } else {
0086       anal()->addErrorCode(sistrip::unexpectedExtraInfo_);
0087     }
0088   }
0089 }
0090
0091 // -----------------------------------------------------------------------------
0092 //
0093 void PedestalsAlgorithm::analyse() {
0094   if (!anal()) {
0095     edm::LogWarning(mlCommissioning_) << "[PedestalsAlgorithm::" << __func__ << "]"
0096                                       << " NULL pointer to base Analysis object!";
0097     return;
0098   }
0099
0100   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0101   PedestalsAnalysis* anal = dynamic_cast<PedestalsAnalysis*>(tmp);
0102   if (!anal) {
0103     edm::LogWarning(mlCommissioning_) << "[PedestalsAlgorithm::" << __func__ << "]"
0104                                       << " NULL pointer to derived Analysis object!";
0105     return;
0106   }
0107
0108   if (!hPeds_.first) {
0109     anal->addErrorCode(sistrip::nullPtr_);
0110     return;
0111   }
0112
0113   if (!hNoise_.first) {
0114     anal->addErrorCode(sistrip::nullPtr_);
0115     return;
0116   }
0117
0118   TProfile* peds_histo = dynamic_cast<TProfile*>(hPeds_.first);
0119   TProfile* noise_histo = dynamic_cast<TProfile*>(hNoise_.first);
0120
0121   if (!peds_histo) {
0122     anal->addErrorCode(sistrip::nullPtr_);
0123     return;
0124   }
0125
0126   if (!noise_histo) {
0127     anal->addErrorCode(sistrip::nullPtr_);
0128     return;
0129   }
0130
0131   if (peds_histo->GetNbinsX() != 256) {
0132     anal->addErrorCode(sistrip::numberOfBins_);
0133     return;
0134   }
0135
0136   if (noise_histo->GetNbinsX() != 256) {
0137     anal->addErrorCode(sistrip::numberOfBins_);
0138     return;
0139   }
0140
0141   // Iterate through APVs
0142   for (uint16_t iapv = 0; iapv < 2; iapv++) {
0143     // Used to calc mean and rms for peds and noise
0144     float p_sum = 0., p_sum2 = 0., p_max = -1. * sistrip::invalid_, p_min = sistrip::invalid_;
0145     float n_sum = 0., n_sum2 = 0., n_max = -1. * sistrip::invalid_, n_min = sistrip::invalid_;
0146     float r_sum = 0., r_sum2 = 0., r_max = -1. * sistrip::invalid_, r_min = sistrip::invalid_;
0147
0148     // Iterate through strips of APV
0149     for (uint16_t istr = 0; istr < 128; istr++) {
0150       uint16_t strip = iapv * 128 + istr;
0151
0152       // Pedestals and raw noise
0153       if (peds_histo) {
0154         if (peds_histo->GetBinEntries(strip + 1)) {
0155           anal->peds_[iapv][istr] = peds_histo->GetBinContent(strip + 1);
0156           p_sum += anal->peds_[iapv][istr];
0157           p_sum2 += (anal->peds_[iapv][istr] * anal->peds_[iapv][istr]);
0158           if (anal->peds_[iapv][istr] > p_max) {
0159             p_max = anal->peds_[iapv][istr];
0160           }
0161           if (anal->peds_[iapv][istr] < p_min) {
0162             p_min = anal->peds_[iapv][istr];
0163           }
0164
0165           anal->raw_[iapv][istr] = peds_histo->GetBinError(strip + 1);
0166           r_sum += anal->raw_[iapv][istr];
0167           r_sum2 += (anal->raw_[iapv][istr] * anal->raw_[iapv][istr]);
0168           if (anal->raw_[iapv][istr] > r_max) {
0169             r_max = anal->raw_[iapv][istr];
0170           }
0171           if (anal->raw_[iapv][istr] < r_min) {
0172             r_min = anal->raw_[iapv][istr];
0173           }
0174         }
0175       }
0176
0177       // Noise
0178       if (noise_histo) {
0179         if (noise_histo->GetBinEntries(strip + 1)) {
0180           anal->noise_[iapv][istr] = noise_histo->GetBinContent(strip + 1);
0181           n_sum += anal->noise_[iapv][istr];
0182           n_sum2 += (anal->noise_[iapv][istr] * anal->noise_[iapv][istr]);
0183           if (anal->noise_[iapv][istr] > n_max) {
0184             n_max = anal->noise_[iapv][istr];
0185           }
0186           if (anal->noise_[iapv][istr] < n_min) {
0187             n_min = anal->noise_[iapv][istr];
0188           }
0189         }
0190       }
0191
0192     }  // strip loop
0193
0194     // Calc mean and rms for peds
0195     if (!anal->peds_[iapv].empty()) {
0196       p_sum /= static_cast<float>(anal->peds_[iapv].size());
0197       p_sum2 /= static_cast<float>(anal->peds_[iapv].size());
0198       anal->pedsMean_[iapv] = p_sum;
0199       anal->pedsSpread_[iapv] = sqrt(fabs(p_sum2 - p_sum * p_sum));
0200     }
0201
0202     // Calc mean and rms for noise
0203     if (!anal->noise_[iapv].empty()) {
0204       n_sum /= static_cast<float>(anal->noise_[iapv].size());
0205       n_sum2 /= static_cast<float>(anal->noise_[iapv].size());
0206       anal->noiseMean_[iapv] = n_sum;
0207       anal->noiseSpread_[iapv] = sqrt(fabs(n_sum2 - n_sum * n_sum));
0208     }
0209
0210     // Calc mean and rms for raw noise
0211     if (!anal->raw_[iapv].empty()) {
0212       r_sum /= static_cast<float>(anal->raw_[iapv].size());
0213       r_sum2 /= static_cast<float>(anal->raw_[iapv].size());
0214       anal->rawMean_[iapv] = r_sum;
0215       anal->rawSpread_[iapv] = sqrt(fabs(r_sum2 - r_sum * r_sum));
0216     }
0217
0218     // Set max and min values for peds, noise and raw noise
0219     if (p_max > -1. * sistrip::maximum_) {
0220       anal->pedsMax_[iapv] = p_max;
0221     }
0222     if (p_min < 1. * sistrip::maximum_) {
0223       anal->pedsMin_[iapv] = p_min;
0224     }
0225     if (n_max > -1. * sistrip::maximum_) {
0226       anal->noiseMax_[iapv] = n_max;
0227     }
0228     if (n_min < 1. * sistrip::maximum_) {
0229       anal->noiseMin_[iapv] = n_min;
0230     }
0231     if (r_max > -1. * sistrip::maximum_) {
0232       anal->rawMax_[iapv] = r_max;
0233     }
0234     if (r_min < 1. * sistrip::maximum_) {
0235       anal->rawMin_[iapv] = r_min;
0236     }
0237
0238     // Set dead and noisy strips
0239     for (uint16_t istr = 0; istr < 128; istr++) {
0240       if (anal->noiseMin_[iapv] > sistrip::maximum_ || anal->noiseMax_[iapv] > sistrip::maximum_) {
0241         continue;
0242       }
0243       if (anal->noise_[iapv][istr] < (anal->noiseMean_[iapv] - deadStripMax_ * anal->noiseSpread_[iapv])) {
0244         anal->dead_[iapv].push_back(istr);
0245       } else if (anal->noise_[iapv][istr] > (anal->noiseMean_[iapv] + noisyStripMin_ * anal->noiseSpread_[iapv])) {
0246         anal->noisy_[iapv].push_back(istr);
0247       }
0248     }
0249
0250   }  // apv loop
0251 }