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

Updates to precision timing integration in reconstruction and miniaod #20338

Merged
merged 51 commits into from
Sep 23, 2017
Merged
Show file tree
Hide file tree
Changes from 47 commits
Commits
Show all changes
51 commits
Select commit Hold shift + click to select a range
604ef34
SimTracker/TrackAssociation: don't propagate NaNs to the ValueMap
gpetruc Dec 6, 2016
96b03b9
use timing for primary vertex assignment and sorting
bendavid Dec 6, 2016
22bda77
remove debug printout
bendavid Dec 7, 2016
d3c2bd3
Use the time out of PFCandidates when running the PF sorter, and rear…
gpetruc Dec 7, 2016
091885f
RecoParticleFlow/PFClusterProducer: propagate time to cluster if poss…
gpetruc Dec 8, 2016
e8b4b67
Import GSF track time in PF (phase2)
gpetruc Dec 8, 2016
7d80e04
Combine GSF and Cluster time in PFEGammaAlgo
gpetruc Dec 8, 2016
cdd3a86
Fix input tag for GSF track value map
gpetruc Dec 9, 2016
7f47d0a
Fix PFEGammaAlgo: read track timeError using timeError(), not time()
gpetruc Dec 9, 2016
b555e21
PFClusterTimeAssigner to propagate also the time uncertainty
gpetruc Dec 9, 2016
eaa2695
SimPFProducer: use also PFCluster time for electrons
gpetruc Dec 9, 2016
8956001
neglect cluster time for electrons for now (PFCluster time is buggy i…
bendavid Jan 19, 2017
89a315f
protect against pathological case in PrimaryVertexAssignment (since f…
bendavid Jan 22, 2017
4a37aac
split TrackTimeValueMap producer into seperate instances for generalT…
bendavid Jan 24, 2017
d12c8a9
add random number source for gsfTrackTimeValueMap
bendavid Jan 24, 2017
5107342
store hard process position and time in miniaod (4 floats per event)
bendavid Mar 3, 2017
e4e879e
add eta bounds to trackTimeValueMapProducer
lgray Nov 11, 2016
ee49778
impose realistic acceptance cuts for track timing maps by default
bendavid Mar 14, 2017
d91029e
make track-time importer that operates on all tracks after complete i…
lgray Mar 15, 2017
378e0e2
migrate also the gsf track time assignment in pf to the new importer
bendavid Mar 22, 2017
a2fe468
remove unused vector
bendavid Mar 22, 2017
989dcce
properly assign no valid timestamp to tracks out of acceptance, and a…
bendavid Mar 22, 2017
6b92f6b
fix incorrect weighting of tracks from time uncertainty
bendavid Mar 31, 2017
91645db
consistently use twice the beamspot width as the uncertainty for trac…
bendavid Mar 31, 2017
15b1308
remove leftover use of timing since it has moved elsewhere
bendavid Aug 30, 2017
e34878e
fix logic error in pv assignment
bendavid Aug 30, 2017
9ea1aff
remove extra pt cut for 4d vertexing
bendavid Aug 30, 2017
b51de43
update 4d vertexing sequence, remove redundant 1d collection and add …
bendavid Aug 30, 2017
f21f833
revert PV association to standard vertices for timing scenario
bendavid Aug 31, 2017
2e207f8
add missing parameter
bendavid Aug 31, 2017
e02bf9d
modify PrimaryVertexAssignment logic to preserve efficiency for track…
bendavid Mar 4, 2017
aeb1777
preserve default behaviour for vertex assignment in miniaod
bendavid Aug 31, 2017
4b6fb7a
Packed timing in packed candidates
gpetruc Dec 6, 2016
a4f4a18
In phase2_timing scenario, propagate timing to packed pf candidates
gpetruc Dec 6, 2016
e23adf9
Fix copy constructor of PackedCandidate
gpetruc Dec 6, 2016
253c0fa
Unit test for time packing (with optional verbose printing out of ran…
gpetruc Dec 6, 2016
813454e
Fix initializer order and checksum for PackedCandidates
bendavid Aug 31, 2017
bd5df20
fix for time assignment of displaced tracks
bendavid Aug 25, 2017
e33b8ad
fix units
bendavid Aug 25, 2017
92e7a03
remove debug printout
bendavid Aug 25, 2017
68d2027
don't use cluster time also in PFAlgo, since it's not reliable yet
bendavid Sep 1, 2017
caa0dcb
fix parameter set name in vertex config
bendavid Aug 31, 2017
86878cb
code checks
bendavid Sep 12, 2017
1b346e3
restore pt>0.7 GeV cut in 4d vertexing pending further performance tu…
bendavid Sep 15, 2017
8bfd57a
minor code cleanup
bendavid Sep 15, 2017
7ca1d5a
Document time encoding in PackedCandidate, fix some corner cases
gpetruc Sep 15, 2017
bf86758
remove spurious size_t ipv=0 argument in PackedCandidate::time()
gpetruc Sep 15, 2017
2e974d5
minor cleanup and variable name rationalization
bendavid Sep 19, 2017
e4cf917
remove currently unused code for pf candidate time assignment from pf…
bendavid Sep 19, 2017
96ee090
do not reorder elements when assigning time
lgray Sep 21, 2017
213e0b3
Merge pull request #4 from lgray/timingupdate93x
bendavid Sep 21, 2017
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
31 changes: 26 additions & 5 deletions CommonTools/RecoAlgos/interface/PrimaryVertexAssignment.h
Original file line number Diff line number Diff line change
Expand Up @@ -23,65 +23,86 @@ class PrimaryVertexAssignment {
maxDzSigForPrimaryAssignment_(iConfig.getParameter<double>("maxDzSigForPrimaryAssignment")),
maxDzForPrimaryAssignment_(iConfig.getParameter<double>("maxDzForPrimaryAssignment")),
maxDzErrorForPrimaryAssignment_(iConfig.getParameter<double>("maxDzErrorForPrimaryAssignment")),
maxDtSigForPrimaryAssignment_(iConfig.getParameter<double>("maxDtSigForPrimaryAssignment")),
maxJetDeltaR_(iConfig.getParameter<double>("maxJetDeltaR")),
minJetPt_(iConfig.getParameter<double>("minJetPt")),
maxDistanceToJetAxis_(iConfig.getParameter<double>("maxDistanceToJetAxis")),
maxDzForJetAxisAssigment_(iConfig.getParameter<double>("maxDzForJetAxisAssigment")),
maxDxyForJetAxisAssigment_(iConfig.getParameter<double>("maxDxyForJetAxisAssigment")),
maxDxySigForNotReconstructedPrimary_(iConfig.getParameter<double>("maxDxySigForNotReconstructedPrimary")),
maxDxyForNotReconstructedPrimary_(iConfig.getParameter<double>("maxDxyForNotReconstructedPrimary"))
maxDxyForNotReconstructedPrimary_(iConfig.getParameter<double>("maxDxyForNotReconstructedPrimary")),
useTiming_(iConfig.getParameter<bool>("useTiming")),
preferHighRanked_(iConfig.getParameter<bool>("preferHighRanked"))
{}

~PrimaryVertexAssignment(){}

std::pair<int,PrimaryVertexAssignment::Quality> chargedHadronVertex(const reco::VertexCollection& vertices,
const reco::TrackRef& trackRef,
const reco::Track * track,
float trackTime,
float trackTimeResolution, // <0 if timing not available for this object
const edm::View<reco::Candidate> & jets,
const TransientTrackBuilder & builder) const;

std::pair<int,PrimaryVertexAssignment::Quality> chargedHadronVertex(const reco::VertexCollection& vertices,
const reco::TrackRef& trackRef,
float trackTime,
float trackTimeResolution, // <0 if timing not available for this object
const edm::View<reco::Candidate> & jets,
const TransientTrackBuilder & builder) const
{
return chargedHadronVertex(vertices,trackRef,&(*trackRef),jets,builder);
return chargedHadronVertex(vertices,trackRef,&(*trackRef),trackTime,trackTimeResolution,jets,builder);
}

std::pair<int,PrimaryVertexAssignment::Quality> chargedHadronVertex( const reco::VertexCollection& vertices,
const reco::PFCandidate& pfcand,
const edm::View<reco::Candidate>& jets,
const TransientTrackBuilder& builder) const {
float time = 0, timeResolution = -1;
if (useTiming_ && pfcand.isTimeValid()) {
time = pfcand.time(); timeResolution = pfcand.timeError();
}
if(pfcand.gsfTrackRef().isNull())
{
if(pfcand.trackRef().isNull())
return std::pair<int,PrimaryVertexAssignment::Quality>(-1,PrimaryVertexAssignment::Unassigned);
else
return chargedHadronVertex(vertices,pfcand.trackRef(),jets,builder);
return chargedHadronVertex(vertices,pfcand.trackRef(),time,timeResolution,jets,builder);
}
return chargedHadronVertex(vertices,reco::TrackRef(),&(*pfcand.gsfTrackRef()),jets,builder);
return chargedHadronVertex(vertices,reco::TrackRef(),&(*pfcand.gsfTrackRef()),time,timeResolution,jets,builder);
}
std::pair<int,PrimaryVertexAssignment::Quality> chargedHadronVertex( const reco::VertexCollection& vertices,
const reco::RecoChargedRefCandidate& chcand,
const edm::ValueMap<float> *trackTimeTag,
const edm::ValueMap<float> *trackTimeResoTag,
const edm::View<reco::Candidate>& jets,
const TransientTrackBuilder& builder) const {
float time = 0, timeResolution = -1;
if (useTiming_) {
time = (*trackTimeTag)[chcand.track()];
timeResolution = (*trackTimeResoTag)[chcand.track()];
}
if(chcand.track().isNull())
return std::pair<int,PrimaryVertexAssignment::Quality>(-1,PrimaryVertexAssignment::Unassigned);
return chargedHadronVertex(vertices,chcand.track(),jets,builder);
return chargedHadronVertex(vertices,chcand.track(),time,timeResolution,jets,builder);
}


private :
double maxDzSigForPrimaryAssignment_;
double maxDzForPrimaryAssignment_;
double maxDzErrorForPrimaryAssignment_;
double maxDtSigForPrimaryAssignment_;
double maxJetDeltaR_;
double minJetPt_;
double maxDistanceToJetAxis_;
double maxDzForJetAxisAssigment_;
double maxDxyForJetAxisAssigment_;
double maxDxySigForNotReconstructedPrimary_;
double maxDxyForNotReconstructedPrimary_;
bool useTiming_;
bool preferHighRanked_;
};

#endif
70 changes: 66 additions & 4 deletions CommonTools/RecoAlgos/plugins/PrimaryVertexSorter.h
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@

*/


template <class ParticlesCollection>

class PrimaryVertexSorter : public edm::stream::EDProducer<> {
Expand Down Expand Up @@ -55,13 +54,21 @@ class PrimaryVertexSorter : public edm::stream::EDProducer<> {
/// vertices
edm::EDGetTokenT<reco::VertexCollection> tokenVertices_;
edm::EDGetTokenT<edm::View<reco::Candidate> > tokenJets_;
edm::EDGetTokenT<edm::ValueMap<float> > tokenTrackTimeTag_;
edm::EDGetTokenT<edm::ValueMap<float> > tokenTrackTimeResoTag_;

bool produceOriginalMapping_;
bool produceSortedVertices_;
bool producePFPileUp_;
bool producePFNoPileUp_;
int qualityCut_;
bool useMET_;
bool useTiming_;

void doConsumesForTiming(const edm::ParameterSet &iConfig) ;
bool needsProductsForTiming() ;
std::pair<int,PrimaryVertexAssignment::Quality> runAlgo( const reco::VertexCollection& vertices, const typename ParticlesCollection::value_type & pf, const edm::ValueMap<float> *trackTimeTag,
const edm::ValueMap<float> *trackTimeResoTag, const edm::View<reco::Candidate>& jets, const TransientTrackBuilder& builder) ;
};


Expand Down Expand Up @@ -89,7 +96,8 @@ PrimaryVertexSorter<ParticlesCollection>::PrimaryVertexSorter(const edm::Paramet
producePFPileUp_(iConfig.getParameter<bool>("producePileUpCollection")),
producePFNoPileUp_(iConfig.getParameter<bool>("produceNoPileUpCollection")),
qualityCut_(iConfig.getParameter<int>("qualityForPrimary")),
useMET_(iConfig.getParameter<bool>("usePVMET"))
useMET_(iConfig.getParameter<bool>("usePVMET")),
useTiming_(iConfig.getParameterSet("assignment").getParameter<bool>("useTiming"))
{

using namespace std;
Expand Down Expand Up @@ -122,6 +130,7 @@ using namespace reco;
produces< PFCollection> ("NoPileUp");
}

if (useTiming_) doConsumesForTiming(iConfig);

}

Expand All @@ -148,7 +157,20 @@ using namespace reco;

Handle<ParticlesCollection> particlesHandle;
iEvent.getByToken( tokenCandidates_, particlesHandle);


Handle<edm::ValueMap<float> > trackTimeTagHandle;
Handle<edm::ValueMap<float> > trackTimeResoTagHandle;

const edm::ValueMap<float> *trackTimeTag = 0;
const edm::ValueMap<float> *trackTimeResoTag = 0;
if (useTiming_ && needsProductsForTiming()) {
iEvent.getByToken(tokenTrackTimeTag_, trackTimeTagHandle);
iEvent.getByToken(tokenTrackTimeResoTag_, trackTimeResoTagHandle);

trackTimeTag = trackTimeTagHandle.product();
trackTimeResoTag = trackTimeResoTagHandle.product();
}

ParticlesCollection particles = *particlesHandle.product();
std::vector<int> pfToPVVector;
std::vector<PrimaryVertexAssignment::Quality> pfToPVQualityVector;
Expand All @@ -160,7 +182,7 @@ using namespace reco;
std::vector<float> vertexScore(vertices->size());

for(auto const & pf : particles) {
std::pair<int,PrimaryVertexAssignment::Quality> vtxWithQuality=assignmentAlgo_.chargedHadronVertex(*vertices,pf,*jets,*builder);
std::pair<int,PrimaryVertexAssignment::Quality> vtxWithQuality = runAlgo(*vertices,pf,trackTimeTag,trackTimeResoTag,*jets,*builder);
pfToPVVector.push_back(vtxWithQuality.first);
pfToPVQualityVector.push_back(vtxWithQuality.second);
}
Expand Down Expand Up @@ -297,4 +319,44 @@ using namespace reco;
}


template<>
void PrimaryVertexSorter<std::vector<reco::RecoChargedRefCandidate>>::doConsumesForTiming(const edm::ParameterSet &iConfig)
{
tokenTrackTimeTag_ = consumes<edm::ValueMap<float> > (iConfig.getParameter<edm::InputTag>("trackTimeTag"));
tokenTrackTimeResoTag_ = consumes<edm::ValueMap<float> > (iConfig.getParameter<edm::InputTag>("trackTimeResoTag"));
}

template<>
void PrimaryVertexSorter<std::vector<reco::PFCandidate>>::doConsumesForTiming(const edm::ParameterSet &iConfig)
{
}

template<>
bool PrimaryVertexSorter<std::vector<reco::RecoChargedRefCandidate>>::needsProductsForTiming()
{
return true;
}

template<>
bool PrimaryVertexSorter<std::vector<reco::PFCandidate>>::needsProductsForTiming()
{
return false;
}

template<>
std::pair<int,PrimaryVertexAssignment::Quality>
PrimaryVertexSorter<std::vector<reco::RecoChargedRefCandidate>>::runAlgo( const reco::VertexCollection& vertices, const reco::RecoChargedRefCandidate & pf, const edm::ValueMap<float> *trackTimeTag,
const edm::ValueMap<float> *trackTimeResoTag, const edm::View<reco::Candidate>& jets, const TransientTrackBuilder& builder)
{
return assignmentAlgo_.chargedHadronVertex( vertices, pf, trackTimeTag, trackTimeResoTag, jets, builder);
}

template<>
std::pair<int,PrimaryVertexAssignment::Quality>
PrimaryVertexSorter<std::vector<reco::PFCandidate>>::runAlgo( const reco::VertexCollection& vertices, const reco::PFCandidate & pf, const edm::ValueMap<float> *trackTimeTag,
const edm::ValueMap<float> *trackTimeResoTag, const edm::View<reco::Candidate>& jets, const TransientTrackBuilder& builder)
{
return assignmentAlgo_.chargedHadronVertex( vertices, pf, jets, builder);
}

#endif
4 changes: 3 additions & 1 deletion CommonTools/RecoAlgos/python/sortedPFPrimaryVertices_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
maxDzSigForPrimaryAssignment = cms.double(5.0), # in AND with next
maxDzForPrimaryAssignment = cms.double(0.1), # in AND with prev
maxDzErrorForPrimaryAssignment = cms.double(0.05), # in AND with prev, tracks with uncertainty above 500um cannot tell us which pv they come from

maxDtSigForPrimaryAssignment = cms.double(3.0), # *FIXME* this parameter needs to be double checked before useTiming is switched back on
# cuts used to recover b-tracks if they are closed to jet axis
maxJetDeltaR = cms.double(0.5),
minJetPt = cms.double(25),
Expand All @@ -17,6 +17,8 @@
#cuts used to identify primary tracks compatible with beamspot
maxDxySigForNotReconstructedPrimary = cms.double(2), #in AND with next
maxDxyForNotReconstructedPrimary = cms.double(0.01), #in AND with prev
useTiming = cms.bool(False),
preferHighRanked = cms.bool(False)
),
particles = cms.InputTag("particleFlow"),
vertices= cms.InputTag("offlinePrimaryVertices"),
Expand Down
6 changes: 5 additions & 1 deletion CommonTools/RecoAlgos/python/sortedPrimaryVertices_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
maxDzSigForPrimaryAssignment = cms.double(5.0), # in AND with next
maxDzForPrimaryAssignment = cms.double(0.1), # in AND with prev
maxDzErrorForPrimaryAssignment = cms.double(0.05), # in AND with prev, tracks with uncertainty above 500um cannot tell us which pv they come from

maxDtSigForPrimaryAssignment = cms.double(4.0),
# cuts used to recover b-tracks if they are closed to jet axis
maxJetDeltaR = cms.double(0.5),
minJetPt = cms.double(25),
Expand All @@ -17,8 +17,12 @@
#cuts used to identify primary tracks compatible with beamspot
maxDxySigForNotReconstructedPrimary = cms.double(2), #in AND with next
maxDxyForNotReconstructedPrimary = cms.double(0.01), #in AND with prev
useTiming = cms.bool(False),
preferHighRanked = cms.bool(False),
),
particles = cms.InputTag("trackRefsForJets"),
trackTimeTag = cms.InputTag(""),
trackTimeResoTag = cms.InputTag(""),
vertices= cms.InputTag("offlinePrimaryVertices"),
# Jets= cms.InputTag("ak4PFJets"),
jets= cms.InputTag("ak4CaloJetsForTrk"),
Expand Down
100 changes: 76 additions & 24 deletions CommonTools/RecoAlgos/src/PrimaryVertexAssignment.cc
Original file line number Diff line number Diff line change
Expand Up @@ -4,50 +4,102 @@
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "TrackingTools/IPTools/interface/IPTools.h"

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


std::pair<int,PrimaryVertexAssignment::Quality>
PrimaryVertexAssignment::chargedHadronVertex( const reco::VertexCollection& vertices,
const reco::TrackRef& trackRef,
const reco::Track* track,
float time,
float timeReso, // <0 if timing not available for this object
const edm::View<reco::Candidate>& jets,
const TransientTrackBuilder& builder) const {

int iVertex = -1;
size_t index=0;
typedef reco::VertexCollection::const_iterator IV;
typedef reco::Vertex::trackRef_iterator IT;

int iVertex = -1;
size_t index=0;
float bestweight=0;
for( auto const & vtx : vertices) {
float w = vtx.trackWeight(trackRef);
if (w > bestweight){
bestweight=w;
iVertex=index;
}
bestweight=w;
iVertex=index;
}
index++;
}

bool useTime = useTiming_;
if (edm::isNotFinite(time) || timeReso<1e-6) {
useTime = false;
time = 0.;
timeReso = -1.;
}

if (preferHighRanked_) {
for(IV iv=vertices.begin(); iv!=vertices.end(); ++iv) {
int ivtx = iv - vertices.begin();
if (iVertex == ivtx) return std::pair<int,PrimaryVertexAssignment::Quality>(ivtx,PrimaryVertexAssignment::UsedInFit);

double dz = std::abs(track->dz(iv->position()));
double dt = std::abs(time-iv->t());

bool useTimeVtx = useTime && iv->tError()>0.;

if ((dz < maxDzForPrimaryAssignment_ or dz/track->dzError() < maxDzSigForPrimaryAssignment_ ) and (!useTimeVtx or dt/timeReso < maxDtSigForPrimaryAssignment_)) {
Copy link
Contributor

Choose a reason for hiding this comment

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

why is the vertex error not a part of the significance computation here (it is in other places in this file)?

Copy link
Contributor

Choose a reason for hiding this comment

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

the vertex supposedly has a time resolution at this point.
Shouldn't it be used as well in the significance computation?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The vertex time error is only checked if it's >0 (ie to check that the vertex has valid timing info.)

The numerical value is intentionally not used in order to avoid preferring to assign tracks to vertices with large time uncertainty (and the logic is coherent for z and t).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

(where the dz significance is used later in the code including the vertex error, it will not cause poor vertices to be favoured, because it is only checked against the closest vertex)

Copy link
Contributor

Choose a reason for hiding this comment

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

the reason to add the vtx z error was to "regularize" the selection for tracks with very small dzError. IIRC, values down to 5 um are possible which can simply kill a good track with respect to still rather reasonable vertices.

return std::pair<int,PrimaryVertexAssignment::Quality>(ivtx,PrimaryVertexAssignment::PrimaryDz);
}
}
}


if(iVertex >= 0 ) return std::pair<int,PrimaryVertexAssignment::Quality>(iVertex,PrimaryVertexAssignment::UsedInFit);

double dzmin = 1e99;
int vtxIdMinDz = -1;
for(IV iv=vertices.begin(); iv!=vertices.end(); ++iv) {
double dz = std::abs(track->dz(iv->position()));
if(dz<dzmin) {
dzmin = dz;
vtxIdMinDz = iv-vertices.begin();
}
}

double distmin = std::numeric_limits<double>::max();
double dzmin = std::numeric_limits<double>::max();
double dtmin = std::numeric_limits<double>::max();
int vtxIdMinDist = -1;
Copy link
Contributor

Choose a reason for hiding this comment

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

vtxIdMinSignif may be a more intuitive name.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This is vestigial from the pre-timing version of the algorithm. Can we leave it for future cleanup?

Copy link
Contributor

Choose a reason for hiding this comment

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

OK

for(IV iv=vertices.begin(); iv!=vertices.end(); ++iv) {
double dz = std::abs(track->dz(iv->position()));
double dt = std::abs(time-iv->t());

double dzsig = dz/track->dzError();
double dist = dzsig*dzsig;
Copy link
Contributor

Choose a reason for hiding this comment

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

why is the vertex uncertainty not a part of the computation?
Is it because of more frequent possible assignment to poorly reconstructed vertices?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, that is the reason.

Actually I do recall now that the vertex uncertainty IS used in once place later in the file. Let me double check what the reason was.


bool useTimeVtx = useTime && iv->tError()>0.;
if (useTime && !useTimeVtx) {
dt = timeReso;
}

if (useTime) {
double dtsig = dt/timeReso;

dist += dtsig*dtsig;
}
if(dist<distmin) {
distmin = dist;
dzmin = dz;
dtmin = dt;
vtxIdMinDist = iv-vertices.begin();
}
}

// first use "closest in Z" with tight cuts (targetting primary particles)
float dzE=sqrt(track->dzError()*track->dzError()+vertices[vtxIdMinDz].covariance(2,2));
if(vtxIdMinDz>=0 and (dzmin < maxDzForPrimaryAssignment_ and dzmin/dzE < maxDzSigForPrimaryAssignment_ and track->dzError()<maxDzErrorForPrimaryAssignment_))
const float add_cov = vtxIdMinDist >= 0 ? vertices[vtxIdMinDist].covariance(2,2) : 0.f;
const float dzE=sqrt(track->dzError()*track->dzError()+add_cov);
if(vtxIdMinDist>=0 and
(dzmin < maxDzForPrimaryAssignment_ and dzmin/dzE < maxDzSigForPrimaryAssignment_ and track->dzError()<maxDzErrorForPrimaryAssignment_) and
Copy link
Contributor Author

Choose a reason for hiding this comment

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

There is in fact a small issue which has been introduced here when forward porting this code to 93x and resolving merge conflicts with 9d480b8

It would not affect the non-timing vertices, but could have a small impact on the sorting of the 4d vertices. I'll commit a fix momentarily.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sorry, that's not true, I was misreading/misremembering what this code does.

Everything which is here is fine, with indeed a small inconsistency between the handling of z and t, since the vertex error is used here for z, but not for t. I would prefer to leave it like this for now, also because we need to study more closely the pulls on the time errors of the 4d vertices.

(!useTime or dtmin/timeReso < maxDtSigForPrimaryAssignment_) )
{
iVertex=vtxIdMinDz;
iVertex=vtxIdMinDist;
}

if(iVertex >= 0 ) return std::pair<int,PrimaryVertexAssignment::Quality>(iVertex,PrimaryVertexAssignment::PrimaryDz);



// if track not assigned yet, it could be a b-decay secondary , use jet axis dist criterion
// first find the closest jet within maxJetDeltaR_
int jetIdx = -1;
Expand Down Expand Up @@ -94,14 +146,14 @@ PrimaryVertexAssignment::chargedHadronVertex( const reco::VertexCollection& vert

// if the track is not compatible with other PVs but is compatible with the BeamSpot, we may simply have not reco'ed the PV!
// we still point it to the closest in Z, but flag it as possible orphan-primary
if(std::abs(track->dxy(vertices[0].position()))<maxDxyForNotReconstructedPrimary_ && std::abs(track->dxy(vertices[0].position())/track->dxyError())<maxDxySigForNotReconstructedPrimary_)
return std::pair<int,PrimaryVertexAssignment::Quality>(vtxIdMinDz,PrimaryVertexAssignment::NotReconstructedPrimary);
if(!vertices.empty() && std::abs(track->dxy(vertices[0].position()))<maxDxyForNotReconstructedPrimary_ && std::abs(track->dxy(vertices[0].position())/track->dxyError())<maxDxySigForNotReconstructedPrimary_)
return std::pair<int,PrimaryVertexAssignment::Quality>(vtxIdMinDist,PrimaryVertexAssignment::NotReconstructedPrimary);

//FIXME: here we could better handle V0s and NucInt

// all other tracks could be non-B secondaries and we just attach them with closest Z
if(vtxIdMinDz>=0)
return std::pair<int,PrimaryVertexAssignment::Quality>(vtxIdMinDz,PrimaryVertexAssignment::OtherDz);
if(vtxIdMinDist>=0)
return std::pair<int,PrimaryVertexAssignment::Quality>(vtxIdMinDist,PrimaryVertexAssignment::OtherDz);
//If for some reason even the dz failed (when?) we consider the track not assigned
return std::pair<int,PrimaryVertexAssignment::Quality>(-1,PrimaryVertexAssignment::Unassigned);
}
Expand Down
Loading