TkDetLayers/src/TIBRing.cc

0001 #include "TIBRing.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0006 #include "TrackingTools/DetLayers/interface/CylinderBuilderFromDet.h"
0007 #include "TrackingTools/DetLayers/interface/simple_stat.h"
0008 #include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossing2OrderLocal.h"
0009 #include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
0010 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0011 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0012
0013 #include "LayerCrossingSide.h"
0014 #include "DetGroupMerger.h"
0015 #include "CompatibleDetToGroupAdder.h"
0016
0017 using namespace std;
0018
0019 typedef GeometricSearchDet::DetWithState DetWithState;
0020
0021 TIBRing::TIBRing(vector<const GeomDet*>& theGeomDets)
0022     : GeometricSearchDet(true), theDets(theGeomDets.begin(), theGeomDets.end()) {
0023   //checkRadius( first, last);
0024   //sort( theDets.begin(), theDets.end(), DetLessPhi());
0025   //checkPeriodicity( theDets.begin(), theDets.end());
0026
0027   theBinFinder = BinFinderType(theDets.front()->surface().position().phi(), theDets.size());
0028
0029   theCylinder = CylinderBuilderFromDet()(theDets.begin(), theDets.end());
0030
0031   computeHelicity();
0032
0033   LogDebug("TkDetLayers") << "==== DEBUG TIBRing =====";
0034   LogDebug("TkDetLayers") << "radius, thickness, lenght: " << theCylinder->radius() << " , "
0035                           << theCylinder->bounds().thickness() << " , " << theCylinder->bounds().length();
0036
0037   for (vector<const GeomDet*>::const_iterator i = theDets.begin(); i != theDets.end(); i++) {
0038     LogDebug("TkDetLayers") << "Ring's Det pos z,perp,eta,phi: " << (**i).position().z() << " , "
0039                             << (**i).position().perp() << " , " << (**i).position().eta() << " , "
0040                             << (**i).position().phi();
0041   }
0042   LogDebug("TkDetLayers") << "==== end DEBUG TIBRing =====";
0043 }
0044
0045 const vector<const GeometricSearchDet*>& TIBRing::components() const {
0046   throw DetLayerException("TIBRing doesn't have GeometricSearchDet components");
0047 }
0048
0049 void TIBRing::checkRadius(vector<const GeomDet*>::const_iterator first, vector<const GeomDet*>::const_iterator last) {
0050   // check radius range
0051   float rMin = 10000.;
0052   float rMax = 0.;
0053   for (vector<const GeomDet*>::const_iterator i = first; i != last; i++) {
0054     float r = (**i).surface().position().perp();
0055     if (r < rMin)
0056       rMin = r;
0057     if (r > rMax)
0058       rMax = r;
0059   }
0060   if (rMax - rMin > 0.1)
0061     throw DetLayerException("TIBRing construction failed: detectors not at constant radius");
0062 }
0063
0064 void TIBRing::checkPeriodicity(vector<const GeomDet*>::const_iterator first,
0065                                vector<const GeomDet*>::const_iterator last) {
0066   vector<double> adj_diff(last - first - 1);
0067   for (int i = 0; i < static_cast<int>(adj_diff.size()); i++) {
0068     vector<const GeomDet*>::const_iterator curent = first + i;
0069     adj_diff[i] = (**(curent + 1)).surface().position().phi() - (**curent).surface().position().phi();
0070   }
0071   double step = stat_mean(adj_diff);
0072   double phi_step = 2. * Geom::pi() / (last - first);
0073
0074   if (std::abs(step - phi_step) / phi_step > 0.01) {
0075     int ndets = last - first;
0076     edm::LogError("TkDetLayers") << "TIBRing Warning: not periodic. ndets=" << ndets;
0077     for (int j = 0; j < ndets; j++) {
0078       edm::LogError("TkDetLayers") << "Dets(r,phi): (" << theDets[j]->surface().position().perp() << ","
0079                                    << theDets[j]->surface().position().phi() << ") ";
0080     }
0081     throw DetLayerException("Error: TIBRing is not periodic");
0082   }
0083 }
0084
0085 void TIBRing::computeHelicity() {
0086   const GeomDet& det = *theDets.front();
0087   GlobalVector radial = det.surface().position() - GlobalPoint(0, 0, 0);
0088   GlobalVector normal = det.surface().toGlobal(LocalVector(0, 0, 1));
0089   if (normal.dot(radial) <= 0)
0090     normal *= -1;
0091   //   edm::LogInfo(TkDetLayers) << "BarrelDetRing::computeHelicity: phi(normal) " << normal.phi()
0092   //        << " phi(radial) " << radial.phi() ;
0093   if (Geom::phiLess(normal.phi(), radial.phi())) {
0094     theHelicity = 1;  // smaller phi angles mean "inner" group
0095   } else {
0096     theHelicity = 0;  // smaller phi angles mean "outer" group
0097   }
0098 }
0099
0100 TIBRing::~TIBRing() {}
0101
0102 pair<bool, TrajectoryStateOnSurface> TIBRing::compatible(const TrajectoryStateOnSurface& ts,
0103                                                          const Propagator&,
0104                                                          const MeasurementEstimator&) const {
0105   edm::LogError("TkDetLayers") << "temporary dummy implementation of TIBRing::compatible()!!";
0106   return pair<bool, TrajectoryStateOnSurface>();
0107 }
0108
0109 void TIBRing::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0110                                      const Propagator& prop,
0111                                      const MeasurementEstimator& est,
0112                                      vector<DetGroup>& result) const {
0113   vector<DetGroup> closestResult;
0114   auto crossings = computeCrossings(tsos, prop.propagationDirection());
0115   if (not crossings)
0116     return;
0117
0118   typedef CompatibleDetToGroupAdder Adder;
0119   Adder::add(*theDets[theBinFinder.binIndex(crossings->closestIndex)], tsos, prop, est, closestResult);
0120
0121   if (closestResult.empty()) {
0122     Adder::add(*theDets[theBinFinder.binIndex(crossings->nextIndex)], tsos, prop, est, result);
0123     return;
0124   }
0125
0126   DetGroupElement closestGel(closestResult.front().front());
0127   float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0128
0129   float detWidth = closestGel.det()->surface().bounds().width();
0130   if (crossings->nextDistance < detWidth + window) {
0131     vector<DetGroup> nextResult;
0132     if (Adder::add(*theDets[theBinFinder.binIndex(crossings->nextIndex)], tsos, prop, est, nextResult)) {
0133       int crossingSide = LayerCrossingSide::barrelSide(tsos, prop);
0134       if (crossings->closestIndex < crossings->nextIndex) {
0135         DetGroupMerger::orderAndMergeTwoLevels(
0136             std::move(closestResult), std::move(nextResult), result, theHelicity, crossingSide);
0137       } else {
0138         DetGroupMerger::orderAndMergeTwoLevels(
0139             std::move(nextResult), std::move(closestResult), result, theHelicity, crossingSide);
0140       }
0141     } else {
0142       result.swap(closestResult);
0143     }
0144   } else {
0145     result.swap(closestResult);
0146   }
0147
0148   // only loop over neighbors (other than closest and next) if window is BIG
0149   if (window > 0.5f * detWidth) {
0150     searchNeighbors(tsos, prop, est, *crossings, window, result);
0151   }
0152 }
0153
0154 void TIBRing::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0155                               const Propagator& prop,
0156                               const MeasurementEstimator& est,
0157                               const SubRingCrossings& crossings,
0158                               float window,
0159                               vector<DetGroup>& result) const {
0160   typedef CompatibleDetToGroupAdder Adder;
0161   int crossingSide = LayerCrossingSide::barrelSide(tsos, prop);
0162   typedef DetGroupMerger Merger;
0163
0164   int negStart = std::min(crossings.closestIndex, crossings.nextIndex) - 1;
0165   int posStart = std::max(crossings.closestIndex, crossings.nextIndex) + 1;
0166
0167   int quarter = theDets.size() / 4;
0168   for (int idet = negStart; idet >= negStart - quarter + 1; idet--) {
0169     const GeomDet* neighbor = theDets[theBinFinder.binIndex(idet)];
0170     // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0171     // maybe also add shallow crossing angle test here???
0172     vector<DetGroup> tmp1;
0173     if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0174       break;
0175     vector<DetGroup> tmp2;
0176     tmp2.swap(result);
0177     vector<DetGroup> newResult;
0178     Merger::orderAndMergeTwoLevels(std::move(tmp1), std::move(tmp2), newResult, theHelicity, crossingSide);
0179     result.swap(newResult);
0180   }
0181   for (int idet = posStart; idet < posStart + quarter - 1; idet++) {
0182     const GeomDet* neighbor = theDets[theBinFinder.binIndex(idet)];
0183     // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0184     // maybe also add shallow crossing angle test here???
0185     vector<DetGroup> tmp1;
0186     if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0187       break;
0188     vector<DetGroup> tmp2;
0189     tmp2.swap(result);
0190     vector<DetGroup> newResult;
0191     Merger::orderAndMergeTwoLevels(std::move(tmp2), std::move(tmp1), newResult, theHelicity, crossingSide);
0192     result.swap(newResult);
0193   }
0194 }
0195
0196 std::optional<TIBRing::SubRingCrossings> TIBRing::computeCrossings(const TrajectoryStateOnSurface& startingState,
0197                                                                    PropagationDirection propDir) const {
0198   typedef HelixBarrelPlaneCrossing2OrderLocal Crossing;
0199   typedef MeasurementEstimator::Local2DVector Local2DVector;
0200
0201   GlobalPoint startPos(startingState.globalPosition());
0202   GlobalVector startDir(startingState.globalMomentum());
0203   auto rho = startingState.transverseCurvature();
0204
0205   HelixBarrelCylinderCrossing cylCrossing(
0206       startPos, startDir, rho, propDir, specificSurface(), HelixBarrelCylinderCrossing::bestSol);
0207
0208   if (!cylCrossing.hasSolution())
0209     return std::nullopt;
0210
0211   GlobalPoint cylPoint(cylCrossing.position());
0212   GlobalVector cylDir(cylCrossing.direction());
0213   int closestIndex = theBinFinder.binIndex(cylPoint.barePhi());
0214
0215   const Plane& closestPlane(theDets[closestIndex]->surface());
0216
0217   LocalPoint closestPos = Crossing::positionOnly(cylPoint, cylDir, rho, closestPlane);
0218   float closestDist = closestPos.x();  // use fact that local X perp to global Z
0219
0220   //int next = cylPoint.phi() - closestPlane.position().phi() > 0 ? closest+1 : closest-1;
0221   int nextIndex =
0222       Geom::phiLess(closestPlane.position().barePhi(), cylPoint.barePhi()) ? closestIndex + 1 : closestIndex - 1;
0223
0224   const Plane& nextPlane(theDets[theBinFinder.binIndex(nextIndex)]->surface());
0225   LocalPoint nextPos = Crossing::positionOnly(cylPoint, cylDir, rho, nextPlane);
0226   float nextDist = nextPos.x();
0227
0228   if (std::abs(closestDist) < std::abs(nextDist)) {
0229     return SubRingCrossings(closestIndex, nextIndex, nextDist);
0230   } else {
0231     return SubRingCrossings(nextIndex, closestIndex, closestDist);
0232   }
0233 }
0234
0235 float TIBRing::computeWindowSize(const GeomDet* det,
0236                                  const TrajectoryStateOnSurface& tsos,
0237                                  const MeasurementEstimator& est) const {
0238   return est.maximalLocalDisplacement(tsos, det->surface()).x();
0239 }