TkDetLayers/src/Phase1PixelBlade.cc

0001 #include "Phase1PixelBlade.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "BladeShapeBuilderFromDet.h"
0006 #include "LayerCrossingSide.h"
0007 #include "DetGroupMerger.h"
0008 #include "CompatibleDetToGroupAdder.h"
0009 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
0010
0011 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0012 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0013 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0014
0015 using namespace std;
0016
0017 typedef GeometricSearchDet::DetWithState DetWithState;
0018
0019 Phase1PixelBlade::~Phase1PixelBlade() {}
0020
0021 Phase1PixelBlade::Phase1PixelBlade(vector<const GeomDet*>& frontDets, vector<const GeomDet*>& backDets)
0022     : GeometricSearchDet(true),
0023       theFrontDets(frontDets),
0024       theBackDets(backDets),
0025       front_radius_range_(std::make_pair(0, 0)),
0026       back_radius_range_(std::make_pair(0, 0)) {
0027   theDets.assign(theFrontDets.begin(), theFrontDets.end());
0028   theDets.insert(theDets.end(), theBackDets.begin(), theBackDets.end());
0029
0030   theDiskSector = BladeShapeBuilderFromDet::build(theDets);
0031   theFrontDiskSector = BladeShapeBuilderFromDet::build(theFrontDets);
0032   theBackDiskSector = BladeShapeBuilderFromDet::build(theBackDets);
0033
0034   //--------- DEBUG INFO --------------
0035   LogDebug("TkDetLayers") << "DEBUG INFO for Phase1PixelBlade";
0036   LogDebug("TkDetLayers") << "Blade z, perp, innerRadius, outerR[disk, front, back]: " << this->position().z() << " , "
0037                           << this->position().perp() << " , (" << theDiskSector->innerRadius() << " , "
0038                           << theDiskSector->outerRadius() << "), (" << theFrontDiskSector->innerRadius() << " , "
0039                           << theFrontDiskSector->outerRadius() << "), (" << theBackDiskSector->innerRadius() << " , "
0040                           << theBackDiskSector->outerRadius() << ")" << std::endl;
0041
0042   for (vector<const GeomDet*>::const_iterator it = theFrontDets.begin(); it != theFrontDets.end(); it++) {
0043     LogDebug("TkDetLayers") << "frontDet phi,z,r: " << (*it)->position().phi() << " , " << (*it)->position().z()
0044                             << " , " << (*it)->position().perp() << " , "
0045                             << " rmin: " << (*it)->surface().rSpan().first
0046                             << " rmax: " << (*it)->surface().rSpan().second << std::endl;
0047   }
0048
0049   for (vector<const GeomDet*>::const_iterator it = theBackDets.begin(); it != theBackDets.end(); it++) {
0050     LogDebug("TkDetLayers") << "backDet phi,z,r: " << (*it)->position().phi() << " , " << (*it)->position().z() << " , "
0051                             << (*it)->position().perp() << " , "
0052                             << " rmin: " << (*it)->surface().rSpan().first
0053                             << " rmax: " << (*it)->surface().rSpan().second << std::endl;
0054   }
0055 }
0056
0057 const vector<const GeometricSearchDet*>& Phase1PixelBlade::components() const {
0058   throw DetLayerException("TOBRod doesn't have GeometricSearchDet components");
0059 }
0060
0061 pair<bool, TrajectoryStateOnSurface> Phase1PixelBlade::compatible(const TrajectoryStateOnSurface& ts,
0062                                                                   const Propagator&,
0063                                                                   const MeasurementEstimator&) const {
0064   edm::LogError("TkDetLayers") << "temporary dummy implementation of Phase1PixelBlade::compatible()!!";
0065   return pair<bool, TrajectoryStateOnSurface>();
0066 }
0067
0068 void Phase1PixelBlade::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0069                                               const Propagator& prop,
0070                                               const MeasurementEstimator& est,
0071                                               std::vector<DetGroup>& result) const {
0072   SubLayerCrossings crossings;
0073
0074   crossings = computeCrossings(tsos, prop.propagationDirection());
0075   if (!crossings.isValid())
0076     return;
0077
0078   vector<DetGroup> closestResult;
0079   addClosest(tsos, prop, est, crossings.closest(), closestResult);
0080
0081   if (closestResult.empty()) {
0082     vector<DetGroup> nextResult;
0083     addClosest(tsos, prop, est, crossings.other(), nextResult);
0084     if (nextResult.empty())
0085       return;
0086
0087     //DetGroupElement nextGel( nextResult.front().front());
0088     //int crossingSide = LayerCrossingSide().endcapSide( nextGel.trajectoryState(), prop);
0089
0090     DetGroupMerger::orderAndMergeTwoLevels(std::move(closestResult),
0091                                            std::move(nextResult),
0092                                            result,
0093                                            0,
0094                                            0);  //fixme gc patched for SLHC - already correctly sorted for SLHC
0095     //crossings.closestIndex(), crossingSide);
0096   } else {
0097     DetGroupElement closestGel(closestResult.front().front());
0098     float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0099
0100     searchNeighbors(tsos, prop, est, crossings.closest(), window, closestResult, false);
0101
0102     vector<DetGroup> nextResult;
0103     searchNeighbors(tsos, prop, est, crossings.other(), window, nextResult, true);
0104
0105     //int crossingSide = LayerCrossingSide().endcapSide( closestGel.trajectoryState(), prop);
0106     DetGroupMerger::orderAndMergeTwoLevels(std::move(closestResult),
0107                                            std::move(nextResult),
0108                                            result,
0109                                            0,
0110                                            0);  //fixme gc patched for SLHC - already correctly sorted for SLHC
0111     //crossings.closestIndex(), crossingSide);
0112   }
0113 }
0114
0115 SubLayerCrossings Phase1PixelBlade::computeCrossings(const TrajectoryStateOnSurface& startingState,
0116                                                      PropagationDirection propDir) const {
0117   HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0118   HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0119   double rho(startingState.transverseCurvature());
0120   HelixArbitraryPlaneCrossing crossing(startPos, startDir, rho, propDir);
0121
0122   pair<bool, double> innerPath = crossing.pathLength(*theFrontDiskSector);
0123   if (!innerPath.first)
0124     return SubLayerCrossings();
0125
0126   GlobalPoint gInnerPoint(crossing.position(innerPath.second));
0127   //Code for use of binfinder
0128   //int innerIndex = theInnerBinFinder.binIndex(gInnerPoint.perp());
0129   //float innerDist = std::abs( theInnerBinFinder.binPosition(innerIndex) - gInnerPoint.z());
0130
0131   //  int innerIndex = findBin(gInnerPoint.perp(),0);
0132   int innerIndex = findBin2(gInnerPoint, 0);
0133
0134   //fixme gc patched for SLHC - force order here to be in z
0135   //float innerDist = std::abs( findPosition(innerIndex,0).perp() - gInnerPoint.perp());
0136   float innerDist = (startingState.globalPosition() - gInnerPoint).mag();
0137   SubLayerCrossing innerSLC(0, innerIndex, gInnerPoint);
0138
0139   pair<bool, double> outerPath = crossing.pathLength(*theBackDiskSector);
0140   if (!outerPath.first)
0141     return SubLayerCrossings();
0142
0143   GlobalPoint gOuterPoint(crossing.position(outerPath.second));
0144   //Code for use of binfinder
0145   //int outerIndex = theOuterBinFinder.binIndex(gOuterPoint.perp());
0146   //float outerDist = std::abs( theOuterBinFinder.binPosition(outerIndex) - gOuterPoint.perp());
0147   //  int outerIndex  = findBin(gOuterPoint.perp(),1);
0148   int outerIndex = findBin2(gOuterPoint, 1);
0149
0150   //fixme gc patched for SLHC - force order here to be in z
0151   //float outerDist = std::abs( findPosition(outerIndex,1).perp() - gOuterPoint.perp());
0152   float outerDist = (startingState.globalPosition() - gOuterPoint).mag();
0153   SubLayerCrossing outerSLC(1, outerIndex, gOuterPoint);
0154
0155   if (innerDist < outerDist) {
0156     return SubLayerCrossings(innerSLC, outerSLC, 0);
0157   } else {
0158     return SubLayerCrossings(outerSLC, innerSLC, 1);
0159   }
0160 }
0161
0162 bool Phase1PixelBlade::addClosest(const TrajectoryStateOnSurface& tsos,
0163                                   const Propagator& prop,
0164                                   const MeasurementEstimator& est,
0165                                   const SubLayerCrossing& crossing,
0166                                   vector<DetGroup>& result) const {
0167   const vector<const GeomDet*>& sBlade(subBlade(crossing.subLayerIndex()));
0168
0169   return CompatibleDetToGroupAdder().add(*sBlade[crossing.closestDetIndex()], tsos, prop, est, result);
0170 }
0171
0172 float Phase1PixelBlade::computeWindowSize(const GeomDet* det,
0173                                           const TrajectoryStateOnSurface& tsos,
0174                                           const MeasurementEstimator& est) const {
0175   return est.maximalLocalDisplacement(tsos, det->surface()).x();
0176 }
0177
0178 void Phase1PixelBlade::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0179                                        const Propagator& prop,
0180                                        const MeasurementEstimator& est,
0181                                        const SubLayerCrossing& crossing,
0182                                        float window,
0183                                        vector<DetGroup>& result,
0184                                        bool checkClosest) const {
0185   const GlobalPoint& gCrossingPos = crossing.position();
0186   const vector<const GeomDet*>& sBlade(subBlade(crossing.subLayerIndex()));
0187
0188   int closestIndex = crossing.closestDetIndex();
0189   int negStartIndex = closestIndex - 1;
0190   int posStartIndex = closestIndex + 1;
0191
0192   if (checkClosest) {  // must decide if the closest is on the neg or pos side
0193     if (gCrossingPos.perp() < sBlade[closestIndex]->surface().position().perp()) {
0194       posStartIndex = closestIndex;
0195     } else {
0196       negStartIndex = closestIndex;
0197     }
0198   }
0199
0200   typedef CompatibleDetToGroupAdder Adder;
0201   for (int idet = negStartIndex; idet >= 0; idet--) {
0202     if (!overlap(gCrossingPos, *sBlade[idet], window))
0203       break;
0204     if (!Adder::add(*sBlade[idet], tsos, prop, est, result))
0205       break;
0206   }
0207   for (int idet = posStartIndex; idet < static_cast<int>(sBlade.size()); idet++) {
0208     if (!overlap(gCrossingPos, *sBlade[idet], window))
0209       break;
0210     if (!Adder::add(*sBlade[idet], tsos, prop, est, result))
0211       break;
0212   }
0213 }
0214
0215 bool Phase1PixelBlade::overlap(const GlobalPoint& crossPoint, const GeomDet& det, float window) const {
0216   // check if the z window around TSOS overlaps with the detector theDet (with a 1% margin added)
0217
0218   //   const float tolerance = 0.1;
0219   const float relativeMargin = 1.01;
0220
0221   LocalPoint localCrossPoint(det.surface().toLocal(crossPoint));
0222   //   if (std::abs(localCrossPoint.z()) > tolerance) {
0223   //     edm::LogInfo(TkDetLayers) << "Phase1PixelBlade::overlap calculation assumes point on surface, but it is off by "
0224   //     << localCrossPoint.z() ;
0225   //   }
0226
0227   float localX = localCrossPoint.x();
0228   float detHalfLength = det.surface().bounds().length() / 2.;
0229
0230   //   edm::LogInfo(TkDetLayers) << "Phase1PixelBlade::overlap: Det at " << det.position() << " hit at " << localY
0231   //        << " Window " << window << " halflength "  << detHalfLength ;
0232
0233   if ((std::abs(localX) - window) < relativeMargin * detHalfLength) {  // FIXME: margin hard-wired!
0234     return true;
0235   } else {
0236     return false;
0237   }
0238 }
0239
0240 int Phase1PixelBlade::findBin(float R, int diskSectorIndex) const {
0241   vector<const GeomDet*> localDets = diskSectorIndex == 0 ? theFrontDets : theBackDets;
0242
0243   int theBin = 0;
0244   float rDiff = std::abs(R - localDets.front()->surface().position().perp());
0245   for (vector<const GeomDet*>::const_iterator i = localDets.begin(); i != localDets.end(); i++) {
0246     float testDiff = std::abs(R - (**i).surface().position().perp());
0247     if (testDiff < rDiff) {
0248       rDiff = testDiff;
0249       theBin = i - localDets.begin();
0250     }
0251   }
0252   return theBin;
0253 }
0254
0255 int Phase1PixelBlade::findBin2(GlobalPoint thispoint, int diskSectorIndex) const {
0256   const vector<const GeomDet*>& localDets = diskSectorIndex == 0 ? theFrontDets : theBackDets;
0257
0258   int theBin = 0;
0259   float sDiff = (thispoint - localDets.front()->surface().position()).mag();
0260
0261   for (vector<const GeomDet*>::const_iterator i = localDets.begin(); i != localDets.end(); i++) {
0262     float testDiff = (thispoint - (**i).surface().position()).mag();
0263     if (testDiff < sDiff) {
0264       sDiff = testDiff;
0265       theBin = i - localDets.begin();
0266     }
0267   }
0268   return theBin;
0269 }
0270
0271 GlobalPoint Phase1PixelBlade::findPosition(int index, int diskSectorType) const {
0272   vector<const GeomDet*> diskSector = diskSectorType == 0 ? theFrontDets : theBackDets;
0273   return (diskSector[index])->surface().position();
0274 }
0275
0276 std::pair<float, float> Phase1PixelBlade::computeRadiusRanges(const std::vector<const GeomDet*>& current_dets) {
0277   typedef Surface::PositionType::BasicVectorType Vector;
0278   Vector posSum(0, 0, 0);
0279   for (auto i : current_dets)
0280     posSum += (*i).surface().position().basicVector();
0281
0282   Surface::PositionType meanPos(0., 0., posSum.z() / float(current_dets.size()));
0283
0284   // temporary plane - for the computation of bounds
0285   const Plane& tmpplane = current_dets.front()->surface();
0286
0287   GlobalVector xAxis;
0288   GlobalVector yAxis;
0289   GlobalVector zAxis;
0290
0291   GlobalVector planeXAxis = tmpplane.toGlobal(LocalVector(1, 0, 0));
0292   const GlobalPoint& planePosition = tmpplane.position();
0293
0294   if (planePosition.x() * planeXAxis.x() + planePosition.y() * planeXAxis.y() > 0.) {
0295     yAxis = planeXAxis;
0296   } else {
0297     yAxis = -planeXAxis;
0298   }
0299
0300   GlobalVector planeZAxis = tmpplane.toGlobal(LocalVector(0, 0, 1));
0301   if (planeZAxis.z() * planePosition.z() > 0.) {
0302     zAxis = planeZAxis;
0303   } else {
0304     zAxis = -planeZAxis;
0305   }
0306
0307   xAxis = yAxis.cross(zAxis);
0308
0309   Surface::RotationType rotation = Surface::RotationType(xAxis, yAxis);
0310   Plane plane(meanPos, rotation);
0311
0312   Surface::PositionType tmpPos = current_dets.front()->surface().position();
0313
0314   float rmin(plane.toLocal(tmpPos).perp());
0315   float rmax(plane.toLocal(tmpPos).perp());
0316
0317   for (auto it : current_dets) {
0318     vector<GlobalPoint> corners = BoundingBox().corners(it->specificSurface());
0319     for (const auto& i : corners) {
0320       float r = plane.toLocal(i).perp();
0321       rmin = min(rmin, r);
0322       rmax = max(rmax, r);
0323     }
0324   }
0325
0326   return std::make_pair(rmin, rmax);
0327 }