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

Do not export missing hits, avoid extra material effect application in backward-fit #376

Merged
merged 2 commits into from
Oct 27, 2021
Merged

Do not export missing hits, avoid extra material effect application in backward-fit #376

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
649a526
For clone-engine keep candidate state at last valid hit.
osschar Oct 26, 2021
ee825c6
Merge branch 'devel' into no_missing_hits_und_no_mateff_bkfit
mmasciov Oct 27, 2021
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
20 changes: 13 additions & 7 deletions mkFit/HitStructures.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -445,27 +445,33 @@ EventOfHits::EventOfHits(const TrackerInfo &trk_inf) :
// TrackCand
//==============================================================================

Track TrackCand::exportTrack() const
Track TrackCand::exportTrack(bool remove_missing_hits) const
{
// printf("TrackCand::exportTrack label=%5d, total_hits=%2d, overlaps=%2d -- n_seed_hits=%d,prod_type=%d\n",
// label(), nTotalHits(), nOverlapHits_, getNSeedHits(), (int)prodType());

Track res(*this);
res.resizeHits(nTotalHits(), nFoundHits());
res.resizeHits(remove_missing_hits ? nFoundHits() : nTotalHits(), nFoundHits());
res.setNOverlapHits(nOverlapHits());

int nh = nTotalHits();
int ch = lastHitIdx_;
int good_hits_pos = nFoundHits();
while (--nh >= 0)
{
HoTNode& hot_node = m_comb_candidate->m_hots[ch];

res.setHitIdxAtPos(nh, hot_node.m_hot);

if (remove_missing_hits)
{
if (hot_node.m_hot.index >= 0)
res.setHitIdxAtPos(--good_hits_pos, hot_node.m_hot);
}
else
{
res.setHitIdxAtPos(nh, hot_node.m_hot);
}
// printf(" nh=%2d, ch=%d, idx=%d lyr=%d prev_idx=%d\n",
// nh, ch, hot_node.m_hot.index, hot_node.m_hot.layer, hot_node.m_prev_idx);

ch = hot_node.m_prev_idx;
ch = hot_node.m_prev_idx;
}

return res;
Expand Down
2 changes: 1 addition & 1 deletion mkFit/HitStructures.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -528,7 +528,7 @@ class TrackCand : public TrackBase

void incOverlapCount() { ++nOverlapHits_; }

Track exportTrack() const;
Track exportTrack(bool remove_missing_hits=false) const;

void resetShortTrack() { score_ = getScoreWorstPossible(); m_comb_candidate = nullptr; }

Expand Down
10 changes: 6 additions & 4 deletions mkFit/MkBase.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,8 @@ class MkBase
Err[iP], Par[iP], N_proc, pf);
}

void PropagateTracksToHitR(const MPlexHV& par, const int N_proc, const PropagationFlags pf)
void PropagateTracksToHitR(const MPlexHV& par, const int N_proc, const PropagationFlags pf,
const MPlexQI *noMatEffPtr=nullptr)
{
MPlexQF msRad;
#pragma omp simd
Expand All @@ -55,7 +56,7 @@ class MkBase
}

propagateHelixToRMPlex(Err[iC], Par[iC], Chg, msRad,
Err[iP], Par[iP], N_proc, pf);
Err[iP], Par[iP], N_proc, pf, noMatEffPtr);
}

//----------------------------------------------------------------------------
Expand All @@ -73,7 +74,8 @@ class MkBase
Err[iP], Par[iP], N_proc, pf);
}

void PropagateTracksToHitZ(const MPlexHV& par, const int N_proc, const PropagationFlags pf)
void PropagateTracksToHitZ(const MPlexHV& par, const int N_proc, const PropagationFlags pf,
const MPlexQI *noMatEffPtr=nullptr)
{
MPlexQF msZ;
#pragma omp simd
Expand All @@ -83,7 +85,7 @@ class MkBase
}

propagateHelixToZMPlex(Err[iC], Par[iC], Chg, msZ,
Err[iP], Par[iP], N_proc, pf);
Err[iP], Par[iP], N_proc, pf, noMatEffPtr);
}

void PropagateTracksToPCAZ(const int N_proc, const PropagationFlags pf)
Expand Down
13 changes: 8 additions & 5 deletions mkFit/MkBuilder.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -303,14 +303,14 @@ int MkBuilder::filter_comb_cands(std::function<filter_track_cand_foo> filter)
return n_removed;
}

void MkBuilder::select_best_comb_cands(bool clear_m_tracks)
void MkBuilder::select_best_comb_cands(bool clear_m_tracks, bool remove_missing_hits)
{
if (clear_m_tracks)
m_tracks.clear();
export_best_comb_cands(m_tracks);
export_best_comb_cands(m_tracks, remove_missing_hits);
}

