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   SubRingCrossings crossings;
0115   crossings = computeCrossings(tsos, prop.propagationDirection());
0116   if (!crossings.isValid_)
0117     return;
0118
0119   typedef CompatibleDetToGroupAdder Adder;
0120   Adder::add(*theDets[theBinFinder.binIndex(crossings.closestIndex)], tsos, prop, est, closestResult);
0121
0122   if (closestResult.empty()) {
0123     Adder::add(*theDets[theBinFinder.binIndex(crossings.nextIndex)], tsos, prop, est, result);
0124     return;
0125   }
0126
0127   DetGroupElement closestGel(closestResult.front().front());
0128   float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0129
0130   float detWidth = closestGel.det()->surface().bounds().width();
0131   if (crossings.nextDistance < detWidth + window) {
0132     vector<DetGroup> nextResult;
0133     if (Adder::add(*theDets[theBinFinder.binIndex(crossings.nextIndex)], tsos, prop, est, nextResult)) {
0134       int crossingSide = LayerCrossingSide::barrelSide(tsos, prop);
0135       if (crossings.closestIndex < crossings.nextIndex) {
0136         DetGroupMerger::orderAndMergeTwoLevels(
0137             std::move(closestResult), std::move(nextResult), result, theHelicity, crossingSide);
0138       } else {
0139         DetGroupMerger::orderAndMergeTwoLevels(
0140             std::move(nextResult), std::move(closestResult), result, theHelicity, crossingSide);
0141       }
0142     } else {
0143       result.swap(closestResult);
0144     }
0145   } else {
0146     result.swap(closestResult);
0147   }
0148
0149   // only loop over neighbors (other than closest and next) if window is BIG
0150   if (window > 0.5f * detWidth) {
0151     searchNeighbors(tsos, prop, est, crossings, window, result);
0152   }
0153 }
0154
0155 void TIBRing::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0156                               const Propagator& prop,
0157                               const MeasurementEstimator& est,
0158                               const SubRingCrossings& crossings,
0159                               float window,
0160                               vector<DetGroup>& result) const {
0161   typedef CompatibleDetToGroupAdder Adder;
0162   int crossingSide = LayerCrossingSide::barrelSide(tsos, prop);
0163   typedef DetGroupMerger Merger;
0164
0165   int negStart = std::min(crossings.closestIndex, crossings.nextIndex) - 1;
0166   int posStart = std::max(crossings.closestIndex, crossings.nextIndex) + 1;
0167
0168   int quarter = theDets.size() / 4;
0169   for (int idet = negStart; idet >= negStart - quarter + 1; idet--) {
0170     const GeomDet* neighbor = theDets[theBinFinder.binIndex(idet)];
0171     // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0172     // maybe also add shallow crossing angle test here???
0173     vector<DetGroup> tmp1;
0174     if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0175       break;
0176     vector<DetGroup> tmp2;
0177     tmp2.swap(result);
0178     vector<DetGroup> newResult;
0179     Merger::orderAndMergeTwoLevels(std::move(tmp1), std::move(tmp2), newResult, theHelicity, crossingSide);
0180     result.swap(newResult);
0181   }
0182   for (int idet = posStart; idet < posStart + quarter - 1; idet++) {
0183     const GeomDet* neighbor = theDets[theBinFinder.binIndex(idet)];
0184     // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0185     // maybe also add shallow crossing angle test here???
0186     vector<DetGroup> tmp1;
0187     if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0188       break;
0189     vector<DetGroup> tmp2;
0190     tmp2.swap(result);
0191     vector<DetGroup> newResult;
0192     Merger::orderAndMergeTwoLevels(std::move(tmp2), std::move(tmp1), newResult, theHelicity, crossingSide);
0193     result.swap(newResult);
0194   }
0195 }
0196
0197 TIBRing::SubRingCrossings TIBRing::computeCrossings(const TrajectoryStateOnSurface& startingState,
0198                                                     PropagationDirection propDir) const {
0199   typedef HelixBarrelPlaneCrossing2OrderLocal Crossing;
0200   typedef MeasurementEstimator::Local2DVector Local2DVector;
0201
0202   GlobalPoint startPos(startingState.globalPosition());
0203   GlobalVector startDir(startingState.globalMomentum());
0204   auto rho = startingState.transverseCurvature();
0205
0206   HelixBarrelCylinderCrossing cylCrossing(
0207       startPos, startDir, rho, propDir, specificSurface(), HelixBarrelCylinderCrossing::bestSol);
0208
0209   if (!cylCrossing.hasSolution())
0210     return SubRingCrossings();
0211
0212   GlobalPoint cylPoint(cylCrossing.position());
0213   GlobalVector cylDir(cylCrossing.direction());
0214   int closestIndex = theBinFinder.binIndex(cylPoint.barePhi());
0215
0216   const Plane& closestPlane(theDets[closestIndex]->surface());
0217
0218   LocalPoint closestPos = Crossing::positionOnly(cylPoint, cylDir, rho, closestPlane);
0219   float closestDist = closestPos.x();  // use fact that local X perp to global Z
0220
0221   //int next = cylPoint.phi() - closestPlane.position().phi() > 0 ? closest+1 : closest-1;
0222   int nextIndex =
0223       Geom::phiLess(closestPlane.position().barePhi(), cylPoint.barePhi()) ? closestIndex + 1 : closestIndex - 1;
0224
0225   const Plane& nextPlane(theDets[theBinFinder.binIndex(nextIndex)]->surface());
0226   LocalPoint nextPos = Crossing::positionOnly(cylPoint, cylDir, rho, nextPlane);
0227   float nextDist = nextPos.x();
0228
0229   if (std::abs(closestDist) < std::abs(nextDist)) {
0230     return SubRingCrossings(closestIndex, nextIndex, nextDist);
0231   } else {
0232     return SubRingCrossings(nextIndex, closestIndex, closestDist);
0233   }
0234 }
0235
0236 float TIBRing::computeWindowSize(const GeomDet* det,
0237                                  const TrajectoryStateOnSurface& tsos,
0238                                  const MeasurementEstimator& est) const {
0239   return est.maximalLocalDisplacement(tsos, det->surface()).x();
0240 }