SiStripCommissioningAnalysis/src/ApvLatencyAlgorithm.cc

0001 #include "DQM/SiStripCommissioningAnalysis/interface/ApvLatencyAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/ApvLatencyAnalysis.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 <iostream>
0008 #include <cmath>
0009
0010 using namespace sistrip;
0011
0012 // ----------------------------------------------------------------------------
0013 //
0014 ApvLatencyAlgorithm::ApvLatencyAlgorithm(const edm::ParameterSet& pset, ApvLatencyAnalysis* const anal)
0015     : CommissioningAlgorithm(anal), histo_(nullptr, "") {
0016   ;
0017 }
0018
0019 // ----------------------------------------------------------------------------
0020 //
0021 void ApvLatencyAlgorithm::extract(const std::vector<TH1*>& histos) {
0022   if (!anal()) {
0023     edm::LogWarning(mlCommissioning_) << "[ApvLatencyAlgorithm::" << __func__ << "]"
0024                                       << " NULL pointer to Analysis object!";
0025     return;
0026   }
0027
0028   // Check
0029   if (histos.size() != 1) {
0030     anal()->addErrorCode(sistrip::numberOfHistos_);
0031   }
0032
0033   // Extract FED key from histo title
0034   if (!histos.empty()) {
0035     anal()->fedKey(extractFedKey(histos.front()));
0036   }
0037
0038   // Extract histograms
0039   std::vector<TH1*>::const_iterator ihis = histos.begin();
0040   for (; ihis != histos.end(); ihis++) {
0041     // Check for NULL pointer
0042     if (!(*ihis)) {
0043       continue;
0044     }
0045
0046     // Check name
0047     SiStripHistoTitle title((*ihis)->GetName());
0048     if (title.runType() != sistrip::APV_LATENCY) {
0049       anal()->addErrorCode(sistrip::unexpectedTask_);
0050       continue;
0051     }
0052
0053     // Extract timing histo
0054     histo_.first = *ihis;
0055     histo_.second = (*ihis)->GetName();
0056   }
0057 }
0058
0059 // ----------------------------------------------------------------------------
0060 //
0061 void ApvLatencyAlgorithm::analyse() {
0062   if (!anal()) {
0063     edm::LogWarning(mlCommissioning_) << "[ApvLatencyAlgorithm::" << __func__ << "]"
0064                                       << " NULL pointer to base Analysis object!";
0065     return;
0066   }
0067
0068   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0069   ApvLatencyAnalysis* anal = dynamic_cast<ApvLatencyAnalysis*>(tmp);
0070   if (!anal) {
0071     edm::LogWarning(mlCommissioning_) << "[ApvLatencyAlgorithm::" << __func__ << "]"
0072                                       << " NULL pointer to derived Analysis object!";
0073     return;
0074   }
0075
0076   // was in deprecated()
0077
0078   std::vector<const TProfile*> histos;
0079   std::vector<unsigned short> monitorables;
0080
0081   // was in analysis()
0082
0083   histos.clear();
0084   histos.push_back(const_cast<const TProfile*>(dynamic_cast<TProfile*>(histo_.first)));
0085   if (!histos[0]) {
0086     anal->addErrorCode(sistrip::nullPtr_);
0087     return;
0088   }
0089
0090   monitorables.clear();
0091
0092   //LogDebug("Commissioning|Algorithm") << "[ApvLatencyAlgorithm::analysis]";
0093
0094   //extract root histogram
0095   //check
0096   if (histos.size() != 1) {
0097     //     edm::LogWarning("Commissioning|Algorithm") << "[ApvLatencyAlgorithm::analysis]: Requires \"const std::vector<const TH1F*>& \" argument to have size 1. Actual size: " << histos.size() << ". Monitorables set to 0.";
0098     monitorables.push_back(0);
0099     return;
0100   }
0101   const TProfile* histo = histos[0];
0102
0103   //monitorable
0104   unsigned short latency;
0105
0106   std::vector<unsigned short> binContent;
0107   binContent.reserve((unsigned short)histo->GetNbinsX());
0108   binContent.resize((unsigned short)histo->GetNbinsX(), 0);
0109
0110   for (unsigned short k = 0; k < (unsigned short)histo->GetNbinsX(); k++) {  // k is bin number
0111
0112     //fill std::vector with histogram contents
0113     binContent.push_back((unsigned int)(histo->GetBinContent(k)));
0114   }
0115
0116   //calculate median
0117
0118   sort(binContent.begin(), binContent.end());
0119
0120   //calculate mean and mean2 of the readout within cutoffs
0121
0122   float meanNoise = 0.;  //M.W method
0123   float mean2Noise = 0.;
0124
0125   for (unsigned short k = (unsigned short)(binContent.size() * .1); k < (unsigned short)(binContent.size() * .9); k++) {
0126     meanNoise += binContent[k];
0127     mean2Noise += binContent[k] * binContent[k];
0128     ;
0129   }
0130
0131   meanNoise = meanNoise * binContent.size() * 0.8;
0132   mean2Noise = mean2Noise * binContent.size() * 0.8;
0133   float sigmaNoise = sqrt(fabs(meanNoise * meanNoise - mean2Noise));
0134
0135   //loop to look for signal > 5* sigma_noise
0136   unsigned short count = 0;
0137   unsigned short maxlatency = 0;
0138   unsigned int maxhits = 0;
0139
0140   for (unsigned short k = 1; k < ((unsigned short)histo->GetNbinsX() + 1); k++) {  // k is bin number
0141     if (histo->GetBinContent((Int_t)k) > maxhits)
0142       maxlatency = k - 1;
0143     if ((float)histo->GetBinContent((Int_t)k) > (meanNoise + 5 * sigmaNoise)) {
0144       latency = k - 1;
0145       count++;
0146     }
0147   }
0148
0149   if (!count) {
0150     //   LogDebug("Commissioning|Algorithm") << "[ApvLatencyAlgorithm::analysis]: Warning: no signal found > mean + 5*sigma(noise). Returning latency of highest number of recorded hits.";
0151     latency = maxlatency;
0152   }
0153
0154   if (count > 1) {
0155     //    LogDebug("Commissioning|Algorithm") << "[ApvLatencyAlgorithm::analysis]: Warning: more than one signal found > mean + 5*sigma(noise). Returning latency of highest number of recorded hits.";
0156     latency = maxlatency;
0157   }
0158
0159   //set monitorables
0160   monitorables.clear();
0161   monitorables.push_back(latency);
0162
0163   anal->latency_ = monitorables[0];
0164 }