Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc

0001 #include <memory>
0002
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include "DataFormats/L1TParticleFlow/interface/PFCluster.h"
0010 #include "L1Trigger/Phase2L1ParticleFlow/src/corrector.h"
0011 #include "L1Trigger/Phase2L1ParticleFlow/interface/ParametricResolution.h"
0012 #include "L1Trigger/Phase2L1ParticleFlow/interface/HGC3DClusterEgID.h"
0013 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0014 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0015
0016 namespace l1tpf {
0017   class PFClusterProducerFromHGC3DClusters : public edm::stream::EDProducer<> {
0018   public:
0019     explicit PFClusterProducerFromHGC3DClusters(const edm::ParameterSet &);
0020     ~PFClusterProducerFromHGC3DClusters() override {}
0021
0022   private:
0023     enum class UseEmInterp { No, EmOnly, AllKeepHad, AllKeepTot };
0024
0025     edm::EDGetTokenT<l1t::HGCalMulticlusterBxCollection> src_;
0026     UseEmInterp scenario_;
0027     bool emOnly_;
0028     double etCut_;
0029     StringCutObjectSelector<l1t::HGCalMulticluster> preEmId_;
0030     l1tpf::HGC3DClusterEgID emVsPionID_, emVsPUID_;
0031     bool hasEmId_;
0032     l1tpf::corrector corrector_;
0033     l1tpf::ParametricResolution resol_;
0034
0035     void produce(edm::Event &, const edm::EventSetup &) override;
0036
0037   };  // class
0038 }  // namespace l1tpf
0039
0040 l1tpf::PFClusterProducerFromHGC3DClusters::PFClusterProducerFromHGC3DClusters(const edm::ParameterSet &iConfig)
0041     : src_(consumes<l1t::HGCalMulticlusterBxCollection>(iConfig.getParameter<edm::InputTag>("src"))),
0042       scenario_(UseEmInterp::No),
0043       emOnly_(iConfig.getParameter<bool>("emOnly")),
0044       etCut_(iConfig.getParameter<double>("etMin")),
0045       preEmId_(iConfig.getParameter<std::string>("preEmId")),
0046       emVsPionID_(iConfig.getParameter<edm::ParameterSet>("emVsPionID")),
0047       emVsPUID_(iConfig.getParameter<edm::ParameterSet>("emVsPUID")),
0048       hasEmId_((iConfig.existsAs<std::string>("preEmId") && !iConfig.getParameter<std::string>("preEmId").empty()) ||
0049                !emVsPionID_.method().empty()),
0050       corrector_(iConfig.getParameter<std::string>("corrector"),
0051                  emOnly_ || iConfig.getParameter<std::string>("corrector").empty()
0052                      ? -1
0053                      : iConfig.getParameter<double>("correctorEmfMax")),
0054       resol_(iConfig.getParameter<edm::ParameterSet>("resol")) {
0055   if (!emVsPionID_.method().empty()) {
0056     emVsPionID_.prepareTMVA();
0057   }
0058   if (!emVsPUID_.method().empty()) {
0059     emVsPUID_.prepareTMVA();
0060   }
0061
0062   produces<l1t::PFClusterCollection>();
0063   produces<l1t::PFClusterCollection>("egamma");
0064   if (hasEmId_) {
0065     produces<l1t::PFClusterCollection>("em");
0066     produces<l1t::PFClusterCollection>("had");
0067   }
0068
0069   std::string scenario = iConfig.getParameter<std::string>("useEMInterpretation");
0070   if (scenario == "emOnly") {
0071     scenario_ = UseEmInterp::EmOnly;
0072   } else if (scenario == "allKeepHad") {
0073     scenario_ = UseEmInterp::AllKeepHad;
0074     if (emOnly_) {
0075       throw cms::Exception("Configuration", "Unsupported emOnly = True when useEMInterpretation is " + scenario);
0076     }
0077   } else if (scenario == "allKeepTot") {
0078     scenario_ = UseEmInterp::AllKeepTot;
0079     if (emOnly_) {
0080       throw cms::Exception("Configuration", "Unsupported emOnly = True when useEMInterpretation is " + scenario);
0081     }
0082   } else if (scenario != "no") {
0083     throw cms::Exception("Configuration", "Unsupported useEMInterpretation scenario " + scenario);
0084   }
0085 }
0086
0087 void l1tpf::PFClusterProducerFromHGC3DClusters::produce(edm::Event &iEvent, const edm::EventSetup &) {
0088   auto out = std::make_unique<l1t::PFClusterCollection>();
0089   auto outEgamma = std::make_unique<l1t::PFClusterCollection>();
0090   std::unique_ptr<l1t::PFClusterCollection> outEm, outHad;
0091   if (hasEmId_) {
0092     outEm = std::make_unique<l1t::PFClusterCollection>();
0093     outHad = std::make_unique<l1t::PFClusterCollection>();
0094   }
0095   edm::Handle<l1t::HGCalMulticlusterBxCollection> multiclusters;
0096   iEvent.getByToken(src_, multiclusters);
0097
0098   for (auto it = multiclusters->begin(0), ed = multiclusters->end(0); it != ed; ++it) {
0099     float pt = it->pt(), hoe = it->hOverE();
0100     bool isEM = hasEmId_ ? preEmId_(*it) : emOnly_;
0101     if (emOnly_) {
0102       if (hoe == -1)
0103         continue;
0104       pt /= (1 + hoe);
0105       hoe = 0;
0106     }
0107     if (pt <= etCut_)
0108       continue;
0109
0110     // this block below is to support the older EG emulators, and is not used in newer ones
0111     if (it->hwQual()) {  // this is the EG ID shipped with the HGC TPs
0112       // we use the EM interpretation of the cluster energy
0113       l1t::PFCluster egcluster(
0114           it->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM), it->eta(), it->phi(), hoe, false);
0115       egcluster.setHwQual(it->hwQual());
0116       egcluster.addConstituent(edm::Ptr<l1t::L1Candidate>(multiclusters, multiclusters->key(it)));
0117       outEgamma->push_back(egcluster);
0118     }
0119
0120     l1t::PFCluster cluster(pt, it->eta(), it->phi(), hoe);
0121     if (scenario_ == UseEmInterp::EmOnly) {  // for emID objs, use EM interp as pT and set H = 0
0122       if (isEM) {
0123         float pt_new = it->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0124         float hoe_new = 0.;
0125         cluster = l1t::PFCluster(pt_new, it->eta(), it->phi(), hoe_new, /*isEM=*/isEM);
0126       }
0127     } else if (scenario_ == UseEmInterp::AllKeepHad) {  // for all objs, replace EM part with EM interp, preserve H
0128       float had_old = pt - cluster.emEt();
0129       //float em_old = cluster.emEt();
0130       float em_new = it->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0131       float pt_new = had_old + em_new;
0132       float hoe_new = em_new > 0 ? (had_old / em_new) : -1;
0133       cluster = l1t::PFCluster(pt_new, it->eta(), it->phi(), hoe_new, /*isEM=*/isEM);
0134       //printf("Scenario %d: pt %7.2f eta %+5.3f em %7.2f, EMI %7.2f, h/e % 8.3f --> pt %7.2f, em %7.2f, h/e % 8.3f\n",
0135       //        2, pt, it->eta(), em_old, em_new, hoe, cluster.pt(), cluster.emEt(), cluster.hOverE());
0136     } else if (scenario_ == UseEmInterp::AllKeepTot) {  // for all objs, replace EM part with EM interp, preserve pT
0137       //float em_old = cluster.emEt();
0138       float em_new = it->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0139       float hoe_new = em_new > 0 ? (it->pt() / em_new - 1) : -1;
0140       cluster = l1t::PFCluster(it->pt(), it->eta(), it->phi(), hoe_new, /*isEM=*/isEM);
0141       //printf("Scenario %d: pt %7.2f eta %+5.3f em %7.2f, EMI %7.2f, h/e % 8.3f --> pt %7.2f, em %7.2f, h/e % 8.3f\n",
0142       //        3, pt, it->eta(), em_old, em_new, hoe, cluster.pt(), cluster.emEt(), cluster.hOverE());
0143     }
0144
0145     if (!emVsPUID_.method().empty()) {
0146       if (!emVsPUID_.passID(*it, cluster)) {
0147         continue;
0148       }
0149     }
0150     if (!emOnly_ && !emVsPionID_.method().empty()) {
0151       isEM = emVsPionID_.passID(*it, cluster);
0152     }
0153     cluster.setHwQual((isEM ? 1 : 0) + (it->hwQual() << 1));
0154
0155     if (corrector_.valid())
0156       corrector_.correctPt(cluster);
0157     cluster.setPtError(resol_(cluster.pt(), std::abs(cluster.eta())));
0158
0159     out->push_back(cluster);
0160     out->back().addConstituent(edm::Ptr<l1t::L1Candidate>(multiclusters, multiclusters->key(it)));
0161     if (hasEmId_) {
0162       (isEM ? outEm : outHad)->push_back(out->back());
0163     }
0164   }
0165
0166   iEvent.put(std::move(out));
0167   iEvent.put(std::move(outEgamma), "egamma");
0168   if (hasEmId_) {
0169     iEvent.put(std::move(outEm), "em");
0170     iEvent.put(std::move(outHad), "had");
0171   }
0172 }
0173 using l1tpf::PFClusterProducerFromHGC3DClusters;
0174 DEFINE_FWK_MODULE(PFClusterProducerFromHGC3DClusters);