Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Switch FastCircleFit to use Eigen, generalize FastCircleFit and RZLine interfaces #15260

Merged
merged 7 commits into from Jul 29, 2016
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions CommonTools/Utils/interface/DynArray.h
Expand Up @@ -44,6 +44,7 @@ public :
unsigned int size() const { return s;}
bool empty() const { return 0==s;}

T const * data() const { return a; }
T const &front() const { return a[0];}
T const & back() const { return a[s-1];}

Expand Down
7 changes: 2 additions & 5 deletions RecoPixelVertexing/PixelLowPtUtilities/src/TrackFitter.cc
Expand Up @@ -28,7 +28,7 @@

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

#include "RecoPixelVertexing/PixelTrackFitting/src/RZLine.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h"
#include "RecoPixelVertexing/PixelTrackFitting/src/CircleFromThreePoints.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackBuilder.h"

Expand Down Expand Up @@ -134,10 +134,7 @@ reco::Track* TrackFitter::run
if(nhits > 2)
{
RZLine rzLine(points,errors,isBarrel);
float cotTheta, intercept, covss, covii, covsi;
rzLine.fit(cotTheta, intercept, covss, covii, covsi);
chi2 = rzLine.chi2(cotTheta, intercept);
//FIXME: check which intercept to use!
chi2 = rzLine.chi2();
}

// build pixel track
Expand Down
120 changes: 120 additions & 0 deletions RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h
@@ -0,0 +1,120 @@
#ifndef PixelTrackFitting_RZLine_H
#define PixelTrackFitting_RZLine_H

#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
#include "CommonTools/Utils/interface/DynArray.h"

#include "CommonTools/Statistics/interface/LinearFit.h"

#include <vector>

class RZLine {
public:
struct ErrZ2_tag {};

/**
* Constructor for containers of GlobalPoint, GlobalError, and bool
*
* @tparam P Container of GlobalPoint
* @tparam E Container of GlobalError
* @tparam B Container of bool
*
* Container can be e.g. std::vector, std::array, or DynArray.
*
* Although for std::array use this constructor could be specialized
* to use std::array instead of DynArray for temporary storage.
*/
template <typename P, typename E, typename B>
RZLine(const P& points, const E& errors, const B& isBarrel) {
const size_t n = points.size();
declareDynArray(float, n, r);
declareDynArray(float, n, z);
declareDynArray(float, n, errZ2);
for(size_t i=0; i<n; ++i) {
const GlobalPoint& p = points[i];
r[i] = p.perp();
z[i] = p.z();
}

float simpleCot2 = sqr( (z[n-1]-z[0])/ (r[n-1]-r[0]) );
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a division by zero possibility here if (r[n-1]-r[0]) == 0? Or would that be prevented from happening? Or would it take too long to check for the zero value?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Note that this code is exactly as before) I guess it could, in principle, happen that the first and last hit would have exactly the same r, if all hits come from FPix/TID/TEC. But I'd expect this situation to be very unlikely, especially given the constraints of seeds to point towards beamspot (even if loosely on strip-triplet steps).

for (size_t i=0; i<n; ++i) {
errZ2[i] = (isBarrel[i]) ? errors[i].czz() :
errors[i].rerr(points[i]) * simpleCot2;
}

calculate(r, z, errZ2);
}

/**
* Constructor for std::vector of r, z, and z standard deviation
*/
RZLine(const std::vector<float> & r,
const std::vector<float> & z,
const std::vector<float> & errZ) {
const size_t n = errZ.size();
declareDynArray(float, n, errZ2);
for(size_t i=0; i<n; ++i) errZ2[i] = sqr(errZ[i]);
calculate(r, z, errZ2);
}

/**
* Constructor for std::array of r, z, and z standard deviation
*/
template <size_t N>
RZLine(const std::array<float, N>& r,
const std::array<float, N>& z,
const std::array<float, N>& errZ) {
std::array<float, N> errZ2;
for(size_t i=0; i<N; ++i) errZ2[i] = sqr(errZ[i]);
calculate(r, z, errZ2);
}

/**
* Constructor for container of r, z, and z variance
*
* @tparam T Container of float
*
* Container can be e.g. std::vector, std::array, or DynArray.
*
* The ErrZ2_tag parameter is used to distinguish this constructor
* from other 3-parameter constructors.
*
* Passing variance is useful in cases where it is already available
* to avoid making a square of a square root.
*/
template <typename T>
RZLine (const T& r, const T& z, const T& errZ2, ErrZ2_tag) {
calculate(r, z, errZ2);
}

