Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-05 02:13:32

0001 /* Software to emulate the hardware 2-layer jet-finding algorithm (fixed-point). *Layers 1 and 2*
0002  *
0003  * Created based on Rishi Patel's L1TrackJetProducer.cc file.
0004  *      Authors: Samuel Edwin Leigh, Tyler Wu
0005  *               Rutgers, the State University of New Jersey
0006  *      Modifications: Georgios Karathanasis
0007  *                     georgios.karathanasis@cern.ch, CU Boulder
0008  *
0009  * Last update: 19-9-2022 (by GK)     
0010  */
0011 
0012 //Holds data from tracks, converted from their integer versions.
0013 
0014 // system include files
0015 
0016 #include "DataFormats/Common/interface/Ref.h"
0017 #include "DataFormats/L1TCorrelator/interface/TkJet.h"
0018 #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h"
0019 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0020 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
0021 #include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h"
0022 #include "DataFormats/L1Trigger/interface/TkJetWord.h"
0023 #include "DataFormats/L1Trigger/interface/VertexWord.h"
0024 // user include files
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/StreamID.h"
0032 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0033 #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
0034 #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h"
0035 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0036 
0037 #include "DataFormats/L1Trigger/interface/Vertex.h"
0038 
0039 #include <memory>
0040 #include <iostream>
0041 #include <fstream>
0042 #include <cmath>
0043 #include <cstdlib>
0044 #include <string>
0045 #include <cstdlib>
0046 #include "TH1D.h"
0047 #include "TH2D.h"
0048 #include <TMath.h>
0049 #include <ap_int.h>
0050 #include "L1Trigger/L1TTrackMatch/interface/L1TrackJetEmulationProducer.h"
0051 
0052 using namespace std;
0053 using namespace edm;
0054 class L1TrackJetEmulationProducer : public stream::EDProducer<> {
0055 public:
0056   explicit L1TrackJetEmulationProducer(const ParameterSet &);
0057   ~L1TrackJetEmulationProducer() override;
0058   typedef TTTrack<Ref_Phase2TrackerDigi_> L1TTTrackType;
0059   typedef vector<L1TTTrackType> L1TTTrackCollectionType;
0060 
0061   static void fillDescriptions(ConfigurationDescriptions &descriptions);
0062   bool trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2);
0063   void L2_cluster(vector<Ptr<L1TTTrackType>> L1TrkPtrs_, TrackJetEmulationMaxZBin &mzb);
0064   virtual TrackJetEmulationEtaPhiBin *L1_cluster(TrackJetEmulationEtaPhiBin *phislice);
0065 
0066 private:
0067   void beginStream(StreamID) override;
0068   void produce(Event &, const EventSetup &) override;
0069   void endStream() override;
0070 
0071   // ----------member data ---------------------------
0072 
0073   vector<Ptr<L1TTTrackType>> L1TrkPtrs_;
0074   vector<int> zBinCount_;
0075   vector<int> tdtrk_;
0076   const float MaxDzTrackPV;
0077   const float trkZMax_;
0078   const float trkPtMax_;
0079   const float trkPtMin_;
0080   const float trkEtaMax_;
0081   const float trkChi2dofMax_;
0082   const float trkBendChi2Max_;
0083   const int trkNPSStubMin_;
0084   const double minTrkJetpT_;
0085   int etaBins_;
0086   int phiBins_;
0087   int zBins_;
0088   float d0CutNStubs4_;
0089   float d0CutNStubs5_;
0090   int lowpTJetMinTrackMultiplicity_;
0091   float lowpTJetMinpT_;
0092   int highpTJetMinTrackMultiplicity_;
0093   float highpTJetMinpT_;
0094   bool displaced_;
0095   float nStubs4DisplacedChi2_;
0096   float nStubs5DisplacedChi2_;
0097   float nStubs4DisplacedBend_;
0098   float nStubs5DisplacedBend_;
0099   int nDisplacedTracks_;
0100 
0101   float PVz;
0102   z0_intern zStep_;
0103   glbeta_intern etaStep_;
0104   glbphi_intern phiStep_;
0105 
0106   const EDGetTokenT<vector<TTTrack<Ref_Phase2TrackerDigi_>>> trackToken_;
0107   const EDGetTokenT<l1t::VertexWordCollection> PVtxToken_;
0108   ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0109 };
0110 
0111 L1TrackJetEmulationProducer::L1TrackJetEmulationProducer(const ParameterSet &iConfig)
0112     : MaxDzTrackPV((float)iConfig.getParameter<double>("MaxDzTrackPV")),
0113       trkZMax_((float)iConfig.getParameter<double>("trk_zMax")),
0114       trkPtMax_((float)iConfig.getParameter<double>("trk_ptMax")),
0115       trkPtMin_((float)iConfig.getParameter<double>("trk_ptMin")),
0116       trkEtaMax_((float)iConfig.getParameter<double>("trk_etaMax")),
0117       trkChi2dofMax_((float)iConfig.getParameter<double>("trk_chi2dofMax")),
0118       trkBendChi2Max_((float)iConfig.getParameter<double>("trk_bendChi2Max")),
0119       trkNPSStubMin_((int)iConfig.getParameter<int>("trk_nPSStubMin")),
0120       minTrkJetpT_(iConfig.getParameter<double>("minTrkJetpT")),
0121       etaBins_((int)iConfig.getParameter<int>("etaBins")),
0122       phiBins_((int)iConfig.getParameter<int>("phiBins")),
0123       zBins_((int)iConfig.getParameter<int>("zBins")),
0124       d0CutNStubs4_((float)iConfig.getParameter<double>("d0_cutNStubs4")),
0125       d0CutNStubs5_((float)iConfig.getParameter<double>("d0_cutNStubs5")),
0126       lowpTJetMinTrackMultiplicity_((int)iConfig.getParameter<int>("lowpTJetMinTrackMultiplicity")),
0127       lowpTJetMinpT_((float)iConfig.getParameter<double>("lowpTJetMinpT")),
0128       highpTJetMinTrackMultiplicity_((int)iConfig.getParameter<int>("highpTJetMinTrackMultiplicity")),
0129       highpTJetMinpT_((float)iConfig.getParameter<double>("highpTJetMinpT")),
0130       displaced_(iConfig.getParameter<bool>("displaced")),
0131       nStubs4DisplacedChi2_((float)iConfig.getParameter<double>("nStubs4DisplacedChi2")),
0132       nStubs5DisplacedChi2_((float)iConfig.getParameter<double>("nStubs5DisplacedChi2")),
0133       nStubs4DisplacedBend_((float)iConfig.getParameter<double>("nStubs4Displacedbend")),
0134       nStubs5DisplacedBend_((float)iConfig.getParameter<double>("nStubs5Displacedbend")),
0135       nDisplacedTracks_((int)iConfig.getParameter<int>("nDisplacedTracks")),
0136       trackToken_(consumes<vector<TTTrack<Ref_Phase2TrackerDigi_>>>(iConfig.getParameter<InputTag>("L1TrackInputTag"))),
0137       PVtxToken_(consumes<l1t::VertexWordCollection>(iConfig.getParameter<InputTag>("VertexInputTag"))),
0138       tTopoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>(edm::ESInputTag("", ""))) {
0139   zStep_ = convert::makeZ0(2.0 * trkZMax_ / (zBins_ + 1));
0140   etaStep_ = convert::makeGlbEta(2.0 * trkEtaMax_ / etaBins_);  //etaStep is the width of an etabin
0141   phiStep_ = convert::makeGlbPhi(2.0 * M_PI / phiBins_);        ////phiStep is the width of a phibin
0142 
0143   if (displaced_)
0144     produces<l1t::TkJetWordCollection>("L1TrackJetsExtended");
0145   else
0146     produces<l1t::TkJetWordCollection>("L1TrackJets");
0147 }
0148 
0149 L1TrackJetEmulationProducer::~L1TrackJetEmulationProducer() {}
0150 
0151 void L1TrackJetEmulationProducer::produce(Event &iEvent, const EventSetup &iSetup) {
0152   unique_ptr<l1t::TkJetWordCollection> L1L1TrackJetProducer(new l1t::TkJetWordCollection);
0153 
0154   // For TTStubs
0155   const TrackerTopology &tTopo = iSetup.getData(tTopoToken_);
0156 
0157   edm::Handle<vector<TTTrack<Ref_Phase2TrackerDigi_>>> TTTrackHandle;
0158   iEvent.getByToken(trackToken_, TTTrackHandle);
0159   vector<TTTrack<Ref_Phase2TrackerDigi_>>::const_iterator iterL1Track;
0160 
0161   edm::Handle<l1t::VertexWordCollection> PVtx;
0162   iEvent.getByToken(PVtxToken_, PVtx);
0163   float PVz = (PVtx->at(0)).z0();
0164 
0165   L1TrkPtrs_.clear();
0166   zBinCount_.clear();
0167   tdtrk_.clear();
0168   unsigned int this_l1track = 0;
0169   for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) {
0170     edm::Ptr<L1TTTrackType> trkPtr(TTTrackHandle, this_l1track);
0171     this_l1track++;
0172     float trk_pt = trkPtr->momentum().perp();
0173     int trk_nstubs = (int)trkPtr->getStubRefs().size();
0174     float trk_chi2dof = trkPtr->chi2Red();
0175     float trk_d0 = trkPtr->d0();
0176     float trk_bendchi2 = trkPtr->stubPtConsistency();
0177     float trk_z0 = trkPtr->z0();
0178 
0179     int trk_nPS = 0;
0180     for (int istub = 0; istub < trk_nstubs; istub++) {  // loop over the stubs
0181       DetId detId(trkPtr->getStubRefs().at(istub)->getDetId());
0182       if (detId.det() == DetId::Detector::Tracker) {
0183         if ((detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3) ||
0184             (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9))
0185           trk_nPS++;
0186       }
0187     }
0188 
0189     if (trk_nPS < trkNPSStubMin_)
0190       continue;
0191     if (!trackQualityCuts(trk_nstubs, trk_chi2dof, trk_bendchi2))
0192       continue;
0193     if (std::abs(trk_z0 - PVz) > MaxDzTrackPV)
0194       continue;
0195     if (std::abs(trk_z0) > trkZMax_)
0196       continue;
0197     if (std::abs(trkPtr->momentum().eta()) > trkEtaMax_)
0198       continue;
0199     if (trk_pt < trkPtMin_)
0200       continue;
0201     L1TrkPtrs_.push_back(trkPtr);
0202     zBinCount_.push_back(0);
0203 
0204     if ((std::abs(trk_d0) > d0CutNStubs5_ && trk_nstubs >= 5 && d0CutNStubs5_ >= 0) ||
0205         (trk_nstubs == 4 && std::abs(trk_d0) > d0CutNStubs4_ && d0CutNStubs4_ >= 0))
0206       tdtrk_.push_back(1);  //displaced track
0207     else
0208       tdtrk_.push_back(0);  // not displaced track
0209   }
0210 
0211   if (!L1TrkPtrs_.empty()) {
0212     TrackJetEmulationMaxZBin mzb;
0213 
0214     L2_cluster(L1TrkPtrs_, mzb);
0215     if (mzb.clusters != nullptr) {
0216       for (int j = 0; j < mzb.nclust; ++j) {
0217         //FILL Two Layer Jets for Jet Collection
0218         if (mzb.clusters[j].pTtot < convert::makePtFromFloat(trkPtMin_))
0219           continue;  //protects against reading bad memory
0220         if (mzb.clusters[j].ntracks < 1)
0221           continue;
0222         if (mzb.clusters[j].ntracks > 5000)
0223           continue;
0224         l1t::TkJetWord::glbeta_t jetEta = mzb.clusters[j].eta * convert::ETA_LSB_POW;
0225         l1t::TkJetWord::glbphi_t jetPhi = mzb.clusters[j].phi * convert::PHI_LSB_POW;
0226         l1t::TkJetWord::z0_t jetZ0 = mzb.zbincenter * convert::Z0_LSB_POW;
0227         l1t::TkJetWord::pt_t jetPt = mzb.clusters[j].pTtot;
0228         l1t::TkJetWord::nt_t totalntracks_ = mzb.clusters[j].ntracks;
0229         l1t::TkJetWord::nx_t totalxtracks_ = mzb.clusters[j].nxtracks;
0230         l1t::TkJetWord::tkjetunassigned_t unassigned_ = 0;
0231 
0232         l1t::TkJetWord trkJet(jetPt, jetEta, jetPhi, jetZ0, totalntracks_, totalxtracks_, unassigned_);
0233         //trkJet.setDispCounters(DispCounters);
0234         L1L1TrackJetProducer->push_back(trkJet);
0235       }
0236     } else if (mzb.clusters == nullptr) {
0237       edm::LogWarning("L1TrackJetEmulationProducer") << "mzb.clusters Not Assigned!\n";
0238     }
0239 
0240     if (displaced_)
0241       iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended");
0242     else
0243       iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets");
0244     delete[] mzb.clusters;
0245   } else if (L1TrkPtrs_.empty()) {
0246     edm::LogWarning("L1TrackJetEmulationProducer") << "L1TrkPtrs Not Assigned!\n";
0247     if (displaced_)
0248       iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended");
0249     else
0250       iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets");
0251   }
0252 }
0253 
0254 void L1TrackJetEmulationProducer::L2_cluster(vector<Ptr<L1TTTrackType>> L1TrkPtrs_, TrackJetEmulationMaxZBin &mzb) {
0255   enum TrackBitWidths {
0256     kEtaSize = 16,    // Width of z-position (40cm / 0.1)
0257     kEtaMagSize = 3,  // Width of z-position magnitude (signed)
0258     kPtSize = 14,     // Width of pt
0259     kPtMagSize = 9,   // Width of pt magnitude (unsigned)
0260     kPhiSize = 12,
0261     kPhiMagSize = 1,
0262   };
0263 
0264   const int nz = zBins_ + 1;
0265   TrackJetEmulationMaxZBin all_zBins[nz];
0266   static TrackJetEmulationMaxZBin mzbtemp;
0267   for (int z = 0; z < nz; ++z)
0268     all_zBins[z] = mzbtemp;
0269 
0270   z0_intern zmin = convert::makeZ0(-1.0 * trkZMax_);
0271   z0_intern zmax = zmin + 2 * zStep_;
0272 
0273   TrackJetEmulationEtaPhiBin epbins[phiBins_][etaBins_];  // create grid of phiBins
0274   glbphi_intern phi = convert::makeGlbPhi(-1.0 * M_PI);
0275   glbeta_intern eta;
0276   glbeta_intern etamin, etamax, phimin, phimax;
0277   for (int i = 0; i < phiBins_; ++i) {
0278     eta = convert::makeGlbEta(-1.0 * trkEtaMax_);
0279     for (int j = 0; j < etaBins_; ++j) {
0280       phimin = phi;
0281       phimax = phi + phiStep_;
0282       etamin = eta;
0283       eta = eta + etaStep_;
0284       etamax = eta;
0285       epbins[i][j].phi = (phimin + phimax) / 2;
0286       epbins[i][j].eta = (etamin + etamax) / 2;
0287       epbins[i][j].pTtot = 0;
0288       epbins[i][j].ntracks = 0;
0289       epbins[i][j].nxtracks = 0;
0290     }  // for each etabin
0291     phi = phi + phiStep_;
0292   }  // for each phibin (finished creating epbins)
0293 
0294   mzb = all_zBins[0];
0295   mzb.ht = 0;
0296   int ntracks = L1TrkPtrs_.size();
0297   // uninitalized arrays
0298   TrackJetEmulationEtaPhiBin *L1clusters[phiBins_];
0299   TrackJetEmulationEtaPhiBin L2cluster[ntracks];
0300 
0301   for (int zbin = 0; zbin < zBins_; ++zbin) {
0302     for (int i = 0; i < phiBins_; ++i) {  //First initialize pT, numtracks, used to 0 (or false)
0303       for (int j = 0; j < etaBins_; ++j) {
0304         epbins[i][j].pTtot = 0;
0305         epbins[i][j].used = false;
0306         epbins[i][j].ntracks = 0;
0307         epbins[i][j].nxtracks = 0;
0308       }  //for each etabin
0309       L1clusters[i] = epbins[i];
0310     }  //for each phibin
0311 
0312     for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) {
0313       ap_ufixed<TrackBitWidths::kPtSize, TrackBitWidths::kPtMagSize, AP_RND_CONV, AP_SAT> inputTrkPt = 0;
0314       inputTrkPt.V = L1TrkPtrs_[k]->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1,
0315                                                    TTTrack_TrackWord::TrackBitLocations::kRinvLSB);
0316       pt_intern trkpt = inputTrkPt;
0317 
0318       ap_fixed<TrackBitWidths::kEtaSize, TrackBitWidths::kEtaMagSize, AP_RND_CONV, AP_SAT> trketainput = 0;
0319       trketainput.V = L1TrkPtrs_[k]->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kTanlMSB,
0320                                                     TTTrack_TrackWord::TrackBitLocations::kTanlLSB);
0321       ap_ufixed<64 + ETA_EXTRABITS, 32 + ETA_EXTRABITS> eta_conv =
0322           1.0 / convert::ETA_LSB;  //conversion factor from input eta format to internal format
0323       glbeta_intern trketa = eta_conv * trketainput;
0324 
0325       int Sector = L1TrkPtrs_[k]->phiSector();
0326       glbphi_intern trkphiSector = 0;
0327       if (Sector < 5) {
0328         trkphiSector = Sector * convert::makeGlbPhi(2.0 * M_PI / 9.0);
0329       } else {
0330         trkphiSector =
0331             convert::makeGlbPhi(-1.0 * M_PI + M_PI / 9.0) + (Sector - 5) * convert::makeGlbPhi(2.0 * M_PI / 9.0);
0332       }
0333 
0334       ap_fixed<TrackBitWidths::kPhiSize, TrackBitWidths::kPhiMagSize, AP_RND_CONV, AP_SAT> trkphiinput = 0;
0335       trkphiinput.V = L1TrkPtrs_[k]->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kPhiMSB,
0336                                                     TTTrack_TrackWord::TrackBitLocations::kPhiLSB);
0337       ap_ufixed<64 + PHI_EXTRABITS, 32 + PHI_EXTRABITS> phi_conv =
0338           TTTrack_TrackWord::stepPhi0 / convert::PHI_LSB *
0339           (1 << (TrackBitWidths::kPhiSize - TrackBitWidths::kPhiMagSize));
0340       //phi_conv is a conversion factor from input phi format to the internal format
0341       glbeta_intern localphi = phi_conv * trkphiinput;
0342       glbeta_intern trkphi = localphi + trkphiSector;
0343 
0344       ap_int<TTTrack_TrackWord::TrackBitWidths::kZ0Size> inputTrkZ0 = L1TrkPtrs_[k]->getTrackWord()(
0345           TTTrack_TrackWord::TrackBitLocations::kZ0MSB, TTTrack_TrackWord::TrackBitLocations::kZ0LSB);
0346       ap_ufixed<32 + Z0_EXTRABITS, 1 + Z0_EXTRABITS> z0_conv =
0347           TTTrack_TrackWord::stepZ0 / convert::Z0_LSB;  //conversion factor from input z format to internal format
0348       z0_intern trkZ = z0_conv * inputTrkZ0;
0349 
0350       for (int i = 0; i < phiBins_; ++i) {
0351         for (int j = 0; j < etaBins_; ++j) {
0352           L2cluster[k] = epbins[i][j];
0353           if ((zmin <= trkZ && zmax >= trkZ) &&
0354               ((epbins[i][j].eta - etaStep_ / 2 < trketa && epbins[i][j].eta + etaStep_ / 2 >= trketa) &&
0355                epbins[i][j].phi - phiStep_ / 2 < trkphi && epbins[i][j].phi + phiStep_ / 2 >= trkphi &&
0356                (zBinCount_[k] != 2))) {
0357             zBinCount_.at(k) = zBinCount_.at(k) + 1;
0358 
0359             if (trkpt < convert::makePtFromFloat(trkPtMax_))
0360               epbins[i][j].pTtot += trkpt;
0361             else
0362               epbins[i][j].pTtot += convert::makePtFromFloat(trkPtMax_);
0363             ++epbins[i][j].ntracks;
0364             //x-bit is currently not used in firmware, so we leave nxtracks = 0 for now
0365           }  // if right bin
0366         }    // for each phibin: j loop
0367       }      // for each phibin: i loop
0368     }        // end loop over tracks
0369 
0370     for (int phislice = 0; phislice < phiBins_; ++phislice) {
0371       L1clusters[phislice] = L1_cluster(epbins[phislice]);
0372       for (int ind = 0; L1clusters[phislice][ind].pTtot != 0; ++ind) {
0373         L1clusters[phislice][ind].used = false;
0374       }
0375     }
0376 
0377     //Create clusters array to hold output cluster data for Layer2; can't have more clusters than tracks.
0378     //Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used.
0379     pt_intern hipT = 0;
0380     int nclust = 0;
0381     int phibin = 0;
0382     int imax = -1;
0383     int index1;  //index of clusters array for each phislice
0384     pt_intern E1 = 0;
0385     pt_intern E0 = 0;
0386     pt_intern E2 = 0;
0387     l1t::TkJetWord::nt_t ntrk1, ntrk2;
0388     l1t::TkJetWord::nx_t nxtrk1, nxtrk2;
0389     int used1, used2, used3, used4;
0390 
0391     for (phibin = 0; phibin < phiBins_; ++phibin) {  //Find eta-phibin with highest pT
0392       while (true) {
0393         hipT = 0;
0394         for (index1 = 0; L1clusters[phibin][index1].pTtot > 0; ++index1) {
0395           if (!L1clusters[phibin][index1].used && L1clusters[phibin][index1].pTtot >= hipT) {
0396             hipT = L1clusters[phibin][index1].pTtot;
0397             imax = index1;
0398           }
0399         }  // for each index within the phibin
0400 
0401         if (hipT == 0)
0402           break;    // If highest pT is 0, all bins are used
0403         E0 = hipT;  // E0 is pT of first phibin of the cluster
0404         E1 = 0;
0405         E2 = 0;
0406         ntrk1 = 0;
0407         ntrk2 = 0;
0408         nxtrk1 = 0;
0409         nxtrk2 = 0;
0410         L2cluster[nclust] = L1clusters[phibin][imax];
0411 
0412         L1clusters[phibin][imax].used = true;
0413         // Add pT of upper neighbor
0414         // E1 is pT of the middle phibin (should be highest pT)
0415         if (phibin != phiBins_ - 1) {
0416           used1 = -1;
0417           used2 = -1;
0418           for (index1 = 0; L1clusters[phibin + 1][index1].pTtot != 0; ++index1) {
0419             if (L1clusters[phibin + 1][index1].used)
0420               continue;
0421             if (L1clusters[phibin + 1][index1].eta - L1clusters[phibin][imax].eta <= 3 * etaStep_ / 2 &&
0422                 L1clusters[phibin][imax].eta - L1clusters[phibin + 1][index1].eta <= 3 * etaStep_ / 2) {
0423               E1 += L1clusters[phibin + 1][index1].pTtot;
0424               ntrk1 += L1clusters[phibin + 1][index1].ntracks;
0425               nxtrk1 += L1clusters[phibin + 1][index1].nxtracks;
0426               if (used1 < 0)
0427                 used1 = index1;
0428               else
0429                 used2 = index1;
0430             }  // if cluster is within one phibin
0431           }    // for each cluster in above phibin
0432 
0433           if (E1 < E0) {  // if E1 isn't higher, E0 and E1 are their own cluster
0434             L2cluster[nclust].pTtot += E1;
0435             L2cluster[nclust].ntracks += ntrk1;
0436             L2cluster[nclust].nxtracks += nxtrk1;
0437             if (used1 >= 0)
0438               L1clusters[phibin + 1][used1].used = true;
0439             if (used2 >= 0)
0440               L1clusters[phibin + 1][used2].used = true;
0441             nclust++;
0442             continue;
0443           }
0444 
0445           if (phibin != phiBins_ - 2) {  // E2 will be the pT of the third phibin (should be lower than E1)
0446             used3 = -1;
0447             used4 = -1;
0448             for (index1 = 0; L1clusters[phibin + 2][index1].pTtot != 0; ++index1) {
0449               if (L1clusters[phibin + 2][index1].used)
0450                 continue;
0451               if (L1clusters[phibin + 2][index1].eta - L1clusters[phibin][imax].eta <= 3 * etaStep_ / 2 &&
0452                   L1clusters[phibin][imax].eta - L1clusters[phibin + 2][index1].eta <= 3 * etaStep_ / 2) {
0453                 E2 += L1clusters[phibin + 2][index1].pTtot;
0454                 ntrk2 += L1clusters[phibin + 2][index1].ntracks;
0455                 nxtrk2 += L1clusters[phibin + 2][index1].nxtracks;
0456                 if (used3 < 0)
0457                   used3 = index1;
0458                 else
0459                   used4 = index1;
0460               }
0461             }
0462             // if indeed E2 < E1, add E1 and E2 to E0, they're all a cluster together
0463             // otherwise, E0 is its own cluster
0464             if (E2 < E1) {
0465               L2cluster[nclust].pTtot += E1 + E2;
0466               L2cluster[nclust].ntracks += ntrk1 + ntrk2;
0467               L2cluster[nclust].nxtracks += nxtrk1 + nxtrk2;
0468               L2cluster[nclust].phi = L1clusters[phibin + 1][used1].phi;
0469               if (used1 >= 0)
0470                 L1clusters[phibin + 1][used1].used = true;
0471               if (used2 >= 0)
0472                 L1clusters[phibin + 1][used2].used = true;
0473               if (used3 >= 0)
0474                 L1clusters[phibin + 2][used3].used = true;
0475               if (used4 >= 0)
0476                 L1clusters[phibin + 2][used4].used = true;
0477             }
0478             nclust++;
0479             continue;
0480           }  // end Not phiBins-2
0481           else {
0482             L2cluster[nclust].pTtot += E1;
0483             L2cluster[nclust].ntracks += ntrk1;
0484             L2cluster[nclust].nxtracks += nxtrk1;
0485             L2cluster[nclust].phi = L1clusters[phibin + 1][used1].phi;
0486             if (used1 >= 0)
0487               L1clusters[phibin + 1][used1].used = true;
0488             if (used2 >= 0)
0489               L1clusters[phibin + 1][used2].used = true;
0490             nclust++;
0491             continue;
0492           }
0493         }       //End not last phibin(23)
0494         else {  //if it is phibin 23
0495           L1clusters[phibin][imax].used = true;
0496           nclust++;
0497         }
0498       }  // while hipT not 0
0499     }    // for each phibin
0500 
0501     for (phibin = 0; phibin < phiBins_; ++phibin)
0502       delete[] L1clusters[phibin];
0503 
0504     // Now merge clusters, if necessary
0505     for (int m = 0; m < nclust - 1; ++m) {
0506       for (int n = m + 1; n < nclust; ++n) {
0507         if (L2cluster[n].eta == L2cluster[m].eta &&
0508             ((L2cluster[n].phi - L2cluster[m].phi < 3 * phiStep_ / 2 &&
0509               L2cluster[m].phi - L2cluster[n].phi < 3 * phiStep_ / 2) ||
0510              (L2cluster[n].phi - L2cluster[m].phi > 2 * convert::makeGlbPhi(M_PI) - phiStep_ ||
0511               L2cluster[m].phi - L2cluster[n].phi > 2 * convert::makeGlbPhi(M_PI) - phiStep_))) {
0512           if (L2cluster[n].pTtot > L2cluster[m].pTtot) {
0513             L2cluster[m].phi = L2cluster[n].phi;
0514           }
0515           L2cluster[m].pTtot += L2cluster[n].pTtot;
0516           L2cluster[m].ntracks += L2cluster[n].ntracks;
0517           L2cluster[m].nxtracks += L2cluster[n].nxtracks;
0518           for (int m1 = n; m1 < nclust - 1; ++m1) {
0519             L2cluster[m1] = L2cluster[m1 + 1];
0520           }
0521           nclust--;
0522           m = -1;
0523           break;  //?????
0524         }         // end if clusters neighbor in eta
0525       }
0526     }  // end for (m) loop
0527     // sum up all pTs in this zbin to find ht
0528 
0529     pt_intern ht = 0;
0530     for (int k = 0; k < nclust; ++k) {
0531       if (L2cluster[k].pTtot > convert::makePtFromFloat(lowpTJetMinpT_) &&
0532           L2cluster[k].ntracks < lowpTJetMinTrackMultiplicity_)
0533         continue;
0534       if (L2cluster[k].pTtot > convert::makePtFromFloat(highpTJetMinpT_) &&
0535           L2cluster[k].ntracks < highpTJetMinTrackMultiplicity_)
0536         continue;
0537       if (L2cluster[k].pTtot > convert::makePtFromFloat(minTrkJetpT_))
0538         ht += L2cluster[k].pTtot;
0539     }
0540     // if ht is larger than previous max, this is the new vertex zbin
0541     all_zBins[zbin].znum = zbin;
0542     all_zBins[zbin].clusters = new TrackJetEmulationEtaPhiBin[nclust];
0543     all_zBins[zbin].nclust = nclust;
0544     all_zBins[zbin].zbincenter = (zmin + zmax) / 2;
0545     for (int k = 0; k < nclust; ++k) {
0546       all_zBins[zbin].clusters[k].phi = L2cluster[k].phi;
0547       all_zBins[zbin].clusters[k].eta = L2cluster[k].eta;
0548       all_zBins[zbin].clusters[k].pTtot = L2cluster[k].pTtot;
0549       all_zBins[zbin].clusters[k].ntracks = L2cluster[k].ntracks;
0550       all_zBins[zbin].clusters[k].nxtracks = L2cluster[k].nxtracks;
0551     }
0552     all_zBins[zbin].ht = ht;
0553     if (ht >= mzb.ht) {
0554       mzb = all_zBins[zbin];
0555       mzb.zbincenter = (zmin + zmax) / 2;
0556     }
0557     // Prepare for next zbin!
0558     zmin = zmin + zStep_;
0559     zmax = zmax + zStep_;
0560   }  // for each zbin
0561 
0562   for (int zbin = 0; zbin < zBins_; ++zbin) {
0563     if (zbin == mzb.znum) {
0564       continue;
0565     }
0566     delete[] all_zBins[zbin].clusters;
0567   }
0568 }
0569 
0570 TrackJetEmulationEtaPhiBin *L1TrackJetEmulationProducer::L1_cluster(TrackJetEmulationEtaPhiBin *phislice) {
0571   TrackJetEmulationEtaPhiBin *clusters = new TrackJetEmulationEtaPhiBin[etaBins_ / 2];
0572   for (int etabin = 0; etabin < etaBins_ / 2; ++etabin) {
0573     clusters[etabin].pTtot = 0;
0574     clusters[etabin].ntracks = 0;
0575     clusters[etabin].nxtracks = 0;
0576     clusters[etabin].phi = 0;
0577     clusters[etabin].eta = 0;
0578     clusters[etabin].used = false;
0579   }
0580 
0581   if (clusters == nullptr)
0582     edm::LogWarning("L1TrackJetEmulationProducer") << "Clusters memory not assigned!\n";
0583 
0584   // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used
0585   pt_intern my_pt, left_pt, right_pt, right2pt;
0586   int nclust = 0;
0587   right2pt = 0;
0588   for (int etabin = 0; etabin < etaBins_; ++etabin) {
0589     // assign values for my pT and neighbors' pT
0590     if (phislice[etabin].used)
0591       continue;
0592     my_pt = phislice[etabin].pTtot;
0593     if (etabin > 0 && !phislice[etabin - 1].used) {
0594       left_pt = phislice[etabin - 1].pTtot;
0595     } else
0596       left_pt = 0;
0597     if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) {
0598       right_pt = phislice[etabin + 1].pTtot;
0599       if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) {
0600         right2pt = phislice[etabin + 2].pTtot;
0601       } else
0602         right2pt = 0;
0603     } else
0604       right_pt = 0;
0605 
0606     // if I'm not a cluster, move on
0607     if (my_pt < left_pt || my_pt <= right_pt) {
0608       // if unused pT in the left neighbor, spit it out as a cluster
0609       if (left_pt > 0) {
0610         clusters[nclust] = phislice[etabin - 1];
0611         phislice[etabin - 1].used = true;
0612         nclust++;
0613       }
0614       continue;
0615     }
0616 
0617     // I guess I'm a cluster-- should I use my right neighbor?
0618     // Note: left neighbor will definitely be used because if it
0619     //       didn't belong to me it would have been used already
0620     clusters[nclust] = phislice[etabin];
0621     phislice[etabin].used = true;
0622     if (left_pt > 0) {
0623       clusters[nclust].pTtot += left_pt;
0624       clusters[nclust].ntracks += phislice[etabin - 1].ntracks;
0625       clusters[nclust].nxtracks += phislice[etabin - 1].nxtracks;
0626     }
0627     if (my_pt >= right2pt && right_pt > 0) {
0628       clusters[nclust].pTtot += right_pt;
0629       clusters[nclust].ntracks += phislice[etabin + 1].ntracks;
0630       clusters[nclust].nxtracks += phislice[etabin + 1].nxtracks;
0631       phislice[etabin + 1].used = true;
0632     }
0633     nclust++;
0634   }  // for each etabin
0635 
0636   // Now merge clusters, if necessary
0637   for (int m = 0; m < nclust - 1; ++m) {
0638     if (clusters[m + 1].eta - clusters[m].eta < 3 * etaStep_ / 2 &&
0639         clusters[m].eta - clusters[m + 1].eta < 3 * etaStep_ / 2) {
0640       if (clusters[m + 1].pTtot > clusters[m].pTtot) {
0641         clusters[m].eta = clusters[m + 1].eta;
0642       }
0643       clusters[m].pTtot += clusters[m + 1].pTtot;
0644       clusters[m].ntracks += clusters[m + 1].ntracks;  // Previous version didn't add tracks when merging
0645       clusters[m].nxtracks += clusters[m + 1].nxtracks;
0646       for (int m1 = m + 1; m1 < nclust - 1; ++m1)
0647         clusters[m1] = clusters[m1 + 1];
0648       nclust--;
0649       m = -1;
0650     }  // end if clusters neighbor in eta
0651   }    // end for (m) loop
0652 
0653   for (int i = nclust; i < etaBins_ / 2; ++i)  // zero out remaining unused clusters
0654     clusters[i].pTtot = 0;
0655   return clusters;
0656 }
0657 
0658 void L1TrackJetEmulationProducer::beginStream(StreamID) {}
0659 
0660 void L1TrackJetEmulationProducer::endStream() {}
0661 
0662 bool L1TrackJetEmulationProducer::trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2) {
0663   bool PassQuality = false;
0664   if (trk_bendchi2 < trkBendChi2Max_ && trk_chi2 < trkChi2dofMax_ && trk_nstub >= 4 && !displaced_)
0665     PassQuality = true;
0666   if (displaced_ && trk_bendchi2 < nStubs4DisplacedBend_ && trk_chi2 < nStubs4DisplacedChi2_ && trk_nstub == 4)
0667     PassQuality = true;
0668   if (displaced_ && trk_bendchi2 < nStubs5DisplacedBend_ && trk_chi2 < nStubs5DisplacedChi2_ && trk_nstub > 4)
0669     PassQuality = true;
0670   return PassQuality;
0671 }
0672 
0673 void L1TrackJetEmulationProducer::fillDescriptions(ConfigurationDescriptions &descriptions) {
0674   //The following says we do not know what parameters are allowed so do no validation
0675   // Please change this to state exactly what you do use, even if it is no parameters
0676   ParameterSetDescription desc;
0677   desc.setUnknown();
0678   descriptions.addDefault(desc);
0679 }
0680 
0681 //define this as a plug-in
0682 DEFINE_FWK_MODULE(L1TrackJetEmulationProducer);