VectorHitBuilderAlgorithm Class Reference

#include <VectorHitBuilderAlgorithm.h>

Inheritance diagram for VectorHitBuilderAlgorithm:

Classes
struct	CurvatureAndPhi

Public Member Functions
VectorHit	buildVectorHit (const StackGeomDet *stack, Phase2TrackerCluster1DRef lower, Phase2TrackerCluster1DRef upper) const override
void	buildVectorHits (VectorHitCollection &vhAcc, VectorHitCollection &vhRej, DetId detIdStack, const StackGeomDet *stack, edm::Handle< edmNew::DetSetVector< Phase2TrackerCluster1D >> clusters, const Detset &DSVinner, const Detset &DSVouter, const std::vector< bool > &phase2OTClustersToSkip=std::vector< bool >()) const override
bool	checkClustersCompatibility (Local3DPoint &posinner, Local3DPoint &posouter, LocalError &errinner, LocalError &errouter) const
bool	checkClustersCompatibilityBeforeBuilding (edm::Handle< edmNew::DetSetVector< Phase2TrackerCluster1D >> clusters, const Detset &theLowerDetSet, const Detset &theUpperDetSet) const
CurvatureAndPhi	curvatureANDphi (Global3DPoint gPositionLower, Global3DPoint gPositionUpper, GlobalError gErrorLower, GlobalError gErrorUpper) const
void	fit (float x[2], float y[2], float sigy[2], Local3DPoint &pos, Local3DVector &dir, AlgebraicSymMatrix22 &covMatrix, double &chi2) const
void	fit2Dzx (const Local3DPoint lpCI, const Local3DPoint lpCO, const LocalError leCI, const LocalError leCO, Local3DPoint &pos, Local3DVector &dir, AlgebraicSymMatrix22 &covMatrix, double &chi2) const
void	fit2Dzy (const Local3DPoint lpCI, const Local3DPoint lpCO, const LocalError leCI, const LocalError leCO, Local3DPoint &pos, Local3DVector &dir, AlgebraicSymMatrix22 &covMatrix, double &chi2) const
void	run (edm::Handle< edmNew::DetSetVector< Phase2TrackerCluster1D >> clusters, VectorHitCollection &vhAcc, VectorHitCollection &vhRej, edmNew::DetSetVector< Phase2TrackerCluster1D > &clustersAcc, edmNew::DetSetVector< Phase2TrackerCluster1D > &clustersRej) const override
	VectorHitBuilderAlgorithm (const edm::ParameterSet &conf, const TrackerGeometry *tkGeomProd, const TrackerTopology *tkTopoProd, const ClusterParameterEstimator< Phase2TrackerCluster1D > *cpeProd)
	~VectorHitBuilderAlgorithm () override=default
Public Member Functions inherited from VectorHitBuilderAlgorithmBase
double	computeParallaxCorrection (const PixelGeomDetUnit *, const Point3DBase< float, LocalTag > &, const PixelGeomDetUnit *, const Point3DBase< float, LocalTag > &) const
void	printCluster (const GeomDet *geomDetUnit, const Phase2TrackerCluster1D *cluster) const
void	printClusters (const edmNew::DetSetVector< Phase2TrackerCluster1D > &clusters) const
	VectorHitBuilderAlgorithmBase (const edm::ParameterSet &, const TrackerGeometry *, const TrackerTopology *, const ClusterParameterEstimator< Phase2TrackerCluster1D > *)
virtual	~VectorHitBuilderAlgorithmBase ()

Additional Inherited Members
Public Types inherited from VectorHitBuilderAlgorithmBase
typedef Detset::const_iterator	const_iterator
typedef edmNew::DetSet< Phase2TrackerCluster1D >	Detset
typedef edmNew::DetSetVector< VectorHit >	output_t
typedef edm::Ref< edmNew::DetSetVector< Phase2TrackerCluster1D >, Phase2TrackerCluster1D >	Phase2TrackerCluster1DRef
typedef std::pair< StackGeomDet, std::vector< Phase2TrackerCluster1D > >	StackClusters
Public Attributes inherited from VectorHitBuilderAlgorithmBase
std::vector< double >	barrelCut_
const ClusterParameterEstimator< Phase2TrackerCluster1D > *	cpe_
std::vector< double >	endcapCut_
unsigned int	nMaxVHforeachStack_
const TrackerGeometry *	tkGeom_
const TrackerTopology *	tkTopo_

Detailed Description

Definition at line 17 of file VectorHitBuilderAlgorithm.h.

Constructor & Destructor Documentation

VectorHitBuilderAlgorithm()

