TkTrackingRegions/plugins/AreaSeededTrackingRegionsBuilder.cc

0001 #include "AreaSeededTrackingRegionsBuilder.h"
0002
0003 #include "FWCore/Utilities/interface/InputTag.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "DataFormats/Math/interface/PtEtaPhiMass.h"
0006 #include "DataFormats/Math/interface/deltaPhi.h"
0007 #include "DataFormats/Math/interface/normalizedPhi.h"
0008 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0009
0010 #include <array>
0011 #include <limits>
0012
0013 namespace {
0014   float perp2(const std::array<float, 2>& a) { return a[0] * a[0] + a[1] * a[1]; }
0015 }  // namespace
0016
0017 AreaSeededTrackingRegionsBuilder::AreaSeededTrackingRegionsBuilder(const edm::ParameterSet& regPSet,
0018                                                                    edm::ConsumesCollector& iC)
0019     : candidates_(regPSet, iC), token_field(iC.esConsumes()) {
0020   m_extraPhi = regPSet.getParameter<double>("extraPhi");
0021   m_extraEta = regPSet.getParameter<double>("extraEta");
0022
0023   // RectangularEtaPhiTrackingRegion parameters:
0024   m_ptMin = regPSet.getParameter<double>("ptMin");
0025   m_originRadius = regPSet.getParameter<double>("originRadius");
0026   m_precise = regPSet.getParameter<bool>("precise");
0027   m_whereToUseMeasurementTracker = RectangularEtaPhiTrackingRegion::stringToUseMeasurementTracker(
0028       regPSet.getParameter<std::string>("whereToUseMeasurementTracker"));
0029   if (m_whereToUseMeasurementTracker != RectangularEtaPhiTrackingRegion::UseMeasurementTracker::kNever) {
0030     token_measurementTracker =
0031         iC.consumes<MeasurementTrackerEvent>(regPSet.getParameter<edm::InputTag>("measurementTrackerName"));
0032   }
0033   m_searchOpt = regPSet.getParameter<bool>("searchOpt");
0034   if (m_precise) {
0035     token_msmaker = iC.esConsumes();
0036   }
0037 }
0038
0039 void AreaSeededTrackingRegionsBuilder::fillDescriptions(edm::ParameterSetDescription& desc) {
0040   desc.add<double>("extraPhi", 0.);
0041   desc.add<double>("extraEta", 0.);
0042
0043   desc.add<double>("ptMin", 0.9);
0044   desc.add<double>("originRadius", 0.2);
0045   desc.add<bool>("precise", true);
0046
0047   desc.add<std::string>("whereToUseMeasurementTracker", "Never");
0048   desc.add<edm::InputTag>("measurementTrackerName", edm::InputTag(""));
0049   TrackingSeedCandidates::fillDescriptions(desc);
0050   desc.add<bool>("searchOpt", false);
0051 }
0052
0053 AreaSeededTrackingRegionsBuilder::Builder AreaSeededTrackingRegionsBuilder::beginEvent(
0054     const edm::Event& e, const edm::EventSetup& es) const {
0055   const auto& field = es.getData(token_field);
0056   const MultipleScatteringParametrisationMaker* msmaker = nullptr;
0057   if (m_precise) {
0058     msmaker = &es.getData(token_msmaker);
0059   }
0060   auto builder = Builder(this, &field, msmaker);
0061
0062   if (!token_measurementTracker.isUninitialized()) {
0063     edm::Handle<MeasurementTrackerEvent> hmte;
0064     e.getByToken(token_measurementTracker, hmte);
0065     builder.setMeasurementTracker(hmte.product());
0066   }
0067   builder.setCandidates((candidates_.objects(e)));
0068   return builder;
0069 }
0070
0071 std::vector<std::unique_ptr<TrackingRegion> > AreaSeededTrackingRegionsBuilder::Builder::regions(
0072     const Origins& origins, const std::vector<Area>& areas) const {
0073   std::vector<std::unique_ptr<TrackingRegion> > result;
0074
0075   // create tracking regions in directions of the points of interest
0076   int n_regions = 0;
0077   for (const auto& origin : origins) {
0078     auto reg = region(origin, areas);
0079     if (!reg)
0080       continue;
0081     result.push_back(std::move(reg));
0082     ++n_regions;
0083   }
0084   LogDebug("AreaSeededTrackingRegionsBuilder") << "produced " << n_regions << " regions";
0085
0086   return result;
0087 }
0088
0089 std::unique_ptr<TrackingRegion> AreaSeededTrackingRegionsBuilder::Builder::region(
0090     const Origin& origin, const std::vector<Area>& areas) const {
0091   return regionImpl(origin, areas);
0092 }
0093 std::unique_ptr<TrackingRegion> AreaSeededTrackingRegionsBuilder::Builder::region(
0094     const Origin& origin, const edm::VecArray<Area, 2>& areas) const {
0095   return regionImpl(origin, areas);
0096 }
0097
0098 template <typename T>
0099 std::unique_ptr<TrackingRegion> AreaSeededTrackingRegionsBuilder::Builder::regionImpl(const Origin& origin,
0100                                                                                       const T& areas) const {
0101   float minEta = std::numeric_limits<float>::max(), maxEta = std::numeric_limits<float>::lowest();
0102   float minPhi = std::numeric_limits<float>::max(), maxPhi = std::numeric_limits<float>::lowest();
0103
0104   const auto& orig = origin.first;
0105
0106   LogDebug("AreaSeededTrackingRegionsProducer") << "Origin x,y,z " << orig.x() << "," << orig.y() << "," << orig.z();
0107
0108   auto vecFromOrigPlusRadius = [&](const std::array<float, 2>& vec) {
0109     const auto invlen = 1.f / std::sqrt(perp2(vec));
0110     const auto tmp = (1.f - m_conf->m_originRadius * invlen);
0111     return std::array<float, 2>{{vec[0] * tmp, vec[1] * tmp}};
0112   };
0113   auto tangentVec = [&](const std::array<float, 2>& vec, int sign) {
0114     const auto d = std::sqrt(perp2(vec));
0115     const auto r = m_conf->m_originRadius;
0116     const auto tmp = r / std::sqrt(d * d - r * r);
0117     // sign+ for counterclockwise, sign- for clockwise
0118     const auto orthvec = sign > 0 ? std::array<float, 2>{{-vec[1] * tmp, vec[0] * tmp}}
0119                                   : std::array<float, 2>{{vec[1] * tmp, -vec[0] * tmp}};
0120     return std::array<float, 2>{{vec[0] - orthvec[0], vec[1] - orthvec[1]}};
0121   };
0122
0123   for (const auto& area : areas) {
0124     // straight line assumption is conservative, accounding for
0125     // low-pT bending would only tighten the eta-phi window
0126     LogTrace("AreaSeededTrackingRegionsBuilder")
0127         << " area x,y points " << area.x_rmin_phimin() << "," << area.y_rmin_phimin() << " "  // go around
0128         << area.x_rmin_phimax() << "," << area.y_rmin_phimax() << " " << area.x_rmax_phimax() << ","
0129         << area.y_rmax_phimax() << " " << area.x_rmax_phimin() << "," << area.y_rmax_phimin() << " "
0130         << "z " << area.zmin() << "," << area.zmax();
0131
0132     // some common variables
0133     const float x_rmin_phimin = area.x_rmin_phimin() - orig.x();
0134     const float x_rmin_phimax = area.x_rmin_phimax() - orig.x();
0135     const float y_rmin_phimin = area.y_rmin_phimin() - orig.y();
0136     const float y_rmin_phimax = area.y_rmin_phimax() - orig.y();
0137
0138     const std::array<float, 2> p_rmin_phimin = {{x_rmin_phimin, y_rmin_phimin}};
0139     const std::array<float, 2> p_rmin_phimax = {{x_rmin_phimax, y_rmin_phimax}};
0140     const std::array<float, 2> p_rmax_phimin = {{area.x_rmax_phimin() - orig.x(), area.y_rmax_phimin() - orig.y()}};
0141     const std::array<float, 2> p_rmax_phimax = {{area.x_rmax_phimax() - orig.x(), area.y_rmax_phimax() - orig.y()}};
0142
0143     // eta
0144     {
0145       // find which of the two rmin points is closer to the origin
0146       //
0147       // closest point for p_rmin, farthest point for p_rmax
0148       const std::array<float, 2> p_rmin = perp2(p_rmin_phimin) < perp2(p_rmin_phimax) ? p_rmin_phimin : p_rmin_phimax;
0149       const std::array<float, 2> p_rmax = perp2(p_rmax_phimin) > perp2(p_rmax_phimax) ? p_rmax_phimin : p_rmax_phimax;
0150
0151       // then calculate the xy vector from the point closest to p_rmin on
0152       // the (origin,originRadius) circle to the p_rmin
0153       // this will maximize the eta window
0154       const auto p_min = vecFromOrigPlusRadius(p_rmin);
0155       const auto p_max = vecFromOrigPlusRadius(p_rmax);
0156
0157       // then we calculate the maximal eta window
0158       const auto etamin = std::min(etaFromXYZ(p_min[0], p_min[1], area.zmin() - (orig.z() + origin.second)),
0159                                    etaFromXYZ(p_max[0], p_max[1], area.zmin() - (orig.z() + origin.second)));
0160       const auto etamax = std::max(etaFromXYZ(p_min[0], p_min[1], area.zmax() - (orig.z() - origin.second)),
0161                                    etaFromXYZ(p_max[0], p_max[1], area.zmax() - (orig.z() - origin.second)));
0162
0163       LogTrace("AreaSeededTrackingRegionBuilder") << "  eta min,max " << etamin << "," << etamax;
0164
0165       minEta = std::min(minEta, etamin);
0166       maxEta = std::max(maxEta, etamax);
0167     }
0168
0169     // phi
0170     {
0171       // For phi we construct the tangent lines of (origin,
0172       // originRadius) that go though each of the 4 points (in xy
0173       // plane) of the area. Easiest is to make a vector orthogonal to
0174       // origin->point vector which has a length of
0175       //
0176       // length = r*d/std::sqrt(d*d-r*r)
0177       //
0178       // where r is the originRadius and d is the distance from origin
0179       // to the point (to derive draw the situation and start with
0180       // definitions of sin/cos of one of the angles of the
0181       // right-angled triangle.
0182
0183       // But we start with a "reference phi" so that we can easily
0184       // decide which phi is the largest/smallest without having too
0185       // much of headache with the wrapping. The reference is in
0186       // principle defined by the averages of y&x phimin/phimax at
0187       // rmin, but the '0.5f*' factor is omitted here to reduce
0188       // computations.
0189       const auto phi_ref = std::atan2(y_rmin_phimin + y_rmin_phimax, x_rmin_phimin + x_rmin_phimax);
0190
0191       // for maximum phi we need the orthogonal vector to the left
0192       const auto tan_rmin_phimax = tangentVec(p_rmin_phimax, +1);
0193       const auto tan_rmax_phimax = tangentVec(p_rmax_phimax, +1);
0194       const auto phi_rmin_phimax = std::atan2(tan_rmin_phimax[1], tan_rmin_phimax[0]);
0195       const auto phi_rmax_phimax = std::atan2(tan_rmax_phimax[1], tan_rmax_phimax[0]);
0196
0197       auto phimax =
0198           std::abs(reco::deltaPhi(phi_rmin_phimax, phi_ref)) > std::abs(reco::deltaPhi(phi_rmax_phimax, phi_ref))
0199               ? phi_rmin_phimax
0200               : phi_rmax_phimax;
0201
0202       LogTrace("AreaSeededTrackingRegionBuilder")
0203           << "   rmin_phimax vec " << p_rmin_phimax[0] << "," << p_rmin_phimax[1] << " tangent " << tan_rmin_phimax[0]
0204           << "," << tan_rmin_phimax[1] << " phi " << phi_rmin_phimax << "\n"
0205           << "   rmax_phimax vec " << p_rmax_phimax[0] << "," << p_rmax_phimax[1] << " tangent " << tan_rmax_phimax[0]
0206           << "," << tan_rmax_phimax[1] << " phi " << phi_rmax_phimax << "\n"
0207           << "   phimax " << phimax;
0208
0209       // for minimum phi we need the orthogonal vector to the right
0210       const auto tan_rmin_phimin = tangentVec(p_rmin_phimin, -1);
0211       const auto tan_rmax_phimin = tangentVec(p_rmax_phimin, -1);
0212       const auto phi_rmin_phimin = std::atan2(tan_rmin_phimin[1], tan_rmin_phimin[0]);
0213       const auto phi_rmax_phimin = std::atan2(tan_rmax_phimin[1], tan_rmax_phimin[0]);
0214
0215       auto phimin =
0216           std::abs(reco::deltaPhi(phi_rmin_phimin, phi_ref)) > std::abs(reco::deltaPhi(phi_rmax_phimin, phi_ref))
0217               ? phi_rmin_phimin
0218               : phi_rmax_phimin;
0219
0220       LogTrace("AreaSeededTrackingRegionBuilder")
0221           << "   rmin_phimin vec " << p_rmin_phimin[0] << "," << p_rmin_phimin[1] << " tangent " << tan_rmin_phimin[0]
0222           << "," << tan_rmin_phimin[1] << " phi " << phi_rmin_phimin << "\n"
0223           << "   rmax_phimin vec " << p_rmax_phimin[0] << "," << p_rmax_phimin[1] << " tangent " << tan_rmax_phimin[0]
0224           << "," << tan_rmax_phimin[1] << " phi " << phi_rmax_phimin << "\n"
0225           << "   phimin " << phimin;
0226
0227       // wrapped around, need to decide which one to wrap
0228       if (phimax < phimin) {
0229         if (phimax < 0)
0230           phimax += 2 * M_PI;
0231         else
0232           phimin -= 2 * M_PI;
0233       }
0234
0235       LogTrace("AreaSeededTrackingRegionBuilder") << "  phi min,max " << phimin << "," << phimax;
0236
0237       minPhi = std::min(minPhi, phimin);
0238       maxPhi = std::max(maxPhi, phimax);
0239     }
0240
0241     LogTrace("AreaSeededTrackingRegionBuilder")
0242         << "  windows after this area  eta " << minEta << "," << maxEta << " phi " << minPhi << "," << maxPhi;
0243   }
0244
0245   const auto meanEta = (minEta + maxEta) / 2.f;
0246   const auto meanPhi = (minPhi + maxPhi) / 2.f;
0247   const auto dEta = maxEta - meanEta + m_conf->m_extraEta;
0248   const auto dPhi = maxPhi - meanPhi + m_conf->m_extraPhi;
0249
0250   auto useCandidates = false;
0251   if (candidates.first)
0252     useCandidates = true;
0253
0254   if (useCandidates) {
0255     // If we have objects used for seeding, loop over objects and find overlap with the found region. Return overlaps as tracking regions to use
0256     for (const auto& object : *candidates.first) {
0257       float dEta_Cand = candidates.second.first;
0258       float dPhi_Cand = candidates.second.second;
0259       float eta_Cand = object.eta();
0260       float phi_Cand = object.phi();
0261       float dEta_Cand_Point = std::abs(eta_Cand - meanEta);
0262       float dPhi_Cand_Point = std::abs(deltaPhi(phi_Cand, meanPhi));
0263
0264       if (dEta_Cand_Point > (dEta_Cand + dEta) || dPhi_Cand_Point > (dPhi_Cand + dPhi))
0265         continue;
0266
0267       float etaMin_RoI = std::max(eta_Cand - dEta_Cand, meanEta - dEta);
0268       float etaMax_RoI = std::min(eta_Cand + dEta_Cand, meanEta + dEta);
0269
0270       float phi_Cand_minus = normalizedPhi(phi_Cand - dPhi_Cand);
0271       float phi_Point_minus = normalizedPhi(meanPhi - dPhi);
0272       float phi_Cand_plus = normalizedPhi(phi_Cand + dPhi_Cand);
0273       float phi_Point_plus = normalizedPhi(meanPhi + dPhi);
0274
0275       float phiMin_RoI = deltaPhi(phi_Cand_minus, phi_Point_minus) > 0. ? phi_Cand_minus : phi_Point_minus;
0276       float phiMax_RoI = deltaPhi(phi_Cand_plus, phi_Point_plus) < 0. ? phi_Cand_plus : phi_Point_plus;
0277
0278       const auto meanEtaTemp = (etaMin_RoI + etaMax_RoI) / 2.f;
0279       auto meanPhiTemp = (phiMin_RoI + phiMax_RoI) / 2.f;
0280       if (phiMax_RoI < phiMin_RoI)
0281         meanPhiTemp -= M_PI;
0282       meanPhiTemp = normalizedPhi(meanPhiTemp);
0283
0284       const auto dPhiTemp = deltaPhi(phiMax_RoI, meanPhiTemp);
0285       const auto dEtaTemp = etaMax_RoI - meanEtaTemp;
0286
0287       const auto x = std::cos(meanPhiTemp);
0288       const auto y = std::sin(meanPhiTemp);
0289       const auto z = 1. / std::tan(2.f * std::atan(std::exp(-meanEtaTemp)));
0290
0291       LogTrace("AreaSeededTrackingRegionsBuilder")
0292           << "Direction x,y,z " << x << "," << y << "," << z << " eta,phi " << meanEtaTemp << "," << meanPhiTemp
0293           << " window eta " << (meanEtaTemp - dEtaTemp) << "," << (meanEtaTemp + dEtaTemp) << " phi "
0294           << (meanPhiTemp - dPhiTemp) << "," << (meanPhiTemp + dPhiTemp);
0295
0296       return std::make_unique<RectangularEtaPhiTrackingRegion>(GlobalVector(x, y, z),
0297                                                                origin.first,  // GlobalPoint
0298                                                                m_conf->m_ptMin,
0299                                                                m_conf->m_originRadius,
0300                                                                origin.second,
0301                                                                dEtaTemp,
0302                                                                dPhiTemp,
0303                                                                *m_field,
0304                                                                m_msmaker,
0305                                                                m_conf->m_precise,
0306                                                                m_conf->m_whereToUseMeasurementTracker,
0307                                                                m_measurementTracker,
0308                                                                m_conf->m_searchOpt);
0309     }
0310     // Have to retun nullptr here to ensure that we always return something
0311     return nullptr;
0312
0313   } else {
0314     const auto x = std::cos(meanPhi);
0315     const auto y = std::sin(meanPhi);
0316     const auto z = 1. / std::tan(2.f * std::atan(std::exp(-meanEta)));
0317
0318     LogTrace("AreaSeededTrackingRegionsBuilder")
0319         << "Direction x,y,z " << x << "," << y << "," << z << " eta,phi " << meanEta << "," << meanPhi << " window eta "
0320         << (meanEta - dEta) << "," << (meanEta + dEta) << " phi " << (meanPhi - dPhi) << "," << (meanPhi + dPhi);
0321
0322     return std::make_unique<RectangularEtaPhiTrackingRegion>(GlobalVector(x, y, z),
0323                                                              origin.first,  // GlobalPoint
0324                                                              m_conf->m_ptMin,
0325                                                              m_conf->m_originRadius,
0326                                                              origin.second,
0327                                                              dEta,
0328                                                              dPhi,
0329                                                              *m_field,
0330                                                              m_msmaker,
0331                                                              m_conf->m_precise,
0332                                                              m_conf->m_whereToUseMeasurementTracker,
0333                                                              m_measurementTracker,
0334                                                              m_conf->m_searchOpt);
0335   }
0336 }