DAClusterizerInZ_vect.h

Go to the documentation of this file.

#ifndef DAClusterizerInZ_vect_h
#define DAClusterizerInZ_vect_h

#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 DAClusterizerInZ_vect  final : public TrackClusterizerInZ {

public:
  // Internal data structure to 
  struct track_t {
    
    void AddItem( double new_z, double new_dz2, const reco::TransientTrack* new_tt, double new_pi   )
    {
      z.push_back( new_z );
      dz2.push_back( new_dz2 );
      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();
      _dz2 = &dz2.front();
      _Z_sum = &Z_sum.front();
      _pi = &pi.front();
    }
    
    double * __restrict__ _z; // z-coordinate at point of closest approach to the beamline
    double * __restrict__  _dz2; // square of the error of z(pca)
    
    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> dz2; // square of the error of z(pca)
    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> 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> se;
    std::vector<double> swE;
    
    
    unsigned int GetSize() const
    {
      return z.size();
    }
    
    void AddItem( double new_z, double new_pk   )
    {
      z.push_back( new_z);
      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);
      se.push_back( 0.0);
      swE.push_back( 0.0);
      
      ExtractRaw();
    }
    
    void InsertItem( unsigned int i, double new_z, double new_pk   )
    {
      z.insert(z.begin() + i, new_z);
      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 );
      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 );
      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);
      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] << " pk = " << _pk[i] << std::endl;
    }
    }
    
    // has to be called everytime the items are modified
    void ExtractRaw()
    {
      _z = &z.front();
      _pk = &pk.front();
      
      _ei = &ei.front();
      _sw = &sw.front();
      _swz = &swz.front();
      _se = &se.front();
      _swE = &swE.front();
      _ei_cache = &ei_cache.front();
      
    }
    
    double * __restrict__ _z;
    double * __restrict__ _pk;
    
    double * __restrict__ _ei_cache;
    double * __restrict__ _ei;
    double * __restrict__ _sw;
    double * __restrict__ _swz;
    double * __restrict__ _se;
    double * __restrict__ _swE;
    
  };
  
  DAClusterizerInZ_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 betapurge_;

};


//#ifndef DAClusterizerInZ_new_h
#endif