Skip to content
Snippets Groups Projects
Commit 32dc97be authored by Andre Gunther's avatar Andre Gunther :island:
Browse files

Merge branch 'wh_updatepvfinderforsmog' into 'master' 

Update of TBLV for smog reconstruction

Closes Moore#693

See merge request !3781
parents b45c912d 16c60ba2
No related branches found
No related tags found
1 merge request!3781Update of TBLV for smog reconstruction
Pipeline #7105327 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
Tr/PatPV/src/TrackBeamLineVertexFinderSoA.cpp
+ 249
− 142
View file @ 32dc97be
......@@ -94,17 +94,24 @@ public:
private:
Gaudi::Property<uint16_t> m_minNumTracksPerVertex{this, "MinNumTracksPerVertex", defaultMinNumTracks,
"Min number of tracks per PV"};
Gaudi::Property<float> m_minTotalTrackWeightPerVertex{this, "MinTotalTrackWeightPerVertex", 0.,
"Min sum of Tukey weights per PV"};
Gaudi::Property<float> m_zmin{this, "MinZ", defaultMinZ, "Min z position of vertex seed"};
Gaudi::Property<float> m_zmax{this, "MaxZ", defaultMaxZ, "Max z position of vertex seed"};
Gaudi::Property<float> m_zminIR{this, "MinZIR", -300., "Min z position of vertex seed for Interaction Region"};
Gaudi::Property<float> m_zmaxIR{this, "MaxZIR", +300., "Max z position of vertex seed for Interaction Region"};
Gaudi::Property<float> m_dz{this, "ZBinSize", defaultZBinSize, "Z histogram bin size"};
Gaudi::Property<float> m_maxTrackZ0Err{this, "MaxTrackZ0Err", 1.5 * Gaudi::Units::mm,
Gaudi::Property<float> m_maxTrackZ0Err{this, "MaxTrackZ0Err", 10.0 * Gaudi::Units::mm,
"Maximum z0-error for adding track to histo"};
Gaudi::Property<float> m_minDensity{this, "MinDensity", 0. / Gaudi::Units::mm,
"Minimal density at cluster peak (inverse resolution)"};
Gaudi::Property<float> m_minDipDensity{this, "MinDipDensity", 3.0 / Gaudi::Units::mm,
Gaudi::Property<float> m_maxTrackZ0ErrIR{this, "MaxTrackZ0ErrIR", 1.5 * Gaudi::Units::mm,
"Maximum z0-error for adding track to histo for Interaction Region"};
Gaudi::Property<float> m_Z0ErrMax{this, "Z0ErrMax", 1.0 * Gaudi::Units::mm, "Maximum for z0-error binning"};
Gaudi::Property<float> m_minBinContent{this, "MinBinContent", 0.01, "Minimal integral in bin for cluster boundary"};
Gaudi::Property<float> m_minDipSize{this, "MinDipSize", 1.0, "Minimal track integral on one side of histogram peak"};
Gaudi::Property<float> m_minDipDensity{this, "MinDipDensity", 0.0 / Gaudi::Units::mm,
"Minimal depth of a dip to split cluster (inverse resolution)"};
Gaudi::Property<float> m_minTracksInSeed{this, "MinTrackIntegralInSeed", 2.5};
Gaudi::Property<float> m_maxVertexRho{this, "BeamSpotRCut", 0.3 * Gaudi::Units::mm,
Gaudi::Property<float> m_minTracksInSeed{this, "MinTrackIntegralInSeed", 2.9};
Gaudi::Property<float> m_maxVertexRho{this, "BeamSpotRCut", 0.3 * Gaudi::Units::mm,
"Maximum distance of vertex to beam line"};
Gaudi::Property<uint16_t> m_maxFitIter{this, "MaxFitIter", defaultMaxFitIter,
"Maximum number of iterations for vertex fit"};
......@@ -116,19 +123,19 @@ private:
"Limit on change in chi2 to determine if vertex fit has converged"};
Gaudi::Property<float> m_minVertexZSeparationChi2{this, "MinVertexZSeparationChi2", defaultMinVertexZSeparationChi2};
Gaudi::Property<float> m_minVertexZSeparation{this, "MinVertexZSeparation", defaultMinVertexZSeparation};
Gaudi::Property<float> m_maxTrackBLChi2{this, "MaxTrackBLChi2", 10,
Gaudi::Property<float> m_maxTrackBLChi2{this, "MaxTrackBLChi2", 25.,
"Maximum chi2 of track to beam line contributing to seeds"};
Gaudi::Property<float> m_beamLineOffsetX{this, "BeamLineOffsetX", 0.0f, "X correction applied to beamline position"};
Gaudi::Property<float> m_beamLineOffsetY{this, "BeamLineOffsetY", 0.0f, "Y correction applied to beamline position"};
static constexpr uint16_t Nztemplatebins = 16;
static constexpr uint16_t Nztemplates = 32;
static constexpr uint16_t TemplateNormalization = 128;
uint16_t m_ztemplates[2 * Nztemplates][Nztemplatebins]; // odd and even, see note
static constexpr int Nztemplatebins = 16;
static constexpr int Nztemplates = 32;
static constexpr unsigned TemplateNormalization = 128;
using ZHistoType = uint32_t;
ZHistoType m_ztemplates[2 * Nztemplates][Nztemplatebins]; // odd and even, see note
mutable Gaudi::Accumulators::SummingCounter<unsigned int> m_nbPVsCounter{this, "Nb PVs"};
mutable Gaudi::Accumulators::BinomialCounter<unsigned int> m_frFailedFits{this, "Fraction of failed PV fits"};
mutable Gaudi::Accumulators::BinomialCounter<unsigned int> m_frUnconvergedFits{this,
"Fraction of unconverged PV fits"};
#ifdef TIMINGHISTOGRAMMING
AIDA::IProfile1D* m_timeperstepPr{nullptr};
AIDA::IProfile1D* m_timevsntrksPr{nullptr};
......@@ -157,6 +164,14 @@ StatusCode TrackBeamLineVertexFinderSoA::initialize() {
inputLocation<LHCb::Conditions::InteractionRegion>() );
} );
// Print a message if the number of bins in the histogram is too small for
// stack allocation (not sure it makes a real difference).
const int nZBins = ( m_zmax - m_zmin ) / m_dz;
if ( nZBins < defaultNZBins )
info() << "With this configuration, the number of bins in the zpoca histogram (" << nZBins
<< ") is smaller than used for static allocation (defaultNZB=" << defaultNZBins
<< ") . Please increase the latter in the C++." << endmsg;
// Fill the odd templates first
const double sqrthalf = std::sqrt( 0.5 );
{
......@@ -164,15 +179,21 @@ StatusCode TrackBeamLineVertexFinderSoA::initialize() {
const double zmaxodd = m_dz * ( Nztemplatebins / 2 - 1 + 0.5 );
const double zminodd = -zmaxodd;
for ( int itemplate = 0; itemplate < Nztemplates; ++itemplate ) {
const double sigmaz = m_maxTrackZ0Err * double( itemplate + 1 ) / Nztemplates;
const double sigmaz = m_Z0ErrMax * double( itemplate + 1 ) / Nztemplates;
double integral = 0.5 * std::erf( sqrthalf * zminodd / sigmaz );
const double offset = ( 1 + 2 * integral ) / ( Nztemplatebins - 1 );
for ( int ibin = 0; ibin < Nztemplatebins - 1; ++ibin ) {
double thisintegral = 0.5 * std::erf( sqrthalf * ( zminodd + ( ibin + 1 ) * m_dz ) / sigmaz );
double bincontent = thisintegral - integral;
m_ztemplates[2 * itemplate + 1][ibin] = int( bincontent * TemplateNormalization );
double bincontent = thisintegral - integral + offset;
m_ztemplates[2 * itemplate + 1][ibin] = std::lround( bincontent * TemplateNormalization );
integral = thisintegral;
}
// The last bin is empty because it isn't use for odd templates
m_ztemplates[2 * itemplate + 1][Nztemplatebins - 1] = 0;
debug() << "Odd template integral: " << itemplate << " "
<< std::accumulate( m_ztemplates[2 * itemplate + 1], m_ztemplates[2 * ( itemplate ) + 2], 0 ) /
float( TemplateNormalization )
<< endmsg;
}
}
......@@ -182,14 +203,19 @@ StatusCode TrackBeamLineVertexFinderSoA::initialize() {
const double zmaxeven = m_dz * Nztemplatebins / 2;
const double zmineven = -zmaxeven;
for ( int itemplate = 0; itemplate < Nztemplates; ++itemplate ) {
const double sigmaz = m_maxTrackZ0Err * double( itemplate + 1 ) / Nztemplates;
const double sigmaz = m_Z0ErrMax * double( itemplate + 1 ) / Nztemplates;
double integral = 0.5 * std::erf( sqrthalf * zmineven / sigmaz );
const double offset = ( 1 + 2 * integral ) / Nztemplatebins;
for ( int ibin = 0; ibin < Nztemplatebins; ++ibin ) {
double thisintegral = 0.5 * std::erf( sqrthalf * ( zmineven + ( ibin + 1 ) * m_dz ) / sigmaz );
double bincontent = thisintegral - integral;
m_ztemplates[2 * itemplate][ibin] = int( bincontent * TemplateNormalization );
double bincontent = thisintegral - integral + offset;
m_ztemplates[2 * itemplate][ibin] = std::lround( bincontent * TemplateNormalization );
integral = thisintegral;
}
debug() << "Even template integral: " << itemplate << " "
<< std::accumulate( m_ztemplates[2 * itemplate], m_ztemplates[2 * itemplate + 1], 0 ) /
float( TemplateNormalization )
<< endmsg;
}
}
......@@ -233,6 +259,70 @@ namespace {
uint16_t ntracks{0};
};
// Function that extracts clusters from a container of extrema. It
// finds the two largest maxima, then tries to partition on the
// smallest minimum in between. Selector is a function that takes a
// maximum, a minimum and another maximum. If the partitioning is
// successful, the container of extrema is split, and the function
// called on each container. If it is not successful, the second
// maximum, and all extrema in between the first and second maximum,
// are removed, and the function is called on the remaining
// container. The recursion stops of the container has three
// entries: that constitutes a cluster. At least one cluster is
// returned.
template <typename ExtremaContainer, typename ClusterContainer, typename DipSelector>
void clusterize( const ExtremaContainer& extrema, ClusterContainer& clusters, DipSelector dipselector ) {
// all odd extrema are maxima. find the two largest
auto const N = extrema.size();
assert( N % 2 == 1 );
if ( N < 3 ) {
return;
} else if ( N == 3 ) {
clusters.emplace_back( extrema[0].index, extrema[2].index, extrema[1].index );
} else {
// find the two largest extrema
unsigned largest = 1;
unsigned nextlargest = 3;
if ( extrema[largest].value < extrema[nextlargest].value ) std::swap( largest, nextlargest );
for ( unsigned i = 5; i < N; i += 2 ) {
if ( extrema[i].value > extrema[largest].value ) {
nextlargest = largest;
largest = i;
} else if ( extrema[i].value > extrema[nextlargest].value ) {
nextlargest = i;
}
}
// find the smallest minimum in between
const int delta = largest < nextlargest ? 1 : -1;
int smallest = largest + delta;
for ( unsigned i = largest + 3 * delta; i != nextlargest + delta; i += 2 * delta ) {
if ( extrema[smallest].value > extrema[i].value ) smallest = i;
}
// can we partition at this point?
if ( dipselector( extrema, std::min( largest, nextlargest ), smallest, std::max( largest, nextlargest ) ) ) {
// create two new extrema containers and recurcively call this function on them.
// Note that they overlap one element, which is okay.
ExtremaContainer c1{extrema.begin(), extrema.begin() + smallest + 1};
ExtremaContainer c2{extrema.begin() + smallest, extrema.end()};
clusterize( c1, clusters, dipselector );
clusterize( c2, clusters, dipselector );
} else {
// we cannot partition between the two maximu.
// create a new extrema container in which the second maximum and everything up to the first maximum are
// removed. then try again.
ExtremaContainer c1;
if ( largest < nextlargest ) {
c1.insert( c1.end(), extrema.begin(), extrema.begin() + largest + 1 );
c1.insert( c1.end(), extrema.begin() + nextlargest + 1, extrema.end() );
} else {
c1.insert( c1.end(), extrema.begin(), extrema.begin() + nextlargest );
c1.insert( c1.end(), extrema.begin() + largest, extrema.end() );
}
clusterize( c1, clusters, dipselector );
}
}
}
// inline std::ostream& operator<<( std::ostream& os, Cluster const& c ) {
// os << "[" << c.izfirst << ", " << c.izlast << ", " << c.izmax << "]";
// return os;
......@@ -354,50 +444,78 @@ PrimaryVertexContainer TrackBeamLineVertexFinderSoA::
// NOTE/TODO: As discussed in Rec#122, this could (and perhaps should) be drawn from an algorithm-local pool instead
// of the event-local pool used here. Alternatively, Wouter suggests that specific example could just be
// std::array, but this doesn't invalidate the algorithm-local idea!
const int nZBins = ( m_zmax - m_zmin ) / m_dz;
boost::container::small_vector<uint16_t, defaultNZBins> zhisto( nZBins, 0 );
const float halfwindow = ( Nztemplatebins / 2 + 1 ) * m_dz;
const float zmin = m_zmin + halfwindow;
const float zmax = m_zmax - halfwindow;
const int nZBins = ( m_zmax - m_zmin ) / m_dz;
boost::container::small_vector<ZHistoType, defaultNZBins + 1> zhisto( nZBins, 0 );
{
const float halfwindow = ( Nztemplatebins / 2 + 1 ) * m_dz;
const float zmin = m_zmin + halfwindow;
const float zmax = m_zmax - halfwindow;
const float zminIR = m_zminIR.value();
const float zmaxIR = m_zmaxIR.value();
int nacceptedtracks{0};
for ( const auto pvtrack : pvseedtracks.simd() ) {
auto const loop_mask = pvtrack.loop_mask();
const auto zbeam = pvtrack.get<PVTrackTag::z>();
const auto zerr = pvtrack.get<PVTrackTag::zerr>();
const auto blchi2 = pvtrack.get<PVTrackTag::blchi2>();
// bin in which z0 is, in floating point
auto const jbin =
int_v{4 * ( zbeam - m_zmin.value() ) / m_dz.value()}; // we need factor 4 to determine odd or even
int_v const dbin = jbin & 3; // & 3 is the same as % 4
auto const minbin =
to_std_array( ( jbin >> 2 ) - Nztemplatebins / 2 + select( dbin == 0, int_v( 0 ), int_v( 1 ) ) );
auto const oddtemplate = to_std_array( select( ( dbin == 0 ) || ( dbin == 3 ), int_v( 0 ), int_v( 1 ) ) );
//// make sure the template fits. make sure first and last bin
//// remain empty. there will be lot's of funny effects at edge of
//// histogram but we do not care.
//// compute the bealine chi2. we make a cut for tracks contributing to seeding
//// compute the beamline chi2. we make a cut for tracks contributing to seeding
// (bl_x, bl_y)^T * W * (bl_x, bl_y) for diagonal W
auto const zerrbin = int_v{Nztemplates / m_maxTrackZ0Err.value() * zerr};
auto const mask = to_std_array( int_v( float_v( zerrbin ) < Nztemplates && blchi2 < m_maxTrackBLChi2.value() &&
zmin < zbeam && zbeam < zmax && loop_mask ) );
auto const zerrbin_arr = to_std_array( zerrbin );
/// add selected entries to the histogram.
auto const mask =
to_std_array( int_v( loop_mask && blchi2 < m_maxTrackBLChi2.value() && zmin < zbeam && zbeam < zmax &&
zerr < m_maxTrackZ0Err.value() &&
( zminIR > zbeam || zmaxIR < zbeam || zerr < m_maxTrackZ0ErrIR.value() ) ) );
auto const zerrrebinfactor = int_v( 1 + zerr / m_Z0ErrMax.value() );
auto const zerrscale = float_v{1.0} / zerrrebinfactor;
// auto const zerrbin = int_v(Nztemplates / m_Z0ErrMax.value() * (zerr/zerrrebinfactor)) ;
auto const zerrbin = int_v( ( Nztemplates / m_Z0ErrMax.value() ) * ( zerr * zerrscale ) );
// bin in which z0 is, in floating point
auto const jbin =
int_v{4 * ( zbeam - m_zmin.value() ) / m_dz.value()}; // we need factor 4 to determine odd or even
int_v const dbin = jbin & 3; // & 3 is the same as % 4
auto const minbin = to_std_array( ( jbin >> 2 ) - zerrrebinfactor * ( Nztemplatebins / 2 ) +
select( dbin == 0, int_v( 0 ), int_v( 1 ) ) );
auto const oddtemplate = select( ( dbin == 0 ) || ( dbin == 3 ), int_v( 0 ), int_v( 1 ) );
auto const templatebin = to_std_array( 2 * zerrbin + oddtemplate );
auto const zerrscale_arr = to_std_array( zerrscale );
auto const zerrrebinfactor_arr = to_std_array( zerrrebinfactor );
for ( std::size_t simd_idx = 0; simd_idx < simd::size; ++simd_idx ) {
if ( mask[simd_idx] ) {
for ( int j = 0; j < Nztemplatebins; ++j ) {
zhisto[j + minbin[simd_idx]] += m_ztemplates[2 * zerrbin_arr[simd_idx] + oddtemplate[simd_idx]][j];
++nacceptedtracks;
// minbin is the bin in which the templates start (already corrected for scale)
const int beginbin = std::max( 0, -minbin[simd_idx] );
const int endbin = std::min( Nztemplatebins * zerrrebinfactor_arr[simd_idx], nZBins - minbin[simd_idx] - 1 );
// if there were no error scaling, this run would run from 0 to Nztemplatebins=16. now it runs from 0 to
// Nztemplatebins*rebinfactor. this has become much more timeconsuming perhaps because the compiler no longer
// unrolls the loop.
for ( int j = beginbin; j < endbin; ++j ) {
const auto weight =
m_ztemplates[templatebin[simd_idx]][j / zerrrebinfactor_arr[simd_idx]] * zerrscale_arr[simd_idx];
zhisto[j + minbin[simd_idx]] += std::lround( weight );
}
}
}
}
if ( msgLevel( MSG::DEBUG ) )
debug() << "Histogram integral: "
<< std::accumulate( zhisto.begin(), zhisto.end(), 0 ) / float( TemplateNormalization ) << " "
<< pvseedtracks.size() << " " << nacceptedtracks << endmsg;
}
#ifdef TIMINGHISTOGRAMMING
timer[2].stop();
timer[3].start();
#endif
// info() << "Start of cluster search" << endmsg ;
////
////// Step 3: perform a peak search in the histogram. This used to be
////// very simple but the logic needed to find 'significant dips' made
......@@ -406,31 +524,35 @@ PrimaryVertexContainer TrackBeamLineVertexFinderSoA::
////
boost::container::small_vector<Cluster, 32> clusters;
{
// step A: make 'ProtoClusters'
// Step B: for each such ProtoClusters
// - find the significant extrema (an odd number, start with a minimum. you can always achieve this by adding a
// zero bin at the beginning)
// an extremum is a bin-index, plus the integral till that point, plus the content of the bin
// - find the highest extremum and
// - try and partition at the lowest minimum besides it
// - if that doesn't work, try the other extremum
// - if that doesn't work, accept as cluster
// Step A: make 'proto-clusters': these are subsequent bins with non-zero content and an integral above the
// threshold.
const uint32_t minTracksInSeed = m_minTracksInSeed * TemplateNormalization;
const uint32_t mindip = m_minDipDensity * m_dz * TemplateNormalization; // need to invent something
const uint32_t minpeak = m_minDensity * m_dz * TemplateNormalization;
using BinIndex = uint16_t;
// FIXME: how dangerous is it to use a fixed capacity vector?
// Step A: make 'ProtoClusters': these are subsequent bins with
// non-zero content and an integral above the threshold.
//
// Step B: for each such ProtoClusters create a list of 'extrema':
// these are points where the derivative changes sign. an
// extremum is a bin-index, plus the integral till that point,
// plus the content of the bin.
//
// Step C: using the extrema, find 'partition points': a partition
// point is the smallest minimum between the two largest
// maxima. if the partition point is valid, then use it, otherwise
// delete the smallest of the two maximum and try with the next
// pair. If no more partition points are found, accept the
// cluster.
// Step A: make 'proto-clusters'
const unsigned minTracksInSeed = m_minTracksInSeed * TemplateNormalization;
const unsigned mindipsize = m_minDipSize * TemplateNormalization;
const unsigned mindip = m_minDipDensity * m_dz * TemplateNormalization; // need to invent something
const unsigned threshold = std::max( int( m_minBinContent * TemplateNormalization ), 1 );
using BinIndex = unsigned;
boost::container::small_vector<BinIndex, 64> clusteredges;
{
bool prevempty = true;
uint32_t integral = zhisto[0];
for ( BinIndex i = 1; i < zhisto.size(); ++i ) {
integral += zhisto[i];
bool empty = zhisto[i] < 1;
bool empty = !( zhisto[i] >= threshold );
if ( empty != prevempty ) {
if ( empty ) {
if ( integral >= minTracksInSeed )
......@@ -445,91 +567,78 @@ PrimaryVertexContainer TrackBeamLineVertexFinderSoA::
}
}
}
// Step B: turn these into clusters. There can be more than one cluster per proto-cluster.
// Step B: find the extrema. All even extrema are maxima.
const size_t Nproto = clusteredges.size() / 2;
for ( uint16_t i = 0; i < Nproto; ++i ) {
const BinIndex ibegin = clusteredges[i * 2];
const BinIndex iend = clusteredges[i * 2 + 1];
// find the extrema. all this complicated logic is to be able to
// split close seeds.
boost::container::small_vector<Extremum, 32> extrema;
{
bool rising = true;
uint32_t integral = zhisto[ibegin];
extrema.emplace_back( ibegin, zhisto[ibegin], integral );
for ( uint16_t i = ibegin; i < iend; ++i ) {
const auto value = zhisto[i];
if ( value != zhisto[i + 1] ) {
bool stillrising = zhisto[i + 1] > value;
if ( stillrising != rising ) { // found a local extremum
const uint8_t n = extrema.size();
const int16_t deltavalue = value - extrema.back().value;
if ( rising ) { // found a local maximum
if ( n % 2 == 1 ) { // the last one was a minimum. check that this maximum is high enough. we always
// accept the first maximum.
if ( ( n == 1 && value >= minpeak ) || deltavalue >= int( mindip ) )
extrema.emplace_back( i, value, integral + value / 2 );
} else { // the last one was a maximum, but apparently there were no good minima in between
if ( deltavalue > 0 ) {
extrema.pop_back();
extrema.emplace_back( i, value, integral + value / 2 );
}
}
} else { // found a local minimum
if ( n % 2 == 0 ) { // the last one was a maximum. check that this minimum is small enough
if ( -1 * deltavalue >= int( mindip ) ) extrema.emplace_back( i, value, integral + 0.5f * value );
} else { // the last one was a minimum, but apparently there were no good maxima in between
if ( deltavalue < 0 ) {
extrema.pop_back();
extrema.emplace_back( i, value, integral + value / 2 );
}
}
}
}
rising = stillrising;
for ( uint16_t j = 0; j < Nproto; ++j ) {
const auto ibegin = clusteredges[j * 2];
const auto iend = clusteredges[j * 2 + 1];
// Find the extrema. All odd extrema are minima. All even extrema are maxima.
boost::container::small_vector<Extremum, 64> extrema;
uint32_t integral = zhisto[ibegin];
extrema.emplace_back( ibegin, 0, integral / 2 );
for ( BinIndex i = ibegin; i < iend; ++i ) {
const auto thisvalue = zhisto[i];
const auto thisintegral = integral + thisvalue / 2;
integral += thisvalue;
// Four cases:
// A. last extremum was a maximum and we are larger: replace the last extremum
// B. last extremum was a mimimum and we are smaller: replace the last extremum
// C. last extremum was a maximum and we are smaller: accept if sufficiently smaller.
// D. last extremum was a minimum and we are larger: accept if sufficiently larger.
if ( extrema.size() % 2 == 0 ) { // case A or C
if ( thisvalue > extrema.back().value ) { // case A
extrema.back().index = i;
extrema.back().value = thisvalue;
extrema.back().integral = thisintegral;
} else { // case C
if ( thisvalue < zhisto[i + 1] && thisvalue + threshold < extrema.back().value )
extrema.emplace_back( i, thisvalue, thisintegral );
}
integral += value;
}
if ( extrema.size() % 2 == 1 ) { // last was minimum. this one should replace it
extrema.pop_back();
}
extrema.emplace_back( iend, zhisto[iend], integral );
}
// FIXME: temporary logic check
// if( extrema.size()%2==0 ) {
// warning() << "ERROR: even number of extrema found." << extrema.size() << endmsg ;
//}
if ( extrema.size() >= 3 ) {
// now partition on extrema
const auto N = extrema.size();
boost::container::small_vector<Cluster, 16> subclusters;
if ( N > 3 ) {
for ( uint32_t i = 1; i < N / 2 + 1; ++i ) {
if ( extrema[2 * i].integral - extrema[2 * i - 2].integral > minTracksInSeed ) {
subclusters.emplace_back( extrema[2 * i - 2].index, extrema[2 * i].index, extrema[2 * i - 1].index );
}
} else { // Case B or D
if ( thisvalue <= extrema.back().value ) { // case B
extrema.back().index = i;
extrema.back().value = thisvalue;
extrema.back().integral = thisintegral;
} else { // Case D
if ( thisvalue > zhisto[i + 1] && thisvalue > extrema.back().value + threshold )
extrema.emplace_back( i, thisvalue, thisintegral );
}
}
if ( subclusters.empty() ) {
clusters.emplace_back( extrema.front().index, extrema.back().index, extrema[1].index );
} else {
// adjust the limit of the first and last to extend to the entire protocluster
subclusters.front().izfirst = ibegin;
subclusters.back().izlast = iend;
clusters.insert( std::end( clusters ), std::begin( subclusters ), std::end( subclusters ) );
}
}
// replace the last minimum, or add a minimum
const auto N = extrema.size();
if ( N % 2 == 1 ) extrema.pop_back();
extrema.emplace_back( iend, 0, integral );
// Step C: call the clusterization
clusterize( extrema, clusters,
[mindip, minTracksInSeed, mindipsize]( const auto& v, int posmax1, int posmin, int posmax2 ) {
// this functions returns true if the minimum between these maxima is an acceptable partition point.
// check that there are enough tracks on either side of the minimum
if ( v[posmin].integral - v.front().integral < minTracksInSeed ) return false;
if ( v.back().integral - v[posmin].integral < minTracksInSeed ) return false;
// check that the maxima are large enough
if ( v[posmax1].value < v[posmin].value + mindip ) return false;
if ( v[posmax2].value < v[posmin].value + mindip ) return false;
if ( mindipsize > 0 ) {
if ( ( v[posmin].integral - v[posmax1].integral ) <
v[posmin].value * ( v[posmin].index - v[posmax1].index ) + mindipsize )
return false;
if ( ( v[posmax2].integral - v[posmin].integral ) <
v[posmin].value * ( v[posmax2].index - v[posmin].index ) + mindipsize )
return false;
}
return true;
} );
}
if ( msgLevel( MSG::DEBUG ) ) debug() << "Number of clusters: " << clusters.size() << endmsg;
}
// std::cout << "Number of clusters: " << clusters.size() << std::endl ;
#ifdef TIMINGHISTOGRAMMING
timer[3].stop();
#endif
// FIXME: we set up a lot of navigation below which is difficult to
// fix if a PV is removed in the final selection step. Luckily this
// happens so rarely that we can choose a lazy solution: if a PV is
......@@ -587,27 +696,26 @@ PrimaryVertexContainer TrackBeamLineVertexFinderSoA::
clusteringready = smallestcluster->ntracks >= m_minNumTracksPerVertex;
if ( !clusteringready ) clusters.erase( smallestcluster );
}
#ifdef TIMINGHISTOGRAMMING
timer[4].stop();
#endif
}
//// Step 5: Seed the primary vertices and assign tracks
// Now we need to construct the output
#ifdef TIMINGHISTOGRAMMING
timer[4].stop();
#endif
if ( !clusters.empty() ) {
#ifdef TIMINGHISTOGRAMMING
timer[5].start();
#endif
// std::vector<Vertex, LHCb::Allocators::EventLocal<Vertex>> vertices{memResource};
// First initialize the vertex parameters. I found that this funny
// weighted 'maximum' is better than most other inexpensive
// solutions.
auto zClusterMean = [this, &zhisto]( auto izmax ) -> float {
const uint16_t* b = zhisto.data() + izmax;
int d1 = *b - *( b - 1 );
int d2 = *b - *( b + 1 );
float idz = d1 + d2 > 0 ? ( 0.5f * ( d1 - d2 ) ) / ( d1 + d2 ) : 0.0f;
const auto* b = zhisto.data() + izmax;
int d1 = *b - *( b - 1 );
int d2 = *b - *( b + 1 );
float idz = d1 + d2 > 0 ? ( 0.5f * ( d1 - d2 ) ) / ( d1 + d2 ) : 0.0f;
return m_zmin + m_dz * ( izmax + idz + 0.5f );
};
......@@ -647,13 +755,12 @@ PrimaryVertexContainer TrackBeamLineVertexFinderSoA::
// Remove any vertices that do not pass the selection. We assign their tracks to the next or previous vertex,
// depending on which one is closer.
applyBeamlineSelection( output, m_minNumTracksPerVertex, m_maxVertexRho, beamlineX, beamlineY, beamlineTx,
beamlineTy, fitconfig );
applyBeamlineSelection( output, m_minNumTracksPerVertex, m_minTotalTrackWeightPerVertex, m_maxVertexRho, beamlineX,
beamlineY, beamlineTx, beamlineTy, fitconfig );
#ifdef TIMINGHISTOGRAMMING
timer[7].stop();
#endif
}
// Set the keys of the PVs. Unfortunately, as long as PrimaryVertex
// derives from KeyedObject, we can only set the key once, so we
// really need to do this at the end. Perhaps we should get rid of
......
Tr/PatPV/src/TrackPVFinderUtils.h
+ 7
− 6
View file @ 32dc97be
......@@ -26,8 +26,8 @@ namespace LHCb::TrackPVFinder {
using int_v = simd::int_v;
/// Some default values of parameters used by different PV algorithms
constexpr float defaultMinZ = -300 * Gaudi::Units::mm;
constexpr float defaultMaxZ = +300 * Gaudi::Units::mm;
constexpr float defaultMinZ = -600 * Gaudi::Units::mm;
constexpr float defaultMaxZ = +600 * Gaudi::Units::mm;
constexpr float defaultMaxDeltaChi2 = 12;
constexpr float defaultMaxDeltaZConverged = 0.001 * Gaudi::Units::mm;
constexpr float defaultMaxDeltaChi2Converged = 0.01;
......@@ -148,15 +148,16 @@ namespace LHCb::TrackPVFinder {
}
inline void applyBeamlineSelection( PrimaryVertexContainer& pvcontainer, const size_t minNumTracksPerVertex,
const double maxVertexRho, const double beamlineX, const double beamlineY,
const double beamlineTx, const double beamlineTy,
const AdaptiveFitConfig fitconfig ) {
const double minTotalTrackWeightPerVertex, const double maxVertexRho,
const double beamlineX, const double beamlineY, const double beamlineTx,
const double beamlineTy, const AdaptiveFitConfig fitconfig ) {
applyPVSelection( pvcontainer, fitconfig, [=]( const auto& vertex ) {
const auto beamlinedx = vertex.position().x() - ( beamlineX + beamlineTx * vertex.position().z() );
const auto beamlinedy = vertex.position().y() - ( beamlineY + beamlineTy * vertex.position().z() );
const auto beamlinerho2 = sqr( beamlinedx ) + sqr( beamlinedy );
const auto maxVertexRho2 = sqr( maxVertexRho );
return vertex.nDoF() > 0 && vertex.nTracks() >= int( minNumTracksPerVertex ) && beamlinerho2 < maxVertexRho2;
return vertex.nDoF() > 0 && vertex.nTracks() >= int( minNumTracksPerVertex ) &&
vertex.sumOfWeights() >= minTotalTrackWeightPerVertex && beamlinerho2 < maxVertexRho2;
} );
}
......
Tr/PatPV/src/TrackUnbiasedPVFinderSoA.cpp
+ 4
− 2
View file @ 32dc97be
......@@ -65,6 +65,8 @@ public:
private:
Gaudi::Property<uint32_t> m_minNumTracksPerVertex{this, "MinNumTracksPerVertex", 5};
Gaudi::Property<uint16_t> m_minTotalTrackWeightPerVertex{this, "MinTotalTrackWeightPerVertex", 4.,
"Min sum of Tukey weights per PV"};
Gaudi::Property<float> m_zmin{this, "MinZ", defaultMinZ, "Min z position of vertex seed"};
Gaudi::Property<float> m_zmax{this, "MaxZ", defaultMaxZ, "Max z position of vertex seed"};
Gaudi::Property<float> m_maxTrackZ0Err{this, "MaxTrackZ0Err", 1.5 * Gaudi::Units::mm,
......@@ -688,8 +690,8 @@ PrimaryVertexContainer TrackUnbiasedPVFinderSoA::operator()( const EventContext&
mergeCloseVertices( output, m_minVertexZSeparation, m_minVertexZSeparationChi2, fitconfig );
// Remove vertices that are not compatible with the beamline
applyBeamlineSelection( output, m_minNumTracksPerVertex, m_maxVertexRho, beamlineX, beamlineY, beamlineTx,
beamlineTy, fitconfig );
applyBeamlineSelection( output, m_minNumTracksPerVertex, m_minTotalTrackWeightPerVertex, m_maxVertexRho, beamlineX,
beamlineY, beamlineTx, beamlineTy, fitconfig );
#ifdef TIMINGHISTOGRAMMING
timer[6].stop();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment