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.

Sign up for GitHub

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

Jump to bottom

Avoid Nans in primary vertex #5490

Merged
merged 7 commits into from Sep 26, 2014
Merged

Avoid Nans in primary vertex #5490

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
c975167
robistify DA
VinInn Sep 20, 2014
4371659
compute pt once
VinInn Sep 21, 2014
0dd30dc
optimize sorting, add warning
VinInn Sep 21, 2014
d233608
avoid to trash memory
VinInn Sep 21, 2014
9dae5d1
clean config
VinInn Sep 23, 2014
44dec68
remove g
VinInn Sep 23, 2014
23ba733
avoid trigonometry
VinInn Sep 23, 2014
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
90 changes: 40 additions & 50 deletions RecoVertex/AdaptiveVertexFit/src/AdaptiveVertexFitter.cc
View file
Open in desktop
Expand Up @@ -18,12 +18,16 @@ using namespace edm;
using namespace std;

namespace {
struct CompareTwoTracks {
int operator() ( const reco::TransientTrack & a, const reco::TransientTrack & b ) {
return ( a.impactPointState().globalMomentum().perp() >
b.impactPointState().globalMomentum().perp() ) ;
};
};
void sortTracksByPt(std::vector<reco::TransientTrack> & cont) {
auto s = cont.size();
float pt2[s]; int ind[s]; int i=0;
for (auto const & tk : cont) { ind[i]=i; pt2[++i] = tk.impactPointState().globalMomentum().perp2();}
Copy link
Contributor

Choose a reason for hiding this comment

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

this should have been i++, or, more straightforward and less error prone, pt2[i] = ...; ++i;
Valgrind reports errors here.
I suspect this is where our vertex/btag discrepancies originate.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

oops, is this in 72?

Copy link
Contributor

Choose a reason for hiding this comment

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

@davidlt @ktf
Did this get picked up by ASAN (I'm guessing we are running it) or by valgrind (which I guess doesn't run regularly) ?

Copy link
Contributor

Choose a reason for hiding this comment

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

I think this is only in 73X

Copy link
Contributor

Choose a reason for hiding this comment

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

@slava77 currently I don't run ASAN builds. I discussed with @ktf recently of doing so in some future. Yet, I might do it manually in some days if you want.

std::sort(ind,ind+s, [&](int i, int j){return pt2[i]>pt2[j];} );
std::vector<reco::TransientTrack> tmp; tmp.reserve(s);
for (auto i=0U; i<s; ++i) tmp.emplace_back(std::move( cont[ind[i]] ) );
cont.swap(tmp);
}

// typedef ReferenceCountingPointer<VertexTrack<5> > RefCountedVertexTrack;
typedef AdaptiveVertexFitter::RefCountedVertexTrack RefCountedVertexTrack;

Expand All @@ -37,11 +41,8 @@ namespace {
ret(2,2)=initialError;
return ret;
}
GlobalError fitError()
{
static const GlobalError err( initFitError() );
return err;
}

GlobalError const fitError = initFitError();

AlgebraicSymMatrix33 initLinePointError() {
// that's how we model the error of the linearization point.
Expand All @@ -51,28 +52,24 @@ namespace {
// ret(0,0)=1e-7; ret(1,1)=1e-7; ret(2,2)=1e-7;
return ret;
}
GlobalError linPointError()
{
static const GlobalError err( initLinePointError() );
return err;
}

struct DistanceToRefPoint {
DistanceToRefPoint ( const GlobalPoint & ref ) : theRef ( ref ) {}
GlobalError const linPointError = initLinePointError();

bool operator() ( const RefCountedVertexTrack & v1, const RefCountedVertexTrack & v2 )
{
return ( distance ( v1 ) < distance ( v2 ) );
}

float distance ( const RefCountedVertexTrack & v1 )
{
return ( v1->linearizedTrack()->track().initialFreeState().position() - theRef ).mag2();
}

private:
GlobalPoint theRef;
};

void
sortByDistanceToRefPoint (std::vector<RefCountedVertexTrack> & cont, const GlobalPoint ref ) {
auto s = cont.size();
float d2[s]; int ind[s]; int i=0;
for (auto const & tk : cont) { ind[i]=i; d2[++i] = (tk->linearizedTrack()->track().initialFreeState().position() - ref ).mag2();}
std::sort(ind,ind+s, [&](int i, int j){return d2[i]<d2[j];} );
std::vector<RefCountedVertexTrack> tmp; tmp.reserve(s);
for (auto i=0U; i<s; ++i) tmp.emplace_back(std::move( cont[ind[i]] ) );
cont.swap(tmp);
}



#ifdef STORE_WEIGHTS
//NOTE: This is not thread safe
Expand Down Expand Up @@ -163,14 +160,14 @@ AdaptiveVertexFitter::vertex(const vector<reco::TransientTrack> & unstracks) con
return CachingVertex<5>(); // return invalid vertex
};
vector < reco::TransientTrack > tracks = unstracks;
sort ( tracks.begin(), tracks.end(), CompareTwoTracks() );
sortTracksByPt( tracks);
// Linearization Point
GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
// Initial vertex seed, with a very large error matrix
VertexState lseed (linP, linPointError() );
VertexState lseed (linP, linPointError );
vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, lseed);

