Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:33:15 

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 #include "TGraphAsymmErrors.h"
0022 
0023 //
0024 // class declaration
0025 //
0026 
0027 class PFJetDQMPostProcessor : public DQMEDHarvester {
0028 public:
0029   explicit PFJetDQMPostProcessor(const edm::ParameterSet&);
0030   ~PFJetDQMPostProcessor() override;
0031 
0032 private:
0033   void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0034   void fitResponse(TH1F* hreso,
0035                    TH1F* h_genjet_pt,
0036                    int ptbinlow,
0037                    int ietahigh,
0038                    double recoptcut,
0039                    double& resp,
0040                    double& resp_err,
0041                    double& reso,
0042                    double& reso_err);
0043   double getRespUnc(double width, double width_err, double mean, double mean_err);
0044 
0045   std::vector<std::string> jetResponseDir;
0046   std::string genjetDir;
0047   std::string offsetDir;
0048   std::vector<double> ptBins;
0049   std::vector<double> etaBins;
0050 
0051   double recoptcut;
0052 
0053   bool debug = false;
0054 
0055   std::string seta(double eta) {
0056     std::string seta = std::to_string(int(eta * 10.));
0057     if (seta.length() < 2)
0058       seta = "0" + seta;
0059     return seta;
0060   }
0061 
0062   std::string spt(double ptmin, double ptmax) {
0063     std::string spt = std::to_string(int(ptmin)) + "_" + std::to_string(int(ptmax));
0064     return spt;
0065   }
0066 };
0067 
0068 // Some switches
0069 //
0070 // constuctors and destructor
0071 //
0072 PFJetDQMPostProcessor::PFJetDQMPostProcessor(const edm::ParameterSet& iConfig) {
0073   jetResponseDir = iConfig.getParameter<std::vector<std::string>>("jetResponseDir");
0074   genjetDir = iConfig.getParameter<std::string>("genjetDir");
0075   offsetDir = iConfig.getParameter<std::string>("offsetDir");
0076   ptBins = iConfig.getParameter<std::vector<double>>("ptBins");
0077   etaBins = iConfig.getParameter<std::vector<double>>("etaBins");
0078   recoptcut = iConfig.getParameter<double>("recoPtCut");
0079 }
0080 
0081 PFJetDQMPostProcessor::~PFJetDQMPostProcessor() {}
0082 
0083 // ------------ method called right after a run ends ------------
0084 void PFJetDQMPostProcessor::dqmEndJob(DQMStore::IBooker& ibook_, DQMStore::IGetter& iget_) {
0085   for (unsigned int idir = 0; idir < jetResponseDir.size(); idir++) {
0086     iget_.setCurrentFolder(genjetDir);
0087     std::vector<std::string> sME_genjets = iget_.getMEs();
0088     std::for_each(sME_genjets.begin(), sME_genjets.end(), [&](auto& s) { s.insert(0, genjetDir); });
0089 
0090     iget_.setCurrentFolder(offsetDir);
0091     std::vector<std::string> sME_offset = iget_.getMEs();
0092     std::for_each(sME_offset.begin(), sME_offset.end(), [&](auto& s) { s.insert(0, offsetDir); });
0093 
0094     iget_.setCurrentFolder(jetResponseDir[idir]);
0095     std::vector<std::string> sME_response = iget_.getMEs();
0096     std::for_each(sME_response.begin(), sME_response.end(), [&](auto& s) { s.insert(0, jetResponseDir[idir]); });
0097 
0098     iget_.setCurrentFolder(jetResponseDir[idir]);
0099 
0100     double ptBinsArray[ptBins.size()];
0101     unsigned int nPtBins = ptBins.size() - 1;
0102     std::copy(ptBins.begin(), ptBins.end(), ptBinsArray);
0103     //for(unsigned int ipt = 0; ipt < ptBins.size(); ++ipt) std::cout << ptBins[ipt] << std::endl;
0104 
0105     std::string stitle;
0106     char ctitle[50];
0107     std::vector<MonitorElement*> vME_presponse;
0108     std::vector<MonitorElement*> vME_presponse_mean;
0109     std::vector<MonitorElement*> vME_presponse_median;
0110     std::vector<MonitorElement*> vME_preso;
0111     std::vector<MonitorElement*> vME_preso_rms;
0112     std::vector<MonitorElement*> vME_efficiency;
0113     std::vector<MonitorElement*> vME_purity;
0114     std::vector<MonitorElement*> vME_ratePUJet;
0115 
0116     MonitorElement* me;
0117     TH1F* h_resp;
0118     TH1F *h_genjet_pt, *h_genjet_matched_pt = nullptr;
0119     TH1F *h_recojet_pt = nullptr, *h_recojet_matched_pt = nullptr, *h_recojet_unmatched_pt = nullptr;
0120 
0121     stitle = offsetDir + "mu";
0122     std::vector<std::string>::const_iterator it = std::find(sME_offset.begin(), sME_offset.end(), stitle);
0123     if (it == sME_offset.end())
0124       continue;
0125     me = iget_.get(stitle);
0126     int nEvents = ((TH1F*)me->getTH1F())->GetEntries();
0127     if (nEvents == 0)
0128       continue;
0129     iget_.setCurrentFolder(jetResponseDir[idir]);
0130 
0131     bool isNoJEC = (jetResponseDir[idir].find("noJEC") != std::string::npos);
0132     bool isJEC = false;
0133     if (!isNoJEC)
0134       isJEC = (jetResponseDir[idir].find("JEC") != std::string::npos);
0135 
0136     //
0137     // Response distributions
0138     //
0139     for (unsigned int ieta = 1; ieta < etaBins.size(); ++ieta) {
0140       stitle = genjetDir + "genjet_pt" + "_eta" + seta(etaBins[ieta]);
0141 
0142       std::vector<std::string>::const_iterator it = std::find(sME_genjets.begin(), sME_genjets.end(), stitle);
0143       if (it == sME_genjets.end())
0144         continue;
0145       me = iget_.get(stitle);
0146       h_genjet_pt = (TH1F*)me->getTH1F();
0147 
0148       if (isNoJEC) {
0149         // getting the histogram for matched gen jets
0150         stitle = jetResponseDir[idir] + "genjet_pt" + "_eta" + seta(etaBins[ieta]) + "_matched";
0151         me = iget_.get(stitle);
0152         h_genjet_matched_pt = (TH1F*)me->getTH1F();
0153 
0154         /*// getting the histogram for unmatched gen jets
0155         stitle = jetResponseDir[idir] + "genjet_pt" + "_eta" + seta(etaBins[ieta]) + "_unmatched";
0156         me = iget_.get(stitle);
0157         h_genjet_unmatched_pt = (TH1F*)me->getTH1F();*/
0158       }
0159       if (isJEC) {
0160         // getting the histogram for reco jets
0161         stitle = jetResponseDir[idir] + "recojet_pt" + "_eta" + seta(etaBins[ieta]);
0162         me = iget_.get(stitle);
0163         h_recojet_pt = (TH1F*)me->getTH1F();
0164 
0165         // getting the histogram for matched reco jets
0166         stitle = jetResponseDir[idir] + "recojet_pt" + "_eta" + seta(etaBins[ieta]) + "_matched";
0167         me = iget_.get(stitle);
0168         h_recojet_matched_pt = (TH1F*)me->getTH1F();
0169 
0170         // getting the histogram for unmatched reco jets
0171         stitle = jetResponseDir[idir] + "recojet_pt" + "_eta" + seta(etaBins[ieta]) + "_unmatched";
0172         me = iget_.get(stitle);
0173         h_recojet_unmatched_pt = (TH1F*)me->getTH1F();
0174       }
0175 
0176       stitle = "presponse_eta" + seta(etaBins[ieta]);
0177       // adding "Raw" to the title of raw jet response histograms
0178       if (isNoJEC) {
0179         sprintf(ctitle, "Raw Jet pT response, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0180       } else {
0181         sprintf(ctitle, "Jet pT response, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0182       }
0183       TH1F* h_presponse = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0184 
0185       stitle = "presponse_eta" + seta(etaBins[ieta]) + "_mean";
0186       // adding "Raw" to the title of raw jet response histograms
0187       if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
0188         sprintf(ctitle, "Raw Jet pT response using Mean, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0189       } else {
0190         sprintf(ctitle, "Jet pT response using Mean, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0191       }
0192       TH1F* h_presponse_mean = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0193 
0194       stitle = "presponse_eta" + seta(etaBins[ieta]) + "_median";
0195       // adding "Raw" to the title of raw jet response histograms
0196       if (isNoJEC) {
0197         sprintf(ctitle, "Raw Jet pT response using Med., %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0198       } else {
0199         sprintf(ctitle, "Jet pT response using Med., %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0200       }
0201       TH1F* h_presponse_median = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0202 
0203       stitle = "preso_eta" + seta(etaBins[ieta]);
0204       if (isNoJEC) {
0205         sprintf(ctitle, "Raw Jet pT resolution, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0206       } else {
0207         sprintf(ctitle, "Jet pT resolution, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0208       }
0209       TH1F* h_preso = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0210 
0211       stitle = "preso_eta" + seta(etaBins[ieta]) + "_rms";
0212       if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
0213         sprintf(ctitle, "Raw Jet pT resolution using RMS, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0214       } else {
0215         sprintf(ctitle, "Jet pT resolution using RMS, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0216       }
0217       TH1F* h_preso_rms = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0218 
0219       //Booking histogram for Jet Efficiency vs pT
0220       stitle = "efficiency_eta" + seta(etaBins[ieta]);
0221       sprintf(ctitle, "Efficiency, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0222       TH1F* h_efficiency = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0223 
0224       //Booking histogram for Jet Purity vs pT
0225       stitle = "purity_eta" + seta(etaBins[ieta]);
0226       sprintf(ctitle, "Purity, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0227       TH1F* h_purity = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0228 
0229       //Booking histogram for #PU jets vs pT
0230       stitle = "ratePUJet_eta" + seta(etaBins[ieta]);
0231       sprintf(ctitle, "PU Jet Rate, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
0232       TH1F* h_ratePUJet = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
0233 
0234       if (isNoJEC) {
0235         h_efficiency->Divide(h_genjet_matched_pt, h_genjet_pt, 1, 1, "B");
0236       }
0237       if (isJEC) {
0238         h_purity->Divide(h_recojet_matched_pt, h_recojet_pt, 1, 1, "B");
0239         h_ratePUJet = (TH1F*)h_recojet_unmatched_pt->Clone();
0240         h_ratePUJet->SetName("h_ratePUJet");
0241         h_ratePUJet->Scale(1. / double(nEvents));
0242       }
0243 
0244       for (unsigned int ipt = 0; ipt < ptBins.size() - 1; ++ipt) {
0245         stitle = jetResponseDir[idir] + "reso_dist_" + spt(ptBins[ipt], ptBins[ipt + 1]) + "_eta" + seta(etaBins[ieta]);
0246         std::vector<std::string>::const_iterator it = std::find(sME_response.begin(), sME_response.end(), stitle);
0247         if (it == sME_response.end())
0248           continue;
0249         me = iget_.get(stitle);
0250         h_resp = (TH1F*)me->getTH1F();
0251 
0252         // Fit-based
0253         double resp = 1.0, resp_err = 0.0, reso = 0.0, reso_err = 0.0;
0254         fitResponse(h_resp, h_genjet_pt, ipt, ieta, recoptcut, resp, resp_err, reso, reso_err);
0255 
0256         h_presponse->SetBinContent(ipt + 1, resp);
0257         h_presponse->SetBinError(ipt + 1, resp_err);
0258         h_preso->SetBinContent(ipt + 1, reso);
0259         h_preso->SetBinError(ipt + 1, reso_err);
0260 
0261         // RMS-based
0262         double std = h_resp->GetStdDev();
0263         double std_error = h_resp->GetStdDevError();
0264 
0265         // Scale each bin with mean response
0266         double mean = 1.0;
0267         double mean_error = 0.0;
0268         double err = 0.0;
0269         if (h_resp->GetMean() > 0) {
0270           mean = h_resp->GetMean();
0271           mean_error = h_resp->GetMeanError();
0272 
0273           // Scale resolution by response.
0274           std /= mean;
0275           std_error /= mean;
0276 
0277           err = getRespUnc(std, std_error, mean, mean_error);
0278         }
0279 
0280         h_preso_rms->SetBinContent(ipt + 1, std);
0281         h_preso_rms->SetBinError(ipt + 1, err);
0282 
0283         // Mean-based
0284         h_presponse_mean->SetBinContent(ipt + 1, h_resp->GetMean());
0285         h_presponse_mean->SetBinError(ipt + 1, h_resp->GetMeanError());
0286 
0287         // Median-based
0288         if (h_resp->GetEntries() > 0) {
0289           int numBins = h_resp->GetXaxis()->GetNbins();
0290           Double_t x1[numBins];
0291           Double_t y1[numBins];
0292           for (int i = 0; i < numBins; i++) {
0293             x1[i] = h_resp->GetBinCenter(i + 1);
0294             y1[i] = h_resp->GetBinContent(i + 1) > 0 ? h_resp->GetBinContent(i + 1) : 0.0;
0295           }
0296           const auto x = x1, y = y1;
0297           double median = TMath::Median(numBins, x, y);
0298           h_presponse_median->SetBinContent(ipt + 1, median);
0299           h_presponse_median->SetBinError(ipt + 1, 1.2533 * (h_resp->GetMeanError()));
0300         }
0301 
0302       }  // ipt
0303 
0304       if (isNoJEC) {
0305         stitle = "efficiency_eta" + seta(etaBins[ieta]);
0306         me = ibook_.book1D(stitle.c_str(), h_efficiency);
0307         vME_efficiency.push_back(me);
0308       }
0309 
0310       if (isJEC) {
0311         stitle = "purity_eta" + seta(etaBins[ieta]);
0312         me = ibook_.book1D(stitle.c_str(), h_purity);
0313         vME_purity.push_back(me);
0314 
0315         stitle = "ratePUJet_eta" + seta(etaBins[ieta]);
0316         me = ibook_.book1D(stitle.c_str(), h_ratePUJet);
0317         vME_ratePUJet.push_back(me);
0318       }
0319 
0320       stitle = "presponse_eta" + seta(etaBins[ieta]);
0321       me = ibook_.book1D(stitle.c_str(), h_presponse);
0322       vME_presponse.push_back(me);
0323 
0324       stitle = "presponse_eta" + seta(etaBins[ieta]) + "_mean";
0325       me = ibook_.book1D(stitle.c_str(), h_presponse_mean);
0326       vME_presponse_mean.push_back(me);
0327 
0328       stitle = "presponse_eta" + seta(etaBins[ieta]) + "_median";
0329       me = ibook_.book1D(stitle.c_str(), h_presponse_median);
0330       vME_presponse_median.push_back(me);
0331 
0332       stitle = "preso_eta" + seta(etaBins[ieta]);
0333       me = ibook_.book1D(stitle.c_str(), h_preso);
0334       vME_preso.push_back(me);
0335 
0336       stitle = "preso_eta" + seta(etaBins[ieta]) + "_rms";
0337       me = ibook_.book1D(stitle.c_str(), h_preso_rms);
0338       vME_preso_rms.push_back(me);
0339 
0340     }  // ieta
0341 
0342     //
0343     // Checks
0344     //
0345     if (debug) {
0346       if (isNoJEC) {
0347         for (std::vector<MonitorElement*>::const_iterator i = vME_efficiency.begin(); i != vME_efficiency.end(); ++i)
0348           (*i)->getTH1F()->Print();
0349       }
0350       if (isJEC) {
0351         for (std::vector<MonitorElement*>::const_iterator i = vME_purity.begin(); i != vME_purity.end(); ++i)
0352           (*i)->getTH1F()->Print();
0353         for (std::vector<MonitorElement*>::const_iterator i = vME_ratePUJet.begin(); i != vME_ratePUJet.end(); ++i)
0354           (*i)->getTH1F()->Print();
0355       }
0356 
0357       for (std::vector<MonitorElement*>::const_iterator i = vME_presponse.begin(); i != vME_presponse.end(); ++i)
0358         (*i)->getTH1F()->Print();
0359       for (std::vector<MonitorElement*>::const_iterator i = vME_presponse_mean.begin(); i != vME_presponse_mean.end();
0360            ++i)
0361         (*i)->getTH1F()->Print();
0362       for (std::vector<MonitorElement*>::const_iterator i = vME_presponse_median.begin();
0363            i != vME_presponse_median.end();
0364            ++i)
0365         (*i)->getTH1F()->Print();
0366       for (std::vector<MonitorElement*>::const_iterator i = vME_preso.begin(); i != vME_preso.end(); ++i)
0367         (*i)->getTH1F()->Print();
0368       for (std::vector<MonitorElement*>::const_iterator i = vME_preso_rms.begin(); i != vME_preso_rms.end(); ++i)
0369         (*i)->getTH1F()->Print();
0370     }
0371   }
0372 }
0373 
0374 void PFJetDQMPostProcessor::fitResponse(TH1F* hreso,
0375                                         TH1F* h_genjet_pt,
0376                                         int ptbinlow,
0377                                         int ietahigh,
0378                                         double recoptcut,
0379                                         double& resp,
0380                                         double& resp_err,
0381                                         double& reso,
0382                                         double& reso_err) {
0383   // This 'smartfitter' is converted from the original Python smart_fit() -function
0384   // implemented in test/helperFunctions.py. See that file for more commentary.
0385   // Juska 23 May 2019
0386 
0387   // Only do plots if needed for debugging
0388   // NOTE a directory called 'debug' must exist in the working directory
0389   //
0390   bool doPlots = false;
0391 
0392   double ptlow = ptBins[ptbinlow];
0393   double pthigh = ptBins[ptbinlow + 1];
0394 
0395   // Take range by Mikko's advice: -1.5 and + 1.5 * RMS width
0396 
0397   double rmswidth = hreso->GetStdDev();
0398   double rmsmean = hreso->GetMean();
0399   double fitlow = rmsmean - 1.5 * rmswidth;
0400   fitlow = TMath::Max(recoptcut / ptlow, fitlow);
0401   double fithigh = rmsmean + 1.5 * rmswidth;
0402 
0403   TF1* fg = new TF1("mygaus", "gaus", fitlow, fithigh);
0404   TF1* fg2 = new TF1("fg2", "TMath::Gaus(x,[0],[1],true)*[2]", fitlow, fithigh);
0405 
0406   hreso->Fit("mygaus", "RQN");
0407 
0408   fg2->SetParameter(0, fg->GetParameter(1));
0409   fg2->SetParameter(1, fg->GetParameter(2));
0410 
0411   // Extract ngenjet in the current pT bin from the genjet histo
0412   float ngenjet = h_genjet_pt->GetBinContent(ptbinlow + 1);
0413 
0414   // Here the fit is forced to take the area of ngenjets.
0415   // The area is further normalized for the response histogram x-axis length
0416   // (3) and number of bins (100)
0417   fg2->FixParameter(2, ngenjet * 3. / 100.);
0418 
0419   hreso->Fit("fg2", "RQN");
0420 
0421   fitlow = fg2->GetParameter(0) - 1.5 * fg2->GetParameter(1);
0422   fitlow = TMath::Max(recoptcut / ptlow, fitlow);
0423   fithigh = fg2->GetParameter(0) + 1.5 * fg2->GetParameter(1);
0424 
0425   fg2->SetRange(fitlow, fithigh);
0426 
0427   hreso->Fit("fg2", "RQ");
0428 
0429   fg->SetRange(0, 3);
0430   fg2->SetRange(0, 3);
0431   fg->SetLineWidth(2);
0432   fg2->SetLineColor(kGreen + 2);
0433 
0434   hreso->GetXaxis()->SetRangeUser(0, 2);
0435 
0436   // Save plots to a subdirectory if asked (debug-directory must exist!)
0437   if (doPlots & (hreso->GetEntries() > 0)) {
0438     TCanvas* cfit = new TCanvas(Form("respofit_%i", int(ptlow)), "respofit", 600, 600);
0439     hreso->Draw("ehist");
0440     fg->Draw("same");
0441     fg2->Draw("same");
0442     cfit->SaveAs(
0443         Form("debug/respo_smartfit_%04d_%i_eta%s.pdf", (int)ptlow, (int)pthigh, seta(etaBins[ietahigh]).c_str()));
0444   }
0445 
0446   resp = fg2->GetParameter(0);
0447   resp_err = fg2->GetParError(0);
0448 
0449   // Scale resolution by response. Avoid division by zero.
0450   if (0 == resp)
0451     reso = 0;
0452   else
0453     reso = fg2->GetParameter(1) / resp;
0454   reso_err = fg2->GetParError(1) / resp;
0455 
0456   reso_err = getRespUnc(reso, reso_err, resp, resp_err);
0457 }
0458 
0459 // Calculate resolution uncertainty
0460 double PFJetDQMPostProcessor::getRespUnc(double width, double width_err, double mean, double mean_err) {
0461   if (0 == width || 0 == mean)
0462     return 0;
0463   return TMath::Sqrt(pow(width_err, 2) / pow(width, 2) + pow(mean_err, 2) / pow(mean, 2)) * width;
0464 }
0465 
0466 #include "FWCore/Framework/interface/MakerMacros.h"
0467 DEFINE_FWK_MODULE(PFJetDQMPostProcessor);
[ 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!