Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:29:15

0001 #ifndef RecoVertex_PrimaryVertexProducer_DAClusterizerInZT_vect_h
0002 #define RecoVertex_PrimaryVertexProducer_DAClusterizerInZT_vect_h
0003 
0004 /**\class DAClusterizerInZT_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 "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include <vector>
0018 #include "DataFormats/Math/interface/Error.h"
0019 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0020 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0021 
0022 #include <memory>
0023 
0024 //#define USEVTXDT2
0025 
0026 class DAClusterizerInZT_vect final : public TrackClusterizerInZ {
0027 public:
0028   // internal data structure for tracks
0029   struct track_t {
0030     std::vector<double> zpca_vec;                  // z-coordinate at point of closest approach to the beamline
0031     std::vector<double> tpca_vec;                  // t-coordinate at point of closest approach to the beamline
0032     std::vector<double> dz2_vec;                   // square of the error of z(pca)
0033     std::vector<double> dt2_vec;                   // square of the error of t(pca)
0034     std::vector<double> sum_Z_vec;                 // track contribution to the partition function, Z
0035     std::vector<double> tkwt_vec;                  // track weight, close to 1.0 for most tracks
0036     std::vector<unsigned int> kmin;                // index of the first cluster within zrange
0037     std::vector<unsigned int> kmax;                // 1 + index of the last cluster within zrange
0038     std::vector<const reco::TransientTrack *> tt;  // a pointer to the Transient Track
0039 
0040     double osumtkwt;  // 1. / (sum of all track weights)
0041 
0042     void addItem(double new_zpca,
0043                  double new_tpca,
0044                  double new_dz2,
0045                  double new_dt2,
0046                  const reco::TransientTrack *new_tt,
0047                  double new_tkwt) {
0048       zpca_vec.push_back(new_zpca);
0049       tpca_vec.push_back(new_tpca);
0050       dz2_vec.push_back(new_dz2);
0051       dt2_vec.push_back(new_dt2);
0052       tt.push_back(new_tt);
0053       tkwt_vec.push_back(new_tkwt);
0054       sum_Z_vec.push_back(1.0);
0055       kmin.push_back(0);
0056       kmax.push_back(0);
0057     }
0058 
0059     void insertItem(unsigned int i,
0060                     double new_zpca,
0061                     double new_tpca,
0062                     double new_dz2,
0063                     double new_dt2,
0064                     const reco::TransientTrack *new_tt,
0065                     double new_tkwt) {
0066       zpca_vec.insert(zpca_vec.begin() + i, new_zpca);
0067       tpca_vec.insert(tpca_vec.begin() + i, new_tpca);
0068       dz2_vec.insert(dz2_vec.begin() + i, new_dz2);
0069       dt2_vec.insert(dt2_vec.begin() + i, new_dt2);
0070       tt.insert(tt.begin() + i, new_tt);
0071       tkwt_vec.insert(tkwt_vec.begin() + i, new_tkwt);
0072       sum_Z_vec.insert(sum_Z_vec.begin() + i, 1.0);
0073       kmin.insert(kmin.begin() + i, 0);
0074       kmax.insert(kmax.begin() + i, 0);
0075     }
0076 
0077     unsigned int getSize() const { return zpca_vec.size(); }
0078 
0079     // has to be called everytime the items are modified
0080     void extractRaw() {
0081       zpca = &zpca_vec.front();
0082       tpca = &tpca_vec.front();
0083       dz2 = &dz2_vec.front();
0084       dt2 = &dt2_vec.front();
0085       tkwt = &tkwt_vec.front();
0086       sum_Z = &sum_Z_vec.front();
0087     }
0088 
0089     // pointers to the first element of vectors, needed for vectorized code
0090     double *__restrict__ zpca;
0091     double *__restrict__ tpca;
0092     double *__restrict__ dz2;
0093     double *__restrict__ dt2;
0094     double *__restrict__ tkwt;
0095     double *__restrict__ sum_Z;
0096   };
0097 
0098   // internal data structure for clusters
0099   struct vertex_t {
0100     std::vector<double> zvtx_vec;  // z coordinate
0101     std::vector<double> tvtx_vec;  // t coordinate
0102     std::vector<double> rho_vec;   //  vertex "mass" for mass-constrained clustering
0103 #ifdef USEVTXDT2
0104     std::vector<double> dt2_vec;   // only used with vertex time uncertainties
0105     std::vector<double> sumw_vec;  // only used with vertex time uncertainties
0106 #endif
0107     // --- temporary numbers, used during update
0108     std::vector<double> exp_arg_vec;
0109     std::vector<double> exp_vec;
0110     std::vector<double> sw_vec;
0111     std::vector<double> swz_vec;
0112     std::vector<double> swt_vec;
0113     std::vector<double> se_vec;
0114     std::vector<double> nuz_vec;
0115     std::vector<double> nut_vec;
0116     std::vector<double> szz_vec;
0117     std::vector<double> stt_vec;
0118     std::vector<double> szt_vec;
0119 
0120     unsigned int getSize() const { return zvtx_vec.size(); }
0121 
0122     void addItem(double new_zvtx, double new_tvtx, double new_rho) {
0123       zvtx_vec.push_back(new_zvtx);
0124       tvtx_vec.push_back(new_tvtx);
0125       rho_vec.push_back(new_rho);
0126       exp_arg_vec.push_back(0.0);
0127       exp_vec.push_back(0.0);
0128       swz_vec.push_back(0.0);
0129       swt_vec.push_back(0.0);
0130       se_vec.push_back(0.0);
0131       nuz_vec.push_back(0.0);
0132       nut_vec.push_back(0.0);
0133       szz_vec.push_back(0.0);
0134       stt_vec.push_back(0.0);
0135       szt_vec.push_back(0.0);
0136 #ifdef USEVTXDT2
0137       dt2_vec.push_back(0.0);
0138       sumw_vec.push_back(0.0);
0139 #endif
0140 
0141       extractRaw();
0142     }
0143 
0144     void insertItem(unsigned int k, double new_zvtx, double new_tvtx, double new_rho, track_t &tks) {
0145       zvtx_vec.insert(zvtx_vec.begin() + k, new_zvtx);
0146       tvtx_vec.insert(tvtx_vec.begin() + k, new_tvtx);
0147       rho_vec.insert(rho_vec.begin() + k, new_rho);
0148       exp_arg_vec.insert(exp_arg_vec.begin() + k, 0.0);
0149       exp_vec.insert(exp_vec.begin() + k, 0.0);
0150       swz_vec.insert(swz_vec.begin() + k, 0.0);
0151       swt_vec.insert(swt_vec.begin() + k, 0.0);
0152       se_vec.insert(se_vec.begin() + k, 0.0);
0153       nuz_vec.insert(nuz_vec.begin() + k, 0.0);
0154       nut_vec.insert(nut_vec.begin() + k, 0.0);
0155       szz_vec.insert(szz_vec.begin() + k, 0.0);
0156       stt_vec.insert(stt_vec.begin() + k, 0.0);
0157       szt_vec.insert(szt_vec.begin() + k, 0.0);
0158 #ifdef USEVTXDT2
0159       dt2_vec.insert(dt2_vec.begin() + k, 0.0);
0160       sumw_vec.insert(sumw_vec.begin() + k, 0.0);
0161 #endif
0162 
0163       // adjust vertex lists of tracks
0164       for (unsigned int i = 0; i < tks.getSize(); i++) {
0165         if (tks.kmin[i] > k) {
0166           tks.kmin[i]++;
0167         }
0168         if ((tks.kmax[i] >= k) || (tks.kmax[i] == tks.kmin[i])) {
0169           tks.kmax[i]++;
0170         }
0171       }
0172 
0173       extractRaw();
0174     }
0175 
0176     void removeItem(unsigned int k, track_t &tks) {
0177       zvtx_vec.erase(zvtx_vec.begin() + k);
0178       tvtx_vec.erase(tvtx_vec.begin() + k);
0179       rho_vec.erase(rho_vec.begin() + k);
0180       exp_arg_vec.erase(exp_arg_vec.begin() + k);
0181       exp_vec.erase(exp_vec.begin() + k);
0182       swz_vec.erase(swz_vec.begin() + k);
0183       swt_vec.erase(swt_vec.begin() + k);
0184       se_vec.erase(se_vec.begin() + k);
0185       nuz_vec.erase(nuz_vec.begin() + k);
0186       nut_vec.erase(nut_vec.begin() + k);
0187       szz_vec.erase(szz_vec.begin() + k);
0188       stt_vec.erase(stt_vec.begin() + k);
0189       szt_vec.erase(szt_vec.begin() + k);
0190 #ifdef USEVTXDT2
0191       dt2_vec.erase(dt2_vec.begin() + k);
0192       sumw_vec.erase(sumw_vec.begin() + k);
0193 #endif
0194 
0195       // adjust vertex lists of tracks
0196       for (unsigned int i = 0; i < tks.getSize(); i++) {
0197         if (tks.kmax[i] > k) {
0198           tks.kmax[i]--;
0199         }
0200         if ((tks.kmin[i] > k) || (((tks.kmax[i] < (tks.kmin[i] + 1)) && (tks.kmin[i] > 0)))) {
0201           tks.kmin[i]--;
0202         }
0203       }
0204 
0205       extractRaw();
0206     }
0207 
0208     unsigned int insertOrdered(double zvtx, double tvtx, double rho, track_t &tks) {
0209       // insert a new cluster according to it's z-position, return the index at which it was inserted
0210 
0211       unsigned int k = 0;
0212       for (; k < getSize(); k++) {
0213         if (zvtx < zvtx_vec[k])
0214           break;
0215       }
0216       insertItem(k, zvtx, tvtx, rho, tks);
0217       return k;
0218     }
0219 
0220     // pointers to the first element of vectors, needed for vectorized code
0221     double *__restrict__ zvtx;
0222     double *__restrict__ tvtx;
0223     double *__restrict__ rho;
0224     double *__restrict__ exp_arg;
0225     double *__restrict__ exp;
0226     double *__restrict__ swt;
0227     double *__restrict__ swz;
0228     double *__restrict__ se;
0229     double *__restrict__ nuz;
0230     double *__restrict__ nut;
0231     double *__restrict__ szz;
0232     double *__restrict__ stt;
0233     double *__restrict__ szt;
0234 #ifdef USEVTXDT2
0235     double *__restrict__ dt2;
0236     double *__restrict__ sumw;
0237 #endif
0238 
0239     // has to be called everytime the items are modified
0240     void extractRaw() {
0241       zvtx = &zvtx_vec.front();
0242       tvtx = &tvtx_vec.front();
0243       rho = &rho_vec.front();
0244       exp_arg = &exp_arg_vec.front();
0245       exp = &exp_vec.front();
0246       swz = &swz_vec.front();
0247       swt = &swt_vec.front();
0248       se = &se_vec.front();
0249       nuz = &nuz_vec.front();
0250       nut = &nut_vec.front();
0251       szz = &szz_vec.front();
0252       stt = &stt_vec.front();
0253       szt = &szt_vec.front();
0254 #ifdef USEVTXDT2
0255       dt2 = &dt2_vec.front();
0256       sumw = &sumw_vec.front();
0257 #endif
0258     }
0259   };
0260 
0261   DAClusterizerInZT_vect(const edm::ParameterSet &conf);
0262 
0263   static void fillPSetDescription(edm::ParameterSetDescription &desc);
0264 
0265   std::vector<std::vector<reco::TransientTrack> > clusterize(
0266       const std::vector<reco::TransientTrack> &tracks) const override;
0267 
0268   std::vector<TransientVertex> vertices(const std::vector<reco::TransientTrack> &tracks) const override;
0269 
0270   track_t fill(const std::vector<reco::TransientTrack> &tracks) const;
0271 
0272   void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const;
0273 
0274   void clear_vtx_range(track_t &gtracks, vertex_t &gvertices) const;
0275 
0276   unsigned int thermalize(
0277       double beta, track_t &gtracks, vertex_t &gvertices, const double delta_max, const double rho0 = 0.) const;
0278 
0279   double update(
0280       double beta, track_t &gtracks, vertex_t &gvertices, const double rho0 = 0, const bool updateTc = false) const;
0281 
0282   void dump(const double beta, const vertex_t &y, const track_t &tks, const int verbosity = 0) const;
0283   bool zorder(vertex_t &y) const;
0284   bool find_nearest(double z, double t, vertex_t &y, unsigned int &k_min, double dz, double dt) const;
0285   bool merge(vertex_t &, track_t &, double &beta) const;
0286   bool purge(vertex_t &, track_t &, double &, const double) const;
0287   bool split(const double beta, track_t &t, vertex_t &y, double threshold = 1.) const;
0288 
0289   double beta0(const double betamax, track_t const &tks, vertex_t const &y) const;
0290 
0291   double get_Tc(const vertex_t &y, int k) const;
0292   void verify(const vertex_t &v, const track_t &tks, unsigned int nv = 999999, unsigned int nt = 999999) const;
0293 
0294 private:
0295   double zdumpcenter_;
0296   double zdumpwidth_;
0297 
0298   double vertexSize_;
0299   double vertexSizeTime_;
0300   unsigned int maxIterations_;
0301   double coolingFactor_;
0302   double betamax_;
0303   double betastop_;
0304   double dzCutOff_;
0305   double d0CutOff_;
0306   double dtCutOff_;
0307   double t0Max_;
0308 
0309   double mintrkweight_;
0310   double uniquetrkweight_;
0311   double uniquetrkminp_;
0312   double zmerge_;
0313   double tmerge_;
0314   double betapurge_;
0315 
0316   unsigned int convergence_mode_;
0317   double delta_highT_;
0318   double delta_lowT_;
0319 
0320   double sel_zrange_;
0321   const double zrange_min_ = 0.1;  // smallest z-range to be included in a tracks cluster list
0322 };
0323 
0324 //#ifndef DAClusterizerInZT_vect_h
0325 #endif