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.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 new file mode 100644 index 0000000000000..363e9acdd3686 --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -0,0 +1,240 @@ +#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" + +#include + +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 * 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 * Z_sum_; // Z[i] for DA clustering + + 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 { + + 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(); + } + + 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); + 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; + } + } + + std::vector z; // z coordinate + std::vector t; // t coordinate + std::vector pk; // vertex weight for "constrained" clustering + + double * z_; + double * t_; + double * pk_; + + double * ei_cache_; + double * ei_; + double * sw_; + double * swz_; + double * swt_; + double * se_; + 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); + + std::vector > + clusterize(const std::vector & tracks) const override; + + 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_; + double vertexSizeTime_; + int maxIterations_; + double coolingFactor_; + double betamax_; + double betastop_; + double dzCutOff_; + double d0CutOff_; + double dtCutOff_; + 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..e58958cbc818c 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h @@ -34,12 +34,13 @@ #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" #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" @@ -54,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/plugins/PrimaryVertexProducer.cc b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc index 3aef3af31403d..6599da9825bfe 100644 --- a/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc +++ b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc @@ -49,9 +49,8 @@ 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")); + } else if( clusteringAlgorithm=="DA2D_vect" ) { + theTrackClusterizer = new DAClusterizerInZT_vect(conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); f4D = true; } diff --git a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py index 10afac03823ca..1625308e1b260 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.01), # ~ 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="< +#include +#include +#include +#include "FWCore/Utilities/interface/isFinite.h" +#include "vdt/vdtMath.h" + +using namespace std; + +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 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"); + vertexSizeTime_ = conf.getParameter ("vertexSizeTime"); + coolingFactor_ = conf.getParameter ("coolingFactor"); + useTc_=true; + if(coolingFactor_<0){ + coolingFactor_=-coolingFactor_; + useTc_=false; + } + 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"); + +#ifdef VI_DEBUG + 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 = " << 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: 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; + } +#endif + + + if (minT == 0) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tmin" << minT + << " reset do default " << 1. / betamax_; + } else { + betamax_ = 1. / minT; + } + + + if ((purgeT > minT) || (purgeT == 0)) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tpurge" << purgeT + << " set to " << minT; + purgeT = minT; + } + betapurge_ = 1./purgeT; + + + if ((stopT > purgeT) || (stopT == 0)) { + edm::LogWarning("DAClusterizerinZT_vectorized") << "DAClusterizerInZT: invalid Tstop" << stopT + << " set to " << max(1., purgeT); + stopT = max(1., purgeT) ; + } + betastop_ = 1./stopT; + +} + + +namespace { + inline double local_exp( double const& inp) { + return vdt::fast_exp( inp ); + } + + 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]); + } + +} + +//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( const auto& tk : tracks ) { + if (!tk.isValid()) continue; + double t_pi=1.; + double t_z = tk.stateAtBeamLine().trackStateAtPCA().position().z(); + double t_t = tk.timeExt(); + if (std::fabs(t_z) > 1000.) continue; + auto const & t_mom = tk.stateAtBeamLine().trackStateAtPCA().momentum(); + // get the beam-spot + reco::BeamSpot beamspot = tk.stateAtBeamLine().beamSpot(); + double t_dz2 = + 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(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 + 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, &tk, t_pi); + } + tks.extractRaw(); + +#ifdef VI_DEBUG + if (verbose_) { + std::cout << "Track count " << tks.getSize() << std::endl; + } +#endif + + 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 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 ); + } + }; + + 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_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; + } + }; + + + 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, 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.) { + 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.empty()) 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){ +#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; + }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 { + 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; + + std::vector inverse_zsums(nt), arg_cache(nt), eik_cache(nt); + 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; + } + + for (unsigned int k = 0; k < nv; ++k) { + + 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]; + parg_cache[i] = botrack_dz2 * ( mult_resz * mult_resz ) + botrack_dt2 * ( mult_rest * mult_rest ); + } + local_exp_list(parg_cache, peik_cache, nt); + for (unsigned int i = 0; i < nt; ++i) { + const double p = y.pk_[k] * peik_cache[i] * pinverse_zsums[i]; + sump += p; + nUnique += ( ( p > pcut ) & ( tks.pi_[i] > 0 ) ); + } + + if ((nUnique < 2) && (sump < sumpmin)) { + sumpmin = sump; + k0 = k; + } + + } + + 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 { + 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) + +#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) { + 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-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)); + } + } + if (critical.empty()) 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] ) * 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; + 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.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]; + } + +#ifdef VI_DEBUG + 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; + } + } + } +#endif + + // 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; + +#ifdef VI_DEBUG + if (verbose_) { + std::cout << "Before Split "<< std::endl; + y.DebugOut(); + } +#endif + + 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)) ) ) ) + { + // 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(); + +#ifdef VI_DEBUG + if (verbose_) { + std::cout << "After split " << std::endl; + y.DebugOut(); + } +#endif +} + + + +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); +#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) && + (niter++ < maxIterations_)) {} + + // annealing loop, stop when T1/minT) + + 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 +#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; + 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);} +#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; + } + }else{ + // merge collapsed clusters + 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=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 + } + + 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_ ){ + 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) > 2.5e-7 * y.getSize() ) && (niter++ < maxIterations_)) { + } + } + +#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_ ){ + beta = min( beta/coolingFactor_, betastop_); + niter =0; + 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); + + // 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 { + +#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.empty()) { + 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{ +#ifdef VI_DEBUG + if(verbose_){ + std::cout << " one track cluster at " << k->position().z() << " suppressed" << std::endl; + } +#endif + } + 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 { +#ifdef VI_DEBUG + 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; + } +#endif +} diff --git a/Validation/RecoTrack/python/TrackValidation_cff.py b/Validation/RecoTrack/python/TrackValidation_cff.py index 5c98bdf9bae9b..4606f78bd360e 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_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') ) 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 )