File indexing completed on 2023-03-17 11:22:51
0001 #include "PixelForwardLayerPhase1.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
0006 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0007
0008 #include "TrackingTools/DetLayers/interface/simple_stat.h"
0009 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing2Order.h"
0010 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0011 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0012
0013 #include "LayerCrossingSide.h"
0014 #include "DetGroupMerger.h"
0015 #include "CompatibleDetToGroupAdder.h"
0016
0017 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0018
0019 #include <algorithm>
0020 #include <climits>
0021
0022 using namespace std;
0023
0024 typedef GeometricSearchDet::DetWithState DetWithState;
0025
0026 PixelForwardLayerPhase1::PixelForwardLayerPhase1(vector<const Phase1PixelBlade*>& blades)
0027 : ForwardDetLayer(true), theComps(blades.begin(), blades.end()) {
0028
0029
0030
0031
0032
0033
0034
0035 float theRmin = (*(theComps.begin()))->surface().position().perp();
0036 float theRmax = theRmin;
0037 for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0038 theRmin = std::min(theRmin, (*it)->surface().position().perp());
0039 theRmax = std::max(theRmax, (*it)->surface().position().perp());
0040 }
0041
0042
0043 float split_inner_outer_radius = 10.;
0044 _num_innerpanels = 0;
0045 for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0046 if ((**it).surface().position().perp() <= split_inner_outer_radius)
0047 ++_num_innerpanels;
0048 }
0049 _num_outerpanels = theComps.size() - _num_innerpanels;
0050 edm::LogInfo("TkDetLayers") << " Rmin, Rmax, R_average = " << theRmin << ", " << theRmax << ", "
0051 << split_inner_outer_radius << std::endl
0052 << " num inner, outer disks = " << _num_innerpanels << ", " << _num_outerpanels
0053 << std::endl;
0054
0055 for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0056 theBasicComps.insert(theBasicComps.end(), (**it).basicComponents().begin(), (**it).basicComponents().end());
0057 }
0058
0059
0060
0061 setSurface(computeSurface());
0062
0063 theBinFinder_inner = BinFinderType(theComps.front()->surface().position().phi(), _num_innerpanels);
0064 theBinFinder_outer = BinFinderType(theComps[_num_innerpanels]->surface().position().phi(), _num_outerpanels);
0065
0066
0067 LogDebug("TkDetLayers") << "DEBUG INFO for PixelForwardLayerPhase1"
0068 << "\n"
0069 << "Num of inner and outer panels = " << _num_innerpanels << ", " << _num_outerpanels << "\n"
0070 << "Based on phi separation for inner: " << theComps.front()->surface().position().phi()
0071 << " and on phi separation for outer: "
0072 << theComps[_num_innerpanels]->surface().position().phi()
0073 << "PixelForwardLayerPhase1.surfcace.phi(): " << std::endl
0074 << this->surface().position().phi() << "\n"
0075 << "PixelForwardLayerPhase1.surfcace.z(): " << this->surface().position().z() << "\n"
0076 << "PixelForwardLayerPhase1.surfcace.innerR(): " << this->specificSurface().innerRadius()
0077 << "\n"
0078 << "PixelForwardLayerPhase1.surfcace.outerR(): " << this->specificSurface().outerRadius();
0079
0080 for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0081 LogDebug("TkDetLayers") << "blades phi,z,r: " << (*it)->surface().position().phi() << " , "
0082 << (*it)->surface().position().z() << " , " << (*it)->surface().position().perp();
0083
0084
0085
0086
0087 }
0088
0089 }
0090
0091 PixelForwardLayerPhase1::~PixelForwardLayerPhase1() {
0092 vector<const GeometricSearchDet*>::const_iterator i;
0093 for (i = theComps.begin(); i != theComps.end(); i++) {
0094 delete *i;
0095 }
0096 }
0097
0098 void PixelForwardLayerPhase1::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0099 const Propagator& prop,
0100 const MeasurementEstimator& est,
0101 std::vector<DetGroup>& result) const {
0102 std::vector<DetGroup> closestResult_inner;
0103 std::vector<DetGroup> closestResult_outer;
0104 std::vector<DetGroup> nextResult_inner;
0105 std::vector<DetGroup> nextResult_outer;
0106 std::vector<DetGroup> result_inner;
0107 std::vector<DetGroup> result_outer;
0108 int frontindex_inner = 0;
0109 int frontindex_outer = 0;
0110 SubTurbineCrossings crossings_inner;
0111 SubTurbineCrossings crossings_outer;
0112
0113 if (_num_innerpanels > 0)
0114 crossings_inner = computeCrossings(tsos, prop.propagationDirection(), true);
0115 crossings_outer = computeCrossings(tsos, prop.propagationDirection(), false);
0116
0117
0118
0119
0120 if (!crossings_outer.isValid) {
0121 edm::LogInfo("TkDetLayers")
0122 << "outer computeCrossings returns invalid in PixelForwardLayerPhase1::groupedCompatibleDets:";
0123 return;
0124 }
0125
0126 typedef CompatibleDetToGroupAdder Adder;
0127 if (crossings_inner.isValid) {
0128 Adder::add(
0129 *theComps[theBinFinder_inner.binIndex(crossings_inner.closestIndex)], tsos, prop, est, closestResult_inner);
0130
0131 if (closestResult_inner.empty()) {
0132 Adder::add(*theComps[theBinFinder_inner.binIndex(crossings_inner.nextIndex)], tsos, prop, est, result_inner);
0133 frontindex_inner = crossings_inner.nextIndex;
0134 } else {
0135 if (Adder::add(
0136 *theComps[theBinFinder_inner.binIndex(crossings_inner.nextIndex)], tsos, prop, est, nextResult_inner)) {
0137 int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0138 int theHelicity = computeHelicity(theComps[theBinFinder_inner.binIndex(crossings_inner.closestIndex)],
0139 theComps[theBinFinder_inner.binIndex(crossings_inner.nextIndex)]);
0140 std::vector<DetGroup> tmp99 = closestResult_inner;
0141 DetGroupMerger::orderAndMergeTwoLevels(
0142 std::move(tmp99), std::move(nextResult_inner), result_inner, theHelicity, crossingSide);
0143 if (theHelicity == crossingSide)
0144 frontindex_inner = crossings_inner.closestIndex;
0145 else
0146 frontindex_inner = crossings_inner.nextIndex;
0147 } else {
0148 result_inner.swap(closestResult_inner);
0149 frontindex_inner = crossings_inner.closestIndex;
0150 }
0151 }
0152 }
0153 if (!closestResult_inner.empty()) {
0154 DetGroupElement closestGel(closestResult_inner.front().front());
0155 float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0156 searchNeighbors(tsos, prop, est, crossings_inner, window, result_inner, true);
0157 }
0158
0159
0160
0161
0162
0163
0164 Adder::add(*theComps[(theBinFinder_outer.binIndex(crossings_outer.closestIndex)) + _num_innerpanels],
0165 tsos,
0166 prop,
0167 est,
0168 closestResult_outer);
0169
0170 if (closestResult_outer.empty()) {
0171 Adder::add(*theComps[theBinFinder_outer.binIndex(crossings_outer.nextIndex) + _num_innerpanels],
0172 tsos,
0173 prop,
0174 est,
0175 result_outer);
0176 frontindex_outer = crossings_outer.nextIndex;
0177 } else {
0178 if (Adder::add(*theComps[theBinFinder_outer.binIndex(crossings_outer.nextIndex) + _num_innerpanels],
0179 tsos,
0180 prop,
0181 est,
0182 nextResult_outer)) {
0183 int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0184 int theHelicity =
0185 computeHelicity(theComps[theBinFinder_outer.binIndex(crossings_outer.closestIndex) + _num_innerpanels],
0186 theComps[theBinFinder_outer.binIndex(crossings_outer.nextIndex) + _num_innerpanels]);
0187 std::vector<DetGroup> tmp99 = closestResult_outer;
0188 DetGroupMerger::orderAndMergeTwoLevels(
0189 std::move(tmp99), std::move(nextResult_outer), result_outer, theHelicity, crossingSide);
0190 if (theHelicity == crossingSide)
0191 frontindex_outer = crossings_outer.closestIndex;
0192 else
0193 frontindex_outer = crossings_outer.nextIndex;
0194 } else {
0195 result_outer.swap(closestResult_outer);
0196 frontindex_outer = crossings_outer.closestIndex;
0197 }
0198 }
0199 if (!closestResult_outer.empty()) {
0200 DetGroupElement closestGel(closestResult_outer.front().front());
0201 float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0202 searchNeighbors(tsos, prop, est, crossings_outer, window, result_outer, false);
0203 }
0204
0205 if (result_inner.empty() && result_outer.empty())
0206 return;
0207 if (result_inner.empty())
0208 result.swap(result_outer);
0209 else if (result_outer.empty())
0210 result.swap(result_inner);
0211 else {
0212 int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0213 int theHelicity = computeHelicity(theComps[theBinFinder_inner.binIndex(frontindex_inner)],
0214 theComps[theBinFinder_outer.binIndex(frontindex_outer) + _num_innerpanels]);
0215 DetGroupMerger::orderAndMergeTwoLevels(
0216 std::move(result_inner), std::move(result_outer), result, theHelicity, crossingSide);
0217 }
0218 }
0219
0220 void PixelForwardLayerPhase1::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0221 const Propagator& prop,
0222 const MeasurementEstimator& est,
0223 const SubTurbineCrossings& crossings,
0224 float window,
0225 vector<DetGroup>& result,
0226 bool innerDisk) const {
0227 typedef CompatibleDetToGroupAdder Adder;
0228 int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0229 typedef DetGroupMerger Merger;
0230
0231 int negStart = min(crossings.closestIndex, crossings.nextIndex) - 1;
0232 int posStart = max(crossings.closestIndex, crossings.nextIndex) + 1;
0233
0234 int quarter = (innerDisk ? _num_innerpanels / 4 : _num_outerpanels / 4);
0235
0236 for (int idet = negStart; idet >= negStart - quarter + 1; idet--) {
0237 std::vector<DetGroup> tmp1;
0238 std::vector<DetGroup> newResult;
0239 if (innerDisk) {
0240 const GeometricSearchDet* neighbor = theComps[theBinFinder_inner.binIndex(idet)];
0241
0242
0243 if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0244 break;
0245 int theHelicity =
0246 computeHelicity(theComps[theBinFinder_inner.binIndex(idet)], theComps[theBinFinder_inner.binIndex(idet + 1)]);
0247 std::vector<DetGroup> tmp2;
0248 tmp2.swap(result);
0249 Merger::orderAndMergeTwoLevels(std::move(tmp1), std::move(tmp2), newResult, theHelicity, crossingSide);
0250 } else {
0251 const GeometricSearchDet* neighbor = theComps[(theBinFinder_outer.binIndex(idet)) + _num_innerpanels];
0252
0253
0254 if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0255 break;
0256 int theHelicity = computeHelicity(theComps[(theBinFinder_outer.binIndex(idet)) + _num_innerpanels],
0257 theComps[(theBinFinder_outer.binIndex(idet + 1)) + _num_innerpanels]);
0258 std::vector<DetGroup> tmp2;
0259 tmp2.swap(result);
0260 Merger::orderAndMergeTwoLevels(std::move(tmp1), std::move(tmp2), newResult, theHelicity, crossingSide);
0261 }
0262 result.swap(newResult);
0263 }
0264 for (int idet = posStart; idet < posStart + quarter - 1; idet++) {
0265 std::vector<DetGroup> tmp1;
0266 std::vector<DetGroup> newResult;
0267 if (innerDisk) {
0268 const GeometricSearchDet* neighbor = theComps[theBinFinder_inner.binIndex(idet)];
0269
0270
0271 if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0272 break;
0273 int theHelicity =
0274 computeHelicity(theComps[theBinFinder_inner.binIndex(idet - 1)], theComps[theBinFinder_inner.binIndex(idet)]);
0275 std::vector<DetGroup> tmp2;
0276 tmp2.swap(result);
0277 Merger::orderAndMergeTwoLevels(std::move(tmp2), std::move(tmp1), newResult, theHelicity, crossingSide);
0278 } else {
0279 const GeometricSearchDet* neighbor = theComps[(theBinFinder_outer.binIndex(idet)) + _num_innerpanels];
0280
0281
0282 if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0283 break;
0284 int theHelicity = computeHelicity(theComps[(theBinFinder_outer.binIndex(idet - 1)) + _num_innerpanels],
0285 theComps[(theBinFinder_outer.binIndex(idet)) + _num_innerpanels]);
0286 std::vector<DetGroup> tmp2;
0287 tmp2.swap(result);
0288 Merger::orderAndMergeTwoLevels(std::move(tmp2), std::move(tmp1), newResult, theHelicity, crossingSide);
0289 }
0290 result.swap(newResult);
0291 }
0292 }
0293
0294 int PixelForwardLayerPhase1::computeHelicity(const GeometricSearchDet* firstBlade,
0295 const GeometricSearchDet* secondBlade) {
0296 return std::abs(firstBlade->position().z()) < std::abs(secondBlade->position().z()) ? 0 : 1;
0297 }
0298
0299 PixelForwardLayerPhase1::SubTurbineCrossings PixelForwardLayerPhase1::computeCrossings(
0300 const TrajectoryStateOnSurface& startingState, PropagationDirection propDir, bool innerDisk) const {
0301 typedef MeasurementEstimator::Local2DVector Local2DVector;
0302
0303 HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0304 HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0305 auto rho = startingState.transverseCurvature();
0306
0307 HelixArbitraryPlaneCrossing turbineCrossing(startPos, startDir, rho, propDir);
0308
0309 pair<bool, double> thePath = turbineCrossing.pathLength(specificSurface());
0310
0311 if (!thePath.first) {
0312 return SubTurbineCrossings();
0313 }
0314
0315 HelixPlaneCrossing::PositionType turbinePoint(turbineCrossing.position(thePath.second));
0316 HelixPlaneCrossing::DirectionType turbineDir(turbineCrossing.direction(thePath.second));
0317 int closestIndex = 0;
0318 int nextIndex = 0;
0319 if (innerDisk)
0320 closestIndex = theBinFinder_inner.binIndex(turbinePoint.barePhi());
0321 else
0322 closestIndex = theBinFinder_outer.binIndex(turbinePoint.barePhi());
0323
0324 HelixArbitraryPlaneCrossing2Order theBladeCrossing(turbinePoint, turbineDir, rho);
0325
0326 float closestDist = 0;
0327 if (innerDisk) {
0328 const BoundPlane& closestPlane(static_cast<const BoundPlane&>(theComps[closestIndex]->surface()));
0329
0330 pair<bool, double> theClosestBladePath = theBladeCrossing.pathLength(closestPlane);
0331 LocalPoint closestPos = closestPlane.toLocal(GlobalPoint(theBladeCrossing.position(theClosestBladePath.second)));
0332 closestDist = closestPos.x();
0333
0334 nextIndex = Geom::phiLess(closestPlane.phi(), turbinePoint.barePhi()) ? closestIndex + 1 : closestIndex - 1;
0335 } else {
0336 const BoundPlane& closestPlane(
0337 static_cast<const BoundPlane&>(theComps[closestIndex + _num_innerpanels]->surface()));
0338
0339 pair<bool, double> theClosestBladePath = theBladeCrossing.pathLength(closestPlane);
0340 LocalPoint closestPos = closestPlane.toLocal(GlobalPoint(theBladeCrossing.position(theClosestBladePath.second)));
0341 closestDist = closestPos.x();
0342
0343 nextIndex = Geom::phiLess(closestPlane.phi(), turbinePoint.barePhi()) ? closestIndex + 1 : closestIndex - 1;
0344 }
0345
0346 float nextDist = 0;
0347 if (innerDisk) {
0348 const BoundPlane& nextPlane(
0349 static_cast<const BoundPlane&>(theComps[theBinFinder_inner.binIndex(nextIndex)]->surface()));
0350 pair<bool, double> theNextBladePath = theBladeCrossing.pathLength(nextPlane);
0351 LocalPoint nextPos = nextPlane.toLocal(GlobalPoint(theBladeCrossing.position(theNextBladePath.second)));
0352 nextDist = nextPos.x();
0353 } else {
0354 const BoundPlane& nextPlane(
0355 static_cast<const BoundPlane&>(theComps[theBinFinder_outer.binIndex(nextIndex) + _num_innerpanels]->surface()));
0356 pair<bool, double> theNextBladePath = theBladeCrossing.pathLength(nextPlane);
0357 LocalPoint nextPos = nextPlane.toLocal(GlobalPoint(theBladeCrossing.position(theNextBladePath.second)));
0358 nextDist = nextPos.x();
0359 }
0360 if (std::abs(closestDist) < std::abs(nextDist)) {
0361 return SubTurbineCrossings(closestIndex, nextIndex, nextDist);
0362 } else {
0363 return SubTurbineCrossings(nextIndex, closestIndex, closestDist);
0364 }
0365 }
0366
0367 float PixelForwardLayerPhase1::computeWindowSize(const GeomDet* det,
0368 const TrajectoryStateOnSurface& tsos,
0369 const MeasurementEstimator& est) {
0370 return est.maximalLocalDisplacement(tsos, det->surface()).x();
0371 }