SiStripCommissioningAnalysis/src/FedTimingAlgorithm.cc

0001 #include "DQM/SiStripCommissioningAnalysis/interface/FedTimingAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/FedTimingAnalysis.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 FedTimingAlgorithm::FedTimingAlgorithm(const edm::ParameterSet& pset, FedTimingAnalysis* const anal)
0017     : CommissioningAlgorithm(anal), histo_(nullptr, "") {
0018   ;
0019 }
0020
0021 // ----------------------------------------------------------------------------
0022 //
0023 void FedTimingAlgorithm::extract(const std::vector<TH1*>& histos) {
0024   if (!anal()) {
0025     edm::LogWarning(mlCommissioning_) << "[FedTimingAlgorithm::" << __func__ << "]"
0026                                       << " NULL pointer to Analysis object!";
0027     return;
0028   }
0029
0030   // Check number of histograms
0031   if (histos.size() != 1) {
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 name
0049     SiStripHistoTitle title((*ihis)->GetName());
0050     if (title.runType() != sistrip::FED_TIMING) {
0051       anal()->addErrorCode(sistrip::unexpectedTask_);
0052       continue;
0053     }
0054
0055     // Extract timing histo
0056     histo_.first = *ihis;
0057     histo_.second = (*ihis)->GetName();
0058   }
0059 }
0060
0061 // ----------------------------------------------------------------------------
0062 //
0063 void FedTimingAlgorithm::analyse() {
0064   if (!anal()) {
0065     edm::LogWarning(mlCommissioning_) << "[FedTimingAlgorithm::" << __func__ << "]"
0066                                       << " NULL pointer to base Analysis object!";
0067     return;
0068   }
0069
0070   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0071   FedTimingAnalysis* anal = dynamic_cast<FedTimingAnalysis*>(tmp);
0072   if (!anal) {
0073     edm::LogWarning(mlCommissioning_) << "[FedTimingAlgorithm::" << __func__ << "]"
0074                                       << " NULL pointer to derived Analysis object!";
0075     return;
0076   }
0077
0078   if (!histo_.first) {
0079     anal->addErrorCode(sistrip::nullPtr_);
0080     return;
0081   }
0082
0083   // Transfer histogram contents/errors/stats to containers
0084   uint16_t non_zero = 0;
0085   float max = -1.e9;
0086   float min = 1.e9;
0087   uint16_t nbins = static_cast<uint16_t>(histo_.first->GetNbinsX());
0088   std::vector<float> bin_contents;
0089   std::vector<float> bin_errors;
0090   std::vector<float> bin_entries;
0091   bin_contents.reserve(nbins);
0092   bin_errors.reserve(nbins);
0093   bin_entries.reserve(nbins);
0094   for (uint16_t ibin = 0; ibin < nbins; ibin++) {
0095     bin_contents.push_back(histo_.first->GetBinContent(ibin + 1));
0096     bin_errors.push_back(histo_.first->GetBinError(ibin + 1));
0097     //bin_entries.push_back( histo_.first->GetBinEntries(ibin+1) );
0098     if (bin_entries[ibin]) {
0099       if (bin_contents[ibin] > max) {
0100         max = bin_contents[ibin];
0101       }
0102       if (bin_contents[ibin] < min) {
0103         min = bin_contents[ibin];
0104       }
0105       non_zero++;
0106     }
0107   }
0108
0109   //LogTrace(mlCommissioning_) << " Number of bins with non-zero entries: " << non_zero;
0110   if (bin_contents.size() < 100) {
0111     anal->addErrorCode(sistrip::numberOfBins_);
0112     return;
0113   }
0114
0115   // Calculate range (max-min) and threshold level (range/2)
0116   float range = max - min;
0117   float threshold = min + range / 2.;
0118   if (range < 50.) {
0119     anal->addErrorCode(sistrip::smallDataRange_);
0120     return;
0121   }
0122   //LogTrace(mlCommissioning_) << " ADC samples: max/min/range/threshold: "
0123   //<< max << "/" << min << "/" << range << "/" << threshold;
0124
0125   // Associate samples with either "tick mark" or "baseline"
0126   std::vector<float> tick;
0127   std::vector<float> base;
0128   for (uint16_t ibin = 0; ibin < nbins; ibin++) {
0129     if (bin_entries[ibin]) {
0130       if (bin_contents[ibin] < threshold) {
0131         base.push_back(bin_contents[ibin]);
0132       } else {
0133         tick.push_back(bin_contents[ibin]);
0134       }
0135     }
0136   }
0137   //LogTrace(mlCommissioning_) << " Number of 'tick mark' samples: " << tick.size()
0138   //<< " Number of 'baseline' samples: " << base.size();
0139
0140   // Find median level of tick mark and baseline
0141   float tickmark = 0.;
0142   float baseline = 0.;
0143   sort(tick.begin(), tick.end());
0144   sort(base.begin(), base.end());
0145   if (!tick.empty()) {
0146     tickmark = tick[tick.size() % 2 ? tick.size() / 2 : tick.size() / 2];
0147   }
0148   if (!base.empty()) {
0149     baseline = base[base.size() % 2 ? base.size() / 2 : base.size() / 2];
0150   }
0151   //LogTrace(mlCommissioning_) << " Tick mark level: " << tickmark << " Baseline level: " << baseline
0152   //<< " Range: " << (tickmark-baseline);
0153   if ((tickmark - baseline) < 50.) {
0154     anal->addErrorCode(sistrip::smallDataRange_);
0155     return;
0156   }
0157
0158   // Find rms spread in "baseline" samples
0159   float mean = 0.;
0160   float mean2 = 0.;
0161   for (uint16_t ibin = 0; ibin < base.size(); ibin++) {
0162     mean += base[ibin];
0163     mean2 += base[ibin] * base[ibin];
0164   }
0165   if (!base.empty()) {
0166     mean = mean / base.size();
0167     mean2 = mean2 / base.size();
0168   } else {
0169     mean = 0.;
0170     mean2 = 0.;
0171   }
0172   float baseline_rms = 0.;
0173   if (mean2 > mean * mean) {
0174     baseline_rms = sqrt(mean2 - mean * mean);
0175   } else {
0176     baseline_rms = 0.;
0177   }
0178   //LogTrace(mlCommissioning_) << " Spread in baseline samples: " << baseline_rms;
0179
0180   // Find rising edges (derivative across two bins > range/2)
0181   std::map<uint16_t, float> edges;
0182   for (uint16_t ibin = 1; ibin < nbins - 1; ibin++) {
0183     if (bin_entries[ibin + 1] && bin_entries[ibin - 1]) {
0184       float derivative = bin_contents[ibin + 1] - bin_contents[ibin - 1];
0185       if (derivative > 5. * baseline_rms) {
0186         edges[ibin] = derivative;
0187         //LogTrace(mlCommissioning_) << " Found edge #" << edges.size() << " at bin " << ibin
0188         //<< " and with derivative " << derivative;
0189       }
0190     }
0191   }
0192
0193   // Iterate through "edges" std::map
0194   bool found = false;
0195   uint16_t deriv_bin = sistrip::invalid_;
0196   float max_deriv = -1. * sistrip::invalid_;
0197   std::map<uint16_t, float>::iterator iter = edges.begin();
0198   while (!found && iter != edges.end()) {
0199     // Iterate through 50 subsequent samples
0200     bool valid = true;
0201     for (uint16_t ii = 0; ii < 50; ii++) {
0202       uint16_t bin = iter->first + ii;
0203
0204       // Calc local derivative
0205       float temp_deriv = 0;
0206       if (static_cast<uint32_t>(bin) < 1 || static_cast<uint32_t>(bin + 1) >= nbins) {
0207         continue;
0208       }
0209       temp_deriv = bin_contents[bin + 1] - bin_contents[bin - 1];
0210
0211       // Store max derivative
0212       if (temp_deriv > max_deriv) {
0213         max_deriv = temp_deriv;
0214         deriv_bin = bin;
0215       }
0216
0217       // Check if samples following edge are all "high"
0218       if (ii > 10 && ii < 40 && bin_entries[bin] && bin_contents[bin] < baseline + 5 * baseline_rms) {
0219         valid = false;
0220       }
0221     }
0222
0223     // Break from loop if tick mark found
0224     if (valid) {
0225       found = true;
0226     } else {
0227       max_deriv = -1. * sistrip::invalid_;
0228       deriv_bin = sistrip::invalid_;
0229       edges.erase(iter);
0230     }
0231
0232     iter++;
0233   }
0234
0235   // Set monitorables (but not PLL coarse and fine here)
0236   if (!edges.empty()) {
0237     anal->time_ = deriv_bin;
0238     anal->error_ = 0.;
0239     anal->base_ = baseline;
0240     anal->peak_ = tickmark;
0241     anal->height_ = tickmark - baseline;
0242   } else {
0243     anal->addErrorCode(sistrip::missingTickMark_);
0244     anal->base_ = baseline;
0245     anal->peak_ = tickmark;
0246     anal->height_ = tickmark - baseline;
0247   }
0248 }