float cotTheta() const { return cotTheta_; }
float intercept() const { return intercept_; }
float covss() const { return covss_; }
float covii() const { return covii_; }
float covsi() const { return covsi_; }

float chi2() const { return chi2_; }

private:
template <typename R, typename Z, typename E>
void calculate(const R& r, const Z& z, const E& errZ2) {
const size_t n = r.size();
linearFit(r.data(), z.data(), n, errZ2.data(), cotTheta_, intercept_, covss_, covii_, covsi_);
chi2_ = 0.f;
for(size_t i=0; i<n; ++i) {
chi2_ += sqr( ((z[i]-intercept_) - cotTheta_*r[i]) ) / errZ2[i];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How safe is this division? Could unusual circumstances cause errZ2 to contain some zero values? Is it worth some possible performance cost to protect against division by zero?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Note that this code is exactly as before) In practice errZ2 can be zero only if some TrackingRecHit has zero czz or zero rerr. I guess those would be a sign of something going wrong elsewhere.

}
}

template <typename T>
T sqr(T t) { return t*t; }

float cotTheta_;
float intercept_;
float covss_;
float covii_;
float covsi_;
float chi2_;
};
#endif
Expand Up @@ -27,7 +27,7 @@

#include "ConformalMappingFit.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackBuilder.h"
#include "RZLine.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h"

