EgammaElectronAlgos/src/PixelHitMatcher.cc

0001 #include "RecoEgamma/EgammaElectronAlgos/interface/utils.h"
0002 #include "RecoEgamma/EgammaElectronAlgos/interface/PixelHitMatcher.h"
0003
0004 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
0005 #include "DataFormats/GeometryCommonDetAlgo/interface/PerpendicularBoundPlaneBuilder.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0008 #include "TrackingTools/TrajectoryState/interface/ftsFromVertexToPoint.h"
0009 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0011
0012 using namespace reco;
0013 using namespace std;
0014
0015 bool PixelHitMatcher::ForwardMeasurementEstimator::operator()(const GlobalPoint &vprim,
0016                                                               const TrajectoryStateOnSurface &absolute_ts,
0017                                                               const GlobalPoint &absolute_gp,
0018                                                               int charge) const {
0019   GlobalVector ts = absolute_ts.globalParameters().position() - vprim;
0020   GlobalVector gp = absolute_gp - vprim;
0021
0022   float rDiff = gp.perp() - ts.perp();
0023   float rMin = theRMin;
0024   float rMax = theRMax;
0025   float myZ = gp.z();
0026   if ((std::abs(myZ) > 70.f) && (std::abs(myZ) < 170.f)) {
0027     rMin = theRMinI;
0028     rMax = theRMaxI;
0029   }
0030
0031   if ((rDiff >= rMax) | (rDiff <= rMin))
0032     return false;
0033
0034   float phiDiff = -charge * normalizedPhi(gp.barePhi() - ts.barePhi());
0035
0036   return (phiDiff < thePhiMax) & (phiDiff > thePhiMin);
0037 }
0038
0039 bool PixelHitMatcher::BarrelMeasurementEstimator::operator()(const GlobalPoint &vprim,
0040                                                              const TrajectoryStateOnSurface &absolute_ts,
0041                                                              const GlobalPoint &absolute_gp,
0042                                                              int charge) const {
0043   GlobalVector ts = absolute_ts.globalParameters().position() - vprim;
0044   GlobalVector gp = absolute_gp - vprim;
0045
0046   float myZ = gp.z();
0047   float zDiff = myZ - ts.z();
0048   float myZmax = theZMax;
0049   float myZmin = theZMin;
0050   if ((std::abs(myZ) < 30.f) && (gp.perp() > 8.f)) {
0051     myZmax = 0.09f;
0052     myZmin = -0.09f;
0053   }
0054
0055   if ((zDiff >= myZmax) | (zDiff <= myZmin))
0056     return false;
0057
0058   float phiDiff = -charge * normalizedPhi(gp.barePhi() - ts.barePhi());
0059
0060   return (phiDiff < thePhiMax) & (phiDiff > thePhiMin);
0061 }
0062
0063 PixelHitMatcher::PixelHitMatcher(float phi1min,
0064                                  float phi1max,
0065                                  float phi2minB,
0066                                  float phi2maxB,
0067                                  float phi2minF,
0068                                  float phi2maxF,
0069                                  float z2maxB,
0070                                  float r2maxF,
0071                                  float rMaxI,
0072                                  bool useRecoVertex)
0073     :  //zmin1 and zmax1 are dummy at this moment, set from beamspot later
0074       meas1stBLayer{phi1min, phi1max, 0., 0.},
0075       meas2ndBLayer{phi2minB, phi2maxB, -z2maxB, z2maxB},
0076       meas1stFLayer{phi1min, phi1max, 0., 0., -rMaxI, rMaxI},
0077       meas2ndFLayer{phi2minF, phi2maxF, -r2maxF, r2maxF, -rMaxI, rMaxI},
0078       useRecoVertex_(useRecoVertex) {}
0079
0080 void PixelHitMatcher::set1stLayer(float dummyphi1min, float dummyphi1max) {
0081   meas1stBLayer.thePhiMin = dummyphi1min;
0082   meas1stBLayer.thePhiMax = dummyphi1max;
0083   meas1stFLayer.thePhiMin = dummyphi1min;
0084   meas1stFLayer.thePhiMax = dummyphi1max;
0085 }
0086
0087 void PixelHitMatcher::set1stLayerZRange(float zmin1, float zmax1) {
0088   meas1stBLayer.theZMin = zmin1;
0089   meas1stBLayer.theZMax = zmax1;
0090   meas1stFLayer.theRMin = zmin1;
0091   meas1stFLayer.theRMax = zmax1;
0092 }
0093
0094 void PixelHitMatcher::set2ndLayer(float dummyphi2minB, float dummyphi2maxB, float dummyphi2minF, float dummyphi2maxF) {
0095   meas2ndBLayer.thePhiMin = dummyphi2minB;
0096   meas2ndBLayer.thePhiMax = dummyphi2maxB;
0097   meas2ndFLayer.thePhiMin = dummyphi2minF;
0098   meas2ndFLayer.thePhiMax = dummyphi2maxF;
0099 }
0100
0101 void PixelHitMatcher::setES(MagneticField const &magField, TrackerGeometry const &trackerGeometry) {
0102   theMagField = &magField;
0103   theTrackerGeometry = &trackerGeometry;
0104   constexpr float mass = .000511;  // electron propagation
0105   prop1stLayer = std::make_unique<PropagatorWithMaterial>(oppositeToMomentum, mass, theMagField);
0106   prop2ndLayer = std::make_unique<PropagatorWithMaterial>(alongMomentum, mass, theMagField);
0107 }
0108
0109 std::vector<SeedWithInfo> PixelHitMatcher::operator()(const std::vector<const TrajectorySeedCollection *> &seedsV,
0110                                                       const GlobalPoint &xmeas,
0111                                                       const GlobalPoint &vprim,
0112                                                       float energy,
0113                                                       int charge) const {
0114   auto xmeas_r = xmeas.perp();
0115
0116   const float phicut = std::cos(2.5);
0117
0118   auto fts = trackingTools::ftsFromVertexToPoint(*theMagField, xmeas, vprim, energy, charge);
0119   PerpendicularBoundPlaneBuilder bpb;
0120   TrajectoryStateOnSurface tsos(fts, *bpb(fts.position(), fts.momentum()));
0121
0122   std::vector<SeedWithInfo> result;
0123   unsigned int allSeedsSize = 0;
0124   for (auto const sc : seedsV)
0125     allSeedsSize += sc->size();
0126
0127   IntGlobalPointPairUnorderedMap<TrajectoryStateOnSurface> mapTsos2Fast(allSeedsSize);
0128
0129   auto ndets = theTrackerGeometry->dets().size();
0130
0131   int iTsos[ndets];
0132   for (auto &i : iTsos)
0133     i = -1;
0134   std::vector<TrajectoryStateOnSurface> vTsos;
0135   vTsos.reserve(allSeedsSize);
0136
0137   std::vector<GlobalPoint> hitGpMap;
0138   for (const auto seeds : seedsV) {
0139     for (const auto &seed : *seeds) {
0140       hitGpMap.clear();
0141       if (seed.nHits() > 9) {
0142         edm::LogWarning("GsfElectronAlgo|UnexpectedSeed") << "We cannot deal with seeds having more than 9 hits.";
0143         continue;
0144       }
0145
0146       auto const &hits = seed.recHits();
0147       // cache the global points
0148
0149       for (auto const &hit : hits) {
0150         hitGpMap.emplace_back(hit.globalPosition());
0151       }
0152
0153       //iterate on the hits
0154       auto he = hits.end() - 1;
0155       for (auto it1 = hits.begin(); it1 < he; ++it1) {
0156         if (!it1->isValid())
0157           continue;
0158         auto idx1 = std::distance(hits.begin(), it1);
0159         const DetId id1 = it1->geographicalId();
0160         const GeomDet *geomdet1 = it1->det();
0161
0162         auto ix1 = geomdet1->gdetIndex();
0163
0164         /*  VI: this generates regression (other cut is just in phi). in my opinion it is safe and makes sense
0165             auto away = geomdet1->position().basicVector().dot(xmeas.basicVector()) <0;
0166             if (away) continue;
0167         */
0168
0169         const GlobalPoint &hit1Pos = hitGpMap[idx1];
0170         auto dt = hit1Pos.x() * xmeas.x() + hit1Pos.y() * xmeas.y();
0171         if (dt < 0)
0172           continue;
0173         if (dt < phicut * (xmeas_r * hit1Pos.perp()))
0174           continue;
0175
0176         if (iTsos[ix1] < 0) {
0177           iTsos[ix1] = vTsos.size();
0178           vTsos.push_back(prop1stLayer->propagate(tsos, geomdet1->surface()));
0179         }
0180         auto tsos1 = &vTsos[iTsos[ix1]];
0181
0182         if (!tsos1->isValid())
0183           continue;
0184         bool est = id1.subdetId() % 2 ? meas1stBLayer(vprim, *tsos1, hit1Pos, charge)
0185                                       : meas1stFLayer(vprim, *tsos1, hit1Pos, charge);
0186         if (!est)
0187           continue;
0188         EleRelPointPair pp1(hit1Pos, tsos1->globalParameters().position(), vprim);
0189         const math::XYZPoint relHit1Pos(hit1Pos - vprim), relTSOSPos(tsos1->globalParameters().position() - vprim);
0190         const int subDet1 = id1.subdetId();
0191         const float dRz1 = (id1.subdetId() % 2 ? pp1.dZ() : pp1.dPerp());
0192         const float dPhi1 = pp1.dPhi();
0193         // setup our vertex
0194         double zVertex;
0195         if (!useRecoVertex_) {
0196           // we don't know the z vertex position, get it from linear extrapolation
0197           // compute the z vertex from the cluster point and the found pixel hit
0198           const double pxHit1z = hit1Pos.z();
0199           const double pxHit1x = hit1Pos.x();
0200           const double pxHit1y = hit1Pos.y();
0201           const double r1diff =
0202               std::sqrt((pxHit1x - vprim.x()) * (pxHit1x - vprim.x()) + (pxHit1y - vprim.y()) * (pxHit1y - vprim.y()));
0203           const double r2diff =
0204               std::sqrt((xmeas.x() - pxHit1x) * (xmeas.x() - pxHit1x) + (xmeas.y() - pxHit1y) * (xmeas.y() - pxHit1y));
0205           zVertex = pxHit1z - r1diff * (xmeas.z() - pxHit1z) / r2diff;
0206         } else {
0207           // here use rather the reco vertex z position
0208           zVertex = vprim.z();
0209         }
0210         GlobalPoint vertex(vprim.x(), vprim.y(), zVertex);
0211         auto fts2 = trackingTools::ftsFromVertexToPoint(*theMagField, hit1Pos, vertex, energy, charge);
0212         // now find the matching hit
0213         for (auto it2 = it1 + 1; it2 != hits.end(); ++it2) {
0214           if (!it2->isValid())
0215             continue;
0216           auto idx2 = std::distance(hits.begin(), it2);
0217           const DetId id2 = it2->geographicalId();
0218           const GeomDet *geomdet2 = it2->det();
0219           const auto det_key = std::make_pair(geomdet2->gdetIndex(), hit1Pos);
0220           const TrajectoryStateOnSurface *tsos2;
0221           auto tsos2_itr = mapTsos2Fast.find(det_key);
0222           if (tsos2_itr != mapTsos2Fast.end()) {
0223             tsos2 = &(tsos2_itr->second);
0224           } else {
0225             auto empl_result = mapTsos2Fast.emplace(det_key, prop2ndLayer->propagate(fts2, geomdet2->surface()));
0226             tsos2 = &(empl_result.first->second);
0227           }
0228           if (!tsos2->isValid())
0229             continue;
0230           const GlobalPoint &hit2Pos = hitGpMap[idx2];
0231           bool est2 = id2.subdetId() % 2 ? meas2ndBLayer(vertex, *tsos2, hit2Pos, charge)
0232                                          : meas2ndFLayer(vertex, *tsos2, hit2Pos, charge);
0233           if (est2) {
0234             EleRelPointPair pp2(hit2Pos, tsos2->globalParameters().position(), vertex);
0235             const int subDet2 = id2.subdetId();
0236             const float dRz2 = (subDet2 % 2 == 1) ? pp2.dZ() : pp2.dPerp();
0237             const float dPhi2 = pp2.dPhi();
0238             const unsigned char hitsMask = (1 << idx1) | (1 << idx2);
0239             result.push_back({seed, hitsMask, subDet2, dRz2, dPhi2, subDet1, dRz1, dPhi1});
0240           }
0241         }  // inner loop on hits
0242       }  // outer loop on hits
0243     }  // loop on seeds
0244   }  //loop on vector of seeds
0245
0246   return result;
0247 }