ForwardDetLayer.cc

Go to the documentation of this file.

#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" 
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
#include "DataFormats/GeometrySurface/interface/BoundingBox.h"

using namespace std;


ForwardDetLayer::~ForwardDetLayer() {
}


//--- Extension of the interface
void ForwardDetLayer::setSurface( BoundDisk* cp) { 
  theDisk = cp;
}

bool ForwardDetLayer::contains(const Local3DPoint& p) const {
  return surface().bounds().inside(p);
}


void ForwardDetLayer::initialize() {
  setSurface( computeSurface());
}


BoundDisk* ForwardDetLayer::computeSurface() {
  LogDebug("DetLayers") << "ForwaLayer::computeSurface callded" ;
  vector<const GeomDet*> comps= basicComponents();

  vector<const GeomDet*>::const_iterator ifirst = comps.begin();
  vector<const GeomDet*>::const_iterator ilast  = comps.end();

  // Find extension in R
  float theRmin = components().front()->position().perp(); 
  float theRmax = theRmin;
  float theZmin = components().back()->position().z();
  float theZmax = theZmin;
  for ( vector<const GeomDet*>::const_iterator deti = ifirst;
    deti != ilast; deti++) {
    vector<GlobalPoint> corners = 
      BoundingBox().corners( dynamic_cast<const Plane&>((**deti).surface()));
    for (vector<GlobalPoint>::const_iterator ic = corners.begin();
     ic != corners.end(); ic++) {
      float r = ic->perp();
      LogDebug("DetLayers") << "corner.perp(): " << r ;
      float z = ic->z();
      theRmin = min( theRmin, r);
      theRmax = max( theRmax, r);
      theZmin = min( theZmin, z);
      theZmax = max( theZmax, z);
    }

    // in addition to the corners we have to check the middle of the 
    // det +/- length/2
    // , since the min (max) radius for typical fw dets is reached there

    float rdet  = (**deti).position().perp();
    float len   = (**deti).surface().bounds().length();
    float width = (**deti).surface().bounds().width();

    GlobalVector xAxis = (**deti).toGlobal(LocalVector(1,0,0));
    GlobalVector yAxis = (**deti).toGlobal(LocalVector(0,1,0));
    GlobalVector perpDir = GlobalVector( (**deti).position() - GlobalPoint(0,0,(**deti).position().z()) );

    double xAxisCos = xAxis.unit().dot(perpDir.unit());
    double yAxisCos = yAxis.unit().dot(perpDir.unit());

    LogDebug("DetLayers") << "in ForwardDetLayer::computeSurface(),xAxisCos,yAxisCos: " << xAxisCos << " , " << yAxisCos ;
    LogDebug("DetLayers") << "det pos.perp,length,width: " 
              << rdet << " , " 
              << len  << " , "
              << width ;

    if( fabs(xAxisCos) > fabs(yAxisCos) ) {
      theRmin = min( theRmin, rdet-width/2.F);
      theRmax = max( theRmax, rdet+width/2.F);
    }else{
      theRmin = min( theRmin, rdet-len/2.F);
      theRmax = max( theRmax, rdet+len/2.F);
    }
  }


  LogDebug("DetLayers") << "creating SimpleDiskBounds with r range" << theRmin << " " 
            << theRmax << " and z range " << theZmin << " " << theZmax ;

  // By default the forward layers are positioned around the z axis of the
  // global frame, and the axes of their local frame coincide with 
  // those of the global grame (z along the disk axis)
  float zPos = (theZmax+theZmin)/2.;
  PositionType pos(0.,0.,zPos);
  RotationType rot;

  return new BoundDisk( pos, rot, 
            new SimpleDiskBounds( theRmin, theRmax, 
                              theZmin-zPos, theZmax-zPos));
}  


pair<bool, TrajectoryStateOnSurface>
ForwardDetLayer::compatible( const TrajectoryStateOnSurface& ts, 
                 const Propagator& prop, 
                 const MeasurementEstimator&) const
{
  if unlikely(theDisk == nullptr)  edm::LogError("DetLayers") 
    << "ERROR: BarrelDetLayer::compatible() is used before the layer surface is initialized" ;
  // throw an exception? which one?

  TrajectoryStateOnSurface myState = prop.propagate( ts, specificSurface());
  if unlikely( !myState.isValid()) return make_pair( false, myState);


  // check z;  (are we sure?????)
  // auto z = myState.localPosition().z();
  // if ( z<bounds().minZ | z>bounds().maxZ) return make_pair( false, myState);

  // check r
  auto r2 =  myState.localPosition().perp2();
  if ( (r2 > rmin()*rmin()) & (r2< rmax()*rmax()) ) return make_pair( true, myState);


  // take into account the thickness of the layer
  float deltaR = 0.5f*bounds().thickness() *
    myState.localDirection().perp()/std::abs(myState.localDirection().z());

  // and take into account the error on the predicted state
  const float nSigma = 3.;
  if (myState.hasError()) {
    LocalError err = myState.localError().positionError();
    // ignore correlation for the moment...
    deltaR += nSigma * sqrt(err.xx() + err.yy());
  }

  // check r again;
  auto ri2 = std::max(rmin()-deltaR,0.f); ri2*=ri2;
  auto ro2 = rmax()+deltaR; ro2*=ro2;
  return make_pair( (r2>ri2) & (r2<ro2), myState);
}