HLTrigger/plugins/HLTGenValSource.cc

0001 //********************************************************************************
0002 //
0003 //  Description:
0004 //    Producing and filling histograms for generator-level HLT path efficiency histograms. Harvested by a HLTGenValClient.
0005 //
0006 // Implementation:
0007 //   Histograms for objects of a certain type are created for multiple paths chosen by the user: for all objects,
0008 //   and for objects passing filters in the path, determined by deltaR matching.
0009 //   Each HLTGenValSource can handle a single object type and any number of paths (and filters in them)
0010 //
0011 //  Author: Finn Labe, UHH, Jul. 2022
0012 //
0013 //********************************************************************************
0014
0015 // system include files
0016 #include <memory>
0017 #include <chrono>
0018 #include <ctime>
0019
0020 // user include files
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Utilities/interface/InputTag.h"
0026
0027 // including GenParticles
0028 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0029
0030 // icnluding GenMET
0031 #include "DataFormats/METReco/interface/GenMETCollection.h"
0032 #include "DataFormats/METReco/interface/GenMET.h"
0033
0034 // includes needed for histogram creation
0035 #include "FWCore/ServiceRegistry/interface/Service.h"
0036 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0037
0038 // FunctionDefs
0039 #include "DQMOffline/Trigger/interface/FunctionDefs.h"
0040
0041 // includes of histogram collection class
0042 #include "Validation/HLTrigger/interface/HLTGenValHistCollPath.h"
0043
0044 // DQMStore
0045 #include "DQMServices/Core/interface/DQMStore.h"
0046 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0047
0048 // trigger Event
0049 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0050 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0051
0052 // object that can be a GenJet, GenParticle or energy sum
0053 #include "Validation/HLTrigger/interface/HLTGenValObject.h"
0054
0055 class HLTGenValSource : public DQMEDAnalyzer {
0056 public:
0057   explicit HLTGenValSource(const edm::ParameterSet&);
0058   ~HLTGenValSource() override = default;
0059   HLTGenValSource(const HLTGenValSource&) = delete;
0060   HLTGenValSource& operator=(const HLTGenValSource&) = delete;
0061
0062   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0063
0064 private:
0065   void analyze(const edm::Event&, const edm::EventSetup&) override;
0066   void bookHistograms(DQMStore::IBooker&, edm::Run const& run, edm::EventSetup const& c) override;
0067   void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
0068
0069   // functions to get correct object collection for chosen object type
0070   std::vector<HLTGenValObject> getObjectCollection(const edm::Event&);
0071   std::vector<HLTGenValObject> getGenParticles(const edm::Event&);
0072   reco::GenParticle getLastCopyPreFSR(reco::GenParticle part);
0073   reco::GenParticle getLastCopy(reco::GenParticle part);
0074
0075   // ----------member data ---------------------------
0076
0077   // tokens to get collections
0078   const edm::EDGetTokenT<reco::GenParticleCollection> genParticleToken_;
0079   const edm::EDGetTokenT<reco::GenMETCollection> genMETToken_;
0080   const edm::EDGetTokenT<reco::GenJetCollection> ak4genJetToken_;
0081   const edm::EDGetTokenT<reco::GenJetCollection> ak8genJetToken_;
0082   const edm::EDGetTokenT<trigger::TriggerEvent> trigEventToken_;
0083
0084   // config strings/Psets
0085   std::string objType_;
0086   std::string dirName_;
0087   std::vector<edm::ParameterSet> histConfigs_;
0088   std::vector<edm::ParameterSet> histConfigs2D_;
0089   std::vector<edm::ParameterSet> binnings_;
0090   std::string hltProcessName_;
0091
0092   // constructing the info string, which will be written to the output file for display of information in the GUI
0093   // the string will have a JSON formating, thus starting here with the opening bracket, which will be close directly before saving to the root file
0094   std::string infoString_ = "{";
0095
0096   // histogram collection per path
0097   std::vector<HLTGenValHistCollPath> collectionPath_;
0098
0099   // HLT config provider/getter
0100   HLTConfigProvider hltConfig_;
0101
0102   // some miscellaneous member variables
0103   std::vector<std::string> hltPathsToCheck_;
0104   std::vector<std::string> hltPaths;
0105   std::vector<std::string> hltPathSpecificCuts;
0106   double dR2limit_;
0107   bool doOnlyLastFilter_;
0108 };
0109
0110 HLTGenValSource::HLTGenValSource(const edm::ParameterSet& iConfig)
0111     : genParticleToken_(consumes<reco::GenParticleCollection>(
0112           iConfig.getParameterSet("inputCollections").getParameter<edm::InputTag>("genParticles"))),
0113       genMETToken_(consumes<reco::GenMETCollection>(
0114           iConfig.getParameterSet("inputCollections").getParameter<edm::InputTag>("genMET"))),
0115       ak4genJetToken_(consumes<reco::GenJetCollection>(
0116           iConfig.getParameterSet("inputCollections").getParameter<edm::InputTag>("ak4GenJets"))),
0117       ak8genJetToken_(consumes<reco::GenJetCollection>(
0118           iConfig.getParameterSet("inputCollections").getParameter<edm::InputTag>("ak8GenJets"))),
0119       trigEventToken_(consumes<trigger::TriggerEvent>(
0120           iConfig.getParameterSet("inputCollections").getParameter<edm::InputTag>("TrigEvent"))) {
0121   // getting the histogram configurations
0122   histConfigs_ = iConfig.getParameterSetVector("histConfigs");
0123   histConfigs2D_ = iConfig.getParameterSetVector("histConfigs2D");
0124   binnings_ = iConfig.getParameterSetVector("binnings");
0125
0126   // getting all other configurations
0127   dirName_ = iConfig.getParameter<std::string>("dqmDirName");
0128   objType_ = iConfig.getParameter<std::string>("objType");
0129   dR2limit_ = iConfig.getParameter<double>("dR2limit");
0130   doOnlyLastFilter_ = iConfig.getParameter<bool>("doOnlyLastFilter");
0131   hltProcessName_ = iConfig.getParameter<std::string>("hltProcessName");
0132   hltPathsToCheck_ = iConfig.getParameter<std::vector<std::string>>("hltPathsToCheck");
0133 }
0134
0135 void HLTGenValSource::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0136   // writing general information to info JSON
0137   auto t = std::time(nullptr);
0138   auto tm = *std::localtime(&t);
0139
0140   // date and time of running this
0141   std::ostringstream timeStringStream;
0142   timeStringStream << std::put_time(&tm, "%d-%m-%Y %H-%M-%S");
0143   auto timeString = timeStringStream.str();
0144   infoString_ += "\"date & time\":\"" + timeString + "\",";
0145
0146   // CMSSW version
0147   std::string cmsswVersion = std::getenv("CMSSW_VERSION");
0148   infoString_ += std::string("\"CMSSW release\":\"") + cmsswVersion + "\",";
0149
0150   // Initialize hltConfig, for cross-checking whether chosen paths exist
0151   bool changedConfig;
0152   if (!hltConfig_.init(iRun, iSetup, hltProcessName_, changedConfig)) {
0153     edm::LogError("HLTGenValSource") << "Initialization of HLTConfigProvider failed!";
0154     return;
0155   }
0156
0157   // global tag
0158   infoString_ += std::string("\"global tag\":\"") + hltConfig_.globalTag() + "\",";
0159
0160   // confDB table name
0161   infoString_ += std::string("\"HLT ConfDB table\":\"") + hltConfig_.tableName() + "\",";
0162
0163   // Get the set of trigger paths we want to make plots for
0164   std::vector<std::string> notFoundPaths;
0165   for (auto const& pathToCheck : hltPathsToCheck_) {
0166     // It is possible to add additional requirements to each path, seperated by a colon from the path name
0167     // these additional requirements are split from the path name here
0168     std::string cleanedPathToCheck;
0169     std::string pathSpecificCuts = "";
0170     if (pathToCheck.find(':') != std::string::npos) {
0171       // splitting the string
0172       std::stringstream hltPathToCheckInputStream(pathToCheck);
0173       std::string hltPathToCheckInputSegment;
0174       std::vector<std::string> hltPathToCheckInputSeglist;
0175       while (std::getline(hltPathToCheckInputStream, hltPathToCheckInputSegment, ':')) {
0176         hltPathToCheckInputSeglist.push_back(hltPathToCheckInputSegment);
0177       }
0178
0179       // here, exactly two parts are expected
0180       if (hltPathToCheckInputSeglist.size() != 2)
0181         throw cms::Exception("InputError")
0182             << "Path string can not be properly split into path and cuts: please use exactly one colon!.\n";
0183
0184       // the first part is the name of the path
0185       cleanedPathToCheck = hltPathToCheckInputSeglist.at(0);
0186
0187       // second part are the cuts, to be parsed later
0188       pathSpecificCuts = hltPathToCheckInputSeglist.at(1);
0189
0190     } else {
0191       cleanedPathToCheck = pathToCheck;
0192     }
0193
0194     bool pathfound = false;
0195     for (auto const& pathFromConfig : hltConfig_.triggerNames()) {
0196       if (pathFromConfig.find(cleanedPathToCheck) != std::string::npos) {
0197         hltPaths.push_back(pathFromConfig);
0198
0199         // in case the path was added twice, we'll add a tag automatically
0200         int count = std::count(hltPaths.begin(), hltPaths.end(), pathFromConfig);
0201         if (count > 1) {
0202           pathSpecificCuts += std::string(",autotag=v") + std::to_string(count);
0203         }
0204         hltPathSpecificCuts.push_back(pathSpecificCuts);
0205         pathfound = true;
0206       }
0207     }
0208     if (!pathfound)
0209       notFoundPaths.push_back(cleanedPathToCheck);
0210   }
0211   if (!notFoundPaths.empty()) {
0212     // error handling in case some paths do not exist
0213     std::string notFoundPathsMessage = "";
0214     for (const auto& path : notFoundPaths)
0215       notFoundPathsMessage += "- " + path + "\n";
0216     edm::LogError("HLTGenValSource") << "The following paths could not be found and will not be used: \n"
0217                                      << notFoundPathsMessage << std::endl;
0218   }
0219
0220   // before creating the collections for each path, we'll store the needed configurations in a pset
0221   // we'll copy this base multiple times and add the respective path
0222   // most of these options are not needed in the pathColl, but in the filterColls created in the pathColl
0223   edm::ParameterSet pathCollConfig;
0224   pathCollConfig.addParameter<std::string>("objType", objType_);
0225   pathCollConfig.addParameter<double>("dR2limit", dR2limit_);
0226   pathCollConfig.addParameter<bool>("doOnlyLastFilter", doOnlyLastFilter_);
0227   pathCollConfig.addParameter<std::string>("hltProcessName", hltProcessName_);
0228
0229   // creating a histogram collection for each path
0230   for (const auto& path : hltPaths) {
0231     edm::ParameterSet pathCollConfigStep = pathCollConfig;
0232     pathCollConfigStep.addParameter<std::string>("triggerPath", path);
0233     collectionPath_.emplace_back(HLTGenValHistCollPath(pathCollConfigStep, hltConfig_));
0234   }
0235 }
0236
0237 // ------------ method called for each event  ------------
0238 void HLTGenValSource::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0239   // creating the collection of HLTGenValObjects
0240   const std::vector<HLTGenValObject> objects = getObjectCollection(iEvent);
0241
0242   // init triggerEvent, which is always needed
0243   edm::Handle<trigger::TriggerEvent> triggerEvent;
0244   iEvent.getByToken(trigEventToken_, triggerEvent);
0245
0246   // loop over all objects and fill hists
0247   for (const auto& object : objects) {
0248     for (auto& collection_path : collectionPath_) {
0249       collection_path.fillHists(object, triggerEvent);
0250     }
0251   }
0252 }
0253
0254 // ------------ method called once each job just before starting event loop  ------------
0255 void HLTGenValSource::bookHistograms(DQMStore::IBooker& iBooker, const edm::Run& run, const edm::EventSetup& setup) {
0256   iBooker.setCurrentFolder(dirName_);
0257
0258   if (infoString_.back() == ',')
0259     infoString_.pop_back();
0260   infoString_ += "}";  // adding the closing bracked to the JSON string
0261   iBooker.bookString("HLTGenValInfo", infoString_);
0262
0263   // booking all histograms
0264   for (long unsigned int i = 0; i < collectionPath_.size(); i++) {
0265     std::vector<edm::ParameterSet> histConfigs = histConfigs_;
0266     for (auto& histConfig : histConfigs) {
0267       histConfig.addParameter<std::string>("pathSpecificCuts", hltPathSpecificCuts.at(i));
0268       histConfig.addParameter<std::vector<edm::ParameterSet>>("binnings",
0269                                                               binnings_);  // passing along the user-defined binnings
0270     }
0271
0272     collectionPath_.at(i).bookHists(iBooker, histConfigs, histConfigs2D_);
0273   }
0274 }
0275
0276 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0277 void HLTGenValSource::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0278   edm::ParameterSetDescription desc;
0279
0280   // basic parameter strings
0281   desc.add<std::string>(
0282       "objType");  // this deliberately has no default, as this is the main thing the user needs to chose
0283   desc.add<std::vector<std::string>>(
0284       "hltPathsToCheck");  // this for the moment also has no default: maybe there can be some way to handle this later?
0285   desc.add<std::string>("dqmDirName", "HLTGenVal");
0286   desc.add<std::string>("hltProcessName", "HLT");
0287   desc.add<double>("dR2limit", 0.1);
0288   desc.add<bool>("doOnlyLastFilter", false);
0289
0290   // input collections, a PSet
0291   edm::ParameterSetDescription inputCollections;
0292   inputCollections.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
0293   inputCollections.add<edm::InputTag>("genMET", edm::InputTag("genMetTrue"));
0294   inputCollections.add<edm::InputTag>("ak4GenJets", edm::InputTag("ak4GenJets"));
0295   inputCollections.add<edm::InputTag>("ak8GenJets", edm::InputTag("ak8GenJets"));
0296   inputCollections.add<edm::InputTag>("TrigEvent", edm::InputTag("hltTriggerSummaryAOD"));
0297   desc.add<edm::ParameterSetDescription>("inputCollections", inputCollections);
0298
0299   // hist descriptors, which are a vector of PSets
0300
0301   // defining single histConfig
0302   // this is generally without default, but a default set of histConfigs is specified below
0303   edm::ParameterSetDescription histConfig;
0304   histConfig.add<std::string>("vsVar");
0305   histConfig.add<std::vector<double>>("binLowEdges");
0306   histConfig.addVPSet(
0307       "rangeCuts", VarRangeCut<HLTGenValObject>::makePSetDescription(), std::vector<edm::ParameterSet>());
0308
0309   // default set of histConfigs
0310   std::vector<edm::ParameterSet> histConfigDefaults;
0311
0312   edm::ParameterSet histConfigDefault0;
0313   histConfigDefault0.addParameter<std::string>("vsVar", "pt");
0314   std::vector<double> defaultPtBinning{0,  5,  10, 12.5, 15,  17.5, 20,  22.5, 25,  30,  35, 40,
0315                                        45, 50, 60, 80,   100, 150,  200, 250,  300, 350, 400};
0316   histConfigDefault0.addParameter<std::vector<double>>("binLowEdges", defaultPtBinning);
0317   histConfigDefaults.push_back(histConfigDefault0);
0318
0319   edm::ParameterSet histConfigDefault1;
0320   histConfigDefault1.addParameter<std::string>("vsVar", "eta");
0321   std::vector<double> defaultetaBinning{-10, -8, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 8, 10};
0322   histConfigDefault1.addParameter<std::vector<double>>("binLowEdges", defaultetaBinning);
0323   histConfigDefaults.push_back(histConfigDefault1);
0324
0325   desc.addVPSet("histConfigs", histConfig, histConfigDefaults);
0326
0327   // defining single histConfig2D
0328   edm::ParameterSetDescription histConfig2D;
0329   histConfig2D.add<std::string>("vsVarX");
0330   histConfig2D.add<std::string>("vsVarY");
0331   histConfig2D.add<std::vector<double>>("binLowEdgesX");
0332   histConfig2D.add<std::vector<double>>("binLowEdgesY");
0333
0334   // default set of histConfigs
0335   std::vector<edm::ParameterSet> histConfigDefaults2D;
0336
0337   edm::ParameterSet histConfigDefault2D0;
0338   histConfigDefault2D0.addParameter<std::string>("vsVarX", "pt");
0339   histConfigDefault2D0.addParameter<std::string>("vsVarY", "eta");
0340   histConfigDefault2D0.addParameter<std::vector<double>>("binLowEdgesX", defaultPtBinning);
0341   histConfigDefault2D0.addParameter<std::vector<double>>("binLowEdgesY", defaultetaBinning);
0342   histConfigDefaults2D.push_back(histConfigDefault2D0);
0343
0344   desc.addVPSet("histConfigs2D", histConfig2D, histConfigDefaults2D);
0345
0346   // binnings, which are vectors of PSets
0347   // there are no default for this
0348   edm::ParameterSetDescription binningConfig;
0349   binningConfig.add<std::string>("name");
0350   binningConfig.add<std::string>("vsVar");
0351   binningConfig.add<std::vector<double>>("binLowEdges");
0352
0353   // this by default is empty
0354   std::vector<edm::ParameterSet> binningConfigDefaults;
0355
0356   desc.addVPSet("binnings", binningConfig, binningConfigDefaults);
0357
0358   descriptions.addDefault(desc);
0359 }
0360
0361 // this method handles the different object types and collections that can be used for efficiency calculation
0362 std::vector<HLTGenValObject> HLTGenValSource::getObjectCollection(const edm::Event& iEvent) {
0363   std::vector<HLTGenValObject> objects;  // the vector of objects to be filled
0364
0365   // handle object type
0366   std::vector<std::string> implementedGenParticles = {"ele", "pho", "mu", "tau"};
0367   if (std::find(implementedGenParticles.begin(), implementedGenParticles.end(), objType_) !=
0368       implementedGenParticles.end()) {
0369     objects = getGenParticles(iEvent);
0370   } else if (objType_ == "AK4jet") {  // ak4 jets, using the ak4GenJets collection
0371     const auto& genJets = iEvent.getHandle(ak4genJetToken_);
0372     for (size_t i = 0; i < genJets->size(); i++) {
0373       const reco::GenJet p = (*genJets)[i];
0374       objects.emplace_back(p);
0375     }
0376   } else if (objType_ == "AK8jet") {  // ak8 jets, using the ak8GenJets collection
0377     const auto& genJets = iEvent.getHandle(ak8genJetToken_);
0378     for (size_t i = 0; i < genJets->size(); i++) {
0379       const reco::GenJet p = (*genJets)[i];
0380       objects.emplace_back(p);
0381     }
0382   } else if (objType_ == "AK4HT") {  // ak4-based HT, using the ak4GenJets collection
0383     const auto& genJets = iEvent.getHandle(ak4genJetToken_);
0384     if (!genJets->empty()) {
0385       double HTsum = 0.;
0386       for (const auto& genJet : *genJets) {
0387         if (genJet.pt() > 30 && std::abs(genJet.eta()) < 2.5)
0388           HTsum += genJet.pt();
0389       }
0390       if (HTsum > 0)
0391         objects.emplace_back(reco::Candidate::PolarLorentzVector(HTsum, 0, 0, 0));
0392     }
0393   } else if (objType_ == "AK8HT") {  // ak8-based HT, using the ak8GenJets collection
0394     const auto& genJets = iEvent.getHandle(ak8genJetToken_);
0395     if (!genJets->empty()) {
0396       double HTsum = 0.;
0397       for (const auto& genJet : *genJets) {
0398         if (genJet.pt() > 200 && std::abs(genJet.eta()) < 2.5)
0399           HTsum += genJet.pt();
0400       }
0401       if (HTsum > 0)
0402         objects.emplace_back(reco::Candidate::PolarLorentzVector(HTsum, 0, 0, 0));
0403     }
0404   } else if (objType_ == "MET") {  // MET, using genMET
0405     const auto& genMET = iEvent.getHandle(genMETToken_);
0406     if (!genMET->empty()) {
0407       auto genMETpt = (*genMET)[0].pt();
0408       objects.emplace_back(reco::Candidate::PolarLorentzVector(genMETpt, 0, 0, 0));
0409     }
0410   } else
0411     throw cms::Exception("InputError") << "Generator-level validation is not available for type " << objType_ << ".\n"
0412                                        << "Please check for a potential spelling error.\n";
0413
0414   return objects;
0415 }
0416
0417 // in case of GenParticles, a subset of the entire collection needs to be chosen
0418 std::vector<HLTGenValObject> HLTGenValSource::getGenParticles(const edm::Event& iEvent) {
0419   std::vector<HLTGenValObject> objects;  // vector to be filled
0420
0421   const auto& genParticles = iEvent.getHandle(genParticleToken_);  // getting all GenParticles
0422
0423   // we need to ge the ID corresponding to the desired GenParticle type
0424   int pdgID = -1;  // setting to -1 should not be needed, but prevents the compiler warning :)
0425   if (objType_ == "ele")
0426     pdgID = 11;
0427   else if (objType_ == "pho")
0428     pdgID = 22;
0429   else if (objType_ == "mu")
0430     pdgID = 13;
0431   else if (objType_ == "tau")
0432     pdgID = 15;
0433
0434   // main loop over GenParticles
0435   for (size_t i = 0; i < genParticles->size(); ++i) {
0436     const reco::GenParticle p = (*genParticles)[i];
0437
0438     // only GenParticles with correct ID
0439     if (std::abs(p.pdgId()) != pdgID)
0440       continue;
0441
0442     // checking if particle comes from "hard process"
0443     if (p.isHardProcess()) {
0444       // depending on the particle type, last particle before or after FSR is chosen
0445       if ((objType_ == "ele") || (objType_ == "pho"))
0446         objects.emplace_back(getLastCopyPreFSR(p));
0447       else if ((objType_ == "mu") || (objType_ == "tau"))
0448         objects.emplace_back(getLastCopy(p));
0449     }
0450   }
0451
0452   return objects;
0453 }
0454
0455 // function returning the last GenParticle in a decay chain before FSR
0456 reco::GenParticle HLTGenValSource::getLastCopyPreFSR(reco::GenParticle part) {
0457   const auto& daughters = part.daughterRefVector();
0458   if (daughters.size() == 1 && daughters.at(0)->pdgId() == part.pdgId())
0459     return getLastCopyPreFSR(*daughters.at(0).get());  // recursion, whooo
0460   else
0461     return part;
0462 }
0463
0464 // function returning the last GenParticle in a decay chain
0465 reco::GenParticle HLTGenValSource::getLastCopy(reco::GenParticle part) {
0466   for (const auto& daughter : part.daughterRefVector()) {
0467     if (daughter->pdgId() == part.pdgId())
0468       return getLastCopy(*daughter.get());
0469   }
0470   return part;
0471 }
0472
0473 //define this as a framework plug-in
0474 DEFINE_FWK_MODULE(HLTGenValSource);