Skip to content

Commit

Permalink
address Adriano's comments
Browse files Browse the repository at this point in the history
  • Loading branch information
ccaillol committed Sep 15, 2022
1 parent 4fc043d commit 626df8f
Show file tree
Hide file tree
Showing 3 changed files with 236 additions and 282 deletions.
246 changes: 231 additions & 15 deletions L1Trigger/L1TTrackMatch/interface/L1TrackJetProducer.h
View file
Original file line number Diff line number Diff line change
@@ -1,21 +1,16 @@
#ifndef L1Trigger_L1TTrackMatch_L1TrackJetProducer_HH
#define L1Trigger_L1TTrackMatch_L1TrackJetProducer_HH
#pragma once
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>

//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<unsigned int> trackidx;
Expand All @@ -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<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
};
#endif

inline std::vector<EtaPhiBin> L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) {
std::vector<EtaPhiBin> 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<unsigned int> 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<EtaPhiBin> L2_clustering(std::vector<std::vector<EtaPhiBin>> &L1clusters,
int phiBins_,
float phiStep_,
float etaStep_) {
std::vector<EtaPhiBin> 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<unsigned int> trkidx1;
std::vector<unsigned int> 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<unsigned int> 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<pt1, merge both clusters
Fill_L2Cluster(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1);
for (unsigned int iused : used_already)
L1clusters[phibin + 1][iused].used = true;
continue;
}
// if phi = next to last bin there is no "next to next"
if (phibin == phiBins_ - 2) {
Fill_L2Cluster(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1);
clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi;
for (unsigned int iused : used_already)
L1clusters[phibin + 1][iused].used = true;
continue;
}
std::vector<int> 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<unsigned int> 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;
}
20 changes: 5 additions & 15 deletions L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <TMath.h>
Expand Down Expand Up @@ -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
Expand All @@ -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;
Expand Down Expand Up @@ -227,18 +229,6 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) {
L1clusters.reserve(phiBins_);
std::vector<EtaPhiBin> 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];
Expand Down
Loading

0 comments on commit 626df8f

Please sign in to comment.