Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 // Original Author:  Rishi Patel
0002 // Modifications:    George Karathanasis, georgios.karathanasis@cern.ch, CU Boulder
0003 //         Created:  Wed, 01 Aug 2018 14:01:41 GMT
0004 //         Latest update: Nov 2021 (by GK)
0005 //
0006 // Track jets are clustered in a two-layer process, first by clustering in phi,
0007 // then by clustering in eta
0008 // Introduction to object (p10-13):
0009 // https://indico.cern.ch/event/791517/contributions/3341650/attachments/1818736/2973771/TrackBasedAlgos_L1TMadrid_MacDonald.pdf
0010 
0011 // system include files
0012 
0013 #include "DataFormats/Common/interface/Ref.h"
0014 #include "DataFormats/L1TCorrelator/interface/TkJet.h"
0015 #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h"
0016 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0017 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
0018 // user include files
0019 #include "FWCore/Framework/interface/Frameworkfwd.h"
0020 #include "FWCore/Framework/interface/stream/EDProducer.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Utilities/interface/StreamID.h"
0026 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0027 #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
0028 #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h"
0029 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0030 
0031 #include "DataFormats/L1Trigger/interface/Vertex.h"
0032 
0033 #include "L1Trigger/L1TTrackMatch/interface/L1Clustering.h"
0034 
0035 #include "TH1D.h"
0036 #include "TH2D.h"
0037 #include <TMath.h>
0038 #include <cmath>
0039 #include <cstdlib>
0040 #include <fstream>
0041 #include <iostream>
0042 #include <memory>
0043 #include <string>
0044 
0045 using namespace std;
0046 using namespace edm;
0047 using namespace l1t;
0048 
0049 class L1TrackJetProducer : public stream::EDProducer<> {
0050 public:
0051   explicit L1TrackJetProducer(const ParameterSet &);
0052   ~L1TrackJetProducer() override;
0053   typedef TTTrack<Ref_Phase2TrackerDigi_> L1TTTrackType;
0054   typedef vector<L1TTTrackType> L1TTTrackCollectionType;
0055 
0056   static void fillDescriptions(ConfigurationDescriptions &descriptions);
0057   bool trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2);
0058 
0059 private:
0060   void beginStream(StreamID) override;
0061   void produce(Event &, const EventSetup &) override;
0062   void endStream() override;
0063 
0064   // ----------member data ---------------------------
0065 
0066   vector<Ptr<L1TTTrackType>> L1TrkPtrs_;
0067   vector<int> tdtrk_;
0068   const float trkZMax_;
0069   const float trkPtMax_;
0070   const float trkPtMin_;
0071   const float trkEtaMax_;
0072   const float nStubs4PromptChi2_;
0073   const float nStubs5PromptChi2_;
0074   const float nStubs4PromptBend_;
0075   const float nStubs5PromptBend_;
0076   const int trkNPSStubMin_;
0077   const int lowpTJetMinTrackMultiplicity_;
0078   const float lowpTJetThreshold_;
0079   const int highpTJetMinTrackMultiplicity_;
0080   const float highpTJetThreshold_;
0081   const int zBins_;
0082   const int etaBins_;
0083   const int phiBins_;
0084   const double minTrkJetpT_;
0085   float zStep_;
0086   float etaStep_;
0087   float phiStep_;
0088   const bool displaced_;
0089   const float d0CutNStubs4_;
0090   const float d0CutNStubs5_;
0091   const float nStubs4DisplacedChi2_;
0092   const float nStubs5DisplacedChi2_;
0093   const float nStubs4DisplacedBend_;
0094   const float nStubs5DisplacedBend_;
0095   const int nDisplacedTracks_;
0096   const float dzPVTrk_;
0097 
0098   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0099   const EDGetTokenT<vector<TTTrack<Ref_Phase2TrackerDigi_>>> trackToken_;
0100   const edm::EDGetTokenT<std::vector<l1t::Vertex>> PVtxToken_;
0101 };
0102 
0103 L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig)
0104     : trkZMax_((float)iConfig.getParameter<double>("trk_zMax")),
0105       trkPtMax_((float)iConfig.getParameter<double>("trk_ptMax")),
0106       trkPtMin_((float)iConfig.getParameter<double>("trk_ptMin")),
0107       trkEtaMax_((float)iConfig.getParameter<double>("trk_etaMax")),
0108       nStubs4PromptChi2_((float)iConfig.getParameter<double>("nStubs4PromptChi2")),
0109       nStubs5PromptChi2_((float)iConfig.getParameter<double>("nStubs5PromptChi2")),
0110       nStubs4PromptBend_((float)iConfig.getParameter<double>("nStubs4PromptBend")),
0111       nStubs5PromptBend_((float)iConfig.getParameter<double>("nStubs5PromptBend")),
0112       trkNPSStubMin_((int)iConfig.getParameter<int>("trk_nPSStubMin")),
0113       lowpTJetMinTrackMultiplicity_((int)iConfig.getParameter<int>("lowpTJetMinTrackMultiplicity")),
0114       lowpTJetThreshold_((float)iConfig.getParameter<double>("lowpTJetThreshold")),
0115       highpTJetMinTrackMultiplicity_((int)iConfig.getParameter<int>("highpTJetMinTrackMultiplicity")),
0116       highpTJetThreshold_((float)iConfig.getParameter<double>("highpTJetThreshold")),
0117       zBins_((int)iConfig.getParameter<int>("zBins")),
0118       etaBins_((int)iConfig.getParameter<int>("etaBins")),
0119       phiBins_((int)iConfig.getParameter<int>("phiBins")),
0120       minTrkJetpT_(iConfig.getParameter<double>("minTrkJetpT")),
0121       displaced_(iConfig.getParameter<bool>("displaced")),
0122       d0CutNStubs4_((float)iConfig.getParameter<double>("d0_cutNStubs4")),
0123       d0CutNStubs5_((float)iConfig.getParameter<double>("d0_cutNStubs5")),
0124       nStubs4DisplacedChi2_((float)iConfig.getParameter<double>("nStubs4DisplacedChi2")),
0125       nStubs5DisplacedChi2_((float)iConfig.getParameter<double>("nStubs5DisplacedChi2")),
0126       nStubs4DisplacedBend_((float)iConfig.getParameter<double>("nStubs4DisplacedBend")),
0127       nStubs5DisplacedBend_((float)iConfig.getParameter<double>("nStubs5DisplacedBend")),
0128       nDisplacedTracks_((int)iConfig.getParameter<int>("nDisplacedTracks")),
0129       dzPVTrk_((float)iConfig.getParameter<double>("MaxDzTrackPV")),
0130       tTopoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>(edm::ESInputTag("", ""))),
0131       trackToken_(consumes<vector<TTTrack<Ref_Phase2TrackerDigi_>>>(iConfig.getParameter<InputTag>("L1TrackInputTag"))),
0132       PVtxToken_(consumes<vector<l1t::Vertex>>(iConfig.getParameter<InputTag>("L1PVertexCollection"))) {
0133   zStep_ = 2.0 * trkZMax_ / (zBins_ + 1);  // added +1 in denom
0134   etaStep_ = 2.0 * trkEtaMax_ / etaBins_;  //etaStep is the width of an etabin
0135   phiStep_ = 2 * M_PI / phiBins_;          ////phiStep is the width of a phibin
0136 
0137   if (displaced_)
0138     produces<TkJetCollection>("L1TrackJetsExtended");
0139   else
0140     produces<TkJetCollection>("L1TrackJets");
0141 }
0142 
0143 L1TrackJetProducer::~L1TrackJetProducer() {}
0144 
0145 void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) {
0146   unique_ptr<TkJetCollection> L1L1TrackJetProducer(new TkJetCollection);
0147 
0148   // Read inputs
0149   const TrackerTopology &tTopo = iSetup.getData(tTopoToken_);
0150 
0151   edm::Handle<vector<TTTrack<Ref_Phase2TrackerDigi_>>> TTTrackHandle;
0152   iEvent.getByToken(trackToken_, TTTrackHandle);
0153 
0154   edm::Handle<std::vector<l1t::Vertex>> PVtx;
0155   iEvent.getByToken(PVtxToken_, PVtx);
0156   float PVz = (PVtx->at(0)).z0();
0157 
0158   L1TrkPtrs_.clear();
0159   tdtrk_.clear();
0160 
0161   for (unsigned int this_l1track = 0; this_l1track < TTTrackHandle->size(); this_l1track++) {
0162     edm::Ptr<L1TTTrackType> trkPtr(TTTrackHandle, this_l1track);
0163     float trk_pt = trkPtr->momentum().perp();
0164     int trk_nstubs = (int)trkPtr->getStubRefs().size();
0165     float trk_chi2dof = trkPtr->chi2Red();
0166     float trk_d0 = trkPtr->d0();
0167     float trk_bendchi2 = trkPtr->stubPtConsistency();
0168     float trk_z0 = trkPtr->z0();
0169 
0170     int trk_nPS = 0;
0171     for (int istub = 0; istub < trk_nstubs; istub++) {  // loop over the stubs
0172       DetId detId(trkPtr->getStubRefs().at(istub)->getDetId());
0173       if (detId.det() == DetId::Detector::Tracker) {
0174         if ((detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3) ||
0175             (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9))
0176           trk_nPS++;
0177       }
0178     }
0179     // select tracks
0180     if (trk_nPS < trkNPSStubMin_)
0181       continue;
0182     if (!trackQualityCuts(trk_nstubs, trk_chi2dof, trk_bendchi2))
0183       continue;
0184     if (std::abs(PVz - trk_z0) > dzPVTrk_ && dzPVTrk_ > 0)
0185       continue;
0186     if (std::abs(trk_z0) > trkZMax_)
0187       continue;
0188     if (std::abs(trkPtr->momentum().eta()) > trkEtaMax_)
0189       continue;
0190     if (trk_pt < trkPtMin_)
0191       continue;
0192     L1TrkPtrs_.push_back(trkPtr);
0193 
0194     if ((std::abs(trk_d0) > d0CutNStubs5_ && trk_nstubs >= 5 && d0CutNStubs5_ >= 0) ||
0195         (trk_nstubs == 4 && std::abs(trk_d0) > d0CutNStubs4_ && d0CutNStubs4_ >= 0))
0196       tdtrk_.push_back(1);  //displaced track
0197     else
0198       tdtrk_.push_back(0);  // not displaced track
0199   }
0200 
0201   if (!L1TrkPtrs_.empty()) {
0202     MaxZBin mzb;
0203     mzb.znum = -1;
0204     mzb.nclust = -1;
0205     mzb.ht = -1;
0206     EtaPhiBin epbins_default[phiBins_][etaBins_];  // create grid of phiBins
0207 
0208     float phi = -1.0 * M_PI;
0209     for (int i = 0; i < phiBins_; ++i) {
0210       float eta = -1.0 * trkEtaMax_;
0211       for (int j = 0; j < etaBins_; ++j) {
0212         epbins_default[i][j].phi = (phi + (phi + phiStep_)) / 2.0;  // phimin=phi; phimax= phimin+step
0213         epbins_default[i][j].eta = (eta + (eta + etaStep_)) / 2.0;  // phimin=phi; phimax phimin+step
0214         eta += etaStep_;
0215       }  // for each etabin
0216       phi += phiStep_;
0217     }  // for each phibin (finished creating epbins)
0218 
0219     std::vector<float> zmins, zmaxs;
0220     for (int zbin = 0; zbin < zBins_; zbin++) {
0221       zmins.push_back(-1.0 * trkZMax_ + zStep_ * zbin);
0222       zmaxs.push_back(-1.0 * trkZMax_ + zStep_ * zbin + zStep_ * 2);
0223     }
0224 
0225     // create vectors that hold data
0226     std::vector<std::vector<EtaPhiBin>> L1clusters;
0227     L1clusters.reserve(phiBins_);
0228     std::vector<EtaPhiBin> L2clusters;
0229 
0230     for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) {
0231       // initialize matrices
0232       float zmin = zmins[zbin];
0233       float zmax = zmaxs[zbin];
0234       EtaPhiBin epbins[phiBins_][etaBins_];
0235       std::copy(&epbins_default[0][0], &epbins_default[0][0] + phiBins_ * etaBins_, &epbins[0][0]);
0236 
0237       //clear containers
0238       L1clusters.clear();
0239       L2clusters.clear();
0240 
0241       // fill grid
0242       for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) {
0243         float trkZ = L1TrkPtrs_[k]->z0();
0244         if (zmax < trkZ)
0245           continue;
0246         if (zmin > trkZ)
0247           continue;
0248         if (zbin == 0 && zmin == trkZ)
0249           continue;
0250         float trkpt = L1TrkPtrs_[k]->momentum().perp();
0251         float trketa = L1TrkPtrs_[k]->momentum().eta();
0252         float trkphi = L1TrkPtrs_[k]->momentum().phi();
0253         for (int i = 0; i < phiBins_; ++i) {
0254           for (int j = 0; j < etaBins_; ++j) {
0255             float eta_min = epbins[i][j].eta - etaStep_ / 2.0;  //eta min
0256             float eta_max = epbins[i][j].eta + etaStep_ / 2.0;  //eta max
0257             float phi_min = epbins[i][j].phi - phiStep_ / 2.0;  //phi min
0258             float phi_max = epbins[i][j].phi + phiStep_ / 2.0;  //phi max
0259 
0260             if ((trketa < eta_min) || (trketa > eta_max) || (trkphi < phi_min) || (trkphi > phi_max))
0261               continue;
0262 
0263             if (trkpt < trkPtMax_)
0264               epbins[i][j].pTtot += trkpt;
0265             else
0266               epbins[i][j].pTtot += trkPtMax_;
0267             epbins[i][j].numtdtrks += tdtrk_[k];
0268             epbins[i][j].trackidx.push_back(k);
0269             ++epbins[i][j].numtracks;
0270           }  // for each phibin
0271         }    // for each phibin
0272       }      //end loop over tracks
0273 
0274       // first layer clustering - in eta for all phi bins
0275       for (int phibin = 0; phibin < phiBins_; ++phibin) {
0276         L1clusters.push_back(L1_clustering(epbins[phibin], etaBins_, etaStep_));
0277       }
0278 
0279       //second layer clustering in phi for all eta clusters
0280       L2clusters = L2_clustering(L1clusters, phiBins_, phiStep_, etaStep_);
0281 
0282       // sum all cluster pTs in this zbin to find max
0283       float sum_pt = 0;
0284       for (unsigned int k = 0; k < L2clusters.size(); ++k) {
0285         if (L2clusters[k].pTtot > lowpTJetThreshold_ && L2clusters[k].numtracks < lowpTJetMinTrackMultiplicity_)
0286           continue;
0287         if (L2clusters[k].pTtot > highpTJetThreshold_ && L2clusters[k].numtracks < highpTJetMinTrackMultiplicity_)
0288           continue;
0289         if (L2clusters[k].pTtot > minTrkJetpT_)
0290           sum_pt += L2clusters[k].pTtot;
0291       }
0292 
0293       if (sum_pt < mzb.ht)
0294         continue;
0295       // if ht is larger than previous max, this is the new vertex zbin
0296       mzb.ht = sum_pt;
0297       mzb.znum = zbin;
0298       mzb.clusters = L2clusters;
0299       mzb.nclust = L2clusters.size();
0300       mzb.zbincenter = (zmin + zmax) / 2.0;
0301     }  //zbin loop
0302 
0303     vector<Ptr<L1TTTrackType>> L1TrackAssocJet;
0304     for (unsigned int j = 0; j < mzb.clusters.size(); ++j) {
0305       float jetEta = mzb.clusters[j].eta;
0306       float jetPhi = mzb.clusters[j].phi;
0307       float jetPt = mzb.clusters[j].pTtot;
0308       float jetPx = jetPt * cos(jetPhi);
0309       float jetPy = jetPt * sin(jetPhi);
0310       float jetPz = jetPt * sinh(jetEta);
0311       float jetP = jetPt * cosh(jetEta);
0312       int totalDisptrk = mzb.clusters[j].numtdtrks;
0313       bool isDispJet = (totalDisptrk > nDisplacedTracks_ || totalDisptrk == nDisplacedTracks_);
0314 
0315       math::XYZTLorentzVector jetP4(jetPx, jetPy, jetPz, jetP);
0316       L1TrackAssocJet.clear();
0317       for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++)
0318         L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]);
0319 
0320       TkJet trkJet(jetP4, L1TrackAssocJet, mzb.zbincenter, mzb.clusters[j].numtracks, 0, totalDisptrk, 0, isDispJet);
0321 
0322       if (!L1TrackAssocJet.empty())
0323         L1L1TrackJetProducer->push_back(trkJet);
0324     }
0325 
0326     if (displaced_)
0327       iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended");
0328     else
0329       iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets");
0330   }
0331 }
0332 
0333 void L1TrackJetProducer::beginStream(StreamID) {}
0334 
0335 void L1TrackJetProducer::endStream() {}
0336 
0337 bool L1TrackJetProducer::trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2) {
0338   bool PassQuality = false;
0339   if (!displaced_) {
0340     if (trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ &&
0341         trk_chi2 < nStubs4PromptChi2_)  // 4 stubs are the lowest track quality and have different cuts
0342       PassQuality = true;
0343     if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ &&
0344         trk_chi2 < nStubs5PromptChi2_)  // above 4 stubs diffent selection imposed (genrally looser)
0345       PassQuality = true;
0346   } else {
0347     if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ &&
0348         trk_chi2 < nStubs4DisplacedChi2_)  // 4 stubs are the lowest track quality and have different cuts
0349       PassQuality = true;
0350     if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ &&
0351         trk_chi2 < nStubs5DisplacedChi2_)  // above 4 stubs diffent selection imposed (genrally looser)
0352       PassQuality = true;
0353   }
0354   return PassQuality;
0355 }
0356 
0357 void L1TrackJetProducer::fillDescriptions(ConfigurationDescriptions &descriptions) {
0358   // l1tTrackJets
0359   edm::ParameterSetDescription desc;
0360   desc.add<edm::InputTag>("L1TrackInputTag", edm::InputTag("l1tTTTracksFromTrackletEmulation", "Level1TTTracks"));
0361   desc.add<edm::InputTag>("L1PVertexCollection", edm::InputTag("l1tVertexProducer", "l1vertices"));
0362   desc.add<double>("MaxDzTrackPV", 1.0);
0363   desc.add<double>("trk_zMax", 15.0);
0364   desc.add<double>("trk_ptMax", 200.0);
0365   desc.add<double>("trk_ptMin", 3.0);
0366   desc.add<double>("trk_etaMax", 2.4);
0367   desc.add<double>("nStubs4PromptChi2", 5.0);
0368   desc.add<double>("nStubs4PromptBend", 1.7);
0369   desc.add<double>("nStubs5PromptChi2", 2.75);
0370   desc.add<double>("nStubs5PromptBend", 3.5);
0371   desc.add<int>("trk_nPSStubMin", -1);
0372   desc.add<double>("minTrkJetpT", -1.0);
0373   desc.add<int>("etaBins", 24);
0374   desc.add<int>("phiBins", 27);
0375   desc.add<int>("zBins", 1);
0376   desc.add<double>("d0_cutNStubs4", -1);
0377   desc.add<double>("d0_cutNStubs5", -1);
0378   desc.add<int>("lowpTJetMinTrackMultiplicity", 2);
0379   desc.add<double>("lowpTJetThreshold", 50.0);
0380   desc.add<int>("highpTJetMinTrackMultiplicity", 3);
0381   desc.add<double>("highpTJetThreshold", 100.0);
0382   desc.add<bool>("displaced", false);
0383   desc.add<double>("nStubs4DisplacedChi2", 5.0);
0384   desc.add<double>("nStubs4DisplacedBend", 1.7);
0385   desc.add<double>("nStubs5DisplacedChi2", 2.75);
0386   desc.add<double>("nStubs5DisplacedBend", 3.5);
0387   desc.add<int>("nDisplacedTracks", 2);
0388   descriptions.add("l1tTrackJets", desc);
0389 }
0390 
0391 //define this as a plug-in
0392 DEFINE_FWK_MODULE(L1TrackJetProducer);