File indexing completed on 2021-02-14 13:11:27
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
0031 if (histos.size() != 1) {
0032 anal()->addErrorCode(sistrip::numberOfHistos_);
0033 }
0034
0035
0036 if (!histos.empty()) {
0037 anal()->fedKey(extractFedKey(histos.front()));
0038 }
0039
0040
0041 std::vector<TH1*>::const_iterator ihis = histos.begin();
0042 for (; ihis != histos.end(); ihis++) {
0043
0044 if (!(*ihis)) {
0045 continue;
0046 }
0047
0048
0049 SiStripHistoTitle title((*ihis)->GetName());
0050 if (title.runType() != sistrip::APV_TIMING) {
0051 anal()->addErrorCode(sistrip::unexpectedTask_);
0052 continue;
0053 }
0054
0055
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
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
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
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
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
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
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
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
0200 bool valid = true;
0201 for (uint16_t ii = 0; ii < 50; ii++) {
0202 uint16_t bin = iter->first + ii;
0203
0204
0205 float temp = 0.;
0206 if (static_cast<uint32_t>(bin) < 1 || static_cast<uint32_t>(bin + 1) >= nbins) {
0207 valid = false;
0208 anal->addErrorCode(sistrip::incompletePlateau_);
0209 continue;
0210 }
0211 temp = bin_contents[bin + 1] - bin_contents[bin - 1];
0212
0213
0214 if (temp > max_derivative) {
0215 max_derivative = temp;
0216 max_derivative_bin = bin;
0217 }
0218
0219
0220 if (ii > 10 && ii < 40 && bin_entries[bin] && bin_contents[bin] < baseline + 5. * baseline_rms) {
0221 valid = false;
0222 }
0223 }
0224
0225
0226 if (valid) {
0227 found = true;
0228 }
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239 iter++;
0240 }
0241
0242 if (!found) {
0243
0244 max_derivative_bin = sistrip::invalid_;
0245 max_derivative = -1. * sistrip::invalid_;
0246
0247
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
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
0275 bool valid_r = true;
0276 int lowpointcount_r = 0;
0277 const int lowpointallow_r = 25;
0278 for (uint16_t ii_r = 0; ii_r < 50; ii_r++) {
0279 uint16_t bin_r = iter_r->first + ii_r;
0280
0281
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;
0285 anal->addErrorCode(sistrip::incompletePlateau_);
0286 continue;
0287 }
0288 temp_r = bin_contents[bin_r + 1] - bin_contents[bin_r - 1];
0289
0290
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
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
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
0316 anal->addErrorCode(sistrip::rejectedCandidate_);
0317 }
0318
0319 iter_r++;
0320 }
0321 }
0322
0323
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 }