File indexing completed on 2023-10-25 09:55:56
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/interface/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 };
0038 }
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
0111 if (it->hwQual()) {
0112
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) {
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);
0126 }
0127 } else if (scenario_ == UseEmInterp::AllKeepHad) {
0128 float had_old = pt - cluster.emEt();
0129
0130 float em_new = it->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0131 float pt_new = had_old + em_new;
0132
0133 float hoe_new = em_new > 0 ? (had_old / em_new) : 999;
0134 cluster = l1t::PFCluster(pt_new, it->eta(), it->phi(), hoe_new, isEM);
0135
0136
0137 } else if (scenario_ == UseEmInterp::AllKeepTot) {
0138
0139 float em_new = it->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM);
0140 float hoe_new = em_new > 0 ? (it->pt() / em_new - 1) : -1;
0141 cluster = l1t::PFCluster(it->pt(), it->eta(), it->phi(), hoe_new, isEM);
0142
0143
0144 }
0145
0146 if (!emVsPUID_.method().empty()) {
0147 if (!emVsPUID_.passID(*it, cluster)) {
0148 continue;
0149 }
0150 }
0151 if (!emOnly_ && !emVsPionID_.method().empty()) {
0152 isEM = emVsPionID_.passID(*it, cluster);
0153 }
0154 cluster.setHwQual((isEM ? 1 : 0) + (it->hwQual() << 1));
0155
0156 if (corrector_.valid())
0157 corrector_.correctPt(cluster);
0158 cluster.setPtError(resol_(cluster.pt(), std::abs(cluster.eta())));
0159
0160
0161 cluster.setAbsZBarycenter(fabs(it->zBarycenter()));
0162 cluster.setSigmaRR(it->sigmaRRTot());
0163
0164 out->push_back(cluster);
0165 out->back().addConstituent(edm::Ptr<l1t::L1Candidate>(multiclusters, multiclusters->key(it)));
0166 if (hasEmId_) {
0167 (isEM ? outEm : outHad)->push_back(out->back());
0168 }
0169 }
0170
0171 iEvent.put(std::move(out));
0172 iEvent.put(std::move(outEgamma), "egamma");
0173 if (hasEmId_) {
0174 iEvent.put(std::move(outEm), "em");
0175 iEvent.put(std::move(outHad), "had");
0176 }
0177 }
0178 using l1tpf::PFClusterProducerFromHGC3DClusters;
0179 DEFINE_FWK_MODULE(PFClusterProducerFromHGC3DClusters);