Skip to content

Update HLTPixelThrustFilter to consider saturated pixels #50167

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cmsbuild merged 3 commits into from
Mar 7, 2026
Merged

Update HLTPixelThrustFilter to consider saturated pixels #50167

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
bf1d445
Updated HLTPixelThrustFilter to select saturated pixels
Mar 4, 2026
915148b
Count pixels before reserving
Mar 5, 2026
56faab0
Add upper limit on number of pixels
Mar 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 18 additions & 1 deletion HLTrigger/special/plugins/HLTPixelThrustFilter.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@ class HLTPixelThrustFilter : public edm::global::EDFilter<> {
private:
const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryRcdToken_;
const edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster> > inputToken_;
const bool useOnlySaturatedPixels_;
const unsigned int min_nPixels_, min_nSatPixels_;
const double min_thrust_; // minimum thrust
const double max_thrust_; // maximum thrust
};
Expand All @@ -36,12 +38,18 @@ class HLTPixelThrustFilter : public edm::global::EDFilter<> {
HLTPixelThrustFilter::HLTPixelThrustFilter(const edm::ParameterSet& config)
: trackerGeometryRcdToken_(esConsumes()),
inputToken_(consumes<edmNew::DetSetVector<SiPixelCluster> >(config.getParameter<edm::InputTag>("inputTag"))),
useOnlySaturatedPixels_(config.getParameter<bool>("useOnlySaturatedPixels")),
min_nPixels_(config.getParameter<unsigned int>("minNPixels")),
min_nSatPixels_(config.getParameter<unsigned int>("minNSaturatedPixels")),
min_thrust_(config.getParameter<double>("minThrust")),
max_thrust_(config.getParameter<double>("maxThrust")) {}

void HLTPixelThrustFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("inputTag", edm::InputTag("hltSiPixelClusters"));
desc.add<bool>("useOnlySaturatedPixels", false);
desc.add<unsigned int>("minNPixels", 2);
desc.add<unsigned int>("minNSaturatedPixels", 0);
desc.add<double>("minThrust", 0);
desc.add<double>("maxThrust", 0);
descriptions.add("hltPixelThrustFilter", desc);
Expand All @@ -56,15 +64,24 @@ bool HLTPixelThrustFilter::filter(edm::StreamID, edm::Event& event, edm::EventSe
auto const& clusters = event.get(inputToken_);
auto const& trackerGeo = iSetup.getData(trackerGeometryRcdToken_);

size_t nSatPixels(0);
std::vector<reco::LeafCandidate> vec;
vec.reserve(clusters.size() * 1000);
Copy link
Copy Markdown
Contributor

@mmusich mmusich Mar 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

where's the 1000 coming from?

Copy link
Copy Markdown
Contributor Author

@stahlleiton stahlleiton Mar 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I checked locally with HIForward data and UPC MC, and saw that there were up to around 250 clusters in a detset, so to be conservative I multiplied the number of detsets by 1000 to reserve the vector. If prefered, I could either remove the reserve or count the total number clusters before reserving.

Copy link
Copy Markdown
Contributor

@mmusich mmusich Mar 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure how stable that number is going to be at the very end of the lifetime of the detector.

Copy link
Copy Markdown
Contributor Author

@stahlleiton stahlleiton Mar 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok, in 915148b now it counts the total before reserving

for (auto DSViter = clusters.begin(); DSViter != clusters.end(); DSViter++) {
auto const& pgdu = static_cast<const PixelGeomDetUnit*>(trackerGeo.idToDetUnit(DSViter->detId()));
for (auto const& cluster : *DSViter) {
if (cluster.isSaturated())
nSatPixels += 1;
else if (useOnlySaturatedPixels_)
continue;
auto const& pos = pgdu->surface().toGlobal(pgdu->specificTopology().localPosition({cluster.x(), cluster.y()}));
auto const mag = std::sqrt(pos.x() * pos.x() + pos.y() * pos.y());
vec.emplace_back(0, reco::Particle::LorentzVector(pos.x() / mag, pos.y() / mag, 0, 0));
if (mag > 0)
vec.emplace_back(0, reco::Particle::LorentzVector(pos.x() / mag, pos.y() / mag, 0, 0));
}
}
if (vec.size() < min_nPixels_ || nSatPixels < min_nSatPixels_)
return false;
auto const thrust = Thrust(vec.begin(), vec.end()).thrust();

bool accept = (thrust >= min_thrust_);
Expand Down
1 change: 1 addition & 0 deletions PhysicsTools/CandUtils/BuildFile.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
<use name="DataFormats/Candidate"/>
<use name="roothistmatrix"/>
<use name="vdt_headers"/>
<export>
<lib name="1"/>
</export>
5 changes: 4 additions & 1 deletion PhysicsTools/CandUtils/interface/Thrust.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,20 +33,23 @@
*/
#include "DataFormats/Math/interface/Vector3D.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include <vdt/vdtMath.h>
#include <vector>

class Thrust {
public:
/// spatial vector
typedef math::XYZVector Vector;
static constexpr double PI = M_PI, PI2 = 2. * PI, PI_2 = PI / 2., PI_4 = PI / 4.;
/// constructor from first and last iterators
template <typename const_iterator>
Thrust(const_iterator begin, const_iterator end) : thrust_(0), axis_(0, 0, 0), pSum_(0), n_(end - begin), p_(n_) {
if (n_ == 0)
return;
std::vector<const reco::Candidate *> cands;
cands.reserve(n_);
for (const_iterator i = begin; i != end; ++i) {
cands.push_back(&*i);
cands.emplace_back(&*i);
}
init(cands);
}
Expand Down
86 changes: 42 additions & 44 deletions PhysicsTools/CandUtils/src/Thrust.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,8 @@
#include "PhysicsTools/CandUtils/interface/Thrust.h"
#include <cmath>
using namespace reco;
const double pi = M_PI, pi2 = 2 * pi, pi_2 = pi / 2, pi_4 = pi / 4;

void Thrust::init(const std::vector<const Candidate*>& cands) {
void Thrust::init(const std::vector<const reco::Candidate*>& cands) {
int i = 0;
for (std::vector<const Candidate*>::const_iterator t = cands.begin(); t != cands.end(); ++t, ++i)
for (std::vector<const reco::Candidate*>::const_iterator t = cands.begin(); t != cands.end(); ++t, ++i)
pSum_ += (p_[i] = (*t)->momentum()).r();
axis_ = axis(finalAxis(initialAxis()));
if (axis_.z() < 0)
Expand All @@ -14,21 +11,23 @@ void Thrust::init(const std::vector<const Candidate*>& cands) {
}

Thrust::ThetaPhi Thrust::initialAxis() const {
static const int nSegsTheta = 10, nSegsPhi = 10, nSegs = nSegsTheta * nSegsPhi;
int i, j;
constexpr int nSegsTheta = 10, nSegsPhi = 10, nSegs = nSegsTheta * nSegsPhi;
constexpr double fPhi = PI2 / nSegsPhi;
constexpr double fTheta = PI / (nSegsTheta - 1);
double thr[nSegs], max = 0;
int indI = 0, indJ = 0, index = -1;
for (i = 0; i < nSegsTheta; ++i) {
double z = cos(pi * i / (nSegsTheta - 1));
double r = sqrt(1 - z * z);
for (j = 0; j < nSegsPhi; ++j) {
double phi = pi2 * j / nSegsPhi;
thr[i * nSegsPhi + j] = thrust(Vector(r * cos(phi), r * sin(phi), z));
if (thr[i * nSegsPhi + j] > max) {
index = i * nSegsPhi + j;
int indI = 0, indJ = 0;
for (int i = 0; i < nSegsTheta; ++i) {
double z = vdt::fast_cos(i * fTheta);
double r = std::sqrt(1 - z * z);
for (int j = 0; j < nSegsPhi; ++j) {
double sinPhi, cosPhi;
vdt::fast_sincos(j * fPhi, sinPhi, cosPhi);
double val = thrust(Vector(r * cosPhi, r * sinPhi, z));
thr[i * nSegsPhi + j] = val;
if (val > max) {
indI = i;
indJ = j;
max = thr[index];
max = val;
}
}
}
Expand Down Expand Up @@ -61,11 +60,11 @@ Thrust::ThetaPhi Thrust::initialAxis() const {
if (a != 0)
maxThetaInd = -b / (2 * a);
}
return ThetaPhi(pi * (maxThetaInd + indI) / (nSegsTheta - 1), pi2 * (maxPhiInd + indJ) / nSegsPhi);
return ThetaPhi((maxThetaInd + indI) * fTheta, (maxPhiInd + indJ) * fPhi);
}

Thrust::ThetaPhi Thrust::finalAxis(ThetaPhi best) const {
static const double epsilon = 0.0001;
constexpr double epsilon = 0.0001;
double maxChange1 = 0.0, maxChange2 = 0.0, a = 0.0, b = 0.0, c = 0.0, thr = 0.0;
int mandCt = 3, maxCt = 1000;
bool done;
Expand All @@ -74,12 +73,12 @@ Thrust::ThetaPhi Thrust::finalAxis(ThetaPhi best) const {
maxChange1 = 10 * (b < 0 ? -1 : 1);
if (a != 0)
maxChange1 = -b / (2 * a);
while (fabs(maxChange1 * epsilon) > pi_4)
while (std::abs(maxChange1 * epsilon) > PI_4)
maxChange1 /= 2;
if (maxChange1 == 0 && (best.theta == 0 || best.theta == pi)) {
best.phi += pi_2;
if (best.phi > pi2)
best.phi -= pi2;
if (maxChange1 == 0 && (best.theta == 0 || best.theta == PI)) {
best.phi += PI_2;
if (best.phi > PI2)
best.phi -= PI2;
parabola(a, b, c, axis(best), axis(best.theta + epsilon, best.phi), axis(best.theta - epsilon, best.phi));
maxChange1 = 10 * (b < 0 ? -1 : 1);
if (a != 0)
Expand All @@ -93,23 +92,23 @@ Thrust::ThetaPhi Thrust::finalAxis(ThetaPhi best) const {
} while (thr < c);

best.theta += maxChange1 * epsilon;
if (best.theta > pi) {
best.theta = pi - (best.theta - pi);
best.phi += pi;
if (best.phi > pi2)
best.phi -= pi2;
if (best.theta > PI) {
best.theta = PI - (best.theta - PI);
best.phi += PI;
if (best.phi > PI2)
best.phi -= PI2;
}
if (best.theta < 0) {
best.theta *= -1;
best.phi += pi;
if (best.phi > pi2)
best.phi -= pi2;
best.phi += PI;
if (best.phi > PI2)
best.phi -= PI2;
}
parabola(a, b, c, axis(best), axis(best.theta, best.phi + epsilon), axis(best.theta, best.phi - epsilon));
maxChange2 = 10 * (b < 0 ? -1 : 1);
if (a != 0)
maxChange2 = -b / (2 * a);
while (fabs(maxChange2 * epsilon) > pi_4) {
while (std::abs(maxChange2 * epsilon) > PI_4) {
maxChange2 /= 2;
}
do {
Expand All @@ -119,14 +118,14 @@ Thrust::ThetaPhi Thrust::finalAxis(ThetaPhi best) const {
maxChange2 /= 2;
} while (thr < c);
best.phi += maxChange2 * epsilon;
if (best.phi > pi2)
best.phi -= pi2;
if (best.phi > PI2)
best.phi -= PI2;
if (best.phi < 0)
best.phi += pi2;
best.phi += PI2;
if (mandCt > 0)
mandCt--;
maxCt--;
done = (fabs(maxChange1) > 1 || fabs(maxChange2) > 1 || mandCt) && (maxCt > 0);
done = (std::abs(maxChange1) > 1 || std::abs(maxChange2) > 1 || mandCt) && (maxCt > 0);
} while (done);

return best;
Expand All @@ -140,16 +139,15 @@ void Thrust::parabola(double& a, double& b, double& c, const Vector& a1, const V
}

Thrust::Vector Thrust::axis(double theta, double phi) const {
double theSin = sin(theta);
return Vector(theSin * cos(phi), theSin * sin(phi), cos(theta));
double sinTheta, cosTheta, sinPhi, cosPhi;
vdt::fast_sincos(theta, sinTheta, cosTheta);
vdt::fast_sincos(phi, sinPhi, cosPhi);
return Vector(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta);
}

double Thrust::thrust(const Vector& axis) const {
double result = 0;
double sum = 0;
for (unsigned int i = 0; i < n_; ++i)
sum += fabs(axis.Dot(p_[i]));
if (pSum_ > 0)
result = sum / pSum_;
return result;
sum += std::abs(axis.Dot(p_[i]));
return pSum_ > 0 ? sum / pSum_ : 0;
}