File indexing completed on 2021-02-14 13:11:28
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
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::FED_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 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
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
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
0110 if (bin_contents.size() < 100) {
0111 anal->addErrorCode(sistrip::numberOfBins_);
0112 return;
0113 }
0114
0115
0116 float range = max - min;
0117 float threshold = min + range / 2.;
0118 if (range < 50.) {
0119 anal->addErrorCode(sistrip::smallDataRange_);
0120 return;
0121 }
0122
0123
0124
0125
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
0138
0139
0140
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
0152
0153 if ((tickmark - baseline) < 50.) {
0154 anal->addErrorCode(sistrip::smallDataRange_);
0155 return;
0156 }
0157
0158
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
0179
0180
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
0188
0189 }
0190 }
0191 }
0192
0193
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
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_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
0212 if (temp_deriv > max_deriv) {
0213 max_deriv = temp_deriv;
0214 deriv_bin = bin;
0215 }
0216
0217
0218 if (ii > 10 && ii < 40 && bin_entries[bin] && bin_contents[bin] < baseline + 5 * baseline_rms) {
0219 valid = false;
0220 }
0221 }
0222
0223
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
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 }