Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:57:29 

0001 //
0002 // Revision 1.19  2010/07/20 02:58:27  wmtan
0003 // Add missing #include files
0004 //
0005 // Revision 1.18  2010/04/02 20:48:12  wittich
0006 // updates to scale entries by received number of FU's
0007 //
0008 // Revision 1.17  2010/03/17 20:56:19  wittich
0009 // Check for good updates based on mergeCount values
0010 // add code for rates normalized per FU
0011 //
0012 // Revision 1.16  2010/03/16 22:19:19  wittich
0013 // updates for per-LS normalization for variable
0014 // number of FU's sending information back to the clients.
0015 //
0016 // Revision 1.15  2010/02/11 23:55:18  wittich
0017 // - adapt to shorter Lumi Section length
0018 // - fix bug in how history of counts was filled
0019 //
0020 // Revision 1.14  2010/02/02 11:44:20  wittich
0021 // more diagnostics for online scalers
0022 //
0023 // Revision 1.13  2009/12/15 20:41:16  wittich
0024 // better hlt scalers client
0025 //
0026 // Revision 1.12  2009/11/22 13:33:05  puigh
0027 // clean beginJob
0028 //
0029 // Revision 1.11  2009/11/07 03:35:05  lorenzo
0030 // fixed binning
0031 //
0032 // Revision 1.10  2009/11/04 06:00:05  lorenzo
0033 // changed folders
0034 //
0035 
0036 ///// THIS ONLY WILL WORK WITH THE DATA FROM THE FU'S
0037 
0038 #include <cassert>
0039 #include <numeric>
0040 
0041 #include "DQM/TrigXMonitorClient/interface/HLTScalersClient.h"
0042 
0043 #include "FWCore/Framework/interface/Event.h"
0044 #include "FWCore/Framework/interface/LuminosityBlock.h"
0045 #include "FWCore/Framework/interface/Run.h"
0046 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0047 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0048 #include "FWCore/ServiceRegistry/interface/Service.h"
0049 #include "FWCore/Utilities/interface/isFinite.h"
0050 
0051 #include "DQMServices/Core/interface/DQMStore.h"
0052 
0053 using edm::LogInfo;
0054 using edm::LogWarning;
0055 
0056 // I am not sure this is right at more than 10%
0057 #define SECS_PER_LUMI_SECTION 23.31
0058 const int kPerHisto = 20;
0059 
0060 /// Constructors
0061 HLTScalersClient::HLTScalersClient(const edm::ParameterSet &ps)
0062     : dbe_(nullptr),
0063       nLumi_(0),
0064       currentRate_(nullptr),
0065       currentLumiBlockNumber_(0),
0066       first_(true),
0067       missingPathNames_(true),
0068       folderName_(ps.getUntrackedParameter<std::string>("dqmFolder", "HLT/HLScalers_EvF")),
0069       kRateIntegWindow_(ps.getUntrackedParameter<unsigned int>("rateIntegWindow", 3)),
0070       processName_(ps.getParameter<std::string>("processName")),
0071       ignores_(),
0072       debug_(ps.getUntrackedParameter<bool>("debugDump", false)),
0073       maxFU_(ps.getUntrackedParameter<unsigned int>("maxFU", false)),
0074       recentOverallCountsPerLS_(kRateIntegWindow_),
0075       recentNormedOverallCountsPerLS_(2) {
0076   LogDebug("HLTScalersClient") << "constructor";
0077   if (debug_) {
0078     textfile_.open("debug.txt");
0079     if (!textfile_) {
0080       std::cout << "constructor: can't open text file" << std::endl;
0081     }
0082   }
0083   // get back-end interface
0084   usesResource("DQMStore");
0085   dbe_ = edm::Service<DQMStore>().operator->();
0086   dbe_->setCurrentFolder(folderName_);
0087 
0088   std::string rawdir(folderName_ + "/raw");
0089   dbe_->setCurrentFolder(rawdir);
0090   hltRate_ = dbe_->book1D("hltRate", "Overall HLT Accept rate vs LS", MAX_LUMI_SEG_HLT, -0.5, MAX_LUMI_SEG_HLT - 0.5);
0091   dbe_->setCurrentFolder(folderName_);
0092   hltNormRate_ = dbe_->book1D(
0093       "hltRateNorm", "Overall HLT Accept rate vs LS, scaled", MAX_LUMI_SEG_HLT, -0.5, MAX_LUMI_SEG_HLT - 0.5);
0094   hltCount_ = dbe_->book1D("hltCount", "Overall HLT Counts vs LS", MAX_LUMI_SEG_HLT, -0.5, MAX_LUMI_SEG_HLT - 0.5);
0095   //   hltCountN_= dbe_->book1D("hltCountN", "Overall HLT Counts per LS vs LS",
0096   //              MAX_LUMI_SEG_HLT, -0.5, MAX_LUMI_SEG_HLT-0.5);
0097   mergeCount_ =
0098       dbe_->book1D("mergeCount", "Number of merge counts vs LS", MAX_LUMI_SEG_HLT, -0.5, MAX_LUMI_SEG_HLT - 0.5);
0099 
0100   updates_ = dbe_->book1D("updates", "Status of Updates", 2, 0, 2);
0101   updates_->setBinLabel(1, "Good Updates");
0102   updates_->setBinLabel(2, "Incomplete Updates");
0103 
0104 }  // end constructor
0105 
0106 /// BeginJob
0107 void HLTScalersClient::beginJob(void) {
0108   LogDebug("HLTScalersClient") << "beingJob";
0109   if (dbe_) {
0110     dbe_->setCurrentFolder(folderName_);
0111   }
0112   first_ = true;
0113   missingPathNames_ = true;
0114 }
0115 
0116 /// beginRun
0117 void HLTScalersClient::beginRun(const edm::Run &run, const edm::EventSetup &c) {
0118   missingPathNames_ = true;
0119   first_ = true;
0120   LogDebug("HLTScalersClient") << "beginRun, run " << run.id();
0121 
0122 }  // beginRun
0123 
0124 /// EndRun
0125 void HLTScalersClient::endRun(const edm::Run &run, const edm::EventSetup &c) {
0126   missingPathNames_ = true;
0127   first_ = true;
0128 }
0129 
0130 /// End LumiBlock
0131 /// DQM Client Diagnostic should be performed here
0132 void HLTScalersClient::endLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &c) {
0133   nLumi_ = lumiSeg.id().luminosityBlock();
0134 
0135   // get raw data
0136   std::string scalHisto = folderName_ + "/raw/hltScalers";
0137   MonitorElement *scalers = dbe_->get(scalHisto);
0138   if (scalers == nullptr) {
0139     LogDebug("HLTScalersClient") << "cannot get hlt scalers histogram, "
0140                                  << "bailing out.";
0141     if (debug_)
0142       std::cout << "No scalers ? Looking for " << scalHisto << std::endl;
0143     return;
0144   }
0145 
0146   int npaths = scalers->getNbinsX();
0147   if (npaths > MAX_PATHS)
0148     npaths = MAX_PATHS;  // HARD CODE FOR NOW
0149   LogDebug("HLTScalersClient") << "I see " << npaths << " paths. ";
0150 
0151   // set the bin labels on the first go-through
0152   // I need to do this here because we don't have the paths yet
0153   // on begin-run. I should do this in a less ugly way (see FV?)
0154   if (first_) {
0155     std::string rawdir(folderName_ + "/raw");
0156 
0157     LogDebug("HLTScalersClient") << "Setting up paths on first endLumiBlock " << npaths;
0158     dbe_->setCurrentFolder(rawdir);
0159     currentRate_ = dbe_->book1D("cur_rate", "current lumi section rate per path", npaths, -0.5, npaths - 0.5);
0160     currentNormRate_ =
0161         dbe_->book1D("cur_rate_norm", "current norm. lumi section rate per path", npaths, -0.5, npaths - 0.5);
0162     recentPathCountsPerLS_.reserve(npaths);
0163     recentNormedPathCountsPerLS_.reserve(npaths);
0164     char rates_subfolder[256];
0165     snprintf(rates_subfolder, 256, "%s/RateHistory", folderName_.c_str());
0166     char counts_subfolder[256];
0167     snprintf(counts_subfolder, 256, "%s/CountHistory", folderName_.c_str());
0168 
0169     hltCurrentRate_.reserve(npaths);
0170     rateHistories_.reserve(npaths);
0171     countHistories_.reserve(npaths);
0172     hltCurrentNormRate_.reserve(npaths);
0173     rateNormHistories_.reserve(npaths);
0174 
0175     dbe_->setCurrentFolder(folderName_);  // these belong in top-level
0176     for (int i = 0; i < npaths; ++i) {
0177       dbe_->setCurrentFolder(std::string(rates_subfolder) + "/raw");
0178 
0179       char name[256];
0180       snprintf(name, 256, "raw_rate_p%03d", i);
0181       //     LogDebug("HLTScalersClient") << "name " << i << " is " << name ;
0182       rateHistories_.push_back(dbe_->book1D(name, name, MAX_LUMI_SEG_HLT, -0.5, MAX_LUMI_SEG_HLT - 0.5));
0183       snprintf(name, 256, "norm_rate_p%03d", i);
0184       dbe_->setCurrentFolder(rates_subfolder);
0185       rateNormHistories_.push_back(dbe_->book1D(name, name, MAX_LUMI_SEG_HLT, -0.5, MAX_LUMI_SEG_HLT - 0.5));
0186       dbe_->setCurrentFolder(counts_subfolder);
0187       snprintf(name, 256, "counts_p%03d", i);
0188       countHistories_.push_back(dbe_->book1D(name, name, MAX_LUMI_SEG_HLT, -0.5, MAX_LUMI_SEG_HLT - 0.5));
0189       // prefill the data structures
0190       recentPathCountsPerLS_.push_back(CountLSFifo_t(kRateIntegWindow_));
0191       recentNormedPathCountsPerLS_.push_back(CountLSFifo_t(2));
0192     }
0193     dbe_->setCurrentFolder(folderName_);
0194 
0195     // split hlt scalers up into groups of 20
0196     const int maxlen = 64;
0197     char metitle[maxlen];                     // histo name
0198     char mename[maxlen];                      // ME name
0199     int numHistos = int(npaths / kPerHisto);  // this hasta be w/o remainders
0200 
0201     int remainder = npaths % kPerHisto;
0202     if (remainder)
0203       numHistos += 1;
0204 
0205     for (int k = 0; k < numHistos; k++) {
0206       int npath_low = kPerHisto * k;
0207       int npath_high = kPerHisto * (k + 1) - 1;
0208       snprintf(mename, maxlen, "hltScalers_%0d", k);
0209       snprintf(metitle, maxlen, "HLT scalers - Paths %d to %d", npath_low, npath_high);
0210       dbe_->setCurrentFolder(rawdir);
0211       hltCurrentRate_.push_back(dbe_->book1D(mename, metitle, kPerHisto, -0.5 + npath_low, npath_high + 0.5));
0212       dbe_->setCurrentFolder(folderName_);  // these belong in top-level
0213       snprintf(mename, maxlen, "hltScalersNorm_%0d", k);
0214       snprintf(metitle, maxlen, "HLT Rate (scaled) - Paths %d to %d", npath_low, npath_high);
0215       hltCurrentNormRate_.push_back(dbe_->book1D(mename, metitle, kPerHisto, -0.5 + npath_low, npath_high + 0.5));
0216     }
0217 
0218     first_ = false;
0219   }
0220 
0221   if (missingPathNames_) {
0222     // first try the scalers histogram and see if the bin names are set
0223     for (int i = 0; i < npaths; ++i) {
0224       // needs to be i+1 since root bins start at 1
0225       const char *name = scalers->getTH1()->GetXaxis()->GetBinLabel(i + 1);
0226       if (name && (strlen(name) > 0)) {
0227         if (debug_) {
0228           std::cout << "path " << i << " name is " << name << std::endl;
0229         }
0230         int whichHisto = i / kPerHisto;
0231         int whichBin = i % kPerHisto + 1;
0232         char pname[256];
0233         hltCurrentRate_[whichHisto]->setBinLabel(whichBin, name);
0234         hltCurrentNormRate_[whichHisto]->setBinLabel(whichBin, name);
0235         snprintf(pname, 256, "Rate - path %s (Path # %03d)", name, i);
0236         rateHistories_[i]->setTitle(pname);
0237         rateNormHistories_[i]->setTitle(pname);
0238         snprintf(pname, 256, "Counts - path %s (Path # %03d)", name, i);
0239         countHistories_[i]->setTitle(pname);
0240 
0241         currentRate_->setBinLabel(i + 1, name);
0242         currentNormRate_->setBinLabel(i + 1, name);
0243 
0244         missingPathNames_ = false;
0245       }
0246     }
0247   }
0248   // MEGA-HACK
0249   if (missingPathNames_) {
0250     // if that didn't work we load 'em from a text file. damn straight.
0251     int ipath = 1;
0252     std::ifstream names("names.dat");
0253     if (!names) {
0254       if (debug_) {
0255         std::ostringstream msg;
0256         msg << "open of "
0257             << "names.dat";
0258         perror(msg.str().c_str());
0259       }
0260     } else {  // open succeeded
0261       missingPathNames_ = false;
0262       std::string line;
0263       while (!names.eof()) {
0264         getline(names, line);
0265         std::istringstream fnames(line);
0266         std::string label;
0267         int bin;
0268         if (fnames.str().find('#') == 0)  // skip comment lines
0269           continue;
0270         if (fnames >> bin >> label) {
0271           if (debug_) {
0272             std::cout << bin << "--" << label << "(" << ipath << ")" << std::endl;
0273           }
0274           currentRate_->setBinLabel(ipath, label);
0275           currentNormRate_->setBinLabel(ipath, label);
0276           countHistories_[ipath - 1]->setTitle(label);
0277           rateHistories_[ipath - 1]->setTitle(label);
0278           rateNormHistories_[ipath - 1]->setTitle(label);
0279           int whichHisto = (ipath - 1) / kPerHisto;
0280           int whichBin = (ipath - 1) % kPerHisto + 1;
0281           hltCurrentRate_[whichHisto]->setBinLabel(whichBin, label);
0282           hltCurrentNormRate_[whichHisto]->setBinLabel(whichBin, label);
0283           ++ipath;
0284           if (ipath > npaths)
0285             break;
0286         }  //
0287       }    // loop lines
0288     }      // open ok
0289   }
0290 
0291   // END SETUP
0292 
0293   std::string nLumiHisto(folderName_ + "/nLumiBlock");
0294   MonitorElement *nLumi = dbe_->get(nLumiHisto);
0295   if (nLumi == nullptr) {
0296     nLumiHisto = folderName_ + "/raw/nLumiBlock";
0297     nLumi = dbe_->get(nLumiHisto);
0298   }
0299   int testval = (nLumi != nullptr ? nLumi->getIntValue() : -1);
0300   LogDebug("HLTScalersClient") << "Lumi Block from DQM: " << testval << ", local is " << nLumi_;
0301   int nL = (nLumi != nullptr ? nLumi->getIntValue() : nLumi_);
0302   if (nL > MAX_LUMI_SEG_HLT) {
0303     LogDebug("HLTScalersClient") << "Too many Lumi segments, " << nL << " is greater than MAX_LUMI_SEG_HLT,"
0304                                  << " wrapping to " << (nL % MAX_LUMI_SEG_HLT);
0305     // nL = MAX_LUMI_SEG_HLT;
0306     nL = nL % MAX_LUMI_SEG_HLT;
0307   }
0308 
0309   // merging counts
0310   double num_fu = -1.0;
0311   std::string mergeName(folderName_ + "/raw/hltMerge");
0312   MonitorElement *merge = dbe_->get(mergeName);
0313   if (merge != nullptr) {
0314     num_fu = merge->getBinContent(1);
0315     if (debug_) {
0316       std::cout << "Number of received entries: " << num_fu << std::endl;
0317     }
0318     mergeCount_->Fill(nL, num_fu);
0319   }
0320   // end
0321 
0322   // evaluate the data
0323   // loop over current counts
0324 
0325   // this gets filled for all times. Mainly for debugging.
0326   for (int i = 1; i <= npaths; ++i) {  // bins start at 1
0327     double current_count = scalers->getBinContent(i);
0328     // nb that this will now _overwrite_ some settings
0329     countHistories_[i - 1]->setBinContent(nL, current_count);  // good or bad
0330   }
0331 
0332   std::string overallScalerName(folderName_ + "/raw/hltOverallScaler");
0333   MonitorElement *hltScaler = dbe_->get(overallScalerName);
0334   if (hltScaler != nullptr) {
0335     double current_count = hltScaler->getBinContent(1);
0336     hltCount_->setBinContent(nL, current_count);
0337     recentOverallCountsPerLS_.update(CountLS_t(nL, current_count));
0338     std::pair<double, double> sl = getSlope_(recentOverallCountsPerLS_);
0339     double slope = sl.first;
0340     double slope_err = sl.second;
0341     if (slope > 0) {
0342       hltRate_->setBinContent(nL, slope);
0343       if (!edm::isNotFinite(slope_err) && (slope_err >= 0))
0344         hltRate_->setBinError(nL, slope_err);
0345     }
0346   }  // found  histo
0347 
0348   if (num_fu >= 0.95 * maxFU_) {
0349     if (num_fu > maxFU_) {
0350       maxFU_ = num_fu;
0351       if (debug_)
0352         std::cout << "maxFU is now " << maxFU_ << std::endl;
0353     }
0354     // good data
0355     for (int i = 1; i <= npaths; ++i) {  // bins start at 1
0356       double current_count = scalers->getBinContent(i);
0357       // DEBUG
0358       if (!recentPathCountsPerLS_[i - 1].empty() && debug_)
0359         std::cout << i << "\t-> good one: new => cnt, ls = " << current_count << ", " << nL
0360                   << ", old = " << recentPathCountsPerLS_[i - 1].back().second << "\t"
0361                   << recentPathCountsPerLS_[i - 1].back().first << std::endl;
0362       // END DEBUG
0363       recentPathCountsPerLS_[i - 1].update(CountLS_t(nL, current_count));
0364 
0365       // NB: we do not fill a new entry in the rate histo if we can't
0366       // calculate it
0367       std::pair<double, double> sl = getSlope_(recentPathCountsPerLS_[i - 1]);
0368       double slope = sl.first;
0369       double slope_err = sl.second;
0370       // rateHistories_[i-1]->Fill(nL,slope);
0371       if (slope > 0) {
0372         rateHistories_[i - 1]->setBinContent(nL, slope);
0373         // set the current rate(s)
0374         hltCurrentRate_[(i - 1) / kPerHisto]->setBinContent(i % kPerHisto, slope);
0375         currentRate_->setBinContent(i, slope);
0376         if (!edm::isNotFinite(slope_err) && (slope_err >= 0)) {
0377           currentRate_->setBinError(i, slope_err);
0378           hltCurrentRate_[(i - 1) / kPerHisto]->setBinError(i % kPerHisto, slope_err);
0379           rateHistories_[i - 1]->setBinError(nL, slope_err);
0380         }
0381       }
0382       //// HACK - normalized path rates
0383       //// END HACK
0384 
0385     }  // loop over paths
0386 
0387     // ---------------------------- overall rate, absolute counts
0388     std::string overallScalerName(folderName_ + "/raw/hltOverallScaler");
0389     MonitorElement *hltScaler = dbe_->get(overallScalerName);
0390     if (hltScaler != nullptr) {
0391       double current_count = hltScaler->getBinContent(1);
0392       hltCount_->setBinContent(nL, current_count);
0393       recentOverallCountsPerLS_.update(CountLS_t(nL, current_count));
0394       std::pair<double, double> sl = getSlope_(recentOverallCountsPerLS_);
0395       double slope = sl.first;
0396       double slope_err = sl.second;
0397       if (slope >= 0) {
0398         hltRate_->setBinContent(nL, slope);
0399         if (!edm::isNotFinite(slope_err) && (slope_err >= 0))
0400           hltRate_->setBinError(nL, slope_err);
0401       }
0402     }                   // found  histo
0403     updates_->Fill(0);  // good
0404   }                     // check on number of FU's - good data
0405   else {
0406     updates_->Fill(1);  // missing updates
0407   }
0408 
0409   // PW DEBUG
0410   if (debug_) {
0411     textfile_ << nL << "\t" << npaths << "\t";
0412     for (int i = 0; i < npaths; ++i) {
0413       textfile_ << scalers->getBinContent(i) << " ";
0414     }
0415     textfile_ << std::endl;
0416   }
0417   // end DEBUG
0418 
0419 #ifdef LATER
0420   // ------ overall rate normalized - all data
0421   overallScalerName = std::string(folderName_ + "/raw/hltOverallScalerN");
0422   hltScaler = dbe_->get(overallScalerName);
0423   if (hltScaler != 0) {
0424     double cnt = hltScaler->getBinContent(1);
0425     //     hltCountN_->setBinContent(nL,cnt);
0426     if (debug_) {
0427       std::cout << "Overall Norm: new => cnt, ls = " << cnt << ", " << nL << ", num_fu = " << num_fu << std::endl;
0428     }
0429     recentNormedOverallCountsPerLS_.update(CountLS_t(nL, cnt / num_fu));
0430     cnt = recentNormedOverallCountsPerLS_.getCount(nL);  // for dupes/partials
0431     double slope = cnt / num_fu / SECS_PER_LUMI_SECTION;
0432     if (debug_) {
0433       std::cout << "Normalized slope = " << slope << std::endl;
0434     }
0435     if (slope > 0)
0436       hltNormRate_->setBinContent(nL, slope);
0437   }
0438   //
0439   std::string scalHistoNorm = folderName_ + "/raw/hltScalersN";
0440   MonitorElement *scalersN = dbe_->get(scalHistoNorm);
0441   if (scalersN) {
0442     for (int i = 0; i < npaths; ++i) {
0443       double cnt = scalersN->getBinContent(i);
0444       double slope = cnt / num_fu / SECS_PER_LUMI_SECTION;
0445       if (slope > 0) {
0446         rateNormHistories_[i - 1]->setBinContent(nL, slope);
0447         // set the current rate(s)
0448         hltCurrentNormRate_[(i - 1) / kPerHisto]->setBinContent(i % kPerHisto, slope);
0449         currentNormRate_->setBinContent(i, slope);
0450       }
0451     }
0452   }
0453 #else  // NOT LATER
0454   // ------ overall rate normalized - all data
0455   overallScalerName = std::string(folderName_ + "/raw/hltOverallScaler");
0456   hltScaler = dbe_->get(overallScalerName);
0457   if (hltScaler != nullptr) {
0458     double cnt = hltScaler->getBinContent(1);
0459     //     hltCountN_->setBinContent(nL,cnt);
0460     float sf = num_fu / maxFU_;
0461     if (debug_) {
0462       std::cout << "Overall Norm: new => cnt, ls = " << cnt << ", " << nL << ", num_fu = " << num_fu << ", sf = " << sf
0463                 << std::endl;
0464     }
0465     recentNormedOverallCountsPerLS_.update(CountLS_t(nL, cnt / sf));
0466     cnt = recentNormedOverallCountsPerLS_.getCount(nL);  // for dupes/partials
0467     std::pair<double, double> sl = getSlope_(recentNormedOverallCountsPerLS_);
0468     double slope = sl.first;
0469     double slope_err = sl.second;
0470     if (debug_) {
0471       std::cout << "Normalized slope = " << slope << std::endl;
0472     }
0473     if (slope > 0) {
0474       hltNormRate_->setBinContent(nL, slope);
0475       if (cnt > 0)
0476         slope_err = slope * sqrt(2. / num_fu + 2. / cnt);
0477       if (!edm::isNotFinite(slope_err) && (slope_err >= 0))
0478         hltNormRate_->setBinError(nL, slope_err);
0479     }
0480   }
0481   //
0482   std::string scalHistoNorm = folderName_ + "/raw/hltScalers";
0483   MonitorElement *scalersN = dbe_->get(scalHistoNorm);
0484   if (scalersN) {
0485     double sf = num_fu / maxFU_;
0486     for (int i = 1; i <= npaths; ++i) {
0487       double cnt = scalersN->getBinContent(i);
0488       recentNormedPathCountsPerLS_[i - 1].update(CountLS_t(nL, cnt / sf));
0489       std::pair<double, double> sl = getSlope_(recentNormedPathCountsPerLS_[i - 1]);
0490       double slope = sl.first;
0491       double slope_err = sl.second;
0492       if (slope >= 0) {
0493         rateNormHistories_[i - 1]->setBinContent(nL, slope);
0494         // set the current rate(s)
0495         hltCurrentNormRate_[(i - 1) / kPerHisto]->setBinContent(i % kPerHisto, slope);
0496         currentNormRate_->setBinContent(i, slope);
0497         if (slope_err <= 0 && cnt > 0) {
0498           // ignores error on prev point, so scale by sqrt(2)
0499           slope_err = slope * sqrt(2. / num_fu + 2. / cnt);
0500           if (debug_) {
0501             std::cout << "Slope err " << i << " = " << slope_err << std::endl;
0502           }
0503         }
0504         if (!edm::isNotFinite(slope_err) && (slope_err >= 0)) {
0505           rateNormHistories_[i - 1]->setBinError(nL, slope_err);
0506           // set the current rate(s)
0507           hltCurrentNormRate_[(i - 1) / kPerHisto]->setBinError(i % kPerHisto, slope_err);
0508           currentNormRate_->setBinError(i, slope_err);
0509         }
0510       }
0511     }
0512   }
0513 
0514 #endif  // LATER
0515   //
0516 }
0517 
0518 // unused
0519 void HLTScalersClient::analyze(const edm::Event &e, const edm::EventSetup &c) {
0520   // nothing to do here
0521 }
0522 
0523 // this is probably overkill. Do a least-squares fit to get slope
0524 // note that the data is in units of counts, ls number
0525 // but we return a value in Hz...
0526 std::pair<double, double> HLTScalersClient::getSlope_(const HLTScalersClient::CountLSFifo_t &points) {
0527   double slope, sigma_m;
0528   if (points.size() < points.targetSize()) {
0529     return std::pair<double, double>(-1, -1);
0530   }
0531   // just do a delta if we just want two bins
0532   else if (points.size() == 2) {
0533     // just do diff and be done with it
0534     double delta_ls = points.front().first - points.back().first;
0535     double delta_cnt = points.front().second - points.back().second;
0536     slope = delta_cnt / delta_ls;
0537     sigma_m = -1;
0538   } else {  // do a fit
0539     double xy = 0;
0540     double x = 0;
0541     double xsq = 0;
0542     double y = 0;
0543     double n = double(points.size());
0544     for (auto i(points.begin()); i != points.end(); ++i) {
0545       if (debug_)
0546         std::cout << "x = " << i->first << ", y = " << i->second << std::endl;
0547       xy += i->first * i->second;
0548       x += i->first;
0549       xsq += i->first * i->first;
0550       y += i->second;
0551     }
0552     slope = (n * xy - x * y) / (n * xsq - x * x);
0553 
0554     // now get the uncertainty on the slope. Need intercept for this.
0555     double intercept = (xsq * y - xy * x) / (n * xsq - x * x);
0556     double sigma_ysq = 0;
0557     for (auto i(points.begin()); i != points.end(); ++i) {
0558       sigma_ysq += pow((i->second - slope * i->first - intercept), 2.);
0559     }
0560     //     if ( debug_ )
0561     //       std::cout << "chi^2 = " << sigma_ysq << std::endl;
0562     sigma_ysq *= 1. / (n - 2.);
0563 
0564     sigma_m = sqrt(n * sigma_ysq / (n * xsq - x * x));
0565   }
0566 
0567   //  std::cout << "Slope is " << slope << std::endl;
0568   slope /= SECS_PER_LUMI_SECTION;
0569   if (sigma_m > 0)
0570     sigma_m /= SECS_PER_LUMI_SECTION;
0571   //   if ( debug_ ) {
0572   //     std::cout << "Slope = " << slope << " +- " << sigma_m
0573   //          << std::endl;
0574   //   std::cout << "intercept is " << intercept
0575   //        << std::endl;
0576   //  }
0577 
0578   return std::pair<double, double>(slope, sigma_m);
0579 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!