JetMET/src/JetMETDQMPostProcessor.cc

0001 // Migrated to use DQMEDHarvester by: Jyothsna Rani Komaragiri, Oct 2014
0002
0003 #include "HLTriggerOffline/JetMET/interface/JetMETDQMPostProcessor.h"
0004
0005 #include "DQMServices/Core/interface/DQMStore.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008
0009 #include <cmath>
0010 #include <cstring>
0011 #include <fstream>
0012 #include <iomanip>
0013 #include <iostream>
0014
0015 JetMETDQMPostProcessor::JetMETDQMPostProcessor(const edm::ParameterSet &pset) {
0016   subDir_ = pset.getUntrackedParameter<std::string>("subDir");
0017   patternJetTrg_ = pset.getUntrackedParameter<std::string>("PatternJetTrg", "");
0018   patternMetTrg_ = pset.getUntrackedParameter<std::string>("PatternMetTrg", "");
0019 }
0020
0021 void JetMETDQMPostProcessor::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0022   //////////////////////////////////
0023   // setup DQM stor               //
0024   //////////////////////////////////
0025
0026   bool isJetDir = false;
0027   bool isMetDir = false;
0028
0029   TPRegexp patternJet(patternJetTrg_);
0030   TPRegexp patternMet(patternMetTrg_);
0031
0032   // go to the directory to be processed
0033   if (igetter.dirExists(subDir_))
0034     ibooker.cd(subDir_);
0035   else {
0036     edm::LogWarning("JetMETDQMPostProcessor") << "cannot find directory: " << subDir_ << " , skipping";
0037     return;
0038   }
0039
0040   std::vector<std::string> subdirectories = igetter.getSubdirs();
0041   for (std::vector<std::string>::iterator dir = subdirectories.begin(); dir != subdirectories.end(); dir++) {
0042     ibooker.cd(*dir);
0043
0044     isJetDir = false;
0045     isMetDir = false;
0046
0047     if (TString(*dir).Contains(patternJet))
0048       isJetDir = true;
0049     if (TString(*dir).Contains(patternMet))
0050       isMetDir = true;
0051
0052     if (isMetDir) {
0053       // std::cout << "JetMETDQMPostProcessor - Met paths: " << ibooker.pwd() <<
0054       // " " << *dir << std::endl;
0055
0056       // GenMET
0057       dividehistos(ibooker,
0058                    igetter,
0059                    "_meGenMETTrgMC",
0060                    "_meGenMET",
0061                    "_meTurnOngMET",
0062                    "Gen Missing ET",
0063                    "Gen Missing ET Turn-On RelVal");
0064       dividehistos(ibooker,
0065                    igetter,
0066                    "_meGenMETTrg",
0067                    "_meGenMETTrgLow",
0068                    "_meTurnOngMETLow",
0069                    "Gen Missing ETLow",
0070                    "Gen Missing ET Turn-On Data");
0071
0072       // HLTMET
0073       dividehistos(ibooker,
0074                    igetter,
0075                    "_meHLTMETTrgMC",
0076                    "_meHLTMET",
0077                    "_meTurnOnhMET",
0078                    "HLT Missing ET",
0079                    "HLT Missing ET Turn-On RelVal");
0080       dividehistos(ibooker,
0081                    igetter,
0082                    "_meHLTMETTrg",
0083                    "_meHLTMETTrgLow",
0084                    "_meTurnOnhMETLow",
0085                    "HLT Missing ETLow",
0086                    "HLT Missing ET Turn-On Data");
0087     }
0088
0089     if (isJetDir) {
0090       // std::cout << "JetMETDQMPostProcessor - Jet paths: " << ibooker.pwd() <<
0091       // " " << *dir << std::endl;
0092
0093       // GenJets
0094       dividehistos(ibooker,
0095                    igetter,
0096                    "_meGenJetPtTrgMC",
0097                    "_meGenJetPt",
0098                    "_meTurnOngJetPt",
0099                    "Gen Jet Pt",
0100                    "Gen Jet Pt Turn-On RelVal");
0101       dividehistos(ibooker,
0102                    igetter,
0103                    "_meGenJetPtTrg",
0104                    "_meGenJetPtTrgLow",
0105                    "_meTurnOngJetPtLow",
0106                    "Gen Jet PtLow",
0107                    "Gen Jet Pt Turn-On Data");
0108       dividehistos(ibooker,
0109                    igetter,
0110                    "_meGenJetEtaTrgMC",
0111                    "_meGenJetEta",
0112                    "_meTurnOngJetEta",
0113                    "Gen Jet Eta",
0114                    "Gen Jet Eta Turn-On RelVal");
0115       dividehistos(ibooker,
0116                    igetter,
0117                    "_meGenJetEtaTrg",
0118                    "_meGenJetEtaTrgLow",
0119                    "_meTurnOngJetEtaLow",
0120                    "Gen Jet EtaLow",
0121                    "Gen Jet Eta Turn-On Data");
0122       dividehistos(ibooker,
0123                    igetter,
0124                    "_meGenJetPhiTrgMC",
0125                    "_meGenJetPhi",
0126                    "_meTurnOngJetPhi",
0127                    "Gen Jet Phi",
0128                    "Gen Jet Phi Turn-On RelVal");
0129       dividehistos(ibooker,
0130                    igetter,
0131                    "_meGenJetPhiTrg",
0132                    "_meGenJetPhiTrgLow",
0133                    "_meTurnOngJetPhiLow",
0134                    "Gen Jet PhiLow",
0135                    "Gen Jet Phi Turn-On Data");
0136
0137       // HLTJets
0138       dividehistos(ibooker,
0139                    igetter,
0140                    "_meHLTJetPtTrgMC",
0141                    "_meHLTJetPt",
0142                    "_meTurnOnhJetPt",
0143                    "HLT Jet Pt",
0144                    "HLT Jet Pt Turn-On RelVal");
0145       dividehistos(ibooker,
0146                    igetter,
0147                    "_meHLTJetPtTrg",
0148                    "_meHLTJetPtTrgLow",
0149                    "_meTurnOnhJetPtLow",
0150                    "HLT Jet PtLow",
0151                    "HLT Jet Pt Turn-On Data");
0152       dividehistos(ibooker,
0153                    igetter,
0154                    "_meHLTJetEtaTrgMC",
0155                    "_meHLTJetEta",
0156                    "_meTurnOnhJetEta",
0157                    "HLT Jet Eta",
0158                    "HLT Jet Eta Turn-On RelVal");
0159       dividehistos(ibooker,
0160                    igetter,
0161                    "_meHLTJetEtaTrg",
0162                    "_meHLTJetEtaTrgLow",
0163                    "_meTurnOnhJetEtaLow",
0164                    "HLT Jet EtaLow",
0165                    "HLT Jet Eta Turn-On Data");
0166       dividehistos(ibooker,
0167                    igetter,
0168                    "_meHLTJetPhiTrgMC",
0169                    "_meHLTJetPhi",
0170                    "_meTurnOnhJetPhi",
0171                    "HLT Jet Phi",
0172                    "HLT Jet Phi Turn-On RelVal");
0173       dividehistos(ibooker,
0174                    igetter,
0175                    "_meHLTJetPhiTrg",
0176                    "_meHLTJetPhiTrgLow",
0177                    "_meTurnOnhJetPhiLow",
0178                    "HLT Jet PhiLow",
0179                    "HLT Jet Phi Turn-On Data");
0180     }
0181
0182     ibooker.goUp();
0183   }
0184 }
0185
0186 //----------------------------------------------------------------------
0187 TProfile *JetMETDQMPostProcessor::dividehistos(DQMStore::IBooker &ibooker,
0188                                                DQMStore::IGetter &igetter,
0189                                                const std::string &numName,
0190                                                const std::string &denomName,
0191                                                const std::string &outName,
0192                                                const std::string &label,
0193                                                const std::string &titel) {
0194   // ibooker.pwd();
0195   // std::cout << "In dividehistos: " << ibooker.pwd() << std::endl;
0196
0197   // std::cout << numName <<std::endl;
0198   TH1F *num = getHistogram(ibooker, igetter, ibooker.pwd() + "/" + numName);
0199
0200   // std::cout << denomName << std::endl;
0201   TH1F *denom = getHistogram(ibooker, igetter, ibooker.pwd() + "/" + denomName);
0202
0203   if (num == nullptr)
0204     edm::LogWarning("JetMETDQMPostProcessor")
0205         << "numerator histogram " << ibooker.pwd() + "/" + numName << " does not exist";
0206   if (denom == nullptr)
0207     edm::LogWarning("JetMETDQMPostProcessor")
0208         << "denominator histogram " << ibooker.pwd() + "/" + denomName << " does not exist";
0209
0210   // Check if histograms actually exist
0211   if (!num || !denom)
0212     return nullptr;
0213
0214   MonitorElement *meOut = ibooker.bookProfile(
0215       outName, titel, num->GetXaxis()->GetNbins(), num->GetXaxis()->GetXmin(), num->GetXaxis()->GetXmax(), 0., 1.2);
0216   meOut->setEfficiencyFlag();
0217   TProfile *out = meOut->getTProfile();
0218   out->GetXaxis()->SetTitle(label.c_str());
0219   out->SetYTitle("Efficiency");
0220   out->SetOption("PE");
0221   out->SetLineColor(2);
0222   out->SetLineWidth(2);
0223   out->SetMarkerStyle(20);
0224   out->SetMarkerSize(0.8);
0225   out->SetStats(kFALSE);
0226
0227   for (int i = 1; i <= num->GetNbinsX(); i++) {
0228     double e, low, high;
0229     Efficiency((int)num->GetBinContent(i), (int)denom->GetBinContent(i), 0.683, e, low, high);
0230     double err = e - low > high - e ? e - low : high - e;
0231     // here is the trick to store info in TProfile:
0232     out->SetBinContent(i, e);
0233     out->SetBinEntries(i, 1);
0234     out->SetBinError(i, sqrt(e * e + err * err));
0235   }
0236   return out;
0237 }
0238
0239 //----------------------------------------------------------------------
0240 TH1F *JetMETDQMPostProcessor::getHistogram(DQMStore::IBooker &ibooker,
0241                                            DQMStore::IGetter &igetter,
0242                                            const std::string &histoPath) {
0243   ibooker.pwd();
0244   MonitorElement *monElement = igetter.get(histoPath);
0245   if (monElement != nullptr)
0246     return monElement->getTH1F();
0247   else
0248     return nullptr;
0249 }
0250
0251 //----------------------------------------------------------------------
0252 void JetMETDQMPostProcessor::Efficiency(
0253     int passing, int total, double level, double &mode, double &lowerBound, double &upperBound) {
0254   // protection
0255   if (total == 0) {
0256     mode = 0.5;
0257     lowerBound = 0;
0258     upperBound = 1;
0259     return;
0260   }
0261   mode = passing / ((double)total);
0262
0263   // see http://root.cern.ch/root/html/TEfficiency.html#compare
0264   lowerBound = TEfficiency::Wilson(total, passing, level, false);
0265   upperBound = TEfficiency::Wilson(total, passing, level, true);
0266 }
0267
0268 //----------------------------------------------------------------------
0269
0270 DEFINE_FWK_MODULE(JetMETDQMPostProcessor);