RecoAlgos/src/PrimaryVertexSorting.cc

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) {  // otherwise, what is the point to cluster zeroes
0053         fjInputs_.push_back(fastjet::PseudoJet(c->px() * scale, c->py() * scale, c->pz() * scale, c->p4().E() * scale));
0054         //      fjInputs_.back().set_user_index(i);
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;  //leave some epsilon value to sort vertices with unknown pt
0072   return sumPt2;
0073 }