CMS 3D CMS Logo

BtoCharmDecayVertexMerger.cc
Go to the documentation of this file.
1 // this code is partially taken from BCandidateProducer of L. Wehrli
2 
3 // first revised prototype version for CMSSW 29.Aug.2012
4 
10 
20 
22 
23 #include <algorithm>
24 
27  const std::vector<reco::CandidatePtr> &tracks = vtx.daughterPtrVector();
28 
29  for (std::vector<reco::CandidatePtr>::const_iterator track = tracks.begin(); track != tracks.end(); ++track) {
30  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> vec;
31  vec.SetPx((*track)->px());
32  vec.SetPy((*track)->py());
33  vec.SetPz((*track)->pz());
35  sum += vec;
36  }
37  return sum;
38 }
39 
40 template <typename VTX>
42 public:
44 
45  void produce(edm::Event &event, const edm::EventSetup &es) override;
46 
48 
49 private:
53  // double maxFraction;
54  // double minSignificance;
55 
57 
58  // a vertex proxy for sorting using stl
59  struct VertexProxy {
60  VTX vert;
61  double invm;
62  size_t ntracks;
63  };
64 
65  // comparison operator for VertexProxy, used in sorting
66  friend bool operator<(VertexProxy v1, VertexProxy v2) {
67  if (v1.ntracks > 2 && v2.ntracks < 3)
68  return true;
69  if (v1.ntracks < 3 && v2.ntracks > 2)
70  return false;
71  return (v1.invm > v2.invm);
72  }
73 
74  VertexProxy buildVertexProxy(const VTX &vtx);
75 
76  void resolveBtoDchain(std::vector<VertexProxy> &coll, unsigned int i, unsigned int k);
81 };
82 
83 template <typename VTX>
85  : token_primaryVertex(consumes<reco::VertexCollection>(params.getParameter<edm::InputTag>("primaryVertices"))),
86  token_secondaryVertex(consumes<edm::View<VTX>>(params.getParameter<edm::InputTag>("secondaryVertices"))),
87  maxDRForUnique(params.getParameter<double>("maxDRUnique")),
88  maxvecSumIMCUTForUnique(params.getParameter<double>("maxvecSumIMifsmallDRUnique")),
89  minCosPAtomerge(params.getParameter<double>("minCosPAtomerge")),
90  maxPtreltomerge(params.getParameter<double>("maxPtreltomerge"))
91 // maxFraction(params.getParameter<double>("maxFraction")),
92 // minSignificance(params.getParameter<double>("minSignificance"))
93 {
94  produces<std::vector<VTX>>();
95 }
96 //------------------------------------
97 template <typename VTX>
99  using namespace reco;
100  // PV
102  iEvent.getByToken(token_primaryVertex, PVcoll);
103 
104  if (!PVcoll->empty()) {
105  const reco::VertexCollection pvc = *(PVcoll.product());
106  pv = pvc[0];
107 
108  // get the IVF collection
110  iEvent.getByToken(token_secondaryVertex, secondaryVertices);
111 
112  //loop over vertices, fill into own collection for sorting
113  std::vector<VertexProxy> vertexProxyColl;
114  for (typename edm::View<VTX>::const_iterator sv = secondaryVertices->begin(); sv != secondaryVertices->end();
115  ++sv) {
116  vertexProxyColl.push_back(buildVertexProxy(*sv));
117  }
118 
119  // sort the vertices by mass and track multiplicity
120  sort(vertexProxyColl.begin(), vertexProxyColl.end());
121 
122  // loop forward over all vertices
123  for (unsigned int iVtx = 0; iVtx < vertexProxyColl.size(); ++iVtx) {
124  // nested loop backwards (in order to start from light masses)
125  // check all vertices against each other for B->D chain
126  for (unsigned int kVtx = vertexProxyColl.size() - 1; kVtx > iVtx; --kVtx) {
127  // remove D vertices from the collection and add the tracks to the original one
128  resolveBtoDchain(vertexProxyColl, iVtx, kVtx);
129  }
130  }
131 
132  // now create new vertex collection and add to event
133  auto bvertColl = std::make_unique<std::vector<VTX>>();
134  for (typename std::vector<VertexProxy>::const_iterator it = vertexProxyColl.begin(); it != vertexProxyColl.end();
135  ++it)
136  bvertColl->push_back((*it).vert);
137  iEvent.put(std::move(bvertColl));
138  } else {
139  iEvent.put(std::make_unique<std::vector<VTX>>());
140  }
141 }
142 //------------------------------------
143 template <typename VTX>
145  return GlobalVector(v2.x() - v1.x(), v2.y() - v1.y(), v2.z() - v1.z());
146 }
147 
148 template <typename VTX>
151  return GlobalVector(v2.vertex().x() - v1.x(), v2.vertex().y() - v1.y(), v2.vertex().z() - v1.z());
152 }
153 
154 template <typename VTX>
157  return GlobalVector(
158  v2.vertex().x() - v1.vertex().x(), v2.vertex().y() - v1.vertex().y(), v2.vertex().z() - v1.vertex().z());
159 }
160 //-------------- template specializations --------------------
161 
162 // -------------- buildVertexProxy ----------------
163 template <>
165  const reco::Vertex &vtx) {
166  VertexProxy vtxProxy = {vtx, vtx.p4().M(), vtx.tracksSize()};
167  return vtxProxy;
168 }
169 
170 template <>
174  VertexProxy vtxProxy = {vtx, vtxP4(vtx).M(), vtx.numberOfSourceCandidatePtrs()};
175  return vtxProxy;
176 }
177 
178 // -------------- resolveBtoDchain ----------------
179 template <>
181  unsigned int i,
182  unsigned int k) {
183  using namespace reco;
184 
185  GlobalVector momentum1 = GlobalVector(coll[i].vert.p4().X(), coll[i].vert.p4().Y(), coll[i].vert.p4().Z());
186  GlobalVector momentum2 = GlobalVector(coll[k].vert.p4().X(), coll[k].vert.p4().Y(), coll[k].vert.p4().Z());
187 
188  SecondaryVertex sv1(pv, coll[i].vert, momentum1, true);
189  SecondaryVertex sv2(pv, coll[k].vert, momentum2, true);
190 
191  // find out which one is near and far
192  SecondaryVertex svNear = sv1;
193  SecondaryVertex svFar = sv2;
194  GlobalVector momentumNear = momentum1;
195  GlobalVector momentumFar = momentum2;
196 
197  // swap if it is the other way around
198  if (sv1.dist3d().value() >= sv2.dist3d().value()) {
199  svNear = sv2;
200  svFar = sv1;
201  momentumNear = momentum2;
202  momentumFar = momentum1;
203  }
204  GlobalVector nearToFar = flightDirection(svNear, svFar);
205  GlobalVector pvToNear = flightDirection(pv, svNear);
206  GlobalVector pvToFar = flightDirection(pv, svFar);
207 
208  double cosPA = nearToFar.dot(momentumFar) / momentumFar.mag() / nearToFar.mag();
209  double cosa = pvToNear.dot(momentumFar) / pvToNear.mag() / momentumFar.mag();
210  double ptrel = sqrt(1.0 - cosa * cosa) * momentumFar.mag();
211 
212  double vertexDeltaR = Geom::deltaR(pvToNear, pvToFar);
213 
214  // create a set of all tracks from both vertices, avoid double counting by using a std::set<>
215  std::set<reco::TrackRef> trackrefs;
216  // first vertex
217  for (reco::Vertex::trackRef_iterator ti = sv1.tracks_begin(); ti != sv1.tracks_end(); ++ti) {
218  if (sv1.trackWeight(*ti) > 0.5) {
219  reco::TrackRef t = ti->castTo<reco::TrackRef>();
220  trackrefs.insert(t);
221  }
222  }
223  // second vertex
224  for (reco::Vertex::trackRef_iterator ti = sv2.tracks_begin(); ti != sv2.tracks_end(); ++ti) {
225  if (sv2.trackWeight(*ti) > 0.5) {
226  reco::TrackRef t = ti->castTo<reco::TrackRef>();
227  trackrefs.insert(t);
228  }
229  }
230 
231  // now calculate one LorentzVector from the track momenta
232  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> mother;
233  for (std::set<reco::TrackRef>::const_iterator it = trackrefs.begin(); it != trackrefs.end(); ++it) {
234  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> temp(
235  (*it)->px(), (*it)->py(), (*it)->pz(), ParticleMasses::piPlus);
236  mother += temp;
237  }
238 
239  // check if criteria for merging are fulfilled
240  if (vertexDeltaR < maxDRForUnique && mother.M() < maxvecSumIMCUTForUnique && cosPA > minCosPAtomerge &&
241  std::abs(ptrel) < maxPtreltomerge) {
242  // add tracks of second vertex which are missing to first vertex
243  // loop over the second
244  bool bFoundDuplicate = false;
245  for (reco::Vertex::trackRef_iterator ti = sv2.tracks_begin(); ti != sv2.tracks_end(); ++ti) {
246  reco::Vertex::trackRef_iterator it = find(sv1.tracks_begin(), sv1.tracks_end(), *ti);
247  if (it == sv1.tracks_end())
248  coll[i].vert.add(*ti, sv2.refittedTrack(*ti), sv2.trackWeight(*ti));
249  else
250  bFoundDuplicate = true;
251  }
252  // 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.
253  // the weight must be the larger one, otherwise we may have outliers which are not real outliers
254 
255  if (bFoundDuplicate) {
256  // create backup track containers from main vertex
257  std::vector<TrackBaseRef> tracks_;
258  std::vector<Track> refittedTracks_;
259  std::vector<float> weights_;
260  for (reco::Vertex::trackRef_iterator it = coll[i].vert.tracks_begin(); it != coll[i].vert.tracks_end(); ++it) {
261  tracks_.push_back(*it);
262  refittedTracks_.push_back(coll[i].vert.refittedTrack(*it));
263  weights_.push_back(coll[i].vert.trackWeight(*it));
264  }
265  // delete tracks and add all tracks back, and check in which vertex the weight is larger
266  coll[i].vert.removeTracks();
267  std::vector<Track>::iterator it2 = refittedTracks_.begin();
268  std::vector<float>::iterator it3 = weights_.begin();
269  for (reco::Vertex::trackRef_iterator it = tracks_.begin(); it != tracks_.end(); ++it, ++it2, ++it3) {
270  float weight = *it3;
271  float weight2 = sv2.trackWeight(*it);
272  Track refittedTrackWithLargerWeight = *it2;
273  if (weight2 > weight) {
274  weight = weight2;
275  refittedTrackWithLargerWeight = sv2.refittedTrack(*it);
276  }
277  coll[i].vert.add(*it, refittedTrackWithLargerWeight, weight);
278  }
279  }
280 
281  // remove the second vertex from the collection
282  coll.erase(coll.begin() + k);
283  }
284 }
285 
286 template <>
288  unsigned int i,
289  unsigned int k) {
290  using namespace reco;
291 
292  Candidate::LorentzVector vtx1P4 = vtxP4(coll[i].vert);
293  Candidate::LorentzVector vtx2P4 = vtxP4(coll[k].vert);
294  GlobalVector momentum1 = GlobalVector(vtx1P4.X(), vtx1P4.Y(), vtx1P4.Z());
295  GlobalVector momentum2 = GlobalVector(vtx2P4.X(), vtx2P4.Y(), vtx2P4.Z());
296 
297  SecondaryVertex sv1(pv, coll[i].vert, momentum1, true);
298  SecondaryVertex sv2(pv, coll[k].vert, momentum2, true);
299 
300  // find out which one is near and far
301  SecondaryVertex svNear = sv1;
302  SecondaryVertex svFar = sv2;
303  GlobalVector momentumNear = momentum1;
304  GlobalVector momentumFar = momentum2;
305 
306  // swap if it is the other way around
307  if (sv1.dist3d().value() >= sv2.dist3d().value()) {
308  svNear = sv2;
309  svFar = sv1;
310  momentumNear = momentum2;
311  momentumFar = momentum1;
312  }
313  GlobalVector nearToFar = flightDirection(svNear, svFar);
314  GlobalVector pvToNear = flightDirection(pv, svNear);
315  GlobalVector pvToFar = flightDirection(pv, svFar);
316 
317  double cosPA = nearToFar.dot(momentumFar) / momentumFar.mag() / nearToFar.mag();
318  double cosa = pvToNear.dot(momentumFar) / pvToNear.mag() / momentumFar.mag();
319  double ptrel = sqrt(1.0 - cosa * cosa) * momentumFar.mag();
320 
321  double vertexDeltaR = Geom::deltaR(pvToNear, pvToFar);
322 
323  // create a set of all tracks from both vertices, avoid double counting by using a std::set<>
324  std::set<reco::CandidatePtr> trackrefs;
325  // first vertex
326  for (size_t i = 0; i < sv1.numberOfSourceCandidatePtrs(); ++i)
327  trackrefs.insert(sv1.daughterPtr(i));
328  // second vertex
329  for (size_t i = 0; i < sv2.numberOfSourceCandidatePtrs(); ++i)
330  trackrefs.insert(sv2.daughterPtr(i));
331 
332  // now calculate one LorentzVector from the track momenta
334  for (std::set<reco::CandidatePtr>::const_iterator it = trackrefs.begin(); it != trackrefs.end(); ++it) {
335  Candidate::LorentzVector temp((*it)->px(), (*it)->py(), (*it)->pz(), ParticleMasses::piPlus);
336  mother += temp;
337  }
338 
339  // check if criteria for merging are fulfilled
340  if (vertexDeltaR < maxDRForUnique && mother.M() < maxvecSumIMCUTForUnique && cosPA > minCosPAtomerge &&
341  std::abs(ptrel) < maxPtreltomerge) {
342  // add tracks of second vertex which are missing to first vertex
343  // loop over the second
344  const std::vector<reco::CandidatePtr> &tracks1 = sv1.daughterPtrVector();
345  const std::vector<reco::CandidatePtr> &tracks2 = sv2.daughterPtrVector();
346  for (std::vector<reco::CandidatePtr>::const_iterator ti = tracks2.begin(); ti != tracks2.end(); ++ti) {
347  std::vector<reco::CandidatePtr>::const_iterator it = find(tracks1.begin(), tracks1.end(), *ti);
348  if (it == tracks1.end()) {
349  coll[i].vert.addDaughter(*ti);
350  coll[i].vert.setP4((*ti)->p4() + coll[i].vert.p4());
351  }
352  }
353 
354  // remove the second vertex from the collection
355  coll.erase(coll.begin() + k);
356  }
357 }
358 
359 // define specific instances of the templated BtoCharmDecayVertexMergerT class
362 
363 // define plugins
const double piPlus
Definition: ParticleMasses.h:9
reco::TemplatedSecondaryVertex< VTX > SecondaryVertex
double z() const
z coordinate
Definition: Vertex.h:133
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:99
T const * product() const
Definition: Handle.h:70
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
Definition: weight.py:1
void produce(edm::Event &event, const edm::EventSetup &es) override
const Point & vertex() const override
vertex position (overwritten by PF...)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
std::vector< Vertex > VertexCollection
Definition: Vertex.h:31
BtoCharmDecayVertexMergerT< reco::Vertex > BtoCharmDecayVertexMerger
VertexProxy buildVertexProxy(const VTX &vtx)
int iEvent
Definition: GenABIO.cc:224
GlobalVector flightDirection(const reco::Vertex &pv, const reco::Vertex &sv)
friend bool operator<(VertexProxy v1, VertexProxy v2)
T sqrt(T t)
Definition: SSEVec.h:19
T mag() const
Definition: PV3DBase.h:64
def pv(vc)
Definition: MetAnalyzer.py:7
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double x() const
x coordinate
Definition: Vertex.h:129
double y() const
y coordinate
Definition: Vertex.h:131
auto const & tracks
cannot be loose
BtoCharmDecayVertexMergerT(const edm::ParameterSet &params)
void resolveBtoDchain(std::vector< VertexProxy > &coll, unsigned int i, unsigned int k)
GlobalVector flightDirection(const reco::Vertex &v1, const reco::Vertex &v2)
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:36
double value() const
Definition: Measurement1D.h:25
BtoCharmDecayVertexMergerT< reco::VertexCompositePtrCandidate > CandidateBtoCharmDecayVertexMerger
reco::Candidate::LorentzVector vtxP4(const reco::VertexCompositePtrCandidate &vtx)
fixed size matrix
HLT enums.
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:86
edm::EDGetTokenT< edm::View< VTX > > token_secondaryVertex
edm::EDGetTokenT< reco::VertexCollection > token_primaryVertex
def move(src, dest)
Definition: eostools.py:511
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector<TrackRef>
Definition: Vertex.h:38
Definition: event.py:1
Global3DVector GlobalVector
Definition: GlobalVector.h:10