ParticleFlow/src/PFPileUpAlgo.cc

0001 #include "CommonTools/ParticleFlow/interface/PFPileUpAlgo.h"
0002 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0003 #include "DataFormats/VertexReco/interface/Vertex.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 void PFPileUpAlgo::process(const PFCollection& pfCandidates, const reco::VertexCollection& vertices) {
0007   pfCandidatesFromVtx_.clear();
0008   pfCandidatesFromPU_.clear();
0009
0010   for (unsigned i = 0; i < pfCandidates.size(); i++) {
0011     const reco::PFCandidate& cand = *(pfCandidates[i]);
0012
0013     int ivertex;
0014
0015     switch (cand.particleId()) {
0016       case reco::PFCandidate::h:
0017         ivertex = chargedHadronVertex(vertices, cand);
0018         break;
0019       default:
0020         continue;
0021     }
0022
0023     // no associated vertex, or primary vertex
0024     // not pile-up
0025     if ((ivertex == -1 || ivertex == 0 ||
0026          (fNumOfPUVtxsForCharged_ > 0 && !vertices.empty() && ivertex <= int(fNumOfPUVtxsForCharged_) &&
0027           std::abs(cand.vertex().z() - vertices[0].z()) < fDzCutForChargedFromPUVtxs_))) {
0028       if (verbose_)
0029         std::cout << "VTX " << i << " " << *(pfCandidates[i]) << std::endl;
0030       pfCandidatesFromVtx_.push_back(pfCandidates[i]);
0031     } else {
0032       if (verbose_)
0033         std::cout << "PU  " << i << " " << *(pfCandidates[i]) << std::endl;
0034       // associated to a vertex
0035       pfCandidatesFromPU_.push_back(pfCandidates[i]);
0036     }
0037   }
0038 }
0039
0040 int PFPileUpAlgo::chargedHadronVertex(const reco::VertexCollection& vertices, const reco::PFCandidate& pfcand) const {
0041   auto const& track = pfcand.trackRef();
0042   size_t iVertex = 0;
0043   unsigned int index = 0;
0044   unsigned int nFoundVertex = 0;
0045   float bestweight = 0;
0046   for (auto const& vtx : vertices) {
0047     float w = vtx.trackWeight(track);
0048     //select the vertex for which the track has the highest weight
0049     if (w > bestweight) {
0050       bestweight = w;
0051       iVertex = index;
0052       nFoundVertex++;
0053     }
0054     ++index;
0055   }
0056
0057   if (nFoundVertex > 0) {
0058     if (nFoundVertex != 1)
0059       edm::LogWarning("TrackOnTwoVertex") << "a track is shared by at least two verteces. Used to be an assert";
0060     return iVertex;
0061   }
0062   // no vertex found with this track.
0063
0064   // optional: as a secondary solution, associate the closest vertex in z
0065   if (checkClosestZVertex_) {
0066     double dzmin = 10000;
0067     double ztrack = pfcand.vertex().z();
0068     bool foundVertex = false;
0069     index = 0;
0070     for (auto iv = vertices.begin(); iv != vertices.end(); ++iv, ++index) {
0071       double dz = fabs(ztrack - iv->z());
0072       if (dz < dzmin) {
0073         dzmin = dz;
0074         iVertex = index;
0075         foundVertex = true;
0076       }
0077     }
0078
0079     if (foundVertex)
0080       return iVertex;
0081   }
0082
0083   return -1;
0084 }