VertexState seed (linP, fitError() );
VertexState seed (linP, fitError );
return fit(vtContainer, seed, false);
}

Expand All @@ -185,7 +182,7 @@ AdaptiveVertexFitter::vertex(const vector<RefCountedVertexTrack> & tracks) const
};
// Initial vertex seed, with a very small weight matrix
GlobalPoint linP = tracks[0]->linearizedTrack()->linearizationPoint();
VertexState seed (linP, fitError() );
VertexState seed (linP, fitError );
return fit(tracks, seed, false);
}

Expand Down Expand Up @@ -218,9 +215,9 @@ AdaptiveVertexFitter::vertex(const vector<reco::TransientTrack> & tracks,
return CachingVertex<5>(); // return invalid vertex
};
// Initial vertex seed, with a very large error matrix
VertexState seed (linPoint, linPointError() );
VertexState seed (linPoint, linPointError );
vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, seed);
VertexState fitseed (linPoint, fitError() );
VertexState fitseed (linPoint, fitError );
return fit(vtContainer, fitseed, false);
}

Expand All @@ -244,12 +241,12 @@ AdaptiveVertexFitter::vertex(const vector<reco::TransientTrack> & unstracks,
vector<RefCountedVertexTrack> vtContainer;

vector < reco::TransientTrack > tracks = unstracks;
sort ( tracks.begin(), tracks.end(), CompareTwoTracks() );
sortTracksByPt(tracks);

if (tracks.size() > 1) {
// Linearization Point search if there are more than 1 track
GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
VertexState lpState(linP, linPointError() );
VertexState lpState(linP, linPointError );
vtContainer = linearizeTracks(tracks, lpState);
} else {
// otherwise take the beamspot position.
Expand Down Expand Up @@ -437,8 +434,7 @@ AdaptiveVertexFitter::reWeightTracks(

finalTracks.push_back(vTrData);
}
sort ( finalTracks.begin(), finalTracks.end(),
DistanceToRefPoint ( vertex.position() ) );
sortByDistanceToRefPoint( finalTracks, vertex.position() );
// cout << "[AdaptiveVertexFitter] /now reweight" << endl;
return finalTracks;
}
Expand Down Expand Up @@ -534,12 +530,9 @@ AdaptiveVertexFitter::fit( const vector<RefCountedVertexTrack> & tracks,
priorVertexPosition,priorVertexError,initialTracks,0);
}

// vector<RefCountedVertexTrack> globalVTracks = tracks;
std::vector<RefCountedVertexTrack> globalVTracks = tracks;
// sort the tracks, according to distance to seed!
vector<RefCountedVertexTrack> globalVTracks ( tracks.size() );

partial_sort_copy ( tracks.begin(), tracks.end(),
globalVTracks.begin(), globalVTracks.end(), DistanceToRefPoint ( priorSeed.position() ) );
sortByDistanceToRefPoint (globalVTracks, priorSeed.position() );

