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
Avoid Nans in primary vertex #5490
Merged
Merged
Changes from all commits
Commits
Show all changes
7 commits
Select commit
Hold shift + click to select a range
c975167
robistify DA
VinInn 4371659
compute pt once
VinInn 0dd30dc
optimize sorting, add warning
VinInn d233608
avoid to trash memory
VinInn 9dae5d1
clean config
VinInn 44dec68
remove g
VinInn 23ba733
avoid trigonometry
VinInn File filter
Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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();} | ||
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; | ||
|
||
|
@@ -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. | ||
|
@@ -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 | ||
|
@@ -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); | ||
} | ||
|
||
|
@@ -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); | ||
} | ||
|
||
|
@@ -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); | ||
} | ||
|
||
|
@@ -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. | ||
|
@@ -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; | ||
} | ||
|
@@ -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; | ||
|
@@ -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. || | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. probably it's a silly point, There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. If somebody volunteers (the code was there already) |
||
nVertex.position().perp()>120.) | ||
{ | ||
|
@@ -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(); | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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) ?
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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.