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   [[clang::suppress]]
0148   std::string cmsswVersion = std::getenv("CMSSW_VERSION");
0149   infoString_ += std::string("\"CMSSW release\":\"") + cmsswVersion + "\",";
0150
0151   // Initialize hltConfig, for cross-checking whether chosen paths exist
0152   bool changedConfig;
0153   if (!hltConfig_.init(iRun, iSetup, hltProcessName_, changedConfig)) {
0154     edm::LogError("HLTGenValSource") << "Initialization of HLTConfigProvider failed!";
0155     return;
0156   }
0157
0158   // global tag
0159   infoString_ += std::string("\"global tag\":\"") + hltConfig_.globalTag() + "\",";
0160
0161   // confDB table name
0162   infoString_ += std::string("\"HLT ConfDB table\":\"") + hltConfig_.tableName() + "\",";
0163
0164   // Get the set of trigger paths we want to make plots for
0165   std::vector<std::string> notFoundPaths;
0166   for (auto const& pathToCheck : hltPathsToCheck_) {
0167     // It is possible to add additional requirements to each path, seperated by a colon from the path name
0168     // these additional requirements are split from the path name here
0169     std::string cleanedPathToCheck;
0170     std::string pathSpecificCuts = "";
0171     if (pathToCheck.find(':') != std::string::npos) {
0172       // splitting the string
0173       std::stringstream hltPathToCheckInputStream(pathToCheck);
0174       std::string hltPathToCheckInputSegment;
0175       std::vector<std::string> hltPathToCheckInputSeglist;
0176       while (std::getline(hltPathToCheckInputStream, hltPathToCheckInputSegment, ':')) {
0177         hltPathToCheckInputSeglist.push_back(hltPathToCheckInputSegment);
0178       }
0179
0180       // here, exactly two parts are expected
0181       if (hltPathToCheckInputSeglist.size() != 2)
0182         throw cms::Exception("InputError")
0183             << "Path string can not be properly split into path and cuts: please use exactly one colon!.\n";
0184
0185       // the first part is the name of the path
0186       cleanedPathToCheck = hltPathToCheckInputSeglist.at(0);
0187
0188       // second part are the cuts, to be parsed later
0189       pathSpecificCuts = hltPathToCheckInputSeglist.at(1);
0190
0191     } else {
0192       cleanedPathToCheck = pathToCheck;
0193     }
0194
0195     bool pathfound = false;
0196     for (auto const& pathFromConfig : hltConfig_.triggerNames()) {
0197       if (pathFromConfig.find(cleanedPathToCheck) != std::string::npos) {
0198         hltPaths.push_back(pathFromConfig);
0199
0200         // in case the path was added twice, we'll add a tag automatically
0201         int count = std::count(hltPaths.begin(), hltPaths.end(), pathFromConfig);
0202         if (count > 1) {
0203           pathSpecificCuts += std::string(",autotag=v") + std::to_string(count);
0204         }
0205         hltPathSpecificCuts.push_back(pathSpecificCuts);
0206         pathfound = true;
0207       }
0208     }
0209     if (!pathfound)
0210       notFoundPaths.push_back(cleanedPathToCheck);
0211   }
0212   if (!notFoundPaths.empty()) {
0213     // error handling in case some paths do not exist
0214     std::string notFoundPathsMessage = "";
0215     for (const auto& path : notFoundPaths)
0216       notFoundPathsMessage += "- " + path + "\n";
0217     edm::LogError("HLTGenValSource") << "The following paths could not be found and will not be used: \n"
0218                                      << notFoundPathsMessage << std::endl;
0219   }
0220
0221   // before creating the collections for each path, we'll store the needed configurations in a pset
0222   // we'll copy this base multiple times and add the respective path
0223   // most of these options are not needed in the pathColl, but in the filterColls created in the pathColl
0224   edm::ParameterSet pathCollConfig;
0225   pathCollConfig.addParameter<std::string>("objType", objType_);
0226   pathCollConfig.addParameter<double>("dR2limit", dR2limit_);
0227   pathCollConfig.addParameter<bool>("doOnlyLastFilter", doOnlyLastFilter_);
0228   pathCollConfig.addParameter<std::string>("hltProcessName", hltProcessName_);
0229
0230   // creating a histogram collection for each path
0231   for (const auto& path : hltPaths) {
0232     edm::ParameterSet pathCollConfigStep = pathCollConfig;
0233     pathCollConfigStep.addParameter<std::string>("triggerPath", path);
0234     collectionPath_.emplace_back(HLTGenValHistCollPath(pathCollConfigStep, hltConfig_));
0235   }
0236 }
0237
0238 // ------------ method called for each event  ------------
0239 void HLTGenValSource::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0240   // creating the collection of HLTGenValObjects
0241   const std::vector<HLTGenValObject> objects = getObjectCollection(iEvent);
0242
0243   // init triggerEvent, which is always needed
0244   edm::Handle<trigger::TriggerEvent> triggerEvent;
0245   iEvent.getByToken(trigEventToken_, triggerEvent);
0246
0247   // loop over all objects and fill hists
0248   for (const auto& object : objects) {
0249     for (auto& collection_path : collectionPath_) {
0250       collection_path.fillHists(object, triggerEvent);
0251     }
0252   }
0253 }
0254
0255 // ------------ method called once each job just before starting event loop  ------------
0256 void HLTGenValSource::bookHistograms(DQMStore::IBooker& iBooker, const edm::Run& run, const edm::EventSetup& setup) {
0257   iBooker.setCurrentFolder(dirName_);
0258
0259   if (infoString_.back() == ',')
0260     infoString_.pop_back();
0261   infoString_ += "}";  // adding the closing bracked to the JSON string
0262   iBooker.bookString("HLTGenValInfo", infoString_);
0263
0264   // booking all histograms
0265   for (long unsigned int i = 0; i < collectionPath_.size(); i++) {
0266     std::vector<edm::ParameterSet> histConfigs = histConfigs_;
0267     for (auto& histConfig : histConfigs) {
0268       histConfig.addParameter<std::string>("pathSpecificCuts", hltPathSpecificCuts.at(i));
0269       histConfig.addParameter<std::vector<edm::ParameterSet>>("binnings",
0270                                                               binnings_);  // passing along the user-defined binnings
0271     }
0272
0273     collectionPath_.at(i).bookHists(iBooker, histConfigs, histConfigs2D_);
0274   }
0275 }
0276
0277 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0278 void HLTGenValSource::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0279   edm::ParameterSetDescription desc;
0280
0281   // basic parameter strings
0282   desc.add<std::string>(
0283       "objType");  // this deliberately has no default, as this is the main thing the user needs to chose
0284   desc.add<std::vector<std::string>>(
0285       "hltPathsToCheck");  // this for the moment also has no default: maybe there can be some way to handle this later?
0286   desc.add<std::string>("dqmDirName", "HLTGenVal");
0287   desc.add<std::string>("hltProcessName", "HLT");
0288   desc.add<double>("dR2limit", 0.1);
0289   desc.add<bool>("doOnlyLastFilter", false);
0290
0291   // input collections, a PSet
0292   edm::ParameterSetDescription inputCollections;
0293   inputCollections.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
0294   inputCollections.add<edm::InputTag>("genMET", edm::InputTag("genMetTrue"));
0295   inputCollections.add<edm::InputTag>("ak4GenJets", edm::InputTag("ak4GenJets"));
0296   inputCollections.add<edm::InputTag>("ak8GenJets", edm::InputTag("ak8GenJets"));
0297   inputCollections.add<edm::InputTag>("TrigEvent", edm::InputTag("hltTriggerSummaryAOD"));
0298   desc.add<edm::ParameterSetDescription>("inputCollections", inputCollections);
0299
0300   // hist descriptors, which are a vector of PSets
0301
0302   // defining single histConfig
0303   // this is generally without default, but a default set of histConfigs is specified below
0304   edm::ParameterSetDescription histConfig;
0305   histConfig.add<std::string>("vsVar");
0306   histConfig.add<std::vector<double>>("binLowEdges");
0307   histConfig.addVPSet(
0308       "rangeCuts", VarRangeCut<HLTGenValObject>::makePSetDescription(), std::vector<edm::ParameterSet>());
0309
0310   // default set of histConfigs
0311   std::vector<edm::ParameterSet> histConfigDefaults;
0312
0313   edm::ParameterSet histConfigDefault0;
0314   histConfigDefault0.addParameter<std::string>("vsVar", "pt");
0315   std::vector<double> defaultPtBinning{0,  5,  10, 12.5, 15,  17.5, 20,  22.5, 25,  30,  35, 40,
0316                                        45, 50, 60, 80,   100, 150,  200, 250,  300, 350, 400};
0317   histConfigDefault0.addParameter<std::vector<double>>("binLowEdges", defaultPtBinning);
0318   histConfigDefaults.push_back(histConfigDefault0);
0319
0320   edm::ParameterSet histConfigDefault1;
0321   histConfigDefault1.addParameter<std::string>("vsVar", "eta");
0322   std::vector<double> defaultetaBinning{-10, -8, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 8, 10};
0323   histConfigDefault1.addParameter<std::vector<double>>("binLowEdges", defaultetaBinning);
0324   histConfigDefaults.push_back(histConfigDefault1);
0325
0326   desc.addVPSet("histConfigs", histConfig, histConfigDefaults);
0327
0328   // defining single histConfig2D
0329   edm::ParameterSetDescription histConfig2D;
0330   histConfig2D.add<std::string>("vsVarX");
0331   histConfig2D.add<std::string>("vsVarY");
0332   histConfig2D.add<std::vector<double>>("binLowEdgesX");
0333   histConfig2D.add<std::vector<double>>("binLowEdgesY");
0334
0335   // default set of histConfigs
0336   std::vector<edm::ParameterSet> histConfigDefaults2D;
0337
0338   edm::ParameterSet histConfigDefault2D0;
0339   histConfigDefault2D0.addParameter<std::string>("vsVarX", "pt");
0340   histConfigDefault2D0.addParameter<std::string>("vsVarY", "eta");
0341   histConfigDefault2D0.addParameter<std::vector<double>>("binLowEdgesX", defaultPtBinning);
0342   histConfigDefault2D0.addParameter<std::vector<double>>("binLowEdgesY", defaultetaBinning);
0343   histConfigDefaults2D.push_back(histConfigDefault2D0);
0344
0345   desc.addVPSet("histConfigs2D", histConfig2D, histConfigDefaults2D);
0346
0347   // binnings, which are vectors of PSets
0348   // there are no default for this
0349   edm::ParameterSetDescription binningConfig;
0350   binningConfig.add<std::string>("name");
0351   binningConfig.add<std::string>("vsVar");
0352   binningConfig.add<std::vector<double>>("binLowEdges");
0353
0354   // this by default is empty
0355   std::vector<edm::ParameterSet> binningConfigDefaults;
0356
0357   desc.addVPSet("binnings", binningConfig, binningConfigDefaults);
0358
0359   descriptions.addDefault(desc);
0360 }
0361
0362 // this method handles the different object types and collections that can be used for efficiency calculation
0363 std::vector<HLTGenValObject> HLTGenValSource::getObjectCollection(const edm::Event& iEvent) {
0364   std::vector<HLTGenValObject> objects;  // the vector of objects to be filled
0365
0366   // handle object type
0367   std::vector<std::string> implementedGenParticles = {"ele", "pho", "mu", "tau"};
0368   if (std::find(implementedGenParticles.begin(), implementedGenParticles.end(), objType_) !=
0369       implementedGenParticles.end()) {
0370     objects = getGenParticles(iEvent);
0371   } else if (objType_ == "AK4jet") {  // ak4 jets, using the ak4GenJets collection
0372     const auto& genJets = iEvent.getHandle(ak4genJetToken_);
0373     for (size_t i = 0; i < genJets->size(); i++) {
0374       const reco::GenJet p = (*genJets)[i];
0375       objects.emplace_back(p);
0376     }
0377   } else if (objType_ == "AK8jet") {  // ak8 jets, using the ak8GenJets collection
0378     const auto& genJets = iEvent.getHandle(ak8genJetToken_);
0379     for (size_t i = 0; i < genJets->size(); i++) {
0380       const reco::GenJet p = (*genJets)[i];
0381       objects.emplace_back(p);
0382     }
0383   } else if (objType_ == "AK4HT") {  // ak4-based HT, using the ak4GenJets collection
0384     const auto& genJets = iEvent.getHandle(ak4genJetToken_);
0385     if (!genJets->empty()) {
0386       double HTsum = 0.;
0387       for (const auto& genJet : *genJets) {
0388         if (genJet.pt() > 30 && std::abs(genJet.eta()) < 2.5)
0389           HTsum += genJet.pt();
0390       }
0391       if (HTsum > 0)
0392         objects.emplace_back(reco::Candidate::PolarLorentzVector(HTsum, 0, 0, 0));
0393     }
0394   } else if (objType_ == "AK8HT") {  // ak8-based HT, using the ak8GenJets collection
0395     const auto& genJets = iEvent.getHandle(ak8genJetToken_);
0396     if (!genJets->empty()) {
0397       double HTsum = 0.;
0398       for (const auto& genJet : *genJets) {
0399         if (genJet.pt() > 200 && std::abs(genJet.eta()) < 2.5)
0400           HTsum += genJet.pt();
0401       }
0402       if (HTsum > 0)
0403         objects.emplace_back(reco::Candidate::PolarLorentzVector(HTsum, 0, 0, 0));
0404     }
0405   } else if (objType_ == "MET") {  // MET, using genMET
0406     const auto& genMET = iEvent.getHandle(genMETToken_);
0407     if (!genMET->empty()) {
0408       auto genMETpt = (*genMET)[0].pt();
0409       objects.emplace_back(reco::Candidate::PolarLorentzVector(genMETpt, 0, 0, 0));
0410     }
0411   } else
0412     throw cms::Exception("InputError") << "Generator-level validation is not available for type " << objType_ << ".\n"
0413                                        << "Please check for a potential spelling error.\n";
0414
0415   return objects;
0416 }
0417
0418 // in case of GenParticles, a subset of the entire collection needs to be chosen
0419 std::vector<HLTGenValObject> HLTGenValSource::getGenParticles(const edm::Event& iEvent) {
0420   std::vector<HLTGenValObject> objects;  // vector to be filled
0421
0422   const auto& genParticles = iEvent.getHandle(genParticleToken_);  // getting all GenParticles
0423
0424   // we need to ge the ID corresponding to the desired GenParticle type
0425   int pdgID = -1;  // setting to -1 should not be needed, but prevents the compiler warning :)
0426   if (objType_ == "ele")
0427     pdgID = 11;
0428   else if (objType_ == "pho")
0429     pdgID = 22;
0430   else if (objType_ == "mu")
0431     pdgID = 13;
0432   else if (objType_ == "tau")
0433     pdgID = 15;
0434
0435   // main loop over GenParticles
0436   for (size_t i = 0; i < genParticles->size(); ++i) {
0437     const reco::GenParticle p = (*genParticles)[i];
0438
0439     // only GenParticles with correct ID
0440     if (std::abs(p.pdgId()) != pdgID)
0441       continue;
0442
0443     // checking if particle comes from "hard process"
0444     if (p.isHardProcess()) {
0445       // depending on the particle type, last particle before or after FSR is chosen
0446       if ((objType_ == "ele") || (objType_ == "pho"))
0447         objects.emplace_back(getLastCopyPreFSR(p));
0448       else if ((objType_ == "mu") || (objType_ == "tau"))
0449         objects.emplace_back(getLastCopy(p));
0450     }
0451   }
0452
0453   return objects;
0454 }
0455
0456 // function returning the last GenParticle in a decay chain before FSR
0457 reco::GenParticle HLTGenValSource::getLastCopyPreFSR(reco::GenParticle part) {
0458   const auto& daughters = part.daughterRefVector();
0459   if (daughters.size() == 1 && daughters.at(0)->pdgId() == part.pdgId())
0460     return getLastCopyPreFSR(*daughters.at(0).get());  // recursion, whooo
0461   else
0462     return part;
0463 }
0464
0465 // function returning the last GenParticle in a decay chain
0466 reco::GenParticle HLTGenValSource::getLastCopy(reco::GenParticle part) {
0467   for (const auto& daughter : part.daughterRefVector()) {
0468     if (daughter->pdgId() == part.pdgId())
0469       return getLastCopy(*daughter.get());
0470   }
0471   return part;
0472 }
0473
0474 //define this as a framework plug-in
0475 DEFINE_FWK_MODULE(HLTGenValSource);