vectorized 2D simulated annealing vertexing

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)
 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")

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 <vector>
+#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<double> z; // z-coordinate at point of closest approach to the beamline
+    std::vector<double> t; // t-coordinate at point of closest approach to the beamline
+    std::vector<double> dz2; // square of the error of z(pca)
+    std::vector<double> dt2; // square of the error of t(pca)
+    std::vector<double> errsum; // sum of squares of the pca errors
+    std::vector< const reco::TransientTrack* > tt; // a pointer to the Transient Track
+
+    std::vector<double> Z_sum; // Z[i]   for DA clustering
+    std::vector<double> pi; // track weight
+  };
+
+  struct vertex_t {
+    std::vector<double> z; //           z coordinate
+    std::vector<double> t; //           t coordinate
+    std::vector<double> pk; //           vertex weight for "constrained" clustering
+
+    // --- temporary numbers, used during update
+    std::vector<double> ei_cache;
+    std::vector<double> ei;
+    std::vector<double> sw;
+    std::vector<double> swz;
+    std::vector<double> swt;
+    std::vector<double> se;
+    std::vector<double> 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<std::vector<reco::TransientTrack> >
+  clusterize(const std::vector<reco::TransientTrack> & tracks) const;
+
+
+  std::vector<TransientVertex>
+  vertices(const std::vector<reco::TransientTrack> & tracks,
+	   const int verbosity = 0) const ;
+
+  track_t	fill(const std::vector<reco::TransientTrack> & 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

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"

RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc

@@ -53,6 +53,9 @@ PrimaryVertexProducer::PrimaryVertexProducer(const edm::ParameterSet& conf)
   else if( clusteringAlgorithm=="DA2D" ) {
     theTrackClusterizer = new DAClusterizerInZT(conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
     f4D = true;
+  } else if( clusteringAlgorithm=="DA2D_vect" ) {
+    theTrackClusterizer = new DAClusterizerInZT_vect(conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
+    f4D = true;
   }
 
   else{

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 (T<Tmin)
+        zmerge = cms.double(1e-2),        # merge intermediat clusters separated by less than zmerge and tmerge
+        tmerge = cms.double(1e-2),        # merge intermediat clusters separated by less than zmerge and tmerge
+        uniquetrkweight = cms.double(0.8) # require at least two tracks with this weight at T=Tpurge
+        )
+)

RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT.cc

@@ -324,7 +324,7 @@ bool DAClusterizerInZT::split( double beta,
       }
     }
     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;}
+    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); }