Phase1PixelBlade.cc

Go to the documentation of this file.

#include "Phase1PixelBlade.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "BladeShapeBuilderFromDet.h"
#include "LayerCrossingSide.h"
#include "DetGroupMerger.h"
#include "CompatibleDetToGroupAdder.h"
#include "DataFormats/GeometrySurface/interface/BoundingBox.h"

#include "TrackingTools/DetLayers/interface/DetLayerException.h"
#include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"

using namespace std;

typedef GeometricSearchDet::DetWithState DetWithState;

Phase1PixelBlade::~Phase1PixelBlade(){}

Phase1PixelBlade::Phase1PixelBlade(vector<const GeomDet*>& frontDets,
               vector<const GeomDet*>& backDets):
  GeometricSearchDet(true),            
    theFrontDets(frontDets),
    theBackDets(backDets),
    front_radius_range_(std::make_pair(0,0)),
    back_radius_range_(std::make_pair(0,0))
{
  theDets.assign(theFrontDets.begin(),theFrontDets.end());
  theDets.insert(theDets.end(),theBackDets.begin(),theBackDets.end());

  theDiskSector      = BladeShapeBuilderFromDet::build(theDets);
  theFrontDiskSector = BladeShapeBuilderFromDet::build(theFrontDets);
  theBackDiskSector  = BladeShapeBuilderFromDet::build(theBackDets);

  //--------- DEBUG INFO --------------
  LogDebug("TkDetLayers") << "DEBUG INFO for Phase1PixelBlade" ;
  LogDebug("TkDetLayers") << "Blade z, perp, innerRadius, outerR[disk, front, back]: "
              << this->position().z() << " , "
              << this->position().perp() << " , ("
              << theDiskSector->innerRadius() << " , "
              << theDiskSector->outerRadius() << "), ("
              << theFrontDiskSector->innerRadius() << " , "
              << theFrontDiskSector->outerRadius() << "), ("
              << theBackDiskSector->innerRadius() << " , "
              << theBackDiskSector->outerRadius() << ")" << std::endl;

  for(vector<const GeomDet*>::const_iterator it=theFrontDets.begin();
      it!=theFrontDets.end(); it++){
    LogDebug("TkDetLayers") << "frontDet phi,z,r: "
                << (*it)->position().phi() << " , "
                << (*it)->position().z()   << " , "
                << (*it)->position().perp() << " , "
                            << " rmin: " << (*it)->surface().rSpan().first
                            << " rmax: " << (*it)->surface().rSpan().second
                            << std::endl;
  }

  for(vector<const GeomDet*>::const_iterator it=theBackDets.begin();
      it!=theBackDets.end(); it++){
    LogDebug("TkDetLayers") << "backDet phi,z,r: "
                << (*it)->position().phi() << " , "
                << (*it)->position().z()   << " , "
                << (*it)->position().perp() << " , "
                            << " rmin: " << (*it)->surface().rSpan().first
                            << " rmax: " << (*it)->surface().rSpan().second
                            << std::endl;
  }

}


const vector<const GeometricSearchDet*>&
Phase1PixelBlade::components() const{
  throw DetLayerException("TOBRod doesn't have GeometricSearchDet components");
}

pair<bool, TrajectoryStateOnSurface>
Phase1PixelBlade::compatible( const TrajectoryStateOnSurface& ts, const Propagator&,
            const MeasurementEstimator&) const{
  edm::LogError("TkDetLayers") << "temporary dummy implementation of Phase1PixelBlade::compatible()!!" ;
  return pair<bool,TrajectoryStateOnSurface>();
}



void
Phase1PixelBlade::groupedCompatibleDetsV( const TrajectoryStateOnSurface& tsos,
                                    const Propagator& prop,
                                    const MeasurementEstimator& est,
                                    std::vector<DetGroup> & result) const{


  SubLayerCrossings  crossings;

  crossings = computeCrossings( tsos, prop.propagationDirection());
  if(! crossings.isValid()) return;

  vector<DetGroup> closestResult;
  addClosest( tsos, prop, est, crossings.closest(), closestResult);

  if (closestResult.empty()) {
    vector<DetGroup> nextResult;
    addClosest( tsos, prop, est, crossings.other(), nextResult);
    if(nextResult.empty())    return;

    //DetGroupElement nextGel( nextResult.front().front());
    //int crossingSide = LayerCrossingSide().endcapSide( nextGel.trajectoryState(), prop);

    DetGroupMerger::orderAndMergeTwoLevels( std::move(closestResult), std::move(nextResult), result,
                        0,0);//fixme gc patched for SLHC - already correctly sorted for SLHC
    //crossings.closestIndex(), crossingSide);
  }
  else {
    DetGroupElement closestGel( closestResult.front().front());
    float window = computeWindowSize( closestGel.det(), closestGel.trajectoryState(), est);

    searchNeighbors( tsos, prop, est, crossings.closest(), window,
             closestResult, false);

    vector<DetGroup> nextResult;
    searchNeighbors( tsos, prop, est, crossings.other(), window,
             nextResult, true);

    //int crossingSide = LayerCrossingSide().endcapSide( closestGel.trajectoryState(), prop);
    DetGroupMerger::orderAndMergeTwoLevels( std::move(closestResult), std::move(nextResult), result,
                        0,0);//fixme gc patched for SLHC - already correctly sorted for SLHC
    //crossings.closestIndex(), crossingSide);
  }
}

