SecondaryVertex/plugins/BtoCharmDecayVertexMerger.cc

0001 // this code is partially taken from BCandidateProducer of L. Wehrli
0002
0003 // first revised prototype version for CMSSW 29.Aug.2012
0004
0005 #include "FWCore/Framework/interface/stream/EDProducer.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/Utilities/interface/InputTag.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010
0011 #include "DataFormats/BTauReco/interface/ParticleMasses.h"
0012 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0016 #include "DataFormats/VertexReco/interface/Vertex.h"
0017 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0018 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0019 #include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
0020
0021 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0022
0023 #include <algorithm>
0024
0025 reco::Candidate::LorentzVector vtxP4(const reco::VertexCompositePtrCandidate &vtx) {
0026   reco::Candidate::LorentzVector sum;
0027   const std::vector<reco::CandidatePtr> &tracks = vtx.daughterPtrVector();
0028
0029   for (std::vector<reco::CandidatePtr>::const_iterator track = tracks.begin(); track != tracks.end(); ++track) {
0030     ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> vec;
0031     vec.SetPx((*track)->px());
0032     vec.SetPy((*track)->py());
0033     vec.SetPz((*track)->pz());
0034     vec.SetM(reco::ParticleMasses::piPlus);
0035     sum += vec;
0036   }
0037   return sum;
0038 }
0039
0040 template <typename VTX>
0041 class BtoCharmDecayVertexMergerT : public edm::stream::EDProducer<> {
0042 public:
0043   BtoCharmDecayVertexMergerT(const edm::ParameterSet &params);
0044
0045   void produce(edm::Event &event, const edm::EventSetup &es) override;
0046
0047   typedef reco::TemplatedSecondaryVertex<VTX> SecondaryVertex;
0048
0049 private:
0050   edm::EDGetTokenT<reco::VertexCollection> token_primaryVertex;
0051   edm::EDGetTokenT<edm::View<VTX>> token_secondaryVertex;
0052   reco::Vertex pv;
0053   // double                                    maxFraction;
0054   // double                                    minSignificance;
0055
0056   double maxDRForUnique, maxvecSumIMCUTForUnique, minCosPAtomerge, maxPtreltomerge;
0057
0058   // a vertex proxy for sorting using stl
0059   struct VertexProxy {
0060     VTX vert;
0061     double invm;
0062     size_t ntracks;
0063   };
0064
0065   // comparison operator for VertexProxy, used in sorting
0066   friend bool operator<(VertexProxy v1, VertexProxy v2) {
0067     if (v1.ntracks > 2 && v2.ntracks < 3)
0068       return true;
0069     if (v1.ntracks < 3 && v2.ntracks > 2)
0070       return false;
0071     return (v1.invm > v2.invm);
0072   }
0073
0074   VertexProxy buildVertexProxy(const VTX &vtx);
0075
0076   void resolveBtoDchain(std::vector<VertexProxy> &coll, unsigned int i, unsigned int k);
0077   GlobalVector flightDirection(const reco::Vertex &v1, const reco::Vertex &v2);
0078   GlobalVector flightDirection(const reco::Vertex &v1, const reco::VertexCompositePtrCandidate &v2);
0079   GlobalVector flightDirection(const reco::VertexCompositePtrCandidate &v1,
0080                                const reco::VertexCompositePtrCandidate &v2);
0081 };
0082
0083 template <typename VTX>
0084 BtoCharmDecayVertexMergerT<VTX>::BtoCharmDecayVertexMergerT(const edm::ParameterSet &params)
0085     : token_primaryVertex(consumes<reco::VertexCollection>(params.getParameter<edm::InputTag>("primaryVertices"))),
0086       token_secondaryVertex(consumes<edm::View<VTX>>(params.getParameter<edm::InputTag>("secondaryVertices"))),
0087       maxDRForUnique(params.getParameter<double>("maxDRUnique")),
0088       maxvecSumIMCUTForUnique(params.getParameter<double>("maxvecSumIMifsmallDRUnique")),
0089       minCosPAtomerge(params.getParameter<double>("minCosPAtomerge")),
0090       maxPtreltomerge(params.getParameter<double>("maxPtreltomerge"))
0091 // maxFraction(params.getParameter<double>("maxFraction")),
0092 // minSignificance(params.getParameter<double>("minSignificance"))
0093 {
0094   produces<std::vector<VTX>>();
0095 }
0096 //------------------------------------
0097 template <typename VTX>
0098 void BtoCharmDecayVertexMergerT<VTX>::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0099   using namespace reco;
0100   // PV
0101   edm::Handle<reco::VertexCollection> PVcoll;
0102   iEvent.getByToken(token_primaryVertex, PVcoll);
0103
0104   if (!PVcoll->empty()) {
0105     const reco::VertexCollection pvc = *(PVcoll.product());
0106     pv = pvc[0];
0107
0108     // get the IVF collection
0109     edm::Handle<edm::View<VTX>> secondaryVertices;
0110     iEvent.getByToken(token_secondaryVertex, secondaryVertices);
0111
0112     //loop over vertices,  fill into own collection for sorting
0113     std::vector<VertexProxy> vertexProxyColl;
0114     for (typename edm::View<VTX>::const_iterator sv = secondaryVertices->begin(); sv != secondaryVertices->end();
0115          ++sv) {
0116       vertexProxyColl.push_back(buildVertexProxy(*sv));
0117     }
0118
0119     // sort the vertices by mass and track multiplicity
0120     sort(vertexProxyColl.begin(), vertexProxyColl.end());
0121
0122     // loop forward over all vertices
0123     for (unsigned int iVtx = 0; iVtx < vertexProxyColl.size(); ++iVtx) {
0124       // nested loop backwards (in order to start from light masses)
0125       // check all vertices against each other for B->D chain
0126       for (unsigned int kVtx = vertexProxyColl.size() - 1; kVtx > iVtx; --kVtx) {
0127         // remove D vertices from the collection and add the tracks to the original one
0128         resolveBtoDchain(vertexProxyColl, iVtx, kVtx);
0129       }
0130     }
0131
0132     // now create new vertex collection and add to event
0133     auto bvertColl = std::make_unique<std::vector<VTX>>();
0134     for (typename std::vector<VertexProxy>::const_iterator it = vertexProxyColl.begin(); it != vertexProxyColl.end();
0135          ++it)
0136       bvertColl->push_back((*it).vert);
0137     iEvent.put(std::move(bvertColl));
0138   } else {
0139     iEvent.put(std::make_unique<std::vector<VTX>>());
0140   }
0141 }
0142 //------------------------------------
0143 template <typename VTX>
0144 GlobalVector BtoCharmDecayVertexMergerT<VTX>::flightDirection(const reco::Vertex &v1, const reco::Vertex &v2) {
0145   return GlobalVector(v2.x() - v1.x(), v2.y() - v1.y(), v2.z() - v1.z());
0146 }
0147
0148 template <typename VTX>
0149 GlobalVector BtoCharmDecayVertexMergerT<VTX>::flightDirection(const reco::Vertex &v1,
0150                                                               const reco::VertexCompositePtrCandidate &v2) {
0151   return GlobalVector(v2.vertex().x() - v1.x(), v2.vertex().y() - v1.y(), v2.vertex().z() - v1.z());
0152 }
0153
0154 template <typename VTX>
0155 GlobalVector BtoCharmDecayVertexMergerT<VTX>::flightDirection(const reco::VertexCompositePtrCandidate &v1,
0156                                                               const reco::VertexCompositePtrCandidate &v2) {
0157   return GlobalVector(
0158       v2.vertex().x() - v1.vertex().x(), v2.vertex().y() - v1.vertex().y(), v2.vertex().z() - v1.vertex().z());
0159 }
0160 //-------------- template specializations --------------------
0161
0162 // -------------- buildVertexProxy ----------------
0163 template <>
0164 BtoCharmDecayVertexMergerT<reco::Vertex>::VertexProxy BtoCharmDecayVertexMergerT<reco::Vertex>::buildVertexProxy(
0165     const reco::Vertex &vtx) {
0166   VertexProxy vtxProxy = {vtx, vtx.p4().M(), vtx.tracksSize()};
0167   return vtxProxy;
0168 }
0169
0170 template <>
0171 BtoCharmDecayVertexMergerT<reco::VertexCompositePtrCandidate>::VertexProxy
0172 BtoCharmDecayVertexMergerT<reco::VertexCompositePtrCandidate>::buildVertexProxy(
0173     const reco::VertexCompositePtrCandidate &vtx) {
0174   VertexProxy vtxProxy = {vtx, vtxP4(vtx).M(), vtx.numberOfSourceCandidatePtrs()};
0175   return vtxProxy;
0176 }
0177
0178 // -------------- resolveBtoDchain ----------------
0179 template <>
0180 void BtoCharmDecayVertexMergerT<reco::Vertex>::resolveBtoDchain(std::vector<VertexProxy> &coll,
0181                                                                 unsigned int i,
0182                                                                 unsigned int k) {
0183   using namespace reco;
0184
0185   GlobalVector momentum1 = GlobalVector(coll[i].vert.p4().X(), coll[i].vert.p4().Y(), coll[i].vert.p4().Z());
0186   GlobalVector momentum2 = GlobalVector(coll[k].vert.p4().X(), coll[k].vert.p4().Y(), coll[k].vert.p4().Z());
0187
0188   SecondaryVertex sv1(pv, coll[i].vert, momentum1, true);
0189   SecondaryVertex sv2(pv, coll[k].vert, momentum2, true);
0190
0191   // find out which one is near and far
0192   SecondaryVertex svNear = sv1;
0193   SecondaryVertex svFar = sv2;
0194   GlobalVector momentumNear = momentum1;
0195   GlobalVector momentumFar = momentum2;
0196
0197   // swap if it is the other way around
0198   if (sv1.dist3d().value() >= sv2.dist3d().value()) {
0199     svNear = sv2;
0200     svFar = sv1;
0201     momentumNear = momentum2;
0202     momentumFar = momentum1;
0203   }
0204   GlobalVector nearToFar = flightDirection(svNear, svFar);
0205   GlobalVector pvToNear = flightDirection(pv, svNear);
0206   GlobalVector pvToFar = flightDirection(pv, svFar);
0207
0208   double cosPA = nearToFar.dot(momentumFar) / momentumFar.mag() / nearToFar.mag();
0209   double cosa = pvToNear.dot(momentumFar) / pvToNear.mag() / momentumFar.mag();
0210   double ptrel = sqrt(1.0 - cosa * cosa) * momentumFar.mag();
0211
0212   double vertexDeltaR = Geom::deltaR(pvToNear, pvToFar);
0213
0214   // create a set of all tracks from both vertices, avoid double counting by using a std::set<>
0215   std::set<reco::TrackRef> trackrefs;
0216   // first vertex
0217   for (reco::Vertex::trackRef_iterator ti = sv1.tracks_begin(); ti != sv1.tracks_end(); ++ti) {
0218     if (sv1.trackWeight(*ti) > 0.5) {
0219       reco::TrackRef t = ti->castTo<reco::TrackRef>();
0220       trackrefs.insert(t);
0221     }
0222   }
0223   // second vertex
0224   for (reco::Vertex::trackRef_iterator ti = sv2.tracks_begin(); ti != sv2.tracks_end(); ++ti) {
0225     if (sv2.trackWeight(*ti) > 0.5) {
0226       reco::TrackRef t = ti->castTo<reco::TrackRef>();
0227       trackrefs.insert(t);
0228     }
0229   }
0230
0231   // now calculate one LorentzVector from the track momenta
0232   ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> mother;
0233   for (std::set<reco::TrackRef>::const_iterator it = trackrefs.begin(); it != trackrefs.end(); ++it) {
0234     ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> temp(
0235         (*it)->px(), (*it)->py(), (*it)->pz(), ParticleMasses::piPlus);
0236     mother += temp;
0237   }
0238
0239   //   check if criteria for merging are fulfilled
0240   if (vertexDeltaR < maxDRForUnique && mother.M() < maxvecSumIMCUTForUnique && cosPA > minCosPAtomerge &&
0241       std::abs(ptrel) < maxPtreltomerge) {
0242     // add tracks of second vertex which are missing to first vertex
0243     // loop over the second
0244     bool bFoundDuplicate = false;
0245     for (reco::Vertex::trackRef_iterator ti = sv2.tracks_begin(); ti != sv2.tracks_end(); ++ti) {
0246       reco::Vertex::trackRef_iterator it = find(sv1.tracks_begin(), sv1.tracks_end(), *ti);
0247       if (it == sv1.tracks_end())
0248         coll[i].vert.add(*ti, sv2.refittedTrack(*ti), sv2.trackWeight(*ti));
0249       else
0250         bFoundDuplicate = true;
0251     }
0252     // in case a duplicate track is found, need to create the full track list in the vertex from scratch because we need to modify the weight.
0253     // the weight must be the larger one, otherwise we may have outliers which are not real outliers
0254
0255     if (bFoundDuplicate) {
0256       // create backup track containers from main vertex
0257       std::vector<TrackBaseRef> tracks_;
0258       std::vector<Track> refittedTracks_;
0259       std::vector<float> weights_;
0260       for (reco::Vertex::trackRef_iterator it = coll[i].vert.tracks_begin(); it != coll[i].vert.tracks_end(); ++it) {
0261         tracks_.push_back(*it);
0262         refittedTracks_.push_back(coll[i].vert.refittedTrack(*it));
0263         weights_.push_back(coll[i].vert.trackWeight(*it));
0264       }
0265       // delete tracks and add all tracks back, and check in which vertex the weight is larger
0266       coll[i].vert.removeTracks();
0267       std::vector<Track>::iterator it2 = refittedTracks_.begin();
0268       std::vector<float>::iterator it3 = weights_.begin();
0269       for (reco::Vertex::trackRef_iterator it = tracks_.begin(); it != tracks_.end(); ++it, ++it2, ++it3) {
0270         float weight = *it3;
0271         float weight2 = sv2.trackWeight(*it);
0272         Track refittedTrackWithLargerWeight = *it2;
0273         if (weight2 > weight) {
0274           weight = weight2;
0275           refittedTrackWithLargerWeight = sv2.refittedTrack(*it);
0276         }
0277         coll[i].vert.add(*it, refittedTrackWithLargerWeight, weight);
0278       }
0279     }
0280
0281     // remove the second vertex from the collection
0282     coll.erase(coll.begin() + k);
0283   }
0284 }
0285
0286 template <>
0287 void BtoCharmDecayVertexMergerT<reco::VertexCompositePtrCandidate>::resolveBtoDchain(std::vector<VertexProxy> &coll,
0288                                                                                      unsigned int i,
0289                                                                                      unsigned int k) {
0290   using namespace reco;
0291
0292   Candidate::LorentzVector vtx1P4 = vtxP4(coll[i].vert);
0293   Candidate::LorentzVector vtx2P4 = vtxP4(coll[k].vert);
0294   GlobalVector momentum1 = GlobalVector(vtx1P4.X(), vtx1P4.Y(), vtx1P4.Z());
0295   GlobalVector momentum2 = GlobalVector(vtx2P4.X(), vtx2P4.Y(), vtx2P4.Z());
0296
0297   SecondaryVertex sv1(pv, coll[i].vert, momentum1, true);
0298   SecondaryVertex sv2(pv, coll[k].vert, momentum2, true);
0299
0300   // find out which one is near and far
0301   SecondaryVertex svNear = sv1;
0302   SecondaryVertex svFar = sv2;
0303   GlobalVector momentumNear = momentum1;
0304   GlobalVector momentumFar = momentum2;
0305
0306   // swap if it is the other way around
0307   if (sv1.dist3d().value() >= sv2.dist3d().value()) {
0308     svNear = sv2;
0309     svFar = sv1;
0310     momentumNear = momentum2;
0311     momentumFar = momentum1;
0312   }
0313   GlobalVector nearToFar = flightDirection(svNear, svFar);
0314   GlobalVector pvToNear = flightDirection(pv, svNear);
0315   GlobalVector pvToFar = flightDirection(pv, svFar);
0316
0317   double cosPA = nearToFar.dot(momentumFar) / momentumFar.mag() / nearToFar.mag();
0318   double cosa = pvToNear.dot(momentumFar) / pvToNear.mag() / momentumFar.mag();
0319   double ptrel = sqrt(1.0 - cosa * cosa) * momentumFar.mag();
0320
0321   double vertexDeltaR = Geom::deltaR(pvToNear, pvToFar);
0322
0323   // create a set of all tracks from both vertices, avoid double counting by using a std::set<>
0324   std::set<reco::CandidatePtr> trackrefs;
0325   // first vertex
0326   for (size_t i = 0; i < sv1.numberOfSourceCandidatePtrs(); ++i)
0327     trackrefs.insert(sv1.daughterPtr(i));
0328   // second vertex
0329   for (size_t i = 0; i < sv2.numberOfSourceCandidatePtrs(); ++i)
0330     trackrefs.insert(sv2.daughterPtr(i));
0331
0332   // now calculate one LorentzVector from the track momenta
0333   Candidate::LorentzVector mother;
0334   for (std::set<reco::CandidatePtr>::const_iterator it = trackrefs.begin(); it != trackrefs.end(); ++it) {
0335     Candidate::LorentzVector temp((*it)->px(), (*it)->py(), (*it)->pz(), ParticleMasses::piPlus);
0336     mother += temp;
0337   }
0338
0339   //   check if criteria for merging are fulfilled
0340   if (vertexDeltaR < maxDRForUnique && mother.M() < maxvecSumIMCUTForUnique && cosPA > minCosPAtomerge &&
0341       std::abs(ptrel) < maxPtreltomerge) {
0342     // add tracks of second vertex which are missing to first vertex
0343     // loop over the second
0344     const std::vector<reco::CandidatePtr> &tracks1 = sv1.daughterPtrVector();
0345     const std::vector<reco::CandidatePtr> &tracks2 = sv2.daughterPtrVector();
0346     for (std::vector<reco::CandidatePtr>::const_iterator ti = tracks2.begin(); ti != tracks2.end(); ++ti) {
0347       std::vector<reco::CandidatePtr>::const_iterator it = find(tracks1.begin(), tracks1.end(), *ti);
0348       if (it == tracks1.end()) {
0349         coll[i].vert.addDaughter(*ti);
0350         coll[i].vert.setP4((*ti)->p4() + coll[i].vert.p4());
0351       }
0352     }
0353
0354     // remove the second vertex from the collection
0355     coll.erase(coll.begin() + k);
0356   }
0357 }
0358
0359 // define specific instances of the templated BtoCharmDecayVertexMergerT class
0360 typedef BtoCharmDecayVertexMergerT<reco::Vertex> BtoCharmDecayVertexMerger;
0361 typedef BtoCharmDecayVertexMergerT<reco::VertexCompositePtrCandidate> CandidateBtoCharmDecayVertexMerger;
0362
0363 // define plugins
0364 DEFINE_FWK_MODULE(BtoCharmDecayVertexMerger);
0365 DEFINE_FWK_MODULE(CandidateBtoCharmDecayVertexMerger);