Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-03-12 06:22:00 

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