File indexing completed on 2023-10-25 09:36:08
0001 #include "CommonTools/RecoAlgos/interface/PrimaryVertexSorting.h"
0002 #include "DataFormats/VertexReco/interface/Vertex.h"
0003 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0004 #include "DataFormats/Math/interface/deltaR.h"
0005 #include "TrackingTools/IPTools/interface/IPTools.h"
0006 #include <fastjet/internal/base.hh>
0007 #include "fastjet/PseudoJet.hh"
0008 #include "fastjet/JetDefinition.hh"
0009 #include "fastjet/ClusterSequence.hh"
0010 #include "fastjet/Selector.hh"
0011 #include "fastjet/PseudoJet.hh"
0012 #include "FWCore/Utilities/interface/isFinite.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 using namespace fastjet;
0016 using namespace std;
0017
0018 float PrimaryVertexSorting::score(const reco::Vertex &pv,
0019 const std::vector<const reco::Candidate *> &cands,
0020 bool useMet) const {
0021 typedef math::XYZTLorentzVector LorentzVector;
0022 float sumPt2 = 0;
0023 float sumEt = 0;
0024 LorentzVector met;
0025 std::vector<fastjet::PseudoJet> fjInputs_;
0026 fjInputs_.clear();
0027 size_t countScale0 = 0;
0028 for (size_t i = 0; i < cands.size(); i++) {
0029 const reco::Candidate *c = cands[i];
0030 float scale = 1.;
0031 if (c->bestTrack() != nullptr) {
0032 if (c->pt() != 0) {
0033 scale = (c->pt() - c->bestTrack()->ptError()) / c->pt();
0034 }
0035 if (edm::isNotFinite(scale)) {
0036 edm::LogWarning("PrimaryVertexSorting") << "Scaling is NAN ignoring this candidate/track" << std::endl;
0037 scale = 0;
0038 }
0039 if (scale < 0) {
0040 scale = 0;
0041 countScale0++;
0042 }
0043 }
0044
0045 int absId = abs(c->pdgId());
0046 if (absId == 13 or absId == 11) {
0047 float pt = c->pt() * scale;
0048 sumPt2 += pt * pt;
0049 met += c->p4() * scale;
0050 sumEt += c->pt() * scale;
0051 } else {
0052 if (scale != 0) {
0053 fjInputs_.push_back(fastjet::PseudoJet(c->px() * scale, c->py() * scale, c->pz() * scale, c->p4().E() * scale));
0054
0055 }
0056 }
0057 }
0058 fastjet::ClusterSequence sequence(fjInputs_, JetDefinition(antikt_algorithm, 0.4));
0059 auto jets = fastjet::sorted_by_pt(sequence.inclusive_jets(0));
0060 for (const auto &pj : jets) {
0061 auto p4 = LorentzVector(pj.px(), pj.py(), pj.pz(), pj.e());
0062 sumPt2 += (p4.pt() * p4.pt()) * 0.8 * 0.8;
0063 met += p4;
0064 sumEt += p4.pt();
0065 }
0066 float metAbove = met.pt() - 2 * sqrt(sumEt);
0067 if (metAbove > 0 and useMet) {
0068 sumPt2 += metAbove * metAbove;
0069 }
0070 if (countScale0 == cands.size())
0071 sumPt2 = countScale0 * 0.01;
0072 return sumPt2;
0073 }