From 993794f306291cce34d81f1013107710fdb56bf6 Mon Sep 17 00:00:00 2001 From: Georgios Karathanasis Date: Fri, 29 Apr 2022 11:49:24 +0200 Subject: [PATCH 1/8] update trackjet simulation rebase (cherry picked from commit 834a2df7dcc2994cb808dd2018d971faafbc2748) --- .../plugins/L1TrackJetProducer.cc | 738 ++++++------------ .../plugins/L1TrackJetProducer.h | 250 ++++++ .../L1TTrackMatch/python/l1tTrackJets_cfi.py | 69 +- 3 files changed, 524 insertions(+), 533 deletions(-) create mode 100644 L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc index 46f0208e2852b..4d316f0af8196 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc @@ -1,6 +1,7 @@ // Original Author: Rishi Patel // Modifications: George Karathanasis, georgios.karathanasis@cern.ch, CU Boulder // Created: Wed, 01 Aug 2018 14:01:41 GMT +// Latest update: Nov 2021 (by GK) // // Track jets are clustered in a two-layer process, first by clustering in phi, // then by clustering in eta @@ -18,6 +19,7 @@ #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -27,8 +29,10 @@ #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "DataFormats/L1TCorrelator/interface/TkPrimaryVertex.h" #include "DataFormats/L1Trigger/interface/Vertex.h" -#include "L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h" + +#include "L1TrackJetProducer.h" #include "TH1D.h" #include "TH2D.h" #include @@ -42,6 +46,8 @@ using namespace std; using namespace edm; using namespace l1t; + + class L1TrackJetProducer : public stream::EDProducer<> { public: explicit L1TrackJetProducer(const ParameterSet &); @@ -51,8 +57,7 @@ class L1TrackJetProducer : public stream::EDProducer<> { static void fillDescriptions(ConfigurationDescriptions &descriptions); bool trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2); - void L2_cluster(vector > L1TrkPtrs_, vector tdtrk_, MaxZBin &mzb); - virtual EtaPhiBin *L1_cluster(EtaPhiBin *phislice); + private: void beginStream(StreamID) override; @@ -61,73 +66,69 @@ class L1TrackJetProducer : public stream::EDProducer<> { // ----------member data --------------------------- + const EDGetTokenT > > trackToken_; + const edm::EDGetTokenT > PVtxToken_; vector > L1TrkPtrs_; - vector zBinCount_; vector tdtrk_; - const float trkZMax_; - const float trkPtMax_; - const float trkPtMin_; - const float trkEtaMax_; - const float trkChi2dofMax_; - const float trkBendChi2Max_; - const int trkNPSStubMin_; - const double minTrkJetpT_; - const double minJetEtLowPt_; - const double minJetEtHighPt_; + float trkZMax_; + float trkPtMax_; + float trkPtMin_; + float trkEtaMax_; + float nStubs4PromptChi2_; + float nStubs5PromptChi2_; + float nStubs4PromptBend_; + float nStubs5PromptBend_; + int trkNPSStubMin_; + int lowpTJetMinTrackMultiplicity_; + int highpTJetMinTrackMultiplicity_; + int zBins_; int etaBins_; int phiBins_; - int zBins_; + double minTrkJetpT_; + float zStep_; + float etaStep_; + float phiStep_; + bool displaced_; float d0CutNStubs4_; float d0CutNStubs5_; - const int lowpTJetMinTrackMultiplicity_; - const int highpTJetMinTrackMultiplicity_; - bool displaced_; float nStubs4DisplacedChi2_; float nStubs5DisplacedChi2_; float nStubs4DisplacedBend_; float nStubs5DisplacedBend_; int nDisplacedTracks_; float dzPVTrk_; - - const EDGetTokenT > > trackToken_; - const edm::EDGetTokenT > PVtxToken_; - edm::ESGetToken tTopoToken_; - - float zStep_; - float etaStep_; - float phiStep_; }; L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) - : trkZMax_((float)iConfig.getParameter("trk_zMax")), - trkPtMax_((float)iConfig.getParameter("trk_ptMax")), - trkPtMin_((float)iConfig.getParameter("trk_ptMin")), - trkEtaMax_((float)iConfig.getParameter("trk_etaMax")), - trkChi2dofMax_((float)iConfig.getParameter("trk_chi2dofMax")), - trkBendChi2Max_((float)iConfig.getParameter("trk_bendChi2Max")), - trkNPSStubMin_((int)iConfig.getParameter("trk_nPSStubMin")), - minTrkJetpT_(iConfig.getParameter("minTrkJetpT")), - minJetEtLowPt_(iConfig.getParameter("minJetEtLowPt")), - minJetEtHighPt_(iConfig.getParameter("minJetEtHighPt")), - etaBins_((int)iConfig.getParameter("etaBins")), - phiBins_((int)iConfig.getParameter("phiBins")), - zBins_((int)iConfig.getParameter("zBins")), - d0CutNStubs4_((float)iConfig.getParameter("d0_cutNStubs4")), - d0CutNStubs5_((float)iConfig.getParameter("d0_cutNStubs5")), - lowpTJetMinTrackMultiplicity_((int)iConfig.getParameter("lowpTJetMinTrackMultiplicity")), - highpTJetMinTrackMultiplicity_((int)iConfig.getParameter("highpTJetMinTrackMultiplicity")), - displaced_(iConfig.getParameter("displaced")), - nStubs4DisplacedChi2_((float)iConfig.getParameter("nStubs4DisplacedChi2")), - nStubs5DisplacedChi2_((float)iConfig.getParameter("nStubs5DisplacedChi2")), - nStubs4DisplacedBend_((float)iConfig.getParameter("nStubs4Displacedbend")), - nStubs5DisplacedBend_((float)iConfig.getParameter("nStubs5Displacedbend")), - nDisplacedTracks_((int)iConfig.getParameter("nDisplacedTracks")), - dzPVTrk_((float)iConfig.getParameter("MaxDzTrackPV")), - trackToken_( - consumes > >(iConfig.getParameter("L1TrackInputTag"))), - PVtxToken_(consumes >(iConfig.getParameter("L1PVertexCollection"))), - tTopoToken_(esConsumes(edm::ESInputTag("", ""))) { - zStep_ = 2.0 * trkZMax_ / (zBins_ + 1); + : trackToken_( consumes > >(iConfig.getParameter("L1TrackInputTag")) ), + PVtxToken_( consumes< vector >(iConfig.getParameter("L1PVertexCollection")) ) + { + trkZMax_ = (float)iConfig.getParameter("trk_zMax"); + trkPtMax_ = (float)iConfig.getParameter("trk_ptMax"); + trkPtMin_ = (float)iConfig.getParameter("trk_ptMin"); + trkEtaMax_ = (float)iConfig.getParameter("trk_etaMax"); + nStubs4PromptChi2_ = (float)iConfig.getParameter("nStubs4PromptChi2"); + nStubs5PromptChi2_ = (float)iConfig.getParameter("nStubs5PromptChi2"); + nStubs4PromptBend_ = (float)iConfig.getParameter("nStubs4PromptBend"); + nStubs5PromptBend_ = (float)iConfig.getParameter("nStubs5PromptBend"); + trkNPSStubMin_ = (int)iConfig.getParameter("trk_nPSStubMin"); + minTrkJetpT_ = iConfig.getParameter("minTrkJetpT"); + etaBins_ = (int)iConfig.getParameter("etaBins"); + phiBins_ = (int)iConfig.getParameter("phiBins"); + zBins_ = (int)iConfig.getParameter("zBins"); + d0CutNStubs4_ = (float)iConfig.getParameter("d0_cutNStubs4"); + d0CutNStubs5_ = (float)iConfig.getParameter("d0_cutNStubs5"); + lowpTJetMinTrackMultiplicity_ = (int)iConfig.getParameter("lowpTJetMinTrackMultiplicity"); + highpTJetMinTrackMultiplicity_ = (int)iConfig.getParameter("highpTJetMinTrackMultiplicity"); + displaced_ = iConfig.getParameter("displaced"); + nStubs4DisplacedChi2_ = (float)iConfig.getParameter("nStubs4DisplacedChi2"); + nStubs5DisplacedChi2_ = (float)iConfig.getParameter("nStubs5DisplacedChi2"); + nStubs4DisplacedBend_ = (float)iConfig.getParameter("nStubs4DisplacedBend"); + nStubs5DisplacedBend_ = (float)iConfig.getParameter("nStubs5DisplacedBend"); + nDisplacedTracks_ = (int)iConfig.getParameter("nDisplacedTracks"); + dzPVTrk_ = (float)iConfig.getParameter("MaxDzTrackPV"); + + zStep_ = 2.0 * trkZMax_ / (zBins_+1); // added +1 in denom etaStep_ = 2.0 * trkEtaMax_ / etaBins_; //etaStep is the width of an etabin phiStep_ = 2 * M_PI / phiBins_; ////phiStep is the width of a phibin @@ -142,509 +143,236 @@ L1TrackJetProducer::~L1TrackJetProducer() {} void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { unique_ptr L1L1TrackJetProducer(new TkJetCollection); - // For TTStubs - const TrackerTopology &tTopo = iSetup.getData(tTopoToken_); + // Read inputs + ESHandle tTopoHandle; + iSetup.get().get(tTopoHandle); + ESHandle tGeomHandle; + iSetup.get().get(tGeomHandle); + const TrackerTopology *const tTopo = tTopoHandle.product(); edm::Handle > > TTTrackHandle; iEvent.getByToken(trackToken_, TTTrackHandle); - vector >::const_iterator iterL1Track; - edm::Handle > PVtx; + edm::Handle< std::vector > PVtx; iEvent.getByToken(PVtxToken_, PVtx); - float PVz = (PVtx->at(0)).z0(); + float PVz = (PVtx->at(0)).z0(); L1TrkPtrs_.clear(); - zBinCount_.clear(); tdtrk_.clear(); - unsigned int this_l1track = 0; - for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) { + + for (unsigned int this_l1track=0; this_l1tracksize(); this_l1track++){ edm::Ptr trkPtr(TTTrackHandle, this_l1track); - this_l1track++; float trk_pt = trkPtr->momentum().perp(); int trk_nstubs = (int)trkPtr->getStubRefs().size(); float trk_chi2dof = trkPtr->chi2Red(); float trk_d0 = trkPtr->d0(); float trk_bendchi2 = trkPtr->stubPtConsistency(); - float trk_z0 = trkPtr->z0(); - int trk_nPS = 0; + int trk_nPS = 0; for (int istub = 0; istub < trk_nstubs; istub++) { // loop over the stubs DetId detId(trkPtr->getStubRefs().at(istub)->getDetId()); if (detId.det() == DetId::Detector::Tracker) { - if ((detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3) || - (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9)) + if ((detId.subdetId() == StripSubdetector::TOB && tTopo->tobLayer(detId) <= 3) || + (detId.subdetId() == StripSubdetector::TID && tTopo->tidRing(detId) <= 9)) trk_nPS++; } } + // select tracks if (trk_nPS < trkNPSStubMin_) - continue; + continue; if (!trackQualityCuts(trk_nstubs, trk_chi2dof, trk_bendchi2)) continue; - if (std::abs(PVz - trk_z0) > dzPVTrk_) + if ( fabs(PVz-trkPtr->z0()) > dzPVTrk_ && dzPVTrk_>0 ) continue; - if (std::abs(trk_z0) > trkZMax_) + if (fabs(trkPtr->z0()) > trkZMax_) continue; - if (std::abs(trkPtr->momentum().eta()) > trkEtaMax_) + if ( fabs(trkPtr->momentum().eta())>trkEtaMax_) continue; if (trk_pt < trkPtMin_) continue; L1TrkPtrs_.push_back(trkPtr); - zBinCount_.push_back(0); - if ((std::abs(trk_d0) > d0CutNStubs5_ && trk_nstubs >= 5 && d0CutNStubs5_ >= 0) || - (trk_nstubs == 4 && std::abs(trk_d0) > d0CutNStubs4_ && d0CutNStubs4_ >= 0)) - tdtrk_.push_back(1); //displaced track + if ((fabs(trk_d0) > d0CutNStubs5_ && trk_nstubs >= 5 && d0CutNStubs5_>=0) || (trk_nstubs == 4 && fabs(trk_d0) > d0CutNStubs4_ && d0CutNStubs4_>=0)) + tdtrk_.push_back(1); //displaced track else tdtrk_.push_back(0); // not displaced track } if (!L1TrkPtrs_.empty()) { - MaxZBin mzb; - - L2_cluster(L1TrkPtrs_, tdtrk_, mzb); - vector > L1TrackAssocJet; - if (mzb.clusters != nullptr) { - for (int j = 0; j < mzb.nclust; ++j) { - //FILL Two Layer Jets for Jet Collection - if (mzb.clusters[j].pTtot <= trkPtMin_) - continue; //protects against reading bad memory - if (mzb.clusters[j].numtracks < 1) - continue; - if (mzb.clusters[j].numtracks > 5000) - continue; - float jetEta = mzb.clusters[j].eta; - float jetPhi = mzb.clusters[j].phi; - float jetPt = mzb.clusters[j].pTtot; - float jetPx = jetPt * cos(jetPhi); - float jetPy = jetPt * sin(jetPhi); - float jetPz = jetPt * sinh(jetEta); - float jetP = jetPt * cosh(jetEta); - int totalDisptrk = mzb.clusters[j].numtdtrks; - bool isDispJet = totalDisptrk >= nDisplacedTracks_; - math::XYZTLorentzVector jetP4(jetPx, jetPy, jetPz, jetP); - L1TrackAssocJet.clear(); - for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++) { - L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]); - } - - TkJet trkJet(jetP4, L1TrackAssocJet, mzb.zbincenter, mzb.clusters[j].numtracks, 0, totalDisptrk, 0, isDispJet); - if (!L1TrackAssocJet.empty()) { - L1L1TrackJetProducer->push_back(trkJet); - } - } - } - - if (displaced_) - iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended"); - else - iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets"); - delete[] mzb.clusters; + MaxZBin mzb; + mzb.znum=-1; + mzb.nclust=-1; + mzb.ht=-1; + EtaPhiBin epbins_default[phiBins_][etaBins_]; // create grid of phiBins + + float phi = -1.0 * M_PI; + float eta = -1.0 * trkEtaMax_; + for (int i = 0; i < phiBins_; ++i) { + eta = -1.0 * trkEtaMax_; + for (int j = 0; j < etaBins_; ++j) { + epbins_default[i][j].phi = (phi + (phi + phiStep_) ) / 2.0; // phimin=phi; phimax= phimin+step + epbins_default[i][j].eta = (eta + (eta + etaStep_) ) / 2.0; // phimin=phi; phimax phimin+step + eta += etaStep_; + } // for each etabin + phi += phiStep_; + } // for each phibin (finished creating epbins) + + std::vector zmins,zmaxs; + for (int zbin = 0; zbin < zBins_; zbin++) { + zmins.push_back(-1.0 * trkZMax_ + zStep_*zbin); + zmaxs.push_back(-1.0 * trkZMax_ + zStep_*zbin + zStep_*2); } -} -void L1TrackJetProducer::L2_cluster(vector > L1TrkPtrs_, vector tdtrk_, MaxZBin &mzb) { - const int nz = zBins_ + 1; - MaxZBin all_zBins[nz]; - MaxZBin mzbtemp = {}; - for (int z = 0; z < nz; ++z) - all_zBins[z] = mzbtemp; - - float zmin = -1.0 * trkZMax_; - float zmax = zmin + 2 * zStep_; - - EtaPhiBin epbins[phiBins_][etaBins_]; // create grid of phiBins - float phi = -1.0 * M_PI; - float eta; - float etamin, etamax, phimin, phimax; - - for (int i = 0; i < phiBins_; ++i) { - eta = -1.0 * trkEtaMax_; - for (int j = 0; j < etaBins_; ++j) { - phimin = phi; - phimax = phi + phiStep_; - etamin = eta; - eta = eta + etaStep_; - etamax = eta; - epbins[i][j].phi = (phimin + phimax) / 2.0; - epbins[i][j].eta = (etamin + etamax) / 2.0; - } // for each etabin - phi = phi + phiStep_; - } // for each phibin (finished creating epbins) - - mzb = all_zBins[0]; - int ntracks = L1TrkPtrs_.size(); - - // uninitalized arrays - EtaPhiBin *L1clusters[phiBins_]; - EtaPhiBin L2cluster[ntracks]; - - for (int zbin = 0; zbin < zBins_; ++zbin) { - for (int i = 0; i < phiBins_; ++i) { //First initialize pT, numtracks, used to 0 (or false) - for (int j = 0; j < etaBins_; ++j) { - epbins[i][j].pTtot = 0; - epbins[i][j].used = false; - epbins[i][j].numtracks = 0; - epbins[i][j].numttrks = 0; - epbins[i][j].numtdtrks = 0; - epbins[i][j].numttdtrks = 0; - epbins[i][j].trackidx.clear(); - } //for each etabin - L1clusters[i] = epbins[i]; - } //for each phibin - - for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { - float trkpt = L1TrkPtrs_[k]->momentum().perp(); - float trketa = L1TrkPtrs_[k]->momentum().eta(); - float trkphi = L1TrkPtrs_[k]->momentum().phi(); - float trkZ = L1TrkPtrs_[k]->z0(); - - for (int i = 0; i < phiBins_; ++i) { - for (int j = 0; j < etaBins_; ++j) { - L2cluster[k] = epbins[i][j]; - if ((zmin <= trkZ && zmax >= trkZ) && - ((epbins[i][j].eta - etaStep_ / 2.0 < trketa && epbins[i][j].eta + etaStep_ / 2.0 >= trketa) && - epbins[i][j].phi - phiStep_ / 2.0 < trkphi && epbins[i][j].phi + phiStep_ / 2.0 >= trkphi && - (zBinCount_[k] != 2))) { - zBinCount_.at(k) = zBinCount_.at(k) + 1; - if (trkpt < trkPtMax_) - epbins[i][j].pTtot += trkpt; - else - epbins[i][j].pTtot += trkPtMax_; - epbins[i][j].numtdtrks += tdtrk_[k]; - epbins[i][j].trackidx.push_back(k); - ++epbins[i][j].numtracks; - } // if right bin - } // for each phibin: j loop - } // for each phibin: i loop - } // end loop over tracks - - for (int phislice = 0; phislice < phiBins_; ++phislice) { - L1clusters[phislice] = L1_cluster(epbins[phislice]); - if (L1clusters[phislice] != nullptr) { - for (int ind = 0; L1clusters[phislice][ind].pTtot != 0; ++ind) { - L1clusters[phislice][ind].used = false; - } - } + // create vectors that hold data + std::vector> L1clusters; + L1clusters.reserve(phiBins_); + std::vector L2clusters; + + for (int i = 0; i < phiBins_; ++i) { + for (int j = 0; j < etaBins_; ++j) { + epbins_default[i][j].pTtot = 0; + epbins_default[i][j].used = false; + epbins_default[i][j].numtracks = 0; + epbins_default[i][j].numttrks = 0; + epbins_default[i][j].numtdtrks = 0; + epbins_default[i][j].numttdtrks = 0; + epbins_default[i][j].trackidx.clear(); } - - //Create clusters array to hold output cluster data for Layer2; can't have more clusters than tracks. - //Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used. - float hipT = 0; - int nclust = 0; - int phibin = 0; - int imax = -1; - int index1; //index of clusters array for each phislice - float E1 = 0; - float E0 = 0; - float E2 = 0; - int trx1, trx2; - int tdtrk1, tdtrk2; - int used1, used2, used3, used4; - - for (phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT - while (true) { - hipT = 0; - for (index1 = 0; L1clusters[phibin][index1].pTtot > 0; ++index1) { - if (!L1clusters[phibin][index1].used && L1clusters[phibin][index1].pTtot >= hipT) { - hipT = L1clusters[phibin][index1].pTtot; - imax = index1; - } - } // for each index within the phibin - - if (hipT == 0) - break; // If highest pT is 0, all bins are used - E0 = hipT; // E0 is pT of first phibin of the cluster - E1 = 0; - E2 = 0; - trx1 = 0; - trx2 = 0; - tdtrk1 = 0; - tdtrk2 = 0; - std::vector trkidx1; - std::vector trkidx2; - L2cluster[nclust] = L1clusters[phibin][imax]; - L1clusters[phibin][imax].used = true; - // Add pT of upper neighbor - // E1 is pT of the middle phibin (should be highest pT) - if (phibin != phiBins_ - 1) { - used1 = -1; - used2 = -1; - for (index1 = 0; L1clusters[phibin + 1][index1].pTtot != 0; ++index1) { - if (L1clusters[phibin + 1][index1].used) - continue; - - if (std::abs(L1clusters[phibin + 1][index1].eta - L1clusters[phibin][imax].eta) <= 1.5 * etaStep_) { - E1 += L1clusters[phibin + 1][index1].pTtot; - trx1 += L1clusters[phibin + 1][index1].numtracks; - tdtrk1 += L1clusters[phibin + 1][index1].numtdtrks; - for (unsigned int itrk = 0; itrk < L1clusters[phibin + 1][index1].trackidx.size(); itrk++) { - trkidx1.push_back(L1clusters[phibin + 1][index1].trackidx[itrk]); - } - if (used1 < 0) - used1 = index1; - else - used2 = index1; - } // if cluster is within one phibin - } // for each cluster in above phibin - - if (E1 < E0) { // if E1 isn't higher, E0 and E1 are their own cluster - L2cluster[nclust].pTtot += E1; - L2cluster[nclust].numtracks += trx1; - L2cluster[nclust].numtdtrks += tdtrk1; - for (unsigned int itrk = 0; itrk < trkidx1.size(); itrk++) { - L2cluster[nclust].trackidx.push_back(trkidx1[itrk]); - } - - if (used1 >= 0) - L1clusters[phibin + 1][used1].used = true; - if (used2 >= 0) - L1clusters[phibin + 1][used2].used = true; - nclust++; - continue; - } - - if (phibin != phiBins_ - 2) { // E2 will be the pT of the third phibin (should be lower than E1) - used3 = -1; - used4 = -1; - for (index1 = 0; L1clusters[phibin + 2][index1].pTtot != 0; ++index1) { - if (L1clusters[phibin + 2][index1].used) - continue; - if (std::abs(L1clusters[phibin + 2][index1].eta - L1clusters[phibin][imax].eta) <= 1.5 * etaStep_) { - E2 += L1clusters[phibin + 2][index1].pTtot; - trx2 += L1clusters[phibin + 2][index1].numtracks; - tdtrk2 += L1clusters[phibin + 2][index1].numtdtrks; - for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][index1].trackidx.size(); itrk++) - trkidx2.push_back(L1clusters[phibin + 2][index1].trackidx[itrk]); - - if (used3 < 0) - used3 = index1; - else - used4 = index1; - } - } - - // if indeed E2 < E1, add E1 and E2 to E0, they're all a cluster together - // otherwise, E0 is its own cluster - if (E2 < E1) { - L2cluster[nclust].pTtot += E1 + E2; - L2cluster[nclust].numtracks += trx1 + trx2; - L2cluster[nclust].numtdtrks += tdtrk1 + tdtrk2; - L2cluster[nclust].phi = L1clusters[phibin + 1][used1].phi; - for (unsigned int itrk = 0; itrk < trkidx1.size(); itrk++) { - L2cluster[nclust].trackidx.push_back(trkidx1[itrk]); - } - for (unsigned int itrk = 0; itrk < trkidx2.size(); itrk++) { - L2cluster[nclust].trackidx.push_back(trkidx2[itrk]); - } - - if (used1 >= 0) - L1clusters[phibin + 1][used1].used = true; - if (used2 >= 0) - L1clusters[phibin + 1][used2].used = true; - if (used3 >= 0) - L1clusters[phibin + 2][used3].used = true; - if (used4 >= 0) - L1clusters[phibin + 2][used4].used = true; - } - nclust++; + } + + for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) { + // initialize matrices + float zmin = zmins[zbin]; + float zmax = zmaxs[zbin]; + EtaPhiBin epbins[phiBins_][etaBins_]; + std::copy(&epbins_default[0][0], &epbins_default[0][0]+phiBins_*etaBins_,&epbins[0][0]); + + //clear containers + L1clusters.clear(); + L2clusters.clear(); + + // fill grid + for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { + float trkpt = L1TrkPtrs_[k]->momentum().perp(); + float trketa = L1TrkPtrs_[k]->momentum().eta(); + float trkphi = L1TrkPtrs_[k]->momentum().phi(); + float trkZ = L1TrkPtrs_[k]->z0(); + if ( zmax < trkZ) + continue; + if (zbin==0){ + if ( zmin > trkZ ) continue; - } // end Not phiBins-2 - else { - L2cluster[nclust].pTtot += E1; - L2cluster[nclust].numtracks += trx1; - L2cluster[nclust].numtdtrks += tdtrk1; - L2cluster[nclust].phi = L1clusters[phibin + 1][used1].phi; - for (unsigned int itrk = 0; itrk < trkidx1.size(); itrk++) { - L2cluster[nclust].trackidx.push_back(trkidx1[itrk]); - } - - if (used1 >= 0) - L1clusters[phibin + 1][used1].used = true; - if (used2 >= 0) - L1clusters[phibin + 1][used2].used = true; - nclust++; + } else{ + if ( zmin > trkZ || zmin ==trkZ ) continue; - } - } //End not last phibin(23) - else { //if it is phibin 23 - L1clusters[phibin][imax].used = true; - nclust++; - } - } // while hipT not 0 - } // for each phibin - - for (phibin = 0; phibin < phiBins_; ++phibin) - delete[] L1clusters[phibin]; - - // Now merge clusters, if necessary - for (int m = 0; m < nclust - 1; ++m) { - for (int n = m + 1; n < nclust; ++n) { - if (L2cluster[n].eta == L2cluster[m].eta && (std::abs(L2cluster[n].phi - L2cluster[m].phi) < 1.5 * phiStep_ || - std::abs(L2cluster[n].phi - L2cluster[m].phi) > 6.0)) { - if (L2cluster[n].pTtot > L2cluster[m].pTtot) { - L2cluster[m].phi = L2cluster[n].phi; - } - L2cluster[m].pTtot += L2cluster[n].pTtot; - L2cluster[m].numtracks += L2cluster[n].numtracks; - L2cluster[m].numtdtrks += L2cluster[n].numtdtrks; - for (unsigned int itrk = 0; itrk < L2cluster[n].trackidx.size(); itrk++) - L2cluster[m].trackidx.push_back(L2cluster[n].trackidx[itrk]); - - for (int m1 = n; m1 < nclust - 1; ++m1) { - L2cluster[m1] = L2cluster[m1 + 1]; - } - nclust--; - m = -1; - break; //????? - } // end if clusters neighbor in eta - } - } // end for (m) loop - - // sum up all pTs in this zbin to find ht - // Require jet to have at least 2 tracks with a summed eT>50 GeV, or 3 tracks with summed eT>100 GeV - // in order to count toward HT - float ht = 0; - for (int k = 0; k < nclust; ++k) { - if (L2cluster[k].pTtot > minJetEtLowPt_ && L2cluster[k].numtracks < lowpTJetMinTrackMultiplicity_) + } + for (int i = 0; i < phiBins_; ++i) { + for (int j = 0; j < etaBins_; ++j) { + float eta_min=epbins[i][j].eta - etaStep_ / 2.0; //eta min + float eta_max=epbins[i][j].eta + etaStep_ / 2.0; //eta max + float phi_min=epbins[i][j].phi - phiStep_ / 2.0; //phi min + float phi_max=epbins[i][j].phi + phiStep_ / 2.0; //phi max + + if ( ( trketa < eta_min ) || ( trketa > eta_max ) || + ( trkphi < phi_min ) || ( trkphi > phi_max ) ) + continue; + + if (trkpt < trkPtMax_) + epbins[i][j].pTtot += trkpt; + else + epbins[i][j].pTtot += trkPtMax_; + epbins[i][j].numtdtrks += tdtrk_[k]; + epbins[i][j].trackidx.push_back(k); + ++epbins[i][j].numtracks; + } // for each phibin + } // for each phibin + } //end loop over tracks + + // first layer clustering - in eta for all phi bins + for (int phibin = 0; phibin < phiBins_; ++phibin) { + L1clusters.push_back( L1_clustering(epbins[phibin],etaBins_,etaStep_) ); + } + + //second layer clustering in phi for all eta clusters + L2clusters= L2_clustering(L1clusters, phiBins_, phiStep_, etaStep_); + + // sum all cluster pTs in this zbin to find max + float sum_pt = 0; + for (unsigned int k = 0; k < L2clusters.size(); ++k) { + if (L2clusters[k].pTtot > 50 && L2clusters[k].numtracks < lowpTJetMinTrackMultiplicity_) continue; - if (L2cluster[k].pTtot > minJetEtHighPt_ && L2cluster[k].numtracks < highpTJetMinTrackMultiplicity_) + if (L2clusters[k].pTtot > 100 && L2clusters[k].numtracks < highpTJetMinTrackMultiplicity_) continue; - if (L2cluster[k].pTtot > minTrkJetpT_) - ht += L2cluster[k].pTtot; - } - - // if ht is larger than previous max, this is the new vertex zbin - all_zBins[zbin].znum = zbin; - all_zBins[zbin].clusters = new EtaPhiBin[nclust]; - all_zBins[zbin].nclust = nclust; - all_zBins[zbin].zbincenter = (zmin + zmax) / 2.0; - for (int k = 0; k < nclust; ++k) { - all_zBins[zbin].clusters[k].phi = L2cluster[k].phi; - all_zBins[zbin].clusters[k].eta = L2cluster[k].eta; - all_zBins[zbin].clusters[k].pTtot = L2cluster[k].pTtot; - all_zBins[zbin].clusters[k].numtracks = L2cluster[k].numtracks; - all_zBins[zbin].clusters[k].numtdtrks = L2cluster[k].numtdtrks; - for (unsigned int itrk = 0; itrk < L2cluster[k].trackidx.size(); itrk++) - all_zBins[zbin].clusters[k].trackidx.push_back(L2cluster[k].trackidx[itrk]); + if (L2clusters[k].pTtot > minTrkJetpT_) + sum_pt += L2clusters[k].pTtot; + } + + if (sum_pt < mzb.ht) continue; + // if ht is larger than previous max, this is the new vertex zbin + mzb.ht=sum_pt; + mzb.znum = zbin; + mzb.clusters=L2clusters; + mzb.nclust = L2clusters.size(); + mzb.zbincenter = (zmin + zmax) / 2.0; + }//zbin loop + + + vector > L1TrackAssocJet; + for (unsigned int j = 0; j < mzb.clusters.size(); ++j) { + float jetEta = mzb.clusters[j].eta; + float jetPhi = mzb.clusters[j].phi; + float jetPt = mzb.clusters[j].pTtot; + float jetPx = jetPt * cos(jetPhi); + float jetPy = jetPt * sin(jetPhi); + float jetPz = jetPt * sinh(jetEta); + float jetP = jetPt * cosh(jetEta); + int totalDisptrk = mzb.clusters[j].numtdtrks; + bool isDispJet=false; + if (totalDisptrk>nDisplacedTracks_ || totalDisptrk==nDisplacedTracks_) + isDispJet=true; + + math::XYZTLorentzVector jetP4(jetPx, jetPy, jetPz, jetP); + L1TrackAssocJet.clear(); + for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++) + L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]); + + TkJet trkJet(jetP4, L1TrackAssocJet, mzb.zbincenter, + mzb.clusters[j].numtracks, 0, totalDisptrk, 0, isDispJet); + + if (!L1TrackAssocJet.empty()) + L1L1TrackJetProducer->push_back(trkJet); } + - all_zBins[zbin].ht = ht; - if (ht >= mzb.ht) { - mzb = all_zBins[zbin]; - mzb.zbincenter = (zmin + zmax) / 2.0; - } - // Prepare for next zbin! - zmin = zmin + zStep_; - zmax = zmax + zStep_; - } // for each zbin - for (int zbin = 0; zbin < zBins_; ++zbin) { - if (zbin == mzb.znum) - continue; - delete[] all_zBins[zbin].clusters; + if (displaced_) + iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended"); + else + iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets"); } } -EtaPhiBin *L1TrackJetProducer::L1_cluster(EtaPhiBin *phislice) { - EtaPhiBin *clusters = new EtaPhiBin[etaBins_ / 2]; - if (clusters == nullptr) - edm::LogWarning("L1TrackJetProducer") << "Clusters memory not assigned!\n"; - - // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used - float my_pt, left_pt, right_pt, right2pt; - int nclust = 0; - right2pt = 0; - for (int etabin = 0; etabin < etaBins_; ++etabin) { - // assign values for my pT and neighbors' pT - if (phislice[etabin].used) - continue; - my_pt = phislice[etabin].pTtot; - if (etabin > 0 && !phislice[etabin - 1].used) { - left_pt = phislice[etabin - 1].pTtot; - } else - left_pt = 0; - if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) { - right_pt = phislice[etabin + 1].pTtot; - if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) { - right2pt = phislice[etabin + 2].pTtot; - } else - right2pt = 0; - } else - right_pt = 0; - - // if I'm not a cluster, move on - if (my_pt < left_pt || my_pt <= right_pt) { - // if unused pT in the left neighbor, spit it out as a cluster - if (left_pt > 0) { - clusters[nclust] = phislice[etabin - 1]; - phislice[etabin - 1].used = true; - nclust++; - } - continue; - } - - // I guess I'm a cluster-- should I use my right neighbor? - // Note: left neighbor will definitely be used because if it - // didn't belong to me it would have been used already - clusters[nclust] = phislice[etabin]; - phislice[etabin].used = true; - if (left_pt > 0) { - if (clusters != nullptr) { - clusters[nclust].pTtot += left_pt; - clusters[nclust].numtracks += phislice[etabin - 1].numtracks; - clusters[nclust].numtdtrks += phislice[etabin - 1].numtdtrks; - for (unsigned int itrk = 0; itrk < phislice[etabin - 1].trackidx.size(); itrk++) - clusters[nclust].trackidx.push_back(phislice[etabin - 1].trackidx[itrk]); - } - } - if (my_pt >= right2pt && right_pt > 0) { - if (clusters != nullptr) { - clusters[nclust].pTtot += right_pt; - clusters[nclust].numtracks += phislice[etabin + 1].numtracks; - clusters[nclust].numtdtrks += phislice[etabin + 1].numtdtrks; - for (unsigned int itrk = 0; itrk < phislice[etabin + 1].trackidx.size(); itrk++) - clusters[nclust].trackidx.push_back(phislice[etabin + 1].trackidx[itrk]); - phislice[etabin + 1].used = true; - } - } - nclust++; - } // for each etabin - - // Now merge clusters, if necessary - for (int m = 0; m < nclust - 1; ++m) { - if (std::abs(clusters[m + 1].eta - clusters[m].eta) < 1.5 * etaStep_) { - if (clusters[m + 1].pTtot > clusters[m].pTtot) { - clusters[m].eta = clusters[m + 1].eta; - } - clusters[m].pTtot += clusters[m + 1].pTtot; - clusters[m].numtracks += clusters[m + 1].numtracks; // Previous version didn't add tracks when merging - clusters[m].numtdtrks += clusters[m + 1].numtdtrks; - for (unsigned int itrk = 0; itrk < clusters[m + 1].trackidx.size(); itrk++) - clusters[m].trackidx.push_back(clusters[m + 1].trackidx[itrk]); - - for (int m1 = m + 1; m1 < nclust - 1; ++m1) - clusters[m1] = clusters[m1 + 1]; - nclust--; - m = -1; - } // end if clusters neighbor in eta - } // end for (m) loop - - for (int i = nclust; i < etaBins_ / 2; ++i) // zero out remaining unused clusters - clusters[i].pTtot = 0; - return clusters; -} - void L1TrackJetProducer::beginStream(StreamID) {} void L1TrackJetProducer::endStream() {} -bool L1TrackJetProducer::trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2) { +bool L1TrackJetProducer::trackQualityCuts( + int trk_nstub, float trk_chi2, float trk_bendchi2) { bool PassQuality = false; - if (trk_bendchi2 < trkBendChi2Max_ && trk_chi2 < trkChi2dofMax_ && trk_nstub >= 4 && !displaced_) - PassQuality = true; - if (displaced_ && trk_bendchi2 < nStubs4DisplacedBend_ && trk_chi2 < nStubs4DisplacedChi2_ && trk_nstub == 4) - PassQuality = true; - if (displaced_ && trk_bendchi2 < nStubs5DisplacedBend_ && trk_chi2 < nStubs5DisplacedChi2_ && trk_nstub > 4) - PassQuality = true; + if (!displaced_){ + if ( trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ && + trk_chi2 < nStubs4PromptChi2_) + PassQuality = true; + if ( trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ && + trk_chi2 < nStubs5PromptChi2_) + PassQuality = true; + } else{ + if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ && + trk_chi2 < nStubs4DisplacedChi2_ ) + PassQuality = true; + if (trk_nstub >4 && trk_bendchi2 < nStubs5DisplacedBend_ && + trk_chi2 < nStubs5DisplacedChi2_ ) + PassQuality = true; + } return PassQuality; } diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h new file mode 100644 index 0000000000000..841a64c995e3a --- /dev/null +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h @@ -0,0 +1,250 @@ +#pragma once +#include +#include +#include +#include +#include +using namespace std; + +//Each individual box in the eta and phi dimension. +// Also used to store final cluster data for each zbin. +struct EtaPhiBin { + float pTtot; + int numtracks; + int numttrks; + int numtdtrks; + int numttdtrks; + bool used; + float phi; //average phi value (halfway b/t min and max) + float eta; //average eta value + std::vector trackidx; +}; + +//store important information for plots +struct MaxZBin { + int znum; //Numbered from 0 to nzbins (16, 32, or 64) in order + int nclust; //number of clusters in this bin + float zbincenter; + std::vector clusters; //list of all the clusters in this bin + float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored +}; + + +std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_ ) { + + std::vector clusters; + // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used + int nclust = 0; + + // get tracks in eta bins in increasing eta order + for (int etabin = 0; etabin < etaBins_; ++etabin) { + float my_pt=0, previousbin_pt=0;//, nextbin_pt=0, next2bin_pt=0; + float nextbin_pt=0,nextbin2_pt=0; + + // skip (already) used tracks + if (phislice[etabin].used) + continue; + my_pt = phislice[etabin].pTtot; + if (my_pt==0) continue; +//get previous bin pT + if (etabin > 0 && !phislice[etabin - 1].used) + previousbin_pt = phislice[etabin - 1].pTtot; + + // get next bins pt + if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) { + nextbin_pt = phislice[etabin + 1].pTtot; + if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) { + nextbin2_pt = phislice[etabin + 2].pTtot; + } + } + // check if pT of current cluster is higher than neighbors + if (my_pt < previousbin_pt || my_pt <= nextbin_pt) { + // if unused pT in the left neighbor, spit it out as a cluster + if ( previousbin_pt > 0) { + clusters.push_back( phislice[etabin - 1] ); + phislice[etabin - 1].used = true; + nclust++; + } + continue; //if it is not the local max pT skip + } + // here reach only unused local max clusters + clusters.push_back( phislice[etabin]); + phislice[etabin].used = true; //if current bin a cluster + if (previousbin_pt > 0) { + clusters[nclust].pTtot += previousbin_pt; + clusters[nclust].numtracks += phislice[etabin - 1].numtracks; + clusters[nclust].numtdtrks += phislice[etabin - 1].numtdtrks; + for (unsigned int itrk=0; itrk= nextbin2_pt && nextbin_pt > 0) { + clusters[nclust].pTtot += nextbin_pt; + clusters[nclust].numtracks += phislice[etabin + 1].numtracks; + clusters[nclust].numtdtrks += phislice[etabin + 1].numtdtrks; + for (unsigned int itrk=0; itrk clusters[m].pTtot) { + clusters[m].eta = clusters[m + 1].eta; + } + clusters[m].pTtot += clusters[m + 1].pTtot; + clusters[m].numtracks += clusters[m + 1].numtracks; // total ntrk + clusters[m].numtdtrks += clusters[m + 1].numtdtrks; // total ndisp + for (unsigned int itrk=0; itrk trkidx){ + bin.pTtot+=pt; + bin.numtracks+=ntrk; + bin.numtdtrks+=ndtrk; + for (unsigned int itrk=0; itrk= M_PI) + x -= 2*M_PI; + if (x < -1*M_PI) + x += 2*M_PI; + return x; +} + +std::vector L2_clustering(std::vector> &L1clusters, int phiBins_, float phiStep_, float etaStep_ ) { + + std::vector clusters; + for (int phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT + if (L1clusters[phibin].size()==0) continue; + + // sort L1 clusters max -> min + sort(L1clusters[phibin].begin(),L1clusters[phibin].end(),[] ( struct EtaPhiBin &a, struct EtaPhiBin &b){ return a.pTtot>b.pTtot; } ); + for (unsigned int imax = 0; imax trkidx1; + std::vector trkidx2; + clusters.push_back( L1clusters[phibin][imax]); + + L1clusters[phibin][imax].used = true; + + // if we are in the last phi bin, dont check phi+1 phi+2 + if (phibin == phiBins_ -1) + continue; + std::vector used_already; //keep phi+1 clusters that have been used + for (unsigned int icluster = 0; icluster 1.5 * etaStep_) + continue; + pt_next += L1clusters[phibin + 1][icluster].pTtot; + trk1 += L1clusters[phibin + 1][icluster].numtracks; + tdtrk1 += L1clusters[phibin + 1][icluster].numtdtrks; + for (unsigned int itrk=0; itrk used_already2; //keep used clusters in phi+2 + for (unsigned int icluster = 0; icluster 1.5 * etaStep_) + continue; + pt_next2 += L1clusters[phibin + 2][icluster].pTtot; + trk2 += L1clusters[phibin + 2][icluster].numtracks; + tdtrk2 += L1clusters[phibin + 2][icluster].numtdtrks; + for (unsigned int itrk=0; itrk trkidx_both; + trkidx_both.reserve(trkidx1.size()+trkidx2.size()); + trkidx_both.insert(trkidx_both.end(),trkidx1.begin(),trkidx1.end()); + trkidx_both.insert(trkidx_both.end(),trkidx2.begin(),trkidx2.end()); + Fill_L2Cluster( clusters[clusters.size()-1], pt_next+pt_next2, trk1+trk2 , tdtrk1+tdtrk2, trkidx_both); + clusters[clusters.size()-1].phi = L1clusters[phibin + 1][used_already[0]].phi; + for(unsigned int iused: used_already) + L1clusters[phibin + 1][iused].used =true; + for(unsigned int iused: used_already2) + L1clusters[phibin + 2][iused].used =true; + } + } // for each L1 cluster + } // for each phibin + + int nclust=clusters.size(); + + // merge close-by clusters + for (int m = 0; m < nclust - 1; ++m) { + for (int n = m + 1; n < nclust; ++n){ + if (clusters[n].eta != clusters[m].eta) + continue; + if ( fabs(DPhi(clusters[n].phi, clusters[m].phi))>1.5 * phiStep_) + continue; + + if (clusters[n].pTtot > clusters[m].pTtot) + clusters[m].phi = clusters[n].phi; + + clusters[m].pTtot += clusters[n].pTtot; + clusters[m].numtracks += clusters[n].numtracks; + clusters[m].numtdtrks += clusters[n].numtdtrks; + for(unsigned int itrk=0; itrkcut excluded + MaxDzTrackPV = cms.double( 5.0 ), # tracks with dz(trk,PV)>cut excluded trk_zMax = cms.double (15.) , # max track z trk_ptMax = cms.double(200.), # maxi track pT before saturation - trk_ptMin = cms.double(3.0), # min track pt + trk_ptMin = cms.double(3.0), # min track pt trk_etaMax = cms.double(2.4), # max track eta - trk_chi2dofMax=cms.double(40.), # max track chi2/dof for prompt - trk_bendChi2Max=cms.double(40.), # max track bendchi2 for prompt + nStubs4PromptChi2=cms.double(5.0), #Prompt track quality flags for loose/tight + nStubs4PromptBend=cms.double(1.7), + nStubs5PromptChi2=cms.double(2.75), + nStubs5PromptBend=cms.double(3.5), trk_nPSStubMin=cms.int32(-1), # min # PS stubs, -1 means no cut - minTrkJetpT=cms.double(5.), # min track pt to be considered for track jet - etaBins=cms.int32(24), #number of eta bins - phiBins=cms.int32(27), #number of phi bins - zBins=cms.int32(10), #number of z bins - d0_cutNStubs4=cms.double(-1), # -1 excludes nstub=4 from disp tag + minTrkJetpT=cms.double(5.), # min track jet pt to be considered for most energetic zbin finding + etaBins=cms.int32(24), #number of eta bins + phiBins=cms.int32(27), #number of phi bins + zBins=cms.int32(1), #number of z bins + d0_cutNStubs4=cms.double(-1), # -1 excludes nstub=4 from disp tag d0_cutNStubs5=cms.double(0.22), # -1 excludes nstub>4 from disp tag - lowpTJetMinTrackMultiplicity=cms.int32(2), - highpTJetMinTrackMultiplicity=cms.int32(3), - minJetEtLowPt=cms.double(50.0), - minJetEtHighPt=cms.double(100.0), + lowpTJetMinTrackMultiplicity=cms.int32(2), #used only on zbin finding + highpTJetMinTrackMultiplicity=cms.int32(3), #used only on zbin finding displaced=cms.bool(True), #Flag for displaced tracks nStubs4DisplacedChi2=cms.double(3.3), #Disp tracks selection [trk4 from disp tag +# lowpTJetMinTrackMultiplicity=cms.int32(2), #used only on zbin finding +# highpTJetMinTrackMultiplicity=cms.int32(3), #used only on zbin finding +# displaced=cms.bool(True), #Flag for displaced tracks +# nStubs4DisplacedChi2=cms.double(3.3), #Disp tracks selection [trk Date: Tue, 3 May 2022 15:01:58 +0200 Subject: [PATCH 2/8] rebasing update track jets (cherry picked from commit 66b19f9179d6bba6b3590957b5a32a67689c65f5) (cherry picked from commit 4e8b4b038c29c5de21f9b465090e9766d03504a8) --- .../plugins/L1TrackJetProducer.cc | 303 +++++++++--------- .../plugins/L1TrackJetProducer.h | 266 +++++++-------- .../L1TTrackMatch/python/l1tTrackJets_cfi.py | 16 +- 3 files changed, 292 insertions(+), 293 deletions(-) diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc index 4d316f0af8196..b0c2bba6fdeff 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc @@ -29,7 +29,6 @@ #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" -#include "DataFormats/L1TCorrelator/interface/TkPrimaryVertex.h" #include "DataFormats/L1Trigger/interface/Vertex.h" #include "L1TrackJetProducer.h" @@ -47,7 +46,6 @@ using namespace std; using namespace edm; using namespace l1t; - class L1TrackJetProducer : public stream::EDProducer<> { public: explicit L1TrackJetProducer(const ParameterSet &); @@ -58,7 +56,6 @@ class L1TrackJetProducer : public stream::EDProducer<> { static void fillDescriptions(ConfigurationDescriptions &descriptions); bool trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2); - private: void beginStream(StreamID) override; void produce(Event &, const EventSetup &) override; @@ -66,9 +63,9 @@ class L1TrackJetProducer : public stream::EDProducer<> { // ----------member data --------------------------- - const EDGetTokenT > > trackToken_; - const edm::EDGetTokenT > PVtxToken_; - vector > L1TrkPtrs_; + const EDGetTokenT>> trackToken_; + const edm::EDGetTokenT> PVtxToken_; + vector> L1TrkPtrs_; vector tdtrk_; float trkZMax_; float trkPtMax_; @@ -80,7 +77,9 @@ class L1TrackJetProducer : public stream::EDProducer<> { float nStubs5PromptBend_; int trkNPSStubMin_; int lowpTJetMinTrackMultiplicity_; + float lowpTJetThreshold_; int highpTJetMinTrackMultiplicity_; + float highpTJetThreshold_; int zBins_; int etaBins_; int phiBins_; @@ -100,9 +99,8 @@ class L1TrackJetProducer : public stream::EDProducer<> { }; L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) - : trackToken_( consumes > >(iConfig.getParameter("L1TrackInputTag")) ), - PVtxToken_( consumes< vector >(iConfig.getParameter("L1PVertexCollection")) ) - { + : trackToken_(consumes>>(iConfig.getParameter("L1TrackInputTag"))), + PVtxToken_(consumes>(iConfig.getParameter("L1PVertexCollection"))) { trkZMax_ = (float)iConfig.getParameter("trk_zMax"); trkPtMax_ = (float)iConfig.getParameter("trk_ptMax"); trkPtMin_ = (float)iConfig.getParameter("trk_ptMin"); @@ -119,7 +117,9 @@ L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) d0CutNStubs4_ = (float)iConfig.getParameter("d0_cutNStubs4"); d0CutNStubs5_ = (float)iConfig.getParameter("d0_cutNStubs5"); lowpTJetMinTrackMultiplicity_ = (int)iConfig.getParameter("lowpTJetMinTrackMultiplicity"); + lowpTJetThreshold_ = (float)iConfig.getParameter("lowpTJetThreshold"); highpTJetMinTrackMultiplicity_ = (int)iConfig.getParameter("highpTJetMinTrackMultiplicity"); + highpTJetThreshold_ = (float)iConfig.getParameter("highpTJetThreshold"); displaced_ = iConfig.getParameter("displaced"); nStubs4DisplacedChi2_ = (float)iConfig.getParameter("nStubs4DisplacedChi2"); nStubs5DisplacedChi2_ = (float)iConfig.getParameter("nStubs5DisplacedChi2"); @@ -128,7 +128,7 @@ L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) nDisplacedTracks_ = (int)iConfig.getParameter("nDisplacedTracks"); dzPVTrk_ = (float)iConfig.getParameter("MaxDzTrackPV"); - zStep_ = 2.0 * trkZMax_ / (zBins_+1); // added +1 in denom + zStep_ = 2.0 * trkZMax_ / (zBins_ + 1); // added +1 in denom etaStep_ = 2.0 * trkEtaMax_ / etaBins_; //etaStep is the width of an etabin phiStep_ = 2 * M_PI / phiBins_; ////phiStep is the width of a phibin @@ -150,17 +150,17 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { iSetup.get().get(tGeomHandle); const TrackerTopology *const tTopo = tTopoHandle.product(); - edm::Handle > > TTTrackHandle; + edm::Handle>> TTTrackHandle; iEvent.getByToken(trackToken_, TTTrackHandle); - edm::Handle< std::vector > PVtx; + edm::Handle> PVtx; iEvent.getByToken(PVtxToken_, PVtx); - float PVz = (PVtx->at(0)).z0(); + float PVz = (PVtx->at(0)).z0(); L1TrkPtrs_.clear(); tdtrk_.clear(); - for (unsigned int this_l1track=0; this_l1tracksize(); this_l1track++){ + for (unsigned int this_l1track = 0; this_l1track < TTTrackHandle->size(); this_l1track++) { edm::Ptr trkPtr(TTTrackHandle, this_l1track); float trk_pt = trkPtr->momentum().perp(); int trk_nstubs = (int)trkPtr->getStubRefs().size(); @@ -168,7 +168,7 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { float trk_d0 = trkPtr->d0(); float trk_bendchi2 = trkPtr->stubPtConsistency(); - int trk_nPS = 0; + int trk_nPS = 0; for (int istub = 0; istub < trk_nstubs; istub++) { // loop over the stubs DetId detId(trkPtr->getStubRefs().at(istub)->getDetId()); if (detId.det() == DetId::Detector::Tracker) { @@ -179,145 +179,145 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { } // select tracks if (trk_nPS < trkNPSStubMin_) - continue; + continue; if (!trackQualityCuts(trk_nstubs, trk_chi2dof, trk_bendchi2)) continue; - if ( fabs(PVz-trkPtr->z0()) > dzPVTrk_ && dzPVTrk_>0 ) + if (std::abs(PVz - trkPtr->z0()) > dzPVTrk_ && dzPVTrk_ > 0) continue; - if (fabs(trkPtr->z0()) > trkZMax_) + if (std::abs(trkPtr->z0()) > trkZMax_) continue; - if ( fabs(trkPtr->momentum().eta())>trkEtaMax_) + if (std::abs(trkPtr->momentum().eta()) > trkEtaMax_) continue; if (trk_pt < trkPtMin_) continue; L1TrkPtrs_.push_back(trkPtr); - if ((fabs(trk_d0) > d0CutNStubs5_ && trk_nstubs >= 5 && d0CutNStubs5_>=0) || (trk_nstubs == 4 && fabs(trk_d0) > d0CutNStubs4_ && d0CutNStubs4_>=0)) - tdtrk_.push_back(1); //displaced track + if ((std::abs(trk_d0) > d0CutNStubs5_ && trk_nstubs >= 5 && d0CutNStubs5_ >= 0) || + (trk_nstubs == 4 && std::abs(trk_d0) > d0CutNStubs4_ && d0CutNStubs4_ >= 0)) + tdtrk_.push_back(1); //displaced track else tdtrk_.push_back(0); // not displaced track } if (!L1TrkPtrs_.empty()) { - MaxZBin mzb; - mzb.znum=-1; - mzb.nclust=-1; - mzb.ht=-1; - EtaPhiBin epbins_default[phiBins_][etaBins_]; // create grid of phiBins - - float phi = -1.0 * M_PI; - float eta = -1.0 * trkEtaMax_; - for (int i = 0; i < phiBins_; ++i) { - eta = -1.0 * trkEtaMax_; - for (int j = 0; j < etaBins_; ++j) { - epbins_default[i][j].phi = (phi + (phi + phiStep_) ) / 2.0; // phimin=phi; phimax= phimin+step - epbins_default[i][j].eta = (eta + (eta + etaStep_) ) / 2.0; // phimin=phi; phimax phimin+step - eta += etaStep_; - } // for each etabin - phi += phiStep_; - } // for each phibin (finished creating epbins) - - std::vector zmins,zmaxs; - for (int zbin = 0; zbin < zBins_; zbin++) { - zmins.push_back(-1.0 * trkZMax_ + zStep_*zbin); - zmaxs.push_back(-1.0 * trkZMax_ + zStep_*zbin + zStep_*2); - } + MaxZBin mzb; + mzb.znum = -1; + mzb.nclust = -1; + mzb.ht = -1; + EtaPhiBin epbins_default[phiBins_][etaBins_]; // create grid of phiBins + + float phi = -1.0 * M_PI; + float eta = -1.0 * trkEtaMax_; + for (int i = 0; i < phiBins_; ++i) { + eta = -1.0 * trkEtaMax_; + for (int j = 0; j < etaBins_; ++j) { + epbins_default[i][j].phi = (phi + (phi + phiStep_)) / 2.0; // phimin=phi; phimax= phimin+step + epbins_default[i][j].eta = (eta + (eta + etaStep_)) / 2.0; // phimin=phi; phimax phimin+step + eta += etaStep_; + } // for each etabin + phi += phiStep_; + } // for each phibin (finished creating epbins) + + std::vector zmins, zmaxs; + for (int zbin = 0; zbin < zBins_; zbin++) { + zmins.push_back(-1.0 * trkZMax_ + zStep_ * zbin); + zmaxs.push_back(-1.0 * trkZMax_ + zStep_ * zbin + zStep_ * 2); + } - // create vectors that hold data - std::vector> L1clusters; - L1clusters.reserve(phiBins_); - std::vector L2clusters; - - for (int i = 0; i < phiBins_; ++i) { - for (int j = 0; j < etaBins_; ++j) { - epbins_default[i][j].pTtot = 0; - epbins_default[i][j].used = false; - epbins_default[i][j].numtracks = 0; - epbins_default[i][j].numttrks = 0; - epbins_default[i][j].numtdtrks = 0; - epbins_default[i][j].numttdtrks = 0; - epbins_default[i][j].trackidx.clear(); + // create vectors that hold data + std::vector> L1clusters; + L1clusters.reserve(phiBins_); + std::vector L2clusters; + + for (int i = 0; i < phiBins_; ++i) { + for (int j = 0; j < etaBins_; ++j) { + epbins_default[i][j].pTtot = 0; + epbins_default[i][j].used = false; + epbins_default[i][j].numtracks = 0; + epbins_default[i][j].numttrks = 0; + epbins_default[i][j].numtdtrks = 0; + epbins_default[i][j].numttdtrks = 0; + epbins_default[i][j].trackidx.clear(); + } } - } - - for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) { - // initialize matrices - float zmin = zmins[zbin]; - float zmax = zmaxs[zbin]; - EtaPhiBin epbins[phiBins_][etaBins_]; - std::copy(&epbins_default[0][0], &epbins_default[0][0]+phiBins_*etaBins_,&epbins[0][0]); - - //clear containers - L1clusters.clear(); - L2clusters.clear(); - - // fill grid - for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { - float trkpt = L1TrkPtrs_[k]->momentum().perp(); - float trketa = L1TrkPtrs_[k]->momentum().eta(); - float trkphi = L1TrkPtrs_[k]->momentum().phi(); - float trkZ = L1TrkPtrs_[k]->z0(); - if ( zmax < trkZ) + + for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) { + // initialize matrices + float zmin = zmins[zbin]; + float zmax = zmaxs[zbin]; + EtaPhiBin epbins[phiBins_][etaBins_]; + std::copy(&epbins_default[0][0], &epbins_default[0][0] + phiBins_ * etaBins_, &epbins[0][0]); + + //clear containers + L1clusters.clear(); + L2clusters.clear(); + + // fill grid + for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { + float trkpt = L1TrkPtrs_[k]->momentum().perp(); + float trketa = L1TrkPtrs_[k]->momentum().eta(); + float trkphi = L1TrkPtrs_[k]->momentum().phi(); + float trkZ = L1TrkPtrs_[k]->z0(); + if (zmax < trkZ) continue; - if (zbin==0){ - if ( zmin > trkZ ) + if (zbin == 0) { + if (zmin > trkZ) continue; - } else{ - if ( zmin > trkZ || zmin ==trkZ ) + } else { + if (zmin >= trkZ) continue; - } - for (int i = 0; i < phiBins_; ++i) { - for (int j = 0; j < etaBins_; ++j) { - float eta_min=epbins[i][j].eta - etaStep_ / 2.0; //eta min - float eta_max=epbins[i][j].eta + etaStep_ / 2.0; //eta max - float phi_min=epbins[i][j].phi - phiStep_ / 2.0; //phi min - float phi_max=epbins[i][j].phi + phiStep_ / 2.0; //phi max - - if ( ( trketa < eta_min ) || ( trketa > eta_max ) || - ( trkphi < phi_min ) || ( trkphi > phi_max ) ) + } + for (int i = 0; i < phiBins_; ++i) { + for (int j = 0; j < etaBins_; ++j) { + float eta_min = epbins[i][j].eta - etaStep_ / 2.0; //eta min + float eta_max = epbins[i][j].eta + etaStep_ / 2.0; //eta max + float phi_min = epbins[i][j].phi - phiStep_ / 2.0; //phi min + float phi_max = epbins[i][j].phi + phiStep_ / 2.0; //phi max + + if ((trketa < eta_min) || (trketa > eta_max) || (trkphi < phi_min) || (trkphi > phi_max)) continue; - if (trkpt < trkPtMax_) - epbins[i][j].pTtot += trkpt; - else - epbins[i][j].pTtot += trkPtMax_; - epbins[i][j].numtdtrks += tdtrk_[k]; - epbins[i][j].trackidx.push_back(k); - ++epbins[i][j].numtracks; - } // for each phibin - } // for each phibin - } //end loop over tracks - - // first layer clustering - in eta for all phi bins - for (int phibin = 0; phibin < phiBins_; ++phibin) { - L1clusters.push_back( L1_clustering(epbins[phibin],etaBins_,etaStep_) ); - } - - //second layer clustering in phi for all eta clusters - L2clusters= L2_clustering(L1clusters, phiBins_, phiStep_, etaStep_); - - // sum all cluster pTs in this zbin to find max - float sum_pt = 0; - for (unsigned int k = 0; k < L2clusters.size(); ++k) { - if (L2clusters[k].pTtot > 50 && L2clusters[k].numtracks < lowpTJetMinTrackMultiplicity_) - continue; - if (L2clusters[k].pTtot > 100 && L2clusters[k].numtracks < highpTJetMinTrackMultiplicity_) + if (trkpt < trkPtMax_) + epbins[i][j].pTtot += trkpt; + else + epbins[i][j].pTtot += trkPtMax_; + epbins[i][j].numtdtrks += tdtrk_[k]; + epbins[i][j].trackidx.push_back(k); + ++epbins[i][j].numtracks; + } // for each phibin + } // for each phibin + } //end loop over tracks + + // first layer clustering - in eta for all phi bins + for (int phibin = 0; phibin < phiBins_; ++phibin) { + L1clusters.push_back(L1_clustering(epbins[phibin], etaBins_, etaStep_)); + } + + //second layer clustering in phi for all eta clusters + L2clusters = L2_clustering(L1clusters, phiBins_, phiStep_, etaStep_); + + // sum all cluster pTs in this zbin to find max + float sum_pt = 0; + for (unsigned int k = 0; k < L2clusters.size(); ++k) { + if (L2clusters[k].pTtot > lowpTJetThreshold_ && L2clusters[k].numtracks < lowpTJetMinTrackMultiplicity_) + continue; + if (L2clusters[k].pTtot > highpTJetThreshold_ && L2clusters[k].numtracks < highpTJetMinTrackMultiplicity_) + continue; + if (L2clusters[k].pTtot > minTrkJetpT_) + sum_pt += L2clusters[k].pTtot; + } + + if (sum_pt < mzb.ht) continue; - if (L2clusters[k].pTtot > minTrkJetpT_) - sum_pt += L2clusters[k].pTtot; - } - - if (sum_pt < mzb.ht) continue; - // if ht is larger than previous max, this is the new vertex zbin - mzb.ht=sum_pt; - mzb.znum = zbin; - mzb.clusters=L2clusters; - mzb.nclust = L2clusters.size(); - mzb.zbincenter = (zmin + zmax) / 2.0; - }//zbin loop - - - vector > L1TrackAssocJet; + // if ht is larger than previous max, this is the new vertex zbin + mzb.ht = sum_pt; + mzb.znum = zbin; + mzb.clusters = L2clusters; + mzb.nclust = L2clusters.size(); + mzb.zbincenter = (zmin + zmax) / 2.0; + } //zbin loop + + vector> L1TrackAssocJet; for (unsigned int j = 0; j < mzb.clusters.size(); ++j) { float jetEta = mzb.clusters[j].eta; float jetPhi = mzb.clusters[j].phi; @@ -327,22 +327,20 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { float jetPz = jetPt * sinh(jetEta); float jetP = jetPt * cosh(jetEta); int totalDisptrk = mzb.clusters[j].numtdtrks; - bool isDispJet=false; - if (totalDisptrk>nDisplacedTracks_ || totalDisptrk==nDisplacedTracks_) - isDispJet=true; + bool isDispJet = false; + if (totalDisptrk > nDisplacedTracks_ || totalDisptrk == nDisplacedTracks_) + isDispJet = true; math::XYZTLorentzVector jetP4(jetPx, jetPy, jetPz, jetP); L1TrackAssocJet.clear(); for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++) - L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]); + L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]); - TkJet trkJet(jetP4, L1TrackAssocJet, mzb.zbincenter, - mzb.clusters[j].numtracks, 0, totalDisptrk, 0, isDispJet); + TkJet trkJet(jetP4, L1TrackAssocJet, mzb.zbincenter, mzb.clusters[j].numtracks, 0, totalDisptrk, 0, isDispJet); if (!L1TrackAssocJet.empty()) L1L1TrackJetProducer->push_back(trkJet); } - if (displaced_) iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended"); @@ -355,23 +353,18 @@ void L1TrackJetProducer::beginStream(StreamID) {} void L1TrackJetProducer::endStream() {} -bool L1TrackJetProducer::trackQualityCuts( - int trk_nstub, float trk_chi2, float trk_bendchi2) { +bool L1TrackJetProducer::trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2) { bool PassQuality = false; - if (!displaced_){ - if ( trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ && - trk_chi2 < nStubs4PromptChi2_) - PassQuality = true; - if ( trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ && - trk_chi2 < nStubs5PromptChi2_) - PassQuality = true; - } else{ - if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ && - trk_chi2 < nStubs4DisplacedChi2_ ) - PassQuality = true; - if (trk_nstub >4 && trk_bendchi2 < nStubs5DisplacedBend_ && - trk_chi2 < nStubs5DisplacedChi2_ ) - PassQuality = true; + if (!displaced_) { + if (trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ && trk_chi2 < nStubs4PromptChi2_) // 4 stubs are the lowest track quality and have different cuts + PassQuality = true; + if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ && trk_chi2 < nStubs5PromptChi2_) // above 4 stubs diffent selection imposed (genrally looser) + PassQuality = true; + } else { + if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ && trk_chi2 < nStubs4DisplacedChi2_) // 4 stubs are the lowest track quality and have different cuts + PassQuality = true; + if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ && trk_chi2 < nStubs5DisplacedChi2_) // above 4 stubs diffent selection imposed (genrally looser) + PassQuality = true; } return PassQuality; } diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h index 841a64c995e3a..a70e96c443c94 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h @@ -26,27 +26,26 @@ struct MaxZBin { int nclust; //number of clusters in this bin float zbincenter; std::vector clusters; //list of all the clusters in this bin - float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored + float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored }; - -std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_ ) { - +std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) { std::vector clusters; // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used int nclust = 0; // get tracks in eta bins in increasing eta order for (int etabin = 0; etabin < etaBins_; ++etabin) { - float my_pt=0, previousbin_pt=0;//, nextbin_pt=0, next2bin_pt=0; - float nextbin_pt=0,nextbin2_pt=0; + float my_pt = 0, previousbin_pt = 0; //, nextbin_pt=0, next2bin_pt=0; + float nextbin_pt = 0, nextbin2_pt = 0; // skip (already) used tracks if (phislice[etabin].used) continue; my_pt = phislice[etabin].pTtot; - if (my_pt==0) continue; -//get previous bin pT + if (my_pt == 0) + continue; + //get previous bin pT if (etabin > 0 && !phislice[etabin - 1].used) previousbin_pt = phislice[etabin - 1].pTtot; @@ -60,30 +59,30 @@ std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float et // check if pT of current cluster is higher than neighbors if (my_pt < previousbin_pt || my_pt <= nextbin_pt) { // if unused pT in the left neighbor, spit it out as a cluster - if ( previousbin_pt > 0) { - clusters.push_back( phislice[etabin - 1] ); + if (previousbin_pt > 0) { + clusters.push_back(phislice[etabin - 1]); phislice[etabin - 1].used = true; nclust++; } - continue; //if it is not the local max pT skip - } + continue; //if it is not the local max pT skip + } // here reach only unused local max clusters - clusters.push_back( phislice[etabin]); - phislice[etabin].used = true; //if current bin a cluster + clusters.push_back(phislice[etabin]); + phislice[etabin].used = true; //if current bin a cluster if (previousbin_pt > 0) { clusters[nclust].pTtot += previousbin_pt; clusters[nclust].numtracks += phislice[etabin - 1].numtracks; clusters[nclust].numtdtrks += phislice[etabin - 1].numtdtrks; - for (unsigned int itrk=0; itrk= nextbin2_pt && nextbin_pt > 0) { clusters[nclust].pTtot += nextbin_pt; clusters[nclust].numtracks += phislice[etabin + 1].numtracks; clusters[nclust].numtdtrks += phislice[etabin + 1].numtdtrks; - for (unsigned int itrk=0; itrk L1_clustering(EtaPhiBin *phislice, int etaBins_, float et // Merge close-by clusters for (int m = 0; m < nclust - 1; ++m) { - if (fabs(clusters[m + 1].eta - clusters[m].eta) < 1.5 * etaStep_) { + if (std::abs(clusters[m + 1].eta - clusters[m].eta) < 1.5 * etaStep_) { if (clusters[m + 1].pTtot > clusters[m].pTtot) { clusters[m].eta = clusters[m + 1].eta; } clusters[m].pTtot += clusters[m + 1].pTtot; - clusters[m].numtracks += clusters[m + 1].numtracks; // total ntrk - clusters[m].numtdtrks += clusters[m + 1].numtdtrks; // total ndisp - for (unsigned int itrk=0; itrk trkidx){ - bin.pTtot+=pt; - bin.numtracks+=ntrk; - bin.numtdtrks+=ndtrk; - for (unsigned int itrk=0; itrk trkidx) { + bin.pTtot += pt; + bin.numtracks += ntrk; + bin.numtdtrks += ndtrk; + for (unsigned int itrk = 0; itrk < trkidx.size(); itrk++) + bin.trackidx.push_back(trkidx[itrk]); } - -float DPhi(float phi1, float phi2){ - float x = phi1-phi2; - if (x >= M_PI) - x -= 2*M_PI; - if (x < -1*M_PI) - x += 2*M_PI; - return x; +float DPhi(float phi1, float phi2) { + float x = phi1 - phi2; + if (x >= M_PI) + x -= 2 * M_PI; + if (x < -1 * M_PI) + x += 2 * M_PI; + return x; } -std::vector L2_clustering(std::vector> &L1clusters, int phiBins_, float phiStep_, float etaStep_ ) { - +std::vector L2_clustering(std::vector> &L1clusters, + int phiBins_, + float phiStep_, + float etaStep_) { std::vector clusters; for (int phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT - if (L1clusters[phibin].size()==0) continue; - - // sort L1 clusters max -> min - sort(L1clusters[phibin].begin(),L1clusters[phibin].end(),[] ( struct EtaPhiBin &a, struct EtaPhiBin &b){ return a.pTtot>b.pTtot; } ); - for (unsigned int imax = 0; imax trkidx1; - std::vector trkidx2; - clusters.push_back( L1clusters[phibin][imax]); - - L1clusters[phibin][imax].used = true; - - // if we are in the last phi bin, dont check phi+1 phi+2 - if (phibin == phiBins_ -1) - continue; - std::vector used_already; //keep phi+1 clusters that have been used - for (unsigned int icluster = 0; icluster min + sort(L1clusters[phibin].begin(), L1clusters[phibin].end(), [](struct EtaPhiBin &a, struct EtaPhiBin &b) { + return a.pTtot > b.pTtot; + }); + for (unsigned int imax = 0; imax < L1clusters[phibin].size(); ++imax) { + if (L1clusters[phibin][imax].used) + continue; + float pt_current = L1clusters[phibin][imax].pTtot; //current cluster (pt0) + float pt_next = 0; // next phi bin (pt1) + float pt_next2 = 0; // next to next phi bin2 (pt2) + int trk1 = 0; + int trk2 = 0; + int tdtrk1 = 0; + int tdtrk2 = 0; + std::vector trkidx1; + std::vector trkidx2; + clusters.push_back(L1clusters[phibin][imax]); + + L1clusters[phibin][imax].used = true; + + // if we are in the last phi bin, dont check phi+1 phi+2 + if (phibin == phiBins_ - 1) + continue; + std::vector used_already; //keep phi+1 clusters that have been used + for (unsigned int icluster = 0; icluster < L1clusters[phibin + 1].size(); ++icluster) { + if (L1clusters[phibin + 1][icluster].used) continue; - if (fabs(L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) + if (std::abs(L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) continue; - pt_next += L1clusters[phibin + 1][icluster].pTtot; - trk1 += L1clusters[phibin + 1][icluster].numtracks; - tdtrk1 += L1clusters[phibin + 1][icluster].numtdtrks; - for (unsigned int itrk=0; itrk used_already2; //keep used clusters in phi+2 - for (unsigned int icluster = 0; icluster used_already2; //keep used clusters in phi+2 + for (unsigned int icluster = 0; icluster < L1clusters[phibin + 2].size(); ++icluster) { + if (L1clusters[phibin + 2][icluster].used) continue; - if (fabs(L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) + if (std::abs(L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) continue; - pt_next2 += L1clusters[phibin + 2][icluster].pTtot; - trk2 += L1clusters[phibin + 2][icluster].numtracks; - tdtrk2 += L1clusters[phibin + 2][icluster].numtdtrks; - for (unsigned int itrk=0; itrk trkidx_both; - trkidx_both.reserve(trkidx1.size()+trkidx2.size()); - trkidx_both.insert(trkidx_both.end(),trkidx1.begin(),trkidx1.end()); - trkidx_both.insert(trkidx_both.end(),trkidx2.begin(),trkidx2.end()); - Fill_L2Cluster( clusters[clusters.size()-1], pt_next+pt_next2, trk1+trk2 , tdtrk1+tdtrk2, trkidx_both); - clusters[clusters.size()-1].phi = L1clusters[phibin + 1][used_already[0]].phi; - for(unsigned int iused: used_already) - L1clusters[phibin + 1][iused].used =true; - for(unsigned int iused: used_already2) - L1clusters[phibin + 2][iused].used =true; - } - } // for each L1 cluster + pt_next2 += L1clusters[phibin + 2][icluster].pTtot; + trk2 += L1clusters[phibin + 2][icluster].numtracks; + tdtrk2 += L1clusters[phibin + 2][icluster].numtdtrks; + for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][icluster].trackidx.size(); itrk++) + trkidx2.push_back(L1clusters[phibin + 2][icluster].trackidx[itrk]); + used_already2.push_back(icluster); + } + if (pt_next2 < pt_next) { + std::vector trkidx_both; + trkidx_both.reserve(trkidx1.size() + trkidx2.size()); + trkidx_both.insert(trkidx_both.end(), trkidx1.begin(), trkidx1.end()); + trkidx_both.insert(trkidx_both.end(), trkidx2.begin(), trkidx2.end()); + Fill_L2Cluster(clusters[clusters.size() - 1], pt_next + pt_next2, trk1 + trk2, tdtrk1 + tdtrk2, trkidx_both); + clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi; + for (unsigned int iused : used_already) + L1clusters[phibin + 1][iused].used = true; + for (unsigned int iused : used_already2) + L1clusters[phibin + 2][iused].used = true; + } + } // for each L1 cluster } // for each phibin - int nclust=clusters.size(); + int nclust = clusters.size(); // merge close-by clusters for (int m = 0; m < nclust - 1; ++m) { - for (int n = m + 1; n < nclust; ++n){ - if (clusters[n].eta != clusters[m].eta) - continue; - if ( fabs(DPhi(clusters[n].phi, clusters[m].phi))>1.5 * phiStep_) - continue; - + for (int n = m + 1; n < nclust; ++n) { + if (clusters[n].eta != clusters[m].eta) + continue; + if (std::abs(DPhi(clusters[n].phi, clusters[m].phi)) > 1.5 * phiStep_) + continue; + if (clusters[n].pTtot > clusters[m].pTtot) - clusters[m].phi = clusters[n].phi; + clusters[m].phi = clusters[n].phi; clusters[m].pTtot += clusters[n].pTtot; clusters[m].numtracks += clusters[n].numtracks; clusters[m].numtdtrks += clusters[n].numtdtrks; - for(unsigned int itrk=0; itrk4 from disp tag - lowpTJetMinTrackMultiplicity=cms.int32(2), #used only on zbin finding - highpTJetMinTrackMultiplicity=cms.int32(3), #used only on zbin finding + d0_cutNStubs4=cms.double(-1), # -1 excludes nstub=4 from disp tag process + d0_cutNStubs5=cms.double(0.22), # -1 excludes nstub>4 from disp tag process + lowpTJetMinTrackMultiplicity=cms.int32(2), #relevant only when N of z-bins >1; excludes from the HT calculation jets with low number of very energetic tracks; this cut selects the threshold on number of tracks + lowpTJetThreshold=cms.double(50.), # this threshold controls the pT of the jet + highpTJetMinTrackMultiplicity=cms.int32(3), #same as above for a different WP of tracks / pT + highpTJetThreshold=cms.double(100.), displaced=cms.bool(True), #Flag for displaced tracks nStubs4DisplacedChi2=cms.double(3.3), #Disp tracks selection [trk Date: Wed, 1 Jun 2022 16:59:31 +0200 Subject: [PATCH 3/8] changing to ES consumes for topology (cherry picked from commit 5c0a4a56c24cab2edc51eee6d8e80d9d3f72688f) (cherry picked from commit ef80134294fce48bff55af688596a00c31027ec5) --- .../plugins/L1TrackJetProducer.cc | 26 ++++++++++--------- .../plugins/L1TrackJetProducer.h | 16 ++++++------ 2 files changed, 22 insertions(+), 20 deletions(-) diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc index b0c2bba6fdeff..9a46fe3704bd3 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc @@ -63,6 +63,7 @@ class L1TrackJetProducer : public stream::EDProducer<> { // ----------member data --------------------------- + edm::ESGetToken tTopoToken_; const EDGetTokenT>> trackToken_; const edm::EDGetTokenT> PVtxToken_; vector> L1TrkPtrs_; @@ -99,7 +100,8 @@ class L1TrackJetProducer : public stream::EDProducer<> { }; L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) - : trackToken_(consumes>>(iConfig.getParameter("L1TrackInputTag"))), + : tTopoToken_(esConsumes(edm::ESInputTag("", ""))), + trackToken_(consumes>>(iConfig.getParameter("L1TrackInputTag"))), PVtxToken_(consumes>(iConfig.getParameter("L1PVertexCollection"))) { trkZMax_ = (float)iConfig.getParameter("trk_zMax"); trkPtMax_ = (float)iConfig.getParameter("trk_ptMax"); @@ -144,11 +146,7 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { unique_ptr L1L1TrackJetProducer(new TkJetCollection); // Read inputs - ESHandle tTopoHandle; - iSetup.get().get(tTopoHandle); - ESHandle tGeomHandle; - iSetup.get().get(tGeomHandle); - const TrackerTopology *const tTopo = tTopoHandle.product(); + const TrackerTopology& tTopo = iSetup.getData(tTopoToken_); edm::Handle>> TTTrackHandle; iEvent.getByToken(trackToken_, TTTrackHandle); @@ -172,8 +170,8 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { for (int istub = 0; istub < trk_nstubs; istub++) { // loop over the stubs DetId detId(trkPtr->getStubRefs().at(istub)->getDetId()); if (detId.det() == DetId::Detector::Tracker) { - if ((detId.subdetId() == StripSubdetector::TOB && tTopo->tobLayer(detId) <= 3) || - (detId.subdetId() == StripSubdetector::TID && tTopo->tidRing(detId) <= 9)) + if ((detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3) || + (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9)) trk_nPS++; } } @@ -356,14 +354,18 @@ void L1TrackJetProducer::endStream() {} bool L1TrackJetProducer::trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2) { bool PassQuality = false; if (!displaced_) { - if (trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ && trk_chi2 < nStubs4PromptChi2_) // 4 stubs are the lowest track quality and have different cuts + if (trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ && + trk_chi2 < nStubs4PromptChi2_) // 4 stubs are the lowest track quality and have different cuts PassQuality = true; - if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ && trk_chi2 < nStubs5PromptChi2_) // above 4 stubs diffent selection imposed (genrally looser) + if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ && + trk_chi2 < nStubs5PromptChi2_) // above 4 stubs diffent selection imposed (genrally looser) PassQuality = true; } else { - if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ && trk_chi2 < nStubs4DisplacedChi2_) // 4 stubs are the lowest track quality and have different cuts + if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ && + trk_chi2 < nStubs4DisplacedChi2_) // 4 stubs are the lowest track quality and have different cuts PassQuality = true; - if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ && trk_chi2 < nStubs5DisplacedChi2_) // above 4 stubs diffent selection imposed (genrally looser) + if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ && + trk_chi2 < nStubs5DisplacedChi2_) // above 4 stubs diffent selection imposed (genrally looser) PassQuality = true; } return PassQuality; diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h index a70e96c443c94..01b2b505944e8 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h @@ -29,7 +29,7 @@ struct MaxZBin { float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored }; -std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) { +inline std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) { std::vector clusters; // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used int nclust = 0; @@ -117,7 +117,7 @@ std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float et return clusters; } -void Fill_L2Cluster(EtaPhiBin &bin, float pt, int ntrk, int ndtrk, std::vector trkidx) { +inline void Fill_L2Cluster(EtaPhiBin &bin, float pt, int ntrk, int ndtrk, std::vector trkidx) { bin.pTtot += pt; bin.numtracks += ntrk; bin.numtdtrks += ndtrk; @@ -125,7 +125,7 @@ void Fill_L2Cluster(EtaPhiBin &bin, float pt, int ntrk, int ndtrk, std::vector= M_PI) x -= 2 * M_PI; @@ -134,13 +134,13 @@ float DPhi(float phi1, float phi2) { return x; } -std::vector L2_clustering(std::vector> &L1clusters, - int phiBins_, - float phiStep_, - float etaStep_) { +inline std::vector L2_clustering(std::vector> &L1clusters, + int phiBins_, + float phiStep_, + float etaStep_) { std::vector clusters; for (int phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT - if (L1clusters[phibin].size() == 0) + if (L1clusters[phibin].empty()) continue; // sort L1 clusters max -> min From 34106a8cd87ea03142cf52e687883aa97aad94b5 Mon Sep 17 00:00:00 2001 From: ccaillol Date: Sat, 11 Jun 2022 17:53:57 +0200 Subject: [PATCH 4/8] code format (cherry picked from commit 4fc043da37d4520d41222b646ccb0471a85ba172) --- L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc index 9a46fe3704bd3..83c0817efde6f 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc @@ -100,7 +100,7 @@ class L1TrackJetProducer : public stream::EDProducer<> { }; L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) - : tTopoToken_(esConsumes(edm::ESInputTag("", ""))), + : tTopoToken_(esConsumes(edm::ESInputTag("", ""))), trackToken_(consumes>>(iConfig.getParameter("L1TrackInputTag"))), PVtxToken_(consumes>(iConfig.getParameter("L1PVertexCollection"))) { trkZMax_ = (float)iConfig.getParameter("trk_zMax"); @@ -146,7 +146,7 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { unique_ptr L1L1TrackJetProducer(new TkJetCollection); // Read inputs - const TrackerTopology& tTopo = iSetup.getData(tTopoToken_); + const TrackerTopology &tTopo = iSetup.getData(tTopoToken_); edm::Handle>> TTTrackHandle; iEvent.getByToken(trackToken_, TTTrackHandle); From b686f3c58967394b3ea55a98ed2e962f489bdd3e Mon Sep 17 00:00:00 2001 From: ccaillol Date: Thu, 1 Sep 2022 15:34:33 +0200 Subject: [PATCH 5/8] address Adriano's comments (cherry picked from commit 626df8f92723e833e255a81ff9408b95e6cfd7c9) --- .../interface/L1TrackJetProducer.h | 246 +++++++++++++++-- .../plugins/L1TrackJetProducer.cc | 20 +- .../plugins/L1TrackJetProducer.h | 252 ------------------ 3 files changed, 236 insertions(+), 282 deletions(-) delete mode 100644 L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h diff --git a/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h b/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h index 57ac3a39953ae..7de68d459f20b 100644 --- a/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h +++ b/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h @@ -1,21 +1,16 @@ -#ifndef L1Trigger_L1TTrackMatch_L1TrackJetProducer_HH -#define L1Trigger_L1TTrackMatch_L1TrackJetProducer_HH +#pragma once #include #include -#include -#include -#include -#include //Each individual box in the eta and phi dimension. // Also used to store final cluster data for each zbin. struct EtaPhiBin { - float pTtot; - int numtracks; - int numttrks; - int numtdtrks; - int numttdtrks; - bool used; + float pTtot = 0; + int numtracks = 0; + int numttrks = 0; + int numtdtrks = 0; + int numttdtrks = 0; + bool used = false; float phi; //average phi value (halfway b/t min and max) float eta; //average eta value std::vector trackidx; @@ -26,7 +21,228 @@ struct MaxZBin { int znum; //Numbered from 0 to nzbins (16, 32, or 64) in order int nclust; //number of clusters in this bin float zbincenter; - EtaPhiBin *clusters; //list of all the clusters in this bin - float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored + std::vector clusters; //list of all the clusters in this bin + float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored }; -#endif + +inline std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) { + std::vector clusters; + // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used + int nclust = 0; + + // get tracks in eta bins in increasing eta order + for (int etabin = 0; etabin < etaBins_; ++etabin) { + float my_pt = 0, previousbin_pt = 0; //, nextbin_pt=0, next2bin_pt=0; + float nextbin_pt = 0, nextbin2_pt = 0; + + // skip (already) used tracks + if (phislice[etabin].used) + continue; + my_pt = phislice[etabin].pTtot; + if (my_pt == 0) + continue; + //get previous bin pT + if (etabin > 0 && !phislice[etabin - 1].used) + previousbin_pt = phislice[etabin - 1].pTtot; + + // get next bins pt + if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) { + nextbin_pt = phislice[etabin + 1].pTtot; + if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) { + nextbin2_pt = phislice[etabin + 2].pTtot; + } + } + // check if pT of current cluster is higher than neighbors + if (my_pt < previousbin_pt || my_pt <= nextbin_pt) { + // if unused pT in the left neighbor, spit it out as a cluster + if (previousbin_pt > 0) { + clusters.push_back(phislice[etabin - 1]); + phislice[etabin - 1].used = true; + nclust++; + } + continue; //if it is not the local max pT skip + } + // here reach only unused local max clusters + clusters.push_back(phislice[etabin]); + phislice[etabin].used = true; //if current bin a cluster + if (previousbin_pt > 0) { + clusters[nclust].pTtot += previousbin_pt; + clusters[nclust].numtracks += phislice[etabin - 1].numtracks; + clusters[nclust].numtdtrks += phislice[etabin - 1].numtdtrks; + for (unsigned int itrk = 0; itrk < phislice[etabin - 1].trackidx.size(); itrk++) + clusters[nclust].trackidx.push_back(phislice[etabin - 1].trackidx[itrk]); + } + + if (my_pt >= nextbin2_pt && nextbin_pt > 0) { + clusters[nclust].pTtot += nextbin_pt; + clusters[nclust].numtracks += phislice[etabin + 1].numtracks; + clusters[nclust].numtdtrks += phislice[etabin + 1].numtdtrks; + for (unsigned int itrk = 0; itrk < phislice[etabin + 1].trackidx.size(); itrk++) + clusters[nclust].trackidx.push_back(phislice[etabin + 1].trackidx[itrk]); + phislice[etabin + 1].used = true; + } + + nclust++; + + } // for each etabin + + // Merge close-by clusters + for (int m = 0; m < nclust - 1; ++m) { + if (std::abs(clusters[m + 1].eta - clusters[m].eta) < 1.5 * etaStep_) { + if (clusters[m + 1].pTtot > clusters[m].pTtot) { + clusters[m].eta = clusters[m + 1].eta; + } + clusters[m].pTtot += clusters[m + 1].pTtot; + clusters[m].numtracks += clusters[m + 1].numtracks; // total ntrk + clusters[m].numtdtrks += clusters[m + 1].numtdtrks; // total ndisp + for (unsigned int itrk = 0; itrk < clusters[m + 1].trackidx.size(); itrk++) + clusters[m].trackidx.push_back(clusters[m + 1].trackidx[itrk]); + + // if remove the merged cluster - all the others must be closer to 0 + for (int m1 = m + 1; m1 < nclust - 1; ++m1) { + clusters[m1] = clusters[m1 + 1]; + //clusters.erase(clusters.begin()+m1); + } + // clusters[m1] = clusters[m1 + 1]; + clusters.erase(clusters.begin() + nclust); + nclust--; + m = -1; + } // end if clusters neighbor in eta + } // end for (m) loop + + return clusters; +} + +inline void Fill_L2Cluster(EtaPhiBin &bin, float pt, int ntrk, int ndtrk, std::vector trkidx) { + bin.pTtot += pt; + bin.numtracks += ntrk; + bin.numtdtrks += ndtrk; + for (unsigned int itrk = 0; itrk < trkidx.size(); itrk++) + bin.trackidx.push_back(trkidx[itrk]); +} + +inline float DPhi(float phi1, float phi2) { + float x = phi1 - phi2; + if (x >= M_PI) + x -= 2 * M_PI; + if (x < -1 * M_PI) + x += 2 * M_PI; + return x; +} + +inline std::vector L2_clustering(std::vector> &L1clusters, + int phiBins_, + float phiStep_, + float etaStep_) { + std::vector clusters; + for (int phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT + if (L1clusters[phibin].empty()) + continue; + + // sort L1 clusters max -> min + sort(L1clusters[phibin].begin(), L1clusters[phibin].end(), [](struct EtaPhiBin &a, struct EtaPhiBin &b) { + return a.pTtot > b.pTtot; + }); + for (unsigned int imax = 0; imax < L1clusters[phibin].size(); ++imax) { + if (L1clusters[phibin][imax].used) + continue; + float pt_current = L1clusters[phibin][imax].pTtot; //current cluster (pt0) + float pt_next = 0; // next phi bin (pt1) + float pt_next2 = 0; // next to next phi bin2 (pt2) + int trk1 = 0; + int trk2 = 0; + int tdtrk1 = 0; + int tdtrk2 = 0; + std::vector trkidx1; + std::vector trkidx2; + clusters.push_back(L1clusters[phibin][imax]); + + L1clusters[phibin][imax].used = true; + + // if we are in the last phi bin, dont check phi+1 phi+2 + if (phibin == phiBins_ - 1) + continue; + std::vector used_already; //keep phi+1 clusters that have been used + for (unsigned int icluster = 0; icluster < L1clusters[phibin + 1].size(); ++icluster) { + if (L1clusters[phibin + 1][icluster].used) + continue; + if (std::abs(L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) + continue; + pt_next += L1clusters[phibin + 1][icluster].pTtot; + trk1 += L1clusters[phibin + 1][icluster].numtracks; + tdtrk1 += L1clusters[phibin + 1][icluster].numtdtrks; + for (unsigned int itrk = 0; itrk < L1clusters[phibin + 1][icluster].trackidx.size(); itrk++) + trkidx1.push_back(L1clusters[phibin + 1][icluster].trackidx[itrk]); + used_already.push_back(icluster); + } + + if (pt_next < pt_current) { // if pt1 used_already2; //keep used clusters in phi+2 + for (unsigned int icluster = 0; icluster < L1clusters[phibin + 2].size(); ++icluster) { + if (L1clusters[phibin + 2][icluster].used) + continue; + if (std::abs(L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) + continue; + pt_next2 += L1clusters[phibin + 2][icluster].pTtot; + trk2 += L1clusters[phibin + 2][icluster].numtracks; + tdtrk2 += L1clusters[phibin + 2][icluster].numtdtrks; + for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][icluster].trackidx.size(); itrk++) + trkidx2.push_back(L1clusters[phibin + 2][icluster].trackidx[itrk]); + used_already2.push_back(icluster); + } + if (pt_next2 < pt_next) { + std::vector trkidx_both; + trkidx_both.reserve(trkidx1.size() + trkidx2.size()); + trkidx_both.insert(trkidx_both.end(), trkidx1.begin(), trkidx1.end()); + trkidx_both.insert(trkidx_both.end(), trkidx2.begin(), trkidx2.end()); + Fill_L2Cluster(clusters[clusters.size() - 1], pt_next + pt_next2, trk1 + trk2, tdtrk1 + tdtrk2, trkidx_both); + clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi; + for (unsigned int iused : used_already) + L1clusters[phibin + 1][iused].used = true; + for (unsigned int iused : used_already2) + L1clusters[phibin + 2][iused].used = true; + } + } // for each L1 cluster + } // for each phibin + + int nclust = clusters.size(); + + // merge close-by clusters + for (int m = 0; m < nclust - 1; ++m) { + for (int n = m + 1; n < nclust; ++n) { + if (clusters[n].eta != clusters[m].eta) + continue; + if (std::abs(DPhi(clusters[n].phi, clusters[m].phi)) > 1.5 * phiStep_) + continue; + + if (clusters[n].pTtot > clusters[m].pTtot) + clusters[m].phi = clusters[n].phi; + + clusters[m].pTtot += clusters[n].pTtot; + clusters[m].numtracks += clusters[n].numtracks; + clusters[m].numtdtrks += clusters[n].numtdtrks; + for (unsigned int itrk = 0; itrk < clusters[n].trackidx.size(); itrk++) + clusters[m].trackidx.push_back(clusters[n].trackidx[itrk]); + for (int m1 = n; m1 < nclust - 1; ++m1) + clusters[m1] = clusters[m1 + 1]; + clusters.erase(clusters.begin() + nclust); + + nclust--; + m = -1; + } // end of n-loop + } // end of m-loop + return clusters; +} diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc index 83c0817efde6f..7186504bb1f83 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc @@ -31,7 +31,8 @@ #include "DataFormats/L1Trigger/interface/Vertex.h" -#include "L1TrackJetProducer.h" +#include "L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h" + #include "TH1D.h" #include "TH2D.h" #include @@ -165,6 +166,7 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { float trk_chi2dof = trkPtr->chi2Red(); float trk_d0 = trkPtr->d0(); float trk_bendchi2 = trkPtr->stubPtConsistency(); + float trk_z0 = trkPtr->z0(); int trk_nPS = 0; for (int istub = 0; istub < trk_nstubs; istub++) { // loop over the stubs @@ -180,9 +182,9 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { continue; if (!trackQualityCuts(trk_nstubs, trk_chi2dof, trk_bendchi2)) continue; - if (std::abs(PVz - trkPtr->z0()) > dzPVTrk_ && dzPVTrk_ > 0) + if (std::abs(PVz - trk_z0) > dzPVTrk_ && dzPVTrk_ > 0) continue; - if (std::abs(trkPtr->z0()) > trkZMax_) + if (std::abs(trk_z0) > trkZMax_) continue; if (std::abs(trkPtr->momentum().eta()) > trkEtaMax_) continue; @@ -227,18 +229,6 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { L1clusters.reserve(phiBins_); std::vector L2clusters; - for (int i = 0; i < phiBins_; ++i) { - for (int j = 0; j < etaBins_; ++j) { - epbins_default[i][j].pTtot = 0; - epbins_default[i][j].used = false; - epbins_default[i][j].numtracks = 0; - epbins_default[i][j].numttrks = 0; - epbins_default[i][j].numtdtrks = 0; - epbins_default[i][j].numttdtrks = 0; - epbins_default[i][j].trackidx.clear(); - } - } - for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) { // initialize matrices float zmin = zmins[zbin]; diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h deleted file mode 100644 index 01b2b505944e8..0000000000000 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.h +++ /dev/null @@ -1,252 +0,0 @@ -#pragma once -#include -#include -#include -#include -#include -using namespace std; - -//Each individual box in the eta and phi dimension. -// Also used to store final cluster data for each zbin. -struct EtaPhiBin { - float pTtot; - int numtracks; - int numttrks; - int numtdtrks; - int numttdtrks; - bool used; - float phi; //average phi value (halfway b/t min and max) - float eta; //average eta value - std::vector trackidx; -}; - -//store important information for plots -struct MaxZBin { - int znum; //Numbered from 0 to nzbins (16, 32, or 64) in order - int nclust; //number of clusters in this bin - float zbincenter; - std::vector clusters; //list of all the clusters in this bin - float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored -}; - -inline std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) { - std::vector clusters; - // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used - int nclust = 0; - - // get tracks in eta bins in increasing eta order - for (int etabin = 0; etabin < etaBins_; ++etabin) { - float my_pt = 0, previousbin_pt = 0; //, nextbin_pt=0, next2bin_pt=0; - float nextbin_pt = 0, nextbin2_pt = 0; - - // skip (already) used tracks - if (phislice[etabin].used) - continue; - my_pt = phislice[etabin].pTtot; - if (my_pt == 0) - continue; - //get previous bin pT - if (etabin > 0 && !phislice[etabin - 1].used) - previousbin_pt = phislice[etabin - 1].pTtot; - - // get next bins pt - if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) { - nextbin_pt = phislice[etabin + 1].pTtot; - if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) { - nextbin2_pt = phislice[etabin + 2].pTtot; - } - } - // check if pT of current cluster is higher than neighbors - if (my_pt < previousbin_pt || my_pt <= nextbin_pt) { - // if unused pT in the left neighbor, spit it out as a cluster - if (previousbin_pt > 0) { - clusters.push_back(phislice[etabin - 1]); - phislice[etabin - 1].used = true; - nclust++; - } - continue; //if it is not the local max pT skip - } - // here reach only unused local max clusters - clusters.push_back(phislice[etabin]); - phislice[etabin].used = true; //if current bin a cluster - if (previousbin_pt > 0) { - clusters[nclust].pTtot += previousbin_pt; - clusters[nclust].numtracks += phislice[etabin - 1].numtracks; - clusters[nclust].numtdtrks += phislice[etabin - 1].numtdtrks; - for (unsigned int itrk = 0; itrk < phislice[etabin - 1].trackidx.size(); itrk++) - clusters[nclust].trackidx.push_back(phislice[etabin - 1].trackidx[itrk]); - } - - if (my_pt >= nextbin2_pt && nextbin_pt > 0) { - clusters[nclust].pTtot += nextbin_pt; - clusters[nclust].numtracks += phislice[etabin + 1].numtracks; - clusters[nclust].numtdtrks += phislice[etabin + 1].numtdtrks; - for (unsigned int itrk = 0; itrk < phislice[etabin + 1].trackidx.size(); itrk++) - clusters[nclust].trackidx.push_back(phislice[etabin + 1].trackidx[itrk]); - phislice[etabin + 1].used = true; - } - - nclust++; - - } // for each etabin - - // Merge close-by clusters - for (int m = 0; m < nclust - 1; ++m) { - if (std::abs(clusters[m + 1].eta - clusters[m].eta) < 1.5 * etaStep_) { - if (clusters[m + 1].pTtot > clusters[m].pTtot) { - clusters[m].eta = clusters[m + 1].eta; - } - clusters[m].pTtot += clusters[m + 1].pTtot; - clusters[m].numtracks += clusters[m + 1].numtracks; // total ntrk - clusters[m].numtdtrks += clusters[m + 1].numtdtrks; // total ndisp - for (unsigned int itrk = 0; itrk < clusters[m + 1].trackidx.size(); itrk++) - clusters[m].trackidx.push_back(clusters[m + 1].trackidx[itrk]); - - // if remove the merged cluster - all the others must be closer to 0 - for (int m1 = m + 1; m1 < nclust - 1; ++m1) { - clusters[m1] = clusters[m1 + 1]; - //clusters.erase(clusters.begin()+m1); - } - // clusters[m1] = clusters[m1 + 1]; - clusters.erase(clusters.begin() + nclust); - nclust--; - m = -1; - } // end if clusters neighbor in eta - } // end for (m) loop - - return clusters; -} - -inline void Fill_L2Cluster(EtaPhiBin &bin, float pt, int ntrk, int ndtrk, std::vector trkidx) { - bin.pTtot += pt; - bin.numtracks += ntrk; - bin.numtdtrks += ndtrk; - for (unsigned int itrk = 0; itrk < trkidx.size(); itrk++) - bin.trackidx.push_back(trkidx[itrk]); -} - -inline float DPhi(float phi1, float phi2) { - float x = phi1 - phi2; - if (x >= M_PI) - x -= 2 * M_PI; - if (x < -1 * M_PI) - x += 2 * M_PI; - return x; -} - -inline std::vector L2_clustering(std::vector> &L1clusters, - int phiBins_, - float phiStep_, - float etaStep_) { - std::vector clusters; - for (int phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT - if (L1clusters[phibin].empty()) - continue; - - // sort L1 clusters max -> min - sort(L1clusters[phibin].begin(), L1clusters[phibin].end(), [](struct EtaPhiBin &a, struct EtaPhiBin &b) { - return a.pTtot > b.pTtot; - }); - for (unsigned int imax = 0; imax < L1clusters[phibin].size(); ++imax) { - if (L1clusters[phibin][imax].used) - continue; - float pt_current = L1clusters[phibin][imax].pTtot; //current cluster (pt0) - float pt_next = 0; // next phi bin (pt1) - float pt_next2 = 0; // next to next phi bin2 (pt2) - int trk1 = 0; - int trk2 = 0; - int tdtrk1 = 0; - int tdtrk2 = 0; - std::vector trkidx1; - std::vector trkidx2; - clusters.push_back(L1clusters[phibin][imax]); - - L1clusters[phibin][imax].used = true; - - // if we are in the last phi bin, dont check phi+1 phi+2 - if (phibin == phiBins_ - 1) - continue; - std::vector used_already; //keep phi+1 clusters that have been used - for (unsigned int icluster = 0; icluster < L1clusters[phibin + 1].size(); ++icluster) { - if (L1clusters[phibin + 1][icluster].used) - continue; - if (std::abs(L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) - continue; - pt_next += L1clusters[phibin + 1][icluster].pTtot; - trk1 += L1clusters[phibin + 1][icluster].numtracks; - tdtrk1 += L1clusters[phibin + 1][icluster].numtdtrks; - for (unsigned int itrk = 0; itrk < L1clusters[phibin + 1][icluster].trackidx.size(); itrk++) - trkidx1.push_back(L1clusters[phibin + 1][icluster].trackidx[itrk]); - used_already.push_back(icluster); - } - - if (pt_next < pt_current) { // if pt1 used_already2; //keep used clusters in phi+2 - for (unsigned int icluster = 0; icluster < L1clusters[phibin + 2].size(); ++icluster) { - if (L1clusters[phibin + 2][icluster].used) - continue; - if (std::abs(L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) - continue; - pt_next2 += L1clusters[phibin + 2][icluster].pTtot; - trk2 += L1clusters[phibin + 2][icluster].numtracks; - tdtrk2 += L1clusters[phibin + 2][icluster].numtdtrks; - for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][icluster].trackidx.size(); itrk++) - trkidx2.push_back(L1clusters[phibin + 2][icluster].trackidx[itrk]); - used_already2.push_back(icluster); - } - if (pt_next2 < pt_next) { - std::vector trkidx_both; - trkidx_both.reserve(trkidx1.size() + trkidx2.size()); - trkidx_both.insert(trkidx_both.end(), trkidx1.begin(), trkidx1.end()); - trkidx_both.insert(trkidx_both.end(), trkidx2.begin(), trkidx2.end()); - Fill_L2Cluster(clusters[clusters.size() - 1], pt_next + pt_next2, trk1 + trk2, tdtrk1 + tdtrk2, trkidx_both); - clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi; - for (unsigned int iused : used_already) - L1clusters[phibin + 1][iused].used = true; - for (unsigned int iused : used_already2) - L1clusters[phibin + 2][iused].used = true; - } - } // for each L1 cluster - } // for each phibin - - int nclust = clusters.size(); - - // merge close-by clusters - for (int m = 0; m < nclust - 1; ++m) { - for (int n = m + 1; n < nclust; ++n) { - if (clusters[n].eta != clusters[m].eta) - continue; - if (std::abs(DPhi(clusters[n].phi, clusters[m].phi)) > 1.5 * phiStep_) - continue; - - if (clusters[n].pTtot > clusters[m].pTtot) - clusters[m].phi = clusters[n].phi; - - clusters[m].pTtot += clusters[n].pTtot; - clusters[m].numtracks += clusters[n].numtracks; - clusters[m].numtdtrks += clusters[n].numtdtrks; - for (unsigned int itrk = 0; itrk < clusters[n].trackidx.size(); itrk++) - clusters[m].trackidx.push_back(clusters[n].trackidx[itrk]); - for (int m1 = n; m1 < nclust - 1; ++m1) - clusters[m1] = clusters[m1 + 1]; - clusters.erase(clusters.begin() + nclust); - - nclust--; - m = -1; - } // end of n-loop - } // end of m-loop - return clusters; -} From 1098a8f81483378137c08fbe9b52d5ea5e2896a1 Mon Sep 17 00:00:00 2001 From: ccaillol Date: Thu, 15 Sep 2022 12:34:22 +0200 Subject: [PATCH 6/8] Adrianos comments (cherry picked from commit 23cd32cd4107928c1f5f5addd36f31c9c4e300cc) --- .../interface/L1TrackJetProducer.h | 1 + .../plugins/L1TrackJetProducer.cc | 114 +++++++++--------- 2 files changed, 58 insertions(+), 57 deletions(-) diff --git a/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h b/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h index 7de68d459f20b..b8a3ed94703a1 100644 --- a/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h +++ b/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h @@ -1,6 +1,7 @@ #pragma once #include #include +#include //Each individual box in the eta and phi dimension. // Also used to store final cluster data for each zbin. diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc index 7186504bb1f83..d645dffcaa4a7 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc @@ -64,73 +64,73 @@ class L1TrackJetProducer : public stream::EDProducer<> { // ----------member data --------------------------- - edm::ESGetToken tTopoToken_; - const EDGetTokenT>> trackToken_; - const edm::EDGetTokenT> PVtxToken_; vector> L1TrkPtrs_; vector tdtrk_; - float trkZMax_; - float trkPtMax_; - float trkPtMin_; - float trkEtaMax_; - float nStubs4PromptChi2_; - float nStubs5PromptChi2_; - float nStubs4PromptBend_; - float nStubs5PromptBend_; - int trkNPSStubMin_; - int lowpTJetMinTrackMultiplicity_; - float lowpTJetThreshold_; - int highpTJetMinTrackMultiplicity_; - float highpTJetThreshold_; - int zBins_; - int etaBins_; - int phiBins_; - double minTrkJetpT_; + const float trkZMax_; + const float trkPtMax_; + const float trkPtMin_; + const float trkEtaMax_; + const float nStubs4PromptChi2_; + const float nStubs5PromptChi2_; + const float nStubs4PromptBend_; + const float nStubs5PromptBend_; + const int trkNPSStubMin_; + const int lowpTJetMinTrackMultiplicity_; + const float lowpTJetThreshold_; + const int highpTJetMinTrackMultiplicity_; + const float highpTJetThreshold_; + const int zBins_; + const int etaBins_; + const int phiBins_; + const double minTrkJetpT_; float zStep_; float etaStep_; float phiStep_; - bool displaced_; - float d0CutNStubs4_; - float d0CutNStubs5_; - float nStubs4DisplacedChi2_; - float nStubs5DisplacedChi2_; - float nStubs4DisplacedBend_; - float nStubs5DisplacedBend_; - int nDisplacedTracks_; - float dzPVTrk_; + const bool displaced_; + const float d0CutNStubs4_; + const float d0CutNStubs5_; + const float nStubs4DisplacedChi2_; + const float nStubs5DisplacedChi2_; + const float nStubs4DisplacedBend_; + const float nStubs5DisplacedBend_; + const int nDisplacedTracks_; + const float dzPVTrk_; + + edm::ESGetToken tTopoToken_; + const EDGetTokenT>> trackToken_; + const edm::EDGetTokenT> PVtxToken_; }; L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) - : tTopoToken_(esConsumes(edm::ESInputTag("", ""))), + : trkZMax_((float)iConfig.getParameter("trk_zMax")), + trkPtMax_((float)iConfig.getParameter("trk_ptMax")), + trkPtMin_((float)iConfig.getParameter("trk_ptMin")), + trkEtaMax_((float)iConfig.getParameter("trk_etaMax")), + nStubs4PromptChi2_((float)iConfig.getParameter("nStubs4PromptChi2")), + nStubs5PromptChi2_((float)iConfig.getParameter("nStubs5PromptChi2")), + nStubs4PromptBend_((float)iConfig.getParameter("nStubs4PromptBend")), + nStubs5PromptBend_((float)iConfig.getParameter("nStubs5PromptBend")), + trkNPSStubMin_((int)iConfig.getParameter("trk_nPSStubMin")), + lowpTJetMinTrackMultiplicity_((int)iConfig.getParameter("lowpTJetMinTrackMultiplicity")), + lowpTJetThreshold_((float)iConfig.getParameter("lowpTJetThreshold")), + highpTJetMinTrackMultiplicity_((int)iConfig.getParameter("highpTJetMinTrackMultiplicity")), + highpTJetThreshold_((float)iConfig.getParameter("highpTJetThreshold")), + zBins_((int)iConfig.getParameter("zBins")), + etaBins_((int)iConfig.getParameter("etaBins")), + phiBins_((int)iConfig.getParameter("phiBins")), + minTrkJetpT_(iConfig.getParameter("minTrkJetpT")), + displaced_(iConfig.getParameter("displaced")), + d0CutNStubs4_((float)iConfig.getParameter("d0_cutNStubs4")), + d0CutNStubs5_((float)iConfig.getParameter("d0_cutNStubs5")), + nStubs4DisplacedChi2_((float)iConfig.getParameter("nStubs4DisplacedChi2")), + nStubs5DisplacedChi2_((float)iConfig.getParameter("nStubs5DisplacedChi2")), + nStubs4DisplacedBend_((float)iConfig.getParameter("nStubs4DisplacedBend")), + nStubs5DisplacedBend_((float)iConfig.getParameter("nStubs5DisplacedBend")), + nDisplacedTracks_((int)iConfig.getParameter("nDisplacedTracks")), + dzPVTrk_((float)iConfig.getParameter("MaxDzTrackPV")), + tTopoToken_(esConsumes(edm::ESInputTag("", ""))), trackToken_(consumes>>(iConfig.getParameter("L1TrackInputTag"))), PVtxToken_(consumes>(iConfig.getParameter("L1PVertexCollection"))) { - trkZMax_ = (float)iConfig.getParameter("trk_zMax"); - trkPtMax_ = (float)iConfig.getParameter("trk_ptMax"); - trkPtMin_ = (float)iConfig.getParameter("trk_ptMin"); - trkEtaMax_ = (float)iConfig.getParameter("trk_etaMax"); - nStubs4PromptChi2_ = (float)iConfig.getParameter("nStubs4PromptChi2"); - nStubs5PromptChi2_ = (float)iConfig.getParameter("nStubs5PromptChi2"); - nStubs4PromptBend_ = (float)iConfig.getParameter("nStubs4PromptBend"); - nStubs5PromptBend_ = (float)iConfig.getParameter("nStubs5PromptBend"); - trkNPSStubMin_ = (int)iConfig.getParameter("trk_nPSStubMin"); - minTrkJetpT_ = iConfig.getParameter("minTrkJetpT"); - etaBins_ = (int)iConfig.getParameter("etaBins"); - phiBins_ = (int)iConfig.getParameter("phiBins"); - zBins_ = (int)iConfig.getParameter("zBins"); - d0CutNStubs4_ = (float)iConfig.getParameter("d0_cutNStubs4"); - d0CutNStubs5_ = (float)iConfig.getParameter("d0_cutNStubs5"); - lowpTJetMinTrackMultiplicity_ = (int)iConfig.getParameter("lowpTJetMinTrackMultiplicity"); - lowpTJetThreshold_ = (float)iConfig.getParameter("lowpTJetThreshold"); - highpTJetMinTrackMultiplicity_ = (int)iConfig.getParameter("highpTJetMinTrackMultiplicity"); - highpTJetThreshold_ = (float)iConfig.getParameter("highpTJetThreshold"); - displaced_ = iConfig.getParameter("displaced"); - nStubs4DisplacedChi2_ = (float)iConfig.getParameter("nStubs4DisplacedChi2"); - nStubs5DisplacedChi2_ = (float)iConfig.getParameter("nStubs5DisplacedChi2"); - nStubs4DisplacedBend_ = (float)iConfig.getParameter("nStubs4DisplacedBend"); - nStubs5DisplacedBend_ = (float)iConfig.getParameter("nStubs5DisplacedBend"); - nDisplacedTracks_ = (int)iConfig.getParameter("nDisplacedTracks"); - dzPVTrk_ = (float)iConfig.getParameter("MaxDzTrackPV"); - zStep_ = 2.0 * trkZMax_ / (zBins_ + 1); // added +1 in denom etaStep_ = 2.0 * trkEtaMax_ / etaBins_; //etaStep is the width of an etabin phiStep_ = 2 * M_PI / phiBins_; ////phiStep is the width of a phibin From 7f6d1ce2248f3446e16012d3f6d310108275ca55 Mon Sep 17 00:00:00 2001 From: ccaillol Date: Thu, 15 Sep 2022 14:06:21 +0200 Subject: [PATCH 7/8] fix header consistency (cherry picked from commit 2854e1c70cd9100809a06737ff092f6c06fac797) --- L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h | 1 + 1 file changed, 1 insertion(+) diff --git a/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h b/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h index b8a3ed94703a1..0ce88d6472ed6 100644 --- a/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h +++ b/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h @@ -2,6 +2,7 @@ #include #include #include +#include //Each individual box in the eta and phi dimension. // Also used to store final cluster data for each zbin. From 1d7747eb6ddf9a41f12555453edc1ad80ced1521 Mon Sep 17 00:00:00 2001 From: ccaillol Date: Mon, 19 Sep 2022 16:25:54 +0200 Subject: [PATCH 8/8] andrea's comments (cherry picked from commit 7964d2caa2c091ffc7b9b638df2d47e379dfb8b1) --- .../{L1TrackJetProducer.h => L1Clustering.h} | 4 +- .../plugins/L1TrackJetProducer.cc | 63 ++++++++++++------- .../L1TTrackMatch/python/l1tTrackJets_cfi.py | 40 ++++-------- 3 files changed, 55 insertions(+), 52 deletions(-) rename L1Trigger/L1TTrackMatch/interface/{L1TrackJetProducer.h => L1Clustering.h} (98%) diff --git a/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h b/L1Trigger/L1TTrackMatch/interface/L1Clustering.h similarity index 98% rename from L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h rename to L1Trigger/L1TTrackMatch/interface/L1Clustering.h index 0ce88d6472ed6..f97a11e8431ae 100644 --- a/L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h +++ b/L1Trigger/L1TTrackMatch/interface/L1Clustering.h @@ -1,4 +1,5 @@ -#pragma once +#ifndef L1Trigger_L1TTrackMatch_L1Clustering_HH +#define L1Trigger_L1TTrackMatch_L1Clustering_HH #include #include #include @@ -248,3 +249,4 @@ inline std::vector L2_clustering(std::vector> } // end of m-loop return clusters; } +#endif diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc index d645dffcaa4a7..32c40ba9d37a4 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc @@ -19,7 +19,6 @@ #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -31,7 +30,7 @@ #include "DataFormats/L1Trigger/interface/Vertex.h" -#include "L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h" +#include "L1Trigger/L1TTrackMatch/interface/L1Clustering.h" #include "TH1D.h" #include "TH2D.h" @@ -207,9 +206,8 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { EtaPhiBin epbins_default[phiBins_][etaBins_]; // create grid of phiBins float phi = -1.0 * M_PI; - float eta = -1.0 * trkEtaMax_; for (int i = 0; i < phiBins_; ++i) { - eta = -1.0 * trkEtaMax_; + float eta = -1.0 * trkEtaMax_; for (int j = 0; j < etaBins_; ++j) { epbins_default[i][j].phi = (phi + (phi + phiStep_)) / 2.0; // phimin=phi; phimax= phimin+step epbins_default[i][j].eta = (eta + (eta + etaStep_)) / 2.0; // phimin=phi; phimax phimin+step @@ -242,19 +240,16 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { // fill grid for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { - float trkpt = L1TrkPtrs_[k]->momentum().perp(); - float trketa = L1TrkPtrs_[k]->momentum().eta(); - float trkphi = L1TrkPtrs_[k]->momentum().phi(); float trkZ = L1TrkPtrs_[k]->z0(); if (zmax < trkZ) continue; - if (zbin == 0) { - if (zmin > trkZ) - continue; - } else { - if (zmin >= trkZ) - continue; - } + if (zmin > trkZ) + continue; + if (zbin == 0 && zmin == trkZ) + continue; + float trkpt = L1TrkPtrs_[k]->momentum().perp(); + float trketa = L1TrkPtrs_[k]->momentum().eta(); + float trkphi = L1TrkPtrs_[k]->momentum().phi(); for (int i = 0; i < phiBins_; ++i) { for (int j = 0; j < etaBins_; ++j) { float eta_min = epbins[i][j].eta - etaStep_ / 2.0; //eta min @@ -315,9 +310,7 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { float jetPz = jetPt * sinh(jetEta); float jetP = jetPt * cosh(jetEta); int totalDisptrk = mzb.clusters[j].numtdtrks; - bool isDispJet = false; - if (totalDisptrk > nDisplacedTracks_ || totalDisptrk == nDisplacedTracks_) - isDispJet = true; + bool isDispJet = (totalDisptrk > nDisplacedTracks_ || totalDisptrk == nDisplacedTracks_); math::XYZTLorentzVector jetP4(jetPx, jetPy, jetPz, jetP); L1TrackAssocJet.clear(); @@ -362,11 +355,37 @@ bool L1TrackJetProducer::trackQualityCuts(int trk_nstub, float trk_chi2, float t } void L1TrackJetProducer::fillDescriptions(ConfigurationDescriptions &descriptions) { - //The following says we do not know what parameters are allowed so do no validation - // Please change this to state exactly what you do use, even if it is no parameters - ParameterSetDescription desc; - desc.setUnknown(); - descriptions.addDefault(desc); + // l1tTrackJets + edm::ParameterSetDescription desc; + desc.add("L1TrackInputTag", edm::InputTag("l1tTTTracksFromTrackletEmulation", "Level1TTTracks")); + desc.add("L1PVertexCollection", edm::InputTag("l1tVertexProducer", "l1vertices")); + desc.add("MaxDzTrackPV", 1.0); + desc.add("trk_zMax", 15.0); + desc.add("trk_ptMax", 200.0); + desc.add("trk_ptMin", 3.0); + desc.add("trk_etaMax", 2.4); + desc.add("nStubs4PromptChi2", 5.0); + desc.add("nStubs4PromptBend", 1.7); + desc.add("nStubs5PromptChi2", 2.75); + desc.add("nStubs5PromptBend", 3.5); + desc.add("trk_nPSStubMin", -1); + desc.add("minTrkJetpT", -1.0); + desc.add("etaBins", 24); + desc.add("phiBins", 27); + desc.add("zBins", 1); + desc.add("d0_cutNStubs4", -1); + desc.add("d0_cutNStubs5", -1); + desc.add("lowpTJetMinTrackMultiplicity", 2); + desc.add("lowpTJetThreshold", 50.0); + desc.add("highpTJetMinTrackMultiplicity", 3); + desc.add("highpTJetThreshold", 100.0); + desc.add("displaced", false); + desc.add("nStubs4DisplacedChi2", 5.0); + desc.add("nStubs4DisplacedBend", 1.7); + desc.add("nStubs5DisplacedChi2", 2.75); + desc.add("nStubs5DisplacedBend", 3.5); + desc.add("nDisplacedTracks", 2); + descriptions.add("l1tTrackJets", desc); } //define this as a plug-in diff --git a/L1Trigger/L1TTrackMatch/python/l1tTrackJets_cfi.py b/L1Trigger/L1TTrackMatch/python/l1tTrackJets_cfi.py index 53697ab50b2b0..5a10bb03bc44b 100644 --- a/L1Trigger/L1TTrackMatch/python/l1tTrackJets_cfi.py +++ b/L1Trigger/L1TTrackMatch/python/l1tTrackJets_cfi.py @@ -31,35 +31,17 @@ nDisplacedTracks=cms.int32(2) ) -l1tTrackJetsExtended = cms.EDProducer('L1TrackJetProducer', - L1TrackInputTag= cms.InputTag("l1tTTTracksFromExtendedTrackletEmulation", "Level1TTTracks"), - L1PVertexCollection = cms.InputTag("l1tVertexProducer", "l1vertices"), - MaxDzTrackPV = cms.double( 5.0 ), # tracks with dz(trk,PV)>cut excluded - trk_zMax = cms.double (15.) , # max track z - trk_ptMax = cms.double(200.), # maxi track pT before saturation - trk_ptMin = cms.double(3.0), # min track pt - trk_etaMax = cms.double(2.4), # max track eta - nStubs4PromptChi2=cms.double(5.0), #Prompt track quality flags for loose/tight - nStubs4PromptBend=cms.double(1.7), - nStubs5PromptChi2=cms.double(2.75), - nStubs5PromptBend=cms.double(3.5), - trk_nPSStubMin=cms.int32(-1), # min # PS stubs, -1 means no cut - minTrkJetpT=cms.double(5.), # min track jet pt to be considered for most energetic zbin finding - etaBins=cms.int32(24), #number of eta bins - phiBins=cms.int32(27), #number of phi bins - zBins=cms.int32(1), #number of z bins - d0_cutNStubs4=cms.double(-1), # -1 excludes nstub=4 from disp tag process - d0_cutNStubs5=cms.double(0.22), # -1 excludes nstub>4 from disp tag process - lowpTJetMinTrackMultiplicity=cms.int32(2), #relevant only when N of z-bins >1; excludes from the HT calculation jets with low number of very energetic tracks; this cut selects the threshold on number of tracks - lowpTJetThreshold=cms.double(50.), # this threshold controls the pT of the jet - highpTJetMinTrackMultiplicity=cms.int32(3), #same as above for a different WP of tracks / pT - highpTJetThreshold=cms.double(100.), - displaced=cms.bool(True), #Flag for displaced tracks - nStubs4DisplacedChi2=cms.double(3.3), #Disp tracks selection [trkcut excluded + minTrkJetpT= 5., # min track jet pt to be considered for most energetic zbin finding + d0_cutNStubs5= 0.22, # -1 excludes nstub>4 from disp tag process + displaced=True, #Flag for displaced tracks + nStubs4DisplacedChi2= 3.3, #Disp tracks selection [trk