// main loop through all the VTracks
int lpStep = 0; int step = 0;
Expand Down Expand Up @@ -586,11 +579,8 @@ AdaptiveVertexFitter::fit( const vector<RefCountedVertexTrack> & tracks,
if ((**i).weight() > 0.) nVertex = theUpdator->add( fVertex, *i );
else nVertex = fVertex;
if (nVertex.isValid()) {
if ( (**i).weight() >= theWeightThreshold )
{
ns_trks++;
};

if ( (**i).weight() >= theWeightThreshold ) ns_trks++;

if ( fabs ( nVertex.position().z() ) > 10000. ||
Copy link
Contributor

Choose a reason for hiding this comment

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

probably it's a silly point,
but could we get rid of such hard-coded numbers ?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

If somebody volunteers (the code was there already)

nVertex.position().perp()>120.)
{
Expand All @@ -612,7 +602,7 @@ AdaptiveVertexFitter::fit( const vector<RefCountedVertexTrack> & tracks,
<< "The updator returned an invalid vertex when adding track "
<< i-globalVTracks.begin()
<< ".\n Your vertex might just have lost one good track.";
};
}
}
previousPosition = newPosition;
newPosition = fVertex.position();
Expand Down
13 changes: 7 additions & 6 deletions RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc
View file
Open in desktop
Expand Up @@ -155,11 +155,11 @@ PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup


// select tracks
std::vector<reco::TransientTrack> seltks = theTrackFilter->select( t_tks );
std::vector<reco::TransientTrack> && seltks = theTrackFilter->select( t_tks );


// clusterize tracks in Z
std::vector< std::vector<reco::TransientTrack> > clusters = theTrackClusterizer->clusterize(seltks);
std::vector< std::vector<reco::TransientTrack> > && clusters = theTrackClusterizer->clusterize(seltks);
if (fVerbose){std::cout << " clustering returned "<< clusters.size() << " clusters from " << seltks.size() << " selected tracks" <<std::endl;}


Expand All @@ -168,7 +168,7 @@ PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup


std::auto_ptr<reco::VertexCollection> result(new reco::VertexCollection);
reco::VertexCollection vColl;
reco::VertexCollection & vColl = (*result);


std::vector<TransientVertex> pvs;
Expand Down Expand Up @@ -204,8 +204,11 @@ PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup
}


if (clusters.size()>2 && clusters.size() > 2*pvs.size())
edm::LogWarning("PrimaryVertexProducer") << "more than half of candidate vertices lost " << pvs.size() << ' ' << clusters.size();


if (pvs.empty() && seltks.size()>5)
edm::LogWarning("PrimaryVertexProducer") << "no vertex found with " << seltks.size() << " tracks and " << clusters.size() <<" vertex-candidates";

// sort vertices by pt**2 vertex (aka signal vertex tagging)
if(pvs.size()>1){
Expand Down Expand Up @@ -262,8 +265,6 @@ PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup
}
}


*result = vColl;
iEvent.put(result, algorithm->label);
}

Expand Down
6 changes: 1 addition & 5 deletions RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVertices_cfi.py
View file
Open in desktop
Expand Up @@ -17,7 +17,7 @@
),

