SiStripCommissioningAnalysis/src/DaqScopeModeAlgorithm.cc

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