Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-25 22:35:19 

0001 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h"
0002 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0005 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0006 
0007 #include <cmath>
0008 #include <cassert>
0009 #include <limits>
0010 #include <iomanip>
0011 #include "FWCore/Utilities/interface/isFinite.h"
0012 #include "vdt/vdtMath.h"
0013 
0014 using namespace std;
0015 
0016 //#define DEBUG
0017 #ifdef DEBUG
0018 #define DEBUGLEVEL 0
0019 #endif
0020 
0021 DAClusterizerInZT_vect::DAClusterizerInZT_vect(const edm::ParameterSet& conf) {
0022   // hardcoded parameters
0023   maxIterations_ = 1000;
0024   mintrkweight_ = 0.5;
0025 
0026   // configurable debug output
0027 #ifdef DEBUG
0028   zdumpcenter_ = conf.getUntrackedParameter<double>("zdumpcenter", 0.);
0029   zdumpwidth_ = conf.getUntrackedParameter<double>("zdumpwidth", 20.);
0030 #endif
0031 
0032   // configurable parameters
0033   double minT = conf.getParameter<double>("Tmin");
0034   double purgeT = conf.getParameter<double>("Tpurge");
0035   double stopT = conf.getParameter<double>("Tstop");
0036   vertexSize_ = conf.getParameter<double>("vertexSize");
0037   vertexSizeTime_ = conf.getParameter<double>("vertexSizeTime");
0038   coolingFactor_ = conf.getParameter<double>("coolingFactor");
0039   d0CutOff_ = conf.getParameter<double>("d0CutOff");
0040   dzCutOff_ = conf.getParameter<double>("dzCutOff");
0041   dtCutOff_ = conf.getParameter<double>("dtCutOff");
0042   t0Max_ = conf.getParameter<double>("t0Max");
0043   uniquetrkweight_ = conf.getParameter<double>("uniquetrkweight");
0044   uniquetrkminp_ = conf.getParameter<double>("uniquetrkminp");
0045   zmerge_ = conf.getParameter<double>("zmerge");
0046   tmerge_ = conf.getParameter<double>("tmerge");
0047   sel_zrange_ = conf.getParameter<double>("zrange");
0048   convergence_mode_ = conf.getParameter<int>("convergence_mode");
0049   delta_lowT_ = conf.getParameter<double>("delta_lowT");
0050   delta_highT_ = conf.getParameter<double>("delta_highT");
0051 
0052 #ifdef DEBUG
0053   std::cout << "DAClusterizerInZT_vect: mintrkweight = " << mintrkweight_ << std::endl;
0054   std::cout << "DAClusterizerInZT_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl;
0055   std::cout << "DAClusterizerInZT_vect: uniquetrkminp = " << uniquetrkminp_ << std::endl;
0056   std::cout << "DAClusterizerInZT_vect: zmerge = " << zmerge_ << std::endl;
0057   std::cout << "DAClusterizerInZT_vect: tmerge = " << tmerge_ << std::endl;
0058   std::cout << "DAClusterizerInZT_vect: Tmin = " << minT << std::endl;
0059   std::cout << "DAClusterizerInZT_vect: Tpurge = " << purgeT << std::endl;
0060   std::cout << "DAClusterizerInZT_vect: Tstop = " << stopT << std::endl;
0061   std::cout << "DAClusterizerInZT_vect: vertexSize = " << vertexSize_ << std::endl;
0062   std::cout << "DAClusterizerInZT_vect: vertexSizeTime = " << vertexSizeTime_ << std::endl;
0063   std::cout << "DAClusterizerInZT_vect: coolingFactor = " << coolingFactor_ << std::endl;
0064   std::cout << "DAClusterizerInZT_vect: d0CutOff = " << d0CutOff_ << std::endl;
0065   std::cout << "DAClusterizerInZT_vect: dzCutOff = " << dzCutOff_ << std::endl;
0066   std::cout << "DAClusterizerInZT_vect: dtCutoff = " << dtCutOff_ << std::endl;
0067   std::cout << "DAClusterizerInZT_vect: zrange = " << sel_zrange_ << std::endl;
0068   std::cout << "DAClusterizerinZT_vect: convergence mode = " << convergence_mode_ << std::endl;
0069   std::cout << "DAClusterizerinZT_vect: delta_highT = " << delta_highT_ << std::endl;
0070   std::cout << "DAClusterizerinZT_vect: delta_lowT = " << delta_lowT_ << std::endl;
0071   std::cout << "DAClusterizerinZT_vect: DEBUGLEVEL " << DEBUGLEVEL << std::endl;
0072 #endif
0073 
0074   if (convergence_mode_ > 1) {
0075     edm::LogWarning("DAClusterizerinZT_vect")
0076         << "DAClusterizerInZT_vect: invalid convergence_mode " << convergence_mode_ << "  reset to default " << 0;
0077     convergence_mode_ = 0;
0078   }
0079 
0080   if (minT == 0) {
0081     betamax_ = 1.0;
0082     edm::LogWarning("DAClusterizerinZT_vect")
0083         << "DAClusterizerInZT_vect: invalid Tmin " << minT << "  reset to default " << 1. / betamax_;
0084   } else {
0085     betamax_ = 1. / minT;
0086   }
0087 
0088   if ((purgeT > minT) || (purgeT == 0)) {
0089     edm::LogWarning("DAClusterizerinZT_vect")
0090         << "DAClusterizerInZT_vect: invalid Tpurge " << purgeT << "  set to " << minT;
0091     purgeT = minT;
0092   }
0093   betapurge_ = 1. / purgeT;
0094 
0095   if ((stopT > purgeT) || (stopT == 0)) {
0096     edm::LogWarning("DAClusterizerinZT_vect")
0097         << "DAClusterizerInZT_vect: invalid Tstop " << stopT << "  set to  " << max(1., purgeT);
0098     stopT = max(1., purgeT);
0099   }
0100   betastop_ = 1. / stopT;
0101 }
0102 
0103 namespace {
0104   inline double local_exp(double const& inp) { return vdt::fast_exp(inp); }
0105 
0106   inline void local_exp_list_range(double const* __restrict__ arg_inp,
0107                                    double* __restrict__ arg_out,
0108                                    const int kmin,
0109                                    const int kmax) {
0110     for (auto i = kmin; i != kmax; ++i)
0111       arg_out[i] = vdt::fast_exp(arg_inp[i]);
0112   }
0113 
0114 }  // namespace
0115 
0116 void DAClusterizerInZT_vect::verify(const vertex_t& v, const track_t& tks, unsigned int nv, unsigned int nt) const {
0117   if (nv == 999999) {
0118     nv = v.getSize();
0119   } else {
0120     assert(nv == v.getSize());
0121   }
0122 
0123   if (nt == 999999) {
0124     nt = tks.getSize();
0125   } else {
0126     assert(nt == tks.getSize());
0127   }
0128 
0129   // clusters
0130   assert(v.zvtx_vec.size() == nv);
0131   assert(v.tvtx_vec.size() == nv);
0132 #ifdef USEVTXDT2
0133   assert(v.dt2_vec.size() == nv);
0134   assert(v.sumw_vec.size() == nv);
0135 #endif
0136   assert(v.rho_vec.size() == nv);
0137   assert(v.swz_vec.size() == nv);
0138   assert(v.swt_vec.size() == nv);
0139   assert(v.exp_arg_vec.size() == nv);
0140   assert(v.exp_vec.size() == nv);
0141   assert(v.se_vec.size() == nv);
0142   assert(v.nuz_vec.size() == nv);
0143   assert(v.nut_vec.size() == nv);
0144   assert(v.szz_vec.size() == nv);
0145   assert(v.stt_vec.size() == nv);
0146   assert(v.szt_vec.size() == nv);
0147 
0148   assert(v.zvtx == &v.zvtx_vec.front());
0149   assert(v.tvtx == &v.tvtx_vec.front());
0150   assert(v.rho == &v.rho_vec.front());
0151   assert(v.exp_arg == &v.exp_arg_vec.front());
0152   assert(v.swz == &v.swz_vec.front());
0153   assert(v.swt == &v.swt_vec.front());
0154   assert(v.se == &v.se_vec.front());
0155   assert(v.nuz == &v.nuz_vec.front());
0156   assert(v.nut == &v.nut_vec.front());
0157   assert(v.szz == &v.szz_vec.front());
0158   assert(v.stt == &v.stt_vec.front());
0159   assert(v.szt == &v.szt_vec.front());
0160 #ifdef USEVTXDT2
0161   assert(v.sumw == &v.sumw_vec.front());
0162   assert(v.dt2 == &v.dt2_vec.front());
0163 #endif
0164 
0165   for (unsigned int k = 0; k < nv - 1; k++) {
0166     if (v.zvtx[k] <= v.zvtx[k + 1])
0167       continue;
0168     cout << " ZT, cluster z-ordering assertion failure   z[" << k << "] =" << v.zvtx[k] << "    z[" << k + 1
0169          << "] =" << v.zvtx[k + 1] << endl;
0170   }
0171   //for(unsigned int k=0; k< nv-1; k++){
0172   //assert( v.z[k] <= v.z[k+1]);
0173   //}
0174 
0175   // tracks
0176   assert(nt == tks.zpca_vec.size());
0177   assert(nt == tks.tpca_vec.size());
0178   assert(nt == tks.dz2_vec.size());
0179   assert(nt == tks.dt2_vec.size());
0180   assert(nt == tks.tt.size());
0181   assert(nt == tks.tkwt_vec.size());
0182   assert(nt == tks.sum_Z_vec.size());
0183   assert(nt == tks.kmin.size());
0184   assert(nt == tks.kmax.size());
0185 
0186   assert(tks.zpca == &tks.zpca_vec.front());
0187   assert(tks.tpca == &tks.tpca_vec.front());
0188   assert(tks.dz2 == &tks.dz2_vec.front());
0189   assert(tks.dt2 == &tks.dt2_vec.front());
0190   assert(tks.tkwt == &tks.tkwt_vec.front());
0191   assert(tks.sum_Z == &tks.sum_Z_vec.front());
0192 
0193   for (unsigned int i = 0; i < nt; i++) {
0194     if ((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv))
0195       continue;
0196     cout << "track vertex range assertion failure" << i << "/" << nt << "   kmin,kmax=" << tks.kmin[i] << ", "
0197          << tks.kmax[i] << "  nv=" << nv << endl;
0198   }
0199 
0200   for (unsigned int i = 0; i < nt; i++) {
0201     assert((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv));
0202   }
0203 }
0204 
0205 //todo: use r-value possibility of c++11 here
0206 DAClusterizerInZT_vect::track_t DAClusterizerInZT_vect::fill(const vector<reco::TransientTrack>& tracks) const {
0207   // prepare track data for clustering
0208   track_t tks;
0209   double sumtkwt = 0.;
0210 
0211   for (const auto& tk : tracks) {
0212     if (!tk.isValid())
0213       continue;
0214     double t_tkwt = 1.;
0215     double t_z = tk.stateAtBeamLine().trackStateAtPCA().position().z();
0216     double t_t = tk.timeExt();
0217 
0218     if (std::fabs(t_z) > 1000.)
0219       continue;
0220     /*  for comparison with 1d clustering, keep such tracks without timing info, see below
0221     if (std::abs(t_t) > t0Max_)
0222     continue;
0223     */
0224 
0225     auto const& t_mom = tk.stateAtBeamLine().trackStateAtPCA().momentum();
0226     //  get the beam-spot
0227     reco::BeamSpot beamspot = tk.stateAtBeamLine().beamSpot();
0228     double t_dz2 = std::pow(tk.track().dzError(), 2)  // track errror
0229                    + (std::pow(beamspot.BeamWidthX() * t_mom.x(), 2) + std::pow(beamspot.BeamWidthY() * t_mom.y(), 2)) *
0230                          std::pow(t_mom.z(), 2) / std::pow(t_mom.perp2(), 2)  // beam spot width
0231                    + std::pow(vertexSize_, 2);  // intrinsic vertex size, safer for outliers and short lived decays
0232     t_dz2 = 1. / t_dz2;
0233     if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits<double>::min()) {
0234       edm::LogWarning("DAClusterizerinZT_vect") << "rejected track t_dz2 " << t_dz2;
0235       continue;
0236     }
0237 
0238     double t_dt2 =
0239         std::pow(tk.dtErrorExt(), 2.) +
0240         std::pow(vertexSizeTime_, 2.);  // the ~injected~ timing error, need to add a small minimum vertex size in time
0241     if ((tk.dtErrorExt() >= TransientTrackBuilder::defaultInvalidTrackTimeReso) || (std::abs(t_t) > t0Max_)) {
0242       t_dt2 = 0;  // tracks with no time measurement
0243     } else {
0244       t_dt2 = 1. / t_dt2;
0245       if (edm::isNotFinite(t_dt2) || t_dt2 < std::numeric_limits<double>::min()) {
0246         edm::LogWarning("DAClusterizerinZT_vect") << "rejected track t_dt2 " << t_dt2;
0247         continue;
0248       }
0249     }
0250 
0251     if (d0CutOff_ > 0) {
0252       Measurement1D atIP = tk.stateAtBeamLine().transverseImpactParameter();  // error contains beamspot
0253       t_tkwt = 1. / (1. + local_exp(std::pow(atIP.value() / atIP.error(), 2) -
0254                                     std::pow(d0CutOff_, 2)));  // reduce weight for high ip tracks
0255       if (edm::isNotFinite(t_tkwt) || t_tkwt < std::numeric_limits<double>::epsilon()) {
0256         edm::LogWarning("DAClusterizerinZT_vect") << "rejected track t_tkwt " << t_tkwt;
0257         continue;  // usually is > 0.99
0258       }
0259     }
0260 
0261     tks.addItem(t_z, t_t, t_dz2, t_dt2, &tk, t_tkwt);
0262     sumtkwt += t_tkwt;
0263   }
0264 
0265   if (sumtkwt > 0) {
0266     tks.extractRaw();
0267     tks.osumtkwt = 1. / sumtkwt;
0268   } else {
0269     tks.osumtkwt = 0.;
0270   }
0271 
0272 #ifdef DEBUG
0273   if (DEBUGLEVEL > 0) {
0274     std::cout << "Track count (ZT) filled " << tks.getSize() << " initial " << tracks.size() << std::endl;
0275   }
0276 #endif
0277 
0278   return tks;
0279 }
0280 
0281 namespace {
0282   inline double Eik(double t_z, double k_z, double t_dz2, double t_t, double k_t, double t_dt2) {
0283     return std::pow(t_z - k_z, 2) * t_dz2 + std::pow(t_t - k_t, 2) * t_dt2;
0284   }
0285 }  // namespace
0286 
0287 void DAClusterizerInZT_vect::set_vtx_range(double beta, track_t& gtracks, vertex_t& gvertices) const {
0288   const unsigned int nv = gvertices.getSize();
0289   const unsigned int nt = gtracks.getSize();
0290 
0291   if (nv == 0) {
0292     edm::LogWarning("DAClusterizerinZT_vect") << "empty cluster list in set_vtx_range";
0293     return;
0294   }
0295 
0296   for (auto itrack = 0U; itrack < nt; ++itrack) {
0297     double zrange = max(sel_zrange_ / sqrt(beta * gtracks.dz2[itrack]), zrange_min_);
0298 
0299     double zmin = gtracks.zpca[itrack] - zrange;
0300     unsigned int kmin = min(nv - 1, gtracks.kmin[itrack]);
0301     // find the smallest vertex_z that is larger than zmin
0302     if (gvertices.zvtx[kmin] > zmin) {
0303       while ((kmin > 0) && (gvertices.zvtx[kmin - 1] > zmin)) {
0304         kmin--;
0305       }
0306     } else {
0307       while ((kmin < (nv - 1)) && (gvertices.zvtx[kmin] < zmin)) {
0308         kmin++;
0309       }
0310     }
0311 
0312     double zmax = gtracks.zpca[itrack] + zrange;
0313     unsigned int kmax = min(nv - 1, gtracks.kmax[itrack] - 1);
0314     // note: kmax points to the last vertex in the range, while gtracks.kmax points to the entry BEHIND the last vertex
0315     // find the largest vertex_z that is smaller than zmax
0316     if (gvertices.zvtx[kmax] < zmax) {
0317       while ((kmax < (nv - 1)) && (gvertices.zvtx[kmax + 1] < zmax)) {
0318         kmax++;
0319       }
0320     } else {
0321       while ((kmax > 0) && (gvertices.zvtx[kmax] > zmax)) {
0322         kmax--;
0323       }
0324     }
0325 
0326     if (kmin <= kmax) {
0327       gtracks.kmin[itrack] = kmin;
0328       gtracks.kmax[itrack] = kmax + 1;
0329     } else {
0330       gtracks.kmin[itrack] = max(0U, min(kmin, kmax));
0331       gtracks.kmax[itrack] = min(nv, max(kmin, kmax) + 1);
0332     }
0333 
0334 #ifdef DEBUG
0335     if (gtracks.kmin[itrack] >= gtracks.kmax[itrack]) {
0336       cout << "set_vtx_range trk = " << itrack << "  kmin,kmax=" << kmin << "," << kmax
0337            << " gtrack.kmin,kmax = " << gtracks.kmin[itrack] << "," << gtracks.kmax[itrack] << " zrange = " << zrange
0338            << endl;
0339     }
0340 #endif
0341   }
0342 }
0343 
0344 void DAClusterizerInZT_vect::clear_vtx_range(track_t& gtracks, vertex_t& gvertices) const {
0345   const unsigned int nt = gtracks.getSize();
0346   const unsigned int nv = gvertices.getSize();
0347   for (auto itrack = 0U; itrack < nt; ++itrack) {
0348     gtracks.kmin[itrack] = 0;
0349     gtracks.kmax[itrack] = nv;
0350   }
0351 }
0352 
0353 double DAClusterizerInZT_vect::update(
0354     double beta, track_t& gtracks, vertex_t& gvertices, const double rho0, const bool updateTc) const {
0355   //update weights and vertex positions
0356   // returns the maximum of changes of vertex positions
0357   // sums needed for Tc are only updated if updateTC == true
0358 
0359   const unsigned int nt = gtracks.getSize();
0360   const unsigned int nv = gvertices.getSize();
0361   auto osumtkwt = gtracks.osumtkwt;
0362 
0363   double Z_init = 0;
0364   if (rho0 > 0) {
0365     Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);  // cut-off
0366   }
0367 
0368   // define kernels
0369 
0370   auto kernel_calc_exp_arg_range = [beta](const unsigned int itrack,
0371                                           track_t const& tracks,
0372                                           vertex_t const& vertices,
0373                                           const unsigned int kmin,
0374                                           const unsigned int kmax) {
0375     const auto track_z = tracks.zpca[itrack];
0376     const auto track_t = tracks.tpca[itrack];
0377     const auto botrack_dz2 = -beta * tracks.dz2[itrack];
0378 #ifndef USEVTXDT2
0379     const auto botrack_dt2 = -beta * tracks.dt2[itrack];
0380 #endif
0381     // auto-vectorized
0382     for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
0383       const auto mult_resz = track_z - vertices.zvtx[ivertex];
0384       const auto mult_rest = track_t - vertices.tvtx[ivertex];
0385 #ifdef USEVTXDT2
0386       const auto botrack_dt2 =
0387           -beta * tracks.dt2[itrack] * vertices.dt2[ivertex] / (tracks.dt2[itrack] + vertices.dt2[ivertex] + 1.e-10);
0388 #endif
0389       vertices.exp_arg[ivertex] = botrack_dz2 * (mult_resz * mult_resz) + botrack_dt2 * (mult_rest * mult_rest);
0390     }
0391   };
0392 
0393   auto kernel_add_Z_range = [Z_init](
0394                                 vertex_t const& vertices, const unsigned int kmin, const unsigned int kmax) -> double {
0395     double ZTemp = Z_init;
0396     for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
0397       ZTemp += vertices.rho[ivertex] * vertices.exp[ivertex];
0398     }
0399     return ZTemp;
0400   };
0401 
0402   auto kernel_calc_normalization_range = [updateTc](const unsigned int track_num,
0403                                                     track_t& tracks,
0404                                                     vertex_t& vertices,
0405                                                     const unsigned int kmin,
0406                                                     const unsigned int kmax) {
0407     auto tmp_trk_tkwt = tracks.tkwt[track_num];
0408     auto o_trk_sum_Z = 1. / tracks.sum_Z[track_num];
0409     auto o_trk_err_z = tracks.dz2[track_num];
0410     auto o_trk_err_t = tracks.dt2[track_num];
0411     auto tmp_trk_z = tracks.zpca[track_num];
0412     auto tmp_trk_t = tracks.tpca[track_num];
0413     // auto-vectorized
0414     if (updateTc) {
0415 #pragma GCC ivdep
0416       for (unsigned int k = kmin; k < kmax; ++k) {
0417         // parens are important for numerical stability
0418         vertices.se[k] += tmp_trk_tkwt * (vertices.exp[k] * o_trk_sum_Z);
0419         const auto w = tmp_trk_tkwt * (vertices.rho[k] * vertices.exp[k] * o_trk_sum_Z);  // p_{ik}
0420         const auto wz = w * o_trk_err_z;
0421         const auto wt = w * o_trk_err_t;
0422 #ifdef USEVTXDT2
0423         vertices.sumw[k] += w;  // for vtxdt2
0424 #endif
0425         vertices.nuz[k] += wz;
0426         vertices.nut[k] += wt;
0427         vertices.swz[k] += wz * tmp_trk_z;
0428         vertices.swt[k] += wt * tmp_trk_t;
0429         // the following lines are only needed for Tc updates
0430         const auto dsz = (tmp_trk_z - vertices.zvtx[k]) * o_trk_err_z;
0431         const auto dst = (tmp_trk_t - vertices.tvtx[k]) * o_trk_err_t;
0432         vertices.szz[k] += w * dsz * dsz;
0433         vertices.stt[k] += w * dst * dst;
0434         vertices.szt[k] += w * dsz * dst;
0435       }
0436     } else {
0437 #pragma GCC ivdep
0438       for (unsigned int k = kmin; k < kmax; ++k) {
0439         vertices.se[k] += tmp_trk_tkwt * (vertices.exp[k] * o_trk_sum_Z);
0440         const auto w = tmp_trk_tkwt * (vertices.rho[k] * vertices.exp[k] * o_trk_sum_Z);  // p_{ik}
0441         const auto wz = w * o_trk_err_z;
0442         const auto wt = w * o_trk_err_t;
0443 #ifdef USEVTXDT2
0444         vertices.sumw[k] += w;  // for vtxdt2
0445 #endif
0446         vertices.nuz[k] += wz;
0447         vertices.nut[k] += wt;
0448         vertices.swz[k] += wz * tmp_trk_z;
0449         vertices.swt[k] += wt * tmp_trk_t;
0450       }
0451     }
0452   };
0453 
0454   for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
0455     gvertices.se[ivertex] = 0.0;
0456     gvertices.nuz[ivertex] = 0.0;
0457     gvertices.nut[ivertex] = 0.0;
0458     gvertices.swz[ivertex] = 0.0;
0459     gvertices.swt[ivertex] = 0.0;
0460     gvertices.szz[ivertex] = 0.0;
0461     gvertices.stt[ivertex] = 0.0;
0462     gvertices.szt[ivertex] = 0.0;
0463 #ifdef USEVTXDT2
0464     gvertices.sumw[ivertex] = 0.0;
0465 #endif
0466   }
0467 
0468   // loop over tracks
0469   for (auto itrack = 0U; itrack < nt; ++itrack) {
0470     const unsigned int kmin = gtracks.kmin[itrack];
0471     const unsigned int kmax = gtracks.kmax[itrack];
0472 
0473 #ifdef DEBUG
0474     assert((kmin < kmax) && (kmax <= nv));
0475     assert(itrack < gtracks.sum_Z_vec.size());
0476 #endif
0477 
0478     kernel_calc_exp_arg_range(itrack, gtracks, gvertices, kmin, kmax);
0479     local_exp_list_range(gvertices.exp_arg, gvertices.exp, kmin, kmax);
0480     gtracks.sum_Z[itrack] = kernel_add_Z_range(gvertices, kmin, kmax);
0481 
0482     if (edm::isNotFinite(gtracks.sum_Z[itrack]))
0483       gtracks.sum_Z[itrack] = 0.0;
0484 
0485     if (gtracks.sum_Z[itrack] > 1.e-100) {
0486       kernel_calc_normalization_range(itrack, gtracks, gvertices, kmin, kmax);
0487     }
0488   }
0489 
0490   // now update z, t, and rho
0491   auto kernel_calc_zt = [osumtkwt, nv](vertex_t& vertices) -> double {
0492     double delta = 0;
0493 
0494     // does not vectorize
0495     for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
0496       if (vertices.nuz[ivertex] > 0.) {
0497         auto znew = vertices.swz[ivertex] / vertices.nuz[ivertex];
0498         delta = max(std::abs(vertices.zvtx[ivertex] - znew), delta);
0499         vertices.zvtx[ivertex] = znew;
0500       }
0501 
0502       if (vertices.nut[ivertex] > 0.) {
0503         auto tnew = vertices.swt[ivertex] / vertices.nut[ivertex];
0504         //delta = max(std::abs(vertices.t[ ivertex ] - tnew), delta); // FIXME
0505         vertices.tvtx[ivertex] = tnew;
0506 #ifdef USEVTXDT2
0507         vertices.dt2[ivertex] = vertices.nut[ivertex] / vertices.sumw[ivertex];
0508 #endif
0509       } else {
0510         // FIXME
0511         // apparently this cluster has not timing info attached
0512         // this should be taken into account somehow, otherwise we might fail to attach
0513         // new tracks that do have timing information
0514         vertices.tvtx[ivertex] = 0;
0515 #ifdef USEVTXDT2
0516         vertices.dt2[ivertex] = 0;
0517 #endif
0518       }
0519 
0520 #ifdef DEBUG
0521       if ((vertices.nut[ivertex] <= 0.) && (vertices.nuz[ivertex] <= 0.)) {
0522         edm::LogInfo("sumw") << "invalid sum of weights in fit: " << endl;
0523         std::cout << " a cluster melted away ?  "
0524                   << "  zk=" << vertices.zvtx[ivertex] << "  rho=" << vertices.rho[ivertex]
0525                   << " sumw(z,t) =" << vertices.nuz[ivertex] << "," << vertices.nut[ivertex] << endl;
0526         // FIXME: discard this cluster
0527       }
0528 #endif
0529     }
0530 
0531     for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
0532       vertices.rho[ivertex] = vertices.rho[ivertex] * vertices.se[ivertex] * osumtkwt;
0533     }
0534 
0535     return delta;
0536   };
0537 
0538   double delta = kernel_calc_zt(gvertices);
0539 
0540   if (zorder(gvertices)) {
0541     set_vtx_range(beta, gtracks, gvertices);
0542   };
0543 
0544   // return how much the prototypes moved
0545   return delta;
0546 }
0547 
0548 bool DAClusterizerInZT_vect::zorder(vertex_t& y) const {
0549   const unsigned int nv = y.getSize();
0550 
0551 #ifdef DEBUG
0552   assert(y.zvtx_vec.size() == nv);
0553   assert(y.tvtx_vec.size() == nv);
0554   assert(y.rho_vec.size() == nv);
0555 #endif
0556 
0557   if (nv < 2)
0558     return false;
0559 
0560   bool reordering = true;
0561   bool modified = false;
0562 
0563   while (reordering) {
0564     reordering = false;
0565     for (unsigned int k = 0; (k + 1) < nv; k++) {
0566       if (y.zvtx[k + 1] < y.zvtx[k]) {
0567         auto ztemp = y.zvtx[k];
0568         y.zvtx[k] = y.zvtx[k + 1];
0569         y.zvtx[k + 1] = ztemp;
0570         auto ttemp = y.tvtx[k];
0571         y.tvtx[k] = y.tvtx[k + 1];
0572         y.tvtx[k + 1] = ttemp;
0573 #ifdef USEVTXDT2
0574         auto dt2temp = y.dt2[k];
0575         y.dt2[k] = y.dt2[k + 1];
0576         y.dt2[k + 1] = dt2temp;
0577 #endif
0578         auto ptemp = y.rho[k];
0579         y.rho[k] = y.rho[k + 1];
0580         y.rho[k + 1] = ptemp;
0581         reordering = true;
0582       }
0583     }
0584     modified |= reordering;
0585   }
0586 
0587   if (modified) {
0588     y.extractRaw();
0589     return true;
0590   }
0591 
0592   return false;
0593 }
0594 
0595 bool DAClusterizerInZT_vect::find_nearest(
0596     double z, double t, vertex_t& y, unsigned int& k_min, double dz, double dt) const {
0597   // find the cluster nearest to (z,t)
0598   // distance measure is   delta = (delta_z / dz )**2  + (delta_t/ d_t)**2
0599   // assumes that clusters are ordered n z
0600   // return value is false if no neighbour with distance < 1 is found
0601 
0602   unsigned int nv = y.getSize();
0603   if (nv < 2) {
0604     k_min = 0;
0605     return false;
0606   }
0607 
0608   // find nearest in z, binary search later
0609   unsigned int k = 0;
0610   for (unsigned int k0 = 1; k0 < nv; k0++) {
0611     if (std::abs(y.zvtx[k0] - z) < std::abs(y.zvtx[k] - z)) {
0612       k = k0;
0613     }
0614   }
0615 
0616   double delta_min = 1.;
0617 
0618   //search left
0619   unsigned int k1 = k;
0620   while ((k1 > 0) && ((y.zvtx[k] - y.zvtx[--k1]) < dz)) {
0621     auto delta = std::pow((y.zvtx[k] - y.zvtx[k1]) / dz, 2) + std::pow((y.tvtx[k] - y.tvtx[k1]) / dt, 2);
0622     if (delta < delta_min) {
0623       k_min = k1;
0624       delta_min = delta;
0625     }
0626   }
0627 
0628   //search right
0629   k1 = k;
0630   while (((++k1) < nv) && ((y.zvtx[k1] - y.zvtx[k]) < dz)) {
0631     auto delta = std::pow((y.zvtx[k1] - y.zvtx[k]) / dz, 2) + std::pow((y.tvtx[k1] - y.tvtx[k]) / dt, 2);
0632     if (delta < delta_min) {
0633       k_min = k1;
0634       delta_min = delta;
0635     }
0636   }
0637 
0638   return (delta_min < 1.);
0639 }
0640 
0641 unsigned int DAClusterizerInZT_vect::thermalize(
0642     double beta, track_t& tks, vertex_t& v, const double delta_max0, const double rho0) const {
0643   unsigned int niter = 0;
0644   double delta = 0;
0645   double delta_max = delta_lowT_;
0646 
0647   if (convergence_mode_ == 0) {
0648     delta_max = delta_max0;
0649   } else if (convergence_mode_ == 1) {
0650     delta_max = delta_lowT_ / sqrt(std::max(beta, 1.0));
0651   }
0652 
0653   zorder(v);
0654   set_vtx_range(beta, tks, v);
0655 
0656   double delta_sum_range = 0;  // accumulate max(|delta-z|) as a lower bound
0657   std::vector<double> z0 = v.zvtx_vec;
0658 
0659   while (niter++ < maxIterations_) {
0660     delta = update(beta, tks, v, rho0);
0661     delta_sum_range += delta;
0662 
0663     if (delta_sum_range > zrange_min_) {
0664       for (unsigned int k = 0; k < v.getSize(); k++) {
0665         if (std::abs(v.zvtx_vec[k] - z0[k]) > zrange_min_) {
0666           zorder(v);
0667           set_vtx_range(beta, tks, v);
0668           delta_sum_range = 0;
0669           z0 = v.zvtx_vec;
0670           break;
0671         }
0672       }
0673     }
0674 
0675     if (delta < delta_max) {
0676       break;
0677     }
0678   }
0679 
0680 #ifdef DEBUG
0681   if (DEBUGLEVEL > 0) {
0682     std::cout << "DAClusterizerInZT_vect.thermalize niter = " << niter << " at T = " << 1 / beta
0683               << "  nv = " << v.getSize() << std::endl;
0684     if (DEBUGLEVEL > 2)
0685       dump(beta, v, tks, 0);
0686   }
0687 #endif
0688 
0689   return niter;
0690 }
0691 
0692 bool DAClusterizerInZT_vect::merge(vertex_t& y, track_t& tks, double& beta) const {
0693   // merge clusters that collapsed or never separated,
0694   // return true if vertices were merged, false otherwise
0695   const unsigned int nv = y.getSize();
0696 
0697   if (nv < 2)
0698     return false;
0699 
0700   // merge the smallest distance clusters first
0701   unsigned int k1_min = 0, k2_min = 0;
0702   double delta_min = 0;
0703 
0704   for (unsigned int k1 = 0; (k1 + 1) < nv; k1++) {
0705     unsigned int k2 = k1;
0706     while ((++k2 < nv) && (std::fabs(y.zvtx[k2] - y.zvtx[k1]) < zmerge_)) {
0707       auto delta = std::pow((y.zvtx[k2] - y.zvtx[k1]) / zmerge_, 2) + std::pow((y.tvtx[k2] - y.tvtx[k1]) / tmerge_, 2);
0708       if ((delta < delta_min) || (k1_min == k2_min)) {
0709         k1_min = k1;
0710         k2_min = k2;
0711         delta_min = delta;
0712       }
0713     }
0714   }
0715 
0716   if ((k1_min == k2_min) || (delta_min > 1)) {
0717     return false;
0718   }
0719 
0720   double rho = y.rho[k1_min] + y.rho[k2_min];
0721 
0722 #ifdef DEBUG
0723   assert((k1_min < nv) && (k2_min < nv));
0724   if (DEBUGLEVEL > 1) {
0725     std::cout << "merging (" << setw(8) << fixed << setprecision(4) << y.zvtx[k1_min] << ',' << y.tvtx[k1_min]
0726               << ") and (" << y.zvtx[k2_min] << ',' << y.tvtx[k2_min] << ")"
0727               << "  idx=" << k1_min << "," << k2_min << std::endl;
0728   }
0729 #endif
0730 
0731   if (rho > 0) {
0732     y.zvtx[k1_min] = (y.rho[k1_min] * y.zvtx[k1_min] + y.rho[k2_min] * y.zvtx[k2_min]) / rho;
0733     y.tvtx[k1_min] = (y.rho[k1_min] * y.tvtx[k1_min] + y.rho[k2_min] * y.tvtx[k2_min]) / rho;
0734 #ifdef USEVTXDT2
0735     y.dt2[k1_min] = (y.rho[k1_min] * y.dt2[k1_min] + y.rho[k2_min] * y.dt2[k2_min]) / rho;
0736 #endif
0737   } else {
0738     y.zvtx[k1_min] = 0.5 * (y.zvtx[k1_min] + y.zvtx[k2_min]);
0739     y.tvtx[k1_min] = 0.5 * (y.tvtx[k1_min] + y.tvtx[k2_min]);
0740   }
0741   y.rho[k1_min] = rho;
0742   y.removeItem(k2_min, tks);
0743 
0744   zorder(y);
0745   set_vtx_range(beta, tks, y);
0746   y.extractRaw();
0747   return true;
0748 }
0749 
0750 bool DAClusterizerInZT_vect::purge(vertex_t& y, track_t& tks, double& rho0, const double beta) const {
0751   constexpr double eps = 1.e-100;
0752   // eliminate clusters with only one significant/unique track
0753   const unsigned int nv = y.getSize();
0754   const unsigned int nt = tks.getSize();
0755 
0756   if (nv < 2)
0757     return false;
0758 
0759   std::vector<double> sump_v(nv), arg_cache_v(nv), exp_cache_v(nv), pcut_cache_v(nv);
0760   std::vector<int> nUnique_v(nv);
0761   double* __restrict__ parg_cache;
0762   double* __restrict__ pexp_cache;
0763   double* __restrict__ ppcut_cache;
0764   double* __restrict__ psump;
0765   int* __restrict__ pnUnique;
0766   int constexpr nunique_min_ = 2;
0767 
0768   zorder(y);
0769   set_vtx_range(beta, tks, y);
0770 
0771   parg_cache = arg_cache_v.data();
0772   pexp_cache = exp_cache_v.data();
0773   ppcut_cache = pcut_cache_v.data();
0774   psump = sump_v.data();
0775   pnUnique = nUnique_v.data();
0776 
0777   const auto rhoconst = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
0778   for (unsigned int k = 0; k < nv; ++k) {
0779     const double pmax = y.rho[k] / (y.rho[k] + rhoconst);
0780     ppcut_cache[k] = uniquetrkweight_ * pmax;
0781   }
0782 
0783   for (unsigned int i = 0; i < nt; i++) {
0784     const auto invZ = ((tks.sum_Z[i] > eps) && (tks.tkwt[i] > uniquetrkminp_)) ? 1. / tks.sum_Z[i] : 0.;
0785     const auto track_z = tks.zpca[i];
0786     const auto track_t = tks.tpca[i];
0787     const auto botrack_dz2 = -beta * tks.dz2[i];
0788     const auto botrack_dt2 = -beta * tks.dt2[i];
0789     const auto kmin = tks.kmin[i];
0790     const auto kmax = tks.kmax[i];
0791 
0792     for (unsigned int k = kmin; k < kmax; k++) {
0793       const auto mult_resz = track_z - y.zvtx[k];
0794       const auto mult_rest = track_t - y.tvtx[k];
0795       parg_cache[k] = botrack_dz2 * (mult_resz * mult_resz) + botrack_dt2 * (mult_rest * mult_rest);
0796     }
0797 
0798     local_exp_list_range(parg_cache, pexp_cache, kmin, kmax);
0799 
0800     for (unsigned int k = kmin; k < kmax; k++) {
0801       const double p = y.rho[k] * pexp_cache[k] * invZ;
0802       psump[k] += p;
0803       pnUnique[k] += (p > ppcut_cache[k]) ? 1 : 0;
0804     }
0805   }
0806 
0807   double sumpmin = nt;
0808   unsigned int k0 = nv;
0809   for (unsigned k = 0; k < nv; k++) {
0810     if ((pnUnique[k] < nunique_min_) && (psump[k] < sumpmin)) {
0811       sumpmin = psump[k];
0812       k0 = k;
0813     }
0814   }
0815 
0816   if (k0 != nv) {
0817 #ifdef DEBUG
0818     assert(k0 < y.getSize());
0819     if (DEBUGLEVEL > 1) {
0820       std::cout << "eliminating prototype at " << std::setw(10) << std::setprecision(4) << y.zvtx[k0] << ","
0821                 << y.tvtx[k0] << " with sump=" << sumpmin << "  rho*nt =" << y.rho[k0] * nt << endl;
0822     }
0823 #endif
0824     y.removeItem(k0, tks);
0825     set_vtx_range(beta, tks, y);
0826     return true;
0827   } else {
0828     return false;
0829   }
0830 }
0831 
0832 double DAClusterizerInZT_vect::beta0(double betamax, track_t const& tks, vertex_t const& y) const {
0833   double T0 = 0;  // max Tc for beta=0
0834   // estimate critical temperature from beta=0 (T=inf)
0835   const unsigned int nt = tks.getSize();
0836   const unsigned int nv = y.getSize();
0837 
0838   for (unsigned int k = 0; k < nv; k++) {
0839     // vertex fit at T=inf
0840     double sumwz = 0;
0841     double sumwt = 0;
0842     double sumw_z = 0;
0843     double sumw_t = 0;
0844     for (unsigned int i = 0; i < nt; i++) {
0845       double w_z = tks.tkwt[i] * tks.dz2[i];
0846       double w_t = tks.tkwt[i] * tks.dt2[i];
0847       sumwz += w_z * tks.zpca[i];
0848       sumwt += w_t * tks.tpca[i];
0849       sumw_z += w_z;
0850       sumw_t += w_t;
0851     }
0852     y.zvtx[k] = sumwz / sumw_z;
0853     y.tvtx[k] = sumwt / sumw_t;
0854 
0855     // estimate Tc, eventually do this in the same loop
0856     double szz = 0, stt = 0, szt = 0;
0857     double nuz = 0, nut = 0;
0858     for (unsigned int i = 0; i < nt; i++) {
0859       double dz = (tks.zpca[i] - y.zvtx[k]) * tks.dz2[i];
0860       double dt = (tks.tpca[i] - y.tvtx[k]) * tks.dt2[i];
0861       double w = tks.tkwt[i];
0862       szz += w * dz * dz;
0863       stt += w * dt * dt;
0864       szt += w * dz * dt;
0865       nuz += w * tks.dz2[i];
0866       nut += w * tks.dt2[i];
0867     }
0868     double Tz = szz / nuz;
0869     double Tt = 0;
0870     double Tc = 0;
0871     if (nut > 0) {
0872       Tt = stt / nut;
0873       Tc = Tz + Tt + sqrt(pow(Tz - Tt, 2) + 4 * szt * szt / nuz / nut);
0874     } else {
0875       Tc = 2. * Tz;
0876     }
0877 
0878     if (Tc > T0)
0879       T0 = Tc;
0880   }  // vertex loop (normally there should be only one vertex at beta=0)
0881 
0882 #ifdef DEBUG
0883   if (DEBUGLEVEL > 0) {
0884     std::cout << "DAClusterizerInZT_vect.beta0:   Tc = " << T0 << std::endl;
0885     int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
0886     std::cout << "DAClusterizerInZT_vect.beta0:   nstep = " << coolingsteps << std::endl;
0887   }
0888 #endif
0889 
0890   if (T0 > 1. / betamax) {
0891     int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
0892     return betamax * std::pow(coolingFactor_, coolingsteps);
0893   } else {
0894     // ensure at least one annealing step
0895     return betamax * coolingFactor_;
0896   }
0897 }
0898 
0899 double DAClusterizerInZT_vect::get_Tc(const vertex_t& y, int k) const {
0900   double Tz = y.szz[k] / y.nuz[k];  // actually 0.5*Tc(z)
0901   double Tt = 0.;
0902   if (y.nut[k] > 0) {
0903     Tt = y.stt[k] / y.nut[k];
0904     double mx = y.szt[k] / y.nuz[k] * y.szt[k] / y.nut[k];
0905     return Tz + Tt + sqrt(pow(Tz - Tt, 2) + 4 * mx);
0906   }
0907   return 2 * Tz;
0908 }
0909 
0910 bool DAClusterizerInZT_vect::split(const double beta, track_t& tks, vertex_t& y, double threshold) const {
0911   // split only critical vertices (Tc >~ T=1/beta   <==>   beta*Tc>~1)
0912   // an update must have been made just before doing this (same beta, no merging)
0913   // returns true if at least one cluster was split
0914 
0915   constexpr double epsilonz = 1e-3;  // minimum split size z
0916   constexpr double epsilont = 1e-2;  // minimum split size t
0917   unsigned int nv = y.getSize();
0918   const double twoBeta = 2.0 * beta;
0919 
0920   // avoid left-right biases by splitting highest Tc first
0921 
0922   std::vector<std::pair<double, unsigned int> > critical;
0923   for (unsigned int k = 0; k < nv; k++) {
0924     double Tc = get_Tc(y, k);
0925     if (beta * Tc > threshold) {
0926       critical.push_back(make_pair(Tc, k));
0927     }
0928   }
0929   if (critical.empty())
0930     return false;
0931 
0932   std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int> >());
0933 
0934   bool split = false;
0935   const unsigned int nt = tks.getSize();
0936 
0937   for (unsigned int ic = 0; ic < critical.size(); ic++) {
0938     unsigned int k = critical[ic].second;
0939 
0940     // split direction in the (z,t)-plane
0941 
0942     double Mzz = y.nuz[k] - twoBeta * y.szz[k];
0943     double Mtt = y.nut[k] - twoBeta * y.stt[k];
0944     double Mzt = -twoBeta * y.szt[k];
0945     const double twoMzt = 2.0 * Mzt;
0946     double D = sqrt(pow(Mtt - Mzz, 2) + twoMzt * twoMzt);
0947     double q1 = atan2(-Mtt + Mzz + D, -twoMzt);
0948     double l1 = 0.5 * (-Mzz - Mtt + D);
0949     double l2 = 0.5 * (-Mzz - Mtt - D);
0950     if ((std::abs(l1) < 1e-4) && (std::abs(l2) < 1e-4)) {
0951       edm::LogWarning("DAClusterizerInZT_vect") << "warning, bad eigenvalues!  idx=" << k << " z= " << y.zvtx[k]
0952                                                 << " Mzz=" << Mzz << "   Mtt=" << Mtt << "  Mzt=" << Mzt << endl;
0953     }
0954 
0955     double qsplit = q1;
0956     double cq = cos(qsplit);
0957     double sq = sin(qsplit);
0958     if (cq < 0) {
0959       cq = -cq;
0960       sq = -sq;
0961     }  // prefer cq>0 to keep z-ordering
0962 
0963     // estimate subcluster positions and weight
0964     double p1 = 0, z1 = 0, t1 = 0, wz1 = 0, wt1 = 0;
0965     double p2 = 0, z2 = 0, t2 = 0, wz2 = 0, wt2 = 0;
0966     for (unsigned int i = 0; i < nt; ++i) {
0967       if (tks.sum_Z[i] > 1.e-100) {
0968         double lr = (tks.zpca[i] - y.zvtx[k]) * cq + (tks.tpca[i] - y.tvtx[k]) * sq;
0969         // winner-takes-all, usually overestimates splitting
0970         double tl = lr < 0 ? 1. : 0.;
0971         double tr = 1. - tl;
0972 
0973         // soften it, especially at low T
0974         double arg = lr * std::sqrt(beta * (cq * cq * tks.dz2[i] + sq * sq * tks.dt2[i]));
0975         if (std::abs(arg) < 20) {
0976           double t = local_exp(-arg);
0977           tl = t / (t + 1.);
0978           tr = 1 / (t + 1.);
0979         }
0980 
0981         double p = y.rho[k] * tks.tkwt[i] *
0982                    local_exp(-beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i], tks.tpca[i], y.tvtx[k], tks.dt2[i])) /
0983                    tks.sum_Z[i];
0984         double wz = p * tks.dz2[i];
0985         double wt = p * tks.dt2[i];
0986         p1 += p * tl;
0987         z1 += wz * tl * tks.zpca[i];
0988         t1 += wt * tl * tks.tpca[i];
0989         wz1 += wz * tl;
0990         wt1 += wt * tl;
0991         p2 += p * tr;
0992         z2 += wz * tr * tks.zpca[i];
0993         t2 += wt * tr * tks.tpca[i];
0994         wz2 += wz * tr;
0995         wt2 += wt * tr;
0996       }
0997     }
0998 
0999     if (wz1 > 0) {
1000       z1 /= wz1;
1001     } else {
1002       z1 = y.zvtx[k] - epsilonz * cq;
1003       edm::LogWarning("DAClusterizerInZT_vect") << "warning, wz1 = " << scientific << wz1 << endl;
1004     }
1005     if (wt1 > 0) {
1006       t1 /= wt1;
1007     } else {
1008       t1 = y.tvtx[k] - epsilont * sq;
1009       edm::LogWarning("DAClusterizerInZT_vect") << "warning, wt1 = " << scientific << wt1 << endl;
1010     }
1011     if (wz2 > 0) {
1012       z2 /= wz2;
1013     } else {
1014       z2 = y.zvtx[k] + epsilonz * cq;
1015       edm::LogWarning("DAClusterizerInZT_vect") << "warning, wz2 = " << scientific << wz2 << endl;
1016     }
1017     if (wt2 > 0) {
1018       t2 /= wt2;
1019     } else {
1020       t2 = y.tvtx[k] + epsilont * sq;
1021       edm::LogWarning("DAClusterizerInZT_vect") << "warning, wt2 = " << scientific << wt2 << endl;
1022     }
1023 
1024     unsigned int k_min1 = k, k_min2 = k;
1025     constexpr double spliteps = 1e-8;
1026     while (((find_nearest(z1, t1, y, k_min1, epsilonz, epsilont) && (k_min1 != k)) ||
1027             (find_nearest(z2, t2, y, k_min2, epsilonz, epsilont) && (k_min2 != k))) &&
1028            (std::abs(z2 - z1) > spliteps || std::abs(t2 - t1) > spliteps)) {
1029       z1 = 0.5 * (z1 + y.zvtx[k]);
1030       t1 = 0.5 * (t1 + y.tvtx[k]);
1031       z2 = 0.5 * (z2 + y.zvtx[k]);
1032       t2 = 0.5 * (t2 + y.tvtx[k]);
1033     }
1034 
1035 #ifdef DEBUG
1036     assert(k < nv);
1037     if (DEBUGLEVEL > 1) {
1038       if (std::fabs(y.zvtx[k] - zdumpcenter_) < zdumpwidth_) {
1039         std::cout << " T= " << std::setw(10) << std::setprecision(1) << 1. / beta << " Tc= " << critical[ic].first
1040                   << " direction =" << std::setprecision(4) << qsplit << "    splitting (" << std::setw(8) << std::fixed
1041                   << std::setprecision(4) << y.zvtx[k] << "," << y.tvtx[k] << ")"
1042                   << " --> (" << z1 << ',' << t1 << "),(" << z2 << ',' << t2 << ")     [" << p1 << "," << p2 << "]";
1043         if (std::fabs(z2 - z1) > epsilonz || std::fabs(t2 - t1) > epsilont) {
1044           std::cout << std::endl;
1045         } else {
1046           std::cout << "  rejected " << std::endl;
1047         }
1048       }
1049     }
1050 #endif
1051 
1052     if (z1 > z2) {
1053       edm::LogInfo("DAClusterizerInZT") << "warning   swapping z in split  qsplit=" << qsplit << "   cq=" << cq
1054                                         << "  sq=" << sq << endl;
1055       auto ztemp = z1;
1056       auto ttemp = t1;
1057       auto ptemp = p1;
1058       z1 = z2;
1059       t1 = t2;
1060       p1 = p2;
1061       z2 = ztemp;
1062       t2 = ttemp;
1063       p2 = ptemp;
1064     }
1065 
1066     // split if the new subclusters are significantly separated
1067     if (std::fabs(z2 - z1) > epsilonz || std::fabs(t2 - t1) > epsilont) {
1068       split = true;
1069       double rho1 = p1 * y.rho[k] / (p1 + p2);
1070       double rho2 = p2 * y.rho[k] / (p1 + p2);
1071 
1072       // replace the original by (z2,t2)
1073       y.removeItem(k, tks);
1074       unsigned int k2 = y.insertOrdered(z2, t2, rho2, tks);
1075 
1076 #ifdef DEBUG
1077       if (k2 < k) {
1078         std::cout << "unexpected z-ordering in split" << std::endl;
1079       }
1080 #endif
1081 
1082       // adjust pointers if necessary
1083       if (!(k == k2)) {
1084         for (unsigned int jc = ic; jc < critical.size(); jc++) {
1085           if (critical[jc].second > k) {
1086             critical[jc].second--;
1087           }
1088           if (critical[jc].second >= k2) {
1089             critical[jc].second++;
1090           }
1091         }
1092       }
1093 
1094       // insert (z1,t1) where it belongs
1095       unsigned int k1 = y.insertOrdered(z1, t1, rho1, tks);
1096       nv++;
1097 
1098       // adjust remaining pointers
1099       for (unsigned int jc = ic; jc < critical.size(); jc++) {
1100         if (critical[jc].second >= k1) {
1101           critical[jc].second++;
1102         }
1103       }
1104     } else {
1105 #ifdef DEBUG
1106       std::cout << "warning ! split rejected, too small." << endl;
1107 #endif
1108     }
1109   }
1110 
1111   return split;
1112 }
1113 
1114 vector<TransientVertex> DAClusterizerInZT_vect::vertices(const vector<reco::TransientTrack>& tracks) const {
1115   track_t&& tks = fill(tracks);
1116   vector<TransientVertex> clusters;
1117   if (tks.getSize() == 0)
1118     return clusters;
1119   tks.extractRaw();
1120 
1121   unsigned int nt = tks.getSize();
1122   double rho0 = 0.0;  // start with no outlier rejection
1123 
1124   vertex_t y;  // the vertex prototypes
1125 
1126   // initialize:single vertex at infinite temperature
1127   y.addItem(0, 0, 1.0);
1128   clear_vtx_range(tks, y);
1129 
1130   // estimate first critical temperature
1131   double beta = beta0(betamax_, tks, y);
1132 #ifdef DEBUG
1133   if (DEBUGLEVEL > 0)
1134     std::cout << "Beta0 is " << beta << std::endl;
1135 #endif
1136 
1137   thermalize(beta, tks, y, delta_highT_, 0.);
1138 
1139   // annealing loop, stop when T<minT  (i.e. beta>1/minT)
1140 
1141   double betafreeze = betamax_ * sqrt(coolingFactor_);
1142 
1143   while (beta < betafreeze) {
1144     update(beta, tks, y, rho0, true);
1145     while (merge(y, tks, beta)) {
1146       if (zorder(y))
1147         set_vtx_range(beta, tks, y);
1148       update(beta, tks, y, rho0, true);
1149     }
1150     split(beta, tks, y);
1151 
1152     beta = beta / coolingFactor_;
1153     thermalize(beta, tks, y, delta_highT_, 0.);
1154   }
1155 
1156 #ifdef DEBUG
1157   verify(y, tks);
1158 
1159   if (DEBUGLEVEL > 0) {
1160     std::cout << "DAClusterizerInZT_vect::vertices :"
1161               << "merging at T=" << 1 / beta << std::endl;
1162   }
1163 #endif
1164 
1165   zorder(y);
1166   set_vtx_range(beta, tks, y);
1167   // ZT_vect version of merge() does not use Tc, but split() does
1168   update(beta, tks, y, rho0, true);
1169 
1170   while (merge(y, tks, beta)) {
1171     zorder(y);
1172     set_vtx_range(beta, tks, y);
1173     update(beta, tks, y, rho0, true);
1174   }
1175 
1176 #ifdef DEBUG
1177   verify(y, tks);
1178   if (DEBUGLEVEL > 0) {
1179     std::cout << "DAClusterizerInZT_vect::vertices :"
1180               << "splitting/merging at T=" << 1 / beta << std::endl;
1181   }
1182 #endif
1183 
1184   unsigned int ntry = 0;
1185   double threshold = 1.0;
1186   while (split(beta, tks, y, threshold) && (ntry++ < 10)) {
1187     thermalize(beta, tks, y, delta_highT_, 0.);
1188 
1189     while (merge(y, tks, beta)) {
1190       if (zorder(y))
1191         set_vtx_range(beta, tks, y);
1192       update(beta, tks, y, rho0, true);  // updateTc in case we go back to split
1193     }
1194 
1195 #ifdef DEBUG
1196     verify(y, tks);
1197 #endif
1198 
1199     // relax splitting a bit to reduce multiple split-merge cycles of the same cluster
1200     threshold *= 1.1;
1201   }
1202 
1203 #ifdef DEBUG
1204   verify(y, tks);
1205   if (DEBUGLEVEL > 0) {
1206     std::cout << "DAClusterizerInZT_vect::vertices :"
1207               << "turning on outlier rejection at T=" << 1 / beta << std::endl;
1208   }
1209 #endif
1210 
1211   // switch on outlier rejection at T=minT
1212   if (dzCutOff_ > 0) {
1213     rho0 = 1. / nt;
1214     for (unsigned int a = 0; a < 5; a++) {
1215       update(beta, tks, y, a * rho0 / 5);  // adiabatic turn-on
1216     }
1217   }
1218 
1219   thermalize(beta, tks, y, delta_lowT_, rho0);
1220 
1221 #ifdef DEBUG
1222   verify(y, tks);
1223   if (DEBUGLEVEL > 0) {
1224     std::cout << "DAClusterizerInZT_vect::vertices :"
1225               << "merging with outlier rejection at T=" << 1 / beta << std::endl;
1226   }
1227   if (DEBUGLEVEL > 2)
1228     dump(beta, y, tks, 2);
1229 #endif
1230 
1231   // merge again  (some cluster split by outliers collapse here)
1232   zorder(y);
1233   while (merge(y, tks, beta)) {
1234     set_vtx_range(beta, tks, y);
1235     update(beta, tks, y, rho0);
1236   }
1237 
1238 #ifdef DEBUG
1239   if (DEBUGLEVEL > 0) {
1240     std::cout << "DAClusterizerInZT_vect::vertices :"
1241               << "after merging with outlier rejection at T=" << 1 / beta << std::endl;
1242   }
1243   if (DEBUGLEVEL > 1)
1244     dump(beta, y, tks, 2);
1245 #endif
1246 
1247   // go down to the purging temperature (if it is lower than tmin)
1248   while (beta < betapurge_) {
1249     beta = min(beta / coolingFactor_, betapurge_);
1250     thermalize(beta, tks, y, delta_lowT_, rho0);
1251   }
1252 
1253 #ifdef DEBUG
1254   verify(y, tks);
1255   if (DEBUGLEVEL > 0) {
1256     std::cout << "DAClusterizerInZT_vect::vertices :"
1257               << "purging at T=" << 1 / beta << std::endl;
1258   }
1259 #endif
1260 
1261   // eliminate insignificant vertices, this is more restrictive at higher T
1262   while (purge(y, tks, rho0, beta)) {
1263     thermalize(beta, tks, y, delta_lowT_, rho0);
1264     if (zorder(y)) {
1265       set_vtx_range(beta, tks, y);
1266     };
1267   }
1268 
1269 #ifdef DEBUG
1270   verify(y, tks);
1271   if (DEBUGLEVEL > 0) {
1272     std::cout << "DAClusterizerInZT_vect::vertices :"
1273               << "last cooling T=" << 1 / beta << std::endl;
1274   }
1275 #endif
1276 
1277   // optionally cool some more without doing anything, to make the assignment harder
1278   while (beta < betastop_) {
1279     beta = min(beta / coolingFactor_, betastop_);
1280     thermalize(beta, tks, y, delta_lowT_, rho0);
1281     if (zorder(y)) {
1282       set_vtx_range(beta, tks, y);
1283     };
1284   }
1285 
1286 #ifdef DEBUG
1287   verify(y, tks);
1288   if (DEBUGLEVEL > 0) {
1289     std::cout << "DAClusterizerInZT_vect::vertices :"
1290               << "stop cooling at T=" << 1 / beta << std::endl;
1291   }
1292   if (DEBUGLEVEL > 2)
1293     dump(beta, y, tks, 2);
1294 #endif
1295 
1296   // new, merge here and not in "clusterize"
1297   // final merging step
1298   double betadummy = 1;
1299   while (merge(y, tks, betadummy))
1300     ;
1301 
1302   // select significant tracks and use a TransientVertex as a container
1303   zorder(y);
1304   set_vtx_range(beta, tks, y);
1305   const unsigned int nv = y.getSize();
1306   for (unsigned int k = 0; k < nv; k++)
1307     if (edm::isNotFinite(y.rho[k]) || edm::isNotFinite(y.zvtx[k])) {
1308       y.rho[k] = 0;
1309       y.zvtx[k] = 0;
1310     }
1311 
1312   const auto z_sum_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
1313   std::vector<std::vector<unsigned int> > vtx_track_indices(nv);
1314   for (unsigned int i = 0; i < nt; i++) {
1315     const auto kmin = tks.kmin[i];
1316     const auto kmax = tks.kmax[i];
1317     for (auto k = kmin; k < kmax; k++) {
1318       y.exp_arg[k] = -beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i], tks.tpca[i], y.tvtx[k], tks.dt2[i]);
1319     }
1320 
1321     local_exp_list_range(y.exp_arg, y.exp, kmin, kmax);
1322 
1323     tks.sum_Z[i] = z_sum_init;
1324     for (auto k = kmin; k < kmax; k++) {
1325       tks.sum_Z[i] += y.rho[k] * y.exp[k];
1326     }
1327     const double invZ = tks.sum_Z[i] > 1e-100 ? 1. / tks.sum_Z[i] : 0.0;
1328 
1329     for (auto k = kmin; k < kmax; k++) {
1330       double p = y.rho[k] * y.exp[k] * invZ;
1331       if (p > mintrkweight_) {
1332         // assign  track i -> vertex k (hard, mintrkweight_ should be >= 0.5 here
1333         vtx_track_indices[k].push_back(i);
1334         break;
1335       }
1336     }
1337 
1338   }  // track loop
1339 
1340   GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01);
1341   for (unsigned int k = 0; k < nv; k++) {
1342     if (!vtx_track_indices[k].empty()) {
1343       GlobalPoint pos(0, 0, y.zvtx[k]);
1344       vector<reco::TransientTrack> vertexTracks;
1345       for (auto i : vtx_track_indices[k]) {
1346         vertexTracks.push_back(*(tks.tt[i]));
1347       }
1348       TransientVertex v(pos, dummyError, vertexTracks, 0);
1349       clusters.push_back(v);
1350     }
1351   }
1352 
1353   return clusters;
1354 }
1355 
1356 vector<vector<reco::TransientTrack> > DAClusterizerInZT_vect::clusterize(
1357     const vector<reco::TransientTrack>& tracks) const {
1358   vector<vector<reco::TransientTrack> > clusters;
1359   vector<TransientVertex>&& pv = vertices(tracks);
1360 
1361 #ifdef DEBUG
1362   if (DEBUGLEVEL > 0) {
1363     std::cout << "###################################################" << endl;
1364     std::cout << "# vectorized DAClusterizerInZT_vect::clusterize   nt=" << tracks.size() << endl;
1365     std::cout << "# DAClusterizerInZT_vect::clusterize   pv.size=" << pv.size() << endl;
1366     std::cout << "###################################################" << endl;
1367   }
1368 #endif
1369 
1370   if (pv.empty()) {
1371     return clusters;
1372   }
1373 
1374   // fill into clusters, don't merge
1375   for (auto k = pv.begin(); k != pv.end(); k++) {
1376     vector<reco::TransientTrack> aCluster = k->originalTracks();
1377     if (aCluster.size() > 1) {
1378       clusters.push_back(aCluster);
1379     }
1380   }
1381 
1382   return clusters;
1383 }
1384 
1385 void DAClusterizerInZT_vect::dump(const double beta, const vertex_t& y, const track_t& tks, int verbosity) const {
1386 #ifdef DEBUG
1387   const unsigned int nv = y.getSize();
1388   const unsigned int nt = tks.getSize();
1389 
1390   std::vector<unsigned int> iz;
1391   for (unsigned int j = 0; j < nt; j++) {
1392     iz.push_back(j);
1393   }
1394   std::sort(iz.begin(), iz.end(), [tks](unsigned int a, unsigned int b) { return tks.zpca[a] < tks.zpca[b]; });
1395   std::cout << std::endl;
1396   std::cout << "-----DAClusterizerInZT::dump ----" << nv << "  clusters " << std::endl;
1397   string h = "                                                                                 ";
1398   std::cout << h << " k= ";
1399   for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1400     if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1401       std::cout << setw(8) << fixed << ivertex;
1402     }
1403   }
1404   std::cout << endl;
1405 
1406   std::cout << h << " z= ";
1407   std::cout << setprecision(4);
1408   for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1409     if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1410       std::cout << setw(8) << fixed << y.zvtx[ivertex];
1411     }
1412   }
1413   std::cout << endl;
1414 
1415   std::cout << h << " t= ";
1416   std::cout << setprecision(4);
1417   for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1418     if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1419       std::cout << setw(8) << fixed << y.tvtx[ivertex];
1420     }
1421   }
1422   std::cout << endl;
1423 
1424   std::cout << "T=" << setw(15) << 1. / beta << " Tmin =" << setw(10) << 1. / betamax_
1425             << "                                               Tc= ";
1426   for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1427     if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1428       double Tc = get_Tc(y, ivertex);
1429       std::cout << setw(8) << fixed << setprecision(1) << Tc;
1430     }
1431   }
1432   std::cout << endl;
1433 
1434   std::cout << h << "rho= ";
1435   double sumrho = 0;
1436   for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1437     sumrho += y.rho[ivertex];
1438     if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1439       continue;
1440     std::cout << setw(8) << setprecision(4) << fixed << y.rho[ivertex];
1441   }
1442   std::cout << endl;
1443 
1444   std::cout << h << "nt= ";
1445   for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1446     sumrho += y.rho[ivertex];
1447     if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1448       continue;
1449     std::cout << setw(8) << setprecision(1) << fixed << y.rho[ivertex] * nt;
1450   }
1451   std::cout << endl;
1452 
1453   if (verbosity > 0) {
1454     double E = 0, F = 0;
1455     std::cout << endl;
1456     std::cout << "----        z +/- dz          t +/- dt                ip +/-dip       pt    phi  eta   weights  ----"
1457               << endl;
1458     std::cout << setprecision(4);
1459     for (unsigned int i0 = 0; i0 < nt; i0++) {
1460       unsigned int i = iz[i0];
1461       if (tks.sum_Z[i] > 0) {
1462         F -= std::log(tks.sum_Z[i]) / beta;
1463       }
1464       double tz = tks.zpca[i];
1465 
1466       if (std::fabs(tz - zdumpcenter_) > zdumpwidth_)
1467         continue;
1468       std::cout << setw(4) << i << ")" << setw(8) << fixed << setprecision(4) << tz << " +/-" << setw(6)
1469                 << sqrt(1. / tks.dz2[i]);
1470 
1471       if (tks.dt2[i] > 0) {
1472         std::cout << setw(8) << fixed << setprecision(4) << tks.tpca[i] << " +/-" << setw(6) << sqrt(1. / tks.dt2[i]);
1473       } else {
1474         std::cout << "                  ";
1475       }
1476 
1477       if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) {
1478         std::cout << " *";
1479       } else {
1480         std::cout << "  ";
1481       }
1482       if (tks.tt[i]->track().hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) {
1483         std::cout << "+";
1484       } else {
1485         std::cout << "-";
1486       }
1487       std::cout << setw(1)
1488                 << tks.tt[i]
1489                        ->track()
1490                        .hitPattern()
1491                        .pixelBarrelLayersWithMeasurement();  // see DataFormats/TrackReco/interface/HitPattern.h
1492       std::cout << setw(1) << tks.tt[i]->track().hitPattern().pixelEndcapLayersWithMeasurement();
1493       std::cout << setw(1) << hex
1494                 << tks.tt[i]->track().hitPattern().trackerLayersWithMeasurement() -
1495                        tks.tt[i]->track().hitPattern().pixelLayersWithMeasurement()
1496                 << dec;
1497       std::cout << "=" << setw(1) << hex
1498                 << tks.tt[i]->track().hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << dec;
1499 
1500       Measurement1D IP = tks.tt[i]->stateAtBeamLine().transverseImpactParameter();
1501       std::cout << setw(8) << IP.value() << "+/-" << setw(6) << IP.error();
1502       std::cout << " " << setw(6) << setprecision(2) << tks.tt[i]->track().pt() * tks.tt[i]->track().charge();
1503       std::cout << " " << setw(5) << setprecision(2) << tks.tt[i]->track().phi() << " " << setw(5) << setprecision(2)
1504                 << tks.tt[i]->track().eta();
1505 
1506       double sump = 0.;
1507       for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1508         if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1509           continue;
1510 
1511         if ((tks.tkwt[i] > 0) && (tks.sum_Z[i] > 0)) {
1512           double p =
1513               y.rho[ivertex] *
1514               local_exp(-beta *
1515                         Eik(tks.zpca[i], y.zvtx[ivertex], tks.dz2[i], tks.tpca[i], y.tvtx[ivertex], tks.dt2[i])) /
1516               tks.sum_Z[i];
1517 
1518           if ((ivertex >= tks.kmin[i]) && (ivertex < tks.kmax[i])) {
1519             if (p > 0.0001) {
1520               std::cout << setw(8) << setprecision(3) << p;
1521             } else {
1522               std::cout << "    _   ";  // tiny but in the cluster list
1523             }
1524             E += p * Eik(tks.zpca[i], y.zvtx[ivertex], tks.dz2[i], tks.tpca[i], y.tvtx[ivertex], tks.dt2[i]);
1525             sump += p;
1526           } else {
1527             if (p > 0.1) {
1528               std::cout << "XXXXXXXX";  // we have an inconsistency here
1529             } else if (p > 0.0001) {
1530               std::cout << "X" << setw(6) << setprecision(3) << p << "X";
1531             } else {
1532               std::cout << "    .   ";  // not in the cluster list
1533             }
1534           }
1535         } else {
1536           std::cout << "        ";
1537         }
1538       }
1539       std::cout << "  ( " << std::setw(3) << tks.kmin[i] << "," << std::setw(3) << tks.kmax[i] - 1 << " ) ";
1540       std::cout << endl;
1541     }
1542     std::cout << "                                                                                    ";
1543     for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1544       if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1545         std::cout << "   " << setw(3) << ivertex << "  ";
1546       }
1547     }
1548     std::cout << endl;
1549     std::cout << "                                                                                 z= ";
1550     std::cout << setprecision(4);
1551     for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1552       if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1553         std::cout << setw(8) << fixed << y.zvtx[ivertex];
1554       }
1555     }
1556     std::cout << endl;
1557     std::cout << endl
1558               << "T=" << 1 / beta << " E=" << E << " n=" << y.getSize() << "  F= " << F << endl
1559               << "----------" << endl;
1560   }
1561 #endif
1562 }
1563 
1564 void DAClusterizerInZT_vect::fillPSetDescription(edm::ParameterSetDescription& desc) {
1565   DAClusterizerInZ_vect::fillPSetDescription(desc);
1566   desc.add<double>("tmerge", 0.1);            // 4D only
1567   desc.add<double>("dtCutOff", 4.);           // 4D only
1568   desc.add<double>("t0Max", 1.0);             // 4D only
1569   desc.add<double>("vertexSizeTime", 0.008);  // 4D only
1570 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!