Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-11 22:44:12

0001 #ifndef RecoVertex_PrimaryVertexProducer_DAClusterizerInZ_vect_h
0002 #define RecoVertex_PrimaryVertexProducer_DAClusterizerInZ_vect_h
0003 
0004 /**\class DAClusterizerInZ_vect
0005 
0006  Description: separates event tracks into clusters along the beam line
0007 
0008     Version which auto-vectorizes with gcc 4.6 or newer
0009 
0010  */
0011 
0012 #include "RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h"
0013 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include <vector>
0017 #include "DataFormats/Math/interface/Error.h"
0018 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0019 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0020 
0021 class DAClusterizerInZ_vect final : public TrackClusterizerInZ {
0022 public:
0023   static void fillPSetDescription(edm::ParameterSetDescription &desc);
0024 
0025   // internal data structure for tracks
0026   struct track_t {
0027     std::vector<double> zpca_vec;                  // z-coordinate at point of closest approach to the beamline
0028     std::vector<double> dz2_vec;                   // square of the error of z(pca)
0029     std::vector<double> sum_Z_vec;                 // track contribution to the partition function, Z
0030     std::vector<double> tkwt_vec;                  // track weight, close to 1.0 for most tracks
0031     std::vector<unsigned int> kmin;                // index of the first cluster within zrange
0032     std::vector<unsigned int> kmax;                // 1 + index of the last cluster within zrange
0033     std::vector<const reco::TransientTrack *> tt;  // a pointer to the Transient Track
0034 
0035     double osumtkwt;  // 1. / (sum of all track weights)
0036 
0037     void addItemSorted(double new_zpca, double new_dz2, const reco::TransientTrack *new_tt, double new_tkwt) {
0038       // sort tracks with decreasing resolution (note that dz2 = 1/sigma^2)
0039       unsigned int i = 0;
0040       for (i = 0; i < zpca_vec.size(); i++) {
0041         if (new_dz2 > dz2_vec[i])
0042           break;
0043       }
0044       insertItem(i, new_zpca, new_dz2, new_tt, new_tkwt);
0045     }
0046 
0047     void insertItem(
0048         unsigned int i, double new_zpca, double new_dz2, const reco::TransientTrack *new_tt, double new_tkwt) {
0049       zpca_vec.insert(zpca_vec.begin() + i, new_zpca);
0050       dz2_vec.insert(dz2_vec.begin() + i, new_dz2);
0051       tt.insert(tt.begin() + i, new_tt);
0052       tkwt_vec.insert(tkwt_vec.begin() + i, new_tkwt);
0053       sum_Z_vec.insert(sum_Z_vec.begin() + i, 1.0);
0054       kmin.insert(kmin.begin() + i, 0);
0055       kmax.insert(kmax.begin() + i, 0);
0056     }
0057 
0058     unsigned int getSize() const { return zpca_vec.size(); }
0059 
0060     // has to be called everytime the items are modified
0061     void extractRaw() {
0062       zpca = &zpca_vec.front();
0063       dz2 = &dz2_vec.front();
0064       tkwt = &tkwt_vec.front();
0065       sum_Z = &sum_Z_vec.front();
0066     }
0067 
0068     // pointers to the first element of vectors, needed for vectorized code
0069     double *__restrict__ zpca;
0070     double *__restrict__ dz2;
0071     double *__restrict__ tkwt;
0072     double *__restrict__ sum_Z;
0073   };
0074 
0075   // internal data structure for clusters
0076   struct vertex_t {
0077     std::vector<double> zvtx_vec;  // z coordinate
0078     std::vector<double> rho_vec;   // vertex "mass" for mass-constrained clustering
0079     // --- temporary numbers, used during update
0080     std::vector<double> exp_arg_vec;
0081     std::vector<double> exp_vec;
0082     std::vector<double> sw_vec;
0083     std::vector<double> swz_vec;
0084     std::vector<double> se_vec;
0085     std::vector<double> swE_vec;
0086 
0087     unsigned int getSize() const { return zvtx_vec.size(); }
0088 
0089     void addItem(double new_zvtx, double new_rho) {
0090       zvtx_vec.push_back(new_zvtx);
0091       rho_vec.push_back(new_rho);
0092       exp_arg_vec.push_back(0.0);
0093       exp_vec.push_back(0.0);
0094       sw_vec.push_back(0.0);
0095       swz_vec.push_back(0.0);
0096       se_vec.push_back(0.0);
0097       swE_vec.push_back(0.0);
0098 
0099       extractRaw();
0100     }
0101 
0102     void insertItem(unsigned int k, double new_zvtx, double new_rho, track_t &tks) {
0103       zvtx_vec.insert(zvtx_vec.begin() + k, new_zvtx);
0104       rho_vec.insert(rho_vec.begin() + k, new_rho);
0105 
0106       exp_arg_vec.insert(exp_arg_vec.begin() + k, 0.0);
0107       exp_vec.insert(exp_vec.begin() + k, 0.0);
0108       sw_vec.insert(sw_vec.begin() + k, 0.0);
0109       swz_vec.insert(swz_vec.begin() + k, 0.0);
0110       se_vec.insert(se_vec.begin() + k, 0.0);
0111       swE_vec.insert(swE_vec.begin() + k, 0.0);
0112 
0113       // adjust vertex lists of tracks
0114       for (unsigned int i = 0; i < tks.getSize(); i++) {
0115         if (tks.kmin[i] > k) {
0116           tks.kmin[i]++;
0117         }
0118         if ((tks.kmax[i] >= k) || (tks.kmax[i] == tks.kmin[i])) {
0119           tks.kmax[i]++;
0120         }
0121       }
0122 
0123       extractRaw();
0124     }
0125 
0126     void removeItem(unsigned int k, track_t &tks) {
0127       zvtx_vec.erase(zvtx_vec.begin() + k);
0128       rho_vec.erase(rho_vec.begin() + k);
0129 
0130       exp_arg_vec.erase(exp_arg_vec.begin() + k);
0131       exp_vec.erase(exp_vec.begin() + k);
0132       sw_vec.erase(sw_vec.begin() + k);
0133       swz_vec.erase(swz_vec.begin() + k);
0134       se_vec.erase(se_vec.begin() + k);
0135       swE_vec.erase(swE_vec.begin() + k);
0136 
0137       // adjust vertex lists of tracks
0138       for (unsigned int i = 0; i < tks.getSize(); i++) {
0139         if (tks.kmax[i] > k) {
0140           tks.kmax[i]--;
0141         }
0142         if ((tks.kmin[i] > k) || (((tks.kmax[i] < (tks.kmin[i] + 1)) && (tks.kmin[i] > 0)))) {
0143           tks.kmin[i]--;
0144         }
0145       }
0146 
0147       extractRaw();
0148     }
0149 
0150     // pointers to the first element of vectors, needed for vectorized code
0151     double *__restrict__ zvtx;
0152     double *__restrict__ rho;
0153     double *__restrict__ exp_arg;
0154     double *__restrict__ exp;
0155     double *__restrict__ sw;
0156     double *__restrict__ swz;
0157     double *__restrict__ se;
0158     double *__restrict__ swE;
0159 
0160     // has to be called everytime the items are modified
0161     void extractRaw() {
0162       zvtx = &zvtx_vec.front();
0163       rho = &rho_vec.front();
0164       exp = &exp_vec.front();
0165       sw = &sw_vec.front();
0166       swz = &swz_vec.front();
0167       se = &se_vec.front();
0168       swE = &swE_vec.front();
0169       exp_arg = &exp_arg_vec.front();
0170     }
0171   };
0172 
0173   DAClusterizerInZ_vect(const edm::ParameterSet &conf);
0174 
0175   std::vector<std::vector<reco::TransientTrack> > clusterize(
0176       const std::vector<reco::TransientTrack> &tracks) const override;
0177 
0178   std::vector<TransientVertex> vertices(const std::vector<reco::TransientTrack> &tracks) const;
0179 
0180   track_t fill(const std::vector<reco::TransientTrack> &tracks) const;
0181 
0182   void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const;
0183 
0184   void clear_vtx_range(track_t &gtracks, vertex_t &gvertices) const;
0185 
0186   unsigned int thermalize(
0187       double beta, track_t &gtracks, vertex_t &gvertices, const double delta_max, const double rho0 = 0.) const;
0188 
0189   double update(
0190       double beta, track_t &gtracks, vertex_t &gvertices, const double rho0 = 0, const bool updateTc = false) const;
0191 
0192   void dump(
0193       const double beta, const vertex_t &y, const track_t &tks, const int verbosity = 0, const double rho0 = 0.) const;
0194   bool merge(vertex_t &y, track_t &tks, double &beta) const;
0195   bool purge(vertex_t &, track_t &, double &, const double) const;
0196   bool split(const double beta, track_t &t, vertex_t &y, double threshold = 1.) const;
0197 
0198   double beta0(const double betamax, track_t const &tks, vertex_t const &y) const;
0199   void verify(const vertex_t &v, const track_t &tks, unsigned int nv = 999999, unsigned int nt = 999999) const;
0200 
0201 private:
0202   double zdumpcenter_;
0203   double zdumpwidth_;
0204 
0205   double vertexSize_;
0206   unsigned int maxIterations_;
0207   double coolingFactor_;
0208   double betamax_;
0209   double betastop_;
0210   double dzCutOff_;
0211   double d0CutOff_;
0212 
0213   double mintrkweight_;
0214   double uniquetrkweight_;
0215   double uniquetrkminp_;
0216   double zmerge_;
0217   double betapurge_;
0218 
0219   unsigned int convergence_mode_;
0220   double delta_highT_;
0221   double delta_lowT_;
0222 
0223   double sel_zrange_;
0224   const double zrange_min_ = 0.1;  // smallest z-range to be included in a tracks cluster list
0225 };
0226 
0227 //#ifndef DAClusterizerInZ_vect_h
0228 #endif