TkDetUtil.cc

Go to the documentation of this file.

#include "TkDetUtil.h"

#include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "DataFormats/GeometryVector/interface/VectorUtil.h"

namespace tkDetUtil {

  float computeWindowSize( const GeomDet* det, 
               const TrajectoryStateOnSurface& tsos, 
               const MeasurementEstimator& est)
  {
    const Plane& startPlane = det->surface();  
    auto maxDistance =  est.maximalLocalDisplacement( tsos, startPlane);
    return std::copysign(calculatePhiWindow(maxDistance, tsos, startPlane),maxDistance.x()); 
  }


  float 
  calculatePhiWindow( const MeasurementEstimator::Local2DVector & imaxDistance, 
              const TrajectoryStateOnSurface& ts, 
              const Plane& plane)
  {
    MeasurementEstimator::Local2DVector maxDistance(std::abs(imaxDistance.x()),std::abs(imaxDistance.y()));

    constexpr float tolerance=1.e-6;
    LocalPoint start = ts.localPosition();
    //     std::cout << "plane z " << plane.normalVector() << std::endl;
    float dphi=0;
    if LIKELY(std::abs(1.f-std::abs(plane.normalVector().z()))<tolerance) {
      auto ori = plane.toLocal(GlobalPoint(0.,0.,0.));
      auto xc = std::abs(start.x() - ori.x());
      auto yc = std::abs(start.y() - ori.y());
      
      if (yc<maxDistance.y() && xc<maxDistance.x()) return M_PI;

      auto hori = yc>maxDistance.y(); // quadrant 1 (&2), otherwiase quadrant 1&4
      auto y0 = hori ? yc + std::copysign(maxDistance.y(), xc - maxDistance.x()) : xc - maxDistance.x();
      auto x0 = hori ? xc - maxDistance.x() : -yc - maxDistance.y();
      auto y1 = hori ? yc - maxDistance.y() :  xc - maxDistance.x();
      auto x1 = hori ? xc + maxDistance.x() : -yc + maxDistance.y();

      auto sp = (x0*x1+y0*y1)/std::sqrt((x0*x0+y0*y0)*(x1*x1+y1*y1));
      sp = std::min(std::max(sp,-1.f),1.f);
      dphi = std::acos(sp);      

      return dphi;
    }
    

    // generic algo
    float corners[]  =  { plane.toGlobal(LocalPoint( start.x()+maxDistance.x(), start.y()+maxDistance.y() )).barePhi(),
              plane.toGlobal(LocalPoint( start.x()-maxDistance.x(), start.y()+maxDistance.y() )).barePhi(),
              plane.toGlobal(LocalPoint( start.x()-maxDistance.x(), start.y()-maxDistance.y() )).barePhi(),
              plane.toGlobal(LocalPoint( start.x()+maxDistance.x(), start.y()-maxDistance.y() )).barePhi() 
    };
    
    float phimin = corners[0];
    float phimax = phimin;
    for ( int i = 1; i<4; i++) {
      float cPhi = corners[i];
      if ( Geom::phiLess(cPhi, phimin) ) { phimin = cPhi; }
      if ( Geom::phiLess( phimax, cPhi) ) { phimax = cPhi; }
    }
    float phiWindow = phimax - phimin;
    if ( phiWindow < 0.) { phiWindow +=  2.*Geom::pi();}
    // std::cout << "phiWindow " << phiWindow << ' ' << dphi << ' ' << dphi-phiWindow  << std::endl;
    return phiWindow;
  }




}