VectorHitBuilderAlgorithm::VectorHitBuilderAlgorithm ( const edm::ParameterSet &  conf, const TrackerGeometry *  tkGeomProd, const TrackerTopology *  tkTopoProd, const ClusterParameterEstimator< Phase2TrackerCluster1D > *  cpeProd  )

inline

Definition at line 19 of file VectorHitBuilderAlgorithm.h.

      : VectorHitBuilderAlgorithmBase(conf, tkGeomProd, tkTopoProd, cpeProd){};

~VectorHitBuilderAlgorithm()

VectorHitBuilderAlgorithm::~VectorHitBuilderAlgorithm ( )

overridedefault

Member Function Documentation

buildVectorHit()

VectorHit VectorHitBuilderAlgorithm::buildVectorHit ( const StackGeomDet *  stack, Phase2TrackerCluster1DRef  lower, Phase2TrackerCluster1DRef  upper  ) const

overridevirtual

Implements VectorHitBuilderAlgorithmBase.

Definition at line 202 of file VectorHitBuilderAlgorithm.cc.

References checkClustersCompatibility(), VectorHitBuilderAlgorithmBase::cpe_, VectorHitBuilderAlgorithm::CurvatureAndPhi::curvature, curvatureANDphi(), VectorHitBuilderAlgorithm::CurvatureAndPhi::curvatureError, fit2Dzx(), fit2Dzy(), ClusterParameterEstimator< T >::localParameters(), LogTrace, convertSQLiteXML::ok, PV3DBase< T, PVType, FrameType >::perp(), VectorHit::phase2clusterGlobalPos(), VectorHit::phase2clusterGlobalPosErr(), VectorHitBuilderAlgorithm::CurvatureAndPhi::phi, VectorHitBuilderAlgorithmBase::printCluster(), svgfig::stack, GeomDet::surface(), edm::swap(), and Surface::toGlobal().

Referenced by buildVectorHits().

                                                                                           {
  LogTrace("VectorHitBuilderAlgorithm") << "Build VH with: ";
#ifdef EDM_ML_DEBUG
  printCluster(stack->upperDet(), &*upper);
#endif
  const PixelGeomDetUnit* geomDetLower = static_cast<const PixelGeomDetUnit*>(stack->lowerDet());
  const PixelGeomDetUnit* geomDetUpper = static_cast<const PixelGeomDetUnit*>(stack->upperDet());

  auto&& lparamsLower = cpe_->localParameters(*lower, *geomDetLower);  // x, y, z, e2_xx, e2_xy, e2_yy
  Global3DPoint gparamsLower = geomDetLower->surface().toGlobal(lparamsLower.first);
  LogTrace("VectorHitBuilderAlgorithm") << "\t lower global pos: " << gparamsLower;

  auto&& lparamsUpper = cpe_->localParameters(*upper, *geomDetUpper);
  Global3DPoint gparamsUpper = geomDetUpper->surface().toGlobal(lparamsUpper.first);
  LogTrace("VectorHitBuilderAlgorithm") << "\t upper global pos: " << gparamsUpper;

  //local parameters of upper cluster in lower system of reference
  Local3DPoint lparamsUpperInLower = geomDetLower->surface().toLocal(gparamsUpper);

  LogTrace("VectorHitBuilderAlgorithm") << "\t lower global pos: " << gparamsLower;
  LogTrace("VectorHitBuilderAlgorithm") << "\t upper global pos: " << gparamsUpper;

  LogTrace("VectorHitBuilderAlgorithm") << "A:\t lower local pos: " << lparamsLower.first
                                        << " with error: " << lparamsLower.second << std::endl;
  LogTrace("VectorHitBuilderAlgorithm") << "A:\t upper local pos in the lower sof " << lparamsUpperInLower
                                        << " with error: " << lparamsUpper.second << std::endl;

  bool ok =
      checkClustersCompatibility(lparamsLower.first, lparamsUpper.first, lparamsLower.second, lparamsUpper.second);

  if (ok) {
    AlgebraicSymMatrix22 covMat2Dzx;
    double chi22Dzx = 0.0;
    Local3DPoint pos2Dzx;
    Local3DVector dir2Dzx;
    fit2Dzx(lparamsLower.first,
            lparamsUpperInLower,
            lparamsLower.second,
            lparamsUpper.second,
            pos2Dzx,
            dir2Dzx,
            covMat2Dzx,
            chi22Dzx);
    LogTrace("VectorHitBuilderAlgorithm") << "\t  pos2Dzx: " << pos2Dzx;
    LogTrace("VectorHitBuilderAlgorithm") << "\t  dir2Dzx: " << dir2Dzx;
    LogTrace("VectorHitBuilderAlgorithm") << "\t  cov2Dzx: " << covMat2Dzx;
    VectorHit2D vh2Dzx = VectorHit2D(pos2Dzx, dir2Dzx, covMat2Dzx, chi22Dzx);

    AlgebraicSymMatrix22 covMat2Dzy;
    double chi22Dzy = 0.0;
    Local3DPoint pos2Dzy;
    Local3DVector dir2Dzy;
    fit2Dzy(lparamsLower.first,
            lparamsUpperInLower,
            lparamsLower.second,
            lparamsUpper.second,
            pos2Dzy,
            dir2Dzy,
            covMat2Dzy,
            chi22Dzy);
    LogTrace("VectorHitBuilderAlgorithm") << "\t  pos2Dzy: " << pos2Dzy;
    LogTrace("VectorHitBuilderAlgorithm") << "\t  dir2Dzy: " << dir2Dzy;
    LogTrace("VectorHitBuilderAlgorithm") << "\t  cov2Dzy: " << covMat2Dzy;
    VectorHit2D vh2Dzy = VectorHit2D(pos2Dzy, dir2Dzy, covMat2Dzy, chi22Dzy);

    OmniClusterRef lowerOmni(lower);
    OmniClusterRef upperOmni(upper);

    Global3DPoint gPositionLower = VectorHit::phase2clusterGlobalPos(geomDetLower, lower);
    Global3DPoint gPositionUpper = VectorHit::phase2clusterGlobalPos(geomDetUpper, upper);
    GlobalError gErrorLower = VectorHit::phase2clusterGlobalPosErr(geomDetLower);
    GlobalError gErrorUpper = VectorHit::phase2clusterGlobalPosErr(geomDetUpper);

    if (gPositionLower.perp() > gPositionUpper.perp()) {
      std::swap(gPositionLower, gPositionUpper);
      std::swap(gErrorLower, gErrorUpper);
    }

    const CurvatureAndPhi curvatureAndPhi = curvatureANDphi(gPositionLower, gPositionUpper, gErrorLower, gErrorUpper);
    VectorHit vh = VectorHit(*stack,
                             vh2Dzx,
                             vh2Dzy,
                             lowerOmni,
                             upperOmni,
                             curvatureAndPhi.curvature,
                             curvatureAndPhi.curvatureError,
                             curvatureAndPhi.phi);
    return vh;
  }

  return VectorHit();
}

