File indexing completed on 2024-07-28 22:48:48
0001 #include <memory>
0002
0003 #include "RecoParticleFlow/PFTracking/interface/PFDisplacedVertexCandidateFinder.h"
0004
0005 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0006 #include "DataFormats/Math/interface/deltaPhi.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008 #include "DataFormats/VertexReco/interface/Vertex.h"
0009
0010 #include "MagneticField/Engine/interface/MagneticField.h"
0011
0012 using namespace std;
0013 using namespace reco;
0014
0015
0016
0017
0018 PFDisplacedVertexCandidateFinder::PFDisplacedVertexCandidateFinder()
0019 : vertexCandidates_(new PFDisplacedVertexCandidateCollection),
0020 dcaCut_(1000),
0021 primaryVertexCut2_(0.0),
0022 dcaPInnerHitCut2_(1000.0),
0023 vertexCandidatesSize_(50),
0024 theMinimum_(TwoTrackMinimumDistance::SlowMode),
0025 debug_(false),
0026 magField_(nullptr) {}
0027
0028 PFDisplacedVertexCandidateFinder::~PFDisplacedVertexCandidateFinder() {
0029 LogDebug("PFDisplacedVertexCandidateFinder")
0030 << "~PFDisplacedVertexCandidateFinder - number of remaining elements: " << eventTracks_.size();
0031 }
0032
0033 void PFDisplacedVertexCandidateFinder::setPrimaryVertex(edm::Handle<reco::VertexCollection> mainVertexHandle,
0034 edm::Handle<reco::BeamSpot> beamSpotHandle) {
0035 const math::XYZPoint beamSpot = beamSpotHandle.isValid()
0036 ? math::XYZPoint(beamSpotHandle->x0(), beamSpotHandle->y0(), beamSpotHandle->z0())
0037 : math::XYZPoint(0, 0, 0);
0038
0039
0040
0041
0042 pvtx_ = mainVertexHandle.isValid()
0043 ? math::XYZPoint(
0044 mainVertexHandle->begin()->x(), mainVertexHandle->begin()->y(), mainVertexHandle->begin()->z())
0045 : beamSpot;
0046 }
0047
0048
0049
0050 void PFDisplacedVertexCandidateFinder::setInput(const edm::Handle<TrackCollection>& trackh,
0051 const MagneticField* magField) {
0052 magField_ = magField;
0053 trackMask_.clear();
0054 trackMask_.resize(trackh->size());
0055 eventTrackTrajectories_.clear();
0056 eventTrackTrajectories_.resize(trackh->size());
0057
0058 for (unsigned i = 0; i < trackMask_.size(); i++)
0059 trackMask_[i] = true;
0060
0061 eventTracks_.clear();
0062 if (trackh.isValid()) {
0063 for (unsigned i = 0; i < trackh->size(); i++) {
0064 TrackRef tref(trackh, i);
0065 TrackBaseRef tbref(tref);
0066
0067 if (!isSelected(tbref)) {
0068 trackMask_[i] = false;
0069 continue;
0070 }
0071
0072 const Track* trk = tref.get();
0073 eventTracks_.push_back(tbref);
0074 eventTrackTrajectories_[i] = getGlobalTrajectoryParameters(trk);
0075 }
0076 }
0077 track_table_ = edm::soa::PtEtaPhiTable{*trackh};
0078 }
0079
0080
0081
0082 void PFDisplacedVertexCandidateFinder::findDisplacedVertexCandidates() {
0083 LogDebug("PFDisplacedVertexCandidateFinder") << "========= Start Finding Displaced Vertex Candidates =========";
0084
0085
0086 if (vertexCandidates_.get())
0087 vertexCandidates_->clear();
0088 else
0089 vertexCandidates_ = std::make_unique<PFDisplacedVertexCandidateCollection>();
0090
0091 vertexCandidates_->reserve(vertexCandidatesSize_);
0092 for (IE ie = eventTracks_.begin(); ie != eventTracks_.end();) {
0093
0094
0095
0096 PFDisplacedVertexCandidate tempVertexCandidate;
0097
0098 ie = associate(eventTracks_.end(), ie, tempVertexCandidate);
0099 LogDebug("PFDisplacedVertexCandidateFinder") << "elements=" << tempVertexCandidate.elements().size();
0100
0101
0102 if (tempVertexCandidate.isValid()) {
0103 packLinks(tempVertexCandidate);
0104 vertexCandidates_->push_back(tempVertexCandidate);
0105 }
0106 }
0107
0108 LogDebug("PFDisplacedVertexCandidateFinder") << "========= End Finding Displaced Vertex Candidates =========";
0109 }
0110
0111 PFDisplacedVertexCandidateFinder::IE PFDisplacedVertexCandidateFinder::associate(
0112 IE last, IE next, PFDisplacedVertexCandidate& tempVertexCandidate) {
0113 LogDebug("PFDisplacedVertexCandidateFinder") << "== Start the association procedure ==";
0114
0115 if (last != eventTracks_.end()) {
0116 double dist = -1;
0117 GlobalPoint crossingPoint(0, 0, 0);
0118 PFDisplacedVertexCandidate::VertexLinkTest linktest;
0119 link(*last, *next, dist, crossingPoint, linktest);
0120
0121 if (dist < -0.5) {
0122 LogDebug("PFDisplacedVertexCandidateFinder") << "link failed";
0123
0124 return ++next;
0125 } else {
0126
0127 LogDebug("PFDisplacedVertexCandidateFinder")
0128 << "adding element " << (*next).key() << " to PFDisplacedVertexCandidate with "
0129 << tempVertexCandidate.elements().size() << " elements";
0130 tempVertexCandidate.addElement((*next));
0131 trackMask_[(*next).key()] = false;
0132
0133 LogDebug("PFDisplacedVertexCandidateFinder")
0134 << "link parameters "
0135 << " *next = " << (*next).key() << " *last = " << (*last).key() << " dist = " << crossingPoint
0136 << " P.x = " << crossingPoint.x() << " P.y = " << crossingPoint.y() << " P.z = " << crossingPoint.z();
0137 }
0138 } else {
0139
0140 LogDebug("PFDisplacedVertexCandidateFinder") << "adding to block element " << (*next).key();
0141 tempVertexCandidate.addElement((*next));
0142 trackMask_[(*next).key()] = false;
0143 }
0144
0145
0146 #ifdef EDM_ML_DEBUG
0147 for (unsigned i = 0; i < trackMask_.size(); i++)
0148 LogTrace("PFDisplacedVertexCandidateFinder") << " Mask[" << i << "] = " << trackMask_[i];
0149 #endif
0150
0151 for (IE ie = eventTracks_.begin(); ie != eventTracks_.end();) {
0152 if (ie == last || ie == next) {
0153 ++ie;
0154 continue;
0155 }
0156
0157
0158 if (!trackMask_[(*ie).key()]) {
0159 ++ie;
0160 continue;
0161 }
0162
0163 LogDebug("PFDisplacedVertexCandidateFinder") << "calling associate " << (*next).key() << " & " << (*ie).key();
0164 ie = associate(next, ie, tempVertexCandidate);
0165 }
0166
0167 LogDebug("PFDisplacedVertexCandidateFinder") << "**** removing element ";
0168
0169 IE iteratorToNextFreeElement = eventTracks_.erase(next);
0170
0171 LogDebug("PFDisplacedVertexCandidateFinder") << "== End the association procedure ==";
0172
0173 return iteratorToNextFreeElement;
0174 }
0175
0176 void PFDisplacedVertexCandidateFinder::link(const TrackBaseRef& el1,
0177 const TrackBaseRef& el2,
0178 double& dist,
0179 GlobalPoint& crossingPoint,
0180 PFDisplacedVertexCandidate::VertexLinkTest& vertexLinkTest) {
0181 using namespace edm::soa::col;
0182 const auto iel1 = el1.key();
0183 const auto iel2 = el2.key();
0184 const auto pt1 = track_table_.get<Pt>(iel1);
0185 const auto pt2 = track_table_.get<Pt>(iel2);
0186 const auto eta1 = track_table_.get<Eta>(iel1);
0187 const auto eta2 = track_table_.get<Eta>(iel2);
0188 const auto phi1 = track_table_.get<Phi>(iel1);
0189 const auto phi2 = track_table_.get<Phi>(iel2);
0190
0191 if (std::abs(eta1 - eta2) > 1) {
0192 dist = -1;
0193 return;
0194 }
0195 if (pt1 > 2 && pt2 > 2 && std::abs(::deltaPhi(phi1, phi2)) > 1) {
0196 dist = -1;
0197 return;
0198 }
0199
0200 const GlobalTrajectoryParameters& gt1 = eventTrackTrajectories_[iel1];
0201 const GlobalTrajectoryParameters& gt2 = eventTrackTrajectories_[iel2];
0202
0203
0204 theMinimum_.calculate(gt1, gt2);
0205
0206
0207 dist = theMinimum_.distance();
0208 LogDebug("PFDisplacedVertexCandidateFinder") << "link iel1=" << iel1 << " iel2=" << iel2 << " dist=" << dist;
0209 crossingPoint = theMinimum_.crossingPoint();
0210
0211 vertexLinkTest = PFDisplacedVertexCandidate::LINKTEST_DCA;
0212
0213
0214 if (dist > dcaCut_) {
0215 dist = -1;
0216 return;
0217 }
0218
0219
0220 double rho2 = crossingPoint.x() * crossingPoint.x() + crossingPoint.y() * crossingPoint.y();
0221
0222 if (rho2 < primaryVertexCut2_) {
0223 dist = -1;
0224 return;
0225 }
0226
0227
0228
0229 double tob_rho_limit2 = 10000;
0230 double tec_z_limit = 270;
0231
0232 if (rho2 > tob_rho_limit2) {
0233 dist = -1;
0234 return;
0235 }
0236 if (std::abs(crossingPoint.z()) > tec_z_limit) {
0237 dist = -1;
0238 return;
0239 }
0240
0241 return;
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267 }
0268
0269
0270
0271 void PFDisplacedVertexCandidateFinder::packLinks(PFDisplacedVertexCandidate& vertexCandidate) {
0272 const vector<TrackBaseRef>& els = vertexCandidate.elements();
0273
0274
0275 for (unsigned i1 = 0; i1 < els.size(); i1++) {
0276 for (unsigned i2 = i1 + 1; i2 < els.size(); i2++) {
0277
0278 if (i1 == i2)
0279 continue;
0280
0281 double dist = -1;
0282 GlobalPoint crossingPoint(0, 0, 0);
0283 PFDisplacedVertexCandidate::VertexLinkTest linktest;
0284
0285 link(els[i1], els[i2], dist, crossingPoint, linktest);
0286 LogDebug("PFDisplacedVertexCandidateFinder")
0287 << "Setting link between elements " << i1 << " key " << els[i1].key() << " and " << i2 << " key "
0288 << els[i2].key() << " of dist =" << dist << " computed from link test " << linktest;
0289
0290 if (dist > -0.5)
0291 vertexCandidate.setLink(i1, i2, dist, crossingPoint, linktest);
0292 }
0293 }
0294 }
0295
0296
0297
0298
0299 GlobalTrajectoryParameters PFDisplacedVertexCandidateFinder::getGlobalTrajectoryParameters(const Track* track) const {
0300 const GlobalPoint position(track->vx(), track->vy(), track->vz());
0301
0302 const GlobalVector momentum(track->momentum().x(), track->momentum().y(), track->momentum().z());
0303
0304 GlobalTrajectoryParameters gtp(position, momentum, track->charge(), magField_);
0305
0306 return gtp;
0307 }
0308
0309
0310 bool PFDisplacedVertexCandidateFinder::goodPtResolution(const TrackBaseRef& trackref) const {
0311 double nChi2 = trackref->normalizedChi2();
0312 double pt = trackref->pt();
0313 double dpt = trackref->ptError();
0314 double dxy = trackref->dxy(pvtx_);
0315
0316 double pt_error = dpt / pt * 100;
0317 double qoverpError = trackref->qoverpError();
0318
0319 LogDebug("PFDisplacedVertexCandidateFinder")
0320 << " PFDisplacedVertexFinder: PFrecTrack->Track Pt= " << pt << " dPt/Pt = " << pt_error << "% nChi2 = " << nChi2;
0321 if (nChi2 > nChi2_max_ || pt < pt_min_ || qoverpError > qoverpError_max_) {
0322 LogDebug("PFDisplacedVertexCandidateFinder") << " PFBlockAlgo: skip badly measured or low pt track"
0323 << " nChi2_cut = " << 5 << " pt_cut = " << 0.2;
0324 return false;
0325 }
0326
0327 if (std::abs(dxy) < dxy_ && pt < pt_min_prim_)
0328 return false;
0329
0330
0331 return true;
0332 }
0333
0334 ostream& operator<<(std::ostream& out, const PFDisplacedVertexCandidateFinder& a) {
0335 if (!out)
0336 return out;
0337
0338 out << "====== Particle Flow Block Algorithm ======= ";
0339 out << endl;
0340 out << "number of unassociated elements : " << a.eventTracks_.size() << endl;
0341 out << endl;
0342
0343 out << " Tracks selection based on " << std::endl;
0344 out << " pvtx_ = " << a.pvtx_ << std::endl;
0345 out << " std::abs(dxy) < " << a.dxy_ << " and pt < " << a.pt_min_prim_ << std::endl;
0346 out << " nChi2 < " << a.nChi2_max_ << " and pt < " << a.pt_min_ << std::endl;
0347
0348 out << endl;
0349
0350 for (PFDisplacedVertexCandidateFinder::IEC ie = a.eventTracks_.begin(); ie != a.eventTracks_.end(); ie++) {
0351 double pt = (*ie).get()->pt();
0352
0353 math::XYZPoint Pi = (*ie).get()->innerPosition();
0354 math::XYZPoint Po = (*ie).get()->outerPosition();
0355
0356 double innermost_radius = sqrt(Pi.x() * Pi.x() + Pi.y() * Pi.y() + Pi.z() * Pi.z());
0357 double outermost_radius = sqrt(Po.x() * Po.x() + Po.y() * Po.y() + Po.z() * Po.z());
0358 double innermost_rho = sqrt(Pi.x() * Pi.x() + Pi.y() * Pi.y());
0359 double outermost_rho = sqrt(Po.x() * Po.x() + Po.y() * Po.y());
0360
0361 out << "ie = " << (*ie).key() << " pt = " << pt << " innermost hit radius = " << innermost_radius
0362 << " rho = " << innermost_rho << " outermost hit radius = " << outermost_radius << " rho = " << outermost_rho
0363 << endl;
0364 }
0365
0366 const std::unique_ptr<reco::PFDisplacedVertexCandidateCollection>& vertexCandidates = std::move(a.vertexCandidates());
0367
0368 if (!vertexCandidates.get()) {
0369 out << "vertexCandidates already transfered" << endl;
0370 } else {
0371 out << "number of vertexCandidates : " << vertexCandidates->size() << endl;
0372 out << endl;
0373
0374 for (PFDisplacedVertexCandidateFinder::IBC ib = vertexCandidates->begin(); ib != vertexCandidates->end(); ib++)
0375 ib->Dump();
0376 }
0377
0378 return out;
0379 }