SiStripCommissioningAnalysis/src/ApvTimingAlgorithm.cc

0001 #include "DQM/SiStripCommissioningAnalysis/interface/ApvTimingAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/ApvTimingAnalysis.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 ApvTimingAlgorithm::ApvTimingAlgorithm(const edm::ParameterSet& pset, ApvTimingAnalysis* const anal)
0017     : CommissioningAlgorithm(anal), histo_(nullptr, "") {
0018   ;
0019 }
0020
0021 // ----------------------------------------------------------------------------
0022 //
0023 void ApvTimingAlgorithm::extract(const std::vector<TH1*>& histos) {
0024   if (!anal()) {
0025     edm::LogWarning(mlCommissioning_) << "[ApvTimingAlgorithm::" << __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::APV_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 ApvTimingAlgorithm::analyse() {
0064   if (!anal()) {
0065     edm::LogWarning(mlCommissioning_) << "[ApvTimingAlgorithm::" << __func__ << "]"
0066                                       << " NULL pointer to base Analysis object!";
0067     return;
0068   }
0069
0070   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0071   ApvTimingAnalysis* anal = dynamic_cast<ApvTimingAnalysis*>(tmp);
0072   if (!anal) {
0073     edm::LogWarning(mlCommissioning_) << "[ApvTimingAlgorithm::" << __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   TProfile* histo = dynamic_cast<TProfile*>(histo_.first);
0084   if (!histo) {
0085     anal->addErrorCode(sistrip::nullPtr_);
0086     return;
0087   }
0088
0089   // Transfer histogram contents/errors/stats to containers
0090   uint16_t non_zero = 0;
0091   float max = -1. * sistrip::invalid_;
0092   float min = 1. * sistrip::invalid_;
0093   uint16_t nbins = static_cast<uint16_t>(histo->GetNbinsX());
0094   std::vector<float> bin_contents;
0095   std::vector<float> bin_errors;
0096   std::vector<float> bin_entries;
0097   bin_contents.reserve(nbins);
0098   bin_errors.reserve(nbins);
0099   bin_entries.reserve(nbins);
0100   for (uint16_t ibin = 0; ibin < nbins; ibin++) {
0101     bin_contents.push_back(histo->GetBinContent(ibin + 1));
0102     bin_errors.push_back(histo->GetBinError(ibin + 1));
0103     bin_entries.push_back(histo->GetBinEntries(ibin + 1));
0104     if (bin_entries[ibin]) {
0105       if (bin_contents[ibin] > max) {
0106         max = bin_contents[ibin];
0107       }
0108       if (bin_contents[ibin] < min) {
0109         min = bin_contents[ibin];
0110       }
0111       non_zero++;
0112     }
0113   }
0114   if (bin_contents.size() < 100) {
0115     anal->addErrorCode(sistrip::numberOfBins_);
0116     return;
0117   }
0118
0119   // Calculate range (max-min) and threshold level (range/2)
0120   float threshold = min + (max - min) / 2.;
0121   anal->base_ = min;
0122   anal->peak_ = max;
0123   anal->height_ = max - min;
0124   if (max - min < ApvTimingAnalysis::tickMarkHeightThreshold_) {
0125     anal->addErrorCode(sistrip::smallDataRange_);
0126     return;
0127   }
0128
0129   // Associate samples with either "tick mark" or "baseline"
0130   std::vector<float> tick;
0131   std::vector<float> base;
0132   for (uint16_t ibin = 0; ibin < nbins; ibin++) {
0133     if (bin_entries[ibin]) {
0134       if (bin_contents[ibin] < threshold) {
0135         base.push_back(bin_contents[ibin]);
0136       } else {
0137         tick.push_back(bin_contents[ibin]);
0138       }
0139     }
0140   }
0141
0142   // Find median level of tick mark and baseline
0143   float tickmark = 0.;
0144   float baseline = 0.;
0145   sort(tick.begin(), tick.end());
0146   sort(base.begin(), base.end());
0147   if (!tick.empty()) {
0148     tickmark = tick[tick.size() % 2 ? tick.size() / 2 : tick.size() / 2];
0149   }
0150   if (!base.empty()) {
0151     baseline = base[base.size() % 2 ? base.size() / 2 : base.size() / 2];
0152   }
0153   anal->base_ = baseline;
0154   anal->peak_ = tickmark;
0155   anal->height_ = tickmark - baseline;
0156   if (tickmark - baseline < ApvTimingAnalysis::tickMarkHeightThreshold_) {
0157     anal->addErrorCode(sistrip::smallTickMarkHeight_);
0158     return;
0159   }
0160
0161   // Find rms spread in "baseline" samples
0162   float mean = 0.;
0163   float mean2 = 0.;
0164   for (uint16_t ibin = 0; ibin < base.size(); ibin++) {
0165     mean += base[ibin];
0166     mean2 += base[ibin] * base[ibin];
0167   }
0168   if (!base.empty()) {
0169     mean = mean / base.size();
0170     mean2 = mean2 / base.size();
0171   } else {
0172     mean = 0.;
0173     mean2 = 0.;
0174   }
0175   float baseline_rms = sqrt(fabs(mean2 - mean * mean));
0176
0177   // Find rising edges (derivative across two bins > threshold)
0178   std::map<uint16_t, float> edges;
0179   for (uint16_t ibin = 1; ibin < nbins - 1; ibin++) {
0180     if (bin_entries[ibin + 1] && bin_entries[ibin - 1]) {
0181       float derivative = bin_contents[ibin + 1] - bin_contents[ibin - 1];
0182       if (derivative > 3. * baseline_rms) {
0183         edges[ibin] = derivative;
0184       }
0185     }
0186   }
0187   if (edges.empty()) {
0188     anal->addErrorCode(sistrip::noRisingEdges_);
0189     return;
0190   }
0191
0192   // Iterate through "edges" map
0193   uint16_t max_derivative_bin = sistrip::invalid_;
0194   float max_derivative = -1. * sistrip::invalid_;
0195
0196   bool found = false;
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 = 0.;
0206       if (static_cast<uint32_t>(bin) < 1 || static_cast<uint32_t>(bin + 1) >= nbins) {
0207         valid = false;  //@@ require complete plateau is found within histo
0208         anal->addErrorCode(sistrip::incompletePlateau_);
0209         continue;
0210       }
0211       temp = bin_contents[bin + 1] - bin_contents[bin - 1];
0212
0213       // Store max derivative
0214       if (temp > max_derivative) {
0215         max_derivative = temp;
0216         max_derivative_bin = bin;
0217       }
0218
0219       // Check if samples following edge are all "high"
0220       if (ii > 10 && ii < 40 && bin_entries[bin] && bin_contents[bin] < baseline + 5. * baseline_rms) {
0221         valid = false;
0222       }
0223     }
0224
0225     // Break from loop if tick mark found
0226     if (valid) {
0227       found = true;
0228     }
0229
0230     /*
0231     else {
0232       max_derivative = -1.*sistrip::invalid_;
0233       max_derivative_bin = sistrip::invalid_;
0234       //edges.erase(iter);
0235       anal->addErrorCode(sistrip::rejectedCandidate_);
0236     }
0237     */
0238
0239     iter++;  // next candidate
0240   }
0241
0242   if (!found) {  //Try tick mark recovery
0243
0244     max_derivative_bin = sistrip::invalid_;
0245     max_derivative = -1. * sistrip::invalid_;
0246
0247     // Find rising edges_r (derivative_r across five bins > threshold)
0248     std::map<uint16_t, float> edges_r;
0249     for (uint16_t ibin_r = 1; ibin_r < nbins - 1; ibin_r++) {
0250       if (bin_entries[ibin_r + 4] && bin_entries[ibin_r + 3] && bin_entries[ibin_r + 2] && bin_entries[ibin_r + 1] &&
0251           bin_entries[ibin_r] && bin_entries[ibin_r - 1]) {
0252         float derivative_r = bin_contents[ibin_r + 1] - bin_contents[ibin_r - 1];
0253         float derivative_r1 = bin_contents[ibin_r + 1] - bin_contents[ibin_r];
0254         float derivative_r2 = bin_contents[ibin_r + 2] - bin_contents[ibin_r + 1];
0255         float derivative_r3 = bin_contents[ibin_r + 3] - bin_contents[ibin_r + 2];
0256
0257         if (derivative_r > 3. * baseline_rms && derivative_r1 > 1. * baseline_rms &&
0258             derivative_r2 > 1. * baseline_rms && derivative_r3 > 1. * baseline_rms) {
0259           edges_r[ibin_r] = derivative_r;
0260         }
0261       }
0262     }
0263     if (edges_r.empty()) {
0264       anal->addErrorCode(sistrip::noRisingEdges_);
0265       return;
0266     }
0267
0268     // Iterate through "edges_r" map
0269     float max_derivative_r = -1. * sistrip::invalid_;
0270
0271     bool found_r = false;
0272     std::map<uint16_t, float>::iterator iter_r = edges_r.begin();
0273     while (!found_r && iter_r != edges_r.end()) {
0274       // Iterate through 50 subsequent samples
0275       bool valid_r = true;
0276       int lowpointcount_r = 0;
0277       const int lowpointallow_r = 25;  //Number of points allowed to fall below threshhold w/o invalidating tick mark
0278       for (uint16_t ii_r = 0; ii_r < 50; ii_r++) {
0279         uint16_t bin_r = iter_r->first + ii_r;
0280
0281         // Calc local derivative_r
0282         float temp_r = 0.;
0283         if (static_cast<uint32_t>(bin_r) < 1 || static_cast<uint32_t>(bin_r + 1) >= nbins) {
0284           valid_r = false;  //@@ require complete plateau is found_r within histo
0285           anal->addErrorCode(sistrip::incompletePlateau_);
0286           continue;
0287         }
0288         temp_r = bin_contents[bin_r + 1] - bin_contents[bin_r - 1];
0289
0290         // Store max derivative_r
0291         if (temp_r > max_derivative_r && ii_r < 10) {
0292           max_derivative_r = temp_r;
0293           max_derivative = temp_r;
0294           max_derivative_bin = bin_r;
0295         }
0296
0297         // Check if majority of samples following edge are all "high"
0298         if (ii_r > 10 && ii_r < 40 && bin_entries[bin_r] && bin_contents[bin_r] < baseline + 5. * baseline_rms) {
0299           lowpointcount_r++;
0300           if (lowpointcount_r > lowpointallow_r) {
0301             valid_r = false;
0302           }
0303         }
0304       }
0305
0306       // Break from loop if recovery tick mark found
0307       if (valid_r) {
0308         found_r = true;
0309         found = true;
0310         anal->addErrorCode(sistrip::tickMarkRecovered_);
0311       } else {
0312         max_derivative_r = -1. * sistrip::invalid_;
0313         max_derivative = -1. * sistrip::invalid_;
0314         max_derivative_bin = sistrip::invalid_;
0315         //edges_r.erase(iter_r);
0316         anal->addErrorCode(sistrip::rejectedCandidate_);
0317       }
0318
0319       iter_r++;  // next candidate
0320     }
0321   }  //End tick mark recovery
0322
0323   // Record time monitorable and check tick mark height
0324   if (max_derivative_bin <= sistrip::valid_) {
0325     anal->time_ = max_derivative_bin * 25. / 24.;
0326     if (anal->height_ < ApvTimingAnalysis::tickMarkHeightThreshold_) {
0327       anal->addErrorCode(sistrip::tickMarkBelowThresh_);
0328     }
0329   } else {
0330     anal->addErrorCode(sistrip::missingTickMark_);
0331   }
0332 }