From 83bba10c805711238f79837cc327fe4656be5763 Mon Sep 17 00:00:00 2001 From: Lindsey Gray Date: Wed, 26 Jul 2017 19:04:30 +0200 Subject: [PATCH 01/16] vectorized 2D simulated annealing vertexing --- .../Configuration/python/RecoVertex_cff.py | 6 +- .../interface/DAClusterizerInZT_vect.h | 243 +++++ .../interface/PrimaryVertexProducer.h | 2 + .../plugins/PrimaryVertexProducer.cc | 3 + .../python/TkClusParameters_cff.py | 16 + .../src/DAClusterizerInZT.cc | 2 +- .../src/DAClusterizerInZT_vect.cc | 959 ++++++++++++++++++ 7 files changed, 1227 insertions(+), 4 deletions(-) create mode 100644 RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h create mode 100644 RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc diff --git a/RecoVertex/Configuration/python/RecoVertex_cff.py b/RecoVertex/Configuration/python/RecoVertex_cff.py index 0c117617b2a13..578cb10d69d50 100644 --- a/RecoVertex/Configuration/python/RecoVertex_cff.py +++ b/RecoVertex/Configuration/python/RecoVertex_cff.py @@ -33,14 +33,14 @@ ) #timing -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA2DParameters +from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA2D_vectParameters unsortedOfflinePrimaryVertices1D = unsortedOfflinePrimaryVertices.clone() unsortedOfflinePrimaryVertices1D.TkFilterParameters.minPt = cms.double(0.7) offlinePrimaryVertices1D=sortedPrimaryVertices.clone(vertices="unsortedOfflinePrimaryVertices1D", particles="trackRefsForJetsBeforeSorting") offlinePrimaryVertices1DWithBS=sortedPrimaryVertices.clone(vertices="unsortedOfflinePrimaryVertices1D:WithBS", particles="trackRefsForJetsBeforeSorting") -DA2DParameters.TkDAClusParameters.verbose = cms.untracked.bool(False) +DA2D_vectParameters.TkDAClusParameters.verbose = cms.untracked.bool(False) unsortedOfflinePrimaryVertices4D = unsortedOfflinePrimaryVertices.clone( verbose = cms.untracked.bool(False), - TkClusParameters = DA2DParameters ) + TkClusParameters = DA2D_vectParameters ) unsortedOfflinePrimaryVertices4D.TkFilterParameters.minPt = cms.double(0.7) unsortedOfflinePrimaryVertices4D.TrackTimesLabel = cms.InputTag("trackTimeValueMapProducer:generalTracksConfigurableFlatResolutionModel") unsortedOfflinePrimaryVertices4D.TrackTimeResosLabel = cms.InputTag("trackTimeValueMapProducer:generalTracksConfigurableFlatResolutionModelResolution") diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h new file mode 100644 index 0000000000000..196fd88321d97 --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -0,0 +1,243 @@ +#ifndef DAClusterizerInZT_vect_h +#define DAClusterizerInZT_vect_h + +/**\class DAClusterizerInZT_vect + + Description: separates event tracks into clusters along the beam line + + Version which auto-vectorizes with gcc 4.6 or newer + + */ + +#include "RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h" +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include +#include "DataFormats/Math/interface/Error.h" +#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h" +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" + + +class DAClusterizerInZT_vect final : public TrackClusterizerInZ { + +public: + // Internal data structure to + struct track_t { + + void AddItem( double new_z, double new_t, double new_dz2, double new_dt2, const reco::TransientTrack* new_tt, double new_pi ) + { + z.push_back( new_z ); + t.push_back( new_t ); + dz2.push_back( new_dz2 ); + dt2.push_back( new_dt2 ); + errsum.push_back( 1./(1./new_dz2 + 1./new_dt2) ); + tt.push_back( new_tt ); + + pi.push_back( new_pi ); // track weight + Z_sum.push_back( 1.0 ); // Z[i] for DA clustering, initial value as done in ::fill + } + + + + unsigned int GetSize() const + { + return z.size(); + } + + + // has to be called everytime the items are modified + void ExtractRaw() + { + _z = &z.front(); + _t = &t.front(); + _dz2 = &dz2.front(); + _dt2 = &dt2.front(); + _errsum = &errsum.front(); + _Z_sum = &Z_sum.front(); + _pi = &pi.front(); + } + + double * __restrict__ _z; // z-coordinate at point of closest approach to the beamline + double * __restrict__ _t; // t-coordinate at point of closest approach to the beamline + double * __restrict__ _dz2; // square of the error of z(pca) + double * __restrict__ _dt2; // square of the error of t(pca) + double * __restrict__ _errsum; // sum of squares of the pca errors + + double * __restrict__ _Z_sum; // Z[i] for DA clustering + double * __restrict__ _pi; // track weight + + std::vector z; // z-coordinate at point of closest approach to the beamline + std::vector t; // t-coordinate at point of closest approach to the beamline + std::vector dz2; // square of the error of z(pca) + std::vector dt2; // square of the error of t(pca) + std::vector errsum; // sum of squares of the pca errors + std::vector< const reco::TransientTrack* > tt; // a pointer to the Transient Track + + std::vector Z_sum; // Z[i] for DA clustering + std::vector pi; // track weight + }; + + struct vertex_t { + std::vector z; // z coordinate + std::vector t; // t coordinate + std::vector pk; // vertex weight for "constrained" clustering + + // --- temporary numbers, used during update + std::vector ei_cache; + std::vector ei; + std::vector sw; + std::vector swz; + std::vector swt; + std::vector se; + std::vector swE; + + + unsigned int GetSize() const + { + return z.size(); + } + + void AddItem( double new_z, double new_t, double new_pk ) + { + z.push_back( new_z); + t.push_back( new_t); + pk.push_back( new_pk); + + ei_cache.push_back( 0.0 ); + ei.push_back( 0.0 ); + sw.push_back( 0.0 ); + swz.push_back( 0.0); + swt.push_back( 0.0); + se.push_back( 0.0); + swE.push_back( 0.0); + + ExtractRaw(); + } + + void InsertItem( unsigned int i, double new_z, double new_t, double new_pk ) + { + z.insert(z.begin() + i, new_z); + t.insert(t.begin() + i, new_t); + pk.insert(pk.begin() + i, new_pk); + + ei_cache.insert(ei_cache.begin() + i, 0.0 ); + ei.insert( ei.begin() + i, 0.0 ); + sw.insert( sw.begin() + i, 0.0 ); + swz.insert(swz.begin() + i, 0.0 ); + swt.insert(swt.begin() + i, 0.0 ); + se.insert( se.begin() + i, 0.0 ); + swE.insert(swE.begin() + i, 0.0 ); + + ExtractRaw(); + } + + void RemoveItem( unsigned int i ) + { + z.erase( z.begin() + i ); + t.erase( t.begin() + i ); + pk.erase( pk.begin() + i ); + + ei_cache.erase( ei_cache.begin() + i); + ei.erase( ei.begin() + i); + sw.erase( sw.begin() + i); + swz.erase( swz.begin() + i); + swt.erase( swt.begin() + i); + se.erase(se.begin() + i); + swE.erase(swE.begin() + i); + + ExtractRaw(); + } + + void DebugOut() + { + std::cout << "vertex_t size: " << GetSize() << std::endl; + + for ( unsigned int i =0; i < GetSize(); ++ i) + { + std::cout << " z = " << _z[i] << " t = " << _t[i] << " pk = " << _pk[i] << std::endl; + } + } + + // has to be called everytime the items are modified + void ExtractRaw() + { + _z = &z.front(); + _t = &t.front(); + _pk = &pk.front(); + + _ei = &ei.front(); + _sw = &sw.front(); + _swz = &swz.front(); + _swt = &swt.front(); + _se = &se.front(); + _swE = &swE.front(); + _ei_cache = &ei_cache.front(); + + } + + double * __restrict__ _z; + double * __restrict__ _t; + double * __restrict__ _pk; + + double * __restrict__ _ei_cache; + double * __restrict__ _ei; + double * __restrict__ _sw; + double * __restrict__ _swz; + double * __restrict__ _swt; + double * __restrict__ _se; + double * __restrict__ _swE; + + }; + + DAClusterizerInZT_vect(const edm::ParameterSet& conf); + + + std::vector > + clusterize(const std::vector & tracks) const; + + + std::vector + vertices(const std::vector & tracks, + const int verbosity = 0) const ; + + track_t fill(const std::vector & tracks) const; + + double update(double beta, track_t & gtracks, + vertex_t & gvertices, bool useRho0, const double & rho0) const; + + void dump(const double beta, const vertex_t & y, + const track_t & tks, const int verbosity = 0) const; + bool merge(vertex_t & y, double & beta)const; + bool purge(vertex_t &, track_t &, double &, + const double) const; + void splitAll( vertex_t & y) const; + bool split(const double beta, track_t &t, vertex_t & y, double threshold = 1. ) const; + + double beta0(const double betamax, track_t const & tks, vertex_t const & y) const; + + +private: + bool verbose_; + double zdumpcenter_; + double zdumpwidth_; + + double vertexSize_; + int maxIterations_; + double coolingFactor_; + double betamax_; + double betastop_; + double dzCutOff_; + double d0CutOff_; + bool useTc_; + + double mintrkweight_; + double uniquetrkweight_; + double zmerge_; + double tmerge_; + double betapurge_; + +}; + + +//#ifndef DAClusterizerInZT_new_h +#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h index bafb5b6821b7c..e55489e1b8f84 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h @@ -34,6 +34,8 @@ #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFindingBase.h" #include "RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h" #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h" +#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h" + #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h" #include "RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h" diff --git a/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc index 3aef3af31403d..12cb6c8f9310b 100644 --- a/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc +++ b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc @@ -53,6 +53,9 @@ PrimaryVertexProducer::PrimaryVertexProducer(const edm::ParameterSet& conf) else if( clusteringAlgorithm=="DA2D" ) { theTrackClusterizer = new DAClusterizerInZT(conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); f4D = true; + } else if( clusteringAlgorithm=="DA2D_vect" ) { + theTrackClusterizer = new DAClusterizerInZT_vect(conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); + f4D = true; } else{ diff --git a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py index 10afac03823ca..8965c98dac5c4 100644 --- a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py +++ b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py @@ -25,3 +25,19 @@ uniquetrkweight = cms.double(0.8) # require at least two tracks with this weight at T=Tpurge ) ) + +DA2D_vectParameters = cms.PSet( + algorithm = cms.string("DA2D_vect"), + TkDAClusParameters = cms.PSet( + coolingFactor = cms.double(0.6), # moderate annealing speed + Tmin = cms.double(2.0), # end of vertex splitting + Tpurge = cms.double(2.0), # cleaning + Tstop = cms.double(0.5), # end of annealing + vertexSize = cms.double(0.006), # added in quadrature to track-z resolutions + d0CutOff = cms.double(3.), # downweight high IP tracks + dzCutOff = cms.double(3.), # outlier rejection after freeze-out (T0){ z1=z1/w1; t1=t1/w1;} else{ z1=y[ik].z-epsilon; t1=y[ik].t-epsilon; } - if(w2>0){ z2=z2/w2; t2=t2/w2;} else{ z2=y[ik].z+epsilon; t2=y[ik].t+epsilon;} + if(w2>0){ z2=z2/w2; t2=t2/w2;} else{ z2=y[ik].z+epsilon; t2=y[ik].t+epsilon; } // reduce split size if there is not enough room if( ( ik > 0 ) && ( y[ik-1].z>=z1 ) ){ z1=0.5*(y[ik].z+y[ik-1].z); t1=0.5*(y[ik].t+y[ik-1].t); } diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc new file mode 100644 index 0000000000000..c759849a4c30d --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -0,0 +1,959 @@ +#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h" +#include "RecoVertex/VertexPrimitives/interface/VertexException.h" + +#include +#include +#include +#include +#include "FWCore/Utilities/interface/isFinite.h" +#include "vdt/vdtMath.h" + +using namespace std; + +namespace { + constexpr double epsilon = 1.0e-3; + constexpr double vertexSizeTime = 0.008; + constexpr double dtCutOff = 4.0; +} + +DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) { + + // hardcoded parameters + maxIterations_ = 100; + mintrkweight_ = 0.5; // conf.getParameter("mintrkweight"); + + + // configurable debug outptut debug output + verbose_ = conf.getUntrackedParameter ("verbose", false); + zdumpcenter_ = conf.getUntrackedParameter ("zdumpcenter", 0.); + zdumpwidth_ = conf.getUntrackedParameter ("zdumpwidth", 20.); + + // configurable parameters + double Tmin = conf.getParameter ("Tmin"); + double Tpurge = conf.getParameter ("Tpurge"); + double Tstop = conf.getParameter ("Tstop"); + vertexSize_ = conf.getParameter ("vertexSize"); + coolingFactor_ = conf.getParameter ("coolingFactor"); + useTc_=true; + if(coolingFactor_<0){ + coolingFactor_=-coolingFactor_; + useTc_=false; + } + d0CutOff_ = conf.getParameter ("d0CutOff"); + dzCutOff_ = conf.getParameter ("dzCutOff"); + uniquetrkweight_ = conf.getParameter("uniquetrkweight"); + zmerge_ = conf.getParameter("zmerge"); + tmerge_ = conf.getParameter("tmerge"); + + if(verbose_){ + std::cout << "DAClusterizerinZT_vect: mintrkweight = " << mintrkweight_ << std::endl; + std::cout << "DAClusterizerinZT_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl; + std::cout << "DAClusterizerinZT_vect: zmerge = " << zmerge_ << std::endl; + std::cout << "DAClusterizerinZT_vect: tmerge = " << tmerge_ << std::endl; + std::cout << "DAClusterizerinZT_vect: Tmin = " << Tmin << std::endl; + std::cout << "DAClusterizerinZT_vect: Tpurge = " << Tpurge << std::endl; + std::cout << "DAClusterizerinZT_vect: Tstop = " << Tstop << std::endl; + std::cout << "DAClusterizerinZT_vect: vertexSize = " << vertexSize_ << std::endl; + std::cout << "DAClusterizerinZT_vect: coolingFactor = " << coolingFactor_ << std::endl; + std::cout << "DAClusterizerinZT_vect: d0CutOff = " << d0CutOff_ << std::endl; + std::cout << "DAClusterizerinZT_vect: dzCutOff = " << dzCutOff_ << std::endl; + std::cout << "DAClusterizerinZT_vect: dtCutoff = " << dtCutOff << std::endl; + } + + + if (Tmin == 0) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tmin" << Tmin + << " reset do default " << 1. / betamax_; + } else { + betamax_ = 1. / Tmin; + } + + + if ((Tpurge > Tmin) || (Tpurge == 0)) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tpurge" << Tpurge + << " set to " << Tmin; + Tpurge = Tmin; + } + betapurge_ = 1./Tpurge; + + + if ((Tstop > Tpurge) || (Tstop == 0)) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tstop" << Tstop + << " set to " << max(1., Tpurge); + Tstop = max(1., Tpurge) ; + } + betastop_ = 1./Tstop; + +} + + +namespace { + inline double local_exp( double const& inp) { + return vdt::fast_exp( inp ); + } + + inline void local_exp_list( double const * __restrict__ arg_inp, double * __restrict__ arg_out, const int arg_arr_size) { + for(int i=0; i!=arg_arr_size; ++i ) arg_out[i]=vdt::fast_exp(arg_inp[i]); + } + +} + +//todo: use r-value possibility of c++11 here +DAClusterizerInZT_vect::track_t +DAClusterizerInZT_vect::fill(const vector & tracks) const { + + // prepare track data for clustering + track_t tks; + for (auto it = tracks.begin(); it!= tracks.end(); it++){ + if (!(*it).isValid()) continue; + double t_pi=1.; + double t_z = ((*it).stateAtBeamLine().trackStateAtPCA()).position().z(); + double t_t = it->timeExt(); + if (std::fabs(t_z) > 1000.) continue; + auto const & t_mom = (*it).stateAtBeamLine().trackStateAtPCA().momentum(); + // get the beam-spot + reco::BeamSpot beamspot = (it->stateAtBeamLine()).beamSpot(); + double t_dz2 = + std::pow((*it).track().dzError(), 2) // track errror + + (std::pow(beamspot.BeamWidthX()*t_mom.x(),2)+std::pow(beamspot.BeamWidthY()*t_mom.y(),2))*std::pow(t_mom.z(),2)/std::pow(t_mom.perp2(),2) // beam spot width + + std::pow(vertexSize_, 2); // intrinsic vertex size, safer for outliers and short lived decays + t_dz2 = 1./ t_dz2; + double t_dt2 =std::pow((*it).dtErrorExt(),2.) + std::pow(vertexSizeTime,2.); // the ~injected~ timing error, need to add a small minimum vertex size in time + t_dt2 = 1./t_dt2; + if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits::min() ) continue; + if (d0CutOff_ > 0) { + Measurement1D atIP = + (*it).stateAtBeamLine().transverseImpactParameter();// error contains beamspot + t_pi = 1. / (1. + local_exp(std::pow(atIP.value() / atIP.error(), 2) - std::pow(d0CutOff_, 2))); // reduce weight for high ip tracks + if (edm::isNotFinite(t_pi) || t_pi < std::numeric_limits::epsilon()) continue; // usually is > 0.99 + } + LogTrace("DAClusterizerinZT_vectorized") << t_z << ' ' << t_t <<' '<< t_dz2 << ' ' << t_dt2 <<' '<< t_pi; + tks.AddItem(t_z, t_t, t_dz2, t_dt2, &(*it), t_pi); + } + tks.ExtractRaw(); + + if (verbose_) { + std::cout << "Track count " << tks.GetSize() << std::endl; + } + + return tks; +} + + +namespace { + inline + double Eik(double t_z, double k_z, double t_dz2, double t_t, double k_t, double t_dt2) { + return std::pow(t_z - k_z, 2) * t_dz2 + std::pow(t_t - k_t,2) * t_dt2; + } +} + +double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, + vertex_t & gvertices, bool useRho0, const double & rho0) const { + + //update weights and vertex positions + // mass constrained annealing without noise + // returns the squared sum of changes of vertex positions + + const unsigned int nt = gtracks.GetSize(); + const unsigned int nv = gvertices.GetSize(); + + //initialize sums + double sumpi = 0; + + // to return how much the prototype moved + double delta = 0; + + + // intial value of a sum + double Z_init = 0; + // independpent of loop + if ( useRho0 ) + { + Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); // cut-off + } + + // define kernels + auto kernel_calc_exp_arg = [ beta, nv ] ( const unsigned int itrack, + track_t const& tracks, + vertex_t const& vertices ) { + const double track_z = tracks._z[itrack]; + const double track_t = tracks._t[itrack]; + const double botrack_dz2 = -beta*tracks._dz2[itrack]; + const double botrack_dt2 = -beta*tracks._dt2[itrack]; + + // auto-vectorized + for ( unsigned int ivertex = 0; ivertex < nv; ++ivertex) { + const auto mult_resz = track_z - vertices._z[ivertex]; + const auto mult_rest = track_t - vertices._t[ivertex]; + vertices._ei_cache[ivertex] = botrack_dz2 * ( mult_resz * mult_resz ) + botrack_dt2 * ( mult_rest * mult_rest ); + } + }; + + auto kernel_add_Z = [ nv, Z_init ] (vertex_t const& vertices) -> double + { + double ZTemp = Z_init; + for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) { + ZTemp += vertices._pk[ivertex] * vertices._ei[ivertex]; + } + return ZTemp; + }; + + auto kernel_calc_normalization = [ beta, nv ] (const unsigned int track_num, + track_t & tks_vec, + vertex_t & y_vec ) { + auto tmp_trk_pi = tks_vec._pi[track_num]; + auto o_trk_Z_sum = 1./tks_vec._Z_sum[track_num]; + //auto o_trk_dz2 = tks_vec._dz2[track_num]; + //auto o_trk_dt2 = tks_vec._dt2[track_num]; + auto o_trk_err_sum = tks_vec._errsum[track_num]; + auto tmp_trk_z = tks_vec._z[track_num]; + auto tmp_trk_t = tks_vec._t[track_num]; + auto obeta = -1./beta; + + // auto-vectorized + for (unsigned int k = 0; k < nv; ++k) { + y_vec._se[k] += y_vec._ei[k] * (tmp_trk_pi* o_trk_Z_sum); + // parens are important for numerical stability + auto w = tmp_trk_pi * (y_vec._pk[k] * y_vec._ei[k] * o_trk_Z_sum) * o_trk_err_sum; + y_vec._sw[k] += w; + y_vec._swz[k] += w * tmp_trk_z; + y_vec._swt[k] += w * tmp_trk_t; + y_vec._swE[k] += w * y_vec._ei_cache[k]*obeta; + } + }; + + + for (auto ivertex = 0U; ivertex < nv; ++ivertex) { + gvertices._se[ivertex] = 0.0; + gvertices._sw[ivertex] = 0.0; + gvertices._swz[ivertex] = 0.0; + gvertices._swt[ivertex] = 0.0; + gvertices._swE[ivertex] = 0.0; + } + + + + // loop over tracks + for (auto itrack = 0U; itrack < nt; ++itrack) { + kernel_calc_exp_arg(itrack, gtracks, gvertices); + local_exp_list(gvertices._ei_cache, gvertices._ei, nv); + + gtracks._Z_sum[itrack] = kernel_add_Z(gvertices); + if (edm::isNotFinite(gtracks._Z_sum[itrack])) gtracks._Z_sum[itrack] = 0.0; + // used in the next major loop to follow + sumpi += gtracks._pi[itrack]; + + if (gtracks._Z_sum[itrack] > 1.e-100){ + kernel_calc_normalization(itrack, gtracks, gvertices); + } + } + + // now update z and pk + auto kernel_calc_zt = [ sumpi, nv, this, useRho0 ] (vertex_t & vertices ) -> double { + + double delta=0; + // does not vectorizes + for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex ) { + if (vertices._sw[ivertex] > 0) { + auto znew = vertices._swz[ ivertex ] / vertices._sw[ ivertex ]; + // prevents from vectorizing if + delta += std::pow( vertices._z[ ivertex ] - znew, 2 ); + vertices._z[ ivertex ] = znew; + auto tnew = vertices._swt[ ivertex ] / vertices._sw[ ivertex ]; + // prevents from vectorizing if + delta += std::pow( vertices._t[ ivertex ] - tnew, 2 ); + vertices._t[ ivertex ] = tnew; + } +#ifdef VI_DEBUG + else { + edm::LogInfo("sumw") << "invalid sum of weights in fit: " << vertices._sw[ivertex] << endl; + if (this->verbose_) { + std::cout << " a cluster melted away ? pk=" << vertices._pk[ ivertex ] << " sumw=" + << vertices._sw[ivertex] << endl; + } + } +#endif + } + + auto osumpi = 1./sumpi; + for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex ) + vertices._pk[ ivertex ] = vertices._pk[ ivertex ] * vertices._se[ ivertex ] * osumpi; + + return delta; + }; + + delta += kernel_calc_zt(gvertices); + + // return how much the prototypes moved + return delta; +} + + + + + +bool DAClusterizerInZT_vect::merge(vertex_t & y, double & beta)const{ + // merge clusters that collapsed or never separated, + // only merge if the estimated critical temperature of the merged vertex is below the current temperature + // return true if vertices were merged, false otherwise + const unsigned int nv = y.GetSize(); + + if (nv < 2) + return false; + + // merge the smallest distance clusters first + std::vector > critical; + for (unsigned int k = 0; (k + 1) < nv; k++) { + if ( std::fabs(y._z[k + 1] - y._z[k]) < zmerge_ && + std::fabs(y._t[k + 1] - y._t[k]) < tmerge_ ) { + auto dt2 = std::pow(y._t[k + 1] - y._t[k],2); + auto dz2 = std::pow(y._z[k + 1] - y._z[k],2); + critical.push_back( make_pair( dz2 + dt2, k ) ); + } + } + if (critical.size() == 0) return false; + + std::stable_sort(critical.begin(), critical.end(), std::less >() ); + + + for (unsigned int ik=0; ik < critical.size(); ik++){ + unsigned int k = critical[ik].second; + double rho = y._pk[k]+y._pk[k+1]; + double swE = y._swE[k]+y._swE[k+1]-y._pk[k]*y._pk[k+1] / rho*( std::pow(y._z[k+1]-y._z[k],2) + std::pow(y._t[k+1]-y._t[k],2) ); + double Tc = 2*swE / (y._sw[k]+y._sw[k+1]); + + if(Tc*beta < 1){ + if(verbose_){ std::cout << "merging (" << y._z[k + 1] << ',' << y._t[k + 1] << ") and (" << y._z[k] << ',' << y._t[k] << ") Tc = " << Tc << " sw = " << y._sw[k]+y._sw[k+1] < 0){ + y._z[k] = (y._pk[k]*y._z[k] + y._pk[k+1]*y._z[k + 1])/rho; + y._t[k] = (y._pk[k]*y._t[k] + y._pk[k+1]*y._t[k + 1])/rho; + }else{ + y._z[k] = 0.5 * (y._z[k] + y._z[k + 1]); + y._t[k] = 0.5 * (y._t[k] + y._t[k + 1]); + } + y._pk[k] = rho; + y._sw[k] += y._sw[k+1]; + y._swE[k] = swE; + y.RemoveItem(k+1); + return true; + } + } + + return false; +} + + + + +bool +DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const double beta) const { + // eliminate clusters with only one significant/unique track + const unsigned int nv = y.GetSize(); + const unsigned int nt = tks.GetSize(); + + if (nv < 2) + return false; + + double sumpmin = nt; + unsigned int k0 = nv; + + for (unsigned int k = 0; k < nv; k++) { + + int nUnique = 0; + double sump = 0; + + double pmax = y._pk[k] / (y._pk[k] + rho0 * local_exp(-beta * dzCutOff_* dzCutOff_)); + for (unsigned int i = 0; i < nt; i++) { + if (tks._Z_sum[i] > 1.e-100) { + double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i], tks._t[i], y._t[k], tks._dt2[i])) / tks._Z_sum[i]; + sump += p; + if ((p > uniquetrkweight_ * pmax) && (tks._pi[i] > 0)) { + nUnique++; + } + } + } + + if ((nUnique < 2) && (sump < sumpmin)) { + sumpmin = sump; + k0 = k; + } + + } + + if (k0 != nv) { + if (verbose_) { + std::cout << "eliminating prototype at " << std::setw(10) << std::setprecision(4) << y._z[k0] + << " with sump=" << sumpmin + << " rho*nt =" << y._pk[k0]*nt + << endl; + } + y.RemoveItem(k0); + return true; + } else { + return false; + } +} + + + + +double +DAClusterizerInZT_vect::beta0(double betamax, track_t const & tks, vertex_t const & y) const { + + double T0 = 0; // max Tc for beta=0 + // estimate critical temperature from beta=0 (T=inf) + const unsigned int nt = tks.GetSize(); + const unsigned int nv = y.GetSize(); + + for (unsigned int k = 0; k < nv; k++) { + + // vertex fit at T=inf + double sumwz = 0; + double sumwt = 0; + double sumw = 0; + for (unsigned int i = 0; i < nt; i++) { + double w = tks._pi[i] * tks._errsum[i]; + sumwz += w * tks._z[i]; + sumwt += w * tks._t[i]; + sumw += w; + } + y._z[k] = sumwz / sumw; + y._t[k] = sumwt / sumw; + + // estimate Tcrit, eventually do this in the same loop + double a = 0, b = 0; + for (unsigned int i = 0; i < nt; i++) { + double dx = tks._z[i] - y._z[k]; + double dt = tks._t[i] - y._t[k]; + double w = tks._pi[i] * tks._errsum[i]; + a += w * std::pow(dx, 2) * tks._dz2[i] + std::pow(dt,2) * tks._dt2[i]; + b += w; + } + double Tc = 2. * a / b; // the critical temperature of this vertex + if (Tc > T0) T0 = Tc; + }// vertex loop (normally there should be only one vertex at beta=0) + + if(verbose_){ + std::cout << "DAClustrizerInZT_vect.beta0: Tc = " << T0 << std::endl; + int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_)); + std::cout << "DAClustrizerInZT_vect.beta0: nstep = " << coolingsteps << std::endl; + } + + + if (T0 > 1. / betamax) { + return betamax / std::pow(coolingFactor_, int(std::log(T0 * betamax) / std::log(coolingFactor_)) - 1); + } else { + // ensure at least one annealing step + return betamax * coolingFactor_; + } +} + + + +bool +DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, double threshold ) const{ + // split only critical vertices (Tc >~ T=1/beta <==> beta*Tc>~1) + // an update must have been made just before doing this (same beta, no merging) + // returns true if at least one cluster was split + + constexpr double epsilonz=1e-3; // minimum split size z + constexpr double epsilont=1e-3; // minimum split size t + unsigned int nv = y.GetSize(); + + // avoid left-right biases by splitting highest Tc first + + std::vector > critical; + for(unsigned int k=0; k threshold){ + critical.push_back( make_pair(Tc, k)); + } + } + if (critical.size()==0) return false; + + + std::stable_sort(critical.begin(), critical.end(), std::greater >() ); + + + bool split=false; + const unsigned int nt = tks.GetSize(); + + for(unsigned int ic=0; ic 1.e-100) { + // winner-takes-all, usually overestimates splitting + double tl = tks._z[i] < y._z[k] ? 1.: 0.; + double tr = 1. - tl; + + // soften it, especially at low T + double arg = ( (tks._z[i] - y._z[k]) + (tks._t[i] - y._t[k]) ) * sqrt(beta * tks._errsum[i]); + if(std::fabs(arg) < 20){ + double t = local_exp(-arg); + tl = t/(t+1.); + tr = 1/(t+1.); + } + + double p = y._pk[k] * tks._pi[i] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i], + tks._t[i], y._t[k], tks._dt2[i])) / tks._Z_sum[i]; + double w = p*tks._errsum[i]; + p1 += p*tl ; z1 += w*tl*tks._z[i]; t1 += w*tl*tks._t[i]; w1 += w*tl; + p2 += p*tr; z2 += w*tr*tks._z[i]; t2 += w*tr*tks._t[i]; w2 += w*tr; + } + } + + if(w1>0){z1 = z1/w1; t1 = t1/w1;} else {z1=y._z[k]-epsilonz; t1=y._t[k]-epsilont;} + if(w2>0){z2 = z2/w2; t2 = t2/w2;} else {z2=y._z[k]+epsilonz; t2=y._t[k]+epsilont;} + + // reduce split size if there is not enough room + if( ( k > 0 ) && ( z1 < (0.6*y._z[k] + 0.4*y._z[k-1])) ){ z1 = 0.6*y._z[k] + 0.4*y._z[k-1]; } + if( ( k+1 < nv) && ( z2 > (0.6*y._z[k] + 0.4*y._z[k+1])) ){ z2 = 0.6*y._z[k] + 0.4*y._z[k+1]; } + + if(verbose_){ + if (std::fabs(y._z[k] - zdumpcenter_) < zdumpwidth_){ + std::cout << " T= " << std::setw(8) << 1./beta + << " Tc= " << critical[ic].first + << " splitting " << std::fixed << std::setprecision(4) << y._z[k] + << " --> (" << z1 << ',' << t1<< "),(" << z2 << ',' << t2 + << ") [" << p1 << "," << p2 << "]" ; + if (std::fabs(z2-z1) > epsilonz || std::fabs(t2-t1) > epsilont){ + std::cout << std::endl; + }else{ + std::cout << " rejected " << std::endl; + } + } + } + + // split if the new subclusters are significantly separated + if( std::fabs(z2-z1) > epsilonz || std::fabs(t2-t1) > epsilont){ + split = true; + double pk1 = p1*y._pk[k]/(p1+p2); + double pk2 = p2*y._pk[k]/(p1+p2); + y._z[k] = z2; + y._t[k] = t2; + y._pk[k] = pk2; + y.InsertItem(k, z1, t1, pk1); + nv++; + + // adjust remaining pointers + for(unsigned int jc=ic; jc < critical.size(); jc++){ + if (critical[jc].second > k) {critical[jc].second++;} + } + } + } + return split; +} + + + +void DAClusterizerInZT_vect::splitAll( vertex_t & y) const { + + const unsigned int nv = y.GetSize(); + + constexpr double epsilonz = 1e-3; // split all single vertices by 10 um + constexpr double epsilont = 1e-2; // split all single vertices by 10 ps + constexpr double zsep = 2 * epsilonz; // split vertices that are isolated by at least zsep (vertices that haven't collapsed) + constexpr double tsep = 2 * epsilont; + vertex_t y1; + + if (verbose_) { + std::cout << "Before Split "<< std::endl; + y.DebugOut(); + } + + for (unsigned int k = 0; k < nv; k++) { + if ( + ( ( (k == 0) || ( y._z[k - 1] < (y._z[k] - zsep)) ) && + ( ((k + 1) == nv) || ( y._z[k + 1] > (y._z[k] + zsep)) ) ) || + ( ( (k == 0) || ( y._t[k - 1] < (y._t[k] - tsep)) ) && + ( ((k + 1) == nv) || ( y._t[k + 1] > (y._t[k] + tsep)) ) ) ) + { + // isolated prototype, split + double new_z = y.z[k] - epsilonz; + double new_t = y.t[k] - epsilont; + y.z[k] = y.z[k] + epsilonz; + y.t[k] = y.t[k] + epsilont; + + double new_pk = 0.5 * y.pk[k]; + y.pk[k] = 0.5 * y.pk[k]; + + y1.AddItem(new_z, new_t, new_pk); + y1.AddItem(y._z[k], y._t[k], y._pk[k]); + } + else if ( (y1.GetSize() == 0 ) || + (y1._z[y1.GetSize() - 1] < (y._z[k] - zsep) ) || + (y1._t[y1.GetSize() - 1] < (y._t[k] - tsep) )) + { + y1.AddItem(y._z[k], y._t[k], y._pk[k]); + } + else + { + y1._z[y1.GetSize() - 1] = y1._z[y1.GetSize() - 1] - epsilonz; + y1._t[y1.GetSize() - 1] = y1._t[y1.GetSize() - 1] - epsilont; + y._z[k] = y._z[k] + epsilonz; + y._t[k] = y._t[k] + epsilont; + y1.AddItem( y._z[k], y._t[k] , y._pk[k]); + } + }// vertex loop + + y = y1; + y.ExtractRaw(); + + if (verbose_) { + std::cout << "After split " << std::endl; + y.DebugOut(); + } +} + + + +vector +DAClusterizerInZT_vect::vertices(const vector & tracks, const int verbosity) const { + track_t && tks = fill(tracks); + tks.ExtractRaw(); + + unsigned int nt = tks.GetSize(); + double rho0 = 0.0; // start with no outlier rejection + + vector clusters; + if (tks.GetSize() == 0) return clusters; + + vertex_t y; // the vertex prototypes + + // initialize:single vertex at infinite temperature + y.AddItem( 0, 0, 1.0); + + int niter = 0; // number of iterations + + + // estimate first critical temperature + double beta = beta0(betamax_, tks, y); + if ( verbose_) std::cout << "Beta0 is " << beta << std::endl; + + niter = 0; + while ((update(beta, tks, y, false, rho0) > 1.e-6) && + (niter++ < maxIterations_)) {} + + // annealing loop, stop when T1/Tmin) + + double betafreeze = betamax_ * sqrt(coolingFactor_); + + while (beta < betafreeze) { + if(useTc_){ + update(beta, tks,y, false, rho0); + while(merge(y, beta)){update(beta, tks, y, false, rho0);} + split(beta, tks, y); + beta=beta/coolingFactor_; + }else{ + beta=beta/coolingFactor_; + splitAll(y); + } + + // make sure we are not too far from equilibrium before cooling further + niter = 0; + while ((update(beta, tks, y, false, rho0) > 1.e-6) && + (niter++ < maxIterations_)) {} + + if(verbose_){ dump( beta, y, tks, 0); } + } + + + if(useTc_){ + //last round of splitting, make sure no critical clusters are left + if(verbose_){ std::cout << "last spliting at " << 1./beta << std::endl; } + update(beta, tks,y, false, rho0);// make sure Tc is up-to-date + while(merge(y,beta)){update(beta, tks,y, false, rho0);} + unsigned int ntry=0; + double threshold = 1.0; + while( split(beta, tks, y, threshold) && (ntry++<10) ){ + niter=0; + while((update(beta, tks,y, false, rho0)>1.e-6) && (niter++ < maxIterations_)){} + while(merge(y,beta)){update(beta, tks,y, false, rho0);} + if(verbose_){ + std::cout << "after final splitting, try " << ntry << std::endl; + dump(beta, y, tks, 2); + } + // relax splitting a bit to reduce multiple split-merge cycles of the same cluster + threshold *= 1.1; + } + }else{ + // merge collapsed clusters + while(merge(y,beta)){update(beta, tks,y, false, rho0);} + } + + if (verbose_) { + update(beta, tks,y, false, rho0); + std::cout << "dump after 1st merging " << endl; + dump(beta, y, tks, 2); + } + + + // switch on outlier rejection at T=Tmin + if(dzCutOff_ > 0){ + rho0 = 1./nt; + for(unsigned int a=0; a<10; a++){ update(beta, tks, y, true, a*rho0/10);} // adiabatic turn-on + } + + niter=0; + while ((update(beta, tks, y, true, rho0) > 1.e-8) && (niter++ < maxIterations_)) {}; + if (verbose_) { + std::cout << "dump after noise-suppression, rho0=" << rho0 << " niter = " << niter << endl; + dump(beta, y, tks, 2); + } + + // merge again (some cluster split by outliers collapse here) + while (merge(y, beta)) {update(beta, tks, y, true, rho0); } + if (verbose_) { + std::cout << "dump after merging " << endl; + dump(beta, y, tks, 2); + } + + // go down to the purging temperature (if it is lower than tmin) + while( beta < betapurge_ ){ + beta = min( beta/coolingFactor_, betapurge_); + niter = 0; + while ((update(beta, tks, y, false, rho0) > 1.e-8) && (niter++ < maxIterations_)) {} + } + + + // eliminate insigificant vertices, this is more restrictive at higher T + while (purge(y, tks, rho0, beta)) { + niter = 0; + while ((update(beta, tks, y, true, rho0) > 1.e-6) && (niter++ < maxIterations_)) {} + } + + if (verbose_) { + update(beta, tks,y, true, rho0); + std::cout << " after purging " << std:: endl; + dump(beta, y, tks, 2); + } + + // optionally cool some more without doing anything, to make the assignment harder + while( beta < betastop_ ){ + beta = min( beta/coolingFactor_, betastop_); + niter =0; + while ((update(beta, tks, y, true, rho0) > 1.e-8) && (niter++ < maxIterations_)) {} + } + + + if (verbose_) { + std::cout << "Final result, rho0=" << std::scientific << rho0 << endl; + dump(beta, y, tks, 2); + } + + // select significant tracks and use a TransientVertex as a container + GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01); + + // ensure correct normalization of probabilities, should makes double assignment 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; y._z[k]=0;} + + for (unsigned int i = 0; i < nt; i++) // initialize + tks._Z_sum[i] = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); + + // improve vectorization (does not require reduction ....) + for (unsigned int k = 0; k < nv; k++) { + for (unsigned int i = 0; i < nt; i++) + tks._Z_sum[i] += y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],tks._dz2[i], tks._t[i], y._t[k],tks._dt2[i])); + } + + + for (unsigned int k = 0; k < nv; k++) { + GlobalPoint pos(0, 0, y._z[k]); + + vector vertexTracks; + for (unsigned int i = 0; i < nt; i++) { + if (tks._Z_sum[i] > 1e-100) { + + double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i], + tks._t[i], y._t[k], tks._dt2[i] )) / tks._Z_sum[i]; + if ((tks._pi[i] > 0) && (p > mintrkweight_)) { + vertexTracks.push_back(*(tks.tt[i])); + tks._Z_sum[i] = 0; // setting Z=0 excludes double assignment + } + } + } + TransientVertex v(pos, dummyError, vertexTracks, 0); + clusters.push_back(v); + } + + return clusters; + +} + +vector > DAClusterizerInZT_vect::clusterize( + const vector & tracks) const { + + if (verbose_) { + std::cout << "###################################################" << endl; + std::cout << "# vectorized DAClusterizerInZT_vect::clusterize nt=" << tracks.size() << endl; + std::cout << "###################################################" << endl; + } + + vector > clusters; + vector && pv = vertices(tracks); + + if (verbose_) { + std::cout << "# DAClusterizerInZT::clusterize pv.size=" << pv.size() + << endl; + } + if (pv.size() == 0) { + return clusters; + } + + // fill into clusters and merge + vector aCluster = pv.begin()->originalTracks(); + + for (auto k = pv.begin() + 1; k != pv.end(); k++) { + if ( std::abs(k->position().z() - (k - 1)->position().z()) > (2 * vertexSize_)) { + // close a cluster + if (aCluster.size()>1){ + clusters.push_back(aCluster); + }else{ + if(verbose_){ + std::cout << " one track cluster at " << k->position().z() << " suppressed" << std::endl; + } + } + aCluster.clear(); + } + for (unsigned int i = 0; i < k->originalTracks().size(); i++) { + aCluster.push_back(k->originalTracks()[i]); + } + + } + clusters.emplace_back(std::move(aCluster)); + + return clusters; + +} + + + +void DAClusterizerInZT_vect::dump(const double beta, const vertex_t & y, + const track_t & tks, int verbosity) const { + + const unsigned int nv = y.GetSize(); + const unsigned int nt = tks.GetSize(); + + std::vector< unsigned int > iz; + for(unsigned int j=0; j 0) { + F -= std::log(tks._Z_sum[i]) / beta; + } + double tz = tks._z[i]; + + if( std::fabs(tz - zdumpcenter_) > zdumpwidth_) continue; + std::cout << setw(4) << i << ")" << setw(8) << fixed << setprecision(4) + << tz << " +/-" << setw(6) << sqrt(1./tks._dz2[i]); + + if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) { + std::cout << " *"; + } else { + std::cout << " "; + } + if (tks.tt[i]->track().hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) { + std::cout << "+"; + } else { + std::cout << "-"; + } + std::cout << setw(1) + << tks.tt[i]->track().hitPattern().pixelBarrelLayersWithMeasurement(); // see DataFormats/TrackReco/interface/HitPattern.h + std::cout << setw(1) + << tks.tt[i]->track().hitPattern().pixelEndcapLayersWithMeasurement(); + std::cout << setw(1) << hex + << tks.tt[i]->track().hitPattern().trackerLayersWithMeasurement() + - tks.tt[i]->track().hitPattern().pixelLayersWithMeasurement() + << dec; + std::cout << "=" << setw(1) << hex + << tks.tt[i]->track().hitPattern().numberOfHits(reco::HitPattern::MISSING_OUTER_HITS) + << dec; + + Measurement1D IP = + tks.tt[i]->stateAtBeamLine().transverseImpactParameter(); + std::cout << setw(8) << IP.value() << "+/-" << setw(6) << IP.error(); + std::cout << " " << setw(6) << setprecision(2) + << tks.tt[i]->track().pt() * tks.tt[i]->track().charge(); + std::cout << " " << setw(5) << setprecision(2) + << tks.tt[i]->track().phi() << " " << setw(5) + << setprecision(2) << tks.tt[i]->track().eta(); + + double sump = 0.; + for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex) { + if (std::fabs(y._z[ivertex]-zdumpcenter_) > zdumpwidth_) continue; + + if ((tks._pi[i] > 0) && (tks._Z_sum[i] > 0)) { + //double p=pik(beta,tks[i],*k); + double p = y._pk[ivertex] * exp(-beta * Eik(tks._z[i], y._z[ivertex], tks._dz2[i], + tks._t[i], y._t[ivertex], tks._dt2[i] )) / tks._Z_sum[i]; + if (p > 0.0001) { + std::cout << setw(8) << setprecision(3) << p; + } else { + std::cout << " . "; + } + E += p * Eik(tks._z[i], y._z[ivertex], tks._dz2[i], + tks._t[i], y._t[ivertex], tks._dt2[i] ); + sump += p; + } else { + std::cout << " "; + } + } + std::cout << endl; + } + std::cout << endl << "T=" << 1 / beta << " E=" << E << " n=" << y.GetSize() + << " F= " << F << endl << "----------" << endl; + } +} From 9154976ba5aa664511ce377df8b89c324b9e300e Mon Sep 17 00:00:00 2001 From: lgray Date: Tue, 1 Aug 2017 16:23:00 -0500 Subject: [PATCH 02/16] make purging loop use an adaptive tolerance --- .../python/TkClusParameters_cff.py | 10 ++-- .../src/DAClusterizerInZT_vect.cc | 58 +++++++++++-------- 2 files changed, 38 insertions(+), 30 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py index 8965c98dac5c4..32f0b455fd762 100644 --- a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py +++ b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py @@ -5,7 +5,7 @@ TkDAClusParameters = cms.PSet( coolingFactor = cms.double(0.6), # moderate annealing speed Tmin = cms.double(4.), # end of annealing - vertexSize = cms.double(0.01), # ~ resolution / sqrt(Tmin) + vertexSize = cms.double(0.006), # ~ resolution / sqrt(Tmin) d0CutOff = cms.double(3.), # downweight high IP tracks dzCutOff = cms.double(4.) # outlier rejection after freeze-out (T 0) { + //std::cout << " DA2D_vect sw = " << vertices._sw[ ivertex ] << ' ' << vertices._swz[ ivertex ] << ' ' << vertices._swt[ ivertex ] << std::endl; auto znew = vertices._swz[ ivertex ] / vertices._sw[ ivertex ]; // prevents from vectorizing if delta += std::pow( vertices._z[ ivertex ] - znew, 2 ); @@ -278,7 +278,7 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, } auto osumpi = 1./sumpi; - for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex ) + for (unsigned int ivertex = 0; ivertex < nv; ++ivertex ) vertices._pk[ ivertex ] = vertices._pk[ ivertex ] * vertices._se[ ivertex ] * osumpi; return delta; @@ -359,18 +359,18 @@ DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const double sumpmin = nt; unsigned int k0 = nv; - for (unsigned int k = 0; k < nv; k++) { + for (unsigned int k = 0; k < nv; ++k) { int nUnique = 0; double sump = 0; double pmax = y._pk[k] / (y._pk[k] + rho0 * local_exp(-beta * dzCutOff_* dzCutOff_)); - for (unsigned int i = 0; i < nt; i++) { + for (unsigned int i = 0; i < nt; ++i) { if (tks._Z_sum[i] > 1.e-100) { double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i], tks._t[i], y._t[k], tks._dt2[i])) / tks._Z_sum[i]; sump += p; if ((p > uniquetrkweight_ * pmax) && (tks._pi[i] > 0)) { - nUnique++; + ++nUnique; } } } @@ -459,13 +459,13 @@ DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, do // returns true if at least one cluster was split constexpr double epsilonz=1e-3; // minimum split size z - constexpr double epsilont=1e-3; // minimum split size t + constexpr double epsilont=1e-2; // minimum split size t unsigned int nv = y.GetSize(); // avoid left-right biases by splitting highest Tc first std::vector > critical; - for(unsigned int k=0; k threshold){ critical.push_back( make_pair(Tc, k)); @@ -492,29 +492,38 @@ DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, do double tl = tks._z[i] < y._z[k] ? 1.: 0.; double tr = 1. - tl; - // soften it, especially at low T - double arg = ( (tks._z[i] - y._z[k]) + (tks._t[i] - y._t[k]) ) * sqrt(beta * tks._errsum[i]); + // soften it, especially at low T + /* + double arg = ( tks._z[i] - y._z[k] ) * sqrt(beta * tks._dz2[i]); // + std::fabs(tks._t[i] - y._t[k]) + std::cout << arg << std::endl; if(std::fabs(arg) < 20){ double t = local_exp(-arg); tl = t/(t+1.); tr = 1/(t+1.); } + */ double p = y._pk[k] * tks._pi[i] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i], tks._t[i], y._t[k], tks._dt2[i])) / tks._Z_sum[i]; double w = p*tks._errsum[i]; - p1 += p*tl ; z1 += w*tl*tks._z[i]; t1 += w*tl*tks._t[i]; w1 += w*tl; + p1 += p*tl; z1 += w*tl*tks._z[i]; t1 += w*tl*tks._t[i]; w1 += w*tl; p2 += p*tr; z2 += w*tr*tks._z[i]; t2 += w*tr*tks._t[i]; w2 += w*tr; } } if(w1>0){z1 = z1/w1; t1 = t1/w1;} else {z1=y._z[k]-epsilonz; t1=y._t[k]-epsilont;} if(w2>0){z2 = z2/w2; t2 = t2/w2;} else {z2=y._z[k]+epsilonz; t2=y._t[k]+epsilont;} - + // reduce split size if there is not enough room - if( ( k > 0 ) && ( z1 < (0.6*y._z[k] + 0.4*y._z[k-1])) ){ z1 = 0.6*y._z[k] + 0.4*y._z[k-1]; } - if( ( k+1 < nv) && ( z2 > (0.6*y._z[k] + 0.4*y._z[k+1])) ){ z2 = 0.6*y._z[k] + 0.4*y._z[k+1]; } - + if( ( k > 0 ) && ( z1 < (0.6*y._z[k] + 0.4*y._z[k-1])) ){ + z1 = 0.5*y._z[k] + 0.5*y._z[k-1]; + t1 = 0.5*y._t[k] + 0.5*y._t[k-1]; + } + if( ( k+1 < nv) && ( z2 > (0.6*y._z[k] + 0.4*y._z[k+1])) ){ + z2 = 0.5*y._z[k] + 0.5*y._z[k+1]; + t2 = 0.5*y._t[k] + 0.5*y._t[k+1]; + } + if(verbose_){ if (std::fabs(y._z[k] - zdumpcenter_) < zdumpwidth_){ std::cout << " T= " << std::setw(8) << 1./beta @@ -570,9 +579,7 @@ void DAClusterizerInZT_vect::splitAll( vertex_t & y) const { for (unsigned int k = 0; k < nv; k++) { if ( ( ( (k == 0) || ( y._z[k - 1] < (y._z[k] - zsep)) ) && - ( ((k + 1) == nv) || ( y._z[k + 1] > (y._z[k] + zsep)) ) ) || - ( ( (k == 0) || ( y._t[k - 1] < (y._t[k] - tsep)) ) && - ( ((k + 1) == nv) || ( y._t[k + 1] > (y._t[k] + tsep)) ) ) ) + ( ((k + 1) == nv) || ( y._z[k + 1] > (y._z[k] + zsep)) ) ) ) { // isolated prototype, split double new_z = y.z[k] - epsilonz; @@ -725,7 +732,8 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co // eliminate insigificant vertices, this is more restrictive at higher T while (purge(y, tks, rho0, beta)) { niter = 0; - while ((update(beta, tks, y, true, rho0) > 1.e-6) && (niter++ < maxIterations_)) {} + while (( update(beta, tks, y, true, rho0) > 2.5e-7 * y.GetSize() ) && (niter++ < maxIterations_)) { + } } if (verbose_) { From 1663cc162d55aa78c884126296f17b18694304ca Mon Sep 17 00:00:00 2001 From: lgray Date: Wed, 2 Aug 2017 18:00:59 -0500 Subject: [PATCH 03/16] use aligned vectors everywhere --- .../PrimaryVertexProducer/BuildFile.xml | 2 +- .../interface/DAClusterizerInZT_vect.h | 147 +++++++++++++----- .../src/DAClusterizerInZT_vect.cc | 21 +-- 3 files changed, 124 insertions(+), 46 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/BuildFile.xml b/RecoVertex/PrimaryVertexProducer/BuildFile.xml index 7f5f2666c0560..0b584c37b5cd1 100644 --- a/RecoVertex/PrimaryVertexProducer/BuildFile.xml +++ b/RecoVertex/PrimaryVertexProducer/BuildFile.xml @@ -1,4 +1,4 @@ - + diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h index 196fd88321d97..b6991cf08d405 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -17,10 +17,85 @@ #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h" #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include +#include +#include + +template +class AlignmentAllocator { +public: + typedef T value_type; + typedef std::size_t size_type; + typedef std::ptrdiff_t difference_type; + + typedef T * pointer; + typedef const T * const_pointer; + + typedef T & reference; + typedef const T & const_reference; + + public: + inline AlignmentAllocator () throw () { } + + template + inline AlignmentAllocator (const AlignmentAllocator &) throw () { } + + inline ~AlignmentAllocator () throw () { } + + inline pointer adress (reference r) { + return &r; + } + + inline const_pointer adress (const_reference r) const { + return &r; + } + + inline pointer allocate (size_type n) { + void * p = nullptr; + posix_memalign(&p,N,n*sizeof(value_type)); + return (pointer)p; + } + + inline void deallocate (pointer p, size_type) { + free(p); + } + + inline void construct (pointer p, const value_type & wert) { + new (p) value_type (wert); + } + + inline void destroy (pointer p) { + p->~value_type (); + } + + inline size_type max_size () const throw () { + return size_type (-1) / sizeof (value_type); + } + + template + struct rebind { + typedef AlignmentAllocator other; + }; + + bool operator!=(const AlignmentAllocator& other) const { + return !(*this == other); + } + + // Returns true if and only if storage allocated from *this + // can be deallocated from other, and vice versa. + // Always returns true for stateless allocators. + bool operator==(const AlignmentAllocator& other) const { + return true; + } +}; + class DAClusterizerInZT_vect final : public TrackClusterizerInZ { public: + template + using AlignedVector = std::vector >; + // Internal data structure to struct track_t { @@ -57,39 +132,39 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { _pi = &pi.front(); } - double * __restrict__ _z; // z-coordinate at point of closest approach to the beamline - double * __restrict__ _t; // t-coordinate at point of closest approach to the beamline - double * __restrict__ _dz2; // square of the error of z(pca) - double * __restrict__ _dt2; // square of the error of t(pca) - double * __restrict__ _errsum; // sum of squares of the pca errors + double * __restrict__ _z __attribute__ ((aligned (16))); // z-coordinate at point of closest approach to the beamline + double * __restrict__ _t __attribute__ ((aligned (16))); // t-coordinate at point of closest approach to the beamline + double * __restrict__ _dz2 __attribute__ ((aligned (16))); // square of the error of z(pca) + double * __restrict__ _dt2 __attribute__ ((aligned (16))); // square of the error of t(pca) + double * __restrict__ _errsum __attribute__ ((aligned (16))); // sum of squares of the pca errors - double * __restrict__ _Z_sum; // Z[i] for DA clustering - double * __restrict__ _pi; // track weight + double * __restrict__ _Z_sum __attribute__ ((aligned (16))); // Z[i] for DA clustering + double * __restrict__ _pi __attribute__ ((aligned (16))); // track weight - std::vector z; // z-coordinate at point of closest approach to the beamline - std::vector t; // t-coordinate at point of closest approach to the beamline - std::vector dz2; // square of the error of z(pca) - std::vector dt2; // square of the error of t(pca) - std::vector errsum; // sum of squares of the pca errors + AlignedVector z; // z-coordinate at point of closest approach to the beamline + AlignedVector t; // t-coordinate at point of closest approach to the beamline + AlignedVector dz2; // square of the error of z(pca) + AlignedVector dt2; // square of the error of t(pca) + AlignedVector errsum; // sum of squares of the pca errors + AlignedVector Z_sum; // Z[i] for DA clustering + AlignedVector pi; // track weight + std::vector< const reco::TransientTrack* > tt; // a pointer to the Transient Track - - std::vector Z_sum; // Z[i] for DA clustering - std::vector pi; // track weight }; struct vertex_t { - std::vector z; // z coordinate - std::vector t; // t coordinate - std::vector pk; // vertex weight for "constrained" clustering + AlignedVector z; // z coordinate + AlignedVector t; // t coordinate + AlignedVector pk; // vertex weight for "constrained" clustering // --- temporary numbers, used during update - std::vector ei_cache; - std::vector ei; - std::vector sw; - std::vector swz; - std::vector swt; - std::vector se; - std::vector swE; + AlignedVector ei_cache; + AlignedVector ei; + AlignedVector sw; + AlignedVector swz; + AlignedVector swt; + AlignedVector se; + AlignedVector swE; unsigned int GetSize() const @@ -175,17 +250,17 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { } - double * __restrict__ _z; - double * __restrict__ _t; - double * __restrict__ _pk; - - double * __restrict__ _ei_cache; - double * __restrict__ _ei; - double * __restrict__ _sw; - double * __restrict__ _swz; - double * __restrict__ _swt; - double * __restrict__ _se; - double * __restrict__ _swE; + double * __restrict__ _z __attribute__ ((aligned (16))); + double * __restrict__ _t __attribute__ ((aligned (16))); + double * __restrict__ _pk __attribute__ ((aligned (16))); + + double * __restrict__ _ei_cache __attribute__ ((aligned (16))); + double * __restrict__ _ei __attribute__ ((aligned (16))); + double * __restrict__ _sw __attribute__ ((aligned (16))); + double * __restrict__ _swz __attribute__ ((aligned (16))); + double * __restrict__ _swt __attribute__ ((aligned (16))); + double * __restrict__ _se __attribute__ ((aligned (16))); + double * __restrict__ _swE __attribute__ ((aligned (16))); }; diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index 8265e9137c8a2..06a90d1b3c2db 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -94,8 +94,10 @@ namespace { return vdt::fast_exp( inp ); } - inline void local_exp_list( double const * __restrict__ arg_inp, double * __restrict__ arg_out, const int arg_arr_size) { - for(int i=0; i!=arg_arr_size; ++i ) arg_out[i]=vdt::fast_exp(arg_inp[i]); + inline void local_exp_list( double const * __restrict__ arg_inp, + double * __restrict__ arg_out, + const unsigned arg_arr_size) { + for(unsigned i=0; i!=arg_arr_size; ++i ) arg_out[i]=vdt::fast_exp(arg_inp[i]); } } @@ -160,10 +162,10 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, const unsigned int nv = gvertices.GetSize(); //initialize sums - double sumpi = 0; + double sumpi = 0.; // to return how much the prototype moved - double delta = 0; + double delta = 0.; // intial value of a sum @@ -178,6 +180,7 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, auto kernel_calc_exp_arg = [ beta, nv ] ( const unsigned int itrack, track_t const& tracks, vertex_t const& vertices ) { + const auto track_z = tracks._z[itrack]; const auto track_t = tracks._t[itrack]; const auto botrack_dz2 = -beta*tracks._dz2[itrack]; @@ -185,8 +188,8 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, // auto-vectorized for ( unsigned int ivertex = 0; ivertex < nv; ++ivertex) { - const auto mult_resz = track_z - vertices._z[ivertex]; - const auto mult_rest = track_t - vertices._t[ivertex]; + const auto mult_resz = track_z - vertices._z[ivertex]; + const auto mult_rest = track_t - vertices._t[ivertex]; vertices._ei_cache[ivertex] = botrack_dz2 * ( mult_resz * mult_resz ) + botrack_dt2 * ( mult_rest * mult_rest ); } }; @@ -238,13 +241,13 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, for (auto itrack = 0U; itrack < nt; ++itrack) { kernel_calc_exp_arg(itrack, gtracks, gvertices); local_exp_list(gvertices._ei_cache, gvertices._ei, nv); - + gtracks._Z_sum[itrack] = kernel_add_Z(gvertices); if (edm::isNotFinite(gtracks._Z_sum[itrack])) gtracks._Z_sum[itrack] = 0.0; // used in the next major loop to follow sumpi += gtracks._pi[itrack]; - if (gtracks._Z_sum[itrack] > 1.e-100){ + if (gtracks._Z_sum[itrack] > 1.d-100){ kernel_calc_normalization(itrack, gtracks, gvertices); } } @@ -255,7 +258,7 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, double delta=0; // does not vectorizes for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex ) { - if (vertices._sw[ivertex] > 0) { + if (vertices._sw[ivertex] > 0.) { //std::cout << " DA2D_vect sw = " << vertices._sw[ ivertex ] << ' ' << vertices._swz[ ivertex ] << ' ' << vertices._swt[ ivertex ] << std::endl; auto znew = vertices._swz[ ivertex ] / vertices._sw[ ivertex ]; // prevents from vectorizing if From 4c14cb1d631a9c537d25993378b7583f6a2eb7ac Mon Sep 17 00:00:00 2001 From: lgray Date: Thu, 3 Aug 2017 17:56:50 -0500 Subject: [PATCH 04/16] further vectorization, few more seconds down in 200PU --- .../src/DAClusterizerInZT_vect.cc | 53 +++++++++++++------ 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index 06a90d1b3c2db..7302bfdd0973c 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -247,7 +247,7 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, // used in the next major loop to follow sumpi += gtracks._pi[itrack]; - if (gtracks._Z_sum[itrack] > 1.d-100){ + if (gtracks._Z_sum[itrack] > 1.e-100){ kernel_calc_normalization(itrack, gtracks, gvertices); } } @@ -268,13 +268,13 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, // prevents from vectorizing if delta += std::pow( vertices._t[ ivertex ] - tnew, 2 ); vertices._t[ ivertex ] = tnew; - } + } #ifdef VI_DEBUG else { edm::LogInfo("sumw") << "invalid sum of weights in fit: " << vertices._sw[ivertex] << endl; if (this->verbose_) { std::cout << " a cluster melted away ? pk=" << vertices._pk[ ivertex ] << " sumw=" - << vertices._sw[ivertex] << endl; + << vertices._sw[ivertex] << endl; } } #endif @@ -352,30 +352,53 @@ bool DAClusterizerInZT_vect::merge(vertex_t & y, double & beta)const{ bool DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const double beta) const { + constexpr double eps = 1.e-100; // eliminate clusters with only one significant/unique track const unsigned int nv = y.GetSize(); const unsigned int nt = tks.GetSize(); - + if (nv < 2) return false; double sumpmin = nt; unsigned int k0 = nv; + int nUnique = 0; + double sump = 0; + + AlignedVector inverse_zsums(nt), arg_cache(nt), eik_cache(nt); + double * __restrict__ _inverse_zsums __attribute__ ((aligned (16))); + double * __restrict__ _arg_cache __attribute__ ((aligned (16))); + double * __restrict__ _eik_cache __attribute__ ((aligned (16))); + _inverse_zsums = inverse_zsums.data(); + _arg_cache = arg_cache.data(); + _eik_cache = eik_cache.data(); + for(unsigned i = 0; i < nt; ++i) { + inverse_zsums[i] = tks._Z_sum[i] > eps ? 1./tks._Z_sum[i] : 0.0; + } + for (unsigned int k = 0; k < nv; ++k) { - int nUnique = 0; - double sump = 0; - - double pmax = y._pk[k] / (y._pk[k] + rho0 * local_exp(-beta * dzCutOff_* dzCutOff_)); + nUnique = 0; + sump = 0; + + const double pmax = y._pk[k] / (y._pk[k] + rho0 * local_exp(-beta * dzCutOff_* dzCutOff_)); + const double pcut = uniquetrkweight_ * pmax; + for(unsigned i = 0; i < nt; ++i) { + const auto track_z = tks._z[i]; + const auto track_t = tks._t[i]; + const auto botrack_dz2 = -beta*tks._dz2[i]; + const auto botrack_dt2 = -beta*tks._dt2[i]; + + const auto mult_resz = track_z - y._z[k]; + const auto mult_rest = track_t - y._t[k]; + _arg_cache[i] = botrack_dz2 * ( mult_resz * mult_resz ) + botrack_dt2 * ( mult_rest * mult_rest ); + } + local_exp_list(_arg_cache, _eik_cache, nt); for (unsigned int i = 0; i < nt; ++i) { - if (tks._Z_sum[i] > 1.e-100) { - double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i], tks._t[i], y._t[k], tks._dt2[i])) / tks._Z_sum[i]; - sump += p; - if ((p > uniquetrkweight_ * pmax) && (tks._pi[i] > 0)) { - ++nUnique; - } - } + const double p = y._pk[k] * _eik_cache[i] * _inverse_zsums[i]; + sump += p; + nUnique += ( ( p > pcut ) & ( tks._pi[i] > 0 ) ); } if ((nUnique < 2) && (sump < sumpmin)) { From 6793ce53c5be9affb76e07ce84e626c79db146f9 Mon Sep 17 00:00:00 2001 From: lgray Date: Tue, 15 Aug 2017 12:47:03 -0500 Subject: [PATCH 05/16] remove nonvectorized code sync vectorized and non vectorized settings --- .../interface/DAClusterizerInZT.h | 102 --- .../interface/DAClusterizerInZT_vect.h | 145 +--- .../python/TkClusParameters_cff.py | 13 +- .../src/DAClusterizerInZT.cc | 703 ------------------ .../src/DAClusterizerInZT_vect.cc | 55 +- 5 files changed, 84 insertions(+), 934 deletions(-) delete mode 100644 RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT.h delete mode 100644 RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT.cc diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT.h deleted file mode 100644 index 1fdb3ceeed359..0000000000000 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT.h +++ /dev/null @@ -1,102 +0,0 @@ -#ifndef RecoVertex_PrimaryVertexProducer_DAClusterizerInZT_h -#define RecoVertex_PrimaryVertexProducer_DAClusterizerInZT_h - -/**\class DAClusterizerInZT - - Description: separates event tracks into clusters along the beam line - -*/ - -#include "RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h" -#include "TrackingTools/TransientTrack/interface/TransientTrack.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include -#include "DataFormats/Math/interface/Error.h" -#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h" -#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" - -class DAClusterizerInZT : public TrackClusterizerInZ { - -public: - -struct track_t{ - double z; // z-coordinate at point of closest approach to the beamline - double t; // t-coordinate at point of closest approach to the beamline - double dz2; // square of the error of z(pca) - double dtz; // covariance of z-t - double dt2; // square of the error of t(pca) - const reco::TransientTrack* tt; // a pointer to the Transient Track - double zi; // Z[i] for DA clustering - double pi; // track weight -}; - - -struct vertex_t{ - double z; // z coordinate - double t; // t coordinate - double pk; // vertex weight for "constrained" clustering - // --- temporary numbers, used during update - double ei; - double sw; - double swz; - double swt; - double se; - // ---for Tc - double swE; - double tC; -}; - - - - - DAClusterizerInZT(const edm::ParameterSet& conf); - - std::vector< std::vector > - clusterize(const std::vector & tracks)const; - - - std::vector< TransientVertex > - vertices(const std::vector & tracks, const int verbosity=0)const; - - - std::vector fill(const std::vector & tracks)const; - - bool split( double beta, - std::vector & tks, - std::vector & y, - double threshold ) const; - - double update(double beta, - std::vector & tks, - std::vector & y, - const double rho0 = 0.0 )const; - - void dump(const double beta, const std::vector & y, const std::vector & tks, const int verbosity=0) const; - bool merge(std::vector &,int ) const; - bool merge(std::vector &,double & ) const; - bool purge(std::vector &, std::vector & , double &, const double ) const; - - void splitAll( std::vector & y ) const; - - double beta0(const double betamax, - std::vector & tks, - std::vector & y )const; - - double e_ik(const track_t & t, const vertex_t & k)const; - - -private: - bool verbose_; - bool useTc_; - float vertexSize_; - int maxIterations_; - double coolingFactor_; - double logCoolingFactor_; - float betamax_; - float betastop_; - double dzCutOff_; - double d0CutOff_; - double dtCutOff_; // for when the beamspot has time -}; - -#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h index b6991cf08d405..3acc66d96e439 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -18,84 +18,11 @@ #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" #include -#include -#include - -template -class AlignmentAllocator { -public: - typedef T value_type; - typedef std::size_t size_type; - typedef std::ptrdiff_t difference_type; - - typedef T * pointer; - typedef const T * const_pointer; - - typedef T & reference; - typedef const T & const_reference; - - public: - inline AlignmentAllocator () throw () { } - - template - inline AlignmentAllocator (const AlignmentAllocator &) throw () { } - - inline ~AlignmentAllocator () throw () { } - - inline pointer adress (reference r) { - return &r; - } - - inline const_pointer adress (const_reference r) const { - return &r; - } - - inline pointer allocate (size_type n) { - void * p = nullptr; - posix_memalign(&p,N,n*sizeof(value_type)); - return (pointer)p; - } - - inline void deallocate (pointer p, size_type) { - free(p); - } - - inline void construct (pointer p, const value_type & wert) { - new (p) value_type (wert); - } - - inline void destroy (pointer p) { - p->~value_type (); - } - - inline size_type max_size () const throw () { - return size_type (-1) / sizeof (value_type); - } - - template - struct rebind { - typedef AlignmentAllocator other; - }; - - bool operator!=(const AlignmentAllocator& other) const { - return !(*this == other); - } - - // Returns true if and only if storage allocated from *this - // can be deallocated from other, and vice versa. - // Always returns true for stateless allocators. - bool operator==(const AlignmentAllocator& other) const { - return true; - } -}; - class DAClusterizerInZT_vect final : public TrackClusterizerInZ { public: - template - using AlignedVector = std::vector >; - + // Internal data structure to struct track_t { @@ -132,39 +59,39 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { _pi = &pi.front(); } - double * __restrict__ _z __attribute__ ((aligned (16))); // z-coordinate at point of closest approach to the beamline - double * __restrict__ _t __attribute__ ((aligned (16))); // t-coordinate at point of closest approach to the beamline - double * __restrict__ _dz2 __attribute__ ((aligned (16))); // square of the error of z(pca) - double * __restrict__ _dt2 __attribute__ ((aligned (16))); // square of the error of t(pca) - double * __restrict__ _errsum __attribute__ ((aligned (16))); // sum of squares of the pca errors + double * __restrict__ _z; // z-coordinate at point of closest approach to the beamline + double * __restrict__ _t; // t-coordinate at point of closest approach to the beamline + double * __restrict__ _dz2; // square of the error of z(pca) + double * __restrict__ _dt2; // square of the error of t(pca) + double * __restrict__ _errsum; // sum of squares of the pca errors - double * __restrict__ _Z_sum __attribute__ ((aligned (16))); // Z[i] for DA clustering - double * __restrict__ _pi __attribute__ ((aligned (16))); // track weight + double * __restrict__ _Z_sum; // Z[i] for DA clustering + double * __restrict__ _pi; // track weight - AlignedVector z; // z-coordinate at point of closest approach to the beamline - AlignedVector t; // t-coordinate at point of closest approach to the beamline - AlignedVector dz2; // square of the error of z(pca) - AlignedVector dt2; // square of the error of t(pca) - AlignedVector errsum; // sum of squares of the pca errors - AlignedVector Z_sum; // Z[i] for DA clustering - AlignedVector pi; // track weight + std::vector z; // z-coordinate at point of closest approach to the beamline + std::vector t; // t-coordinate at point of closest approach to the beamline + std::vector dz2; // square of the error of z(pca) + std::vector dt2; // square of the error of t(pca) + std::vector errsum; // sum of squares of the pca errors + std::vector Z_sum; // Z[i] for DA clustering + std::vector pi; // track weight std::vector< const reco::TransientTrack* > tt; // a pointer to the Transient Track }; struct vertex_t { - AlignedVector z; // z coordinate - AlignedVector t; // t coordinate - AlignedVector pk; // vertex weight for "constrained" clustering + std::vector z; // z coordinate + std::vector t; // t coordinate + std::vector pk; // vertex weight for "constrained" clustering // --- temporary numbers, used during update - AlignedVector ei_cache; - AlignedVector ei; - AlignedVector sw; - AlignedVector swz; - AlignedVector swt; - AlignedVector se; - AlignedVector swE; + std::vector ei_cache; + std::vector ei; + std::vector sw; + std::vector swz; + std::vector swt; + std::vector se; + std::vector swE; unsigned int GetSize() const @@ -250,17 +177,17 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { } - double * __restrict__ _z __attribute__ ((aligned (16))); - double * __restrict__ _t __attribute__ ((aligned (16))); - double * __restrict__ _pk __attribute__ ((aligned (16))); - - double * __restrict__ _ei_cache __attribute__ ((aligned (16))); - double * __restrict__ _ei __attribute__ ((aligned (16))); - double * __restrict__ _sw __attribute__ ((aligned (16))); - double * __restrict__ _swz __attribute__ ((aligned (16))); - double * __restrict__ _swt __attribute__ ((aligned (16))); - double * __restrict__ _se __attribute__ ((aligned (16))); - double * __restrict__ _swE __attribute__ ((aligned (16))); + double * __restrict__ _z; + double * __restrict__ _t; + double * __restrict__ _pk; + + double * __restrict__ _ei_cache; + double * __restrict__ _ei; + double * __restrict__ _sw; + double * __restrict__ _swz; + double * __restrict__ _swt; + double * __restrict__ _se; + double * __restrict__ _swE; }; diff --git a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py index 32f0b455fd762..0f6775ff8006a 100644 --- a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py +++ b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py @@ -1,16 +1,5 @@ import FWCore.ParameterSet.Config as cms -DA2DParameters = cms.PSet( - algorithm = cms.string("DA2D"), - TkDAClusParameters = cms.PSet( - coolingFactor = cms.double(0.6), # moderate annealing speed - Tmin = cms.double(4.), # end of annealing - vertexSize = cms.double(0.006), # ~ resolution / sqrt(Tmin) - d0CutOff = cms.double(3.), # downweight high IP tracks - dzCutOff = cms.double(4.) # outlier rejection after freeze-out (T -DAClusterizerInZT::fill( const vector & tracks )const{ - // prepare track data for clustering - vector tks; - tks.reserve(tracks.size()); - for(vector::const_iterator it=tracks.begin(); it!=tracks.end(); it++){ - track_t t; - t.pi = 1.; - auto tsPCA = (*it).stateAtBeamLine().trackStateAtPCA(); - t.z = tsPCA.position().z(); - t.t = it->timeExt(); // the time - - if (std::abs(t.z) > 1000.) continue; - auto const & t_mom = tsPCA.momentum(); - // get the beam-spot - reco::BeamSpot beamspot = (*it).stateAtBeamLine().beamSpot(); - t.dz2 = - sqr((*it).track().dzError()) // track errror - + (sqr(beamspot.BeamWidthX()*t_mom.x())+sqr(beamspot.BeamWidthY()*t_mom.y()))*sqr(t_mom.z())/sqr(t_mom.perp2()) // beam spot width - + sqr(vertexSize_); // intrinsic vertex size, safer for outliers and short lived decays - //t.dz2 = 1./ t_dz2; - - t.dtz = 0.; - t.dt2 =sqr((*it).dtErrorExt()) + sqr(vertexSizeTime); // the ~injected~ timing error, need to add a small minimum vertex size in time - - if (d0CutOff_>0){ - Measurement1D IP = (*it).stateAtBeamLine().transverseImpactParameter();// error constains beamspot - t.pi=1./(1.+std::exp(sqr(IP.value()/IP.error()) - sqr(d0CutOff_))); // reduce weight for high ip tracks - }else{ - t.pi=1.; - } - t.tt=&(*it); - t.zi=1.; - if( edm::isFinite(t.pi) && t.pi >= std::numeric_limits::epsilon() ) { - tks.push_back(t); - } - } - - tks.shrink_to_fit(); - return tks; -} - - - - - -double DAClusterizerInZT::e_ik(const track_t & t, const vertex_t &k ) const { - return sqr(t.z-k.z)/t.dz2 + sqr(t.t - k.t)/t.dt2; -} - - - - - - -double DAClusterizerInZT::update( double beta, - vector & tks, - vector & y, - const double rho0 ) const { - // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise - // returns the squared sum of changes of vertex positions - - unsigned int nt=tks.size(); - - //initialize sums - double sumpi = 0.; - for(vector::iterator k=y.begin(); k!=y.end(); k++){ - k->sw = 0.; k->swz = 0.; k->swt = 0.; k->se = 0.; - k->swE = 0.; k->tC=0.; - } - - - // loop over tracks - for(unsigned int i=0; i::iterator k=y.begin(); k!=y.end(); k++){ - k->ei = std::exp(-beta*e_ik(tks[i],*k));// cache exponential for one track at a time - Zi += k->pk * k->ei; - } - tks[i].zi=Zi; - sumpi += tks[i].pi; - - // normalization - if (tks[i].zi>0){ - // accumulate weighted z and weights for vertex update - for(vector::iterator k=y.begin(); k!=y.end(); k++){ - double zratio = k->pk*k->ei/Zi; - - k->se += tks[i].pi*zratio/k->pk; - double w = tks[i].pi * zratio /( tks[i].dz2 + tks[i].dt2 ); - - k->sw += w; - k->swz += w * tks[i].z; - k->swt += w * tks[i].t; - k->swE += w * e_ik(tks[i],*k); - } - } - } // end of track loop - - - // now update z - double delta=0; - for(vector::iterator k=y.begin(); k!=y.end(); k++){ - if ( k->sw > 0){ - const double znew=k->swz/k->sw; - const double tnew=k->swt/k->sw; - delta += sqr(k->z-znew) + sqr(k->t-tnew); - k->z = znew; - k->t = tnew; - k->tC = 2*k->swE/k->sw; - }else{ - edm::LogInfo("sumw") << "invalid sum of weights in fit: " << k->sw << endl; - if(verbose_){cout << " a cluster melted away ? pk=" << k->pk << " sumw=" << k->sw << endl;} - k->tC = (rho0 == 0. ? -1 : 0); - } - if(rho0 == 0.) k->pk = k->pk * k->se / sumpi; - } - - // return how much the prototypes moved - return delta; -} - -bool DAClusterizerInZT::merge(vector & y, int nt)const{ - // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise - - if(y.size()<2) return false; - - for(vector::iterator k=y.begin(); (k+1)!=y.end(); k++){ - if( std::abs( (k+1)->z - k->z ) < epsilon && - std::abs( (k+1)->t - k->t ) < epsilon ){ // with fabs if only called after freeze-out (splitAll() at highter T) - double rho = k->pk + (k+1)->pk; - if(rho>0){ - k->z = ( k->pk * k->z + (k+1)->z * (k+1)->pk)/rho; - k->t = ( k->pk * k->t + (k+1)->t * (k+1)->pk)/rho; - }else{ - k->z = 0.5*(k->z + (k+1)->z); - k->t = 0.5*(k->t + (k+1)->t); - } - k->pk = rho; - - y.erase(k+1); - return true; - } - } - - return false; -} - -bool DAClusterizerInZT::merge(vector & y, double & beta)const{ - // merge clusters that collapsed or never separated, - // only merge if the estimated critical temperature of the merged vertex is below the current temperature - // return true if vertices were merged, false otherwise - if(y.size()<2) return false; - - for(vector::iterator k=y.begin(); (k+1)!=y.end(); k++){ - if ( std::abs((k+1)->z - k->z) < 2*epsilon && - std::abs((k+1)->t - k->t) < 2*epsilon ) { - double rho=k->pk + (k+1)->pk; - double swE=k->swE+(k+1)->swE - k->pk * (k+1)->pk / rho * ( sqr((k+1)->z - k->z) + - sqr((k+1)->t - k->t) ); - double Tc=2*swE/(k->sw+(k+1)->sw); - - if(Tc*beta<1){ - if(rho>0){ - k->z = ( k->pk * k->z + (k+1)->z * (k+1)->pk)/rho; - k->t = ( k->pk * k->t + (k+1)->t * (k+1)->pk)/rho; - }else{ - k->z = 0.5*(k->z + (k+1)->z); - k->t = 0.5*(k->t + (k+1)->t); - } - k->pk = rho; - k->sw += (k+1)->sw; - k->swE = swE; - k->tC = Tc; - y.erase(k+1); - return true; - } - } - } - - return false; -} - -bool DAClusterizerInZT::purge(vector & y, vector & tks, double & rho0, const double beta)const{ - // eliminate clusters with only one significant/unique track - if(y.size()<2) return false; - - unsigned int nt=tks.size(); - double sumpmin=nt; - vector::iterator k0=y.end(); - for(vector::iterator k=y.begin(); k!=y.end(); k++){ - int nUnique=0; - double sump=0; - double pmax=k->pk/(k->pk+rho0*exp(-beta*dzCutOff_*dzCutOff_)); - for(unsigned int i=0; i 0){ - double p = k->pk * std::exp(-beta*e_ik(tks[i],*k)) / tks[i].zi ; - sump+=p; - if( (p > 0.9*pmax) && (tks[i].pi>0) ){ nUnique++; } - } - } - - if((nUnique<2)&&(sumpz << "," << k0->t - << " with sump=" << sumpmin << endl;} - y.erase(k0); - return true; - }else{ - return false; - } -} - -double DAClusterizerInZT::beta0( double betamax, - vector & tks, - vector & y ) const { - - double T0=0; // max Tc for beta=0 - // estimate critical temperature from beta=0 (T=inf) - unsigned int nt=tks.size(); - - for(vector::iterator k=y.begin(); k!=y.end(); k++){ - - // vertex fit at T=inf - double sumwz=0.; - double sumwt=0.; - double sumw=0.; - for(unsigned int i=0; iz = sumwz/sumw; - k->t = sumwt/sumw; - - // estimate Tcrit, eventually do this in the same loop - double a=0, b=0; - for(unsigned int i=0; iz); - double dt = tks[i].t-(k->t); - double w = tks[i].pi/(tks[i].dz2 + tks[i].dt2); - a += w*(sqr(dx)/tks[i].dz2 + sqr(dt)/tks[i].dt2); - b += w; - } - double Tc= 2.*a/b; // the critical temperature of this vertex - if(Tc>T0) T0=Tc; - }// vertex loop (normally there should be only one vertex at beta=0) - - if (T0>1./betamax){ - return betamax/pow(coolingFactor_, int(std::log(T0*betamax)*logCoolingFactor_)-1 ); - }else{ - // ensure at least one annealing step - return betamax/coolingFactor_; - } -} - -bool DAClusterizerInZT::split( double beta, - vector & tks, - vector & y, - double threshold ) const { - // split only critical vertices (Tc >~ T=1/beta <==> beta*Tc>~1) - // an update must have been made just before doing this (same beta, no merging) - // returns true if at least one cluster was split - bool split=false; - - // avoid left-right biases by splitting highest Tc first - - std::vector > critical; - for(unsigned int ik=0; ik 1.){ - critical.push_back( make_pair(y[ik].tC, ik)); - } - } - std::stable_sort(critical.begin(), critical.end(), std::greater >() ); - - for(unsigned int ic=0; ic0){ - //sumpi+=tks[i].pi; - double p=y[ik].pk * exp(-beta*e_ik(tks[i],y[ik])) / tks[i].zi*tks[i].pi; - double w=p/(tks[i].dz2 + tks[i].dt2); - if(tks[i].z < y[ik].z){ - p1+=p; z1+=w*tks[i].z; t1+=w*tks[i].t; w1+=w; - }else{ - p2+=p; z2+=w*tks[i].z; t2+=w*tks[i].t; w2+=w; - } - } - } - if(w1>0){ z1=z1/w1; t1=t1/w1;} else{ z1=y[ik].z-epsilon; t1=y[ik].t-epsilon; } - if(w2>0){ z2=z2/w2; t2=t2/w2;} else{ z2=y[ik].z+epsilon; t2=y[ik].t+epsilon; } - - // reduce split size if there is not enough room - if( ( ik > 0 ) && ( y[ik-1].z>=z1 ) ){ z1=0.5*(y[ik].z+y[ik-1].z); t1=0.5*(y[ik].t+y[ik-1].t); } - if( ( ik+1 < y.size()) && ( y[ik+1].z<=z2 ) ){ z2=0.5*(y[ik].z+y[ik+1].z); t2=0.5*(y[ik].t+y[ik+1].t); } - - // split if the new subclusters are significantly separated - if( std::abs(z2-z1) > epsilon || std::abs(t2-t1) > epsilon ){ - split=true; - vertex_t vnew; - vnew.pk = p1*y[ik].pk/(p1+p2); - y[ik].pk= p2*y[ik].pk/(p1+p2); - vnew.z = z1; - vnew.t = t1; - y[ik].z = z2; - y[ik].t = t2; - y.insert(y.begin()+ik, vnew); - - // adjust remaining pointers - for(unsigned int jc=ic; jcik) {critical[jc].second++;} - } - } - } - - // stable_sort(y.begin(), y.end(), clusterLessZ); - return split; -} - - - - - -void DAClusterizerInZT::splitAll( vector & y ) const { - - constexpr double zsep=2*epsilon; // split vertices that are isolated by at least zsep (vertices that haven't collapsed) - constexpr double tsep=2*epsilon; // check t as well - - vector y1; - - for(vector::iterator k=y.begin(); k!=y.end(); k++){ - if ( ( (k==y.begin())|| (k-1)->z < k->z - zsep) && (((k+1)==y.end() )|| (k+1)->z > k->z + zsep)) { - // isolated prototype, split - vertex_t vnew; - vnew.z = k->z - epsilon; - vnew.t = k->t - epsilon; - (*k).z = k->z + epsilon; - (*k).t = k->t + epsilon; - vnew.pk= 0.5* (*k).pk; - (*k).pk= 0.5* (*k).pk; - y1.push_back(vnew); - y1.push_back(*k); - - }else if( y1.empty() || (y1.back().z < k->z -zsep) || (y1.back().t < k->t - tsep) ){ - y1.push_back(*k); - }else{ - y1.back().z -= epsilon; - y1.back().t -= epsilon; - k->z += epsilon; - k->t += epsilon; - y1.push_back(*k); - } - }// vertex loop - - y=y1; -} - - -DAClusterizerInZT::DAClusterizerInZT(const edm::ParameterSet& conf) : - verbose_(conf.getUntrackedParameter("verbose", false)), - useTc_(true), - vertexSize_(conf.getParameter("vertexSize")), - maxIterations_(100), - coolingFactor_(std::sqrt(conf.getParameter("coolingFactor"))), - betamax_(0.1), - betastop_(1.0), - dzCutOff_(conf.getParameter("dzCutOff")), - d0CutOff_(conf.getParameter("d0CutOff")), - dtCutOff_(dtCutOff) -{ - - double Tmin = conf.getParameter("Tmin")*std::sqrt(2.0);// scale up by sqrt(D=2) - if (Tmin==0){ - edm::LogWarning("DAClusterizerInZT") << "DAClusterizerInZT: invalid Tmin" << Tmin << " reset do default " << 1./betamax_ << endl; - }else{ - betamax_ = 1./Tmin; - } - - // for testing, negative cooling factor: revert to old splitting scheme - if(coolingFactor_<0){ - coolingFactor_=-coolingFactor_; useTc_=false; - } - - logCoolingFactor_ = 1.0/std::log(coolingFactor_); -} - - -void DAClusterizerInZT::dump(const double beta, const vector & y, const vector & tks0, int verbosity)const{ - - // copy and sort for nicer printout - vector tks; - for(vector::const_iterator t=tks0.begin(); t!=tks0.end(); t++){tks.push_back(*t); } - std::stable_sort(tks.begin(), tks.end(), recTrackLessZ1); - - cout << "-----DAClusterizerInZT::dump ----" << endl; - cout << " beta=" << beta << " betamax= " << betamax_ << endl; - cout << " z= "; - cout.precision(4); - for(vector::const_iterator k=y.begin(); k!=y.end(); k++){ - cout << setw(8) << fixed << k->z; - } - cout << endl << " t= "; - for(vector::const_iterator k=y.begin(); k!=y.end(); k++){ - cout << setw(8) << fixed << k->t; - } - cout << endl << "T=" << setw(15) << 1./beta <<" Tc= "; - for(vector::const_iterator k=y.begin(); k!=y.end(); k++){ - cout << setw(8) << fixed << k->tC ; - } - - cout << endl << " pk="; - double sumpk=0; - for(vector::const_iterator k=y.begin(); k!=y.end(); k++){ - cout << setw(8) << setprecision(3) << fixed << k->pk; - sumpk+=k->pk; - } - cout << endl; - - if(verbosity>0){ - double E=0, F=0; - cout << endl; - cout << "---- z +/- dz t +/- dt ip +/-dip pt phi eta weights ----" << endl; - cout.precision(4); - for(unsigned int i=0; i0){ F-=log(tks[i].zi)/beta;} - double tz= tks[i].z; - double tt= tks[i].t; - cout << setw (3)<< i << ")" << setw (8) << fixed << setprecision(4)<< tz << " +/-" << setw (6)<< sqrt(tks[i].dz2) - << setw(8) << fixed << setprecision(4) << tt << " +/-" << setw(6) << std::sqrt(tks[i].dt2) ; - - if(tks[i].tt->track().quality(reco::TrackBase::highPurity)){ cout << " *";}else{cout <<" ";} - if(tks[i].tt->track().hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)){cout <<"+";}else{cout << "-";} - - cout << setw(1) << tks[i].tt->track().hitPattern().pixelBarrelLayersWithMeasurement(); // see DataFormats/TrackReco/interface/HitPattern.h - cout << setw(1) << tks[i].tt->track().hitPattern().pixelEndcapLayersWithMeasurement(); - cout << setw(1) << hex << tks[i].tt->track().hitPattern().trackerLayersWithMeasurement()-tks[i].tt->track().hitPattern().pixelLayersWithMeasurement() <track().pt()*tks[i].tt->track().charge(); - cout << " " << setw(5) << setprecision(2) << tks[i].tt->track().phi() - << " " << setw(5) << setprecision(2) << tks[i].tt->track().eta() ; - - double sump=0.; - for(vector::const_iterator k=y.begin(); k!=y.end(); k++){ - if((tks[i].pi>0)&&(tks[i].zi>0)){ - //double p=pik(beta,tks[i],*k); - double p=k->pk * std::exp(-beta*e_ik(tks[i],*k)) / tks[i].zi; - if( p > 0.0001){ - cout << setw (8) << setprecision(3) << p; - }else{ - cout << " . "; - } - E+=p*e_ik(tks[i],*k); - sump+=p; - }else{ - cout << " "; - } - } - cout << endl; - } - cout << endl << "T=" << 1/beta << " E=" << E << " n="<< y.size() << " F= " << F << endl << "----------" << endl; - } -} - - - - - -vector< TransientVertex > -DAClusterizerInZT::vertices(const vector & tracks, const int verbosity) -const -{ - - vector tks=fill(tracks); - unsigned int nt=tks.size(); - double rho0=0.0; // start with no outlier rejection - - vector< TransientVertex > clusters; - if (tks.empty()) return clusters; - - vector y; // the vertex prototypes - - // initialize:single vertex at infinite temperature - vertex_t vstart; - vstart.z=0.; - vstart.t=0.; - vstart.pk=1.; - y.push_back(vstart); - int niter=0; // number of iterations - - - // estimate first critical temperature - double beta=beta0(betamax_, tks, y); - niter=0; - while((update(beta, tks,y)>1.e-6) && (niter++ < maxIterations_)){ } - - // annealing loop, stop when T1/Tmin) - while(beta1.e-6) && (niter++ < maxIterations_)){ } - - } - - if(useTc_){ - // last round of splitting, make sure no critical clusters are left - update(beta, tks,y); - while(merge(y,beta)){update(beta, tks,y);} - unsigned int ntry=0; - while( split(beta, tks,y,1.) && (ntry++<10) ){ - niter=0; - while((update(beta, tks,y)>1.e-6) && (niter++ < maxIterations_)){} - merge(y,beta); - update(beta, tks,y); - } - }else{ - // merge collapsed clusters - while(merge(y,beta)){update(beta, tks,y);} - if(verbose_ ){ cout << "dump after 1st merging " << endl; dump(beta,y,tks,2);} - } - - - - - // switch on outlier rejection - rho0=1./nt; - for(vector::iterator k=y.begin(); k!=y.end(); k++){ k->pk =1.; } // democratic - niter=0; - while( (update(beta, tks,y,rho0) > 1.e-8) && (niter++ < maxIterations_) ){ } - if(verbose_ ){ cout << "rho0=" << rho0 << " niter=" << niter << endl; dump(beta,y,tks,2);} - - - // merge again (some cluster split by outliers collapse here) - while(merge(y,tks.size())){} - if(verbose_ ){ cout << "dump after 2nd merging " << endl; dump(beta,y,tks,2);} - - - // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices - while(beta<=betastop_){ - while(purge(y,tks,rho0, beta)){ - niter=0; - while((update(beta, tks, y, rho0) > 1.e-6) && (niter++ < maxIterations_)){ } - } - beta/=coolingFactor_; - niter=0; - while((update(beta, tks, y, rho0) > 1.e-6) && (niter++ < maxIterations_)){ } - } - - if(verbose_){ - cout << "Final result, rho0=" << rho0 << endl; - dump(beta,y,tks,2); - } - - // ensure correct normalization of probabilities, should make double assginment reasonably impossible - for(unsigned int i=0; i::iterator k=y.begin(); k!=y.end(); k++){ - tks[i].zi += k->pk * exp(-beta*e_ik(tks[i],*k)); - } - } - - - for(vector::iterator k=y.begin(); k!=y.end(); k++){ - GlobalPoint pos(0, 0, k->z); - double time = k->t; - vector< reco::TransientTrack > vertexTracks; - //double max_track_time_err2 = 0; - double mean = 0.; - double expv_x2 = 0.; - double normw = 0.; - for(unsigned int i=0; i0){ - double p = k->pk * exp(-beta*e_ik(tks[i],*k)) / tks[i].zi; - if( (tks[i].pi>0) && ( p > 0.5 ) ){ - //std::cout << "pushing back " << i << ' ' << tks[i].tt << std::endl; - vertexTracks.push_back(*(tks[i].tt)); tks[i].zi=0; - mean += tks[i].t*invdt*p; - expv_x2 += tks[i].t*tks[i].t*invdt*p; - normw += invdt*p; - } // setting Z=0 excludes double assignment - } - } - mean = mean/normw; - expv_x2 = expv_x2/normw; - const double time_var = expv_x2 - mean*mean; - const double crappy_error_guess = std::sqrt(time_var); - GlobalError dummyErrorWithTime(0, - 0,0, - 0,0,0, - 0,0,0,crappy_error_guess); - TransientVertex v(pos, time, dummyErrorWithTime, vertexTracks, 5); - clusters.push_back(v); - } - - - return clusters; - -} - - - - - -vector< vector > -DAClusterizerInZT::clusterize(const vector & tracks) - const -{ - if(verbose_) { - cout << "###################################################" << endl; - cout << "# DAClusterizerInZT::clusterize nt="< > clusters; - vector< TransientVertex > pv=vertices(tracks); - - if(verbose_){ cout << "# DAClusterizerInZT::clusterize pv.size="< aCluster=pv.begin()->originalTracks(); - - if( verbose_ ) { - std::cout << '\t' << 0; - std::cout << ' ' << pv.begin()->position() << ' ' << pv.begin()->time() << std::endl; - } - - for(vector::iterator k=pv.begin()+1; k!=pv.end(); k++){ - if( verbose_ ) { - std::cout << '\t' << std::distance(pv.begin(),k); - std::cout << ' ' << k->position() << ' ' << k->time() << std::endl; - } - if ( std::abs(k->position().z() - (k-1)->position().z()) > (2*vertexSize_) || - std::abs(k->time() - (k-1)->time()) > 2*vertexSizeTime ) { - // close a cluster - clusters.push_back(aCluster); - aCluster.clear(); - } - for(unsigned int i=0; ioriginalTracks().size(); i++){ - aCluster.push_back( k->originalTracks().at(i)); - } - - } - clusters.push_back(aCluster); - - if(verbose_) { std::cout << "# DAClusterizerInZT::clusterize clusters.size="< ("zdumpwidth", 20.); // configurable parameters - double Tmin = conf.getParameter ("Tmin"); - double Tpurge = conf.getParameter ("Tpurge"); - double Tstop = conf.getParameter ("Tstop"); + double Tmin = conf.getParameter ("Tmin")*std::sqrt(2.0); + double Tpurge = conf.getParameter ("Tpurge")*std::sqrt(2.0); + double Tstop = conf.getParameter ("Tstop")*std::sqrt(2.0); vertexSize_ = conf.getParameter ("vertexSize"); coolingFactor_ = conf.getParameter ("coolingFactor"); useTc_=true; @@ -41,12 +41,14 @@ DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) { coolingFactor_=-coolingFactor_; useTc_=false; } + coolingFactor_ = std::sqrt(coolingFactor_); d0CutOff_ = conf.getParameter ("d0CutOff"); dzCutOff_ = conf.getParameter ("dzCutOff"); uniquetrkweight_ = conf.getParameter("uniquetrkweight"); zmerge_ = conf.getParameter("zmerge"); tmerge_ = conf.getParameter("tmerge"); +#ifdef VI_DEBUG if(verbose_){ std::cout << "DAClusterizerinZT_vect: mintrkweight = " << mintrkweight_ << std::endl; std::cout << "DAClusterizerinZT_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl; @@ -61,6 +63,7 @@ DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) { std::cout << "DAClusterizerinZT_vect: dzCutOff = " << dzCutOff_ << std::endl; std::cout << "DAClusterizerinZT_vect: dtCutoff = " << dtCutOff << std::endl; } +#endif if (Tmin == 0) { @@ -136,9 +139,11 @@ DAClusterizerInZT_vect::fill(const vector & tracks) const } tks.ExtractRaw(); +#ifdef VI_DEBUG if (verbose_) { std::cout << "Track count " << tks.GetSize() << std::endl; } +#endif return tks; } @@ -328,7 +333,9 @@ bool DAClusterizerInZT_vect::merge(vertex_t & y, double & beta)const{ double Tc = 2*swE / (y._sw[k]+y._sw[k+1]); if(Tc*beta < 1){ +#ifdef VI_DEBUG if(verbose_){ std::cout << "merging (" << y._z[k + 1] << ',' << y._t[k + 1] << ") and (" << y._z[k] << ',' << y._t[k] << ") Tc = " << Tc << " sw = " << y._sw[k]+y._sw[k+1] < 0){ y._z[k] = (y._pk[k]*y._z[k] + y._pk[k+1]*y._z[k + 1])/rho; y._t[k] = (y._pk[k]*y._t[k] + y._pk[k+1]*y._t[k + 1])/rho; @@ -366,10 +373,10 @@ DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const int nUnique = 0; double sump = 0; - AlignedVector inverse_zsums(nt), arg_cache(nt), eik_cache(nt); - double * __restrict__ _inverse_zsums __attribute__ ((aligned (16))); - double * __restrict__ _arg_cache __attribute__ ((aligned (16))); - double * __restrict__ _eik_cache __attribute__ ((aligned (16))); + std::vector inverse_zsums(nt), arg_cache(nt), eik_cache(nt); + double * __restrict__ _inverse_zsums; + double * __restrict__ _arg_cache; + double * __restrict__ _eik_cache; _inverse_zsums = inverse_zsums.data(); _arg_cache = arg_cache.data(); _eik_cache = eik_cache.data(); @@ -409,12 +416,14 @@ DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const } if (k0 != nv) { +#ifdef VI_DEBUG if (verbose_) { std::cout << "eliminating prototype at " << std::setw(10) << std::setprecision(4) << y._z[k0] << " with sump=" << sumpmin << " rho*nt =" << y._pk[k0]*nt << endl; } +#endif y.RemoveItem(k0); return true; } else { @@ -461,11 +470,13 @@ DAClusterizerInZT_vect::beta0(double betamax, track_t const & tks, vertex_t con if (Tc > T0) T0 = Tc; }// vertex loop (normally there should be only one vertex at beta=0) +#ifdef VI_DEBUG if(verbose_){ std::cout << "DAClustrizerInZT_vect.beta0: Tc = " << T0 << std::endl; int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_)); std::cout << "DAClustrizerInZT_vect.beta0: nstep = " << coolingsteps << std::endl; } +#endif if (T0 > 1. / betamax) { @@ -550,6 +561,7 @@ DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, do t2 = 0.5*y._t[k] + 0.5*y._t[k+1]; } +#ifdef VI_DEBUG if(verbose_){ if (std::fabs(y._z[k] - zdumpcenter_) < zdumpwidth_){ std::cout << " T= " << std::setw(8) << 1./beta @@ -564,6 +576,7 @@ DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, do } } } +#endif // split if the new subclusters are significantly separated if( std::fabs(z2-z1) > epsilonz || std::fabs(t2-t1) > epsilont){ @@ -597,10 +610,12 @@ void DAClusterizerInZT_vect::splitAll( vertex_t & y) const { constexpr double tsep = 2 * epsilont; vertex_t y1; +#ifdef VI_DEBUG if (verbose_) { std::cout << "Before Split "<< std::endl; y.DebugOut(); } +#endif for (unsigned int k = 0; k < nv; k++) { if ( @@ -638,10 +653,12 @@ void DAClusterizerInZT_vect::splitAll( vertex_t & y) const { y = y1; y.ExtractRaw(); +#ifdef VI_DEBUG if (verbose_) { std::cout << "After split " << std::endl; y.DebugOut(); } +#endif } @@ -667,7 +684,9 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co // estimate first critical temperature double beta = beta0(betamax_, tks, y); +#ifdef VI_DEBUG if ( verbose_) std::cout << "Beta0 is " << beta << std::endl; +#endif niter = 0; while ((update(beta, tks, y, false, rho0) > 1.e-6) && @@ -699,7 +718,9 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co if(useTc_){ //last round of splitting, make sure no critical clusters are left +#ifdef VI_DEBUG if(verbose_){ std::cout << "last spliting at " << 1./beta << std::endl; } +#endif update(beta, tks,y, false, rho0);// make sure Tc is up-to-date while(merge(y,beta)){update(beta, tks,y, false, rho0);} unsigned int ntry=0; @@ -708,10 +729,12 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co niter=0; while((update(beta, tks,y, false, rho0)>1.e-6) && (niter++ < maxIterations_)){} while(merge(y,beta)){update(beta, tks,y, false, rho0);} +#ifdef VI_DEBUG if(verbose_){ std::cout << "after final splitting, try " << ntry << std::endl; dump(beta, y, tks, 2); } +#endif // relax splitting a bit to reduce multiple split-merge cycles of the same cluster threshold *= 1.1; } @@ -720,11 +743,13 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co while(merge(y,beta)){update(beta, tks,y, false, rho0);} } +#ifdef VI_DEBUG if (verbose_) { update(beta, tks,y, false, rho0); std::cout << "dump after 1st merging " << endl; dump(beta, y, tks, 2); } +#endif // switch on outlier rejection at T=Tmin @@ -735,17 +760,21 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co niter=0; while ((update(beta, tks, y, true, rho0) > 1.e-8) && (niter++ < maxIterations_)) {}; +#ifdef VI_DEBUG if (verbose_) { std::cout << "dump after noise-suppression, rho0=" << rho0 << " niter = " << niter << endl; dump(beta, y, tks, 2); } +#endif // merge again (some cluster split by outliers collapse here) while (merge(y, beta)) {update(beta, tks, y, true, rho0); } +#ifdef VI_DEBUG if (verbose_) { std::cout << "dump after merging " << endl; dump(beta, y, tks, 2); } +#endif // go down to the purging temperature (if it is lower than tmin) while( beta < betapurge_ ){ @@ -762,11 +791,13 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co } } +#ifdef VI_DEBUG if (verbose_) { update(beta, tks,y, true, rho0); std::cout << " after purging " << std:: endl; dump(beta, y, tks, 2); } +#endif // optionally cool some more without doing anything, to make the assignment harder while( beta < betastop_ ){ @@ -775,11 +806,12 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co while ((update(beta, tks, y, true, rho0) > 1.e-8) && (niter++ < maxIterations_)) {} } - +#ifdef VI_DEBUG if (verbose_) { std::cout << "Final result, rho0=" << std::scientific << rho0 << endl; dump(beta, y, tks, 2); } +#endif // select significant tracks and use a TransientVertex as a container GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01); @@ -825,19 +857,24 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co vector > DAClusterizerInZT_vect::clusterize( const vector & tracks) const { +#ifdef VI_DEBUG if (verbose_) { std::cout << "###################################################" << endl; std::cout << "# vectorized DAClusterizerInZT_vect::clusterize nt=" << tracks.size() << endl; std::cout << "###################################################" << endl; } +#endif vector > clusters; vector && pv = vertices(tracks); +#ifdef VI_DEBUG if (verbose_) { std::cout << "# DAClusterizerInZT::clusterize pv.size=" << pv.size() << endl; } +#endif + if (pv.size() == 0) { return clusters; } @@ -851,9 +888,11 @@ vector > DAClusterizerInZT_vect::clusterize( if (aCluster.size()>1){ clusters.push_back(aCluster); }else{ +#ifdef VI_DEBUG if(verbose_){ std::cout << " one track cluster at " << k->position().z() << " suppressed" << std::endl; } +#endif } aCluster.clear(); } From 3c27ddab224389fa61409af66a529d739540c351 Mon Sep 17 00:00:00 2001 From: lgray Date: Tue, 15 Aug 2017 13:08:12 -0500 Subject: [PATCH 06/16] oops --- .../PrimaryVertexProducer/interface/PrimaryVertexProducer.h | 1 - .../PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc | 4 ---- 2 files changed, 5 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h index e55489e1b8f84..546d1f312bf8f 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h @@ -41,7 +41,6 @@ #include "RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h" #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h" #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h" -#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT.h" #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h" #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h" //#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h" diff --git a/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc index 12cb6c8f9310b..6599da9825bfe 100644 --- a/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc +++ b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc @@ -49,10 +49,6 @@ PrimaryVertexProducer::PrimaryVertexProducer(const edm::ParameterSet& conf) // provide the vectorized version of the clusterizer, if supported by the build else if(clusteringAlgorithm == "DA_vect") { theTrackClusterizer = new DAClusterizerInZ_vect(conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); - } - else if( clusteringAlgorithm=="DA2D" ) { - theTrackClusterizer = new DAClusterizerInZT(conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); - f4D = true; } else if( clusteringAlgorithm=="DA2D_vect" ) { theTrackClusterizer = new DAClusterizerInZT_vect(conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); f4D = true; From 9e6d6db5d7853e6dc3b5e4f20473a3f6e7df4ac5 Mon Sep 17 00:00:00 2001 From: lgray Date: Thu, 31 Aug 2017 13:02:08 -0500 Subject: [PATCH 07/16] code review --- .../interface/DAClusterizerInZT_vect.h | 72 ++-- .../src/DAClusterizerInZT_vect.cc | 308 +++++++++--------- 2 files changed, 190 insertions(+), 190 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h index 3acc66d96e439..2726f59faed56 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -50,23 +50,23 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { // has to be called everytime the items are modified void ExtractRaw() { - _z = &z.front(); - _t = &t.front(); - _dz2 = &dz2.front(); - _dt2 = &dt2.front(); - _errsum = &errsum.front(); - _Z_sum = &Z_sum.front(); - _pi = &pi.front(); + z_ = &z.front(); + t_ = &t.front(); + dz2_ = &dz2.front(); + dt2_ = &dt2.front(); + errsum_ = &errsum.front(); + Z_sum_ = &Z_sum.front(); + pi_ = &pi.front(); } - double * __restrict__ _z; // z-coordinate at point of closest approach to the beamline - double * __restrict__ _t; // t-coordinate at point of closest approach to the beamline - double * __restrict__ _dz2; // square of the error of z(pca) - double * __restrict__ _dt2; // square of the error of t(pca) - double * __restrict__ _errsum; // sum of squares of the pca errors + double * z_; // z-coordinate at point of closest approach to the beamline + double * t_; // t-coordinate at point of closest approach to the beamline + double * dz2_; // square of the error of z(pca) + double * dt2_; // square of the error of t(pca) + double * errsum_; // sum of squares of the pca errors - double * __restrict__ _Z_sum; // Z[i] for DA clustering - double * __restrict__ _pi; // track weight + double * Z_sum_; // Z[i] for DA clustering + double * pi_; // track weight std::vector z; // z-coordinate at point of closest approach to the beamline std::vector t; // t-coordinate at point of closest approach to the beamline @@ -156,38 +156,38 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { for ( unsigned int i =0; i < GetSize(); ++ i) { - std::cout << " z = " << _z[i] << " t = " << _t[i] << " pk = " << _pk[i] << std::endl; + std::cout << " z = " << z_[i] << " t = " << t_[i] << " pk = " << pk_[i] << std::endl; } } // has to be called everytime the items are modified void ExtractRaw() { - _z = &z.front(); - _t = &t.front(); - _pk = &pk.front(); + z_ = &z.front(); + t_ = &t.front(); + pk_ = &pk.front(); - _ei = &ei.front(); - _sw = &sw.front(); - _swz = &swz.front(); - _swt = &swt.front(); - _se = &se.front(); - _swE = &swE.front(); - _ei_cache = &ei_cache.front(); + ei_ = &ei.front(); + sw_ = &sw.front(); + swz_ = &swz.front(); + swt_ = &swt.front(); + se_ = &se.front(); + swE_ = &swE.front(); + ei_cache_ = &ei_cache.front(); } - double * __restrict__ _z; - double * __restrict__ _t; - double * __restrict__ _pk; - - double * __restrict__ _ei_cache; - double * __restrict__ _ei; - double * __restrict__ _sw; - double * __restrict__ _swz; - double * __restrict__ _swt; - double * __restrict__ _se; - double * __restrict__ _swE; + double * z_; + double * t_; + double * pk_; + + double * ei_cache_; + double * ei_; + double * sw_; + double * swz_; + double * swt_; + double * se_; + double * swE_; }; diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index e8e3b3fd4e4ac..53268fe2a6bae 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -97,8 +97,8 @@ namespace { return vdt::fast_exp( inp ); } - inline void local_exp_list( double const * __restrict__ arg_inp, - double * __restrict__ arg_out, + inline void local_exp_list( double const *arg_inp, + double *arg_out, const unsigned arg_arr_size) { for(unsigned i=0; i!=arg_arr_size; ++i ) arg_out[i]=vdt::fast_exp(arg_inp[i]); } @@ -186,16 +186,16 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, track_t const& tracks, vertex_t const& vertices ) { - const auto track_z = tracks._z[itrack]; - const auto track_t = tracks._t[itrack]; - const auto botrack_dz2 = -beta*tracks._dz2[itrack]; - const auto botrack_dt2 = -beta*tracks._dt2[itrack]; + const auto track_z = tracks.z_[itrack]; + const auto track_t = tracks.t_[itrack]; + const auto botrack_dz2 = -beta*tracks.dz2_[itrack]; + const auto botrack_dt2 = -beta*tracks.dt2_[itrack]; // auto-vectorized for ( unsigned int ivertex = 0; ivertex < nv; ++ivertex) { - const auto mult_resz = track_z - vertices._z[ivertex]; - const auto mult_rest = track_t - vertices._t[ivertex]; - vertices._ei_cache[ivertex] = botrack_dz2 * ( mult_resz * mult_resz ) + botrack_dt2 * ( mult_rest * mult_rest ); + const auto mult_resz = track_z - vertices.z_[ivertex]; + const auto mult_rest = track_t - vertices.t_[ivertex]; + vertices.ei_cache_[ivertex] = botrack_dz2 * ( mult_resz * mult_resz ) + botrack_dt2 * ( mult_rest * mult_rest ); } }; @@ -203,7 +203,7 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, { double ZTemp = Z_init; for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) { - ZTemp += vertices._pk[ivertex] * vertices._ei[ivertex]; + ZTemp += vertices.pk_[ivertex] * vertices.ei_[ivertex]; } return ZTemp; }; @@ -211,48 +211,48 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, auto kernel_calc_normalization = [ beta, nv ] (const unsigned int track_num, track_t & tks_vec, vertex_t & y_vec ) { - auto tmp_trk_pi = tks_vec._pi[track_num]; - auto o_trk_Z_sum = 1./tks_vec._Z_sum[track_num]; - //auto o_trk_dz2 = tks_vec._dz2[track_num]; - //auto o_trk_dt2 = tks_vec._dt2[track_num]; - auto o_trk_err_sum = tks_vec._errsum[track_num]; - auto tmp_trk_z = tks_vec._z[track_num]; - auto tmp_trk_t = tks_vec._t[track_num]; + auto tmp_trk_pi = tks_vec.pi_[track_num]; + auto o_trk_Z_sum = 1./tks_vec.Z_sum_[track_num]; + //auto o_trk_dz2 = tks_vec.dz2_[track_num]; + //auto o_trk_dt2 = tks_vec.dt2_[track_num]; + auto o_trk_err_sum = tks_vec.errsum_[track_num]; + auto tmp_trk_z = tks_vec.z_[track_num]; + auto tmp_trk_t = tks_vec.t_[track_num]; auto obeta = -1./beta; // auto-vectorized for (unsigned int k = 0; k < nv; ++k) { // parens are important for numerical stability - y_vec._se[k] += tmp_trk_pi*( y_vec._ei[k] * o_trk_Z_sum ); - const auto w = tmp_trk_pi * (y_vec._pk[k] * y_vec._ei[k] * o_trk_Z_sum) * o_trk_err_sum; - y_vec._sw[k] += w; - y_vec._swz[k] += w * tmp_trk_z; - y_vec._swt[k] += w * tmp_trk_t; - y_vec._swE[k] += w * y_vec._ei_cache[k]*obeta; + y_vec.se_[k] += tmp_trk_pi*( y_vec.ei_[k] * o_trk_Z_sum ); + const auto w = tmp_trk_pi * (y_vec.pk_[k] * y_vec.ei_[k] * o_trk_Z_sum) * o_trk_err_sum; + y_vec.sw_[k] += w; + y_vec.swz_[k] += w * tmp_trk_z; + y_vec.swt_[k] += w * tmp_trk_t; + y_vec.swE_[k] += w * y_vec.ei_cache_[k]*obeta; } }; for (auto ivertex = 0U; ivertex < nv; ++ivertex) { - gvertices._se[ivertex] = 0.0; - gvertices._sw[ivertex] = 0.0; - gvertices._swz[ivertex] = 0.0; - gvertices._swt[ivertex] = 0.0; - gvertices._swE[ivertex] = 0.0; + gvertices.se_[ivertex] = 0.0; + gvertices.sw_[ivertex] = 0.0; + gvertices.swz_[ivertex] = 0.0; + gvertices.swt_[ivertex] = 0.0; + gvertices.swE_[ivertex] = 0.0; } // loop over tracks for (auto itrack = 0U; itrack < nt; ++itrack) { kernel_calc_exp_arg(itrack, gtracks, gvertices); - local_exp_list(gvertices._ei_cache, gvertices._ei, nv); + local_exp_list(gvertices.ei_cache_, gvertices.ei_, nv); - gtracks._Z_sum[itrack] = kernel_add_Z(gvertices); - if (edm::isNotFinite(gtracks._Z_sum[itrack])) gtracks._Z_sum[itrack] = 0.0; + gtracks.Z_sum_[itrack] = kernel_add_Z(gvertices); + if (edm::isNotFinite(gtracks.Z_sum_[itrack])) gtracks.Z_sum_[itrack] = 0.0; // used in the next major loop to follow - sumpi += gtracks._pi[itrack]; + sumpi += gtracks.pi_[itrack]; - if (gtracks._Z_sum[itrack] > 1.e-100){ + if (gtracks.Z_sum_[itrack] > 1.e-100){ kernel_calc_normalization(itrack, gtracks, gvertices); } } @@ -263,23 +263,23 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, double delta=0; // does not vectorizes for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex ) { - if (vertices._sw[ivertex] > 0.) { - //std::cout << " DA2D_vect sw = " << vertices._sw[ ivertex ] << ' ' << vertices._swz[ ivertex ] << ' ' << vertices._swt[ ivertex ] << std::endl; - auto znew = vertices._swz[ ivertex ] / vertices._sw[ ivertex ]; + if (vertices.sw_[ivertex] > 0.) { + //std::cout << " DA2D_vect sw = " << vertices.sw_[ ivertex ] << ' ' << vertices.swz_[ ivertex ] << ' ' << vertices.swt_[ ivertex ] << std::endl; + auto znew = vertices.swz_[ ivertex ] / vertices.sw_[ ivertex ]; // prevents from vectorizing if - delta += std::pow( vertices._z[ ivertex ] - znew, 2 ); - vertices._z[ ivertex ] = znew; - auto tnew = vertices._swt[ ivertex ] / vertices._sw[ ivertex ]; + delta += std::pow( vertices.z_[ ivertex ] - znew, 2 ); + vertices.z_[ ivertex ] = znew; + auto tnew = vertices.swt_[ ivertex ] / vertices.sw_[ ivertex ]; // prevents from vectorizing if - delta += std::pow( vertices._t[ ivertex ] - tnew, 2 ); - vertices._t[ ivertex ] = tnew; + delta += std::pow( vertices.t_[ ivertex ] - tnew, 2 ); + vertices.t_[ ivertex ] = tnew; } #ifdef VI_DEBUG else { - edm::LogInfo("sumw") << "invalid sum of weights in fit: " << vertices._sw[ivertex] << endl; + edm::LogInfo("sumw") << "invalid sum of weights in fit: " << vertices.sw_[ivertex] << endl; if (this->verbose_) { - std::cout << " a cluster melted away ? pk=" << vertices._pk[ ivertex ] << " sumw=" - << vertices._sw[ivertex] << endl; + std::cout << " a cluster melted away ? pk=" << vertices.pk_[ ivertex ] << " sumw=" + << vertices.sw_[ivertex] << endl; } } #endif @@ -287,7 +287,7 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, auto osumpi = 1./sumpi; for (unsigned int ivertex = 0; ivertex < nv; ++ivertex ) - vertices._pk[ ivertex ] = vertices._pk[ ivertex ] * vertices._se[ ivertex ] * osumpi; + vertices.pk_[ ivertex ] = vertices.pk_[ ivertex ] * vertices.se_[ ivertex ] * osumpi; return delta; }; @@ -314,10 +314,10 @@ bool DAClusterizerInZT_vect::merge(vertex_t & y, double & beta)const{ // merge the smallest distance clusters first std::vector > critical; for (unsigned int k = 0; (k + 1) < nv; k++) { - if ( std::fabs(y._z[k + 1] - y._z[k]) < zmerge_ && - std::fabs(y._t[k + 1] - y._t[k]) < tmerge_ ) { - auto dt2 = std::pow(y._t[k + 1] - y._t[k],2); - auto dz2 = std::pow(y._z[k + 1] - y._z[k],2); + if ( std::fabs(y.z_[k + 1] - y.z_[k]) < zmerge_ && + std::fabs(y.t_[k + 1] - y.t_[k]) < tmerge_ ) { + auto dt2 = std::pow(y.t_[k + 1] - y.t_[k],2); + auto dz2 = std::pow(y.z_[k + 1] - y.z_[k],2); critical.push_back( make_pair( dz2 + dt2, k ) ); } } @@ -328,24 +328,24 @@ bool DAClusterizerInZT_vect::merge(vertex_t & y, double & beta)const{ for (unsigned int ik=0; ik < critical.size(); ik++){ unsigned int k = critical[ik].second; - double rho = y._pk[k]+y._pk[k+1]; - double swE = y._swE[k]+y._swE[k+1]-y._pk[k]*y._pk[k+1] / rho*( std::pow(y._z[k+1]-y._z[k],2) + std::pow(y._t[k+1]-y._t[k],2) ); - double Tc = 2*swE / (y._sw[k]+y._sw[k+1]); + double rho = y.pk_[k]+y.pk_[k+1]; + double swE = y.swE_[k]+y.swE_[k+1]-y.pk_[k]*y.pk_[k+1] / rho*( std::pow(y.z_[k+1]-y.z_[k],2) + std::pow(y.t_[k+1]-y.t_[k],2) ); + double Tc = 2*swE / (y.sw_[k]+y.sw_[k+1]); if(Tc*beta < 1){ #ifdef VI_DEBUG - if(verbose_){ std::cout << "merging (" << y._z[k + 1] << ',' << y._t[k + 1] << ") and (" << y._z[k] << ',' << y._t[k] << ") Tc = " << Tc << " sw = " << y._sw[k]+y._sw[k+1] < (" << z1 << ',' << t1<< "),(" << z2 << ',' << t2 << ") [" << p1 << "," << p2 << "]" ; if (std::fabs(z2-z1) > epsilonz || std::fabs(t2-t1) > epsilont){ @@ -581,11 +581,11 @@ DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, do // split if the new subclusters are significantly separated if( std::fabs(z2-z1) > epsilonz || std::fabs(t2-t1) > epsilont){ split = true; - double pk1 = p1*y._pk[k]/(p1+p2); - double pk2 = p2*y._pk[k]/(p1+p2); - y._z[k] = z2; - y._t[k] = t2; - y._pk[k] = pk2; + double pk1 = p1*y.pk_[k]/(p1+p2); + double pk2 = p2*y.pk_[k]/(p1+p2); + y.z_[k] = z2; + y.t_[k] = t2; + y.pk_[k] = pk2; y.InsertItem(k, z1, t1, pk1); nv++; @@ -619,8 +619,8 @@ void DAClusterizerInZT_vect::splitAll( vertex_t & y) const { for (unsigned int k = 0; k < nv; k++) { if ( - ( ( (k == 0) || ( y._z[k - 1] < (y._z[k] - zsep)) ) && - ( ((k + 1) == nv) || ( y._z[k + 1] > (y._z[k] + zsep)) ) ) ) + ( ( (k == 0) || ( y.z_[k - 1] < (y.z_[k] - zsep)) ) && + ( ((k + 1) == nv) || ( y.z_[k + 1] > (y.z_[k] + zsep)) ) ) ) { // isolated prototype, split double new_z = y.z[k] - epsilonz; @@ -632,21 +632,21 @@ void DAClusterizerInZT_vect::splitAll( vertex_t & y) const { y.pk[k] = 0.5 * y.pk[k]; y1.AddItem(new_z, new_t, new_pk); - y1.AddItem(y._z[k], y._t[k], y._pk[k]); + y1.AddItem(y.z_[k], y.t_[k], y.pk_[k]); } else if ( (y1.GetSize() == 0 ) || - (y1._z[y1.GetSize() - 1] < (y._z[k] - zsep) ) || - (y1._t[y1.GetSize() - 1] < (y._t[k] - tsep) )) + (y1.z_[y1.GetSize() - 1] < (y.z_[k] - zsep) ) || + (y1.t_[y1.GetSize() - 1] < (y.t_[k] - tsep) )) { - y1.AddItem(y._z[k], y._t[k], y._pk[k]); + y1.AddItem(y.z_[k], y.t_[k], y.pk_[k]); } else { - y1._z[y1.GetSize() - 1] = y1._z[y1.GetSize() - 1] - epsilonz; - y1._t[y1.GetSize() - 1] = y1._t[y1.GetSize() - 1] - epsilont; - y._z[k] = y._z[k] + epsilonz; - y._t[k] = y._t[k] + epsilont; - y1.AddItem( y._z[k], y._t[k] , y._pk[k]); + y1.z_[y1.GetSize() - 1] = y1.z_[y1.GetSize() - 1] - epsilonz; + y1.t_[y1.GetSize() - 1] = y1.t_[y1.GetSize() - 1] - epsilont; + y.z_[k] = y.z_[k] + epsilonz; + y.t_[k] = y.t_[k] + epsilont; + y1.AddItem( y.z_[k], y.t_[k] , y.pk_[k]); } }// vertex loop @@ -819,30 +819,30 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co // ensure correct normalization of probabilities, should makes double assignment 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; y._z[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_); + tks.Z_sum_[i] = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); // improve vectorization (does not require reduction ....) for (unsigned int k = 0; k < nv; k++) { for (unsigned int i = 0; i < nt; i++) - tks._Z_sum[i] += y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],tks._dz2[i], tks._t[i], y._t[k],tks._dt2[i])); + tks.Z_sum_[i] += y.pk_[k] * local_exp(-beta * Eik(tks.z_[i], y.z_[k],tks.dz2_[i], tks.t_[i], y.t_[k],tks.dt2_[i])); } for (unsigned int k = 0; k < nv; k++) { - GlobalPoint pos(0, 0, y._z[k]); + GlobalPoint pos(0, 0, y.z_[k]); vector vertexTracks; for (unsigned int i = 0; i < nt; i++) { - if (tks._Z_sum[i] > 1e-100) { + if (tks.Z_sum_[i] > 1e-100) { - double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i], - tks._t[i], y._t[k], tks._dt2[i] )) / tks._Z_sum[i]; - if ((tks._pi[i] > 0) && (p > mintrkweight_)) { + double p = y.pk_[k] * local_exp(-beta * Eik(tks.z_[i], y.z_[k], tks.dz2_[i], + tks.t_[i], y.t_[k], tks.dt2_[i] )) / tks.Z_sum_[i]; + if ((tks.pi_[i] > 0) && (p > mintrkweight_)) { vertexTracks.push_back(*(tks.tt[i])); - tks._Z_sum[i] = 0; // setting Z=0 excludes double assignment + tks.Z_sum_[i] = 0; // setting Z=0 excludes double assignment } } } @@ -917,22 +917,22 @@ void DAClusterizerInZT_vect::dump(const double beta, const vertex_t & y, std::vector< unsigned int > iz; for(unsigned int j=0; jtrack().quality(reco::TrackBase::highPurity)) { std::cout << " *"; @@ -1006,19 +1006,19 @@ void DAClusterizerInZT_vect::dump(const double beta, const vertex_t & y, double sump = 0.; for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex) { - if (std::fabs(y._z[ivertex]-zdumpcenter_) > zdumpwidth_) continue; + if (std::fabs(y.z_[ivertex]-zdumpcenter_) > zdumpwidth_) continue; - if ((tks._pi[i] > 0) && (tks._Z_sum[i] > 0)) { + if ((tks.pi_[i] > 0) && (tks.Z_sum_[i] > 0)) { //double p=pik(beta,tks[i],*k); - double p = y._pk[ivertex] * exp(-beta * Eik(tks._z[i], y._z[ivertex], tks._dz2[i], - tks._t[i], y._t[ivertex], tks._dt2[i] )) / tks._Z_sum[i]; + double p = y.pk_[ivertex] * exp(-beta * Eik(tks.z_[i], y.z_[ivertex], tks.dz2_[i], + tks.t_[i], y.t_[ivertex], tks.dt2_[i] )) / tks.Z_sum_[i]; if (p > 0.0001) { std::cout << setw(8) << setprecision(3) << p; } else { std::cout << " . "; } - E += p * Eik(tks._z[i], y._z[ivertex], tks._dz2[i], - tks._t[i], y._t[ivertex], tks._dt2[i] ); + E += p * Eik(tks.z_[i], y.z_[ivertex], tks.dz2_[i], + tks.t_[i], y.t_[ivertex], tks.dt2_[i] ); sump += p; } else { std::cout << " "; From 5c7814586de7caf3addcbf1eaf1ad98bc9a4d0b5 Mon Sep 17 00:00:00 2001 From: lgray Date: Thu, 31 Aug 2017 13:42:52 -0500 Subject: [PATCH 08/16] code checks --- RecoVertex/PrimaryVertexProducer/BuildFile.xml | 2 +- .../interface/DAClusterizerInZT_vect.h | 2 +- .../PrimaryVertexProducer/interface/PrimaryVertexProducer.h | 6 +++--- .../PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc | 6 +++--- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/BuildFile.xml b/RecoVertex/PrimaryVertexProducer/BuildFile.xml index 0b584c37b5cd1..7f5f2666c0560 100644 --- a/RecoVertex/PrimaryVertexProducer/BuildFile.xml +++ b/RecoVertex/PrimaryVertexProducer/BuildFile.xml @@ -1,4 +1,4 @@ - + diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h index 2726f59faed56..e09227747551e 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -195,7 +195,7 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { std::vector > - clusterize(const std::vector & tracks) const; + clusterize(const std::vector & tracks) const override; std::vector diff --git a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h index 546d1f312bf8f..e58958cbc818c 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h @@ -55,10 +55,10 @@ class PrimaryVertexProducer : public edm::stream::EDProducer<> { public: - explicit PrimaryVertexProducer(const edm::ParameterSet&); - ~PrimaryVertexProducer(); + PrimaryVertexProducer(const edm::ParameterSet&); + ~PrimaryVertexProducer() override; - virtual void produce(edm::Event&, const edm::EventSetup&); + void produce(edm::Event&, const edm::EventSetup&) override; // access to config edm::ParameterSet config() const { return theConfig; } diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index 53268fe2a6bae..57d44ce236533 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -321,7 +321,7 @@ bool DAClusterizerInZT_vect::merge(vertex_t & y, double & beta)const{ critical.push_back( make_pair( dz2 + dt2, k ) ); } } - if (critical.size() == 0) return false; + if (critical.empty()) return false; std::stable_sort(critical.begin(), critical.end(), std::less >() ); @@ -508,7 +508,7 @@ DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, do critical.push_back( make_pair(Tc, k)); } } - if (critical.size()==0) return false; + if (critical.empty()) return false; std::stable_sort(critical.begin(), critical.end(), std::greater >() ); @@ -875,7 +875,7 @@ vector > DAClusterizerInZT_vect::clusterize( } #endif - if (pv.size() == 0) { + if (pv.empty()) { return clusters; } From 17ce31c82b305940fffdb88250a72563a62658af Mon Sep 17 00:00:00 2001 From: lgray Date: Thu, 31 Aug 2017 13:53:48 -0500 Subject: [PATCH 09/16] code review --- .../src/DAClusterizerInZT_vect.cc | 42 +++++++++---------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index 57d44ce236533..d67997ac9ae3c 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -31,9 +31,9 @@ DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) { zdumpwidth_ = conf.getUntrackedParameter ("zdumpwidth", 20.); // configurable parameters - double Tmin = conf.getParameter ("Tmin")*std::sqrt(2.0); - double Tpurge = conf.getParameter ("Tpurge")*std::sqrt(2.0); - double Tstop = conf.getParameter ("Tstop")*std::sqrt(2.0); + double minT = conf.getParameter ("Tmin")*std::sqrt(2.0); + double purgeT = conf.getParameter ("Tpurge")*std::sqrt(2.0); + double stopT = conf.getParameter ("Tstop")*std::sqrt(2.0); vertexSize_ = conf.getParameter ("vertexSize"); coolingFactor_ = conf.getParameter ("coolingFactor"); useTc_=true; @@ -54,9 +54,9 @@ DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) { std::cout << "DAClusterizerinZT_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl; std::cout << "DAClusterizerinZT_vect: zmerge = " << zmerge_ << std::endl; std::cout << "DAClusterizerinZT_vect: tmerge = " << tmerge_ << std::endl; - std::cout << "DAClusterizerinZT_vect: Tmin = " << Tmin << std::endl; - std::cout << "DAClusterizerinZT_vect: Tpurge = " << Tpurge << std::endl; - std::cout << "DAClusterizerinZT_vect: Tstop = " << Tstop << std::endl; + std::cout << "DAClusterizerinZT_vect: Tmin = " << minT << std::endl; + std::cout << "DAClusterizerinZT_vect: Tpurge = " << purgeT << std::endl; + std::cout << "DAClusterizerinZT_vect: Tstop = " << stopT << std::endl; std::cout << "DAClusterizerinZT_vect: vertexSize = " << vertexSize_ << std::endl; std::cout << "DAClusterizerinZT_vect: coolingFactor = " << coolingFactor_ << std::endl; std::cout << "DAClusterizerinZT_vect: d0CutOff = " << d0CutOff_ << std::endl; @@ -66,28 +66,28 @@ DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) { #endif - if (Tmin == 0) { - edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tmin" << Tmin + if (minT == 0) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tmin" << minT << " reset do default " << 1. / betamax_; } else { - betamax_ = 1. / Tmin; + betamax_ = 1. / minT; } - if ((Tpurge > Tmin) || (Tpurge == 0)) { - edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tpurge" << Tpurge - << " set to " << Tmin; - Tpurge = Tmin; + if ((purgeT > minT) || (purgeT == 0)) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tpurge" << purgeT + << " set to " << minT; + purgeT = minT; } - betapurge_ = 1./Tpurge; + betapurge_ = 1./purgeT; - if ((Tstop > Tpurge) || (Tstop == 0)) { - edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tstop" << Tstop - << " set to " << max(1., Tpurge); - Tstop = max(1., Tpurge) ; + if ((stopT > purgeT) || (stopT == 0)) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tstop" << stopT + << " set to " << max(1., purgeT); + stopT = max(1., purgeT) ; } - betastop_ = 1./Tstop; + betastop_ = 1./stopT; } @@ -692,7 +692,7 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co while ((update(beta, tks, y, false, rho0) > 1.e-6) && (niter++ < maxIterations_)) {} - // annealing loop, stop when T1/Tmin) + // annealing loop, stop when T1/minT) double betafreeze = betamax_ * sqrt(coolingFactor_); @@ -752,7 +752,7 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co #endif - // switch on outlier rejection at T=Tmin + // switch on outlier rejection at T=minT if(dzCutOff_ > 0){ rho0 = 1./nt; for(unsigned int a=0; a<10; a++){ update(beta, tks, y, true, a*rho0/10);} // adiabatic turn-on From 7f83a43423170cdb4d87089661f0f35153015905 Mon Sep 17 00:00:00 2001 From: lgray Date: Fri, 1 Sep 2017 11:41:28 -0500 Subject: [PATCH 10/16] further code review --- .../interface/DAClusterizerInZT_vect.h | 90 +++++++++---------- .../src/DAClusterizerInZT_vect.cc | 70 +++++++-------- 2 files changed, 76 insertions(+), 84 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h index e09227747551e..14688c49a3293 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -26,7 +26,7 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { // Internal data structure to struct track_t { - void AddItem( double new_z, double new_t, double new_dz2, double new_dt2, const reco::TransientTrack* new_tt, double new_pi ) + void addItem( double new_z, double new_t, double new_dz2, double new_dt2, const reco::TransientTrack* new_tt, double new_pi ) { z.push_back( new_z ); t.push_back( new_t ); @@ -38,17 +38,14 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { pi.push_back( new_pi ); // track weight Z_sum.push_back( 1.0 ); // Z[i] for DA clustering, initial value as done in ::fill } - - - unsigned int GetSize() const + unsigned int getSize() const { return z.size(); } - // has to be called everytime the items are modified - void ExtractRaw() + void extractRaw() { z_ = &z.front(); t_ = &t.front(); @@ -61,12 +58,12 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { double * z_; // z-coordinate at point of closest approach to the beamline double * t_; // t-coordinate at point of closest approach to the beamline + double * pi_; // track weight + double * dz2_; // square of the error of z(pca) double * dt2_; // square of the error of t(pca) - double * errsum_; // sum of squares of the pca errors - + double * errsum_; // sum of squares of the pca errors double * Z_sum_; // Z[i] for DA clustering - double * pi_; // track weight std::vector z; // z-coordinate at point of closest approach to the beamline std::vector t; // t-coordinate at point of closest approach to the beamline @@ -75,14 +72,13 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { std::vector errsum; // sum of squares of the pca errors std::vector Z_sum; // Z[i] for DA clustering std::vector pi; // track weight - std::vector< const reco::TransientTrack* > tt; // a pointer to the Transient Track }; struct vertex_t { std::vector z; // z coordinate std::vector t; // t coordinate - std::vector pk; // vertex weight for "constrained" clustering + std::vector pk; // vertex weight for "constrained" clustering // --- temporary numbers, used during update std::vector ei_cache; @@ -93,13 +89,7 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { std::vector se; std::vector swE; - - unsigned int GetSize() const - { - return z.size(); - } - - void AddItem( double new_z, double new_t, double new_pk ) + void addItem( double new_z, double new_t, double new_pk ) { z.push_back( new_z); t.push_back( new_t); @@ -113,10 +103,32 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { se.push_back( 0.0); swE.push_back( 0.0); - ExtractRaw(); + extractRaw(); } - void InsertItem( unsigned int i, double new_z, double new_t, double new_pk ) + unsigned int getSize() const + { + return z.size(); + } + + // has to be called everytime the items are modified + void extractRaw() + { + z_ = &z.front(); + t_ = &t.front(); + pk_ = &pk.front(); + + ei_ = &ei.front(); + sw_ = &sw.front(); + swz_ = &swz.front(); + swt_ = &swt.front(); + se_ = &se.front(); + swE_ = &swE.front(); + ei_cache_ = &ei_cache.front(); + + } + + void insertItem( unsigned int i, double new_z, double new_t, double new_pk ) { z.insert(z.begin() + i, new_z); t.insert(t.begin() + i, new_t); @@ -130,10 +142,10 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { se.insert( se.begin() + i, 0.0 ); swE.insert(swE.begin() + i, 0.0 ); - ExtractRaw(); + extractRaw(); } - void RemoveItem( unsigned int i ) + void removeItem( unsigned int i ) { z.erase( z.begin() + i ); t.erase( t.begin() + i ); @@ -147,36 +159,19 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { se.erase(se.begin() + i); swE.erase(swE.begin() + i); - ExtractRaw(); + extractRaw(); } - void DebugOut() + void debugOut() { - std::cout << "vertex_t size: " << GetSize() << std::endl; + std::cout << "vertex_t size: " << getSize() << std::endl; - for ( unsigned int i =0; i < GetSize(); ++ i) + for ( unsigned int i =0; i < getSize(); ++ i) { std::cout << " z = " << z_[i] << " t = " << t_[i] << " pk = " << pk_[i] << std::endl; } } - // has to be called everytime the items are modified - void ExtractRaw() - { - z_ = &z.front(); - t_ = &t.front(); - pk_ = &pk.front(); - - ei_ = &ei.front(); - sw_ = &sw.front(); - swz_ = &swz.front(); - swt_ = &swt.front(); - se_ = &se.front(); - swE_ = &swE.front(); - ei_cache_ = &ei_cache.front(); - - } - double * z_; double * t_; double * pk_; @@ -187,16 +182,13 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { double * swz_; double * swt_; double * se_; - double * swE_; - + double * swE_; }; - DAClusterizerInZT_vect(const edm::ParameterSet& conf); - + DAClusterizerInZT_vect(const edm::ParameterSet& conf); std::vector > - clusterize(const std::vector & tracks) const override; - + clusterize(const std::vector & tracks) const override; std::vector vertices(const std::vector & tracks, diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index d67997ac9ae3c..ddfa5afa053f6 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -135,13 +135,13 @@ DAClusterizerInZT_vect::fill(const vector & tracks) const if (edm::isNotFinite(t_pi) || t_pi < std::numeric_limits::epsilon()) continue; // usually is > 0.99 } LogTrace("DAClusterizerinZT_vectorized") << t_z << ' ' << t_t <<' '<< t_dz2 << ' ' << t_dt2 <<' '<< t_pi; - tks.AddItem(t_z, t_t, t_dz2, t_dt2, &(*it), t_pi); + tks.addItem(t_z, t_t, t_dz2, t_dt2, &(*it), t_pi); } - tks.ExtractRaw(); + tks.extractRaw(); #ifdef VI_DEBUG if (verbose_) { - std::cout << "Track count " << tks.GetSize() << std::endl; + std::cout << "Track count " << tks.getSize() << std::endl; } #endif @@ -163,8 +163,8 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, // mass constrained annealing without noise // returns the squared sum of changes of vertex positions - const unsigned int nt = gtracks.GetSize(); - const unsigned int nv = gvertices.GetSize(); + const unsigned int nt = gtracks.getSize(); + const unsigned int nv = gvertices.getSize(); //initialize sums double sumpi = 0.; @@ -306,7 +306,7 @@ bool DAClusterizerInZT_vect::merge(vertex_t & y, double & beta)const{ // merge clusters that collapsed or never separated, // only merge if the estimated critical temperature of the merged vertex is below the current temperature // return true if vertices were merged, false otherwise - const unsigned int nv = y.GetSize(); + const unsigned int nv = y.getSize(); if (nv < 2) return false; @@ -346,7 +346,7 @@ bool DAClusterizerInZT_vect::merge(vertex_t & y, double & beta)const{ y.pk_[k] = rho; y.sw_[k] += y.sw_[k+1]; y.swE_[k] = swE; - y.RemoveItem(k+1); + y.removeItem(k+1); return true; } } @@ -361,8 +361,8 @@ bool DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const double beta) const { constexpr double eps = 1.e-100; // eliminate clusters with only one significant/unique track - const unsigned int nv = y.GetSize(); - const unsigned int nt = tks.GetSize(); + const unsigned int nv = y.getSize(); + const unsigned int nt = tks.getSize(); if (nv < 2) return false; @@ -424,7 +424,7 @@ DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const << endl; } #endif - y.RemoveItem(k0); + y.removeItem(k0); return true; } else { return false; @@ -439,8 +439,8 @@ DAClusterizerInZT_vect::beta0(double betamax, track_t const & tks, vertex_t con double T0 = 0; // max Tc for beta=0 // estimate critical temperature from beta=0 (T=inf) - const unsigned int nt = tks.GetSize(); - const unsigned int nv = y.GetSize(); + const unsigned int nt = tks.getSize(); + const unsigned int nv = y.getSize(); for (unsigned int k = 0; k < nv; k++) { @@ -497,7 +497,7 @@ DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, do constexpr double epsilonz=1e-3; // minimum split size z constexpr double epsilont=1e-2; // minimum split size t - unsigned int nv = y.GetSize(); + unsigned int nv = y.getSize(); // avoid left-right biases by splitting highest Tc first @@ -515,7 +515,7 @@ DAClusterizerInZT_vect::split(const double beta, track_t &tks, vertex_t & y, do bool split=false; - const unsigned int nt = tks.GetSize(); + const unsigned int nt = tks.getSize(); for(unsigned int ic=0; ic DAClusterizerInZT_vect::vertices(const vector & tracks, const int verbosity) const { track_t && tks = fill(tracks); - tks.ExtractRaw(); + tks.extractRaw(); - unsigned int nt = tks.GetSize(); + unsigned int nt = tks.getSize(); double rho0 = 0.0; // start with no outlier rejection vector clusters; - if (tks.GetSize() == 0) return clusters; + if (tks.getSize() == 0) return clusters; vertex_t y; // the vertex prototypes // initialize:single vertex at infinite temperature - y.AddItem( 0, 0, 1.0); + y.addItem( 0, 0, 1.0); int niter = 0; // number of iterations @@ -787,7 +787,7 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co // eliminate insigificant vertices, this is more restrictive at higher T while (purge(y, tks, rho0, beta)) { niter = 0; - while (( update(beta, tks, y, true, rho0) > 2.5e-7 * y.GetSize() ) && (niter++ < maxIterations_)) { + while (( update(beta, tks, y, true, rho0) > 2.5e-7 * y.getSize() ) && (niter++ < maxIterations_)) { } } @@ -817,7 +817,7 @@ DAClusterizerInZT_vect::vertices(const vector & tracks, co GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01); // ensure correct normalization of probabilities, should makes double assignment reasonably impossible - const unsigned int nv = y.GetSize(); + 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; y.z_[k]=0;} @@ -912,8 +912,8 @@ vector > DAClusterizerInZT_vect::clusterize( void DAClusterizerInZT_vect::dump(const double beta, const vertex_t & y, const track_t & tks, int verbosity) const { - const unsigned int nv = y.GetSize(); - const unsigned int nt = tks.GetSize(); + const unsigned int nv = y.getSize(); + const unsigned int nt = tks.getSize(); std::vector< unsigned int > iz; for(unsigned int j=0; j Date: Fri, 1 Sep 2017 14:15:06 -0500 Subject: [PATCH 11/16] more... --- .../interface/DAClusterizerInZT_vect.h | 29 ++++++++------- .../python/TkClusParameters_cff.py | 2 ++ .../src/DAClusterizerInZT_vect.cc | 36 +++++++++---------- 3 files changed, 35 insertions(+), 32 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h index 14688c49a3293..363e9acdd3686 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -76,18 +76,6 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { }; struct vertex_t { - std::vector z; // z coordinate - std::vector t; // t coordinate - std::vector pk; // vertex weight for "constrained" clustering - - // --- temporary numbers, used during update - std::vector ei_cache; - std::vector ei; - std::vector sw; - std::vector swz; - std::vector swt; - std::vector se; - std::vector swE; void addItem( double new_z, double new_t, double new_pk ) { @@ -172,6 +160,10 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { } } + std::vector z; // z coordinate + std::vector t; // t coordinate + std::vector pk; // vertex weight for "constrained" clustering + double * z_; double * t_; double * pk_; @@ -182,7 +174,16 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { double * swz_; double * swt_; double * se_; - double * swE_; + double * swE_; + + // --- temporary numbers, used during update + std::vector ei_cache; + std::vector ei; + std::vector sw; + std::vector swz; + std::vector swt; + std::vector se; + std::vector swE; }; DAClusterizerInZT_vect(const edm::ParameterSet& conf); @@ -216,12 +217,14 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { double zdumpwidth_; double vertexSize_; + double vertexSizeTime_; int maxIterations_; double coolingFactor_; double betamax_; double betastop_; double dzCutOff_; double d0CutOff_; + double dtCutOff_; bool useTc_; double mintrkweight_; diff --git a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py index 0f6775ff8006a..3f20a56a2b3c0 100644 --- a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py +++ b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py @@ -8,8 +8,10 @@ Tpurge = cms.double(2.0), # cleaning Tstop = cms.double(0.5), # end of annealing vertexSize = cms.double(0.006), # added in quadrature to track-z resolutions + vertexSizeTime = cms.double(0.008), d0CutOff = cms.double(3.), # downweight high IP tracks dzCutOff = cms.double(3.), # outlier rejection after freeze-out (T ("Tpurge")*std::sqrt(2.0); double stopT = conf.getParameter ("Tstop")*std::sqrt(2.0); vertexSize_ = conf.getParameter ("vertexSize"); + vertexSizeTime_ = conf.getParameter ("vertexSizeTime"); coolingFactor_ = conf.getParameter ("coolingFactor"); useTc_=true; if(coolingFactor_<0){ @@ -44,6 +39,7 @@ DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) { coolingFactor_ = std::sqrt(coolingFactor_); d0CutOff_ = conf.getParameter ("d0CutOff"); dzCutOff_ = conf.getParameter ("dzCutOff"); + dtCutOff_ = conf.getParameter ("dtCutOff"); uniquetrkweight_ = conf.getParameter("uniquetrkweight"); zmerge_ = conf.getParameter("zmerge"); tmerge_ = conf.getParameter("tmerge"); @@ -58,10 +54,11 @@ DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) { std::cout << "DAClusterizerinZT_vect: Tpurge = " << purgeT << std::endl; std::cout << "DAClusterizerinZT_vect: Tstop = " << stopT << std::endl; std::cout << "DAClusterizerinZT_vect: vertexSize = " << vertexSize_ << std::endl; + std::cout << "DAClusterizerinZT_vect: vertexSizeTime = " << vertexSizeTime_ << std::endl; std::cout << "DAClusterizerinZT_vect: coolingFactor = " << coolingFactor_ << std::endl; std::cout << "DAClusterizerinZT_vect: d0CutOff = " << d0CutOff_ << std::endl; std::cout << "DAClusterizerinZT_vect: dzCutOff = " << dzCutOff_ << std::endl; - std::cout << "DAClusterizerinZT_vect: dtCutoff = " << dtCutOff << std::endl; + std::cout << "DAClusterizerinZT_vect: dtCutoff = " << dtCutOff_ << std::endl; } #endif @@ -111,31 +108,31 @@ DAClusterizerInZT_vect::fill(const vector & tracks) const // prepare track data for clustering track_t tks; - for (auto it = tracks.begin(); it!= tracks.end(); it++){ - if (!(*it).isValid()) continue; + for( const auto& tk : tracks ) { + if (!tk.isValid()) continue; double t_pi=1.; - double t_z = ((*it).stateAtBeamLine().trackStateAtPCA()).position().z(); - double t_t = it->timeExt(); + double t_z = tk.stateAtBeamLine().trackStateAtPCA().position().z(); + double t_t = tk.timeExt(); if (std::fabs(t_z) > 1000.) continue; - auto const & t_mom = (*it).stateAtBeamLine().trackStateAtPCA().momentum(); + auto const & t_mom = tk.stateAtBeamLine().trackStateAtPCA().momentum(); // get the beam-spot - reco::BeamSpot beamspot = (it->stateAtBeamLine()).beamSpot(); + reco::BeamSpot beamspot = tk.stateAtBeamLine().beamSpot(); double t_dz2 = - std::pow((*it).track().dzError(), 2) // track errror + std::pow(tk.track().dzError(), 2) // track errror + (std::pow(beamspot.BeamWidthX()*t_mom.x(),2)+std::pow(beamspot.BeamWidthY()*t_mom.y(),2))*std::pow(t_mom.z(),2)/std::pow(t_mom.perp2(),2) // beam spot width + std::pow(vertexSize_, 2); // intrinsic vertex size, safer for outliers and short lived decays t_dz2 = 1./ t_dz2; - double t_dt2 =std::pow((*it).dtErrorExt(),2.) + std::pow(vertexSizeTime,2.); // the ~injected~ timing error, need to add a small minimum vertex size in time + double t_dt2 =std::pow(tk.dtErrorExt(),2.) + std::pow(vertexSizeTime_,2.); // the ~injected~ timing error, need to add a small minimum vertex size in time t_dt2 = 1./t_dt2; if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits::min() ) continue; if (d0CutOff_ > 0) { Measurement1D atIP = - (*it).stateAtBeamLine().transverseImpactParameter();// error contains beamspot + tk.stateAtBeamLine().transverseImpactParameter();// error contains beamspot t_pi = 1. / (1. + local_exp(std::pow(atIP.value() / atIP.error(), 2) - std::pow(d0CutOff_, 2))); // reduce weight for high ip tracks if (edm::isNotFinite(t_pi) || t_pi < std::numeric_limits::epsilon()) continue; // usually is > 0.99 } LogTrace("DAClusterizerinZT_vectorized") << t_z << ' ' << t_t <<' '<< t_dz2 << ' ' << t_dt2 <<' '<< t_pi; - tks.addItem(t_z, t_t, t_dz2, t_dt2, &(*it), t_pi); + tks.addItem(t_z, t_t, t_dz2, t_dt2, &it, t_pi); } tks.extractRaw(); @@ -910,8 +907,8 @@ vector > DAClusterizerInZT_vect::clusterize( void DAClusterizerInZT_vect::dump(const double beta, const vertex_t & y, - const track_t & tks, int verbosity) const { - + const track_t & tks, int verbosity) const { +#ifdef VI_DEBUG const unsigned int nv = y.getSize(); const unsigned int nt = tks.getSize(); @@ -1029,4 +1026,5 @@ void DAClusterizerInZT_vect::dump(const double beta, const vertex_t & y, std::cout << endl << "T=" << 1 / beta << " E=" << E << " n=" << y.getSize() << " F= " << F << endl << "----------" << endl; } +#endif } From 0208b2daa0ab492661ae591c59432fa0aaa976b8 Mon Sep 17 00:00:00 2001 From: lgray Date: Fri, 1 Sep 2017 14:54:04 -0500 Subject: [PATCH 12/16] it to tk --- RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index c6a590a51d448..96d0dd359b9a9 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -132,7 +132,7 @@ DAClusterizerInZT_vect::fill(const vector & tracks) const if (edm::isNotFinite(t_pi) || t_pi < std::numeric_limits::epsilon()) continue; // usually is > 0.99 } LogTrace("DAClusterizerinZT_vectorized") << t_z << ' ' << t_t <<' '<< t_dz2 << ' ' << t_dt2 <<' '<< t_pi; - tks.addItem(t_z, t_t, t_dz2, t_dt2, &it, t_pi); + tks.addItem(t_z, t_t, t_dz2, t_dt2, &tk, t_pi); } tks.extractRaw(); From 74d36d802a3a98c45b9114591d16ccadd8f77e93 Mon Sep 17 00:00:00 2001 From: lgray Date: Mon, 4 Sep 2017 15:19:47 -0500 Subject: [PATCH 13/16] put time related parameters in the right PSet --- .../PrimaryVertexProducer/python/TkClusParameters_cff.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py index 3f20a56a2b3c0..1625308e1b260 100644 --- a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py +++ b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py @@ -7,11 +7,9 @@ Tmin = cms.double(2.0), # end of vertex splitting Tpurge = cms.double(2.0), # cleaning Tstop = cms.double(0.5), # end of annealing - vertexSize = cms.double(0.006), # added in quadrature to track-z resolutions - vertexSizeTime = cms.double(0.008), + vertexSize = cms.double(0.006), # added in quadrature to track-z resolutions d0CutOff = cms.double(3.), # downweight high IP tracks - dzCutOff = cms.double(3.), # outlier rejection after freeze-out (T Date: Wed, 6 Sep 2017 14:41:16 -0500 Subject: [PATCH 14/16] some code review and changes to validation for timing --- .../src/DAClusterizerInZT_vect.cc | 24 +++++++++---------- .../RecoTrack/python/TrackValidation_cff.py | 18 ++++++++++++++ 2 files changed, 29 insertions(+), 13 deletions(-) diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index 96d0dd359b9a9..93dfb189602cd 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -125,6 +125,7 @@ DAClusterizerInZT_vect::fill(const vector & tracks) const double t_dt2 =std::pow(tk.dtErrorExt(),2.) + std::pow(vertexSizeTime_,2.); // the ~injected~ timing error, need to add a small minimum vertex size in time t_dt2 = 1./t_dt2; if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits::min() ) continue; + if (edm::isNotFinite(t_dt2) || t_dt2 < std::numeric_limits::min() ) continue; if (d0CutOff_ > 0) { Measurement1D atIP = tk.stateAtBeamLine().transverseImpactParameter();// error contains beamspot @@ -210,8 +211,6 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, vertex_t & y_vec ) { auto tmp_trk_pi = tks_vec.pi_[track_num]; auto o_trk_Z_sum = 1./tks_vec.Z_sum_[track_num]; - //auto o_trk_dz2 = tks_vec.dz2_[track_num]; - //auto o_trk_dt2 = tks_vec.dt2_[track_num]; auto o_trk_err_sum = tks_vec.errsum_[track_num]; auto tmp_trk_z = tks_vec.z_[track_num]; auto tmp_trk_t = tks_vec.t_[track_num]; @@ -254,14 +253,13 @@ double DAClusterizerInZT_vect::update(double beta, track_t & gtracks, } } - // now update z and pk + // now update z, t, and pk auto kernel_calc_zt = [ sumpi, nv, this, useRho0 ] (vertex_t & vertices ) -> double { double delta=0; // does not vectorizes for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex ) { if (vertices.sw_[ivertex] > 0.) { - //std::cout << " DA2D_vect sw = " << vertices.sw_[ ivertex ] << ' ' << vertices.swz_[ ivertex ] << ' ' << vertices.swt_[ ivertex ] << std::endl; auto znew = vertices.swz_[ ivertex ] / vertices.sw_[ ivertex ]; // prevents from vectorizing if delta += std::pow( vertices.z_[ ivertex ] - znew, 2 ); @@ -371,12 +369,12 @@ DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const double sump = 0; std::vector inverse_zsums(nt), arg_cache(nt), eik_cache(nt); - double * inverse_zsums_; - double * arg_cache_; - double * eik_cache_; - inverse_zsums_ = inverse_zsums.data(); - arg_cache_ = arg_cache.data(); - eik_cache_ = eik_cache.data(); + double * pinverse_zsums; + double * parg_cache; + double * peik_cache; + pinverse_zsums = inverse_zsums.data(); + parg_cache = arg_cache.data(); + peik_cache = eik_cache.data(); for(unsigned i = 0; i < nt; ++i) { inverse_zsums[i] = tks.Z_sum_[i] > eps ? 1./tks.Z_sum_[i] : 0.0; } @@ -396,11 +394,11 @@ DAClusterizerInZT_vect::purge(vertex_t & y, track_t & tks, double & rho0, const const auto mult_resz = track_z - y.z_[k]; const auto mult_rest = track_t - y.t_[k]; - arg_cache_[i] = botrack_dz2 * ( mult_resz * mult_resz ) + botrack_dt2 * ( mult_rest * mult_rest ); + parg_cache[i] = botrack_dz2 * ( mult_resz * mult_resz ) + botrack_dt2 * ( mult_rest * mult_rest ); } - local_exp_list(arg_cache_, eik_cache_, nt); + local_exp_list(parg_cache, peik_cache, nt); for (unsigned int i = 0; i < nt; ++i) { - const double p = y.pk_[k] * eik_cache_[i] * inverse_zsums_[i]; + const double p = y.pk_[k] * peik_cache[i] * pinverse_zsums[i]; sump += p; nUnique += ( ( p > pcut ) & ( tks.pi_[i] > 0 ) ); } diff --git a/Validation/RecoTrack/python/TrackValidation_cff.py b/Validation/RecoTrack/python/TrackValidation_cff.py index 5c98bdf9bae9b..6cb28ee566760 100644 --- a/Validation/RecoTrack/python/TrackValidation_cff.py +++ b/Validation/RecoTrack/python/TrackValidation_cff.py @@ -674,3 +674,21 @@ def _getMVASelectors(postfix): tracksValidationTruth + trackValidatorLite ) + +## customization for timing +from Configuration.Eras.Modifier_phase2_timing_cff import phase2_timing +phase2_timing.toModify( generalTracksFromPV, + vertexTag = cms.InputTag('offlinePrimaryVertices4D'), + timesTag = cms.InputTag('trackTimeValueMapProducer:generalTracksConfigurableFlatResolutionModel'), + timeResosTag = cms.InputTag('trackTimeValueMapProducer:generalTracksConfigurableFlatResolutionModelResolution'), + nSigmaDtVertex = cms.double(3) ) +phase2_timing.toModify( trackValidatorStandalone, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +phase2_timing.toModify( trackValidatorFromPVStandalone, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +phase2_timing.toModify( trackValidatorFromPVAllTPStandalone, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +phase2_timing.toModify( trackValidatorConversionStandalone, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +phase2_timing.toModify( trackValidatorGsfTracks, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) From efbe344b8290dadec925c0c5c3e3f9329ecf5e9c Mon Sep 17 00:00:00 2001 From: lgray Date: Thu, 7 Sep 2017 11:53:04 -0500 Subject: [PATCH 15/16] move to timing_layer era --- .../RecoTrack/python/TrackValidation_cff.py | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/Validation/RecoTrack/python/TrackValidation_cff.py b/Validation/RecoTrack/python/TrackValidation_cff.py index 6cb28ee566760..4606f78bd360e 100644 --- a/Validation/RecoTrack/python/TrackValidation_cff.py +++ b/Validation/RecoTrack/python/TrackValidation_cff.py @@ -676,19 +676,19 @@ def _getMVASelectors(postfix): ) ## customization for timing -from Configuration.Eras.Modifier_phase2_timing_cff import phase2_timing -phase2_timing.toModify( generalTracksFromPV, - vertexTag = cms.InputTag('offlinePrimaryVertices4D'), - timesTag = cms.InputTag('trackTimeValueMapProducer:generalTracksConfigurableFlatResolutionModel'), - timeResosTag = cms.InputTag('trackTimeValueMapProducer:generalTracksConfigurableFlatResolutionModelResolution'), - nSigmaDtVertex = cms.double(3) ) -phase2_timing.toModify( trackValidatorStandalone, - label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) -phase2_timing.toModify( trackValidatorFromPVStandalone, - label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) -phase2_timing.toModify( trackValidatorFromPVAllTPStandalone, - label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) -phase2_timing.toModify( trackValidatorConversionStandalone, - label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) -phase2_timing.toModify( trackValidatorGsfTracks, - label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +from Configuration.Eras.Modifier_phase2_timing_layer_cff import phase2_timing_layer +phase2_timing_layer.toModify( generalTracksFromPV, + vertexTag = cms.InputTag('offlinePrimaryVertices4D'), + timesTag = cms.InputTag('trackTimeValueMapProducer:generalTracksConfigurableFlatResolutionModel'), + timeResosTag = cms.InputTag('trackTimeValueMapProducer:generalTracksConfigurableFlatResolutionModelResolution'), + nSigmaDtVertex = cms.double(3) ) +phase2_timing_layer.toModify( trackValidatorStandalone, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +phase2_timing_layer.toModify( trackValidatorFromPVStandalone, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +phase2_timing_layer.toModify( trackValidatorFromPVAllTPStandalone, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +phase2_timing_layer.toModify( trackValidatorConversionStandalone, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) +phase2_timing_layer.toModify( trackValidatorGsfTracks, + label_vertex = cms.untracked.InputTag('offlinePrimaryVertices4D') ) From 163dda924ab37b565c2b2ff011882cbf7b93c968 Mon Sep 17 00:00:00 2001 From: lgray Date: Mon, 11 Sep 2017 12:03:20 -0500 Subject: [PATCH 16/16] 4D vertex validation --- .../PrimaryVertexAnalyzer4PUSlimmed_cfi.py | 23 +++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/Validation/RecoVertex/python/PrimaryVertexAnalyzer4PUSlimmed_cfi.py b/Validation/RecoVertex/python/PrimaryVertexAnalyzer4PUSlimmed_cfi.py index b1235890467f1..d75550870c760 100644 --- a/Validation/RecoVertex/python/PrimaryVertexAnalyzer4PUSlimmed_cfi.py +++ b/Validation/RecoVertex/python/PrimaryVertexAnalyzer4PUSlimmed_cfi.py @@ -24,8 +24,8 @@ vertexRecoCollections = cms.VInputTag("offlinePrimaryVertices", "offlinePrimaryVerticesWithBS", "selectedOfflinePrimaryVertices", - "selectedOfflinePrimaryVerticesWithBS", - ), + "selectedOfflinePrimaryVerticesWithBS" + ), ) vertexAnalysisTrackingOnly = vertexAnalysis.clone( @@ -83,3 +83,22 @@ + pixelVertexAnalysisTrackingOnly ) trackingLowPU.toReplaceWith(vertexAnalysisSequenceTrackingOnly, _vertexAnalysisSequenceTrackingOnly_trackingLowPU) + +from Configuration.Eras.Modifier_phase2_timing_layer_cff import phase2_timing_layer +_vertexRecoCollectionsTiming = cms.VInputTag("offlinePrimaryVertices", + "offlinePrimaryVerticesWithBS", + "selectedOfflinePrimaryVertices", + "selectedOfflinePrimaryVerticesWithBS", + "offlinePrimaryVertices4D", + "selectedOfflinePrimaryVertices4D", + ) +selectedOfflinePrimaryVertices4D = selectedOfflinePrimaryVertices.clone(src = cms.InputTag("offlinePrimaryVertices4D")) + +_vertexAnalysisSelectionTiming = vertexAnalysisSelection.copy() +_vertexAnalysisSelectionTiming += selectedOfflinePrimaryVertices4D + +phase2_timing_layer.toModify( vertexAnalysis, + vertexRecoCollections = _vertexRecoCollectionsTiming + ) +phase2_timing_layer.toReplaceWith( vertexAnalysisSelection, + _vertexAnalysisSelectionTiming )