TkClusParameters = cms.PSet(
algorithm = cms.string("DA"),
algorithm = cms.string("DA_vect"),
TkDAClusParameters = cms.PSet(
coolingFactor = cms.double(0.6), # moderate annealing speed
Tmin = cms.double(4.), # end of annealing
Expand Down Expand Up @@ -47,7 +47,3 @@

)

# These lines should be uncommented only for the gcc4X builds, with X>=6
offlinePrimaryVertices.TkClusParameters.algorithm = cms.string("DA_vect" )
offlinePrimaryVertices.TkClusParameters.TkDAClusParameters.use_vdt = cms.untracked.bool( True )

53 changes: 31 additions & 22 deletions RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ_vect.cc
View file
Open in desktop
Expand Up @@ -7,7 +7,7 @@
#include <cassert>
#include <limits>
#include <iomanip>

#include "FWCore/Utilities/interface/isFinite.h"
#include "vdt/vdtMath.h"

using namespace std;
Expand Down Expand Up @@ -63,10 +63,13 @@ DAClusterizerInZ_vect::track_t
DAClusterizerInZ_vect::fill(const vector<reco::TransientTrack> & tracks) const {

// prepare track data for clustering
LogDebug("DAClusterizerinZ_vectorized") << "input values";
track_t tks;
for (auto it = tracks.begin(); it!= tracks.end(); it++){
double t_pi;
if (!(*it).isValid()) continue;
double t_pi=1.;
double t_z = ((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
if (std::abs(t_z) > 1000.) continue;
Copy link
Contributor

Choose a reason for hiding this comment

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

This version of "1000." in the middle of the code is new. I would write it as a const with a suggestive name.

double phi=((*it).stateAtBeamLine().trackStateAtPCA()).momentum().phi();
double tantheta = std::tan( ((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta());
// get the beam-spot
Expand All @@ -75,13 +78,15 @@ DAClusterizerInZ_vect::fill(const vector<reco::TransientTrack> & tracks) const {
std::pow((*it).track().dzError(), 2) // track errror
+ (std::pow(beamspot.BeamWidthX()*cos(phi),2)+std::pow(beamspot.BeamWidthY()*sin(phi),2))/std::pow(tantheta,2) // beam-width
Copy link
Contributor

Choose a reason for hiding this comment

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

this could be simplified and avoid atan2 calls followed by tan or sin/cos

Copy link
Contributor Author

Choose a reason for hiding this comment

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

sure many things can be simplified, thanks for spotting (and this code was there since ages and it si now essentially unmantained)
this is only at initialization though

Copy link
Contributor

Choose a reason for hiding this comment

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

if there is some dependence on numeric precision, I can imagine that a loop from x,y to angle to x,y is a candidate to not return to the same value.
So, not doing it is more "reproducibility-safe"

Copy link
Contributor Author

Choose a reason for hiding this comment

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

fully agree. I will ask to do the math, fix and validate tomorrow.

+ std::pow(vertexSize_, 2); // intrinsic vertex size, safer for outliers and short lived decays
t_dz2 = 1./ t_dz2;
if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits<double>::min() ) continue;
if (d0CutOff_ > 0) {
Measurement1D IP =
Measurement1D atIP =
(*it).stateAtBeamLine().transverseImpactParameter();// error constains beamspot
t_pi = 1. / (1. + local_exp(std::pow(IP.value() / IP.error(), 2) - std::pow(d0CutOff_, 2))); // reduce weight for high ip tracks
} else {
t_pi = 1.;
}
t_pi = 1. / (1. + local_exp(std::pow(atIP.value() / atIP.error(), 2) - std::pow(d0CutOff_, 2))); // reduce weight for high ip tracks
if (edm::isNotFinite(t_pi) || t_pi < std::numeric_limits<double>::epsilon()) continue; // usually is > 0.99
}
LogTrace("DAClusterizerinZ_vectorized") << t_z <<' '<< t_dz2 <<' '<< t_pi;
tks.AddItem(t_z, t_dz2, &(*it), t_pi);
}
tks.ExtractRaw();
Expand All @@ -97,7 +102,7 @@ DAClusterizerInZ_vect::fill(const vector<reco::TransientTrack> & tracks) const {
namespace {
inline
double Eik(double t_z, double k_z, double t_dz2) {
return std::pow(t_z - k_z, 2) / t_dz2;
return std::pow(t_z - k_z, 2) * t_dz2;
}
}

Expand Down Expand Up @@ -126,7 +131,7 @@ double DAClusterizerInZ_vect::update(double beta, track_t & gtracks,
track_t const& tracks,
vertex_t const& vertices ) {
const double track_z = tracks._z[itrack];
const double botrack_dz2 = -beta/tracks._dz2[itrack];
const double botrack_dz2 = -beta*tracks._dz2[itrack];

// auto-vectorized
for ( unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
Expand All @@ -149,7 +154,7 @@ double DAClusterizerInZ_vect::update(double beta, track_t & gtracks,
vertex_t & y_vec ) {
auto tmp_trk_pi = tks_vec._pi[track_num];
auto o_trk_Z_sum = 1./tks_vec._Z_sum[track_num];
auto o_trk_dz2 = 1./tks_vec._dz2[track_num];
auto o_trk_dz2 = tks_vec._dz2[track_num];
auto tmp_trk_z = tks_vec._z[track_num];
auto obeta = -1./beta;

Expand Down Expand Up @@ -361,7 +366,7 @@ DAClusterizerInZ_vect::beta0(double betamax, track_t const & tks, vertex_t cons
double sumwz = 0;
double sumw = 0;
for (unsigned int i = 0; i < nt; i++) {
double w = tks._pi[i] / tks._dz2[i];
double w = tks._pi[i] * tks._dz2[i];
sumwz += w * tks._z[i];
sumw += w;
}
Expand All @@ -371,8 +376,8 @@ DAClusterizerInZ_vect::beta0(double betamax, track_t const & tks, vertex_t cons
double a = 0, b = 0;
for (unsigned int i = 0; i < nt; i++) {
double dx = tks._z[i] - y._z[k];
double w = tks._pi[i] / tks._dz2[i];
a += w * std::pow(dx, 2) / tks._dz2[i];
double w = tks._pi[i] * tks._dz2[i];
a += w * std::pow(dx, 2) * tks._dz2[i];
b += w;
}
double Tc = 2. * a / b; // the critical temperature of this vertex
Expand Down Expand Up @@ -421,7 +426,7 @@ DAClusterizerInZ_vect::split(const double beta, track_t &tks, vertex_t & y ) co
for(unsigned int i=0; i<nt; i++){
if (tks._Z_sum[i] > 0) {
double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i])) / tks._Z_sum[i];
double w=p/tks._dz2[i];
double w=p*tks._dz2[i];
if(tks._z[i]<y._z[k]){
p1+=p; z1+=w*tks._z[i]; w1+=w;
}else{
Expand Down Expand Up @@ -511,10 +516,10 @@ void DAClusterizerInZ_vect::splitAll( vertex_t & y) const {

vector<TransientVertex>
DAClusterizerInZ_vect::vertices(const vector<reco::TransientTrack> & tracks, const int verbosity) const {
track_t tks = fill(tracks);
track_t && tks = fill(tracks);
tks.ExtractRaw();

unsigned int nt = tracks.size();
unsigned int nt = tks.GetSize();
double rho0 = 0.0; // start with no outlier rejection

vector<TransientVertex> clusters;
Expand Down Expand Up @@ -617,14 +622,19 @@ DAClusterizerInZ_vect::vertices(const vector<reco::TransientTrack> & tracks, con

// ensure correct normalization of probabilities, should makes double assginment reasonably impossible
const unsigned int nv = y.GetSize();

for (unsigned int i = 0; i < nt; i++) {
for (unsigned int k = 0; k < nv; k++)
if ( edm::isNotFinite(y._pk[k]) || edm::isNotFinite(y._z[k]) ) { y._pk[k]=0; y._z[k]=0;}

for (unsigned int i = 0; i < nt; i++) // initialize
tks._Z_sum[i] = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
for (unsigned int k = 0; k < nv; k++) {

// improve vectorization (does not require reduction ....)
for (unsigned int k = 0; k < nv; k++) {
for (unsigned int i = 0; i < nt; i++)
tks._Z_sum[i] += y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],tks._dz2[i]));
}
}


for (unsigned int k = 0; k < nv; k++) {
GlobalPoint pos(0, 0, y._z[k]);

Expand All @@ -635,7 +645,6 @@ DAClusterizerInZ_vect::vertices(const vector<reco::TransientTrack> & tracks, con
double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],
tks._dz2[i])) / tks._Z_sum[i];
if ((tks._pi[i] > 0) && (p > 0.5)) {

vertexTracks.push_back(*(tks.tt[i]));
tks._Z_sum[i] = 0;
} // setting Z=0 excludes double assignment
Expand Down Expand Up @@ -730,7 +739,7 @@ void DAClusterizerInZ_vect::dump(const double beta, const vertex_t & y,
}
double tz = tks._z[i];
LogDebug("DAClusterizerinZ_vectorized") << setw(3) << i << ")" << setw(8) << fixed << setprecision(4)
<< tz << " +/-" << setw(6) << sqrt(tks._dz2[i]);
<< tz << " +/-" << setw(6) << sqrt(1./tks._dz2[i]);

if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) {
LogDebug("DAClusterizerinZ_vectorized") << " *";
Expand Down