buildVectorHits()

void VectorHitBuilderAlgorithm::buildVectorHits ( VectorHitCollection &  vhAcc, VectorHitCollection &  vhRej, DetId  detIdStack, const StackGeomDet *  stack, edm::Handle< edmNew::DetSetVector< Phase2TrackerCluster1D >>  clusters, const Detset &  DSVinner, const Detset &  DSVouter, const std::vector< bool > &  phase2OTClustersToSkip = std::vector<bool>()  ) const

overridevirtual

Implements VectorHitBuilderAlgorithmBase.

Definition at line 79 of file VectorHitBuilderAlgorithm.cc.

References funct::abs(), VectorHitBuilderAlgorithmBase::barrelCut_, edmNew::DetSet< T >::begin(), buildVectorHit(), checkClustersCompatibilityBeforeBuilding(), bsc_activity_cfg::clusters, VectorHitBuilderAlgorithmBase::computeParallaxCorrection(), VectorHitBuilderAlgorithmBase::cpe_, DMR_cfg::cut, edmNew::DetSet< T >::end(), VectorHitBuilderAlgorithmBase::endcapCut_, TrackingRecHit::isValid(), TrackerTopology::layer(), ClusterParameterEstimator< T >::localParameters(), LogDebug, LogTrace, edmNew::makeRefTo(), GeomDetEnumerators::P2OTB, GeomDetEnumerators::P2OTEC, VectorHitBuilderAlgorithmBase::printCluster(), edmNew::DetSetVector< T >::FastFiller::push_back(), svgfig::stack, VectorHitBuilderAlgorithmBase::tkTopo_, and ApeEstimator_cff::width.

