TkDetLayers/src/TIDLayer.cc

0001 #include "TIDLayer.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0006
0007 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0008 #include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
0009
0010 #include <array>
0011 #include "DetGroupMerger.h"
0012
0013 //#include "CommonDet/DetLayout/src/DetLessR.h"
0014
0015 using namespace std;
0016
0017 typedef GeometricSearchDet::DetWithState DetWithState;
0018
0019 namespace {
0020
0021   // 2 0 1
0022   /*
0023 class TIDringLess  {
0024   // z-position ordering of TID rings indices
0025 public:
0026   bool operator()( const pair<vector<DetGroup> const *,int> & a,const pair<vector<DetGroup>const *,int> & b) const {
0027     if(a.second==2) {return true;}
0028     else if(a.second==0) {
0029       if(b.second==2) return false;
0030       if(b.second==1) return true;
0031     }
0032     return false;
0033   };
0034 };
0035 */
0036
0037   // groups in correct order: one may be empty
0038   inline void mergeOutward(std::array<vector<DetGroup>, 3>& groups, std::vector<DetGroup>& result) {
0039     typedef DetGroupMerger Merger;
0040     Merger::orderAndMergeTwoLevels(std::move(groups[0]), std::move(groups[1]), result, 1, 1);
0041     if (!groups[2].empty()) {
0042       std::vector<DetGroup> tmp;
0043       tmp.swap(result);
0044       Merger::orderAndMergeTwoLevels(std::move(tmp), std::move(groups[2]), result, 1, 1);
0045     }
0046   }
0047
0048   inline void mergeInward(std::array<vector<DetGroup>, 3>& groups, std::vector<DetGroup>& result) {
0049     typedef DetGroupMerger Merger;
0050     Merger::orderAndMergeTwoLevels(std::move(groups[2]), std::move(groups[1]), result, 1, 1);
0051     if (!groups[0].empty()) {
0052       std::vector<DetGroup> tmp;
0053       tmp.swap(result);
0054       Merger::orderAndMergeTwoLevels(std::move(tmp), std::move(groups[0]), result, 1, 1);
0055     }
0056   }
0057
0058   void orderAndMergeLevels(const TrajectoryStateOnSurface& tsos,
0059                            const Propagator& prop,
0060                            std::array<vector<DetGroup>, 3>& groups,
0061                            std::vector<DetGroup>& result) {
0062     float zpos = tsos.globalPosition().z();
0063     if (tsos.globalMomentum().z() * zpos > 0) {  // momentum points outwards
0064       if (prop.propagationDirection() == alongMomentum)
0065         mergeOutward(groups, result);
0066       else
0067         mergeInward(groups, result);
0068     } else {  //  momentum points inwards
0069       if (prop.propagationDirection() == oppositeToMomentum)
0070         mergeOutward(groups, result);
0071       else
0072         mergeInward(groups, result);
0073     }
0074   }
0075
0076 }  // namespace
0077
0078 //hopefully is never called!
0079 const std::vector<const GeometricSearchDet*>& TIDLayer::components() const {
0080   if (not theComponents) {
0081     auto temp = std::make_unique<std::vector<const GeometricSearchDet*>>();
0082     temp->reserve(3);
0083     for (auto c : theComps)
0084       temp->push_back(c);
0085     std::vector<const GeometricSearchDet*>* expected = nullptr;
0086     if (theComponents.compare_exchange_strong(expected, temp.get())) {
0087       //this thread set the value
0088       temp.release();
0089     }
0090   }
0091
0092   return *theComponents;
0093 }
0094
0095 void TIDLayer::fillRingPars(int i) {
0096   const BoundDisk& ringDisk = static_cast<const BoundDisk&>(theComps[i]->surface());
0097   float ringMinZ = std::abs(ringDisk.position().z()) - ringDisk.bounds().thickness() / 2.;
0098   float ringMaxZ = std::abs(ringDisk.position().z()) + ringDisk.bounds().thickness() / 2.;
0099   ringPars[i].thetaRingMin = ringDisk.innerRadius() / ringMaxZ;
0100   ringPars[i].thetaRingMax = ringDisk.outerRadius() / ringMinZ;
0101   ringPars[i].theRingR = (ringDisk.innerRadius() + ringDisk.outerRadius()) / 2.;
0102 }
0103
0104 TIDLayer::TIDLayer(vector<const TIDRing*>& rings) : RingedForwardLayer(true), theComponents{nullptr} {
0105   //They should be already R-ordered. TO BE CHECKED!!
0106   //sort( theRings.begin(), theRings.end(), DetLessR());
0107
0108   if (rings.size() != 3)
0109     throw DetLayerException("Number of rings in TID layer is not equal to 3 !!");
0110   setSurface(computeDisk(rings));
0111
0112   for (int i = 0; i != 3; ++i) {
0113     theComps[i] = rings[i];
0114     fillRingPars(i);
0115     theBasicComps.insert(
0116         theBasicComps.end(), (*rings[i]).basicComponents().begin(), (*rings[i]).basicComponents().end());
0117   }
0118
0119   LogDebug("TkDetLayers") << "==== DEBUG TIDLayer =====";
0120   LogDebug("TkDetLayers") << "r,zed pos  , thickness, innerR, outerR: " << this->position().perp() << " , "
0121                           << this->position().z() << " , " << this->specificSurface().bounds().thickness() << " , "
0122                           << this->specificSurface().innerRadius() << " , " << this->specificSurface().outerRadius();
0123 }
0124
0125 BoundDisk* TIDLayer::computeDisk(const vector<const TIDRing*>& rings) const {
0126   float theRmin = rings.front()->specificSurface().innerRadius();
0127   float theRmax = rings.front()->specificSurface().outerRadius();
0128   float theZmin = rings.front()->position().z() - rings.front()->surface().bounds().thickness() / 2;
0129   float theZmax = rings.front()->position().z() + rings.front()->surface().bounds().thickness() / 2;
0130
0131   for (vector<const TIDRing*>::const_iterator i = rings.begin(); i != rings.end(); i++) {
0132     float rmin = (**i).specificSurface().innerRadius();
0133     float rmax = (**i).specificSurface().outerRadius();
0134     float zmin = (**i).position().z() - (**i).surface().bounds().thickness() / 2.;
0135     float zmax = (**i).position().z() + (**i).surface().bounds().thickness() / 2.;
0136     theRmin = min(theRmin, rmin);
0137     theRmax = max(theRmax, rmax);
0138     theZmin = min(theZmin, zmin);
0139     theZmax = max(theZmax, zmax);
0140   }
0141
0142   float zPos = (theZmax + theZmin) / 2.;
0143   PositionType pos(0., 0., zPos);
0144   RotationType rot;
0145
0146   return new BoundDisk(pos, rot, new SimpleDiskBounds(theRmin, theRmax, theZmin - zPos, theZmax - zPos));
0147 }
0148
0149 TIDLayer::~TIDLayer() {
0150   for (auto c : theComps)
0151     delete c;
0152
0153   delete theComponents.load();
0154 }
0155
0156 void TIDLayer::groupedCompatibleDetsV(const TrajectoryStateOnSurface& startingState,
0157                                       const Propagator& prop,
0158                                       const MeasurementEstimator& est,
0159                                       std::vector<DetGroup>& result) const {
0160   std::array<int, 3> const& ringIndices = ringIndicesByCrossingProximity(startingState, prop);
0161   if (ringIndices[0] == -1) {
0162     edm::LogError("TkDetLayers") << "TkRingedForwardLayer::groupedCompatibleDets : error in CrossingProximity";
0163     return;
0164   }
0165
0166   std::array<vector<DetGroup>, 3> groupsAtRingLevel;
0167   //order is ring3,ring1,ring2 i.e. 2 0 1
0168   //                                0 1 2
0169   static constexpr int ringOrder[3]{1, 2, 0};
0170   auto index = [&ringIndices](int i) { return ringOrder[ringIndices[i]]; };
0171
0172   auto& closestResult = groupsAtRingLevel[index(0)];
0173   theComps[ringIndices[0]]->groupedCompatibleDetsV(startingState, prop, est, closestResult);
0174   if (closestResult.empty()) {
0175     theComps[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, result);
0176     return;
0177   }
0178
0179   DetGroupElement closestGel(closestResult.front().front());
0180   float rWindow = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0181
0182   if (!overlapInR(closestGel.trajectoryState(), ringIndices[1], rWindow)) {
0183     result.swap(closestResult);
0184     return;
0185   }
0186
0187   auto& nextResult = groupsAtRingLevel[index(1)];
0188   theComps[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, nextResult);
0189   if (nextResult.empty()) {
0190     result.swap(closestResult);
0191     return;
0192   }
0193
0194   if (!overlapInR(closestGel.trajectoryState(), ringIndices[2], rWindow)) {
0195     //then merge 2 levels & return
0196     orderAndMergeLevels(closestGel.trajectoryState(), prop, groupsAtRingLevel, result);
0197     return;
0198   }
0199
0200   auto& nextNextResult = groupsAtRingLevel[index(2)];
0201   theComps[ringIndices[2]]->groupedCompatibleDetsV(startingState, prop, est, nextNextResult);
0202   if (nextNextResult.empty()) {
0203     // then merge 2 levels and return
0204     orderAndMergeLevels(closestGel.trajectoryState(), prop, groupsAtRingLevel, result);  //
0205     return;
0206   }
0207
0208   // merge 3 level and return merged
0209   orderAndMergeLevels(closestGel.trajectoryState(), prop, groupsAtRingLevel, result);
0210 }
0211
0212 std::array<int, 3> TIDLayer::ringIndicesByCrossingProximity(const TrajectoryStateOnSurface& startingState,
0213                                                             const Propagator& prop) const {
0214   typedef HelixForwardPlaneCrossing Crossing;
0215   typedef MeasurementEstimator::Local2DVector Local2DVector;
0216
0217   HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0218   HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0219   PropagationDirection propDir(prop.propagationDirection());
0220   float rho(startingState.transverseCurvature());
0221
0222   // calculate the crossings with the ring surfaces
0223   // rings are assumed to be sorted in R !
0224
0225   Crossing myXing(startPos, startDir, rho, propDir);
0226
0227   GlobalPoint ringCrossings[3];
0228   // vector<GlobalVector>  ringXDirections;
0229
0230   for (int i = 0; i < 3; i++) {
0231     const BoundDisk& theRing = static_cast<const BoundDisk&>(theComps[i]->surface());
0232     pair<bool, double> pathlen = myXing.pathLength(theRing);
0233     if (pathlen.first) {
0234       ringCrossings[i] = GlobalPoint(myXing.position(pathlen.second));
0235       // ringXDirections.push_back( GlobalVector( myXing.direction(pathlen.second )));
0236     } else {
0237       // TO FIX.... perhaps there is something smarter to do
0238       //throw DetLayerException("trajectory doesn't cross TID rings");
0239       ringCrossings[i] = GlobalPoint(0., 0., 0.);
0240       //  ringXDirections.push_back( GlobalVector( 0.,0.,0.));
0241     }
0242   }
0243
0244   int closestIndex = findClosest(ringCrossings);
0245   int nextIndex = findNextIndex(ringCrossings, closestIndex);
0246   if (closestIndex < 0 || nextIndex < 0)
0247     return std::array<int, 3>{{-1, -1, -1}};
0248   int nextNextIndex = -1;
0249   for (int i = 0; i < 3; i++) {
0250     if (i != closestIndex && i != nextIndex) {
0251       nextNextIndex = i;
0252       break;
0253     }
0254   }
0255
0256   std::array<int, 3> indices{{closestIndex, nextIndex, nextNextIndex}};
0257   return indices;
0258 }
0259
0260 float TIDLayer::computeWindowSize(const GeomDet* det,
0261                                   const TrajectoryStateOnSurface& tsos,
0262                                   const MeasurementEstimator& est) const {
0263   const Plane& startPlane = det->surface();
0264   MeasurementEstimator::Local2DVector maxDistance = est.maximalLocalDisplacement(tsos, startPlane);
0265   return maxDistance.y();
0266 }
0267
0268 int TIDLayer::findClosest(const GlobalPoint ringCrossing[3]) const {
0269   int theBin = 0;
0270   float initialR = ringPars[0].theRingR;
0271   float rDiff = std::abs(ringCrossing[0].perp() - initialR);
0272   for (int i = 1; i < 3; i++) {
0273     float ringR = ringPars[i].theRingR;
0274     float testDiff = std::abs(ringCrossing[i].perp() - ringR);
0275     if (testDiff < rDiff) {
0276       rDiff = testDiff;
0277       theBin = i;
0278     }
0279   }
0280   return theBin;
0281 }
0282
0283 int TIDLayer::findNextIndex(const GlobalPoint ringCrossing[3], int closest) const {
0284   int firstIndexToCheck = (closest != 0) ? 0 : 1;
0285   float initialR = ringPars[firstIndexToCheck].theRingR;
0286   float rDiff = std::abs(ringCrossing[firstIndexToCheck].perp() - initialR);
0287   int theBin = firstIndexToCheck;
0288   for (int i = firstIndexToCheck + 1; i < 3; i++) {
0289     if (i != closest) {
0290       float ringR = ringPars[i].theRingR;
0291       float testDiff = std::abs(ringCrossing[i].perp() - ringR);
0292       if (testDiff < rDiff) {
0293         rDiff = testDiff;
0294         theBin = i;
0295       }
0296     }
0297   }
0298   return theBin;
0299 }
0300
0301 bool TIDLayer::overlapInR(const TrajectoryStateOnSurface& tsos, int index, double ymax) const {
0302   // assume "fixed theta window", i.e. margin in local y = r is changing linearly with z
0303   float tsRadius = tsos.globalPosition().perp();
0304   float thetamin = (max(0., tsRadius - ymax)) / (std::abs(tsos.globalPosition().z()) + 10.f);  // add 10 cm contingency
0305   float thetamax = (tsRadius + ymax) / (std::abs(tsos.globalPosition().z()) - 10.f);
0306
0307   // do the theta regions overlap ?
0308
0309   return !(thetamin > ringPars[index].thetaRingMax || ringPars[index].thetaRingMin > thetamax);
0310 }