DetLayers/src/ForwardDetLayer.cc

0001 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0005 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
0006 #include "FWCore/Utilities/interface/Likely.h"
0007
0008 using namespace std;
0009
0010 ForwardDetLayer::~ForwardDetLayer() {}
0011
0012 //--- Extension of the interface
0013 void ForwardDetLayer::setSurface(BoundDisk* cp) { theDisk = cp; }
0014
0015 bool ForwardDetLayer::contains(const Local3DPoint& p) const { return surface().bounds().inside(p); }
0016
0017 void ForwardDetLayer::initialize() { setSurface(computeSurface()); }
0018
0019 BoundDisk* ForwardDetLayer::computeSurface() {
0020   LogDebug("DetLayers") << "ForwaLayer::computeSurface callded";
0021   vector<const GeomDet*> comps = basicComponents();
0022
0023   vector<const GeomDet*>::const_iterator ifirst = comps.begin();
0024   vector<const GeomDet*>::const_iterator ilast = comps.end();
0025
0026   // Find extension in R
0027   float theRmin = components().front()->position().perp();
0028   float theRmax = theRmin;
0029   float theZmin = components().back()->position().z();
0030   float theZmax = theZmin;
0031   for (vector<const GeomDet*>::const_iterator deti = ifirst; deti != ilast; deti++) {
0032     vector<GlobalPoint> corners = BoundingBox().corners(dynamic_cast<const Plane&>((**deti).surface()));
0033     for (vector<GlobalPoint>::const_iterator ic = corners.begin(); ic != corners.end(); ic++) {
0034       float r = ic->perp();
0035       LogDebug("DetLayers") << "corner.perp(): " << r;
0036       float z = ic->z();
0037       theRmin = min(theRmin, r);
0038       theRmax = max(theRmax, r);
0039       theZmin = min(theZmin, z);
0040       theZmax = max(theZmax, z);
0041     }
0042
0043     // in addition to the corners we have to check the middle of the
0044     // det +/- length/2
0045     // , since the min (max) radius for typical fw dets is reached there
0046
0047     float rdet = (**deti).position().perp();
0048     float len = (**deti).surface().bounds().length();
0049     float width = (**deti).surface().bounds().width();
0050
0051     GlobalVector xAxis = (**deti).toGlobal(LocalVector(1, 0, 0));
0052     GlobalVector yAxis = (**deti).toGlobal(LocalVector(0, 1, 0));
0053     GlobalVector perpDir = GlobalVector((**deti).position() - GlobalPoint(0, 0, (**deti).position().z()));
0054
0055     double xAxisCos = xAxis.unit().dot(perpDir.unit());
0056     double yAxisCos = yAxis.unit().dot(perpDir.unit());
0057
0058     LogDebug("DetLayers") << "in ForwardDetLayer::computeSurface(),xAxisCos,yAxisCos: " << xAxisCos << " , "
0059                           << yAxisCos;
0060     LogDebug("DetLayers") << "det pos.perp,length,width: " << rdet << " , " << len << " , " << width;
0061
0062     if (fabs(xAxisCos) > fabs(yAxisCos)) {
0063       theRmin = min(theRmin, rdet - width / 2.F);
0064       theRmax = max(theRmax, rdet + width / 2.F);
0065     } else {
0066       theRmin = min(theRmin, rdet - len / 2.F);
0067       theRmax = max(theRmax, rdet + len / 2.F);
0068     }
0069   }
0070
0071   LogDebug("DetLayers") << "creating SimpleDiskBounds with r range" << theRmin << " " << theRmax << " and z range "
0072                         << theZmin << " " << theZmax;
0073
0074   // By default the forward layers are positioned around the z axis of the
0075   // global frame, and the axes of their local frame coincide with
0076   // those of the global grame (z along the disk axis)
0077   float zPos = (theZmax + theZmin) / 2.;
0078   PositionType pos(0., 0., zPos);
0079   RotationType rot;
0080
0081   return new BoundDisk(pos, rot, new SimpleDiskBounds(theRmin, theRmax, theZmin - zPos, theZmax - zPos));
0082 }
0083
0084 pair<bool, TrajectoryStateOnSurface> ForwardDetLayer::compatible(const TrajectoryStateOnSurface& ts,
0085                                                                  const Propagator& prop,
0086                                                                  const MeasurementEstimator&) const {
0087   if UNLIKELY (theDisk == nullptr)
0088     edm::LogError("DetLayers") << "ERROR: BarrelDetLayer::compatible() is used before the layer surface is initialized";
0089   // throw an exception? which one?
0090
0091   TrajectoryStateOnSurface myState = prop.propagate(ts, specificSurface());
0092   if UNLIKELY (!myState.isValid())
0093     return make_pair(false, myState);
0094
0095   // check z;  (are we sure?????)
0096   // auto z = myState.localPosition().z();
0097   // if ( z<bounds().minZ | z>bounds().maxZ) return make_pair( false, myState);
0098
0099   // check r
0100   auto r2 = myState.localPosition().perp2();
0101
0102   // Disable bitwise-instead-of-logical warning, see discussion in
0103   // https://github.com/cms-sw/cmssw/issues/39105
0104
0105 #if defined(__clang__) && defined(__has_warning)
0106 #if __has_warning("-Wbitwise-instead-of-logical")
0107 #pragma clang diagnostic push
0108 #pragma clang diagnostic ignored "-Wbitwise-instead-of-logical"
0109 #endif
0110 #endif
0111
0112   if ((r2 > rmin() * rmin()) & (r2 < rmax() * rmax()))
0113     return make_pair(true, myState);
0114
0115 #if defined(__clang__) && defined(__has_warning)
0116 #if __has_warning("-Wbitwise-instead-of-logical")
0117 #pragma clang diagnostic pop
0118 #endif
0119 #endif
0120
0121   // take into account the thickness of the layer
0122   float deltaR = 0.5f * bounds().thickness() * myState.localDirection().perp() / std::abs(myState.localDirection().z());
0123
0124   // and take into account the error on the predicted state
0125   const float nSigma = 3.;
0126   if (myState.hasError()) {
0127     LocalError err = myState.localError().positionError();
0128     // ignore correlation for the moment...
0129     deltaR += nSigma * sqrt(err.xx() + err.yy());
0130   }
0131
0132   // check r again;
0133   auto ri2 = std::max(rmin() - deltaR, 0.f);
0134   ri2 *= ri2;
0135   auto ro2 = rmax() + deltaR;
0136   ro2 *= ro2;
0137   return make_pair((r2 > ri2) & (r2 < ro2), myState);
0138 }