L1TTrackMatch/plugins/L1TrackJetEmulatorProducer.cc

0001 // Original Author: Samuel Edwin Leigh, Tyler Wu,
0002 //                  Rutgers, the State University of New Jersey
0003 //
0004 // Rewritting/Improvements:      George Karathanasis,
0005 //                          georgios.karathanasis@cern.ch, CU Boulder
0006 //                          Claire Savard (claire.savard@colorado.edu)
0007 //
0008 //         Created:  Wed, 01 Aug 2018 14:01:41 GMT
0009 //         Latest update: Nov 2023 (by CS)
0010 //
0011 // Track jets are clustered in a two-layer process, first by clustering in phi,
0012 // then by clustering in eta. The code proceeds as following: putting all tracks// in a grid of eta vs phi space, and then cluster them. Finally we merge the cl
0013 // usters when needed. The code is improved to use the same module between emula
0014 // tion and simulation was also improved, with bug fixes and being faster.
0015 // Introduction to object (p10-13):
0016 // https://indico.cern.ch/event/791517/contributions/3341650/attachments/1818736/2973771/TrackBasedAlgos_L1TMadrid_MacDonald.pdf
0017 // New and improved version: https://indico.cern.ch/event/1203796/contributions/5073056/attachments/2519806/4333006/trackjet_emu.pdf
0018
0019 // L1T include files
0020 #include "DataFormats/L1TCorrelator/interface/TkJet.h"
0021 #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h"
0022 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0023 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
0024 #include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h"
0025 #include "DataFormats/L1Trigger/interface/TkJetWord.h"
0026 #include "DataFormats/L1Trigger/interface/VertexWord.h"
0027
0028 // system include files
0029 #include "DataFormats/Common/interface/Handle.h"
0030 #include "FWCore/Framework/interface/Frameworkfwd.h"
0031 #include "FWCore/Framework/interface/stream/EDProducer.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/EventSetup.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "FWCore/Utilities/interface/StreamID.h"
0037 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0038 #include "DataFormats/Common/interface/Ref.h"
0039 #include "DataFormats/Common/interface/RefVector.h"
0040
0041 //own headers
0042 #include "L1TrackJetClustering.h"
0043
0044 //general
0045 #include <ap_int.h>
0046
0047 using namespace std;
0048 using namespace edm;
0049 using namespace l1t;
0050 using namespace l1ttrackjet;
0051
0052 class L1TrackJetEmulatorProducer : public stream::EDProducer<> {
0053 public:
0054   explicit L1TrackJetEmulatorProducer(const ParameterSet &);
0055   ~L1TrackJetEmulatorProducer() override = default;
0056   typedef TTTrack<Ref_Phase2TrackerDigi_> L1TTTrackType;
0057   typedef vector<L1TTTrackType> L1TTTrackCollectionType;
0058   typedef edm::RefVector<L1TTTrackCollectionType> L1TTTrackRefCollectionType;
0059   static void fillDescriptions(ConfigurationDescriptions &descriptions);
0060
0061 private:
0062   void produce(Event &, const EventSetup &) override;
0063
0064   // ----------member data ---------------------------
0065
0066   std::vector<edm::Ptr<L1TTTrackType>> L1TrkPtrs_;
0067   const float trkZMax_;
0068   const float trkPtMax_;
0069   const float trkEtaMax_;
0070   const int lowpTJetMinTrackMultiplicity_;
0071   const float lowpTJetThreshold_;
0072   const int highpTJetMinTrackMultiplicity_;
0073   const float highpTJetThreshold_;
0074   const int zBins_;
0075   const int etaBins_;
0076   const int phiBins_;
0077   const double minTrkJetpT_;
0078   const bool displaced_;
0079   const float d0CutNStubs4_;
0080   const float d0CutNStubs5_;
0081   const int nDisplacedTracks_;
0082
0083   float zStep_;
0084   glbeta_intern etaStep_;
0085   glbphi_intern phiStep_;
0086
0087   TTTrack_TrackWord trackword;
0088
0089   const EDGetTokenT<L1TTTrackRefCollectionType> trackToken_;
0090 };
0091
0092 //constructor
0093 L1TrackJetEmulatorProducer::L1TrackJetEmulatorProducer(const ParameterSet &iConfig)
0094     : trkZMax_(iConfig.getParameter<double>("trk_zMax")),
0095       trkPtMax_(iConfig.getParameter<double>("trk_ptMax")),
0096       trkEtaMax_(iConfig.getParameter<double>("trk_etaMax")),
0097       lowpTJetMinTrackMultiplicity_(iConfig.getParameter<int>("lowpTJetMinTrackMultiplicity")),
0098       lowpTJetThreshold_(iConfig.getParameter<double>("lowpTJetThreshold")),
0099       highpTJetMinTrackMultiplicity_(iConfig.getParameter<int>("highpTJetMinTrackMultiplicity")),
0100       highpTJetThreshold_(iConfig.getParameter<double>("highpTJetThreshold")),
0101       zBins_(iConfig.getParameter<int>("zBins")),
0102       etaBins_(iConfig.getParameter<int>("etaBins")),
0103       phiBins_(iConfig.getParameter<int>("phiBins")),
0104       minTrkJetpT_(iConfig.getParameter<double>("minTrkJetpT")),
0105       displaced_(iConfig.getParameter<bool>("displaced")),
0106       d0CutNStubs4_(iConfig.getParameter<double>("d0_cutNStubs4")),
0107       d0CutNStubs5_(iConfig.getParameter<double>("d0_cutNStubs5")),
0108       nDisplacedTracks_(iConfig.getParameter<int>("nDisplacedTracks")),
0109       trackToken_(consumes<L1TTTrackRefCollectionType>(iConfig.getParameter<InputTag>("L1TrackInputTag"))) {
0110   zStep_ = 2.0 * trkZMax_ / (zBins_ + 1);                 // added +1 in denom
0111   etaStep_ = glbeta_intern(2.0 * trkEtaMax_ / etaBins_);  //etaStep is the width of an etabin
0112   phiStep_ = DoubleToBit(2.0 * (M_PI) / phiBins_,
0113                          TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit,
0114                          TTTrack_TrackWord::stepPhi0);  ///phiStep is the width of a phibin
0115
0116   if (displaced_)
0117     produces<l1t::TkJetWordCollection>("L1TrackJetsExtended");
0118   else
0119     produces<l1t::TkJetWordCollection>("L1TrackJets");
0120 }
0121
0122 void L1TrackJetEmulatorProducer::produce(Event &iEvent, const EventSetup &iSetup) {
0123   unique_ptr<l1t::TkJetWordCollection> L1TrackJetContainer(new l1t::TkJetWordCollection);
0124
0125   // L1 tracks
0126   edm::Handle<L1TTTrackRefCollectionType> TTTrackHandle;
0127   iEvent.getByToken(trackToken_, TTTrackHandle);
0128
0129   L1TrkPtrs_.clear();
0130   // track selection
0131   for (unsigned int this_l1track = 0; this_l1track < TTTrackHandle->size(); this_l1track++) {
0132     edm::Ptr<L1TTTrackType> trkPtr(TTTrackHandle, this_l1track);
0133     L1TrkPtrs_.push_back(trkPtr);
0134   }
0135
0136   // if no tracks pass selection return empty containers
0137   if (L1TrkPtrs_.empty()) {
0138     if (displaced_)
0139       iEvent.put(std::move(L1TrackJetContainer), "L1TrackJetsExtended");
0140     else
0141       iEvent.put(std::move(L1TrackJetContainer), "L1TrackJets");
0142     return;
0143   }
0144
0145   TrackJetEmulationMaxZBin mzb;
0146   mzb.znum = 0;
0147   mzb.nclust = 0;
0148   mzb.ht = 0;
0149
0150   TrackJetEmulationEtaPhiBin epbins_default[phiBins_][etaBins_];  // create grid of phiBins
0151   glbphi_intern phi = DoubleToBit(
0152       -1.0 * M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0);
0153   for (int i = 0; i < phiBins_; ++i) {
0154     glbeta_intern eta = -1 * trkEtaMax_;
0155     for (int j = 0; j < etaBins_; ++j) {
0156       epbins_default[i][j].phi = (phi + (phi + phiStep_)) / 2;  // phi coord of bin center
0157       epbins_default[i][j].eta = (eta + (eta + etaStep_)) / 2;  // eta coord of bin center
0158       eta += etaStep_;
0159     }  // for each etabin
0160     phi += phiStep_;
0161   }  // for each phibin (finished creating epbins)
0162
0163   // bins for z if multibin option is selected
0164   std::vector<z0_intern> zmins, zmaxs;
0165   for (int zbin = 0; zbin < zBins_; zbin++) {
0166     zmins.push_back(DoubleToBit(
0167         -1.0 * trkZMax_ + zStep_ * zbin, TTTrack_TrackWord::TrackBitWidths::kZ0Size, TTTrack_TrackWord::stepZ0));
0168     zmaxs.push_back(DoubleToBit(-1.0 * trkZMax_ + zStep_ * zbin + 2 * zStep_,
0169                                 TTTrack_TrackWord::TrackBitWidths::kZ0Size,
0170                                 TTTrack_TrackWord::stepZ0));
0171   }
0172
0173   // create vectors that hold clusters
0174   std::vector<std::vector<TrackJetEmulationEtaPhiBin>> L1clusters;
0175   L1clusters.reserve(phiBins_);
0176   std::vector<TrackJetEmulationEtaPhiBin> L2clusters;
0177
0178   // default (empty) grid
0179   for (int i = 0; i < phiBins_; ++i) {
0180     for (int j = 0; j < etaBins_; ++j) {
0181       epbins_default[i][j].pTtot = 0;
0182       epbins_default[i][j].used = false;
0183       epbins_default[i][j].ntracks = 0;
0184       epbins_default[i][j].nxtracks = 0;
0185       epbins_default[i][j].trackidx.clear();
0186     }
0187   }
0188
0189   //Begin Firmware-style clustering
0190   // logic: loop over z bins find tracks in this bin and arrange them in a 2D eta-phi matrix
0191   for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) {
0192     // initialize matrices for every z bin
0193     z0_intern zmin = zmins[zbin];
0194     z0_intern zmax = zmaxs[zbin];
0195
0196     TrackJetEmulationEtaPhiBin epbins[phiBins_][etaBins_];
0197
0198     std::copy(&epbins_default[0][0], &epbins_default[0][0] + phiBins_ * etaBins_, &epbins[0][0]);
0199
0200     L1clusters.clear();
0201     L2clusters.clear();
0202     for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) {
0203       z0_intern trkZ = L1TrkPtrs_[k]->getZ0Word();
0204
0205       if (zmax < trkZ)
0206         continue;
0207       if (zmin > trkZ)
0208         continue;
0209       if (zbin == 0 && zmin == trkZ)
0210         continue;
0211
0212       // Pt
0213       ap_uint<TTTrack_TrackWord::TrackBitWidths::kRinvSize - 1> ptEmulationBits = L1TrkPtrs_[k]->getRinvWord();
0214       pt_intern trkpt;
0215       trkpt.V = ptEmulationBits.range();
0216
0217       // d0
0218       d0_intern abs_trkD0 = L1TrkPtrs_[k]->getD0Word();
0219
0220       // nstubs
0221       int trk_nstubs = (int)L1TrkPtrs_[k]->getStubRefs().size();
0222
0223       // Phi bin
0224       int i = phi_bin_firmwareStyle(L1TrkPtrs_[k]->phiSector(),
0225                                     L1TrkPtrs_[k]->getPhiWord());  //Function defined in L1TrackJetClustering.h
0226
0227       // Eta bin
0228       int j = eta_bin_firmwareStyle(L1TrkPtrs_[k]->getTanlWord());  //Function defined in L1TrackJetClustering.h
0229
0230       //This is a quick fix to eta going outside of scope - also including protection against phi going outside
0231       //of scope as well. The eta index, j, cannot be less than zero or greater than 23 (the number of eta bins
0232       //minus one). The phi index, i, cannot be less than zero or greater than 26 (the number of phi bins
0233       //minus one).
0234       if ((j < 0) || (j > (etaBins_ - 1)) || (i < 0) || (i > (phiBins_ - 1)))
0235         continue;
0236
0237       if (trkpt < pt_intern(trkPtMax_))
0238         epbins[i][j].pTtot += trkpt;
0239       else
0240         epbins[i][j].pTtot += pt_intern(trkPtMax_);
0241       if ((abs_trkD0 >
0242                DoubleToBit(d0CutNStubs5_, TTTrack_TrackWord::TrackBitWidths::kD0Size, TTTrack_TrackWord::stepD0) &&
0243            trk_nstubs >= 5 && d0CutNStubs5_ >= 0) ||
0244           (abs_trkD0 >
0245                DoubleToBit(d0CutNStubs4_, TTTrack_TrackWord::TrackBitWidths::kD0Size, TTTrack_TrackWord::stepD0) &&
0246            trk_nstubs == 4 && d0CutNStubs4_ >= 0))
0247         epbins[i][j].nxtracks += 1;
0248
0249       epbins[i][j].trackidx.push_back(k);
0250       ++epbins[i][j].ntracks;
0251     }
0252     //End Firmware style clustering
0253
0254     // first layer clustering - in eta using grid
0255     for (int phibin = 0; phibin < phiBins_; ++phibin) {
0256       L1clusters.push_back(L1_clustering<TrackJetEmulationEtaPhiBin, pt_intern, glbeta_intern, glbphi_intern>(
0257           epbins[phibin], etaBins_, etaStep_));
0258     }
0259
0260     //second layer clustering in phi for all eta clusters
0261     L2clusters = L2_clustering<TrackJetEmulationEtaPhiBin, pt_intern, glbeta_intern, glbphi_intern>(
0262         L1clusters, phiBins_, phiStep_, etaStep_);
0263
0264     // sum all cluster pTs in this zbin to find max
0265     pt_intern sum_pt = 0;
0266     for (unsigned int k = 0; k < L2clusters.size(); ++k) {
0267       if (L2clusters[k].pTtot > pt_intern(highpTJetThreshold_) && L2clusters[k].ntracks < lowpTJetMinTrackMultiplicity_)
0268         continue;
0269       if (L2clusters[k].pTtot > pt_intern(highpTJetThreshold_) &&
0270           L2clusters[k].ntracks < highpTJetMinTrackMultiplicity_)
0271         continue;
0272
0273       if (L2clusters[k].pTtot > pt_intern(minTrkJetpT_))
0274         sum_pt += L2clusters[k].pTtot;
0275     }
0276     if (sum_pt < mzb.ht)
0277       continue;
0278     // if ht is larger than previous max, this is the new vertex zbin
0279     mzb.ht = sum_pt;
0280     mzb.znum = zbin;
0281     mzb.clusters = L2clusters;
0282     mzb.nclust = L2clusters.size();
0283     mzb.zbincenter = (zmin + zmax) / 2.0;
0284   }  //zbin loop
0285
0286   vector<edm::Ptr<L1TTTrackType>> L1TrackAssocJet;
0287   for (unsigned int j = 0; j < mzb.clusters.size(); ++j) {
0288     l1t::TkJetWord::glbeta_t jetEta = DoubleToBit(double(mzb.clusters[j].eta),
0289                                                   TkJetWord::TkJetBitWidths::kGlbEtaSize,
0290                                                   TkJetWord::MAX_ETA / (1 << TkJetWord::TkJetBitWidths::kGlbEtaSize));
0291     l1t::TkJetWord::glbphi_t jetPhi = DoubleToBit(
0292         BitToDouble(mzb.clusters[j].phi, TTTrack_TrackWord::TrackBitWidths::kPhiSize + 4, TTTrack_TrackWord::stepPhi0),
0293         TkJetWord::TkJetBitWidths::kGlbPhiSize,
0294         (2. * std::abs(M_PI)) / (1 << TkJetWord::TkJetBitWidths::kGlbPhiSize));
0295     l1t::TkJetWord::z0_t jetZ0 = 0;
0296     l1t::TkJetWord::pt_t jetPt = mzb.clusters[j].pTtot;
0297     l1t::TkJetWord::nt_t total_ntracks = mzb.clusters[j].ntracks;
0298     l1t::TkJetWord::nx_t total_disptracks = mzb.clusters[j].nxtracks;
0299     l1t::TkJetWord::dispflag_t dispflag = 0;
0300     l1t::TkJetWord::tkjetunassigned_t unassigned = 0;
0301
0302     if (total_disptracks >= nDisplacedTracks_)
0303       dispflag = 1;
0304     L1TrackAssocJet.clear();
0305     for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++)
0306       L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]);
0307
0308     l1t::TkJetWord trkJet(jetPt, jetEta, jetPhi, jetZ0, total_ntracks, total_disptracks, dispflag, unassigned);
0309
0310     L1TrackJetContainer->push_back(trkJet);
0311   }
0312
0313   std::sort(L1TrackJetContainer->begin(), L1TrackJetContainer->end(), [](auto &a, auto &b) { return a.pt() > b.pt(); });
0314   if (displaced_)
0315     iEvent.put(std::move(L1TrackJetContainer), "L1TrackJetsExtended");
0316   else
0317     iEvent.put(std::move(L1TrackJetContainer), "L1TrackJets");
0318 }
0319
0320 void L1TrackJetEmulatorProducer::fillDescriptions(ConfigurationDescriptions &descriptions) {
0321   //The following says we do not know what parameters are allowed so do no validation
0322   // Please change this to state exactly what you do use, even if it is no parameters
0323   ParameterSetDescription desc;
0324   desc.add<edm::InputTag>("L1TrackInputTag", edm::InputTag("l1tTTTracksFromTrackletEmulation", "Level1TTTracks"));
0325   desc.add<double>("trk_zMax", 15.0);
0326   desc.add<double>("trk_ptMax", 200.0);
0327   desc.add<double>("trk_etaMax", 2.4);
0328   desc.add<double>("minTrkJetpT", -1.0);
0329   desc.add<int>("etaBins", 24);
0330   desc.add<int>("phiBins", 27);
0331   desc.add<int>("zBins", 1);
0332   desc.add<double>("d0_cutNStubs4", -1);
0333   desc.add<double>("d0_cutNStubs5", -1);
0334   desc.add<int>("lowpTJetMinTrackMultiplicity", 2);
0335   desc.add<double>("lowpTJetThreshold", 50.0);
0336   desc.add<int>("highpTJetMinTrackMultiplicity", 3);
0337   desc.add<double>("highpTJetThreshold", 100.0);
0338   desc.add<bool>("displaced", false);
0339   desc.add<int>("nDisplacedTracks", 2);
0340   descriptions.add("l1tTrackJetsEmulator", desc);
0341 }
0342
0343 //define this as a plug-in
0344 DEFINE_FWK_MODULE(L1TrackJetEmulatorProducer);