L1TTrackMatch/plugins/L1TrackJetProducer.cc

0001 // Original simulation Author:  Rishi Patel
0002 //
0003 // Rewritting/improvements:  George Karathanasis,
0004 //                           georgios.karathanasis@cern.ch, CU Boulder
0005 //                           Claire Savard (claire.savard@colorado.edu)
0006 //
0007 //         Created:  Wed, 01 Aug 2018 14:01:41 GMT
0008 //         Latest update: Nov 2023 (by CS)
0009 //
0010 // Track jets are clustered in a two-layer process, first by clustering in phi,
0011 // then by clustering in eta. The code proceeds as following: putting all tracks
0012 // 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 // system include files
0020 #include "DataFormats/Common/interface/Ref.h"
0021 #include "DataFormats/Common/interface/RefVector.h"
0022 #include "DataFormats/L1TCorrelator/interface/TkJet.h"
0023 #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h"
0024 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0025 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
0026 #include "DataFormats/L1Trigger/interface/VertexWord.h"
0027
0028 // user include files
0029 #include "FWCore/Framework/interface/Frameworkfwd.h"
0030 #include "FWCore/Framework/interface/stream/EDProducer.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/Framework/interface/EventSetup.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "FWCore/Utilities/interface/StreamID.h"
0036 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0037
0038 //own headers
0039 #include "L1TrackJetClustering.h"
0040
0041 using namespace std;
0042 using namespace edm;
0043 using namespace l1t;
0044 using namespace l1ttrackjet;
0045
0046 class L1TrackJetProducer : public stream::EDProducer<> {
0047 public:
0048   explicit L1TrackJetProducer(const ParameterSet &);
0049   ~L1TrackJetProducer() override = default;
0050   typedef TTTrack<Ref_Phase2TrackerDigi_> L1TTTrackType;
0051   typedef vector<L1TTTrackType> L1TTTrackCollectionType;
0052   typedef edm::RefVector<L1TTTrackCollectionType> L1TTTrackRefCollectionType;
0053   static void fillDescriptions(ConfigurationDescriptions &descriptions);
0054
0055 private:
0056   void produce(Event &, const EventSetup &) override;
0057
0058   // ----------member data ---------------------------
0059
0060   vector<Ptr<L1TTTrackType>> L1TrkPtrs_;
0061   const float trkZMax_;
0062   const float trkPtMax_;
0063   const float trkEtaMax_;
0064   const int lowpTJetMinTrackMultiplicity_;
0065   const float lowpTJetThreshold_;
0066   const int highpTJetMinTrackMultiplicity_;
0067   const float highpTJetThreshold_;
0068   const int zBins_;
0069   const int etaBins_;
0070   const int phiBins_;
0071   const double minTrkJetpT_;
0072   float zStep_;
0073   float etaStep_;
0074   float phiStep_;
0075   const bool displaced_;
0076   const float d0CutNStubs4_;
0077   const float d0CutNStubs5_;
0078   const int nDisplacedTracks_;
0079
0080   const EDGetTokenT<L1TTTrackRefCollectionType> trackToken_;
0081 };
0082
0083 L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig)
0084     : trkZMax_(iConfig.getParameter<double>("trk_zMax")),
0085       trkPtMax_(iConfig.getParameter<double>("trk_ptMax")),
0086       trkEtaMax_(iConfig.getParameter<double>("trk_etaMax")),
0087       lowpTJetMinTrackMultiplicity_(iConfig.getParameter<int>("lowpTJetMinTrackMultiplicity")),
0088       lowpTJetThreshold_(iConfig.getParameter<double>("lowpTJetThreshold")),
0089       highpTJetMinTrackMultiplicity_(iConfig.getParameter<int>("highpTJetMinTrackMultiplicity")),
0090       highpTJetThreshold_(iConfig.getParameter<double>("highpTJetThreshold")),
0091       zBins_(iConfig.getParameter<int>("zBins")),
0092       etaBins_(iConfig.getParameter<int>("etaBins")),
0093       phiBins_(iConfig.getParameter<int>("phiBins")),
0094       minTrkJetpT_(iConfig.getParameter<double>("minTrkJetpT")),
0095       displaced_(iConfig.getParameter<bool>("displaced")),
0096       d0CutNStubs4_(iConfig.getParameter<double>("d0_cutNStubs4")),
0097       d0CutNStubs5_(iConfig.getParameter<double>("d0_cutNStubs5")),
0098       nDisplacedTracks_(iConfig.getParameter<int>("nDisplacedTracks")),
0099       trackToken_(consumes<L1TTTrackRefCollectionType>(iConfig.getParameter<InputTag>("L1TrackInputTag"))) {
0100   zStep_ = 2.0 * trkZMax_ / (zBins_ + 1);  // added +1 in denom
0101   etaStep_ = 2.0 * trkEtaMax_ / etaBins_;  //etaStep is the width of an etabin
0102   phiStep_ = 2 * M_PI / phiBins_;          ////phiStep is the width of a phibin
0103
0104   if (displaced_)
0105     produces<TkJetCollection>("L1TrackJetsExtended");
0106   else
0107     produces<TkJetCollection>("L1TrackJets");
0108 }
0109
0110 void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) {
0111   unique_ptr<TkJetCollection> L1TrackJetProducer(new TkJetCollection);
0112
0113   // L1 tracks
0114   edm::Handle<L1TTTrackRefCollectionType> TTTrackHandle;
0115   iEvent.getByToken(trackToken_, TTTrackHandle);
0116
0117   L1TrkPtrs_.clear();
0118
0119   // track selection
0120   for (unsigned int this_l1track = 0; this_l1track < TTTrackHandle->size(); this_l1track++) {
0121     edm::Ptr<L1TTTrackType> trkPtr(TTTrackHandle, this_l1track);
0122     L1TrkPtrs_.push_back(trkPtr);
0123   }
0124
0125   // if no tracks pass selection return empty containers
0126   if (L1TrkPtrs_.empty()) {
0127     if (displaced_)
0128       iEvent.put(std::move(L1TrackJetProducer), "L1TrackJetsExtended");
0129     else
0130       iEvent.put(std::move(L1TrackJetProducer), "L1TrackJets");
0131     return;
0132   }
0133
0134   MaxZBin mzb;
0135   mzb.znum = -1;
0136   mzb.nclust = -1;
0137   mzb.ht = -1;
0138
0139   // create 2D grid of eta phi bins
0140   EtaPhiBin epbins_default[phiBins_][etaBins_];
0141   float phi = -1.0 * M_PI;
0142   for (int i = 0; i < phiBins_; ++i) {
0143     float eta = -1.0 * trkEtaMax_;
0144     for (int j = 0; j < etaBins_; ++j) {
0145       epbins_default[i][j].phi = (phi + (phi + phiStep_)) / 2.0;
0146       epbins_default[i][j].eta = (eta + (eta + etaStep_)) / 2.0;
0147       eta += etaStep_;
0148     }  // for each etabin
0149     phi += phiStep_;
0150   }  // for each phibin (finished creating bins)
0151
0152   // create z-bins (might be useful for displaced if we run w/o dz cut)
0153   std::vector<float> zmins, zmaxs;
0154   for (int zbin = 0; zbin < zBins_; zbin++) {
0155     zmins.push_back(-1.0 * trkZMax_ + zStep_ * zbin);
0156     zmaxs.push_back(-1.0 * trkZMax_ + zStep_ * zbin + zStep_ * 2);
0157   }
0158
0159   // create vectors of clusters
0160   std::vector<std::vector<EtaPhiBin>> L1clusters;
0161   L1clusters.reserve(phiBins_);
0162   std::vector<EtaPhiBin> L2clusters;
0163
0164   // default (empty) grid
0165   for (int i = 0; i < phiBins_; ++i) {
0166     for (int j = 0; j < etaBins_; ++j) {
0167       epbins_default[i][j].pTtot = 0;
0168       epbins_default[i][j].used = false;
0169       epbins_default[i][j].ntracks = 0;
0170       epbins_default[i][j].nxtracks = 0;
0171       epbins_default[i][j].trackidx.clear();
0172     }
0173   }
0174
0175   for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) {
0176     // initialize grid
0177     float zmin = zmins[zbin];
0178     float zmax = zmaxs[zbin];
0179     EtaPhiBin epbins[phiBins_][etaBins_];
0180     std::copy(&epbins_default[0][0], &epbins_default[0][0] + phiBins_ * etaBins_, &epbins[0][0]);
0181
0182     //clear cluster containers
0183     L1clusters.clear();
0184     L2clusters.clear();
0185
0186     // fill grid with tracks
0187     for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) {
0188       float trkZ = L1TrkPtrs_[k]->z0();
0189       if (zmax < trkZ)
0190         continue;
0191       if (zmin > trkZ)
0192         continue;
0193       if (zbin == 0 && zmin == trkZ)
0194         continue;
0195       float trkpt = L1TrkPtrs_[k]->momentum().perp();
0196       float trketa = L1TrkPtrs_[k]->momentum().eta();
0197       float trkphi = L1TrkPtrs_[k]->momentum().phi();
0198       float trkd0 = L1TrkPtrs_[k]->d0();
0199       int trknstubs = (int)L1TrkPtrs_[k]->getStubRefs().size();
0200       for (int i = 0; i < phiBins_; ++i) {
0201         for (int j = 0; j < etaBins_; ++j) {
0202           float eta_min = epbins[i][j].eta - etaStep_ / 2.0;  //eta min
0203           float eta_max = epbins[i][j].eta + etaStep_ / 2.0;  //eta max
0204           float phi_min = epbins[i][j].phi - phiStep_ / 2.0;  //phi min
0205           float phi_max = epbins[i][j].phi + phiStep_ / 2.0;  //phi max
0206
0207           if ((trketa < eta_min) || (trketa > eta_max) || (trkphi < phi_min) || (trkphi > phi_max))
0208             continue;
0209
0210           if (trkpt < trkPtMax_)
0211             epbins[i][j].pTtot += trkpt;
0212           else
0213             epbins[i][j].pTtot += trkPtMax_;
0214           if ((std::abs(trkd0) > d0CutNStubs5_ && trknstubs >= 5 && d0CutNStubs5_ >= 0) ||
0215               (trknstubs == 4 && std::abs(trkd0) > d0CutNStubs4_ && d0CutNStubs4_ >= 0))
0216             epbins[i][j].nxtracks += 1;
0217           epbins[i][j].trackidx.push_back(k);
0218           ++epbins[i][j].ntracks;
0219         }  // for each etabin
0220       }    // for each phibin
0221     }      //end loop over tracks
0222
0223     // cluster tracks in eta (first layer) using grid
0224     for (int phibin = 0; phibin < phiBins_; ++phibin) {
0225       L1clusters.push_back(L1_clustering<EtaPhiBin, float, float, float>(epbins[phibin], etaBins_, etaStep_));
0226     }
0227
0228     // second layer clustering in phi for using eta clusters
0229     L2clusters = L2_clustering<EtaPhiBin, float, float, float>(L1clusters, phiBins_, phiStep_, etaStep_);
0230
0231     // sum all cluster pTs in this zbin to find max
0232     float sum_pt = 0;
0233     for (unsigned int k = 0; k < L2clusters.size(); ++k) {
0234       if (L2clusters[k].pTtot > lowpTJetThreshold_ && L2clusters[k].ntracks < lowpTJetMinTrackMultiplicity_)
0235         continue;
0236       if (L2clusters[k].pTtot > highpTJetThreshold_ && L2clusters[k].ntracks < highpTJetMinTrackMultiplicity_)
0237         continue;
0238       if (L2clusters[k].pTtot > minTrkJetpT_)
0239         sum_pt += L2clusters[k].pTtot;
0240     }
0241
0242     if (sum_pt < mzb.ht)
0243       continue;
0244
0245     // if ht is larger than previous max, this is the new vertex zbin
0246     mzb.ht = sum_pt;
0247     mzb.znum = zbin;
0248     mzb.clusters = L2clusters;
0249     mzb.nclust = L2clusters.size();
0250     mzb.zbincenter = (zmin + zmax) / 2.0;
0251   }  //zbin loop
0252
0253   // output
0254   vector<Ptr<L1TTTrackType>> L1TrackAssocJet;
0255   for (unsigned int j = 0; j < mzb.clusters.size(); ++j) {
0256     float jetEta = mzb.clusters[j].eta;
0257     float jetPhi = mzb.clusters[j].phi;
0258     float jetPt = mzb.clusters[j].pTtot;
0259     float jetPx = jetPt * cos(jetPhi);
0260     float jetPy = jetPt * sin(jetPhi);
0261     float jetPz = jetPt * sinh(jetEta);
0262     float jetP = jetPt * cosh(jetEta);
0263     int totalDisptrk = mzb.clusters[j].nxtracks;
0264     bool isDispJet = (totalDisptrk >= nDisplacedTracks_);
0265
0266     math::XYZTLorentzVector jetP4(jetPx, jetPy, jetPz, jetP);
0267     L1TrackAssocJet.clear();
0268     for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++)
0269       L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]);
0270
0271     TkJet trkJet(jetP4, L1TrackAssocJet, mzb.zbincenter, mzb.clusters[j].ntracks, 0, totalDisptrk, 0, isDispJet);
0272
0273     L1TrackJetProducer->push_back(trkJet);
0274   }
0275
0276   std::sort(L1TrackJetProducer->begin(), L1TrackJetProducer->end(), [](auto &a, auto &b) { return a.pt() > b.pt(); });
0277   if (displaced_)
0278     iEvent.put(std::move(L1TrackJetProducer), "L1TrackJetsExtended");
0279   else
0280     iEvent.put(std::move(L1TrackJetProducer), "L1TrackJets");
0281 }
0282
0283 void L1TrackJetProducer::fillDescriptions(ConfigurationDescriptions &descriptions) {
0284   ParameterSetDescription desc;
0285   desc.add<edm::InputTag>(
0286       "L1TrackInputTag", edm::InputTag("l1tTrackVertexAssociationProducerForJets", "Level1TTTracksSelectedAssociated"));
0287   desc.add<double>("trk_zMax", 15.0);
0288   desc.add<double>("trk_ptMax", 200.0);
0289   desc.add<double>("trk_etaMax", 2.4);
0290   desc.add<double>("minTrkJetpT", -1.0);
0291   desc.add<int>("etaBins", 24);
0292   desc.add<int>("phiBins", 27);
0293   desc.add<int>("zBins", 1);
0294   desc.add<double>("d0_cutNStubs4", -1);
0295   desc.add<double>("d0_cutNStubs5", -1);
0296   desc.add<int>("lowpTJetMinTrackMultiplicity", 2);
0297   desc.add<double>("lowpTJetThreshold", 50.0);
0298   desc.add<int>("highpTJetMinTrackMultiplicity", 3);
0299   desc.add<double>("highpTJetThreshold", 100.0);
0300   desc.add<bool>("displaced", false);
0301   desc.add<int>("nDisplacedTracks", 2);
0302   descriptions.add("l1tTrackJets", desc);
0303 }
0304
0305 //define this as a plug-in
0306 DEFINE_FWK_MODULE(L1TrackJetProducer);