TkDetLayers/src/PixelForwardLayer.cc

0001 #include "PixelForwardLayer.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 "DataFormats/GeometryVector/interface/VectorUtil.h"
0010 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing2Order.h"
0011 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0012 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0013
0014 #include "LayerCrossingSide.h"
0015 #include "DetGroupMerger.h"
0016 #include "CompatibleDetToGroupAdder.h"
0017
0018 using namespace std;
0019
0020 typedef GeometricSearchDet::DetWithState DetWithState;
0021
0022 PixelForwardLayer::PixelForwardLayer(vector<const PixelBlade*>& blades)
0023     : ForwardDetLayer(true), theComps(blades.begin(), blades.end()) {
0024   for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0025     theBasicComps.insert(theBasicComps.end(), (**it).basicComponents().begin(), (**it).basicComponents().end());
0026   }
0027
0028   //They should be already phi-ordered. TO BE CHECKED!!
0029   //sort( theBlades.begin(), theBlades.end(), PhiLess());
0030   setSurface(computeSurface());
0031
0032   //Is a "periodic" binFinderInPhi enough?. TO BE CHECKED!!
0033   theBinFinder = BinFinderType(theComps.front()->surface().position().phi(), theComps.size());
0034
0035   //--------- DEBUG INFO --------------
0036   LogDebug("TkDetLayers") << "DEBUG INFO for PixelForwardLayer"
0037                           << "\n"
0038                           << "PixelForwardLayer.surfcace.phi(): " << this->surface().position().phi() << "\n"
0039                           << "PixelForwardLayer.surfcace.z(): " << this->surface().position().z() << "\n"
0040                           << "PixelForwardLayer.surfcace.innerR(): " << this->specificSurface().innerRadius() << "\n"
0041                           << "PixelForwardLayer.surfcace.outerR(): " << this->specificSurface().outerRadius();
0042
0043   for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0044     LogDebug("TkDetLayers") << "blades phi,z,r: " << (*it)->surface().position().phi() << " , "
0045                             << (*it)->surface().position().z() << " , " << (*it)->surface().position().perp();
0046   }
0047   //-----------------------------------
0048 }
0049
0050 PixelForwardLayer::~PixelForwardLayer() {
0051   vector<const GeometricSearchDet*>::const_iterator i;
0052   for (i = theComps.begin(); i != theComps.end(); i++) {
0053     delete *i;
0054   }
0055 }
0056
0057 void PixelForwardLayer::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0058                                                const Propagator& prop,
0059                                                const MeasurementEstimator& est,
0060                                                std::vector<DetGroup>& result) const {
0061   std::vector<DetGroup> closestResult;
0062   SubTurbineCrossings crossings;
0063
0064   crossings = computeCrossings(tsos, prop.propagationDirection());
0065   if (!crossings.isValid) {
0066     //edm::LogInfo("TkDetLayers") << "computeCrossings returns invalid in PixelForwardLayer::groupedCompatibleDets:";
0067     return;
0068   }
0069
0070   typedef CompatibleDetToGroupAdder Adder;
0071   Adder::add(*theComps[theBinFinder.binIndex(crossings.closestIndex)], tsos, prop, est, closestResult);
0072
0073   if (closestResult.empty()) {
0074     Adder::add(*theComps[theBinFinder.binIndex(crossings.nextIndex)], tsos, prop, est, result);
0075     return;
0076   }
0077
0078   DetGroupElement closestGel(closestResult.front().front());
0079   float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0080
0081   //float detWidth = closestGel.det()->surface().bounds().width();
0082   //if (crossings.nextDistance < detWidth + window) {
0083   vector<DetGroup> nextResult;
0084   if (Adder::add(*theComps[theBinFinder.binIndex(crossings.nextIndex)], tsos, prop, est, nextResult)) {
0085     int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0086     int theHelicity = computeHelicity(theComps[theBinFinder.binIndex(crossings.closestIndex)],
0087                                       theComps[theBinFinder.binIndex(crossings.nextIndex)]);
0088     DetGroupMerger::orderAndMergeTwoLevels(
0089         std::move(closestResult), std::move(nextResult), result, theHelicity, crossingSide);
0090   } else {
0091     result.swap(closestResult);
0092   }
0093
0094   /*
0095   }
0096   else {
0097     result.swap(closestResult);
0098   }
0099   */
0100
0101   // --- THIS lines may speed up the reconstruction. But it reduces slightly the efficiency.
0102   // only loop over neighbors (other than closest and next) if window is BIG
0103   //if (window > 0.5*detWidth) {
0104   searchNeighbors(tsos, prop, est, crossings, window, result);
0105   //}
0106 }
0107
0108 void PixelForwardLayer::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0109                                         const Propagator& prop,
0110                                         const MeasurementEstimator& est,
0111                                         const SubTurbineCrossings& crossings,
0112                                         float window,
0113                                         vector<DetGroup>& result) const {
0114   typedef CompatibleDetToGroupAdder Adder;
0115   int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0116   typedef DetGroupMerger Merger;
0117
0118   int negStart = min(crossings.closestIndex, crossings.nextIndex) - 1;
0119   int posStart = max(crossings.closestIndex, crossings.nextIndex) + 1;
0120
0121   int quarter = theComps.size() / 4;
0122
0123   for (int idet = negStart; idet >= negStart - quarter + 1; idet--) {
0124     std::vector<DetGroup> tmp1;
0125     const GeometricSearchDet* neighbor = theComps[theBinFinder.binIndex(idet)];
0126     // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0127     // maybe also add shallow crossing angle test here???
0128     if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0129       break;
0130     int theHelicity = computeHelicity(theComps[theBinFinder.binIndex(idet)], theComps[theBinFinder.binIndex(idet + 1)]);
0131     std::vector<DetGroup> tmp2;
0132     tmp2.swap(result);
0133     std::vector<DetGroup> newResult;
0134     Merger::orderAndMergeTwoLevels(std::move(tmp1), std::move(tmp2), newResult, theHelicity, crossingSide);
0135     result.swap(newResult);
0136   }
0137   for (int idet = posStart; idet < posStart + quarter - 1; idet++) {
0138     vector<DetGroup> tmp1;
0139     const GeometricSearchDet* neighbor = theComps[theBinFinder.binIndex(idet)];
0140     // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0141     // maybe also add shallow crossing angle test here???
0142     if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0143       break;
0144     int theHelicity = computeHelicity(theComps[theBinFinder.binIndex(idet - 1)], theComps[theBinFinder.binIndex(idet)]);
0145     std::vector<DetGroup> tmp2;
0146     tmp2.swap(result);
0147     std::vector<DetGroup> newResult;
0148     Merger::orderAndMergeTwoLevels(std::move(tmp2), std::move(tmp1), newResult, theHelicity, crossingSide);
0149     result.swap(newResult);
0150   }
0151 }
0152
0153 int PixelForwardLayer::computeHelicity(const GeometricSearchDet* firstBlade, const GeometricSearchDet* secondBlade) {
0154   return std::abs(firstBlade->position().z()) < std::abs(secondBlade->position().z()) ? 0 : 1;
0155 }
0156
0157 PixelForwardLayer::SubTurbineCrossings PixelForwardLayer::computeCrossings(
0158     const TrajectoryStateOnSurface& startingState, PropagationDirection propDir) const {
0159   typedef MeasurementEstimator::Local2DVector Local2DVector;
0160
0161   HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0162   HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0163
0164   auto rho = startingState.transverseCurvature();
0165
0166   HelixArbitraryPlaneCrossing turbineCrossing(startPos, startDir, rho, propDir);
0167
0168   pair<bool, double> thePath = turbineCrossing.pathLength(specificSurface());
0169
0170   if (!thePath.first) {
0171     //edm::LogInfo("TkDetLayers") << "ERROR in PixelForwardLayer: disk not crossed by track" ;
0172     return SubTurbineCrossings();
0173   }
0174
0175   HelixPlaneCrossing::PositionType turbinePoint(turbineCrossing.position(thePath.second));
0176   HelixPlaneCrossing::DirectionType turbineDir(turbineCrossing.direction(thePath.second));
0177
0178   int closestIndex = theBinFinder.binIndex(turbinePoint.barePhi());
0179
0180   const Plane& closestPlane(static_cast<const Plane&>(theComps[closestIndex]->surface()));
0181
0182   HelixArbitraryPlaneCrossing2Order theBladeCrossing(turbinePoint, turbineDir, rho);
0183
0184   pair<bool, double> theClosestBladePath = theBladeCrossing.pathLength(closestPlane);
0185   LocalPoint closestPos = closestPlane.toLocal(GlobalPoint(theBladeCrossing.position(theClosestBladePath.second)));
0186
0187   auto closestDist = closestPos.x();  // use fact that local X perp to global Y
0188
0189   //int next = turbinePoint.phi() - closestPlane.position().phi() > 0 ? closest+1 : closest-1;
0190
0191   int nextIndex = Geom::phiLess(closestPlane.phi(), turbinePoint.barePhi()) ? closestIndex + 1 : closestIndex - 1;
0192
0193   const Plane& nextPlane(static_cast<const Plane&>(theComps[theBinFinder.binIndex(nextIndex)]->surface()));
0194
0195   pair<bool, double> theNextBladePath = theBladeCrossing.pathLength(nextPlane);
0196   LocalPoint nextPos = nextPlane.toLocal(GlobalPoint(theBladeCrossing.position(theNextBladePath.second)));
0197
0198   auto nextDist = nextPos.x();
0199
0200   if (std::abs(closestDist) < std::abs(nextDist)) {
0201     return SubTurbineCrossings(closestIndex, nextIndex, nextDist);
0202   } else {
0203     return SubTurbineCrossings(nextIndex, closestIndex, closestDist);
0204   }
0205 }
0206
0207 float PixelForwardLayer::computeWindowSize(const GeomDet* det,
0208                                            const TrajectoryStateOnSurface& tsos,
0209                                            const MeasurementEstimator& est) {
0210   return est.maximalLocalDisplacement(tsos, det->surface()).x();
0211 }