SiStripCommissioningAnalysis/src/NoiseAlgorithm.cc

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