File indexing completed on 2021-02-14 14:27:46
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
0024
0025
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
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
0092
0093 if (Geom::phiLess(normal.phi(), radial.phi())) {
0094 theHelicity = 1;
0095 } else {
0096 theHelicity = 0;
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
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
0172
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
0185
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();
0220
0221
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 }