File indexing completed on 2023-10-25 09:43:14
0001 #include "DQM/SiStripCommissioningAnalysis/interface/SamplingAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/SamplingAnalysis.h"
0003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
0004 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "DQM/SiStripCommissioningAnalysis/interface/SiStripPulseShape.h"
0007 #include "TProfile.h"
0008 #include "TF1.h"
0009 #include <iostream>
0010 #include <sstream>
0011 #include <iomanip>
0012 #include <cmath>
0013
0014 using namespace sistrip;
0015
0016
0017
0018 SamplingAlgorithm::SamplingAlgorithm(const edm::ParameterSet& pset, SamplingAnalysis* const anal, uint32_t latencyCode)
0019 : CommissioningAlgorithm(anal),
0020 histo_(nullptr, ""),
0021 deconv_fitter_(nullptr),
0022 peak_fitterA_(nullptr),
0023 peak_fitterB_(nullptr),
0024 latencyCode_(latencyCode),
0025 samp_(nullptr) {
0026 peak_fitterA_ = new TF1("peak_fitterA", fpeak_convoluted, -4800, 0, 5);
0027 peak_fitterA_->SetNpx(2000);
0028 peak_fitterA_->FixParameter(0, 0);
0029 peak_fitterA_->SetParLimits(1, 0, 4800);
0030 peak_fitterA_->SetParLimits(2, 0, 20);
0031 peak_fitterA_->FixParameter(3, 50);
0032 peak_fitterA_->SetParLimits(4, 0, 75);
0033 peak_fitterA_->SetParameters(0., 1250, 10, 50, 10);
0034
0035 peak_fitterB_ = new TF1("peak_fitterB", fpeak_convoluted, -100, 100, 5);
0036 peak_fitterB_->SetNpx(200);
0037 peak_fitterB_->FixParameter(0, 0);
0038 peak_fitterB_->SetParLimits(1, -100, 100);
0039 peak_fitterB_->SetParLimits(2, 0, 20);
0040 peak_fitterB_->FixParameter(3, 50);
0041 peak_fitterB_->SetParLimits(4, 0, 75);
0042 peak_fitterB_->SetParameters(0., -50, 10, 50, 10);
0043
0044 deconv_fitter_ = new TF1("deconv_fitter", fdeconv_convoluted, -50, 50, 5);
0045 deconv_fitter_->SetNpx(1000);
0046 deconv_fitter_->FixParameter(0, 0);
0047 deconv_fitter_->SetParLimits(1, -50, 50);
0048 deconv_fitter_->SetParLimits(2, 0, 200);
0049 deconv_fitter_->SetParLimits(3, 5, 100);
0050 deconv_fitter_->FixParameter(3, 50);
0051 deconv_fitter_->SetParLimits(4, 0, 75);
0052 deconv_fitter_->SetParameters(0., -2.82, 0.96, 50, 20);
0053 }
0054
0055
0056
0057 void SamplingAlgorithm::extract(const std::vector<TH1*>& histos) {
0058 if (!anal()) {
0059 edm::LogWarning(mlCommissioning_) << "[SamplingAlgorithm::" << __func__ << "]"
0060 << " NULL pointer to Analysis object!";
0061 return;
0062 }
0063
0064 CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0065 samp_ = dynamic_cast<SamplingAnalysis*>(tmp);
0066 if (!samp_) {
0067 edm::LogWarning(mlCommissioning_) << "[SamplingAlgorithm::" << __func__ << "]"
0068 << " NULL pointer to derived Analysis object!";
0069 return;
0070 }
0071
0072
0073 if (histos.size() != 1 && histos.size() != 2) {
0074 samp_->addErrorCode(sistrip::numberOfHistos_);
0075 }
0076
0077
0078 if (!histos.empty()) {
0079 samp_->fedKey(extractFedKey(histos.front()));
0080 }
0081
0082
0083 std::vector<TH1*>::const_iterator ihis = histos.begin();
0084 for (; ihis != histos.end(); ihis++) {
0085
0086 if (!(*ihis)) {
0087 edm::LogWarning(mlCommissioning_) << " NULL pointer to histogram!";
0088 continue;
0089 }
0090
0091
0092 SiStripHistoTitle title((*ihis)->GetName());
0093 if (title.runType() != sistrip::APV_LATENCY && title.runType() != sistrip::FINE_DELAY) {
0094 samp_->addErrorCode(sistrip::unexpectedTask_);
0095 continue;
0096 }
0097
0098 samp_->runType_ = title.runType();
0099
0100
0101 samp_->granularity_ = title.granularity();
0102
0103
0104 if (title.extraInfo().find(sistrip::extrainfo::clusterCharge_) != std::string::npos) {
0105 histo_.first = *ihis;
0106 histo_.second = (*ihis)->GetName();
0107 }
0108 }
0109 }
0110
0111
0112
0113 void SamplingAlgorithm::analyse() {
0114 if (!samp_) {
0115 edm::LogWarning(mlCommissioning_) << "[SamplingAlgorithm::" << __func__ << "]"
0116 << " NULL pointer to derived Analysis object!";
0117 return;
0118 }
0119
0120 TProfile* prof = (TProfile*)(histo_.first);
0121 if (!prof) {
0122 edm::LogWarning(mlCommissioning_) << " NULL pointer to histogram!";
0123 return;
0124 }
0125
0126
0127 prof->SetErrorOption(" ");
0128
0129
0130
0131 for (int i = 0; i < prof->GetNbinsX(); ++i) {
0132 if (prof->GetBinEntries(i) > 0)
0133 prof->SetBinError(i, prof->GetBinError(i) / sqrt(prof->GetBinEntries(i)));
0134 }
0135
0136
0137 pruneProfile(prof);
0138
0139
0140 correctBinning(prof);
0141
0142
0143 correctProfile(prof, samp_->sOnCut_);
0144
0145
0146 if (samp_->runType_ == sistrip::APV_LATENCY) {
0147
0148 float max = prof->GetBinCenter(prof->GetMaximumBin());
0149 float ampl = prof->GetMaximum();
0150 peak_fitterA_->SetParameters(0., 50 - max, ampl / 20., 50, 10);
0151
0152
0153 if (prof->Fit(peak_fitterA_, "Q") == 0)
0154 prof->Fit(peak_fitterA_, "QEM");
0155
0156
0157 samp_->max_ = peak_fitterA_->GetMaximumX();
0158 samp_->error_ = peak_fitterA_->GetParError(1);
0159
0160 } else {
0161
0162
0163 float max = prof->GetBinCenter(prof->GetMaximumBin());
0164 float ampl = prof->GetMaximum();
0165 deconv_fitter_->SetParameters(0., -max, ampl / 10., 50, 20);
0166 peak_fitterB_->SetParameters(0., 50 - max, ampl / 20., 50, 10);
0167 if (latencyCode_ & 0x80) {
0168
0169 if (prof->Fit(deconv_fitter_, "Q") == 0)
0170 prof->Fit(deconv_fitter_, "QEM");
0171
0172 samp_->max_ = deconv_fitter_->GetMaximumX();
0173 samp_->error_ = deconv_fitter_->GetParError(1);
0174 } else {
0175
0176 if (prof->Fit(peak_fitterB_, "Q") == 0)
0177 prof->Fit(peak_fitterB_, "QEM");
0178
0179 samp_->max_ = peak_fitterB_->GetMaximumX();
0180 samp_->error_ = peak_fitterB_->GetParError(1);
0181 }
0182 }
0183 }
0184
0185
0186
0187 void SamplingAlgorithm::pruneProfile(TProfile* profile) const {
0188
0189 uint32_t nbins = profile->GetNbinsX();
0190 uint32_t max = 0;
0191 for (uint32_t bin = 1; bin <= nbins; ++bin)
0192 max = max < profile->GetBinEntries(bin) ? uint32_t(profile->GetBinEntries(bin)) : max;
0193
0194 uint32_t min = max / 10;
0195 for (uint32_t bin = 1; bin <= nbins; ++bin)
0196 if (profile->GetBinEntries(bin) < min) {
0197 profile->SetBinContent(bin, 0.);
0198 profile->SetBinError(bin, 0.);
0199 }
0200 }
0201
0202
0203
0204 void SamplingAlgorithm::correctBinning(TProfile* prof) const {
0205 prof->GetXaxis()->SetLimits(prof->GetXaxis()->GetXmin() - prof->GetBinWidth(1) / 2.,
0206 prof->GetXaxis()->GetXmax() - prof->GetBinWidth(1) / 2.);
0207 }
0208
0209
0210
0211 void SamplingAlgorithm::correctProfile(TProfile* profile, float SoNcut) const {
0212 if (!samp_) {
0213 return;
0214 }
0215 uint32_t nbins = profile->GetNbinsX();
0216 float min = samp_->limit(SoNcut);
0217 for (uint32_t bin = 1; bin <= nbins; ++bin)
0218 if (profile->GetBinContent(bin) < min) {
0219 profile->SetBinContent(bin, 0.);
0220 profile->SetBinError(bin, 0.);
0221 profile->SetBinEntries(bin, 0);
0222 } else {
0223 profile->SetBinContent(
0224 bin, profile->GetBinEntries(bin) * samp_->correctMeasurement(profile->GetBinContent(bin), SoNcut));
0225 }
0226 }