Referenced by run().

                                                                                                     {
  if (checkClustersCompatibilityBeforeBuilding(clusters, theLowerDetSet, theUpperDetSet)) {
    LogDebug("VectorHitBuilderAlgorithm") << "  compatible -> continue ... " << std::endl;
  } else {
    LogTrace("VectorHitBuilderAlgorithm") << "  not compatible, going to the next cluster";
  }

  edmNew::DetSetVector<VectorHit>::FastFiller vh_colAcc(vhAcc, detIdStack);
  edmNew::DetSetVector<VectorHit>::FastFiller vh_colRej(vhRej, detIdStack);

  unsigned int layerStack = tkTopo_->layer(stack->geographicalId());
  if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTB)
    LogDebug("VectorHitBuilderAlgorithm") << " \t is barrel.    " << std::endl;
  if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTEC)
    LogDebug("VectorHitBuilderAlgorithm") << " \t is endcap.    " << std::endl;
  LogDebug("VectorHitBuilderAlgorithm") << " \t layer is : " << layerStack << std::endl;

  float cut = 0.0;
  if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTB)
    cut = barrelCut_.at(layerStack);
  if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTEC)
    cut = endcapCut_.at(layerStack);
  LogDebug("VectorHitBuilderAlgorithm") << " \t the cut is:" << cut << std::endl;

  //only cache local parameters for upper cluster as we loop over lower clusters only once anyway
  std::vector<std::pair<LocalPoint, LocalError>> localParamsUpper;
  std::vector<const PixelGeomDetUnit*> localGDUUpper;

  const PixelGeomDetUnit* gduUpp = dynamic_cast<const PixelGeomDetUnit*>(stack->upperDet());
  for (auto const& clusterUpper : theUpperDetSet) {
    localGDUUpper.push_back(gduUpp);
    localParamsUpper.push_back(cpe_->localParameters(clusterUpper, *gduUpp));
  }
  int upperIterator = 0;
  const PixelGeomDetUnit* gduLow = dynamic_cast<const PixelGeomDetUnit*>(stack->lowerDet());
  for (const_iterator cil = theLowerDetSet.begin(); cil != theLowerDetSet.end(); ++cil) {
    LogDebug("VectorHitBuilderAlgorithm") << " lower clusters " << std::endl;
    Phase2TrackerCluster1DRef cluL = edmNew::makeRefTo(clusters, cil);
#ifdef EDM_ML_DEBUG
    printCluster(stack->lowerDet(), &*cluL);
#endif
    auto&& lparamsLow = cpe_->localParameters(*cluL, *gduLow);
    upperIterator = 0;
    for (const_iterator ciu = theUpperDetSet.begin(); ciu != theUpperDetSet.end(); ++ciu) {
      LogDebug("VectorHitBuilderAlgorithm") << "\t upper clusters " << std::endl;
      Phase2TrackerCluster1DRef cluU = edmNew::makeRefTo(clusters, ciu);
#ifdef EDM_ML_DEBUG
      printCluster(stack->upperDet(), &*cluU);
#endif
      //applying the parallax correction
      double pC = computeParallaxCorrection(
          gduLow, lparamsLow.first, localGDUUpper[upperIterator], localParamsUpper[upperIterator].first);
      LogDebug("VectorHitBuilderAlgorithm") << " \t parallax correction:" << pC << std::endl;
      double lpos_upp_corr = 0.0;
      double lpos_low_corr = 0.0;
      auto const localUpperX = localParamsUpper[upperIterator].first.x();
      if (localUpperX > lparamsLow.first.x()) {
        if (localUpperX > 0) {
          lpos_low_corr = lparamsLow.first.x();
          lpos_upp_corr = localParamsUpper[upperIterator].first.x() - std::abs(pC);
        } else if (localUpperX < 0) {
          lpos_low_corr = lparamsLow.first.x() + std::abs(pC);
          lpos_upp_corr = localUpperX;
        }
      } else if (localUpperX < lparamsLow.first.x()) {
        if (localUpperX > 0) {
          lpos_low_corr = lparamsLow.first.x() - std::abs(pC);
          lpos_upp_corr = localUpperX;
        } else if (localUpperX < 0) {
          lpos_low_corr = lparamsLow.first.x();
          lpos_upp_corr = localUpperX + std::abs(pC);
        }
      } else {
        if (localUpperX > 0) {
          lpos_low_corr = lparamsLow.first.x();
          lpos_upp_corr = localUpperX - std::abs(pC);
        } else if (localUpperX < 0) {
          lpos_low_corr = lparamsLow.first.x();
          lpos_upp_corr = localUpperX + std::abs(pC);
        }
      }

      LogDebug("VectorHitBuilderAlgorithm") << " \t local pos upper corrected (x):" << lpos_upp_corr << std::endl;
      LogDebug("VectorHitBuilderAlgorithm") << " \t local pos lower corrected (x):" << lpos_low_corr << std::endl;

      double width = lpos_low_corr - lpos_upp_corr;
      LogDebug("VectorHitBuilderAlgorithm") << " \t width: " << width << std::endl;

      //old cut: indipendent from layer
      //building my tolerance : 10*sigma
      //double delta = 10.0 * sqrt(lparamsLow.second.xx() + localParamsUpper[upperIterator].second.xx());
      //LogDebug("VectorHitBuilderAlgorithm") << " \t delta: " << delta << std::endl;
      //if( (lpos_upp_corr < lpos_low_corr + delta) &&
      //    (lpos_upp_corr > lpos_low_corr - delta) ){
      //new cut: dependent on layers
      if (std::abs(width) < cut) {
        LogDebug("VectorHitBuilderAlgorithm") << " accepting VH! " << std::endl;
        VectorHit vh = buildVectorHit(stack, cluL, cluU);
        //protection: the VH can also be empty!!
        if (vh.isValid()) {
          vh_colAcc.push_back(vh);
        }

      } else {
        LogDebug("VectorHitBuilderAlgorithm") << " rejecting VH: " << std::endl;
        //storing vh rejected for combinatiorial studies
        VectorHit vh = buildVectorHit(stack, cluL, cluU);
        if (vh.isValid()) {
          vh_colRej.push_back(vh);
        }
      }
      upperIterator += 1;
    }
  }
}

checkClustersCompatibility()

bool VectorHitBuilderAlgorithm::checkClustersCompatibility	(	Local3DPoint &	posinner,
		Local3DPoint &	posouter,
		LocalError &	errinner,
		LocalError &	errouter
	)		const

Definition at line 70 of file VectorHitBuilderAlgorithm.cc.

Referenced by buildVectorHit().

                                                                                       {
  return true;
}

checkClustersCompatibilityBeforeBuilding()

bool VectorHitBuilderAlgorithm::checkClustersCompatibilityBeforeBuilding	(	edm::Handle< edmNew::DetSetVector< Phase2TrackerCluster1D >>	clusters,
		const Detset &	theLowerDetSet,
		const Detset &	theUpperDetSet
	)		const

Definition at line 49 of file VectorHitBuilderAlgorithm.cc.

References edmNew::DetSet< T >::begin(), bsc_activity_cfg::clusters, edmNew::DetSet< T >::end(), LogDebug, edmNew::makeRefTo(), and edmNew::DetSet< T >::size().

Referenced by buildVectorHits().

                                        {
  if (theLowerDetSet.size() == 1 && theUpperDetSet.size() == 1)
    return true;

  //order lower clusters in u
  std::vector<Phase2TrackerCluster1D> lowerClusters;
  lowerClusters.reserve(theLowerDetSet.size());
  if (theLowerDetSet.size() > 1)
    LogDebug("VectorHitBuilderAlgorithm") << " more than 1 lower cluster! " << std::endl;
  if (theUpperDetSet.size() > 1)
    LogDebug("VectorHitBuilderAlgorithm") << " more than 1 upper cluster! " << std::endl;
  for (const_iterator cil = theLowerDetSet.begin(); cil != theLowerDetSet.end(); ++cil) {
    Phase2TrackerCluster1DRef clusterLower = edmNew::makeRefTo(clusters, cil);
    lowerClusters.push_back(*clusterLower);
  }
  return true;
}

curvatureANDphi()

VectorHitBuilderAlgorithm::CurvatureAndPhi VectorHitBuilderAlgorithm::curvatureANDphi	(	Global3DPoint	gPositionLower,
		Global3DPoint	gPositionUpper,
		GlobalError	gErrorLower,
		GlobalError	gErrorUpper
	)		const

Definition at line 366 of file VectorHitBuilderAlgorithm.cc.

References PixelRecoUtilities::curvature(), GlobalErrorBase< T, ErrorWeightType >::cxx(), GlobalErrorBase< T, ErrorWeightType >::cyx(), GlobalErrorBase< T, ErrorWeightType >::cyy(), SiPixelRawToDigiRegional_cfi::deltaPhi, mps_fire::i, M_PI, PV3DBase< T, PVType, FrameType >::perp2(), phi, conifer::pow(), mps_fire::result, mathSSE::sqrt(), groupFilesInBlocks::temp, PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

Referenced by buildVectorHit().

                                                                                                                     {
  VectorHitBuilderAlgorithm::CurvatureAndPhi result;

  float curvature = -999.;
  float errorCurvature = -999.;
  float phi = -999.;

  float h1 = gPositionLower.x() * gPositionUpper.y() - gPositionUpper.x() * gPositionLower.y();

  //determine sign of curvature
  AlgebraicVector2 n1;
  n1[0] = -gPositionLower.y();
  n1[1] = gPositionLower.x();
  AlgebraicVector2 n2;
  n2[0] = gPositionUpper.x() - gPositionLower.x();
  n2[1] = gPositionUpper.y() - gPositionLower.y();

  double n3 = n1[0] * n2[0] + n1[1] * n2[1];
  double signCurv = -copysign(1.0, n3);
  double phi1 = atan2(gPositionUpper.y() - gPositionLower.y(), gPositionUpper.x() - gPositionLower.x());

  double x2Low = pow(gPositionLower.x(), 2);
  double y2Low = pow(gPositionLower.y(), 2);
  double x2Up = pow(gPositionUpper.x(), 2);
  double y2Up = pow(gPositionUpper.y(), 2);

  if (h1 != 0) {
    double h2 = 2 * h1;
    double h2Inf = 1. / (2 * h1);
    double r12 = gPositionLower.perp2();
    double r22 = gPositionUpper.perp2();
    double h3 = pow(n2[0], 2) + pow(n2[1], 2);
    double h4 = -x2Low * gPositionUpper.x() + gPositionLower.x() * x2Up + gPositionLower.x() * y2Up -
                gPositionUpper.x() * y2Low;
    double h5 =
        x2Low * gPositionUpper.y() - x2Up * gPositionLower.y() + y2Low * gPositionUpper.y() - gPositionLower.y() * y2Up;

    //radius of circle
    double invRho2 = (4. * h1 * h1) / (r12 * r22 * h3);
    curvature = sqrt(invRho2);

    //center of circle
    double xcentre = h5 / h2;
    double ycentre = h4 / h2;

    //to compute phi at the cluster points
    double xtg = gPositionLower.y() - ycentre;
    double ytg = -(gPositionLower.x() - xcentre);

    //to compute phi at the origin
    phi = atan2(ytg, xtg);

    AlgebraicROOTObject<4, 4>::Matrix jacobian;

    double denom1 = 1. / sqrt(r12 * r22 * h3);
    double denom2 = 1. / (pow(r12 * r22 * h3, 1.5));
    jacobian[0][0] = 1.0;  // dx1/dx1 dx1/dy1 dx2/dx1 dy2/dx1
    jacobian[1][1] = 1.0;  //dy1/dx1 dy1/dy1 dy2/dx1 dy2/dx1
    jacobian[2][0] =
        -2. * ((h1 * (gPositionLower.x() * r22 * h3 + (gPositionLower.x() - gPositionUpper.x()) * r12 * r22)) * denom2 -
               (gPositionUpper.y()) * denom1);  // dkappa/dx1
    jacobian[2][1] =
        -2. * ((gPositionUpper.x()) * denom1 +
               (h1 * (gPositionLower.y() * r22 * h3 + r12 * r22 * (gPositionLower.y() - gPositionUpper.y()))) *
                   denom2);  // dkappa/dy1
    jacobian[2][2] =
        -2. * ((gPositionLower.y()) * denom1 +
               (h1 * (gPositionUpper.x() * r12 * h3 - (gPositionLower.x() - gPositionUpper.x()) * r12 * r22)) *
                   denom2);  // dkappa/dx2
    jacobian[2][3] =
        -2. * ((h1 * (gPositionUpper.y() * r12 * h3 - r12 * r22 * (gPositionLower.y() - gPositionUpper.y()))) * denom2 -
               (gPositionLower.x()) * denom1);  // dkappa/dy2
    AlgebraicVector2 mVector;
    //to compute phi at the cluster points
    mVector[0] = (gPositionLower.y() - ycentre) * invRho2;   // dphi/dxcentre
    mVector[1] = -(gPositionLower.x() - xcentre) * invRho2;  // dphi/dycentre
    //to compute phi at the origin

    double h22Inv = 1. / pow(h2, 2);

    AlgebraicROOTObject<2, 4>::Matrix kMatrix;
    kMatrix[0][0] =
        2. * ((gPositionLower.x() * gPositionUpper.y()) * h2Inf - (gPositionUpper.y() * h5) * h22Inv);  // dxm/dx1
    kMatrix[0][1] = (2. * gPositionUpper.x() * h5) * h22Inv -
                    (x2Up + y2Up - 2. * gPositionLower.y() * gPositionUpper.y()) * h2Inf;  // dxm/dy1
    kMatrix[0][2] =
        2. * ((gPositionLower.y() * h5) * h22Inv - (gPositionUpper.x() * gPositionLower.y()) * h2Inf);  // dxm/dx2
    kMatrix[0][3] = (x2Low + y2Low - 2. * gPositionUpper.y() * gPositionLower.y()) * h2Inf -
                    (2. * gPositionLower.x() * h5) * h22Inv;  // dxm/dy2
    kMatrix[1][0] = (x2Up - 2. * gPositionLower.x() * gPositionUpper.x() + y2Up) * h2Inf -
                    (2. * gPositionUpper.y() * h4) * h22Inv;  // dym/dx1
    kMatrix[1][1] =
        2. * ((gPositionUpper.x() * h4) * h22Inv - (gPositionUpper.x() * gPositionLower.y()) * h2Inf);  // dym/dy1
    kMatrix[1][2] = (2. * gPositionLower.y() * h4) * h22Inv -
                    (x2Low - 2. * gPositionUpper.x() * gPositionLower.x() + y2Low) * h2Inf;  // dym/dx2
    kMatrix[1][3] =
        2. * (gPositionLower.x() * gPositionUpper.y()) * h2Inf - (gPositionLower.x() * h4) * h22Inv;  // dym/dy2

    AlgebraicVector4 nMatrix = mVector * kMatrix;
    jacobian[3][0] = nMatrix[0];  // dphi/(dx1,dy1,dx2,dy2)
    jacobian[3][1] = nMatrix[1];  // dphi/(dx1,dy1,dx2,dy2)
    jacobian[3][2] = nMatrix[2];  // dphi/(dx1,dy1,dx2,dy2)
    jacobian[3][3] = nMatrix[3];  // dphi/(dx1,dy1,dx2,dy2)

    //assign correct sign to the curvature errors
    if ((signCurv < 0 && curvature > 0) || (signCurv > 0 && curvature < 0)) {
      curvature = -curvature;
      for (int i = 0; i < 4; i++) {
        jacobian[2][i] = -jacobian[2][i];
      }
    }

    // bring phi in the same quadrant as phi1
    if (deltaPhi(phi, phi1) > M_PI / 2.) {
      phi = phi + M_PI;
      if (phi > M_PI)
        phi = phi - 2. * M_PI;
    }

    //computing the curvature error
    AlgebraicVector4 curvatureJacobian;
    for (int i = 0; i < 4; i++) {
      curvatureJacobian[i] = jacobian[2][i];
    }

    AlgebraicROOTObject<4, 4>::Matrix gErrors;

    gErrors[0][0] = gErrorLower.cxx();
    gErrors[0][1] = gErrorLower.cyx();
    gErrors[1][0] = gErrorLower.cyx();
    gErrors[1][1] = gErrorLower.cyy();
    gErrors[2][2] = gErrorUpper.cxx();
    gErrors[2][3] = gErrorUpper.cyx();
    gErrors[3][2] = gErrorUpper.cyx();
    gErrors[3][3] = gErrorUpper.cyy();

    AlgebraicVector4 temp = curvatureJacobian;
    temp = temp * gErrors;
    errorCurvature = temp[0] * curvatureJacobian[0] + temp[1] * curvatureJacobian[1] + temp[2] * curvatureJacobian[2] +
                     temp[3] * curvatureJacobian[3];
  }

  result.curvature = curvature;
  result.curvatureError = errorCurvature;
  result.phi = phi;
  return result;
}

fit()

void VectorHitBuilderAlgorithm::fit	(	float	x[2],
		float	y[2],
		float	sigy[2],
		Local3DPoint &	pos,
		Local3DVector &	dir,
		AlgebraicSymMatrix22 &	covMatrix,
		double &	chi2
	)		const

Definition at line 335 of file VectorHitBuilderAlgorithm.cc.

References nano_mu_local_reco_cff::chi2, DeadROC_duringRun::dir, PVValHelper::dy, dqmiolumiharvest::j, linearFit(), slope, mathSSE::sqrt(), x, and y.

Referenced by fit2Dzx(), fit2Dzy(), and trackingPlots.Iteration::modules().

                                                        {
  float slope = 0.;
  float intercept = 0.;
  float covss = 0.;
  float covii = 0.;
  float covsi = 0.;

  linearFit(x, y, 2, sigy2, slope, intercept, covss, covii, covsi);

  covMatrix[0][0] = covss;  // this is var(dy/dz)
  covMatrix[1][1] = covii;  // this is var(y)
  covMatrix[1][0] = covsi;  // this is cov(dy/dz,y)

  for (unsigned int j = 0; j < 2; j++) {
    const double ypred = intercept + slope * x[j];
    const double dy = (y[j] - ypred) / sqrt(sigy2[j]);
    chi2 += dy * dy;
  }

  pos = Local3DPoint(intercept, 0., 0.);
  //difference in z is the difference of the lowermost and the uppermost cluster z pos
  float slopeZ = x[1] - x[0];
  dir = LocalVector(slope, 0., slopeZ);
}

fit2Dzx()

void VectorHitBuilderAlgorithm::fit2Dzx	(	const Local3DPoint	lpCI,
		const Local3DPoint	lpCO,
		const LocalError	leCI,
		const LocalError	leCO,
		Local3DPoint &	pos,
		Local3DVector &	dir,
		AlgebraicSymMatrix22 &	covMatrix,
		double &	chi2
	)		const

Definition at line 297 of file VectorHitBuilderAlgorithm.cc.

References nano_mu_local_reco_cff::chi2, DeadROC_duringRun::dir, fit(), mathSSE::sqrt(), x, PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), y, and PV3DBase< T, PVType, FrameType >::z().

Referenced by buildVectorHit().

                                                            {
  float x[2] = {lpCI.z(), lpCO.z()};
  float y[2] = {lpCI.x(), lpCO.x()};
  float sqCI = sqrt(leCI.xx());
  float sqCO = sqrt(leCO.xx());
  float sigy2[2] = {sqCI * sqCI, sqCO * sqCO};

  fit(x, y, sigy2, pos, dir, covMatrix, chi2);

  return;
}

fit2Dzy()

void VectorHitBuilderAlgorithm::fit2Dzy	(	const Local3DPoint	lpCI,
		const Local3DPoint	lpCO,
		const LocalError	leCI,
		const LocalError	leCO,
		Local3DPoint &	pos,
		Local3DVector &	dir,
		AlgebraicSymMatrix22 &	covMatrix,
		double &	chi2
	)		const

Definition at line 316 of file VectorHitBuilderAlgorithm.cc.

References nano_mu_local_reco_cff::chi2, DeadROC_duringRun::dir, fit(), mathSSE::sqrt(), x, y, PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by buildVectorHit().

                                                            {
  float x[2] = {lpCI.z(), lpCO.z()};
  float y[2] = {lpCI.y(), lpCO.y()};
  float sqCI = sqrt(leCI.yy());
  float sqCO = sqrt(leCO.yy());
  float sigy2[2] = {sqCI * sqCI, sqCO * sqCO};

  fit(x, y, sigy2, pos, dir, covMatrix, chi2);

  return;
}

run()

void VectorHitBuilderAlgorithm::run ( edm::Handle< edmNew::DetSetVector< Phase2TrackerCluster1D >>  clusters, VectorHitCollection &  vhAcc, VectorHitCollection &  vhRej, edmNew::DetSetVector< Phase2TrackerCluster1D > &  clustersAcc, edmNew::DetSetVector< Phase2TrackerCluster1D > &  clustersRej  ) const

overridevirtual

Implements VectorHitBuilderAlgorithmBase.

Definition at line 8 of file VectorHitBuilderAlgorithm.cc.

References buildVectorHits(), bsc_activity_cfg::clusters, TrackerGeometry::idToDet(), TrackerTopology::isLower(), LogDebug, TrackerTopology::partnerDetId(), DetId::rawId(), TrackerTopology::stack(), VectorHitBuilderAlgorithmBase::tkGeom_, and VectorHitBuilderAlgorithmBase::tkTopo_.

Referenced by VectorHitBuilderEDProducer::run().

                                                                                                   {
  LogDebug("VectorHitBuilderAlgorithm") << "Run VectorHitBuilderAlgorithm ... \n";
  const auto* clustersPhase2Collection = clusters.product();

  //loop over the DetSetVector
  LogDebug("VectorHitBuilderAlgorithm") << "with #clusters : " << clustersPhase2Collection->size() << std::endl;
  for (auto dSViter : *clustersPhase2Collection) {
    unsigned int rawDetId1(dSViter.detId());
    DetId detId1(rawDetId1);
    DetId lowerDetId, upperDetId;
    if (tkTopo_->isLower(detId1)) {
      lowerDetId = detId1;
      upperDetId = tkTopo_->partnerDetId(detId1);
    } else
      continue;

    DetId detIdStack = tkTopo_->stack(detId1);

    //debug
    LogDebug("VectorHitBuilderAlgorithm") << "  DetId stack : " << detIdStack.rawId() << std::endl;
    LogDebug("VectorHitBuilderAlgorithm") << "  DetId lower set of clusters  : " << lowerDetId.rawId();
    LogDebug("VectorHitBuilderAlgorithm") << "  DetId upper set of clusters  : " << upperDetId.rawId() << std::endl;

    const GeomDet* gd;
    const StackGeomDet* stackDet;
    const auto& it_detLower = dSViter;
    const auto& it_detUpper = clustersPhase2Collection->find(upperDetId);

    if (it_detUpper != clustersPhase2Collection->end()) {
      gd = tkGeom_->idToDet(detIdStack);
      stackDet = dynamic_cast<const StackGeomDet*>(gd);
      buildVectorHits(vhAcc, vhRej, detIdStack, stackDet, clusters, it_detLower, *it_detUpper);
    }
  }
  LogDebug("VectorHitBuilderAlgorithm") << "End run VectorHitBuilderAlgorithm ... \n";
}