#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
#include "RecoTracker/TkMSParametrization/interface/LongitudinalBendingCorrection.h"
Expand Down Expand Up @@ -123,15 +123,13 @@ reco::Track* PixelFitterByConformalMappingAndLine::run(
// line fit (R-Z plane)
//
RZLine rzLine(r,z,errZ);
float cottheta, intercept, covss, covii, covsi;
rzLine.fit(cottheta, intercept, covss, covii, covsi);

//
// parameters for track builder
//
Measurement1D zip(intercept, sqrt(covii));
Measurement1D cotTheta(cottheta, sqrt(covss));
float chi2 = parabola.chi2() + rzLine.chi2(cottheta, intercept);
Measurement1D zip(rzLine.intercept(), sqrt(rzLine.covii()));
Measurement1D cotTheta(rzLine.cotTheta(), sqrt(rzLine.covss()));
float chi2 = parabola.chi2() + rzLine.chi2();
int charge = parabola.charge();


Expand Down
Expand Up @@ -26,7 +26,7 @@

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

#include "RZLine.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h"
#include "CircleFromThreePoints.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackBuilder.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackErrorParam.h"
Expand Down Expand Up @@ -161,9 +161,7 @@ reco::Track* PixelFitterByHelixProjections::run(
float chi2 = 0;
if (nhits > 2) {
RZLine rzLine(points,errors,isBarrel);
float cottheta, intercept, covss, covii, covsi;
rzLine.fit(cottheta, intercept, covss, covii, covsi);
chi2 = rzLine.chi2(cottheta, intercept); //FIXME: check which intercept to use!
chi2 = rzLine.chi2();
}

PixelTrackBuilder builder;
Expand Down
55 changes: 0 additions & 55 deletions RecoPixelVertexing/PixelTrackFitting/src/RZLine.cc

This file was deleted.

33 changes: 0 additions & 33 deletions RecoPixelVertexing/PixelTrackFitting/src/RZLine.h

This file was deleted.

Expand Up @@ -109,11 +109,12 @@ CAHitQuadrupletGenerator::findQuadruplets (const TrackingRegion& region, Ordered
const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(es);

// re-used thoughout, need to be vectors because of RZLine interface
std::vector<float> bc_r(4), bc_z(4), bc_errZ(4);

declareDynArray(GlobalPoint, 4, gps);
declareDynArray(GlobalError, 4, ges);
declareDynArray(bool, 4, barrels);
std::array<float, 4> bc_r;
std::array<float, 4> bc_z;
std::array<float, 4> bc_errZ2;
std::array<GlobalPoint, 4> gps;
std::array<GlobalError, 4> ges;
std::array<bool, 4> barrels;

unsigned int numberOfFoundQuadruplets = foundQuadruplets.size();

Expand Down Expand Up @@ -169,18 +170,14 @@ CAHitQuadrupletGenerator::findQuadruplets (const TrackingRegion& region, Ordered
bc_r[i] = sqrt( sqr(point.x() - region.origin().x()) + sqr(point.y() - region.origin().y()) );
bc_r[i] += pixelrecoutilities::LongitudinalBendingCorrection(pt, es)(bc_r[i]);
bc_z[i] = point.z() - region.origin().z();
bc_errZ[i] = (barrels[i]) ? sqrt(error.czz()) : sqrt( error.rerr(point) ) * simpleCot;
bc_errZ2[i] = (barrels[i]) ? error.czz() : error.rerr(point)*sqr(simpleCot);
}
RZLine rzLine(bc_r, bc_z, bc_errZ);
float cottheta, intercept, covss, covii, covsi;
rzLine.fit(cottheta, intercept, covss, covii, covsi);
chi2 = rzLine.chi2(cottheta, intercept);
RZLine rzLine(bc_r, bc_z, bc_errZ2, RZLine::ErrZ2_tag());
chi2 = rzLine.chi2();
} else
{
RZLine rzLine(gps, ges, barrels);
float cottheta, intercept, covss, covii, covsi;
rzLine.fit(cottheta, intercept, covss, covii, covsi);
chi2 = rzLine.chi2(cottheta, intercept);
chi2 = rzLine.chi2();
}
if (edm::isNotFinite(chi2) || chi2 > thisMaxChi2)
{
Expand Down
Expand Up @@ -4,7 +4,7 @@
#include "RecoPixelVertexing/PixelTriplets/interface/HitQuadrupletGenerator.h"
#include "RecoTracker/TkSeedingLayers/interface/SeedComparitorFactory.h"
#include "RecoTracker/TkSeedingLayers/interface/SeedComparitor.h"
#include "RecoPixelVertexing/PixelTrackFitting/src/RZLine.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h"
#include "RecoTracker/TkSeedGenerator/interface/FastCircleFit.h"
#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
#include "RecoTracker/TkMSParametrization/interface/LongitudinalBendingCorrection.h"
Expand Down
Expand Up @@ -7,7 +7,7 @@
#include "RecoPixelVertexing/PixelTriplets/plugins/KDTreeLinkerTools.h"
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"

#include "RecoPixelVertexing/PixelTrackFitting/src/RZLine.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h"
#include "RecoTracker/TkSeedGenerator/interface/FastCircleFit.h"
#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
#include "RecoTracker/TkMSParametrization/interface/LongitudinalBendingCorrection.h"
Expand Down Expand Up @@ -113,8 +113,13 @@ void PixelQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, Orde
const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(es);
const QuantityDependsPtEval extraPhiToleranceEval = extraPhiTolerance.evaluator(es);

// re-used thoughout, need to be vectors because of RZLine interface
std::vector<float> bc_r(4), bc_z(4), bc_errZ(4);
// re-used thoughout
std::array<float, 4> bc_r;
std::array<float, 4> bc_z;
std::array<float, 4> bc_errZ2;
std::array<GlobalPoint, 4> gps;
std::array<GlobalError, 4> ges;
std::array<bool, 4> barrels;

// Loop over triplets
for(const auto& triplet: triplets) {
Expand All @@ -137,9 +142,6 @@ void PixelQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, Orde
return id == PixelSubdetector::PixelBarrel;
};

declareDynArray(GlobalPoint,4, gps);
declareDynArray(GlobalError,4, ges);
declareDynArray(bool,4, barrels);
gps[0] = triplet.inner()->globalPosition();
ges[0] = triplet.inner()->globalPositionError();
barrels[0] = isBarrel(triplet.inner()->geographicalId().subdetId());
Expand Down Expand Up @@ -220,18 +222,14 @@ void PixelQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, Orde
bc_r[i] = sqrt( sqr(point.x()-region.origin().x()) + sqr(point.y()-region.origin().y()) );
bc_r[i] += pixelrecoutilities::LongitudinalBendingCorrection(pt,es)(bc_r[i]);
bc_z[i] = point.z()-region.origin().z();
bc_errZ[i] = (barrels[i]) ? sqrt(error.czz()) : sqrt( error.rerr(point) )*simpleCot;
bc_errZ2[i] = (barrels[i]) ? error.czz() : error.rerr(point)*sqr(simpleCot);
}
RZLine rzLine(bc_r,bc_z,bc_errZ);
float cottheta, intercept, covss, covii, covsi;
rzLine.fit(cottheta, intercept, covss, covii, covsi);
chi2 = rzLine.chi2(cottheta, intercept);
RZLine rzLine(bc_r,bc_z,bc_errZ2, RZLine::ErrZ2_tag());
chi2 = rzLine.chi2();
}
else {
RZLine rzLine(gps, ges, barrels);
float cottheta, intercept, covss, covii, covsi;
rzLine.fit(cottheta, intercept, covss, covii, covsi);
chi2 = rzLine.chi2(cottheta, intercept);
chi2 = rzLine.chi2();
}
if(edm::isNotFinite(chi2) || chi2 > thisMaxChi2) {
continue;
Expand Down