File indexing completed on 2024-04-06 12:28:56
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 }
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
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
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
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
0125
0126 LogTrace("AreaSeededTrackingRegionsBuilder")
0127 << " area x,y points " << area.x_rmin_phimin() << "," << area.y_rmin_phimin() << " "
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
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
0144 {
0145
0146
0147
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
0152
0153
0154 const auto p_min = vecFromOrigPlusRadius(p_rmin);
0155 const auto p_max = vecFromOrigPlusRadius(p_rmax);
0156
0157
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
0170 {
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189 const auto phi_ref = std::atan2(y_rmin_phimin + y_rmin_phimax, x_rmin_phimin + x_rmin_phimax);
0190
0191
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
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
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
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,
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
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,
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 }