File indexing completed on 2024-09-07 04:37:27
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 :
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;
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
0148
0149 for (auto const &hit : hits) {
0150 hitGpMap.emplace_back(hit.globalPosition());
0151 }
0152
0153
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
0165
0166
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
0194 double zVertex;
0195 if (!useRecoVertex_) {
0196
0197
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
0208 zVertex = vprim.z();
0209 }
0210 GlobalPoint vertex(vprim.x(), vprim.y(), zVertex);
0211 auto fts2 = trackingTools::ftsFromVertexToPoint(*theMagField, hit1Pos, vertex, energy, charge);
0212
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 }
0242 }
0243 }
0244 }
0245
0246 return result;
0247 }