Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-29 23:10:47

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0009 #include "DQMServices/Core/interface/DQMStore.h"
0010 
0011 #include "TH1F.h"
0012 #include "TH2F.h"
0013 #include "TH3F.h"
0014 #include "RVersion.h"
0015 
0016 #include "TEfficiency.h"
0017 
0018 using namespace edm;
0019 using namespace std;
0020 
0021 class HeavyFlavorHarvesting : public DQMEDHarvester {
0022 public:
0023   HeavyFlavorHarvesting(const edm::ParameterSet& pset);
0024   ~HeavyFlavorHarvesting() override;
0025   // virtual void endRun(const edm::Run &, const edm::EventSetup &) override;
0026 private:
0027   void calculateEfficiency(const ParameterSet& pset, DQMStore::IBooker&, DQMStore::IGetter&);
0028   void calculateEfficiency1D(TH1* num, TH1* den, string name, DQMStore::IBooker&, DQMStore::IGetter&);
0029   void calculateEfficiency2D(TH2F* num, TH2F* den, string name, DQMStore::IBooker&, DQMStore::IGetter&);
0030 
0031   string myDQMrootFolder;
0032   const VParameterSet efficiencies;
0033 
0034 protected:
0035   void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;  //performed in the endJob
0036 };
0037 
0038 HeavyFlavorHarvesting::HeavyFlavorHarvesting(const edm::ParameterSet& pset)
0039     : myDQMrootFolder(pset.getUntrackedParameter<string>("MyDQMrootFolder")),
0040       efficiencies(pset.getUntrackedParameter<VParameterSet>("Efficiencies")) {}
0041 
0042 void HeavyFlavorHarvesting::dqmEndJob(DQMStore::IBooker& ibooker_, DQMStore::IGetter& igetter_) {
0043   for (VParameterSet::const_iterator pset = efficiencies.begin(); pset != efficiencies.end(); pset++) {
0044     calculateEfficiency(*pset, ibooker_, igetter_);
0045   }
0046 }
0047 
0048 void HeavyFlavorHarvesting::calculateEfficiency(const ParameterSet& pset,
0049                                                 DQMStore::IBooker& ibooker_,
0050                                                 DQMStore::IGetter& igetter_) {
0051   //get hold of numerator and denominator histograms
0052   vector<string> numDenEffMEnames = pset.getUntrackedParameter<vector<string> >("NumDenEffMEnames");
0053   if (numDenEffMEnames.size() != 3) {
0054     LogDebug("HLTriggerOfflineHeavyFlavor") << "NumDenEffMEnames must have three names" << endl;
0055     return;
0056   }
0057   string denMEname = myDQMrootFolder + "/" + numDenEffMEnames[1];
0058   string numMEname = myDQMrootFolder + "/" + numDenEffMEnames[0];
0059   MonitorElement* denME = igetter_.get(denMEname);
0060   MonitorElement* numME = igetter_.get(numMEname);
0061   if (denME == nullptr || numME == nullptr) {
0062     LogDebug("HLTriggerOfflineHeavyFlavor") << "Could not find MEs: " << denMEname << " or " << numMEname << endl;
0063     return;
0064   }
0065   TH1* den = denME->getTH1();
0066   TH1* num = numME->getTH1();
0067   //check compatibility of the histograms
0068   if (den->GetNbinsX() != num->GetNbinsX() || den->GetNbinsY() != num->GetNbinsY() ||
0069       den->GetNbinsZ() != num->GetNbinsZ()) {
0070     LogDebug("HLTriggerOfflineHeavyFlavor")
0071         << "Monitoring elements " << numMEname << " and " << denMEname << "are incompatible" << endl;
0072     return;
0073   }
0074   //figure out the directory and efficiency name
0075   string effName = numDenEffMEnames[2];
0076   string effDir = myDQMrootFolder;
0077   string::size_type slashPos = effName.rfind('/');
0078   if (string::npos != slashPos) {
0079     effDir += "/" + effName.substr(0, slashPos);
0080     effName.erase(0, slashPos + 1);
0081   }
0082   ibooker_.setCurrentFolder(effDir);
0083   //calculate the efficiencies
0084   int dimensions = num->GetDimension();
0085   if (dimensions == 1) {
0086     calculateEfficiency1D(num, den, effName, ibooker_, igetter_);
0087   } else if (dimensions == 2) {
0088     calculateEfficiency2D((TH2F*)num, (TH2F*)den, effName, ibooker_, igetter_);
0089     TH1D* numX = ((TH2F*)num)->ProjectionX();
0090     TH1D* denX = ((TH2F*)den)->ProjectionX();
0091     calculateEfficiency1D(numX, denX, effName + "X", ibooker_, igetter_);
0092     delete numX;
0093     delete denX;
0094     TH1D* numY = ((TH2F*)num)->ProjectionY();
0095     TH1D* denY = ((TH2F*)den)->ProjectionY();
0096     calculateEfficiency1D(numY, denY, effName + "Y", ibooker_, igetter_);
0097     delete numY;
0098     delete denY;
0099   } else {
0100     return;
0101   }
0102 }
0103 
0104 void HeavyFlavorHarvesting::calculateEfficiency1D(
0105     TH1* num, TH1* den, string effName, DQMStore::IBooker& ibooker_, DQMStore::IGetter& igetter_) {
0106   TProfile* eff;
0107   if (num->GetXaxis()->GetXbins()->GetSize() == 0) {
0108     eff = new TProfile(effName.c_str(),
0109                        effName.c_str(),
0110                        num->GetXaxis()->GetNbins(),
0111                        num->GetXaxis()->GetXmin(),
0112                        num->GetXaxis()->GetXmax());
0113   } else {
0114     eff = new TProfile(
0115         effName.c_str(), effName.c_str(), num->GetXaxis()->GetNbins(), num->GetXaxis()->GetXbins()->GetArray());
0116   }
0117   eff->SetTitle(effName.c_str());
0118   eff->SetXTitle(num->GetXaxis()->GetTitle());
0119   eff->SetYTitle("Efficiency");
0120   eff->SetOption("PE");
0121   eff->SetLineColor(2);
0122   eff->SetLineWidth(2);
0123   eff->SetMarkerStyle(20);
0124   eff->SetMarkerSize(0.8);
0125   eff->GetYaxis()->SetRangeUser(-0.001, 1.001);
0126   eff->SetStats(kFALSE);
0127   for (int i = 1; i <= num->GetNbinsX(); i++) {
0128     double e, low, high;
0129     if (int(den->GetBinContent(i)) > 0.)
0130       e = double(num->GetBinContent(i)) / double(den->GetBinContent(i));
0131     else
0132       e = 0.;
0133     low = TEfficiency::Wilson((double)den->GetBinContent(i), (double)num->GetBinContent(i), 0.683, false);
0134     high = TEfficiency::Wilson((double)den->GetBinContent(i), (double)num->GetBinContent(i), 0.683, true);
0135 
0136     double err = e - low > high - e ? e - low : high - e;
0137     //here is the trick to store info in TProfile:
0138     eff->SetBinContent(i, e);
0139     eff->SetBinEntries(i, 1);
0140     eff->SetBinError(i, sqrt(e * e + err * err));
0141   }
0142   ibooker_.bookProfile(effName, eff);
0143   delete eff;
0144 }
0145 
0146 void HeavyFlavorHarvesting::calculateEfficiency2D(
0147     TH2F* num, TH2F* den, string effName, DQMStore::IBooker& ibooker_, DQMStore::IGetter& igetter_) {
0148   TProfile2D* eff;
0149   if (num->GetXaxis()->GetXbins()->GetSize() == 0 && num->GetYaxis()->GetXbins()->GetSize() == 0) {
0150     eff = new TProfile2D(effName.c_str(),
0151                          effName.c_str(),
0152                          num->GetXaxis()->GetNbins(),
0153                          num->GetXaxis()->GetXmin(),
0154                          num->GetXaxis()->GetXmax(),
0155                          num->GetYaxis()->GetNbins(),
0156                          num->GetYaxis()->GetXmin(),
0157                          num->GetYaxis()->GetXmax());
0158   } else if (num->GetXaxis()->GetXbins()->GetSize() != 0 && num->GetYaxis()->GetXbins()->GetSize() == 0) {
0159     eff = new TProfile2D(effName.c_str(),
0160                          effName.c_str(),
0161                          num->GetXaxis()->GetNbins(),
0162                          num->GetXaxis()->GetXbins()->GetArray(),
0163                          num->GetYaxis()->GetNbins(),
0164                          num->GetYaxis()->GetXmin(),
0165                          num->GetYaxis()->GetXmax());
0166   } else if (num->GetXaxis()->GetXbins()->GetSize() == 0 && num->GetYaxis()->GetXbins()->GetSize() != 0) {
0167     eff = new TProfile2D(effName.c_str(),
0168                          effName.c_str(),
0169                          num->GetXaxis()->GetNbins(),
0170                          num->GetXaxis()->GetXmin(),
0171                          num->GetXaxis()->GetXmax(),
0172                          num->GetYaxis()->GetNbins(),
0173                          num->GetYaxis()->GetXbins()->GetArray());
0174   } else {
0175     eff = new TProfile2D(effName.c_str(),
0176                          effName.c_str(),
0177                          num->GetXaxis()->GetNbins(),
0178                          num->GetXaxis()->GetXbins()->GetArray(),
0179                          num->GetYaxis()->GetNbins(),
0180                          num->GetYaxis()->GetXbins()->GetArray());
0181   }
0182   eff->SetTitle(effName.c_str());
0183   eff->SetXTitle(num->GetXaxis()->GetTitle());
0184   eff->SetYTitle(num->GetYaxis()->GetTitle());
0185   eff->SetZTitle("Efficiency");
0186   eff->SetOption("colztexte");
0187   eff->GetZaxis()->SetRangeUser(-0.001, 1.001);
0188   eff->SetStats(kFALSE);
0189   for (int i = 0; i < num->GetSize(); i++) {
0190     double e, low, high;
0191     if (int(den->GetBinContent(i)) > 0.)
0192       e = double(num->GetBinContent(i)) / double(den->GetBinContent(i));
0193     else
0194       e = 0.;
0195     low = TEfficiency::Wilson((double)den->GetBinContent(i), (double)num->GetBinContent(i), 0.683, false);
0196     high = TEfficiency::Wilson((double)den->GetBinContent(i), (double)num->GetBinContent(i), 0.683, true);
0197 
0198     double err = e - low > high - e ? e - low : high - e;
0199     //here is the trick to store info in TProfile:
0200     eff->SetBinContent(i, e);
0201     eff->SetBinEntries(i, 1);
0202     eff->SetBinError(i, sqrt(e * e + err * err));
0203   }
0204   ibooker_.bookProfile2D(effName, eff);
0205   delete eff;
0206 }
0207 
0208 HeavyFlavorHarvesting::~HeavyFlavorHarvesting() {}
0209 
0210 //define this as a plug-in
0211 DEFINE_FWK_MODULE(HeavyFlavorHarvesting);