SiStripCommissioningAnalysis/src/SamplingAlgorithm.cc

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   // Check
0073   if (histos.size() != 1 && histos.size() != 2) {
0074     samp_->addErrorCode(sistrip::numberOfHistos_);
0075   }
0076
0077   // Extract FED key from histo title
0078   if (!histos.empty()) {
0079     samp_->fedKey(extractFedKey(histos.front()));
0080   }
0081
0082   // Extract
0083   std::vector<TH1*>::const_iterator ihis = histos.begin();
0084   for (; ihis != histos.end(); ihis++) {
0085     // Check pointer
0086     if (!(*ihis)) {
0087       edm::LogWarning(mlCommissioning_) << " NULL pointer to histogram!";
0088       continue;
0089     }
0090
0091     // Check name
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     // Set the mode for later fits
0098     samp_->runType_ = title.runType();
0099
0100     // Set the granularity
0101     samp_->granularity_ = title.granularity();
0102
0103     // Extract timing histo
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   // set the right error mode: rms
0127   prof->SetErrorOption(" ");
0128
0129   //that should not be needed, but it seems histos are stored with error option " " and errors "s" in all cases.
0130   //it MUST be removed if the DQM (?) bug is solved
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   // prune the profile
0137   pruneProfile(prof);
0138
0139   // correct for the binning
0140   correctBinning(prof);
0141
0142   // correct for clustering effects
0143   correctProfile(prof, samp_->sOnCut_);
0144
0145   // fit depending on the mode
0146   if (samp_->runType_ == sistrip::APV_LATENCY) {
0147     // initialize  the fit (overal latency)
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     // fit
0153     if (prof->Fit(peak_fitterA_, "Q") == 0)
0154       prof->Fit(peak_fitterA_, "QEM");
0155
0156     // Set monitorables
0157     samp_->max_ = peak_fitterA_->GetMaximumX();
0158     samp_->error_ = peak_fitterA_->GetParError(1);
0159
0160   } else {  // sistrip::FINE_DELAY
0161
0162     // initialize  the fit (overal latency)
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) {  // deconv mode
0168       // fit
0169       if (prof->Fit(deconv_fitter_, "Q") == 0)
0170         prof->Fit(deconv_fitter_, "QEM");
0171       // Set monitorables
0172       samp_->max_ = deconv_fitter_->GetMaximumX();
0173       samp_->error_ = deconv_fitter_->GetParError(1);
0174     } else {  // peak mode
0175       // fit
0176       if (prof->Fit(peak_fitterB_, "Q") == 0)
0177         prof->Fit(peak_fitterB_, "QEM");
0178       // Set monitorables
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   // loop over bins to find the max stat
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   // loop over bins to clean
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 }