RecoParticleFlow/plugins/PFJetDQMPostProcessor.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Validation/RecoParticleFlow
0004 // Class:      PFJetDQMPostProcessor.cc
0005 //
0006 // Main Developer:   "Juska Pekkanen"
0007 // Original Author:  "Kenichi Hatakeyama"
0008 //
0009
0010 #include "FWCore/Framework/interface/LuminosityBlock.h"
0011 #include "FWCore/Framework/interface/Run.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016
0017 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0018 #include "DQMServices/Core/interface/DQMStore.h"
0019
0020 #include "TCanvas.h"
0021
0022 //
0023 // class declaration
0024 //
0025
0026 class PFJetDQMPostProcessor : public DQMEDHarvester {
0027 public:
0028   explicit PFJetDQMPostProcessor(const edm::ParameterSet&);
0029   ~PFJetDQMPostProcessor() override;
0030
0031 private:
0032   void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0033   void fitResponse(TH1F* hreso,
0034                    TH1F* h_genjet_pt,
0035                    int ptbinlow,
0036                    int ietahigh,
0037                    double recoptcut,
0038                    double& resp,
0039                    double& resp_err,
0040                    double& reso,
0041                    double& reso_err);
0042   double getRespUnc(double width, double width_err, double mean, double mean_err);
0043
0044   std::vector<std::string> jetResponseDir;
0045   std::string genjetDir;
0046   std::vector<double> ptBins;
0047   std::vector<double> etaBins;
0048
0049   double recoptcut;
0050
0051   bool debug = false;
0052
0053   std::string seta(double eta) {
0054     std::string seta = std::to_string(int(eta * 10.));
0055     if (seta.length() < 2)
0056       seta = "0" + seta;
0057     return seta;
0058   }
0059
0060   std::string spt(double ptmin, double ptmax) {
0061     std::string spt = std::to_string(int(ptmin)) + "_" + std::to_string(int(ptmax));
0062     return spt;
0063   }
0064 };
0065
0066 // Some switches
0067 //
0068 // constuctors and destructor
0069 //
0070 PFJetDQMPostProcessor::PFJetDQMPostProcessor(const edm::ParameterSet& iConfig) {
0071   jetResponseDir = iConfig.getParameter<std::vector<std::string>>("jetResponseDir");
0072   genjetDir = iConfig.getParameter<std::string>("genjetDir");
0073   ptBins = iConfig.getParameter<std::vector<double>>("ptBins");
0074   etaBins = iConfig.getParameter<std::vector<double>>("etaBins");
0075   recoptcut = iConfig.getParameter<double>("recoPtCut");
0076 }
0077
0078 PFJetDQMPostProcessor::~PFJetDQMPostProcessor() {}
0079
0080 // ------------ method called right after a run ends ------------
0081 void PFJetDQMPostProcessor::dqmEndJob(DQMStore::IBooker& ibook_, DQMStore::IGetter& iget_) {
0082   for (unsigned int idir = 0; idir < jetResponseDir.size(); idir++) {
0083     iget_.setCurrentFolder(genjetDir);
0084     std::vector<std::string> sME_genjets = iget_.getMEs();
0085     std::for_each(sME_genjets.begin(), sME_genjets.end(), [&](auto& s) { s.insert(0, genjetDir); });
0086     //for (unsigned int i=0; i<sME_genjets.size(); i++) std::cout << sME_genjets[i] << std::endl;
0087
0088     iget_.setCurrentFolder(jetResponseDir[idir]);
0089     std::vector<std::string> sME_response = iget_.getMEs();
0090     std::for_each(sME_response.begin(), sME_response.end(), [&](auto& s) { s.insert(0, jetResponseDir[idir]); });
0091     //for (unsigned int i=0; i<sME_response.size(); i++) std::cout << sME_response[i] << std::endl;
0092
0093     iget_.setCurrentFolder(jetResponseDir[idir]);
0094
0095     double ptBinsArray[ptBins.size()];
0096     unsigned int nPtBins = ptBins.size() - 1;
0097     std::copy(ptBins.begin(), ptBins.end(), ptBinsArray);
0098     //for(unsigned int ipt = 0; ipt < ptBins.size(); ++ipt) std::cout << ptBins[ipt] << std::endl;
0099
0100     std::string stitle;
0101     char ctitle[50];
0102     std::vector<MonitorElement*> vME_presponse;
0103     std::vector<MonitorElement*> vME_preso;
0104     std::vector<MonitorElement*> vME_preso_rms;
0105
0106     MonitorElement* me;
0107     TH1F* h_resp;
0108     TH1F* h_genjet_pt;
0109
0110     //
0111     // Response distributions
0112     //
0113     for (unsigned int ieta = 1; ieta < etaBins.size(); ++ieta) {
0114       stitle = genjetDir + "genjet_pt" + "_eta" + seta(etaBins[ieta]);
0115       //std::cout << ieta << " " << stitle << std::endl;
0116
0117       std::vector<std::string>::const_iterator it = std::find(sME_genjets.begin(), sME_genjets.end(), stitle);
0118       if (it == sME_genjets.end())
0119         continue;
0120       me = iget_.get(stitle);
0121       h_genjet_pt = (TH1F*)me->getTH1F();
0122
0123       stitle = "presponse_eta" + seta(etaBins[ieta]);
0124       // adding "Raw" to the title of raw jet response histograms
0125       if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
0126         sprintf(ctitle, "Raw Jet pT response, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0127       } else {
0128         sprintf(ctitle, "Jet pT response, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0129       }
0130       TH1F* h_presponse = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0131
0132       stitle = "preso_eta" + seta(etaBins[ieta]);
0133       if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
0134         sprintf(ctitle, "Raw Jet pT resolution, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0135       } else {
0136         sprintf(ctitle, "Jet pT resolution, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0137       }
0138       TH1F* h_preso = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0139
0140       stitle = "preso_eta" + seta(etaBins[ieta]) + "_rms";
0141       if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
0142         sprintf(ctitle, "Raw Jet pT resolution using RMS, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0143       } else {
0144         sprintf(ctitle, "Jet pT resolution using RMS, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0145       }
0146       TH1F* h_preso_rms = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0147
0148       for (unsigned int ipt = 0; ipt < ptBins.size() - 1; ++ipt) {
0149         stitle = jetResponseDir[idir] + "reso_dist_" + spt(ptBins[ipt], ptBins[ipt + 1]) + "_eta" + seta(etaBins[ieta]);
0150         std::vector<std::string>::const_iterator it = std::find(sME_response.begin(), sME_response.end(), stitle);
0151         if (it == sME_response.end())
0152           continue;
0153         me = iget_.get(stitle);
0154         h_resp = (TH1F*)me->getTH1F();
0155
0156         // Fit-based
0157         double resp = 1.0, resp_err = 0.0, reso = 0.0, reso_err = 0.0;
0158         fitResponse(h_resp, h_genjet_pt, ipt, ieta, recoptcut, resp, resp_err, reso, reso_err);
0159
0160         h_presponse->SetBinContent(ipt + 1, resp);
0161         h_presponse->SetBinError(ipt + 1, resp_err);
0162         h_preso->SetBinContent(ipt + 1, reso);
0163         h_preso->SetBinError(ipt + 1, reso_err);
0164
0165         // RMS-based
0166         double std = h_resp->GetStdDev();
0167         double std_error = h_resp->GetStdDevError();
0168
0169         // Scale each bin with mean response
0170         double mean = 1.0;
0171         double mean_error = 0.0;
0172         double err = 0.0;
0173         if (h_resp->GetMean() > 0) {
0174           mean = h_resp->GetMean();
0175           mean_error = h_resp->GetMeanError();
0176
0177           // Scale resolution by response.
0178           std /= mean;
0179           std_error /= mean;
0180
0181           err = getRespUnc(std, std_error, mean, mean_error);
0182         }
0183
0184         h_preso_rms->SetBinContent(ipt + 1, std);
0185         h_preso_rms->SetBinError(ipt + 1, err);
0186
0187       }  // ipt
0188
0189       stitle = "presponse_eta" + seta(etaBins[ieta]);
0190       me = ibook_.book1D(stitle.c_str(), h_presponse);
0191       vME_presponse.push_back(me);
0192
0193       stitle = "preso_eta" + seta(etaBins[ieta]);
0194       me = ibook_.book1D(stitle.c_str(), h_preso);
0195       vME_preso.push_back(me);
0196
0197       stitle = "preso_eta" + seta(etaBins[ieta]) + "_rms";
0198       me = ibook_.book1D(stitle.c_str(), h_preso_rms);
0199       vME_preso_rms.push_back(me);
0200
0201     }  // ieta
0202
0203     //
0204     // Checks
0205     //
0206     if (debug) {
0207       for (std::vector<MonitorElement*>::const_iterator i = vME_presponse.begin(); i != vME_presponse.end(); ++i)
0208         (*i)->getTH1F()->Print();
0209       for (std::vector<MonitorElement*>::const_iterator i = vME_preso.begin(); i != vME_preso.end(); ++i)
0210         (*i)->getTH1F()->Print();
0211       for (std::vector<MonitorElement*>::const_iterator i = vME_preso_rms.begin(); i != vME_preso_rms.end(); ++i)
0212         (*i)->getTH1F()->Print();
0213     }
0214   }
0215 }
0216
0217 void PFJetDQMPostProcessor::fitResponse(TH1F* hreso,
0218                                         TH1F* h_genjet_pt,
0219                                         int ptbinlow,
0220                                         int ietahigh,
0221                                         double recoptcut,
0222                                         double& resp,
0223                                         double& resp_err,
0224                                         double& reso,
0225                                         double& reso_err) {
0226   // This 'smartfitter' is converted from the original Python smart_fit() -function
0227   // implemented in test/helperFunctions.py. See that file for more commentary.
0228   // Juska 23 May 2019
0229
0230   // Only do plots if needed for debugging
0231   // NOTE a directory called 'debug' must exist in the working directory
0232   //
0233   bool doPlots = false;
0234
0235   double ptlow = ptBins[ptbinlow];
0236   double pthigh = ptBins[ptbinlow + 1];
0237
0238   // Take range by Mikko's advice: -1.5 and + 1.5 * RMS width
0239
0240   double rmswidth = hreso->GetStdDev();
0241   double rmsmean = hreso->GetMean();
0242   double fitlow = rmsmean - 1.5 * rmswidth;
0243   fitlow = TMath::Max(recoptcut / ptlow, fitlow);
0244   double fithigh = rmsmean + 1.5 * rmswidth;
0245
0246   TF1* fg = new TF1("mygaus", "gaus", fitlow, fithigh);
0247   TF1* fg2 = new TF1("fg2", "TMath::Gaus(x,[0],[1],true)*[2]", fitlow, fithigh);
0248
0249   hreso->Fit("mygaus", "RQN");
0250
0251   fg2->SetParameter(0, fg->GetParameter(1));
0252   fg2->SetParameter(1, fg->GetParameter(2));
0253
0254   // Extract ngenjet in the current pT bin from the genjet histo
0255   float ngenjet = h_genjet_pt->GetBinContent(ptbinlow + 1);
0256
0257   // Here the fit is forced to take the area of ngenjets.
0258   // The area is further normalized for the response histogram x-axis length
0259   // (3) and number of bins (100)
0260   fg2->FixParameter(2, ngenjet * 3. / 100.);
0261
0262   hreso->Fit("fg2", "RQN");
0263
0264   fitlow = fg2->GetParameter(0) - 1.5 * fg2->GetParameter(1);
0265   fitlow = TMath::Max(15. / ptlow, fitlow);
0266   fithigh = fg2->GetParameter(0) + 1.5 * fg2->GetParameter(1);
0267
0268   fg2->SetRange(fitlow, fithigh);
0269
0270   hreso->Fit("fg2", "RQ");
0271
0272   fg->SetRange(0, 3);
0273   fg2->SetRange(0, 3);
0274   fg->SetLineWidth(2);
0275   fg2->SetLineColor(kGreen + 2);
0276
0277   hreso->GetXaxis()->SetRangeUser(0, 2);
0278
0279   // Save plots to a subdirectory if asked (debug-directory must exist!)
0280   if (doPlots & (hreso->GetEntries() > 0)) {
0281     TCanvas* cfit = new TCanvas(Form("respofit_%i", int(ptlow)), "respofit", 600, 600);
0282     hreso->Draw("ehist");
0283     fg->Draw("same");
0284     fg2->Draw("same");
0285     cfit->SaveAs(
0286         Form("debug/respo_smartfit_%04d_%i_eta%s.pdf", (int)ptlow, (int)pthigh, seta(etaBins[ietahigh]).c_str()));
0287   }
0288
0289   resp = fg2->GetParameter(0);
0290   resp_err = fg2->GetParError(0);
0291
0292   // Scale resolution by response. Avoid division by zero.
0293   if (0 == resp)
0294     reso = 0;
0295   else
0296     reso = fg2->GetParameter(1) / resp;
0297   reso_err = fg2->GetParError(1) / resp;
0298
0299   reso_err = getRespUnc(reso, reso_err, resp, resp_err);
0300 }
0301
0302 // Calculate resolution uncertainty
0303 double PFJetDQMPostProcessor::getRespUnc(double width, double width_err, double mean, double mean_err) {
0304   if (0 == width || 0 == mean)
0305     return 0;
0306   return TMath::Sqrt(pow(width_err, 2) / pow(width, 2) + pow(mean_err, 2) / pow(mean, 2)) * width;
0307 }
0308
0309 #include "FWCore/Framework/interface/MakerMacros.h"
0310 DEFINE_FWK_MODULE(PFJetDQMPostProcessor);