EgammaPhotonAlgos/src/ConversionHitChecker.cc

0001 #include <iostream>
0002 #include <vector>
0003 #include <memory>
0004 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionHitChecker.h"
0005 // Framework
0006 //
0007
0008 std::pair<uint8_t, Measurement1DFloat> ConversionHitChecker::nHitsBeforeVtx(const reco::TrackExtra &track,
0009                                                                             const reco::Vertex &vtx,
0010                                                                             float sigmaTolerance) const {
0011   // track hits are always inout
0012
0013   GlobalPoint vtxPos(vtx.x(), vtx.y(), vtx.z());
0014
0015   auto const &trajParams = track.trajParams();
0016
0017   //iterate inside out, when distance to vertex starts increasing, we are at the closest hit
0018   // the first (and last, btw) hit is always valid... (apparntly not..., conversion is different????)
0019   TrackingRecHit const *recHit = *track.recHits().begin();
0020   unsigned int closest = 0;
0021   for (auto const &hit : track.recHits()) {
0022     if (hit->isValid()) {
0023       recHit = hit;
0024       break;
0025     }
0026     ++closest;
0027   }
0028   auto globalPosition = recHit->surface()->toGlobal(trajParams[closest].position());
0029   auto distance2 = (vtxPos - globalPosition).mag2();
0030   int nhits = 1;
0031   for (unsigned int h = closest + 1; h < track.recHitsSize(); ++h) {
0032     //check if next valid hit is farther away from vertex than existing closest
0033     auto nextHit = track.recHit(h);
0034     if (!nextHit->isValid())
0035       continue;
0036     globalPosition = nextHit->surface()->toGlobal(trajParams[h].position());
0037     auto nextDistance2 = (vtxPos - globalPosition).mag2();
0038     if (nextDistance2 > distance2)
0039       break;
0040
0041     distance2 = nextDistance2;
0042     ++nhits;
0043     closest = h;
0044   }
0045
0046   //compute signed decaylength significance for closest hit and check if it is before the vertex
0047   //if not then we need to subtract it from the count of hits before the vertex, since it has been implicitly included
0048   recHit = track.recHit(closest).get();
0049   auto momDir = recHit->surface()->toGlobal(trajParams[closest].direction());
0050   globalPosition = recHit->surface()->toGlobal(trajParams[closest].position());
0051   float decayLengthHitToVtx = (vtxPos - globalPosition).dot(momDir);
0052
0053   AlgebraicVector3 j;
0054   j[0] = momDir.x();
0055   j[1] = momDir.y();
0056   j[2] = momDir.z();
0057   float vertexError2 = ROOT::Math::Similarity(j, vtx.covariance());
0058   auto decayLenError = std::sqrt(vertexError2);
0059
0060   Measurement1DFloat decayLength(decayLengthHitToVtx, decayLenError);
0061
0062   if (decayLength.significance() <
0063       sigmaTolerance) {  //decay length is not (significantly) positive, so hit is consistent with the vertex position or late
0064                          //subtract it from wrong hits count
0065     --nhits;
0066   }
0067
0068   return std::pair<unsigned int, Measurement1DFloat>(nhits, decayLength);
0069 }
0070
0071 uint8_t ConversionHitChecker::nSharedHits(const reco::Track &trk1, const reco::Track &trk2) const {
0072   uint8_t nShared = 0;
0073
0074   for (trackingRecHit_iterator iHit1 = trk1.recHitsBegin(); iHit1 != trk1.recHitsEnd(); ++iHit1) {
0075     const TrackingRecHit *hit1 = (*iHit1);
0076     if (hit1->isValid()) {
0077       for (trackingRecHit_iterator iHit2 = trk2.recHitsBegin(); iHit2 != trk2.recHitsEnd(); ++iHit2) {
0078         const TrackingRecHit *hit2 = (*iHit2);
0079         if (hit2->isValid() && hit1->sharesInput(hit2, TrackingRecHit::some)) {
0080           ++nShared;
0081         }
0082       }
0083     }
0084   }
0085
0086   return nShared;
0087 }