SubLayerCrossings
Phase1PixelBlade::computeCrossings( const TrajectoryStateOnSurface& startingState,
                  PropagationDirection propDir) const
{
  HelixPlaneCrossing::PositionType startPos( startingState.globalPosition());
  HelixPlaneCrossing::DirectionType startDir( startingState.globalMomentum());
  double rho( startingState.transverseCurvature());
  HelixArbitraryPlaneCrossing crossing( startPos, startDir, rho, propDir);

  pair<bool,double> innerPath = crossing.pathLength( *theFrontDiskSector);
  if (!innerPath.first) return SubLayerCrossings();

  GlobalPoint gInnerPoint( crossing.position(innerPath.second));
  //Code for use of binfinder
  //int innerIndex = theInnerBinFinder.binIndex(gInnerPoint.perp());
  //float innerDist = std::abs( theInnerBinFinder.binPosition(innerIndex) - gInnerPoint.z());

  //  int innerIndex = findBin(gInnerPoint.perp(),0);
  int innerIndex = findBin2(gInnerPoint,0);

  //fixme gc patched for SLHC - force order here to be in z
  //float innerDist = std::abs( findPosition(innerIndex,0).perp() - gInnerPoint.perp());
  float innerDist = ( startingState.globalPosition() - gInnerPoint).mag();
  SubLayerCrossing innerSLC( 0, innerIndex, gInnerPoint);

  pair<bool,double> outerPath = crossing.pathLength( *theBackDiskSector);
  if (!outerPath.first) return SubLayerCrossings();

  GlobalPoint gOuterPoint( crossing.position(outerPath.second));
  //Code for use of binfinder
  //int outerIndex = theOuterBinFinder.binIndex(gOuterPoint.perp());
  //float outerDist = std::abs( theOuterBinFinder.binPosition(outerIndex) - gOuterPoint.perp());
  //  int outerIndex  = findBin(gOuterPoint.perp(),1);
  int outerIndex  = findBin2(gOuterPoint,1);

  //fixme gc patched for SLHC - force order here to be in z
  //float outerDist = std::abs( findPosition(outerIndex,1).perp() - gOuterPoint.perp());
  float outerDist = ( startingState.globalPosition() - gOuterPoint).mag();
  SubLayerCrossing outerSLC( 1, outerIndex, gOuterPoint);

  if (innerDist < outerDist) {
    return SubLayerCrossings( innerSLC, outerSLC, 0);
  }
  else {
    return SubLayerCrossings( outerSLC, innerSLC, 1);
  }
}



bool
Phase1PixelBlade::addClosest( const TrajectoryStateOnSurface& tsos,
            const Propagator& prop,
            const MeasurementEstimator& est,
            const SubLayerCrossing& crossing,
            vector<DetGroup>& result) const
{

  const vector<const GeomDet*>& sBlade( subBlade( crossing.subLayerIndex()));

  return CompatibleDetToGroupAdder().add( *sBlade[crossing.closestDetIndex()],
                      tsos, prop, est, result);
}


float Phase1PixelBlade::computeWindowSize( const GeomDet* det,
                     const TrajectoryStateOnSurface& tsos,
                     const MeasurementEstimator& est) const
{
  return
      est.maximalLocalDisplacement(tsos, det->surface()).x();
}




void Phase1PixelBlade::searchNeighbors( const TrajectoryStateOnSurface& tsos,
                  const Propagator& prop,
                  const MeasurementEstimator& est,
                  const SubLayerCrossing& crossing,
                  float window,
                  vector<DetGroup>& result,
                  bool checkClosest) const
{
  const GlobalPoint& gCrossingPos = crossing.position();
  const vector<const GeomDet*>& sBlade( subBlade( crossing.subLayerIndex()));

  int closestIndex = crossing.closestDetIndex();
  int negStartIndex = closestIndex-1;
  int posStartIndex = closestIndex+1;

  if (checkClosest) { // must decide if the closest is on the neg or pos side
    if (gCrossingPos.perp() < sBlade[closestIndex]->surface().position().perp()) {
      posStartIndex = closestIndex;
    }
    else {
      negStartIndex = closestIndex;
    }
  }

  typedef CompatibleDetToGroupAdder Adder;
  for (int idet=negStartIndex; idet >= 0; idet--) {
    if (!overlap( gCrossingPos, *sBlade[idet], window)) break;
    if (!Adder::add( *sBlade[idet], tsos, prop, est, result)) break;
  }
  for (int idet=posStartIndex; idet < static_cast<int>(sBlade.size()); idet++) {
    if (!overlap( gCrossingPos, *sBlade[idet], window)) break;
    if (!Adder::add( *sBlade[idet], tsos, prop, est, result)) break;
  }
}