void MkBuilder::export_best_comb_cands(TrackVec &out_vec)
void MkBuilder::export_best_comb_cands(TrackVec &out_vec, bool remove_missing_hits)
{
const EventOfCombCandidates &eoccs = m_event_of_comb_cands;
out_vec.reserve(out_vec.size() + eoccs.m_size);
Expand All @@ -323,7 +323,7 @@ void MkBuilder::export_best_comb_cands(TrackVec &out_vec)
if ( ! eoccs[i].empty())
{
const TrackCand &bcand = eoccs[i].front();
out_vec.emplace_back( bcand.exportTrack() );
out_vec.emplace_back( bcand.exportTrack(remove_missing_hits) );
}
}
}
Expand Down Expand Up @@ -1961,7 +1961,10 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr,
// }

// copy_out the propagated track params, errors only.
mkfndr->CopyOutParErr(eoccs.m_candidates, end - itrack, true);
// Do not, keep cands at last valid hit until actual update,
// this requires change to propagation flags used in MkFinder::UpdateWithLastHit()
// from intra-layer to inter-layer.
// mkfndr->CopyOutParErr(eoccs.m_candidates, end - itrack, true);

dprint("make new candidates");
cloner.begin_iteration();
Expand Down
4 changes: 2 additions & 2 deletions mkFit/MkBuilder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,8 +160,8 @@ class MkBuilder

// XXX filter for rearranging cands that will / will not do backward search.

void select_best_comb_cands(bool clear_m_tracks=false);
void export_best_comb_cands(TrackVec &out_vec);
void select_best_comb_cands(bool clear_m_tracks=false, bool remove_missing_hits=false);
void export_best_comb_cands(TrackVec &out_vec, bool remove_missing_hits=false);
void export_tracks(TrackVec &out_vec);

void CompactifyHitStorageForBestCand(bool remove_seed_hits, int backward_fit_min_hits)
Expand Down
21 changes: 15 additions & 6 deletions mkFit/MkFinder.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1249,8 +1249,10 @@ void MkFinder::UpdateWithLastHit(const LayerOfHits &layer_of_hits, int N_proc,
msPar.CopyIn(i, hit.posArray());
}

// See comment in MkBuilder::find_tracks_in_layer() about intra / inter flags used here
// for propagation to the hit.
(*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar,
Err[iC], Par[iC], N_proc, Config::finding_intra_layer_pflags);
Err[iC], Par[iC], N_proc, Config::finding_inter_layer_pflags);
}


Expand Down Expand Up @@ -1527,6 +1529,7 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits,
// Then we could avoid checking which layers actually do have hits.

MPlexQF tmp_chi2;
MPlexQI no_mat_effs;
float tmp_err[6] = { 666, 0, 666, 0, 0, 666 };
float tmp_pos[3];

Expand All @@ -1542,14 +1545,18 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits,
const Hit *last_hit_ptr[NN];
#endif

int count = 0;
no_mat_effs.SetVal(0);
int done_count = 0;
int here_count = 0;
for (int i = 0; i < N_proc; ++i)
{
while (CurNode[i] >= 0 && HoTNodeArr[ i ][ CurNode[i] ].m_hot.index < 0)
{
CurNode[i] = HoTNodeArr[ i ][ CurNode[i] ].m_prev_idx;
}

if (CurNode[i] < 0) ++done_count;

if (CurNode[i] >= 0 && HoTNodeArr[ i ][ CurNode[i] ].m_hot.layer == layer)
{
// Skip the overlap hits -- if they exist.
Expand All @@ -1569,7 +1576,7 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits,
#endif
msErr.CopyIn(i, hit.errArray());
msPar.CopyIn(i, hit.posArray());
++count;
++here_count;

CurNode[i] = HoTNodeArr[ i ][ CurNode[i] ].m_prev_idx;
}
Expand All @@ -1579,6 +1586,7 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits,
#ifdef DEBUG_BACKWARD_FIT
last_hit_ptr[i] = nullptr;
#endif
no_mat_effs[i] = 1;
tmp_pos[0] = Par[iC](i, 0, 0);
tmp_pos[1] = Par[iC](i, 1, 0);
tmp_pos[2] = Par[iC](i, 2, 0);
Expand All @@ -1587,20 +1595,21 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits,
}
}

if (count == 0) continue;
if (done_count == N_proc) break;
if (here_count == 0) continue;

// ZZZ Could add missing hits here, only if there are any actual matches.

