Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-10-01 01:02:12 

0001 // system include files
0002 #include <memory>
0003 #include <algorithm>
0004 #include <fstream>
0005 #include <cstdio>
0006 
0007 // user include files
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/stream/EDProducer.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ParameterSet/interface/FileInPath.h"
0013 
0014 #include "DataFormats/Common/interface/View.h"
0015 #include "DataFormats/Common/interface/RefToPtr.h"
0016 #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
0017 #include "DataFormats/L1Trigger/interface/Vertex.h"
0018 #include "DataFormats/L1Trigger/interface/VertexWord.h"
0019 
0020 #include "DataFormats/Math/interface/deltaR.h"
0021 
0022 #include "L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/tkinput_ref.h"
0023 #include "L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/muonGmtToL1ct_ref.h"
0024 #include "L1Trigger/Phase2L1ParticleFlow/interface/regionizer/regionizer_base_ref.h"
0025 #include "L1Trigger/Phase2L1ParticleFlow/interface/regionizer/multififo_regionizer_ref.h"
0026 #include "L1Trigger/Phase2L1ParticleFlow/interface/regionizer/tdr_regionizer_ref.h"
0027 #include "L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo2hgc_ref.h"
0028 #include "L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo3_ref.h"
0029 #include "L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo_dummy_ref.h"
0030 #include "L1Trigger/Phase2L1ParticleFlow/interface/puppi/linpuppi_ref.h"
0031 #include "L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h"
0032 #include "L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo_common_ref.h"
0033 #include "L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegsorter_ref.h"
0034 #include "L1Trigger/Phase2L1ParticleFlow/interface/L1TCorrelatorLayer1PatternFileWriter.h"
0035 
0036 #include "DataFormats/L1TCorrelator/interface/TkElectron.h"
0037 #include "DataFormats/L1TCorrelator/interface/TkElectronFwd.h"
0038 #include "DataFormats/L1Trigger/interface/EGamma.h"
0039 #include "DataFormats/L1TCorrelator/interface/TkEm.h"
0040 #include "DataFormats/L1TCorrelator/interface/TkEmFwd.h"
0041 
0042 //--------------------------------------------------------------------------------------------------
0043 class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> {
0044 public:
0045   explicit L1TCorrelatorLayer1Producer(const edm::ParameterSet &);
0046   ~L1TCorrelatorLayer1Producer() override;
0047 
0048 private:
0049   edm::ParameterSet config_;
0050 
0051   bool hasTracks_;
0052   edm::EDGetTokenT<l1t::PFTrackCollection> tkCands_;
0053   float trkPt_;
0054   bool emuTkVtx_;
0055   edm::EDGetTokenT<std::vector<l1t::Vertex>> extTkVtx_;
0056   edm::EDGetTokenT<std::vector<l1t::VertexWord>> tkVtxEmu_;
0057   int nVtx_;
0058 
0059   edm::EDGetTokenT<l1t::SAMuonCollection> muCands_;  // standalone muons
0060 
0061   std::vector<edm::EDGetTokenT<l1t::PFClusterCollection>> emCands_;
0062   std::vector<edm::EDGetTokenT<l1t::PFClusterCollection>> hadCands_;
0063 
0064   float emPtCut_, hadPtCut_;
0065 
0066   l1ct::Event event_;
0067   std::unique_ptr<l1ct::TrackInputEmulator> trackInput_;
0068   std::unique_ptr<l1ct::GMTMuonDecoderEmulator> muonInput_;
0069   std::unique_ptr<l1ct::RegionizerEmulator> regionizer_;
0070   std::unique_ptr<l1ct::PFAlgoEmulatorBase> l1pfalgo_;
0071   std::unique_ptr<l1ct::LinPuppiEmulator> l1pualgo_;
0072   std::unique_ptr<l1ct::PFTkEGAlgoEmulator> l1tkegalgo_;
0073   std::unique_ptr<l1ct::PFTkEGSorterEmulator> l1tkegsorter_;
0074 
0075   bool writeEgSta_;
0076   // Region dump
0077   const std::string regionDumpName_;
0078   bool writeRawHgcalCluster_;
0079   std::fstream fRegionDump_;
0080   const std::vector<edm::ParameterSet> patternWriterConfigs_;
0081   std::vector<std::unique_ptr<L1TCorrelatorLayer1PatternFileWriter>> patternWriters_;
0082 
0083   // region of interest debugging
0084   float debugEta_, debugPhi_, debugR_;
0085 
0086   // these are used to link items back
0087   std::unordered_map<const l1t::PFCluster *, l1t::PFClusterRef> clusterRefMap_;
0088   std::unordered_map<const l1t::PFTrack *, l1t::PFTrackRef> trackRefMap_;
0089   std::unordered_map<const l1t::SAMuon *, l1t::PFCandidate::MuonRef> muonRefMap_;
0090 
0091   // main methods
0092   void beginStream(edm::StreamID) override;
0093   void endStream() override;
0094   void produce(edm::Event &, const edm::EventSetup &) override;
0095   void addUInt(unsigned int value, std::string iLabel, edm::Event &iEvent);
0096 
0097   void initSectorsAndRegions(const edm::ParameterSet &iConfig);
0098   void initEvent(const edm::Event &e);
0099   // add object, tracking references
0100   void addTrack(const l1t::PFTrack &t, l1t::PFTrackRef ref);
0101   void addMuon(const l1t::SAMuon &t, l1t::PFCandidate::MuonRef ref);
0102   void addHadCalo(const l1t::PFCluster &t, l1t::PFClusterRef ref);
0103   void addEmCalo(const l1t::PFCluster &t, l1t::PFClusterRef ref);
0104   // add objects in already-decoded format
0105   void addDecodedTrack(l1ct::DetectorSector<l1ct::TkObjEmu> &sec, const l1t::PFTrack &t);
0106   void addDecodedMuon(l1ct::DetectorSector<l1ct::MuObjEmu> &sec, const l1t::SAMuon &t);
0107   void addDecodedHadCalo(l1ct::DetectorSector<l1ct::HadCaloObjEmu> &sec, const l1t::PFCluster &t);
0108   void addDecodedEmCalo(l1ct::DetectorSector<l1ct::EmCaloObjEmu> &sec, const l1t::PFCluster &t);
0109 
0110   void addRawHgcalCluster(l1ct::DetectorSector<ap_uint<256>> &sec, const l1t::PFCluster &c);
0111 
0112   // fetching outputs
0113   std::unique_ptr<l1t::PFCandidateCollection> fetchHadCalo() const;
0114   std::unique_ptr<l1t::PFCandidateCollection> fetchEmCalo() const;
0115   std::unique_ptr<l1t::PFCandidateCollection> fetchTracks() const;
0116   std::unique_ptr<l1t::PFCandidateCollection> fetchPF() const;
0117   void putPuppi(edm::Event &iEvent) const;
0118 
0119   void putEgStaObjects(edm::Event &iEvent,
0120                        const std::string &egLablel,
0121                        std::vector<edm::Ref<BXVector<l1t::EGamma>>> &egsta_refs);
0122   void putEgObjects(edm::Event &iEvent,
0123                     const bool writeEgSta,
0124                     const std::vector<edm::Ref<BXVector<l1t::EGamma>>> &egsta_refs,
0125                     const std::string &tkEmLabel,
0126                     const std::string &tkEmPerBoardLabel,
0127                     const std::string &tkEleLabel,
0128                     const std::string &tkElePerBoardLabel) const;
0129 
0130   template <typename T>
0131   void setRefs_(l1t::PFCandidate &pf, const T &p) const;
0132 
0133   void doVertexings(std::vector<float> &pvdz) const;
0134   // for multiplicities
0135   enum InputType { caloType = 0, emcaloType = 1, trackType = 2, l1muType = 3 };
0136   static constexpr const char *inputTypeName[l1muType + 1] = {"Calo", "EmCalo", "TK", "Mu"};
0137   std::unique_ptr<std::vector<unsigned>> vecSecInput(InputType i) const;
0138   std::unique_ptr<std::vector<unsigned>> vecRegInput(InputType i) const;
0139   typedef l1ct::OutputRegion::ObjType OutputType;
0140   std::unique_ptr<std::vector<unsigned>> vecOutput(OutputType i, bool usePuppi) const;
0141   std::pair<unsigned int, unsigned int> totAndMax(const std::vector<unsigned> &perRegion) const;
0142 };
0143 
0144 //
0145 // constructors and destructor
0146 //
0147 L1TCorrelatorLayer1Producer::L1TCorrelatorLayer1Producer(const edm::ParameterSet &iConfig)
0148     : config_(iConfig),
0149       hasTracks_(!iConfig.getParameter<edm::InputTag>("tracks").label().empty()),
0150       tkCands_(hasTracks_ ? consumes<l1t::PFTrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"))
0151                           : edm::EDGetTokenT<l1t::PFTrackCollection>()),
0152       trkPt_(iConfig.getParameter<double>("trkPtCut")),
0153       muCands_(consumes<l1t::SAMuonCollection>(iConfig.getParameter<edm::InputTag>("muons"))),
0154       emPtCut_(iConfig.getParameter<double>("emPtCut")),
0155       hadPtCut_(iConfig.getParameter<double>("hadPtCut")),
0156       regionizer_(nullptr),
0157       l1pfalgo_(nullptr),
0158       l1pualgo_(nullptr),
0159       l1tkegalgo_(nullptr),
0160       l1tkegsorter_(nullptr),
0161       regionDumpName_(iConfig.getUntrackedParameter<std::string>("dumpFileName", "")),
0162       writeRawHgcalCluster_(iConfig.getUntrackedParameter<bool>("writeRawHgcalCluster", false)),
0163       patternWriterConfigs_(iConfig.getUntrackedParameter<std::vector<edm::ParameterSet>>(
0164           "patternWriters", std::vector<edm::ParameterSet>())),
0165       debugEta_(iConfig.getUntrackedParameter<double>("debugEta", 0)),
0166       debugPhi_(iConfig.getUntrackedParameter<double>("debugPhi", 0)),
0167       debugR_(iConfig.getUntrackedParameter<double>("debugR", -1)) {
0168   produces<l1t::PFCandidateCollection>("PF");
0169   produces<l1t::PFCandidateCollection>("Puppi");
0170   produces<l1t::PFCandidateRegionalOutput>("PuppiRegional");
0171 
0172   produces<l1t::PFCandidateCollection>("EmCalo");
0173   produces<l1t::PFCandidateCollection>("Calo");
0174   produces<l1t::PFCandidateCollection>("TK");
0175 #if 0  // LATER
0176   produces<l1t::PFCandidateCollection>("TKVtx");
0177 #endif
0178 
0179   for (const auto &tag : iConfig.getParameter<std::vector<edm::InputTag>>("emClusters")) {
0180     emCands_.push_back(consumes<l1t::PFClusterCollection>(tag));
0181   }
0182   for (const auto &tag : iConfig.getParameter<std::vector<edm::InputTag>>("hadClusters")) {
0183     hadCands_.push_back(consumes<l1t::PFClusterCollection>(tag));
0184   }
0185 
0186   if (hasTracks_) {
0187     const std::string &tkInAlgo = iConfig.getParameter<std::string>("trackInputConversionAlgo");
0188     if (tkInAlgo == "Emulator") {
0189       trackInput_ = std::make_unique<l1ct::TrackInputEmulator>(
0190           iConfig.getParameter<edm::ParameterSet>("trackInputConversionParameters"));
0191     } else if (tkInAlgo != "Ideal")
0192       throw cms::Exception("Configuration", "Unsupported trackInputConversionAlgo");
0193   }
0194 
0195   const std::string &muInAlgo = iConfig.getParameter<std::string>("muonInputConversionAlgo");
0196   if (muInAlgo == "Emulator") {
0197     muonInput_ = std::make_unique<l1ct::GMTMuonDecoderEmulator>(
0198         iConfig.getParameter<edm::ParameterSet>("muonInputConversionParameters"));
0199   } else if (muInAlgo != "Ideal")
0200     throw cms::Exception("Configuration", "Unsupported muonInputConversionAlgo");
0201 
0202   const std::string &regalgo = iConfig.getParameter<std::string>("regionizerAlgo");
0203   if (regalgo == "Ideal") {
0204     regionizer_ =
0205         std::make_unique<l1ct::RegionizerEmulator>(iConfig.getParameter<edm::ParameterSet>("regionizerAlgoParameters"));
0206   } else if (regalgo == "Multififo") {
0207     regionizer_ = std::make_unique<l1ct::MultififoRegionizerEmulator>(
0208         iConfig.getParameter<edm::ParameterSet>("regionizerAlgoParameters"));
0209   } else if (regalgo == "TDR") {
0210     regionizer_ = std::make_unique<l1ct::TDRRegionizerEmulator>(
0211         iConfig.getParameter<edm::ParameterSet>("regionizerAlgoParameters"));
0212   } else
0213     throw cms::Exception("Configuration", "Unsupported regionizerAlgo");
0214 
0215   const std::string &algo = iConfig.getParameter<std::string>("pfAlgo");
0216   if (algo == "PFAlgo3") {
0217     l1pfalgo_ = std::make_unique<l1ct::PFAlgo3Emulator>(iConfig.getParameter<edm::ParameterSet>("pfAlgoParameters"));
0218   } else if (algo == "PFAlgo2HGC") {
0219     l1pfalgo_ = std::make_unique<l1ct::PFAlgo2HGCEmulator>(iConfig.getParameter<edm::ParameterSet>("pfAlgoParameters"));
0220   } else if (algo == "PFAlgoDummy") {
0221     l1pfalgo_ =
0222         std::make_unique<l1ct::PFAlgoDummyEmulator>(iConfig.getParameter<edm::ParameterSet>("pfAlgoParameters"));
0223   } else
0224     throw cms::Exception("Configuration", "Unsupported pfAlgo");
0225 
0226   const std::string &pualgo = iConfig.getParameter<std::string>("puAlgo");
0227   if (pualgo == "LinearizedPuppi") {
0228     l1pualgo_ = std::make_unique<l1ct::LinPuppiEmulator>(iConfig.getParameter<edm::ParameterSet>("puAlgoParameters"));
0229   } else
0230     throw cms::Exception("Configuration", "Unsupported puAlgo");
0231 
0232   l1tkegalgo_ = std::make_unique<l1ct::PFTkEGAlgoEmulator>(
0233       l1ct::PFTkEGAlgoEmuConfig(iConfig.getParameter<edm::ParameterSet>("tkEgAlgoParameters")));
0234 
0235   l1tkegsorter_ =
0236       std::make_unique<l1ct::PFTkEGSorterEmulator>(iConfig.getParameter<edm::ParameterSet>("tkEgSorterParameters"));
0237 
0238   if (l1tkegalgo_->writeEgSta())
0239     produces<BXVector<l1t::EGamma>>("L1Eg");
0240   produces<l1t::TkElectronCollection>("L1TkEle");
0241   produces<l1t::TkElectronRegionalOutput>("L1TkElePerBoard");
0242   produces<l1t::TkEmCollection>("L1TkEm");
0243   produces<l1t::TkEmRegionalOutput>("L1TkEmPerBoard");
0244 
0245   emuTkVtx_ = iConfig.getParameter<bool>("vtxCollectionEmulation");
0246   if (emuTkVtx_) {
0247     tkVtxEmu_ = consumes<std::vector<l1t::VertexWord>>(iConfig.getParameter<edm::InputTag>("vtxCollection"));
0248   } else {
0249     extTkVtx_ = consumes<std::vector<l1t::Vertex>>(iConfig.getParameter<edm::InputTag>("vtxCollection"));
0250   }
0251   nVtx_ = iConfig.getParameter<int>("nVtx");
0252 
0253   const char *iprefix[4] = {"totNReg", "maxNReg", "totNSec", "maxNSec"};
0254   for (int i = 0; i <= l1muType; ++i) {
0255     for (int ip = 0; ip < 4; ++ip) {
0256       produces<unsigned int>(std::string(iprefix[ip]) + inputTypeName[i]);
0257     }
0258     produces<std::vector<unsigned>>(std::string("vecNReg") + inputTypeName[i]);
0259     produces<std::vector<unsigned>>(std::string("vecNSec") + inputTypeName[i]);
0260   }
0261   const char *oprefix[4] = {"totNPF", "maxNPF", "totNPuppi", "maxNPuppi"};
0262   for (int i = 0; i < l1ct::OutputRegion::nPFTypes; ++i) {
0263     for (int ip = 0; ip < 4; ++ip) {
0264       produces<unsigned int>(std::string(oprefix[ip]) + l1ct::OutputRegion::objTypeName[i]);
0265     }
0266     produces<std::vector<unsigned>>(std::string("vecNPF") + l1ct::OutputRegion::objTypeName[i]);
0267     produces<std::vector<unsigned>>(std::string("vecNPuppi") + l1ct::OutputRegion::objTypeName[i]);
0268   }
0269 
0270   initSectorsAndRegions(iConfig);
0271 }
0272 
0273 L1TCorrelatorLayer1Producer::~L1TCorrelatorLayer1Producer() {}
0274 
0275 void L1TCorrelatorLayer1Producer::beginStream(edm::StreamID id) {
0276   if (!regionDumpName_.empty()) {
0277     if (id == 0) {
0278       fRegionDump_.open(regionDumpName_.c_str(), std::ios::out | std::ios::binary);
0279     } else {
0280       edm::LogWarning("L1TCorrelatorLayer1Producer")
0281           << "Job running with multiple streams, but dump file will have only events on stream zero.";
0282     }
0283   }
0284   if (!patternWriterConfigs_.empty()) {
0285     if (id == 0) {
0286       for (const auto &pset : patternWriterConfigs_) {
0287         patternWriters_.emplace_back(std::make_unique<L1TCorrelatorLayer1PatternFileWriter>(pset, event_));
0288       }
0289     } else {
0290       edm::LogWarning("L1TCorrelatorLayer1Producer")
0291           << "Job running with multiple streams, but pattern files will be written only on stream zero.";
0292     }
0293   }
0294 }
0295 
0296 void L1TCorrelatorLayer1Producer::endStream() {
0297   for (auto &writer : patternWriters_) {
0298     writer->flush();
0299   }
0300 }
0301 
0302 // ------------ method called to produce the data  ------------
0303 void L1TCorrelatorLayer1Producer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0304   // clear the regions also at the beginning, in case one event didn't complete but the job continues on
0305   initEvent(iEvent);
0306 
0307   /// ------ READ TRACKS ----
0308   if (hasTracks_) {
0309     edm::Handle<l1t::PFTrackCollection> htracks;
0310     iEvent.getByToken(tkCands_, htracks);
0311     const auto &tracks = *htracks;
0312     for (unsigned int itk = 0, ntk = tracks.size(); itk < ntk; ++itk) {
0313       const auto &tk = tracks[itk];
0314       // adding objects to PF
0315       if (debugR_ > 0 && deltaR(tk.eta(), tk.phi(), debugEta_, debugPhi_) > debugR_)
0316         continue;
0317       if (tk.pt() > trkPt_) {
0318         addTrack(tk, l1t::PFTrackRef(htracks, itk));
0319       }
0320     }
0321   }
0322 
0323   /// ------ READ MUONS ----
0324   edm::Handle<l1t::SAMuonCollection> muons;
0325   iEvent.getByToken(muCands_, muons);
0326   for (unsigned int i = 0, n = muons->size(); i < n; ++i) {
0327     const l1t::SAMuon &mu = (*muons)[i];
0328     if (debugR_ > 0 && deltaR(mu.eta(), mu.phi(), debugEta_, debugPhi_) > debugR_)
0329       continue;
0330     addMuon(mu, l1t::PFCandidate::MuonRef(muons, i));
0331   }
0332   // ------ READ CALOS -----
0333   edm::Handle<l1t::PFClusterCollection> caloHandle;
0334   for (const auto &tag : emCands_) {
0335     iEvent.getByToken(tag, caloHandle);
0336     const auto &calos = *caloHandle;
0337     for (unsigned int ic = 0, nc = calos.size(); ic < nc; ++ic) {
0338       const auto &calo = calos[ic];
0339       if (debugR_ > 0 && deltaR(calo.eta(), calo.phi(), debugEta_, debugPhi_) > debugR_)
0340         continue;
0341       if (calo.pt() > emPtCut_)
0342         addEmCalo(calo, l1t::PFClusterRef(caloHandle, ic));
0343     }
0344   }
0345   for (const auto &tag : hadCands_) {
0346     iEvent.getByToken(tag, caloHandle);
0347     const auto &calos = *caloHandle;
0348     for (unsigned int ic = 0, nc = calos.size(); ic < nc; ++ic) {
0349       const auto &calo = calos[ic];
0350       if (debugR_ > 0 && deltaR(calo.eta(), calo.phi(), debugEta_, debugPhi_) > debugR_)
0351         continue;
0352       if (calo.pt() > hadPtCut_)
0353         addHadCalo(calo, l1t::PFClusterRef(caloHandle, ic));
0354     }
0355   }
0356 
0357   regionizer_->run(event_.decoded, event_.pfinputs);
0358 
0359   // First, get a copy of the discretized and corrected inputs, and write them out
0360   iEvent.put(fetchEmCalo(), "EmCalo");
0361   iEvent.put(fetchHadCalo(), "Calo");
0362   iEvent.put(fetchTracks(), "TK");
0363 
0364   // Then do the vertexing, and save it out
0365   std::vector<float> z0s;
0366   std::vector<std::pair<float, float>> ptsums;
0367   float z0 = 0;
0368   double ptsum = 0;
0369   l1t::VertexWord pvwd;
0370   // FIXME: collections seem to be already sorted
0371   if (emuTkVtx_) {
0372     edm::Handle<std::vector<l1t::VertexWord>> vtxEmuHandle;
0373     iEvent.getByToken(tkVtxEmu_, vtxEmuHandle);
0374     for (const auto &vtx : *vtxEmuHandle) {
0375       ptsums.push_back(std::pair<float, float>(vtx.pt(), vtx.z0()));
0376       if (ptsum == 0 || vtx.pt() > ptsum) {
0377         ptsum = vtx.pt();
0378         pvwd = vtx;
0379       }
0380     }
0381   } else {
0382     edm::Handle<std::vector<l1t::Vertex>> vtxHandle;
0383     iEvent.getByToken(extTkVtx_, vtxHandle);
0384     for (const auto &vtx : *vtxHandle) {
0385       ptsums.push_back(std::pair<float, float>(vtx.pt(), vtx.z0()));
0386       if (ptsum == 0 || vtx.pt() > ptsum) {
0387         ptsum = vtx.pt();
0388         z0 = vtx.z0();
0389       }
0390     }
0391     pvwd = l1t::VertexWord(1, z0, 1, ptsum, 1, 1, 1);
0392   }
0393   l1ct::PVObjEmu hwpv;
0394   hwpv.hwZ0 = l1ct::Scales::makeZ0(pvwd.z0());
0395   event_.pvs.push_back(hwpv);
0396   event_.pvs_emu.push_back(pvwd.vertexWord());
0397   //Do a quick histogram vertexing to get multiple vertices (Hack for the taus)
0398   if (nVtx_ > 1) {
0399     std::stable_sort(ptsums.begin(), ptsums.end(), [](const auto &a, const auto &b) { return a.first > b.first; });
0400     for (int i0 = 0; i0 < std::min(int(ptsums.size()), int(nVtx_)); i0++) {
0401       z0s.push_back(ptsums[i0].second);
0402     }
0403     for (unsigned int i = 1; i < z0s.size(); ++i) {
0404       l1ct::PVObjEmu hwpv;
0405       hwpv.hwZ0 = l1ct::Scales::makeZ0(z0s[i]);
0406       event_.pvs.push_back(hwpv);  //Skip emu
0407     }
0408   }
0409   // Then also save the tracks with a vertex cut
0410 #if 0
0411   iEvent.put(l1regions_.fetchTracks(/*ptmin=*/0.0, /*fromPV=*/true), "TKVtx");
0412 #endif
0413 
0414   // Then run PF in each region
0415   event_.out.resize(event_.pfinputs.size());
0416   for (unsigned int ir = 0, nr = event_.pfinputs.size(); ir < nr; ++ir) {
0417     l1pfalgo_->run(event_.pfinputs[ir], event_.out[ir]);
0418     l1pfalgo_->mergeNeutrals(event_.out[ir]);
0419     l1tkegalgo_->run(event_.pfinputs[ir], event_.out[ir]);
0420     l1tkegalgo_->runIso(event_.pfinputs[ir], event_.pvs, event_.out[ir]);
0421   }
0422 
0423   // Then run puppi (regionally)
0424   for (unsigned int ir = 0, nr = event_.pfinputs.size(); ir < nr; ++ir) {
0425     l1pualgo_->run(event_.pfinputs[ir], event_.pvs, event_.out[ir]);
0426     //l1pualgo_->runNeutralsPU(l1region, z0, -1., puGlobals);
0427   }
0428 
0429   // NOTE: This needs to happen before the EG sorting per board so that the EG objects
0430   // get a global reference to the EGSta before being mixed among differente regions
0431   std::vector<edm::Ref<BXVector<l1t::EGamma>>> egsta_refs;
0432   if (l1tkegalgo_->writeEgSta()) {
0433     putEgStaObjects(iEvent, "L1Eg", egsta_refs);
0434   }
0435 
0436   // l1tkegsorter_->setDebug(true);
0437   for (auto &board : event_.board_out) {
0438     l1tkegsorter_->run(event_.pfinputs, event_.out, board.region_index, board.egphoton);
0439     l1tkegsorter_->run(event_.pfinputs, event_.out, board.region_index, board.egelectron);
0440   }
0441 
0442   // save PF into the event
0443   iEvent.put(fetchPF(), "PF");
0444 
0445   // and save puppi
0446   putPuppi(iEvent);
0447 
0448   // save the EG objects
0449   putEgObjects(iEvent, l1tkegalgo_->writeEgSta(), egsta_refs, "L1TkEm", "L1TkEmPerBoard", "L1TkEle", "L1TkElePerBoard");
0450 
0451   // Then go do the multiplicities
0452   for (int i = 0; i <= l1muType; ++i) {
0453     auto vecInputs = vecSecInput(InputType(i));
0454     auto tm = totAndMax(*vecInputs);
0455     addUInt(tm.first, std::string("totNSec") + inputTypeName[i], iEvent);
0456     addUInt(tm.second, std::string("maxNSec") + inputTypeName[i], iEvent);
0457     iEvent.put(std::move(vecInputs), std::string("vecNSec") + inputTypeName[i]);
0458   }
0459   for (int i = 0; i <= l1muType; ++i) {
0460     auto vecInputs = vecRegInput(InputType(i));
0461     auto tm = totAndMax(*vecInputs);
0462     addUInt(tm.first, std::string("totNReg") + inputTypeName[i], iEvent);
0463     addUInt(tm.second, std::string("maxNReg") + inputTypeName[i], iEvent);
0464     iEvent.put(std::move(vecInputs), std::string("vecNReg") + inputTypeName[i]);
0465   }
0466   for (int i = 0; i < l1ct::OutputRegion::nPFTypes; ++i) {
0467     auto vecPF = vecOutput(OutputType(i), false);
0468     auto tmPF = totAndMax(*vecPF);
0469     addUInt(tmPF.first, std::string("totNPF") + l1ct::OutputRegion::objTypeName[i], iEvent);
0470     addUInt(tmPF.second, std::string("maxNPF") + l1ct::OutputRegion::objTypeName[i], iEvent);
0471     iEvent.put(std::move(vecPF), std::string("vecNPF") + l1ct::OutputRegion::objTypeName[i]);
0472     auto vecPuppi = vecOutput(OutputType(i), true);
0473     auto tmPuppi = totAndMax(*vecPuppi);
0474     addUInt(tmPuppi.first, std::string("totNPuppi") + l1ct::OutputRegion::objTypeName[i], iEvent);
0475     addUInt(tmPuppi.second, std::string("maxNPuppi") + l1ct::OutputRegion::objTypeName[i], iEvent);
0476     iEvent.put(std::move(vecPuppi), std::string("vecNPuppi") + l1ct::OutputRegion::objTypeName[i]);
0477   }
0478 
0479   if (fRegionDump_.is_open()) {
0480     event_.write(fRegionDump_);
0481   }
0482   for (auto &writer : patternWriters_) {
0483     writer->write(event_);
0484   }
0485 
0486   // finally clear the regions
0487   event_.clear();
0488 }
0489 
0490 void L1TCorrelatorLayer1Producer::addUInt(unsigned int value, std::string iLabel, edm::Event &iEvent) {
0491   iEvent.put(std::make_unique<unsigned>(value), iLabel);
0492 }
0493 
0494 void L1TCorrelatorLayer1Producer::initSectorsAndRegions(const edm::ParameterSet &iConfig) {
0495   // the track finder geometry is fixed
0496   unsigned int TF_phiSlices = 9;
0497   float TF_phiWidth = 2 * M_PI / TF_phiSlices;
0498   event_.decoded.track.clear();
0499   for (unsigned int ieta = 0, neta = 2; ieta < neta; ++ieta) {
0500     for (unsigned int iphi = 0; iphi < TF_phiSlices; ++iphi) {
0501       float phiCenter = reco::reduceRange(iphi * TF_phiWidth);
0502       event_.decoded.track.emplace_back((ieta ? 0. : -2.5), (ieta ? 2.5 : 0.0), phiCenter, TF_phiWidth);
0503       event_.raw.track.emplace_back((ieta ? 0. : -2.5), (ieta ? 2.5 : 0.0), phiCenter, TF_phiWidth);
0504     }
0505   }
0506 
0507   event_.decoded.emcalo.clear();
0508   event_.decoded.hadcalo.clear();
0509   event_.raw.hgcalcluster.clear();
0510 
0511   for (const edm::ParameterSet &preg : iConfig.getParameter<std::vector<edm::ParameterSet>>("caloSectors")) {
0512     std::vector<double> etaBoundaries = preg.getParameter<std::vector<double>>("etaBoundaries");
0513     if (!std::is_sorted(etaBoundaries.begin(), etaBoundaries.end()))
0514       throw cms::Exception("Configuration", "caloSectors.etaBoundaries not sorted\n");
0515     unsigned int phiSlices = preg.getParameter<uint32_t>("phiSlices");
0516     float phiWidth = 2 * M_PI / phiSlices;
0517     if (phiWidth > 2 * l1ct::Scales::maxAbsPhi())
0518       throw cms::Exception("Configuration", "caloSectors phi range too large for phi_t data type");
0519     double phiZero = preg.getParameter<double>("phiZero");
0520     for (unsigned int ieta = 0, neta = etaBoundaries.size() - 1; ieta < neta; ++ieta) {
0521       float etaWidth = etaBoundaries[ieta + 1] - etaBoundaries[ieta];
0522       if (etaWidth > 2 * l1ct::Scales::maxAbsEta())
0523         throw cms::Exception("Configuration", "caloSectors eta range too large for eta_t data type");
0524       for (unsigned int iphi = 0; iphi < phiSlices; ++iphi) {
0525         float phiCenter = reco::reduceRange(iphi * phiWidth + phiZero);
0526         event_.decoded.hadcalo.emplace_back(etaBoundaries[ieta], etaBoundaries[ieta + 1], phiCenter, phiWidth);
0527         event_.decoded.emcalo.emplace_back(etaBoundaries[ieta], etaBoundaries[ieta + 1], phiCenter, phiWidth);
0528         event_.raw.hgcalcluster.emplace_back(etaBoundaries[ieta], etaBoundaries[ieta + 1], phiCenter, phiWidth);
0529       }
0530     }
0531   }
0532 
0533   event_.decoded.muon.region = l1ct::PFRegionEmu(0., 0.);  // centered at (0,0)
0534   event_.raw.muon.region = l1ct::PFRegionEmu(0., 0.);      // centered at (0,0)
0535 
0536   event_.pfinputs.clear();
0537   for (const edm::ParameterSet &preg : iConfig.getParameter<std::vector<edm::ParameterSet>>("regions")) {
0538     std::vector<double> etaBoundaries = preg.getParameter<std::vector<double>>("etaBoundaries");
0539     if (!std::is_sorted(etaBoundaries.begin(), etaBoundaries.end()))
0540       throw cms::Exception("Configuration", "regions.etaBoundaries not sorted\n");
0541     unsigned int phiSlices = preg.getParameter<uint32_t>("phiSlices");
0542     float etaExtra = preg.getParameter<double>("etaExtra");
0543     float phiExtra = preg.getParameter<double>("phiExtra");
0544     float phiWidth = 2 * M_PI / phiSlices;
0545     for (unsigned int ieta = 0, neta = etaBoundaries.size() - 1; ieta < neta; ++ieta) {
0546       for (unsigned int iphi = 0; iphi < phiSlices; ++iphi) {
0547         float phiCenter = reco::reduceRange(iphi * phiWidth);  //align with L1 TrackFinder phi sector indexing
0548         event_.pfinputs.emplace_back(
0549             etaBoundaries[ieta], etaBoundaries[ieta + 1], phiCenter, phiWidth, etaExtra, phiExtra);
0550       }
0551     }
0552   }
0553 
0554   event_.board_out.clear();
0555   const std::vector<edm::ParameterSet> &board_params = iConfig.getParameter<std::vector<edm::ParameterSet>>("boards");
0556   event_.board_out.resize(board_params.size());
0557   for (unsigned int bidx = 0; bidx < board_params.size(); bidx++) {
0558     event_.board_out[bidx].region_index = board_params[bidx].getParameter<std::vector<unsigned int>>("regions");
0559     float etaBoard = 0.;
0560     float phiBoard = 0.;
0561     for (auto ridx : event_.board_out[bidx].region_index) {
0562       etaBoard += event_.pfinputs[ridx].region.floatEtaCenter();
0563       phiBoard += event_.pfinputs[ridx].region.floatPhiCenter();
0564     }
0565     event_.board_out[bidx].eta = etaBoard / event_.board_out[bidx].region_index.size();
0566     event_.board_out[bidx].phi = phiBoard / event_.board_out[bidx].region_index.size();
0567   }
0568 }
0569 
0570 void L1TCorrelatorLayer1Producer::initEvent(const edm::Event &iEvent) {
0571   event_.clear();
0572   event_.run = iEvent.id().run();
0573   event_.lumi = iEvent.id().luminosityBlock();
0574   event_.event = iEvent.id().event();
0575   clusterRefMap_.clear();
0576   trackRefMap_.clear();
0577   muonRefMap_.clear();
0578 }
0579 
0580 void L1TCorrelatorLayer1Producer::addTrack(const l1t::PFTrack &t, l1t::PFTrackRef ref) {
0581   auto &rawsectors = event_.raw.track;
0582   auto &sectors = event_.decoded.track;
0583   assert(sectors.size() == 18 && rawsectors.size() == 18);
0584   int isec = t.track()->phiSector() + (t.eta() >= 0 ? 9 : 0);
0585   rawsectors[isec].obj.push_back(t.trackWord().getTrackWord());
0586   addDecodedTrack(sectors[isec], t);
0587   trackRefMap_[&t] = ref;
0588 }
0589 void L1TCorrelatorLayer1Producer::addMuon(const l1t::SAMuon &mu, l1t::PFCandidate::MuonRef ref) {
0590   event_.raw.muon.obj.emplace_back(mu.word());
0591   addDecodedMuon(event_.decoded.muon, mu);
0592   muonRefMap_[&mu] = ref;
0593 }
0594 void L1TCorrelatorLayer1Producer::addHadCalo(const l1t::PFCluster &c, l1t::PFClusterRef ref) {
0595   int sidx = 0;
0596   for (auto &sec : event_.decoded.hadcalo) {
0597     if (sec.region.contains(c.eta(), c.phi())) {
0598       addDecodedHadCalo(sec, c);
0599       if (writeRawHgcalCluster_)
0600         addRawHgcalCluster(event_.raw.hgcalcluster[sidx], c);
0601     }
0602     sidx++;
0603   }
0604   clusterRefMap_[&c] = ref;
0605 }
0606 void L1TCorrelatorLayer1Producer::addEmCalo(const l1t::PFCluster &c, l1t::PFClusterRef ref) {
0607   for (auto &sec : event_.decoded.emcalo) {
0608     if (sec.region.contains(c.eta(), c.phi())) {
0609       addDecodedEmCalo(sec, c);
0610     }
0611   }
0612   clusterRefMap_[&c] = ref;
0613 }
0614 
0615 void L1TCorrelatorLayer1Producer::addDecodedTrack(l1ct::DetectorSector<l1ct::TkObjEmu> &sec, const l1t::PFTrack &t) {
0616   std::pair<l1ct::TkObjEmu, bool> tkAndSel;
0617   if (trackInput_) {
0618     tkAndSel = trackInput_->decodeTrack(t.trackWord().getTrackWord(), sec.region);
0619   } else {
0620     tkAndSel.first.hwPt = l1ct::Scales::makePtFromFloat(t.pt());
0621     tkAndSel.first.hwEta =
0622         l1ct::Scales::makeGlbEta(t.caloEta()) -
0623         sec.region.hwEtaCenter;  // important to enforce that the region boundary is on a discrete value
0624     tkAndSel.first.hwPhi = l1ct::Scales::makePhi(sec.region.localPhi(t.caloPhi()));
0625     tkAndSel.first.hwCharge = t.charge() > 0;
0626     tkAndSel.first.hwQuality = t.quality();
0627     tkAndSel.first.hwDEta = l1ct::Scales::makeEta(t.eta() - t.caloEta());
0628     tkAndSel.first.hwDPhi = l1ct::Scales::makePhi(std::abs(reco::deltaPhi(t.phi(), t.caloPhi())));
0629     tkAndSel.first.hwZ0 = l1ct::Scales::makeZ0(t.vertex().Z());
0630     tkAndSel.first.hwDxy = 0;
0631     tkAndSel.second = t.quality() > 0;
0632   }
0633   // CMSSW-only extra info
0634   tkAndSel.first.hwChi2 = round(t.chi2() * 10);
0635   tkAndSel.first.hwStubs = t.nStubs();
0636   tkAndSel.first.simPt = t.pt();
0637   tkAndSel.first.simCaloEta = t.caloEta();
0638   tkAndSel.first.simCaloPhi = t.caloPhi();
0639   tkAndSel.first.simVtxEta = t.eta();
0640   tkAndSel.first.simVtxPhi = t.phi();
0641   tkAndSel.first.simZ0 = t.vertex().Z();
0642   tkAndSel.first.simD0 = t.vertex().Rho();
0643   tkAndSel.first.src = &t;
0644   // If the track fails, we set its pT to zero, so that the decoded tracks are still aligned with the raw tracks
0645   // Downstream, the regionizer will just ignore zero-momentum tracks
0646   if (!tkAndSel.second)
0647     tkAndSel.first.hwPt = 0;
0648   sec.obj.push_back(tkAndSel.first);
0649 }
0650 
0651 void L1TCorrelatorLayer1Producer::addDecodedMuon(l1ct::DetectorSector<l1ct::MuObjEmu> &sec, const l1t::SAMuon &t) {
0652   l1ct::MuObjEmu mu;
0653   if (muonInput_) {
0654     mu = muonInput_->decode(t.word());
0655   } else {
0656     mu.hwPt = l1ct::Scales::makePtFromFloat(t.pt());
0657     mu.hwEta = l1ct::Scales::makeGlbEta(t.eta());  // IMPORTANT: input is in global coordinates!
0658     mu.hwPhi = l1ct::Scales::makeGlbPhi(t.phi());
0659     mu.hwCharge = !t.hwCharge();
0660     mu.hwQuality = t.hwQual() / 2;
0661     mu.hwDEta = 0;
0662     mu.hwDPhi = 0;
0663     mu.hwZ0 = l1ct::Scales::makeZ0(t.vertex().Z());
0664     mu.hwDxy = 0;  // Dxy not defined yet
0665   }
0666   mu.src = &t;
0667   sec.obj.push_back(mu);
0668 }
0669 
0670 void L1TCorrelatorLayer1Producer::addDecodedHadCalo(l1ct::DetectorSector<l1ct::HadCaloObjEmu> &sec,
0671                                                     const l1t::PFCluster &c) {
0672   l1ct::HadCaloObjEmu calo;
0673   calo.hwPt = l1ct::Scales::makePtFromFloat(c.pt());
0674   calo.hwEta = l1ct::Scales::makeGlbEta(c.eta()) -
0675                sec.region.hwEtaCenter;  // important to enforce that the region boundary is on a discrete value
0676   calo.hwPhi = l1ct::Scales::makePhi(sec.region.localPhi(c.phi()));
0677   calo.hwEmPt = l1ct::Scales::makePtFromFloat(c.emEt());
0678   calo.hwEmID = c.hwEmID();
0679   calo.src = &c;
0680   sec.obj.push_back(calo);
0681 }
0682 
0683 void L1TCorrelatorLayer1Producer::addRawHgcalCluster(l1ct::DetectorSector<ap_uint<256>> &sec, const l1t::PFCluster &c) {
0684   ap_uint<256> cwrd = 0;
0685   ap_uint<14> w_pt = round(c.pt() / 0.25);
0686   ap_uint<14> w_empt = round(c.emEt() / 0.25);
0687   constexpr float ETAPHI_LSB = M_PI / 720;
0688   ap_int<9> w_eta = round(sec.region.localEta(c.eta()) / ETAPHI_LSB);
0689   ap_int<9> w_phi = round(sec.region.localPhi(c.phi()) / ETAPHI_LSB);
0690   ap_uint<10> w_qual = c.hwQual();
0691 
0692   cwrd(13, 0) = w_pt;
0693   cwrd(27, 14) = w_empt;
0694   cwrd(72, 64) = w_eta;
0695   cwrd(81, 73) = w_phi;
0696   cwrd(115, 106) = w_qual;
0697 
0698   sec.obj.push_back(cwrd);
0699 }
0700 
0701 void L1TCorrelatorLayer1Producer::addDecodedEmCalo(l1ct::DetectorSector<l1ct::EmCaloObjEmu> &sec,
0702                                                    const l1t::PFCluster &c) {
0703   l1ct::EmCaloObjEmu calo;
0704   calo.hwPt = l1ct::Scales::makePtFromFloat(c.pt());
0705   calo.hwEta = l1ct::Scales::makeGlbEta(c.eta()) -
0706                sec.region.hwEtaCenter;  // important to enforce that the region boundary is on a discrete value
0707   calo.hwPhi = l1ct::Scales::makePhi(sec.region.localPhi(c.phi()));
0708   calo.hwPtErr = l1ct::Scales::makePtFromFloat(c.ptError());
0709   calo.hwEmID = c.hwEmID();
0710   calo.src = &c;
0711   sec.obj.push_back(calo);
0712 }
0713 
0714 template <typename T>
0715 void L1TCorrelatorLayer1Producer::setRefs_(l1t::PFCandidate &pf, const T &p) const {
0716   if (p.srcCluster) {
0717     auto match = clusterRefMap_.find(p.srcCluster);
0718     if (match == clusterRefMap_.end()) {
0719       throw cms::Exception("CorruptData") << "Invalid cluster pointer in PF candidate id " << p.intId() << " pt "
0720                                           << p.floatPt() << " eta " << p.floatEta() << " phi " << p.floatPhi();
0721     }
0722     pf.setPFCluster(match->second);
0723   }
0724   if (p.srcTrack) {
0725     auto match = trackRefMap_.find(p.srcTrack);
0726     if (match == trackRefMap_.end()) {
0727       throw cms::Exception("CorruptData") << "Invalid track pointer in PF candidate id " << p.intId() << " pt "
0728                                           << p.floatPt() << " eta " << p.floatEta() << " phi " << p.floatPhi();
0729     }
0730     pf.setPFTrack(match->second);
0731   }
0732   if (p.srcMu) {
0733     auto match = muonRefMap_.find(p.srcMu);
0734     if (match == muonRefMap_.end()) {
0735       throw cms::Exception("CorruptData") << "Invalid muon pointer in PF candidate id " << p.intId() << " pt "
0736                                           << p.floatPt() << " eta " << p.floatEta() << " phi " << p.floatPhi();
0737     }
0738     pf.setMuon(match->second);
0739   }
0740 }
0741 
0742 template <>
0743 void L1TCorrelatorLayer1Producer::setRefs_<l1ct::PFNeutralObjEmu>(l1t::PFCandidate &pf,
0744                                                                   const l1ct::PFNeutralObjEmu &p) const {
0745   if (p.srcCluster) {
0746     auto match = clusterRefMap_.find(p.srcCluster);
0747     if (match == clusterRefMap_.end()) {
0748       throw cms::Exception("CorruptData") << "Invalid cluster pointer in PF candidate id " << p.intId() << " pt "
0749                                           << p.floatPt() << " eta " << p.floatEta() << " phi " << p.floatPhi();
0750     }
0751     pf.setPFCluster(match->second);
0752   }
0753 }
0754 
0755 template <>
0756 void L1TCorrelatorLayer1Producer::setRefs_<l1ct::HadCaloObjEmu>(l1t::PFCandidate &pf,
0757                                                                 const l1ct::HadCaloObjEmu &p) const {
0758   if (p.src) {
0759     auto match = clusterRefMap_.find(p.src);
0760     if (match == clusterRefMap_.end()) {
0761       throw cms::Exception("CorruptData") << "Invalid cluster pointer in hadcalo candidate  pt " << p.floatPt()
0762                                           << " eta " << p.floatEta() << " phi " << p.floatPhi();
0763     }
0764     pf.setPFCluster(match->second);
0765   }
0766 }
0767 
0768 template <>
0769 void L1TCorrelatorLayer1Producer::setRefs_<l1ct::EmCaloObjEmu>(l1t::PFCandidate &pf,
0770                                                                const l1ct::EmCaloObjEmu &p) const {
0771   if (p.src) {
0772     auto match = clusterRefMap_.find(p.src);
0773     if (match == clusterRefMap_.end()) {
0774       throw cms::Exception("CorruptData") << "Invalid cluster pointer in emcalo candidate  pt " << p.floatPt()
0775                                           << " eta " << p.floatEta() << " phi " << p.floatPhi();
0776     }
0777     pf.setPFCluster(match->second);
0778   }
0779 }
0780 
0781 template <>
0782 void L1TCorrelatorLayer1Producer::setRefs_<l1ct::TkObjEmu>(l1t::PFCandidate &pf, const l1ct::TkObjEmu &p) const {
0783   if (p.src) {
0784     auto match = trackRefMap_.find(p.src);
0785     if (match == trackRefMap_.end()) {
0786       throw cms::Exception("CorruptData") << "Invalid track pointer in track candidate  pt " << p.floatPt() << " eta "
0787                                           << p.floatEta() << " phi " << p.floatPhi();
0788     }
0789     pf.setPFTrack(match->second);
0790   }
0791 }
0792 
0793 std::unique_ptr<l1t::PFCandidateCollection> L1TCorrelatorLayer1Producer::fetchHadCalo() const {
0794   auto ret = std::make_unique<l1t::PFCandidateCollection>();
0795   for (const auto &r : event_.pfinputs) {
0796     const auto &reg = r.region;
0797     for (const auto &p : r.hadcalo) {
0798       if (p.hwPt == 0 || !reg.isFiducial(p))
0799         continue;
0800       reco::Particle::PolarLorentzVector p4(p.floatPt(), reg.floatGlbEtaOf(p), reg.floatGlbPhiOf(p), 0.13f);
0801       l1t::PFCandidate::ParticleType type = p.hwIsEM() ? l1t::PFCandidate::Photon : l1t::PFCandidate::NeutralHadron;
0802       ret->emplace_back(type, 0, p4, 1, p.intPt(), p.intEta(), p.intPhi());
0803       ret->back().setHwEmID(p.hwEmID);
0804       setRefs_(ret->back(), p);
0805     }
0806   }
0807   return ret;
0808 }
0809 std::unique_ptr<l1t::PFCandidateCollection> L1TCorrelatorLayer1Producer::fetchEmCalo() const {
0810   auto ret = std::make_unique<l1t::PFCandidateCollection>();
0811   for (const auto &r : event_.pfinputs) {
0812     const auto &reg = r.region;
0813     for (const auto &p : r.emcalo) {
0814       if (p.hwPt == 0 || !reg.isFiducial(p))
0815         continue;
0816       reco::Particle::PolarLorentzVector p4(p.floatPt(), reg.floatGlbEtaOf(p), reg.floatGlbPhiOf(p), 0.13f);
0817       ret->emplace_back(l1t::PFCandidate::Photon, 0, p4, 1, p.intPt(), p.intEta(), p.intPhi());
0818       ret->back().setHwEmID(p.hwEmID);
0819       setRefs_(ret->back(), p);
0820     }
0821   }
0822   return ret;
0823 }
0824 std::unique_ptr<l1t::PFCandidateCollection> L1TCorrelatorLayer1Producer::fetchTracks() const {
0825   auto ret = std::make_unique<l1t::PFCandidateCollection>();
0826   for (const auto &r : event_.pfinputs) {
0827     const auto &reg = r.region;
0828     for (const auto &p : r.track) {
0829       if (p.hwPt == 0 || !reg.isFiducial(p))
0830         continue;
0831       reco::Particle::PolarLorentzVector p4(
0832           p.floatPt(), reg.floatGlbEta(p.hwVtxEta()), reg.floatGlbPhi(p.hwVtxPhi()), 0.13f);
0833       ret->emplace_back(l1t::PFCandidate::ChargedHadron, p.intCharge(), p4, 1, p.intPt(), p.intEta(), p.intPhi());
0834       setRefs_(ret->back(), p);
0835     }
0836   }
0837   return ret;
0838 }
0839 
0840 std::unique_ptr<l1t::PFCandidateCollection> L1TCorrelatorLayer1Producer::fetchPF() const {
0841   auto ret = std::make_unique<l1t::PFCandidateCollection>();
0842   for (unsigned int ir = 0, nr = event_.pfinputs.size(); ir < nr; ++ir) {
0843     const auto &reg = event_.pfinputs[ir].region;
0844     for (const auto &p : event_.out[ir].pfcharged) {
0845       if (p.hwPt == 0 || !reg.isFiducial(p))
0846         continue;
0847       reco::Particle::PolarLorentzVector p4(
0848           p.floatPt(), reg.floatGlbEta(p.hwVtxEta()), reg.floatGlbPhi(p.hwVtxPhi()), 0.13f);
0849       l1t::PFCandidate::ParticleType type = l1t::PFCandidate::ChargedHadron;
0850       if (p.hwId.isMuon())
0851         type = l1t::PFCandidate::Muon;
0852       else if (p.hwId.isElectron())
0853         type = l1t::PFCandidate::Electron;
0854       ret->emplace_back(type, p.intCharge(), p4, 1, p.intPt(), p.intEta(), p.intPhi());
0855       ret->back().setZ0(p.floatZ0());
0856       ret->back().setDxy(p.floatDxy());
0857       ret->back().setHwZ0(p.hwZ0);
0858       ret->back().setHwDxy(p.hwDxy);
0859       ret->back().setHwTkQuality(p.hwTkQuality);
0860       setRefs_(ret->back(), p);
0861     }
0862     for (const auto &p : event_.out[ir].pfneutral) {
0863       if (p.hwPt == 0 || !reg.isFiducial(p))
0864         continue;
0865       reco::Particle::PolarLorentzVector p4(p.floatPt(), reg.floatGlbEtaOf(p), reg.floatGlbPhiOf(p), 0.13f);
0866       l1t::PFCandidate::ParticleType type =
0867           p.hwId.isPhoton() ? l1t::PFCandidate::Photon : l1t::PFCandidate::NeutralHadron;
0868       ret->emplace_back(type, 0, p4, 1, p.intPt(), p.intEta(), p.intPhi());
0869       ret->back().setHwEmID(p.hwEmID);
0870       setRefs_(ret->back(), p);
0871     }
0872   }
0873   return ret;
0874 }
0875 
0876 void L1TCorrelatorLayer1Producer::putPuppi(edm::Event &iEvent) const {
0877   auto refprod = iEvent.getRefBeforePut<l1t::PFCandidateCollection>("Puppi");
0878   auto coll = std::make_unique<l1t::PFCandidateCollection>();
0879   auto reg = std::make_unique<l1t::PFCandidateRegionalOutput>(refprod);
0880   std::vector<int> nobj;
0881   for (unsigned int ir = 0, nr = event_.pfinputs.size(); ir < nr; ++ir) {
0882     nobj.clear();
0883     for (const auto &p : event_.out[ir].puppi) {
0884       if (p.hwPt == 0)
0885         continue;
0886       // note: Puppi candidates are already in global coordinates & fiducial-only!
0887       l1t::PFCandidate::ParticleType type;
0888       float mass = 0.13f;
0889       if (p.hwId.charged()) {
0890         if (p.hwId.isMuon()) {
0891           type = l1t::PFCandidate::Muon;
0892           mass = 0.105;
0893         } else if (p.hwId.isElectron()) {
0894           type = l1t::PFCandidate::Electron;
0895           mass = 0.005;
0896         } else
0897           type = l1t::PFCandidate::ChargedHadron;
0898       } else {
0899         type = p.hwId.isPhoton() ? l1t::PFCandidate::Photon : l1t::PFCandidate::NeutralHadron;
0900         mass = p.hwId.isPhoton() ? 0.0 : 0.5;
0901       }
0902       reco::Particle::PolarLorentzVector p4(p.floatPt(), p.floatEta(), p.floatPhi(), mass);
0903       coll->emplace_back(type, p.intCharge(), p4, p.floatPuppiW(), p.intPt(), p.intEta(), p.intPhi());
0904       if (p.hwId.charged()) {
0905         coll->back().setZ0(p.floatZ0());
0906         coll->back().setDxy(p.floatDxy());
0907         coll->back().setHwZ0(p.hwZ0());
0908         coll->back().setHwDxy(p.hwDxy());
0909         coll->back().setHwTkQuality(p.hwTkQuality());
0910       } else {
0911         coll->back().setHwPuppiWeight(p.hwPuppiW());
0912         coll->back().setHwEmID(p.hwEmID());
0913       }
0914       coll->back().setEncodedPuppi64(p.pack().to_uint64());
0915       setRefs_(coll->back(), p);
0916       nobj.push_back(coll->size() - 1);
0917     }
0918     reg->addRegion(nobj, event_.pfinputs[ir].region.floatEtaCenter(), event_.pfinputs[ir].region.floatPhiCenter());
0919   }
0920   iEvent.put(std::move(coll), "Puppi");
0921   iEvent.put(std::move(reg), "PuppiRegional");
0922 }
0923 
0924 // NOTE: as a side effect we change the "sta_idx" of TkEle and TkEm objects to an index of the
0925 // vector of refs, for this reason this is not const. We could make this more explicit via arguments
0926 void L1TCorrelatorLayer1Producer::putEgStaObjects(edm::Event &iEvent,
0927                                                   const std::string &egLablel,
0928                                                   std::vector<edm::Ref<BXVector<l1t::EGamma>>> &egsta_refs) {
0929   auto egs = std::make_unique<BXVector<l1t::EGamma>>();
0930   edm::RefProd<BXVector<l1t::EGamma>> ref_egs = iEvent.getRefBeforePut<BXVector<l1t::EGamma>>(egLablel);
0931 
0932   edm::Ref<BXVector<l1t::EGamma>>::key_type idx = 0;
0933   // FIXME: in case more BXes are introduced shuld probably use egs->key(egs->end(bx));
0934 
0935   for (unsigned int ir = 0, nr = event_.pfinputs.size(); ir < nr; ++ir) {
0936     const auto &reg = event_.pfinputs[ir].region;
0937 
0938     std::vector<unsigned int> ref_pos(event_.out[ir].egsta.size());
0939 
0940     // EG standalone objects
0941     for (unsigned int ieg = 0, neg = event_.out[ir].egsta.size(); ieg < neg; ++ieg) {
0942       const auto &p = event_.out[ir].egsta[ieg];
0943       if (p.hwPt == 0 || !reg.isFiducial(p))
0944         continue;
0945       l1t::EGamma eg(
0946           reco::Candidate::PolarLorentzVector(p.floatPt(), reg.floatGlbEta(p.hwEta), reg.floatGlbPhi(p.hwPhi), 0.));
0947       eg.setHwQual(p.hwQual);
0948       egs->push_back(0, eg);
0949       egsta_refs.push_back(edm::Ref<BXVector<l1t::EGamma>>(ref_egs, idx++));
0950       ref_pos[ieg] = egsta_refs.size() - 1;
0951     }
0952 
0953     for (auto &egiso : event_.out[ir].egphoton) {
0954       if (egiso.hwPt == 0)
0955         continue;
0956       egiso.sta_idx = ref_pos[egiso.sta_idx];
0957     }
0958 
0959     for (auto &egele : event_.out[ir].egelectron) {
0960       if (egele.hwPt == 0)
0961         continue;
0962       egele.sta_idx = ref_pos[egele.sta_idx];
0963     }
0964   }
0965 
0966   iEvent.put(std::move(egs), egLablel);
0967 }
0968 
0969 void L1TCorrelatorLayer1Producer::putEgObjects(edm::Event &iEvent,
0970                                                const bool writeEgSta,
0971                                                const std::vector<edm::Ref<BXVector<l1t::EGamma>>> &egsta_refs,
0972                                                const std::string &tkEmLabel,
0973                                                const std::string &tkEmPerBoardLabel,
0974                                                const std::string &tkEleLabel,
0975                                                const std::string &tkElePerBoardLabel) const {
0976   auto tkems = std::make_unique<l1t::TkEmCollection>();
0977   auto tkemRefProd = iEvent.getRefBeforePut<l1t::TkEmCollection>(tkEmLabel);
0978   auto tkemPerBoard = std::make_unique<l1t::TkEmRegionalOutput>(tkemRefProd);
0979   auto tkeles = std::make_unique<l1t::TkElectronCollection>();
0980   auto tkeleRefProd = iEvent.getRefBeforePut<l1t::TkElectronCollection>(tkEleLabel);
0981   auto tkelePerBoard = std::make_unique<l1t::TkElectronRegionalOutput>(tkeleRefProd);
0982 
0983   // TkEG objects are written out after the per-board sorting.
0984   // The mapping to each board is saved into the regionalmap for further (stage-2 consumption)
0985   std::vector<int> nele_obj;
0986   std::vector<int> npho_obj;
0987 
0988   for (const auto &board : event_.board_out) {
0989     npho_obj.clear();
0990     for (const auto &egiso : board.egphoton) {
0991       if (egiso.hwPt == 0)
0992         continue;
0993 
0994       edm::Ref<BXVector<l1t::EGamma>> ref_egsta;
0995       if (writeEgSta) {
0996         ref_egsta = egsta_refs[egiso.sta_idx];
0997       } else {
0998         auto egptr = egiso.srcCluster->constituentsAndFractions()[0].first;
0999         ref_egsta =
1000             edm::Ref<BXVector<l1t::EGamma>>(egptr.id(), dynamic_cast<const l1t::EGamma *>(egptr.get()), egptr.key());
1001       }
1002 
1003       reco::Candidate::PolarLorentzVector mom(egiso.floatPt(), egiso.floatEta(), egiso.floatPhi(), 0.);
1004 
1005       l1t::TkEm tkem(reco::Candidate::LorentzVector(mom),
1006                      ref_egsta,
1007                      egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::TkIso),
1008                      egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::TkIsoPV));
1009       tkem.setHwQual(egiso.hwQual);
1010       tkem.setPFIsol(egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::PfIso));
1011       tkem.setPFIsolPV(egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::PfIsoPV));
1012       tkem.setEgBinaryWord(egiso.pack());
1013       tkems->push_back(tkem);
1014       npho_obj.push_back(tkems->size() - 1);
1015     }
1016     tkemPerBoard->addRegion(npho_obj, board.eta, board.phi);
1017 
1018     nele_obj.clear();
1019     for (const auto &egele : board.egelectron) {
1020       if (egele.hwPt == 0)
1021         continue;
1022 
1023       edm::Ref<BXVector<l1t::EGamma>> ref_egsta;
1024       if (writeEgSta) {
1025         ref_egsta = egsta_refs[egele.sta_idx];
1026       } else {
1027         auto egptr = egele.srcCluster->constituentsAndFractions()[0].first;
1028         ref_egsta =
1029             edm::Ref<BXVector<l1t::EGamma>>(egptr.id(), dynamic_cast<const l1t::EGamma *>(egptr.get()), egptr.key());
1030       }
1031 
1032       reco::Candidate::PolarLorentzVector mom(egele.floatPt(), egele.floatEta(), egele.floatPhi(), 0.);
1033 
1034       l1t::TkElectron tkele(reco::Candidate::LorentzVector(mom),
1035                             ref_egsta,
1036                             edm::refToPtr(egele.srcTrack->track()),
1037                             egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::TkIso));
1038       tkele.setHwQual(egele.hwQual);
1039       tkele.setPFIsol(egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::PfIso));
1040       tkele.setEgBinaryWord(egele.pack());
1041       tkeles->push_back(tkele);
1042       nele_obj.push_back(tkeles->size() - 1);
1043     }
1044     tkelePerBoard->addRegion(nele_obj, board.eta, board.phi);
1045   }
1046 
1047   iEvent.put(std::move(tkems), tkEmLabel);
1048   iEvent.put(std::move(tkemPerBoard), tkEmPerBoardLabel);
1049   iEvent.put(std::move(tkeles), tkEleLabel);
1050   iEvent.put(std::move(tkelePerBoard), tkElePerBoardLabel);
1051 }
1052 
1053 std::unique_ptr<std::vector<unsigned>> L1TCorrelatorLayer1Producer::vecSecInput(InputType t) const {
1054   auto v = std::make_unique<std::vector<unsigned>>();
1055   {
1056     switch (t) {
1057       case caloType:
1058         for (const auto &s : event_.decoded.hadcalo)
1059           v->push_back(s.size());
1060         break;
1061       case emcaloType:
1062         for (const auto &s : event_.decoded.emcalo)
1063           v->push_back(s.size());
1064         break;
1065       case trackType:
1066         for (const auto &s : event_.decoded.track)
1067           v->push_back(s.size());
1068         break;
1069       case l1muType:
1070         v->push_back(event_.decoded.muon.size());
1071         break;
1072     }
1073   }
1074   return v;
1075 }
1076 
1077 std::unique_ptr<std::vector<unsigned>> L1TCorrelatorLayer1Producer::vecRegInput(InputType t) const {
1078   auto v = std::make_unique<std::vector<unsigned>>();
1079   for (const auto &reg : event_.pfinputs) {
1080     switch (t) {
1081       case caloType:
1082         v->push_back(reg.hadcalo.size());
1083         break;
1084       case emcaloType:
1085         v->push_back(reg.emcalo.size());
1086         break;
1087       case trackType:
1088         v->push_back(reg.track.size());
1089         break;
1090       case l1muType:
1091         v->push_back(reg.muon.size());
1092         break;
1093     }
1094   }
1095   return v;
1096 }
1097 
1098 std::unique_ptr<std::vector<unsigned>> L1TCorrelatorLayer1Producer::vecOutput(OutputType i, bool usePuppi) const {
1099   auto v = std::make_unique<std::vector<unsigned>>();
1100   for (const auto &reg : event_.out) {
1101     v->push_back(reg.nObj(i, usePuppi));
1102   }
1103   return v;
1104 }
1105 std::pair<unsigned int, unsigned int> L1TCorrelatorLayer1Producer::totAndMax(
1106     const std::vector<unsigned> &perRegion) const {
1107   unsigned int ntot = 0, nmax = 0;
1108   for (unsigned ni : perRegion) {
1109     ntot += ni;
1110     nmax = std::max(nmax, ni);
1111   }
1112   return std::make_pair(ntot, nmax);
1113 }
1114 //define this as a plug-in
1115 #include "FWCore/Framework/interface/MakerMacros.h"
1116 DEFINE_FWK_MODULE(L1TCorrelatorLayer1Producer);
[ 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!