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
0005
0006
0007
0008
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
0025
0026 class DAClusterizerInZT_vect final : public TrackClusterizerInZ {
0027 public:
0028
0029 struct track_t {
0030 std::vector<double> zpca_vec;
0031 std::vector<double> tpca_vec;
0032 std::vector<double> dz2_vec;
0033 std::vector<double> dt2_vec;
0034 std::vector<double> sum_Z_vec;
0035 std::vector<double> tkwt_vec;
0036 std::vector<unsigned int> kmin;
0037 std::vector<unsigned int> kmax;
0038 std::vector<const reco::TransientTrack *> tt;
0039
0040 double osumtkwt;
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
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
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
0099 struct vertex_t {
0100 std::vector<double> zvtx_vec;
0101 std::vector<double> tvtx_vec;
0102 std::vector<double> rho_vec;
0103 #ifdef USEVTXDT2
0104 std::vector<double> dt2_vec;
0105 std::vector<double> sumw_vec;
0106 #endif
0107
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
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
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
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
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
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 >racks, vertex_t &gvertices) const;
0273
0274 void clear_vtx_range(track_t >racks, vertex_t &gvertices) const;
0275
0276 unsigned int thermalize(
0277 double beta, track_t >racks, vertex_t &gvertices, const double delta_max, const double rho0 = 0.) const;
0278
0279 double update(
0280 double beta, track_t >racks, 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;
0322 };
0323
0324
0325 #endif