if (LI.is_barrel())
{
PropagateTracksToHitR(msPar, N_proc, Config::backward_fit_pflags);
PropagateTracksToHitR(msPar, N_proc, Config::backward_fit_pflags, &no_mat_effs);

kalmanOperation(KFO_Calculate_Chi2 | KFO_Update_Params,
Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc);
}
else
{
PropagateTracksToHitZ(msPar, N_proc, Config::backward_fit_pflags);
PropagateTracksToHitZ(msPar, N_proc, Config::backward_fit_pflags, &no_mat_effs);

kalmanOperationEndcap(KFO_Calculate_Chi2 | KFO_Update_Params,
Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc);
Expand Down
44 changes: 26 additions & 18 deletions mkFit/PropagationMPlex.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -407,8 +407,9 @@ void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, co

void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar,
const MPlexQI &inChg, const MPlexQF& msRad,
MPlexLS &outErr, MPlexLV& outPar,
const int N_proc, const PropagationFlags pflags)
MPlexLS &outErr, MPlexLV& outPar,
const int N_proc, const PropagationFlags pflags,
const MPlexQI *noMatEffPtr)
{
// bool debug = true;

Expand Down Expand Up @@ -448,14 +449,17 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar,
MPlexQF hitsXi;
MPlexQF propSign;
#pragma omp simd
for (int n = 0; n < NN; ++n)
for (int n = 0; n < N_proc; ++n)
{
const int zbin = getZbinME(outPar(n, 2, 0));
const int rbin = getRbinME(msRad (n, 0, 0));

hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations
hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations

if (failFlag(n, 0, 0) || (noMatEffPtr && noMatEffPtr->ConstAt(n, 0, 0))) {
hitsRl(n, 0, 0) = 0.f;
hitsXi(n, 0, 0) = 0.f;
} else {
const int zbin = getZbinME(outPar(n, 2, 0));
const int rbin = getRbinME(msRad (n, 0, 0));
hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations
hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations
}
const float r0 = hipo(inPar(n, 0, 0), inPar(n, 1, 0));
const float r = msRad(n, 0, 0);
propSign(n, 0, 0) = (r>r0 ? 1. : -1.);
Expand Down Expand Up @@ -505,8 +509,9 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar,

void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar,
const MPlexQI &inChg, const MPlexQF& msZ,
MPlexLS &outErr, MPlexLV& outPar,
const int N_proc, const PropagationFlags pflags)
MPlexLS &outErr, MPlexLV& outPar,
const int N_proc, const PropagationFlags pflags,
const MPlexQI *noMatEffPtr)
{
// debug = true;

Expand Down Expand Up @@ -541,14 +546,17 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar,
MPlexQF hitsXi;
MPlexQF propSign;
#pragma omp simd
for (int n = 0; n < NN; ++n)
for (int n = 0; n < N_proc; ++n)
{
const int zbin = getZbinME(msZ(n, 0, 0));
const int rbin = getRbinME(std::hypot(outPar(n, 0, 0), outPar(n, 1, 0)));

hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations
hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations

if (noMatEffPtr && noMatEffPtr->ConstAt(n, 0, 0)) {
hitsRl(n, 0, 0) = 0.f;
hitsXi(n, 0, 0) = 0.f;
} else {
const int zbin = getZbinME(msZ(n, 0, 0));
const int rbin = getRbinME(std::hypot(outPar(n, 0, 0), outPar(n, 1, 0)));
hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations
hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations
}
const float zout = msZ.ConstAt(n, 0, 0);
const float zin = inPar.ConstAt(n, 2, 0);
propSign(n, 0, 0) = (std::abs(zout)>std::abs(zin) ? 1. : -1.);
Expand Down
6 changes: 4 additions & 2 deletions mkFit/PropagationMPlex.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,8 @@ void propagateLineToRMPlex(const MPlexLS &psErr, const MPlexLV& psPar,
void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar,
const MPlexQI &inChg, const MPlexQF& msRad,
MPlexLS &outErr, MPlexLV& outPar,
const int N_proc, const PropagationFlags pflags);
const int N_proc, const PropagationFlags pflags,
const MPlexQI *noMatEffPtr=nullptr);

void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad,
MPlexLV& outPar, MPlexLL& errorProp,
Expand All @@ -44,7 +45,8 @@ void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, const
void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar,
const MPlexQI &inChg, const MPlexQF& msZ,
MPlexLS &outErr, MPlexLV& outPar,
const int N_proc, const PropagationFlags pflags);
const int N_proc, const PropagationFlags pflags,
const MPlexQI *noMatEffPtr=nullptr);

void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ,
MPlexLV& outPar, MPlexLL& errorProp,
Expand Down
2 changes: 1 addition & 1 deletion mkFit/buildtestMPlex.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -715,7 +715,7 @@ void run_OneIteration(const TrackerInfo& trackerInfo, const IterationConfig &itc
builder.filter_comb_cands([&](const TrackCand &t)
{ return StdSeq::qfilter_nan_n_silly(t); });

builder.export_best_comb_cands(out_tracks);
builder.export_best_comb_cands(out_tracks, true);

if (do_remove_duplicates)
{
Expand Down