bool Phase1PixelBlade::overlap( const GlobalPoint& crossPoint, const GeomDet& det, float window) const
{
  // check if the z window around TSOS overlaps with the detector theDet (with a 1% margin added)

  //   const float tolerance = 0.1;
  const float relativeMargin = 1.01;

  LocalPoint localCrossPoint( det.surface().toLocal(crossPoint));
  //   if (std::abs(localCrossPoint.z()) > tolerance) {
  //     edm::LogInfo(TkDetLayers) << "Phase1PixelBlade::overlap calculation assumes point on surface, but it is off by "
  //     << localCrossPoint.z() ;
  //   }

  float localX = localCrossPoint.x();
  float detHalfLength = det.surface().bounds().length()/2.;

  //   edm::LogInfo(TkDetLayers) << "Phase1PixelBlade::overlap: Det at " << det.position() << " hit at " << localY
  //        << " Window " << window << " halflength "  << detHalfLength ;

  if ( ( std::abs(localX)-window) < relativeMargin*detHalfLength ) { // FIXME: margin hard-wired!
    return true;
  } else {
    return false;
  }
}

int
Phase1PixelBlade::findBin( float R,int diskSectorIndex) const
{
  vector<const GeomDet*> localDets = diskSectorIndex==0 ? theFrontDets : theBackDets;


  int theBin = 0;
  float rDiff = std::abs( R - localDets.front()->surface().position().perp());
  for (vector<const GeomDet*>::const_iterator i=localDets.begin(); i !=localDets.end(); i++){
    float testDiff = std::abs( R - (**i).surface().position().perp());
    if ( testDiff < rDiff) {
      rDiff = testDiff;
      theBin = i - localDets.begin();
    }
  }
  return theBin;
}

int
Phase1PixelBlade::findBin2( GlobalPoint thispoint,int diskSectorIndex) const
{
  const vector<const GeomDet*> & localDets = diskSectorIndex==0 ? theFrontDets : theBackDets;



  int theBin = 0;
  float sDiff = (thispoint - localDets.front()->surface().position()).mag();

  for (vector<const GeomDet*>::const_iterator i=localDets.begin(); i !=localDets.end(); i++){
    float testDiff = ( thispoint - (**i).surface().position()).mag();
    if ( testDiff < sDiff) {
      sDiff = testDiff;
      theBin = i - localDets.begin();
    }
  }
  return theBin;
}



GlobalPoint
Phase1PixelBlade::findPosition(int index,int diskSectorType) const
{
  vector<const GeomDet*> diskSector = diskSectorType == 0 ? theFrontDets : theBackDets;
  return (diskSector[index])->surface().position();
}

std::pair<float, float>
Phase1PixelBlade::computeRadiusRanges(const std::vector<const GeomDet*> &current_dets) {

  typedef Surface::PositionType::BasicVectorType Vector;
  Vector posSum(0,0,0);
  for (auto i : current_dets)
    posSum += (*i).surface().position().basicVector();

  Surface::PositionType meanPos( 0.,0.,posSum.z()/float(current_dets.size()));

  // temporary plane - for the computation of bounds
  const Plane& tmpplane = current_dets.front()->surface();

  GlobalVector xAxis;
  GlobalVector yAxis;
  GlobalVector zAxis;

  GlobalVector planeXAxis    = tmpplane.toGlobal(LocalVector(1, 0, 0));
  const GlobalPoint&  planePosition = tmpplane.position();

  if (planePosition.x() * planeXAxis.x()
     + planePosition.y() * planeXAxis.y() > 0.){
    yAxis = planeXAxis;
  } else {
    yAxis = -planeXAxis;
  }

  GlobalVector planeZAxis = tmpplane.toGlobal(LocalVector( 0, 0, 1));
  if (planeZAxis.z() * planePosition.z() > 0.) {
    zAxis = planeZAxis;
  } else {
    zAxis = -planeZAxis;
  }

  xAxis = yAxis.cross(zAxis);

  Surface::RotationType rotation = Surface::RotationType(xAxis, yAxis);
  Plane plane(meanPos, rotation);

  Surface::PositionType tmpPos = current_dets.front()->surface().position();

  float rmin(plane.toLocal(tmpPos).perp());
  float rmax(plane.toLocal(tmpPos).perp());

  for (auto it : current_dets) {
    vector<GlobalPoint> corners = BoundingBox().corners(it->specificSurface());
    for (auto i : corners) {
      float r   = plane.toLocal(i).perp();
      rmin = min(rmin, r);
      rmax = max(rmax, r);
    }
  }

  return std::make_pair(rmin, rmax);
}