TkNavigation/plugins/SimpleNavigableLayer.cc

0001 #include "SimpleNavigableLayer.h"
0002 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
0003
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0007 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0008 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0009 #include "FWCore/Utilities/interface/Likely.h"
0010
0011 #include <set>
0012
0013 using namespace std;
0014
0015 TrajectoryStateOnSurface SimpleNavigableLayer::crossingState(const FreeTrajectoryState& fts,
0016                                                              PropagationDirection dir) const {
0017   //self propagating. step one: go close to the center
0018   GlobalPoint initialPoint = fts.position();
0019   TransverseImpactPointExtrapolator middle;
0020   GlobalPoint center(0, 0, 0);
0021   TSOS propState = middle.extrapolate(fts, center, propagator(dir));
0022   if (!propState.isValid())
0023     return TrajectoryStateOnSurface();
0024
0025   FreeTrajectoryState const& dest = *propState.freeState();
0026   GlobalPoint middlePoint = dest.position();
0027   const float toCloseToEachOther2 = 1e-4 * 1e-4;
0028   if UNLIKELY ((middlePoint - initialPoint).mag2() < toCloseToEachOther2) {
0029     LogDebug("SimpleNavigableLayer")
0030         << "initial state and PCA are identical. Things are bound to fail. Do not add the link.";
0031     return TrajectoryStateOnSurface();
0032   }
0033
0034   /*
0035   std::string dirS;
0036   if (dir==alongMomentum) dirS = "alongMomentum";
0037   else if (dir==oppositeToMomentum) dirS = "oppositeToMomentum";
0038   else dirS = "anyDirection";
0039   */
0040
0041   LogDebug("SimpleNavigableLayer") << "self propagating(" << dir << ") from:\n"
0042                                    << fts << "\n"
0043                                    << dest << "\n"
0044                                    << " and the direction is: " << dir;
0045
0046   //second propagation to go on the other side of the barrel
0047   //propState = propagator(dir).propagate( dest, detLayer()->specificSurface());
0048   propState = propagator(dir).propagate(dest, detLayer()->surface());
0049   if (!propState.isValid())
0050     return TrajectoryStateOnSurface();
0051
0052   const FreeTrajectoryState& dest2 = *propState.freeState();
0053   GlobalPoint finalPoint = dest2.position();
0054   LogDebug("SimpleNavigableLayer") << "second propagation(" << dir << ") to: \n" << dest2;
0055   double finalDot = (middlePoint - initialPoint).basicVector().dot((finalPoint - middlePoint).basicVector());
0056   if UNLIKELY (finalDot < 0) {  // check that before and after are in different side.
0057     LogDebug("SimpleNavigableLayer") << "switch side back: ABORT.";
0058     return TrajectoryStateOnSurface();
0059   }
0060   return propState;
0061 }
0062
0063 bool SimpleNavigableLayer::wellInside(const FreeTrajectoryState& fts,
0064                                       PropagationDirection dir,
0065                                       const BarrelDetLayer* bl,
0066                                       DLC& result) const {
0067   TSOS propState = (bl == detLayer()) ? crossingState(fts, dir) : propagator(dir).propagate(fts, bl->specificSurface());
0068
0069   if (!propState.isValid())
0070     return false;
0071
0072   //if requested check that the layer is crossed on the right side
0073   if (theCheckCrossingSide) {
0074     bool backTobackTransverse =
0075         (fts.position().x() * propState.globalPosition().x() + fts.position().y() * propState.globalPosition().y()) < 0;
0076     bool backToback = propState.globalPosition().basicVector().dot(fts.position().basicVector()) < 0;
0077
0078     if (backTobackTransverse || backToback) {
0079       LogTrace("TkNavigation") << "Crossing over prevented!\nStaring from (x,y,z,r) (" << fts.position().x() << ","
0080                                << fts.position().y() << "," << fts.position().z() << "," << fts.position().perp()
0081                                << ") going to TSOS (x,y,z,r)" << propState.globalPosition().x() << ","
0082                                << propState.globalPosition().y() << "," << propState.globalPosition().z() << ","
0083                                << propState.globalPosition().perp() << ")";
0084       return false;
0085
0086       /*
0087     //we have to check the crossing side only if we are going to something smaller
0088     if (fts.position().perp()>bl->specificSurface().radius() ||
0089     fabs(fts.position().z())>bl->surface().bounds().length()/2. ){
0090     if (propState.globalPosition().basicVector().dot(fts.position().basicVector())<0){
0091     LogTrace("TkNavigation") << "Crossing over prevented!\nStaring from (x,y,z,r) ("
0092     << fts.position().x()<<","<< fts.position().y()<<","<< fts.position().z()<<","<<fts.position().perp()
0093     << ") going to TSOS (x,y,z,r)"
0094     << propState.globalPosition().x()<<","<< propState.globalPosition().y()<<","<< propState.globalPosition().z()<<","<<propState.globalPosition().perp()<<")";;
0095     return false;
0096     }
0097     }
0098     */
0099     }
0100   }
0101
0102   const Bounds& bounds(bl->specificSurface().bounds());
0103   float length = bounds.length() * 0.5f;
0104
0105   // take into account the thickness of the layer
0106   float deltaZ =
0107       0.5f * bounds.thickness() * std::abs(propState.globalDirection().z()) / propState.globalDirection().perp();
0108
0109   // take into account the error on the predicted state
0110   const float nSigma = theEpsilon;  // temporary reuse of epsilon
0111   if (propState.hasError()) {
0112     deltaZ += nSigma * sqrt(fts.cartesianError().position().czz());
0113   }
0114
0115   // cout << "SimpleNavigableLayer BarrelDetLayer deltaZ = " << deltaZ << endl;
0116
0117   float zpos = propState.globalPosition().z();
0118   if (std::abs(zpos) < length + deltaZ)
0119     result.push_back(bl);
0120
0121   return std::abs(zpos) < length - deltaZ;
0122 }
0123
0124 bool SimpleNavigableLayer::wellInside(const FreeTrajectoryState& fts,
0125                                       PropagationDirection dir,
0126                                       const ForwardDetLayer* fl,
0127                                       DLC& result) const {
0128   TSOS propState = propagator(dir).propagate(fts, fl->specificSurface());
0129   if (!propState.isValid())
0130     return false;
0131
0132   if (fl == detLayer()) {
0133     LogDebug("SimpleNavigableLayer") << "self propagating from:\n" << fts << "\n to \n" << *propState.freeState();
0134   }
0135
0136   //if requested avoids crossing over the tracker
0137   if (theCheckCrossingSide) {
0138     bool backTobackTransverse =
0139         (fts.position().x() * propState.globalPosition().x() + fts.position().y() * propState.globalPosition().y()) < 0;
0140     bool backToback = propState.globalPosition().basicVector().dot(fts.position().basicVector()) < 0;
0141
0142     if (backTobackTransverse || backToback) {
0143       LogTrace("TkNavigation") << "Crossing over prevented!\nStaring from (x,y,z,r) (" << fts.position().x() << ","
0144                                << fts.position().y() << "," << fts.position().z() << "," << fts.position().perp()
0145                                << ") going to TSOS (x,y,z,r)" << propState.globalPosition().x() << ","
0146                                << propState.globalPosition().y() << "," << propState.globalPosition().z() << ","
0147                                << propState.globalPosition().perp() << ")";
0148       ;
0149       return false;
0150
0151       //    if (fts.position().z()*propState.globalPosition().z() < 0) return false;
0152     }
0153   }
0154
0155   float rpos = propState.globalPosition().perp();
0156   float innerR = fl->specificSurface().innerRadius();
0157   float outerR = fl->specificSurface().outerRadius();
0158
0159   // take into account the thickness of the layer
0160   float deltaR = 0.5f * fl->surface().bounds().thickness() * propState.localDirection().perp() /
0161                  std::abs(propState.localDirection().z());
0162
0163   // take into account the error on the predicted state
0164   const float nSigma = theEpsilon;
0165   if (propState.hasError()) {
0166     LocalError err = propState.localError().positionError();
0167     // ignore correlation for the moment...
0168     deltaR += nSigma * sqrt(err.xx() + err.yy());
0169   }
0170
0171   // cout << "SimpleNavigableLayer BarrelDetLayer deltaR = " << deltaR << endl;
0172
0173   if (innerR - deltaR < rpos && rpos < outerR + deltaR)
0174     result.push_back(fl);
0175   return (innerR + deltaR < rpos && rpos < outerR - deltaR);
0176 }
0177
0178 bool SimpleNavigableLayer::wellInside(const FreeTrajectoryState& fts,
0179                                       PropagationDirection dir,
0180                                       const DLC& layers,
0181                                       DLC& result) const {
0182   for (auto l : layers) {
0183     if (l->isBarrel()) {
0184       const BarrelDetLayer* bl = reinterpret_cast<const BarrelDetLayer*>(l);
0185       if (wellInside(fts, dir, bl, result))
0186         return true;
0187     } else {
0188       const ForwardDetLayer* fl = reinterpret_cast<const ForwardDetLayer*>(l);
0189       if (wellInside(fts, dir, fl, result))
0190         return true;
0191     }
0192   }
0193   return false;
0194 }
0195
0196 bool SimpleNavigableLayer::wellInside(
0197     const FreeTrajectoryState& fts, PropagationDirection dir, ConstBDLI begin, ConstBDLI end, DLC& result) const {
0198   for (ConstBDLI i = begin; i < end; i++) {
0199     if (wellInside(fts, dir, *i, result))
0200       return true;
0201   }
0202   return false;
0203 }
0204
0205 bool SimpleNavigableLayer::wellInside(
0206     const FreeTrajectoryState& fts, PropagationDirection dir, ConstFDLI begin, ConstFDLI end, DLC& result) const {
0207   for (ConstFDLI i = begin; i < end; i++) {
0208     if (wellInside(fts, dir, *i, result))
0209       return true;
0210   }
0211   return false;
0212 }
0213
0214 std::vector<const DetLayer*> SimpleNavigableLayer::compatibleLayers(const FreeTrajectoryState& fts,
0215                                                                     PropagationDirection timeDirection,
0216                                                                     int& counter) const {
0217   typedef std::vector<const DetLayer*> Lvect;
0218   typedef std::set<const DetLayer*> Lset;
0219
0220   //initiate the first iteration
0221   const Lvect& someLayers = nextLayers(fts, timeDirection);
0222   if (someLayers.empty()) {
0223     LogDebug("SimpleNavigableLayer") << "Number of compatible layers: " << 0;
0224     return someLayers;
0225   }
0226
0227   Lset collect;                     //a container of unique instances. to avoid duplicates
0228   Lset layerToTry, nextLayerToTry;  //set used for iterations
0229   layerToTry.insert(someLayers.begin(), someLayers.end());
0230
0231   while (!layerToTry.empty() && (counter++) <= 150) {
0232     LogDebug("SimpleNavigableLayer") << counter << "] going to check on : " << layerToTry.size() << " next layers.";
0233     //clear this set first, it will be swaped with layerToTry
0234     nextLayerToTry.clear();
0235     for (auto toTry : layerToTry) {
0236       //add the layer you tried.
0237       LogDebug("SimpleNavigableLayer") << counter << "] adding layer with pointer: " << (toTry)
0238                                        << " first detid: " << (toTry)->basicComponents().front()->geographicalId();
0239       if (!collect.insert(toTry).second)
0240         continue;
0241
0242       //find the next layers from it
0243       Lvect&& nextLayers = school->nextLayers(*toTry, fts, timeDirection);
0244       LogDebug("SimpleNavigableLayer") << counter << "] this layer has : " << nextLayers.size() << " next layers.";
0245       nextLayerToTry.insert(nextLayers.begin(), nextLayers.end());
0246     }  // layerToTry
0247     //swap now that you where to go next.
0248     layerToTry.swap(nextLayerToTry);
0249   }
0250   if (counter >= 150) {
0251     edm::LogWarning("SimpleNavigableLayer") << "WARNING: compatibleLayers() more than 150 iterations!!! Bailing out..";
0252     counter = -1;
0253     return Lvect();
0254   }
0255
0256   LogDebug("SimpleNavigableLayer") << "Number of compatible layers: " << collect.size();
0257
0258   return Lvect(collect.begin(), collect.end());
0259 }