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 1 commit
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.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
optimize sorting, add warning
  • Loading branch information
@VinInn
VinInn committed Sep 21, 2014
commit 0dd30dc750b010366c9723c8926b5e8826aff80e
60 changes: 24 additions & 36 deletions RecoVertex/AdaptiveVertexFit/src/AdaptiveVertexFitter.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,14 @@ using namespace std;

namespace {
void sortTracksByPt(std::vector<reco::TransientTrack> & cont) {
auto b = &cont.front();
auto e = b + cont.size();
using tk = reco::TransientTrack const &;
float pt2[e-b];
for (auto p=b; p!=e; ++p) pt2[p-b] = p->impactPointState().globalMomentum().perp2();
std::sort(b,e,[&](tk x, tk y){ return pt2[&x-b]>pt2[&y-b];} );
}
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 Down Expand Up @@ -57,22 +58,16 @@ namespace {



struct DistanceToRefPoint {
DistanceToRefPoint ( const GlobalPoint & ref ) : theRef ( ref ) {}

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);
}



Expand Down Expand Up @@ -439,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 @@ -536,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 @@ -588,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 @@ -614,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
4 changes: 3 additions & 1 deletion RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,7 +206,9 @@ 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
2 changes: 1 addition & 1 deletion RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ_vect.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -623,7 +623,7 @@ 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 k = 0; k < nv; k++)
if ( edm::isNotFinite(y._pk[k]) || edm::isNotFinite(y._z[k]) ) y._pk[k]=0;
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_);
Expand Down