CMS 3D CMS Logo
List of all members | Classes | Public Member Functions | Private Attributes
DAClusterizerInZT_vect Class Referencefinal

#include <DAClusterizerInZT_vect.h>

Inheritance diagram for DAClusterizerInZT_vect:
TrackClusterizerInZ

Classes

struct  track_t
 
struct  vertex_t
 

Public Member Functions

double beta0 (const double betamax, track_t const &tks, vertex_t const &y) const
 
void clear_vtx_range (track_t &gtracks, vertex_t &gvertices) const
 
std::vector< std::vector< reco::TransientTrack > > clusterize (const std::vector< reco::TransientTrack > &tracks) const override
 
 DAClusterizerInZT_vect (const edm::ParameterSet &conf)
 
void dump (const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0) const
 
track_t fill (const std::vector< reco::TransientTrack > &tracks) const
 
bool find_nearest (double z, double t, vertex_t &y, unsigned int &k_min, double dz, double dt) const
 
double get_Tc (const vertex_t &y, int k) const
 
bool merge (vertex_t &, track_t &, double &beta) const
 
bool purge (vertex_t &, track_t &, double &, const double) const
 
void set_vtx_range (double beta, track_t &gtracks, vertex_t &gvertices) const
 
bool split (const double beta, track_t &t, vertex_t &y, double threshold=1.) const
 
unsigned int thermalize (double beta, track_t &gtracks, vertex_t &gvertices, const double delta_max, const double rho0=0.) const
 
double update (double beta, track_t &gtracks, vertex_t &gvertices, const double rho0=0) const
 
void verify (const vertex_t &v, const track_t &tks, unsigned int nv=999999, unsigned int nt=999999) const
 
std::vector< TransientVertex > vertices (const std::vector< reco::TransientTrack > &tracks, const int verbosity=0) const
 
bool zorder (vertex_t &y) const
 
- Public Member Functions inherited from TrackClusterizerInZ
 TrackClusterizerInZ ()
 
 TrackClusterizerInZ (const edm::ParameterSet &conf)
 
virtual ~TrackClusterizerInZ ()
 

Private Attributes

double betamax_
 
double betapurge_
 
double betastop_
 
unsigned int convergence_mode_
 
double coolingFactor_
 
double d0CutOff_
 
double delta_highT_
 
double delta_lowT_
 
double dtCutOff_
 
double dzCutOff_
 
unsigned int maxIterations_
 
double mintrkweight_
 
double sel_zrange_
 
double t0Max_
 
double tmerge_
 
double uniquetrkweight_
 
bool verbose_
 
double vertexSize_
 
double vertexSizeTime_
 
double zdumpcenter_
 
double zdumpwidth_
 
double zmerge_
 
const double zrange_min_ = 0.1
 

Detailed Description

Description: separates event tracks into clusters along the beam line

Version which auto-vectorizes with gcc 4.6 or newer

Definition at line 22 of file DAClusterizerInZT_vect.h.

Constructor & Destructor Documentation

◆ DAClusterizerInZT_vect()

DAClusterizerInZT_vect::DAClusterizerInZT_vect ( const edm::ParameterSet &  conf)

Definition at line 22 of file DAClusterizerInZT_vect.cc.

22  {
23  // hardcoded parameters
24  maxIterations_ = 1000;
25  mintrkweight_ = 0.5;
26 
27  // configurable debug output
28  verbose_ = conf.getUntrackedParameter<bool>("verbose", false);
29  zdumpcenter_ = conf.getUntrackedParameter<double>("zdumpcenter", 0.);
30  zdumpwidth_ = conf.getUntrackedParameter<double>("zdumpwidth", 20.);
31 
32  // configurable parameters
33  double minT = conf.getParameter<double>("Tmin");
34  double purgeT = conf.getParameter<double>("Tpurge");
35  double stopT = conf.getParameter<double>("Tstop");
36  vertexSize_ = conf.getParameter<double>("vertexSize");
37  vertexSizeTime_ = conf.getParameter<double>("vertexSizeTime");
38  coolingFactor_ = conf.getParameter<double>("coolingFactor");
39  d0CutOff_ = conf.getParameter<double>("d0CutOff");
40  dzCutOff_ = conf.getParameter<double>("dzCutOff");
41  dtCutOff_ = conf.getParameter<double>("dtCutOff");
42  t0Max_ = conf.getParameter<double>("t0Max");
43  uniquetrkweight_ = conf.getParameter<double>("uniquetrkweight");
44  zmerge_ = conf.getParameter<double>("zmerge");
45  tmerge_ = conf.getParameter<double>("tmerge");
46 
47  sel_zrange_ = conf.getParameter<double>("zrange");
48  convergence_mode_ = conf.getParameter<int>("convergence_mode");
49  delta_lowT_ = conf.getParameter<double>("delta_lowT");
50  delta_highT_ = conf.getParameter<double>("delta_highT");
51 
52  if (verbose_) {
53  std::cout << "DAClusterizerInZT_vect: mintrkweight = " << mintrkweight_ << std::endl;
54  std::cout << "DAClusterizerInZT_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl;
55  std::cout << "DAClusterizerInZT_vect: zmerge = " << zmerge_ << std::endl;
56  std::cout << "DAClusterizerInZT_vect: tmerge = " << tmerge_ << std::endl;
57  std::cout << "DAClusterizerInZT_vect: Tmin = " << minT << std::endl;
58  std::cout << "DAClusterizerInZT_vect: Tpurge = " << purgeT << std::endl;
59  std::cout << "DAClusterizerInZT_vect: Tstop = " << stopT << std::endl;
60  std::cout << "DAClusterizerInZT_vect: vertexSize = " << vertexSize_ << std::endl;
61  std::cout << "DAClusterizerInZT_vect: vertexSizeTime = " << vertexSizeTime_ << std::endl;
62  std::cout << "DAClusterizerInZT_vect: coolingFactor = " << coolingFactor_ << std::endl;
63  std::cout << "DAClusterizerInZT_vect: d0CutOff = " << d0CutOff_ << std::endl;
64  std::cout << "DAClusterizerInZT_vect: dzCutOff = " << dzCutOff_ << std::endl;
65  std::cout << "DAClusterizerInZT_vect: dtCutoff = " << dtCutOff_ << std::endl;
66  std::cout << "DAClusterizerInZT_vect: zrange = " << sel_zrange_ << std::endl;
67  std::cout << "DAClusterizerinZT_vect: convergence mode = " << convergence_mode_ << std::endl;
68  std::cout << "DAClusterizerinZT_vect: delta_highT = " << delta_highT_ << std::endl;
69  std::cout << "DAClusterizerinZT_vect: delta_lowT = " << delta_lowT_ << std::endl;
70  }
71 #ifdef DEBUG
72  std::cout << "DAClusterizerinZT_vect: DEBUGLEVEL " << DEBUGLEVEL << std::endl;
73 #endif
74 
75  if (convergence_mode_ > 1) {
76  edm::LogWarning("DAClusterizerinZT_vect")
77  << "DAClusterizerInZT_vect: invalid convergence_mode" << convergence_mode_ << " reset to default " << 0;
78  convergence_mode_ = 0;
79  }
80 
81  if (minT == 0) {
82  edm::LogWarning("DAClusterizerinZT_vect")
83  << "DAClusterizerInZT_vect: invalid Tmin" << minT << " reset to default " << 1. / betamax_;
84  } else {
85  betamax_ = 1. / minT;
86  }
87 
88  if ((purgeT > minT) || (purgeT == 0)) {
89  edm::LogWarning("DAClusterizerinZT_vect")
90  << "DAClusterizerInZT_vect: invalid Tpurge" << purgeT << " set to " << minT;
91  purgeT = minT;
92  }
93  betapurge_ = 1. / purgeT;
94 
95  if ((stopT > purgeT) || (stopT == 0)) {
96  edm::LogWarning("DAClusterizerinZT_vect")
97  << "DAClusterizerInZT_vect: invalid Tstop" << stopT << " set to " << max(1., purgeT);
98  stopT = max(1., purgeT);
99  }
100  betastop_ = 1. / stopT;
101 }

References gather_cfg::cout, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), and SiStripPI::max.

Member Function Documentation

◆ beta0()

double DAClusterizerInZT_vect::beta0 ( const double  betamax,
track_t const &  tks,
vertex_t const &  y 
) const

Definition at line 813 of file DAClusterizerInZT_vect.cc.

813  {
814  double T0 = 0; // max Tc for beta=0
815  // estimate critical temperature from beta=0 (T=inf)
816  const unsigned int nt = tks.getSize();
817  const unsigned int nv = y.getSize();
818 
819  for (unsigned int k = 0; k < nv; k++) {
820  // vertex fit at T=inf
821  double sumwz = 0;
822  double sumwt = 0;
823  double sumw_z = 0;
824  double sumw_t = 0;
825  for (unsigned int i = 0; i < nt; i++) {
826  double w_z = tks.pi_ptr[i] * tks.dz2_ptr[i];
827  double w_t = tks.pi_ptr[i] * tks.dt2_ptr[i];
828  sumwz += w_z * tks.z_ptr[i];
829  sumwt += w_t * tks.t_ptr[i];
830  sumw_z += w_z;
831  sumw_t += w_t;
832  }
833  y.z_ptr[k] = sumwz / sumw_z;
834  y.t_ptr[k] = sumwt / sumw_t;
835 
836  // estimate Tc, eventually do this in the same loop
837  double szz = 0, stt = 0, szt = 0;
838  double nuz = 0, nut = 0;
839  for (unsigned int i = 0; i < nt; i++) {
840  double dz = (tks.z_ptr[i] - y.z_ptr[k]) * tks.dz2_ptr[i];
841  double dt = (tks.t_ptr[i] - y.t_ptr[k]) * tks.dt2_ptr[i];
842  double w = tks.pi_ptr[i];
843  szz += w * dz * dz;
844  stt += w * dt * dt;
845  szt += w * dz * dt;
846  nuz += w * tks.dz2_ptr[i];
847  nut += w * tks.dt2_ptr[i];
848  }
849  double Tz = szz / nuz;
850  double Tt = 0;
851  double Tc = 0;
852  if (nut > 0) {
853  Tt = stt / nut;
854  Tc = Tz + Tt + sqrt(pow(Tz - Tt, 2) + 4 * szt * szt / nuz / nut);
855  } else {
856  Tc = 2. * Tz;
857  }
858 
859  if (Tc > T0)
860  T0 = Tc;
861  } // vertex loop (normally there should be only one vertex at beta=0)
862 
863 #ifdef DEBUG
864  if (DEBUGLEVEL > 0) {
865  std::cout << "DAClusterizerInZT_vect.beta0: Tc = " << T0 << std::endl;
866  int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
867  std::cout << "DAClusterizerInZT_vect.beta0: nstep = " << coolingsteps << std::endl;
868  }
869 #endif
870 
871  if (T0 > 1. / betamax) {
872  int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
873  return betamax * std::pow(coolingFactor_, coolingsteps);
874  } else {
875  // ensure at least one annealing step
876  return betamax * coolingFactor_;
877  }
878 }

References gather_cfg::cout, dt, DAClusterizerInZT_vect::track_t::dt2_ptr, PVValHelper::dz, DAClusterizerInZT_vect::track_t::dz2_ptr, DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, createfilelist::int, dqmdumpme::k, dqm-mbProfile::log, nt, DAClusterizerInZT_vect::track_t::pi_ptr, funct::pow(), mathSSE::sqrt(), DAClusterizerInZT_vect::track_t::t_ptr, w, and DAClusterizerInZT_vect::track_t::z_ptr.

◆ clear_vtx_range()

void DAClusterizerInZT_vect::clear_vtx_range ( track_t &  gtracks,
vertex_t &  gvertices 
) const

Definition at line 342 of file DAClusterizerInZT_vect.cc.

342  {
343  const unsigned int nt = gtracks.getSize();
344  const unsigned int nv = gvertices.getSize();
345  for (auto itrack = 0U; itrack < nt; ++itrack) {
346  gtracks.kmin[itrack] = 0;
347  gtracks.kmax[itrack] = nv;
348  }
349 }

References DAClusterizerInZT_vect::track_t::getSize(), DAClusterizerInZT_vect::vertex_t::getSize(), DAClusterizerInZT_vect::track_t::kmax, DAClusterizerInZT_vect::track_t::kmin, nt, and mitigatedMETSequence_cff::U.

◆ clusterize()

vector< vector< reco::TransientTrack > > DAClusterizerInZT_vect::clusterize ( const std::vector< reco::TransientTrack > &  tracks) const
overridevirtual

Implements TrackClusterizerInZ.

Definition at line 1330 of file DAClusterizerInZT_vect.cc.

1331  {
1332  vector<vector<reco::TransientTrack> > clusters;
1333  vector<TransientVertex>&& pv = vertices(tracks);
1334 
1335 #ifdef DEBUG
1336  if (DEBUGLEVEL > 0) {
1337  std::cout << "###################################################" << endl;
1338  std::cout << "# vectorized DAClusterizerInZT_vect::clusterize nt=" << tracks.size() << endl;
1339  std::cout << "# DAClusterizerInZT_vect::clusterize pv.size=" << pv.size() << endl;
1340  std::cout << "###################################################" << endl;
1341  }
1342 #endif
1343 
1344  if (pv.empty()) {
1345  return clusters;
1346  }
1347 
1348  // fill into clusters, don't merge
1349  for (auto k = pv.begin(); k != pv.end(); k++) {
1350  vector<reco::TransientTrack> aCluster = k->originalTracks();
1351  if (aCluster.size() > 1) {
1352  clusters.push_back(aCluster);
1353  }
1354  }
1355 
1356  return clusters;
1357 }

References bsc_activity_cfg::clusters, gather_cfg::cout, dqmdumpme::k, MetAnalyzer::pv(), PDWG_EXOHSCP_cff::tracks, and pwdgSkimBPark_cfi::vertices.

◆ dump()

void DAClusterizerInZT_vect::dump ( const double  beta,
const vertex_t &  y,
const track_t &  tks,
const int  verbosity = 0 
) const

Definition at line 1359 of file DAClusterizerInZT_vect.cc.

1359  {
1360 #ifdef DEBUG
1361  const unsigned int nv = y.getSize();
1362  const unsigned int nt = tks.getSize();
1363 
1364  std::vector<unsigned int> iz;
1365  for (unsigned int j = 0; j < nt; j++) {
1366  iz.push_back(j);
1367  }
1368  std::sort(iz.begin(), iz.end(), [tks](unsigned int a, unsigned int b) { return tks.z_ptr[a] < tks.z_ptr[b]; });
1369  std::cout << std::endl;
1370  std::cout << "-----DAClusterizerInZT::dump ----" << nv << " clusters " << std::endl;
1371  string h = " ";
1372  std::cout << h << " k= ";
1373  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1374  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) < zdumpwidth_) {
1375  std::cout << setw(8) << fixed << ivertex;
1376  }
1377  }
1378  std::cout << endl;
1379 
1380  std::cout << h << " z= ";
1381  std::cout << setprecision(4);
1382  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1383  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) < zdumpwidth_) {
1384  std::cout << setw(8) << fixed << y.z_ptr[ivertex];
1385  }
1386  }
1387  std::cout << endl;
1388 
1389  std::cout << h << " t= ";
1390  std::cout << setprecision(4);
1391  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1392  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) < zdumpwidth_) {
1393  std::cout << setw(8) << fixed << y.t_ptr[ivertex];
1394  }
1395  }
1396  std::cout << endl;
1397 
1398  std::cout << "T=" << setw(15) << 1. / beta << " Tmin =" << setw(10) << 1. / betamax_
1399  << " Tc= ";
1400  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1401  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) < zdumpwidth_) {
1402  double Tc = get_Tc(y, ivertex);
1403  std::cout << setw(8) << fixed << setprecision(1) << Tc;
1404  }
1405  }
1406  std::cout << endl;
1407 
1408  std::cout << h << "pk= ";
1409  double sumpk = 0;
1410  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1411  sumpk += y.pk_ptr[ivertex];
1412  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) > zdumpwidth_)
1413  continue;
1414  std::cout << setw(8) << setprecision(4) << fixed << y.pk_ptr[ivertex];
1415  }
1416  std::cout << endl;
1417 
1418  std::cout << h << "nt= ";
1419  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1420  sumpk += y.pk_ptr[ivertex];
1421  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) > zdumpwidth_)
1422  continue;
1423  std::cout << setw(8) << setprecision(1) << fixed << y.pk_ptr[ivertex] * nt;
1424  }
1425  std::cout << endl;
1426 
1427  if (verbosity > 0) {
1428  double E = 0, F = 0;
1429  std::cout << endl;
1430  std::cout << "---- z +/- dz t +/- dt ip +/-dip pt phi eta weights ----"
1431  << endl;
1432  std::cout << setprecision(4);
1433  for (unsigned int i0 = 0; i0 < nt; i0++) {
1434  unsigned int i = iz[i0];
1435  if (tks.Z_sum_ptr[i] > 0) {
1436  F -= std::log(tks.Z_sum_ptr[i]) / beta;
1437  }
1438  double tz = tks.z_ptr[i];
1439 
1440  if (std::fabs(tz - zdumpcenter_) > zdumpwidth_)
1441  continue;
1442  std::cout << setw(4) << i << ")" << setw(8) << fixed << setprecision(4) << tz << " +/-" << setw(6)
1443  << sqrt(1. / tks.dz2_ptr[i]);
1444 
1445  if (tks.dt2_ptr[i] > 0) {
1446  std::cout << setw(8) << fixed << setprecision(4) << tks.t_ptr[i] << " +/-" << setw(6)
1447  << sqrt(1. / tks.dt2_ptr[i]);
1448  } else {
1449  std::cout << " ";
1450  }
1451 
1452  if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) {
1453  std::cout << " *";
1454  } else {
1455  std::cout << " ";
1456  }
1457  if (tks.tt[i]->track().hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) {
1458  std::cout << "+";
1459  } else {
1460  std::cout << "-";
1461  }
1462  std::cout << setw(1)
1463  << tks.tt[i]
1464  ->track()
1465  .hitPattern()
1466  .pixelBarrelLayersWithMeasurement(); // see DataFormats/TrackReco/interface/HitPattern.h
1467  std::cout << setw(1) << tks.tt[i]->track().hitPattern().pixelEndcapLayersWithMeasurement();
1468  std::cout << setw(1) << hex
1469  << tks.tt[i]->track().hitPattern().trackerLayersWithMeasurement() -
1470  tks.tt[i]->track().hitPattern().pixelLayersWithMeasurement()
1471  << dec;
1472  std::cout << "=" << setw(1) << hex
1473  << tks.tt[i]->track().hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << dec;
1474 
1475  Measurement1D IP = tks.tt[i]->stateAtBeamLine().transverseImpactParameter();
1476  std::cout << setw(8) << IP.value() << "+/-" << setw(6) << IP.error();
1477  std::cout << " " << setw(6) << setprecision(2) << tks.tt[i]->track().pt() * tks.tt[i]->track().charge();
1478  std::cout << " " << setw(5) << setprecision(2) << tks.tt[i]->track().phi() << " " << setw(5) << setprecision(2)
1479  << tks.tt[i]->track().eta();
1480 
1481  double sump = 0.;
1482  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1483  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) > zdumpwidth_)
1484  continue;
1485 
1486  if ((tks.pi_ptr[i] > 0) && (tks.Z_sum_ptr[i] > 0)) {
1487  double p =
1488  y.pk_ptr[ivertex] *
1489  exp(-beta *
1490  Eik(tks.z_ptr[i], y.z_ptr[ivertex], tks.dz2_ptr[i], tks.t_ptr[i], y.t_ptr[ivertex], tks.dt2_ptr[i])) /
1491  tks.Z_sum_ptr[i];
1492 
1493  if ((ivertex >= tks.kmin[i]) && (ivertex < tks.kmax[i])) {
1494  if (p > 0.0001) {
1495  std::cout << setw(8) << setprecision(3) << p;
1496  } else {
1497  std::cout << " _ "; // tiny but in the cluster list
1498  }
1499  E +=
1500  p * Eik(tks.z_ptr[i], y.z_ptr[ivertex], tks.dz2_ptr[i], tks.t_ptr[i], y.t_ptr[ivertex], tks.dt2_ptr[i]);
1501  sump += p;
1502  } else {
1503  if (p > 0.1) {
1504  std::cout << "XXXXXXXX"; // we have an inconsistency here
1505  } else if (p > 0.0001) {
1506  std::cout << "X" << setw(6) << setprecision(3) << p << "X";
1507  } else {
1508  std::cout << " . "; // not in the cluster list
1509  }
1510  }
1511  } else {
1512  std::cout << " ";
1513  }
1514  }
1515  std::cout << " ( " << std::setw(3) << tks.kmin[i] << "," << std::setw(3) << tks.kmax[i] - 1 << " ) ";
1516  std::cout << endl;
1517  }
1518  std::cout << " ";
1519  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1520  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) < zdumpwidth_) {
1521  std::cout << " " << setw(3) << ivertex << " ";
1522  }
1523  }
1524  std::cout << endl;
1525  std::cout << " z= ";
1526  std::cout << setprecision(4);
1527  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1528  if (std::fabs(y.z_ptr[ivertex] - zdumpcenter_) < zdumpwidth_) {
1529  std::cout << setw(8) << fixed << y.z_ptr[ivertex];
1530  }
1531  }
1532  std::cout << endl;
1533  std::cout << endl
1534  << "T=" << 1 / beta << " E=" << E << " n=" << y.getSize() << " F= " << F << endl
1535  << "----------" << endl;
1536  }
1537 #endif
1538 }

References a, b, zMuMuMuonUserData::beta, gather_cfg::cout, TauDecayModes::dec, DAClusterizerInZT_vect::track_t::dt2_ptr, DAClusterizerInZT_vect::track_t::dz2_ptr, JetChargeProducer_cfi::exp, F(), alignBH_cfg::fixed, DAClusterizerInZT_vect::track_t::getSize(), reco::TrackBase::highPurity, mps_fire::i, listHistos::IP, dqmiolumiharvest::j, DAClusterizerInZT_vect::track_t::kmax, DAClusterizerInZT_vect::track_t::kmin, dqm-mbProfile::log, reco::HitPattern::MISSING_OUTER_HITS, nt, AlCaHLTBitMon_ParallelJobs::p, DAClusterizerInZT_vect::track_t::pi_ptr, GeomDetEnumerators::PixelBarrel, mathSSE::sqrt(), DAClusterizerInZT_vect::track_t::t_ptr, DAClusterizerInZT_vect::track_t::tt, HIPAlignmentAlgorithm_cfi::verbosity, DAClusterizerInZT_vect::track_t::z_ptr, and DAClusterizerInZT_vect::track_t::Z_sum_ptr.

◆ fill()

DAClusterizerInZT_vect::track_t DAClusterizerInZT_vect::fill ( const std::vector< reco::TransientTrack > &  tracks) const

Definition at line 214 of file DAClusterizerInZT_vect.cc.

214  {
215  // prepare track data for clustering
216  track_t tks;
217  for (const auto& tk : tracks) {
218  if (!tk.isValid())
219  continue;
220  double t_pi = 1.;
221  double t_z = tk.stateAtBeamLine().trackStateAtPCA().position().z();
222  double t_t = tk.timeExt();
223 
224  if (std::fabs(t_z) > 1000.)
225  continue;
226  /* for comparison with 1d clustering, keep such tracks without timing info, see below
227  if (std::abs(t_t) > t0Max_)
228  continue;
229  */
230 
231  auto const& t_mom = tk.stateAtBeamLine().trackStateAtPCA().momentum();
232  // get the beam-spot
233  reco::BeamSpot beamspot = tk.stateAtBeamLine().beamSpot();
234  double t_dz2 = std::pow(tk.track().dzError(), 2) // track errror
235  + (std::pow(beamspot.BeamWidthX() * t_mom.x(), 2) + std::pow(beamspot.BeamWidthY() * t_mom.y(), 2)) *
236  std::pow(t_mom.z(), 2) / std::pow(t_mom.perp2(), 2) // beam spot width
237  + std::pow(vertexSize_, 2); // intrinsic vertex size, safer for outliers and short lived decays
238  t_dz2 = 1. / t_dz2;
239  if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits<double>::min()) {
240  std::cout << "DAClusterizerInZT_vect.fill rejected track t_dz2 " << t_dz2 << std::endl;
241  continue;
242  }
243 
244  double t_dt2 =
245  std::pow(tk.dtErrorExt(), 2.) +
246  std::pow(vertexSizeTime_, 2.); // the ~injected~ timing error, need to add a small minimum vertex size in time
247  if ((tk.dtErrorExt() > 0.3) || (std::abs(t_t) > t0Max_)) {
248  t_dt2 = 0; // tracks with no time measurement
249  } else {
250  t_dt2 = 1. / t_dt2;
251  if (edm::isNotFinite(t_dt2) || t_dt2 < std::numeric_limits<double>::min()) {
252  std::cout << "DAClusterizerInZT_vect.fill rejected track t_dt2 " << t_dt2 << std::endl;
253  continue;
254  }
255  }
256 
257  if (d0CutOff_ > 0) {
258  Measurement1D atIP = tk.stateAtBeamLine().transverseImpactParameter(); // error contains beamspot
259  t_pi = 1. / (1. + local_exp(std::pow(atIP.value() / atIP.error(), 2) -
260  std::pow(d0CutOff_, 2))); // reduce weight for high ip tracks
261  if (edm::isNotFinite(t_pi) || t_pi < std::numeric_limits<double>::epsilon()) {
262  std::cout << "DAClusterizerInZT_vect.fill rejected track t_pu " << t_pi << std::endl;
263  continue; // usually is > 0.99
264  }
265  }
266 
267  tks.addItem(t_z, t_t, t_dz2, t_dt2, &tk, t_pi);
268  }
269  tks.extractRaw();
270 #ifdef DEBUG
271  if (DEBUGLEVEL > 0) {
272  std::cout << "Track count (ZT) filled " << tks.getSize() << " initial " << tracks.size() << std::endl;
273  }
274 #endif
275 
276  return tks;
277 }

References funct::abs(), DAClusterizerInZT_vect::track_t::addItem(), HLT_2018_cff::atIP, gather_cfg::cout, geometryDiff::epsilon, DAClusterizerInZT_vect::track_t::extractRaw(), DAClusterizerInZT_vect::track_t::getSize(), edm::isNotFinite(), min(), funct::pow(), AlignmentPI::t_z, and PDWG_EXOHSCP_cff::tracks.

◆ find_nearest()

bool DAClusterizerInZT_vect::find_nearest ( double  z,
double  t,
vertex_t &  y,
unsigned int &  k_min,
double  dz,
double  dt 
) const

Definition at line 586 of file DAClusterizerInZT_vect.cc.

587  {
588  // find the cluster nearest to (z,t)
589  // distance measure is delta = (delta_z / dz )**2 + (delta_t/ d_t)**2
590  // assumes that clusters are ordered n z
591  // return value is false if no neighbour with distance < 1 is found
592 
593  unsigned int nv = y.getSize();
594  if (nv < 2) {
595  k_min = 0;
596  return false;
597  }
598 
599  // find nearest in z, binary search later
600  unsigned int k = 0;
601  for (unsigned int k0 = 1; k0 < nv; k0++) {
602  if (std::abs(y.z_ptr[k0] - z) < std::abs(y.z_ptr[k] - z)) {
603  k = k0;
604  }
605  }
606 
607  double delta_min = 1.;
608 
609  //search left
610  unsigned int k1 = k;
611  while ((k1 > 0) && ((y.z[k] - y.z[--k1]) < dz)) {
612  auto delta = std::pow((y.z_ptr[k] - y.z_ptr[k1]) / dz, 2) + std::pow((y.t_ptr[k] - y.t_ptr[k1]) / dt, 2);
613  if (delta < delta_min) {
614  k_min = k1;
615  delta_min = delta;
616  }
617  }
618 
619  //search right
620  k1 = k;
621  while (((++k1) < nv) && ((y.z[k1] - y.z[k]) < dz)) {
622  auto delta = std::pow((y.z_ptr[k1] - y.z_ptr[k]) / dz, 2) + std::pow((y.t_ptr[k1] - y.t_ptr[k]) / dt, 2);
623  if (delta < delta_min) {
624  k_min = k1;
625  delta_min = delta;
626  }
627  }
628 
629  return (delta_min < 1.);
630 }

References funct::abs(), dumpMFGeometry_cfg::delta, dt, PVValHelper::dz, dqmdumpme::k, reco::ParticleMasses::k0, and funct::pow().

◆ get_Tc()

double DAClusterizerInZT_vect::get_Tc ( const vertex_t &  y,
int  k 
) const

Definition at line 880 of file DAClusterizerInZT_vect.cc.

880  {
881  double Tz = y.szz_ptr[k] / y.nuz_ptr[k]; // actually 0.5*Tc(z)
882  double Tt = 0.;
883  if (y.nut_ptr[k] > 0) {
884  Tt = y.stt_ptr[k] / y.nut_ptr[k];
885  double mx = y.szt_ptr[k] / y.nuz_ptr[k] * y.szt_ptr[k] / y.nut_ptr[k];
886  return Tz + Tt + sqrt(pow(Tz - Tt, 2) + 4 * mx);
887  }
888  return 2 * Tz;
889 }

References dqmdumpme::k, funct::pow(), and mathSSE::sqrt().

◆ merge()

bool DAClusterizerInZT_vect::merge ( vertex_t &  y,
track_t &  tks,
double &  beta 
) const

Definition at line 683 of file DAClusterizerInZT_vect.cc.

683  {
684  // merge clusters that collapsed or never separated,
685  // return true if vertices were merged, false otherwise
686  const unsigned int nv = y.getSize();
687 
688  if (nv < 2)
689  return false;
690 
691  // merge the smallest distance clusters first
692  unsigned int k1_min = 0, k2_min = 0;
693  double delta_min = 0;
694 
695  for (unsigned int k1 = 0; (k1 + 1) < nv; k1++) {
696  unsigned int k2 = k1;
697  while ((++k2 < nv) && (std::fabs(y.z[k2] - y.z_ptr[k1]) < zmerge_)) {
698  auto delta =
699  std::pow((y.z_ptr[k2] - y.z_ptr[k1]) / zmerge_, 2) + std::pow((y.t_ptr[k2] - y.t_ptr[k1]) / tmerge_, 2);
700  if ((delta < delta_min) || (k1_min == k2_min)) {
701  k1_min = k1;
702  k2_min = k2;
703  delta_min = delta;
704  }
705  }
706  }
707 
708  if ((k1_min == k2_min) || (delta_min > 1)) {
709  return false;
710  }
711 
712  double rho = y.pk_ptr[k1_min] + y.pk_ptr[k2_min];
713 
714 #ifdef DEBUG
715  assert((k1_min < nv) && (k2_min < nv));
716  if (DEBUGLEVEL > 1) {
717  std::cout << "merging (" << setw(8) << fixed << setprecision(4) << y.z_ptr[k1_min] << ',' << y.t_ptr[k1_min]
718  << ") and (" << y.z_ptr[k2_min] << ',' << y.t_ptr[k2_min] << ")"
719  << " idx=" << k1_min << "," << k2_min << std::endl;
720  }
721 #endif
722 
723  if (rho > 0) {
724  y.z_ptr[k1_min] = (y.pk_ptr[k1_min] * y.z_ptr[k1_min] + y.pk_ptr[k2_min] * y.z_ptr[k2_min]) / rho;
725  y.t_ptr[k1_min] = (y.pk_ptr[k1_min] * y.t_ptr[k1_min] + y.pk_ptr[k2_min] * y.t_ptr[k2_min]) / rho;
726 #ifdef USEVTXDT2
727  y.dt2_ptr[k1_min] = (y.pk_ptr[k1_min] * y.dt2_ptr[k1_min] + y.pk_ptr[k2_min] * y.dt2_ptr[k2_min]) / rho;
728 #endif
729  } else {
730  y.z_ptr[k1_min] = 0.5 * (y.z_ptr[k1_min] + y.z_ptr[k2_min]);
731  y.t_ptr[k1_min] = 0.5 * (y.t_ptr[k1_min] + y.t_ptr[k2_min]);
732  }
733  y.pk_ptr[k1_min] = rho;
734  y.removeItem(k2_min, tks);
735 
736  zorder(y);
737  set_vtx_range(beta, tks, y);
738  y.extractRaw();
739  return true;
740 }

References cms::cuda::assert(), zMuMuMuonUserData::beta, gather_cfg::cout, dumpMFGeometry_cfg::delta, alignBH_cfg::fixed, and funct::pow().

◆ purge()

bool DAClusterizerInZT_vect::purge ( vertex_t &  y,
track_t &  tks,
double &  rho0,
const double  beta 
) const

Definition at line 742 of file DAClusterizerInZT_vect.cc.

742  {
743  constexpr double eps = 1.e-100;
744  // eliminate clusters with only one significant/unique track
745  const unsigned int nv = y.getSize();
746  const unsigned int nt = tks.getSize();
747 
748  if (nv < 2)
749  return false;
750 
751  double sumpmin = nt;
752  unsigned int k0 = nv;
753 
754  int nUnique = 0;
755  double sump = 0;
756 
757  std::vector<double> inverse_zsums(nt), arg_cache(nt), eik_cache(nt);
758  double* pinverse_zsums;
759  double* parg_cache;
760  double* peik_cache;
761  pinverse_zsums = inverse_zsums.data();
762  parg_cache = arg_cache.data();
763  peik_cache = eik_cache.data();
764  for (unsigned i = 0; i < nt; ++i) {
765  inverse_zsums[i] = tks.Z_sum_ptr[i] > eps ? 1. / tks.Z_sum_ptr[i] : 0.0;
766  }
767 
768  for (unsigned int k = 0; k < nv; ++k) {
769  nUnique = 0;
770  sump = 0;
771 
772  const double pmax = y.pk_ptr[k] / (y.pk_ptr[k] + rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_));
773  const double pcut = uniquetrkweight_ * pmax;
774  for (unsigned i = 0; i < nt; ++i) {
775  const auto track_z = tks.z_ptr[i];
776  const auto track_t = tks.t_ptr[i];
777  const auto botrack_dz2 = -beta * tks.dz2_ptr[i];
778  const auto botrack_dt2 = -beta * tks.dt2_ptr[i]; // FIXME usevtxdt2?
779 
780  const auto mult_resz = track_z - y.z_ptr[k];
781  const auto mult_rest = track_t - y.t_ptr[k];
782  parg_cache[i] = botrack_dz2 * (mult_resz * mult_resz) + botrack_dt2 * (mult_rest * mult_rest);
783  }
784  local_exp_list(parg_cache, peik_cache, nt);
785  for (unsigned int i = 0; i < nt; ++i) {
786  const double p = y.pk_ptr[k] * peik_cache[i] * pinverse_zsums[i];
787  sump += p;
788  nUnique += ((p > pcut) & (tks.pi_ptr[i] > 0));
789  }
790 
791  if ((nUnique < 2) && (sump < sumpmin)) {
792  sumpmin = sump;
793  k0 = k;
794  }
795  }
796 
797  if (k0 != nv) {
798 #ifdef DEBUG
799  assert(k0 < y.getSize());
800  if (DEBUGLEVEL > 1) {
801  std::cout << "eliminating prototype at " << std::setw(10) << std::setprecision(4) << y.z_ptr[k0] << ","
802  << y.t_ptr[k0] << " with sump=" << sumpmin << " rho*nt =" << y.pk_ptr[k0] * nt << endl;
803  }
804 #endif
805  y.removeItem(k0, tks);
806  set_vtx_range(beta, tks, y);
807  return true;
808  } else {
809  return false;
810  }
811 }

References cms::cuda::assert(), zMuMuMuonUserData::beta, gather_cfg::cout, DAClusterizerInZT_vect::track_t::dt2_ptr, DAClusterizerInZT_vect::track_t::dz2_ptr, DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, dqmdumpme::k, reco::ParticleMasses::k0, nt, AlCaHLTBitMon_ParallelJobs::p, DAClusterizerInZT_vect::track_t::pi_ptr, DAClusterizerInZT_vect::track_t::t_ptr, DAClusterizerInZT_vect::track_t::z_ptr, and DAClusterizerInZT_vect::track_t::Z_sum_ptr.

◆ set_vtx_range()

void DAClusterizerInZT_vect::set_vtx_range ( double  beta,
track_t &  gtracks,
vertex_t &  gvertices 
) const

Definition at line 285 of file DAClusterizerInZT_vect.cc.

285  {
286  const unsigned int nv = gvertices.getSize();
287  const unsigned int nt = gtracks.getSize();
288 
289  if (nv == 0) {
290  edm::LogWarning("DAClusterizerinZT_vect") << "empty cluster list in set_vtx_range";
291  return;
292  }
293 
294  for (auto itrack = 0U; itrack < nt; ++itrack) {
295  double zrange = max(sel_zrange_ / sqrt(beta * gtracks.dz2[itrack]), zrange_min_);
296 
297  double zmin = gtracks.z_ptr[itrack] - zrange;
298  unsigned int kmin = min(nv - 1, gtracks.kmin[itrack]);
299  // find the smallest vertex_z that is larger than zmin
300  if (gvertices.z_ptr[kmin] > zmin) {
301  while ((kmin > 0) && (gvertices.z_ptr[kmin - 1] > zmin)) {
302  kmin--;
303  }
304  } else {
305  while ((kmin < (nv - 1)) && (gvertices.z_ptr[kmin] < zmin)) {
306  kmin++;
307  }
308  }
309 
310  double zmax = gtracks.z_ptr[itrack] + zrange;
311  unsigned int kmax = min(nv - 1, gtracks.kmax[itrack] - 1);
312  // note: kmax points to the last vertex in the range, while gtracks.kmax points to the entry BEHIND the last vertex
313  // find the largest vertex_z that is smaller than zmax
314  if (gvertices.z_ptr[kmax] < zmax) {
315  while ((kmax < (nv - 1)) && (gvertices.z_ptr[kmax + 1] < zmax)) {
316  kmax++;
317  }
318  } else {
319  while ((kmax > 0) && (gvertices.z_ptr[kmax] > zmax)) {
320  kmax--;
321  }
322  }
323 
324  if (kmin <= kmax) {
325  gtracks.kmin[itrack] = kmin;
326  gtracks.kmax[itrack] = kmax + 1;
327  } else {
328  gtracks.kmin[itrack] = max(0U, min(kmin, kmax));
329  gtracks.kmax[itrack] = min(nv, max(kmin, kmax) + 1);
330  }
331 
332 #ifdef DEBUG
333  if (gtracks.kmin[itrack] >= gtracks.kmax[itrack]) {
334  cout << "set_vtx_range trk = " << itrack << " kmin,kmax=" << kmin << "," << kmax
335  << " gtrack.kmin,kmax = " << gtracks.kmin[itrack] << "," << gtracks.kmax[itrack] << " zrange = " << zrange
336  << endl;
337  }
338 #endif
339  }
340 }

References zMuMuMuonUserData::beta, gather_cfg::cout, DAClusterizerInZT_vect::track_t::dz2, DAClusterizerInZT_vect::track_t::getSize(), DAClusterizerInZT_vect::vertex_t::getSize(), DAClusterizerInZT_vect::track_t::kmax, DAClusterizerInZT_vect::track_t::kmin, SiStripPI::max, min(), nt, mathSSE::sqrt(), mitigatedMETSequence_cff::U, DAClusterizerInZT_vect::track_t::z_ptr, DAClusterizerInZT_vect::vertex_t::z_ptr, SiStripMonitorCluster_cfi::zmax, SiStripMonitorCluster_cfi::zmin, and TkClusParameters_cff::zrange.

◆ split()

bool DAClusterizerInZT_vect::split ( const double  beta,
track_t &  t,
vertex_t &  y,
double  threshold = 1. 
) const

Definition at line 891 of file DAClusterizerInZT_vect.cc.

891  {
892  // split only critical vertices (Tc >~ T=1/beta <==> beta*Tc>~1)
893  // an update must have been made just before doing this (same beta, no merging)
894  // returns true if at least one cluster was split
895 
896  constexpr double epsilonz = 1e-3; // minimum split size z
897  constexpr double epsilont = 1e-2; // minimum split size t
898  unsigned int nv = y.getSize();
899  const double twoBeta = 2.0 * beta;
900 
901  // avoid left-right biases by splitting highest Tc first
902 
903  std::vector<std::pair<double, unsigned int> > critical;
904  for (unsigned int k = 0; k < nv; k++) {
905  double Tc = get_Tc(y, k);
906  if (beta * Tc > threshold) {
907  critical.push_back(make_pair(Tc, k));
908  }
909  }
910  if (critical.empty())
911  return false;
912 
913  std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int> >());
914 
915  bool split = false;
916  const unsigned int nt = tks.getSize();
917 
918  for (unsigned int ic = 0; ic < critical.size(); ic++) {
919  unsigned int k = critical[ic].second;
920 
921  // split direction in the (z,t)-plane
922 
923  double Mzz = y.nuz_ptr[k] - twoBeta * y.szz_ptr[k];
924  double Mtt = y.nut_ptr[k] - twoBeta * y.stt_ptr[k];
925  double Mzt = -twoBeta * y.szt_ptr[k];
926  const double twoMzt = 2.0 * Mzt;
927  double D = sqrt(pow(Mtt - Mzz, 2) + twoMzt * twoMzt);
928  double q1 = atan2(-Mtt + Mzz + D, -twoMzt);
929  double l1 = 0.5 * (-Mzz - Mtt + D);
930  double l2 = 0.5 * (-Mzz - Mtt - D);
931  if ((std::abs(l1) < 1e-4) && (std::abs(l2) < 1e-4)) {
932  edm::LogWarning("DAClusterizerInZT_vect") << "warning, bad eigenvalues! idx=" << k << " z= " << y.z_ptr[k]
933  << " Mzz=" << Mzz << " Mtt=" << Mtt << " Mzt=" << Mzt << endl;
934  }
935 
936  double qsplit = q1;
937  double cq = cos(qsplit);
938  double sq = sin(qsplit);
939  if (cq < 0) {
940  cq = -cq;
941  sq = -sq;
942  } // prefer cq>0 to keep z-ordering
943 
944  // estimate subcluster positions and weight
945  double p1 = 0, z1 = 0, t1 = 0, wz1 = 0, wt1 = 0;
946  double p2 = 0, z2 = 0, t2 = 0, wz2 = 0, wt2 = 0;
947  for (unsigned int i = 0; i < nt; ++i) {
948  if (tks.Z_sum_ptr[i] > 1.e-100) {
949  double lr = (tks.z_ptr[i] - y.z_ptr[k]) * cq + (tks.t[i] - y.t_ptr[k]) * sq;
950  // winner-takes-all, usually overestimates splitting
951  double tl = lr < 0 ? 1. : 0.;
952  double tr = 1. - tl;
953 
954  // soften it, especially at low T
955  double arg = lr * std::sqrt(beta * (cq * cq * tks.dz2_ptr[i] + sq * sq * tks.dt2_ptr[i]));
956  if (std::abs(arg) < 20) {
957  double t = local_exp(-arg);
958  tl = t / (t + 1.);
959  tr = 1 / (t + 1.);
960  }
961 
962  double p =
963  y.pk_ptr[k] * tks.pi_ptr[i] *
964  local_exp(-beta * Eik(tks.z_ptr[i], y.z_ptr[k], tks.dz2_ptr[i], tks.t_ptr[i], y.t_ptr[k], tks.dt2_ptr[i])) /
965  tks.Z_sum_ptr[i];
966  double wz = p * tks.dz2_ptr[i];
967  double wt = p * tks.dt2_ptr[i];
968  p1 += p * tl;
969  z1 += wz * tl * tks.z_ptr[i];
970  t1 += wt * tl * tks.t_ptr[i];
971  wz1 += wz * tl;
972  wt1 += wt * tl;
973  p2 += p * tr;
974  z2 += wz * tr * tks.z_ptr[i];
975  t2 += wt * tr * tks.t_ptr[i];
976  wz2 += wz * tr;
977  wt2 += wt * tr;
978  }
979  }
980 
981  if (wz1 > 0) {
982  z1 /= wz1;
983  } else {
984  z1 = y.z_ptr[k] - epsilonz * cq;
985  edm::LogWarning("DAClusterizerInZT_vect") << "warning, wz1 = " << scientific << wz1 << endl;
986  }
987  if (wt1 > 0) {
988  t1 /= wt1;
989  } else {
990  t1 = y.t_ptr[k] - epsilont * sq;
991  edm::LogWarning("DAClusterizerInZT_vect") << "warning, wt1 = " << scientific << wt1 << endl;
992  }
993  if (wz2 > 0) {
994  z2 /= wz2;
995  } else {
996  z2 = y.z_ptr[k] + epsilonz * cq;
997  edm::LogWarning("DAClusterizerInZT_vect") << "warning, wz2 = " << scientific << wz2 << endl;
998  }
999  if (wt2 > 0) {
1000  t2 /= wt2;
1001  } else {
1002  t2 = y.t_ptr[k] + epsilont * sq;
1003  edm::LogWarning("DAClusterizerInZT_vect") << "warning, wt2 = " << scientific << wt2 << endl;
1004  }
1005 
1006  unsigned int k_min1 = k, k_min2 = k;
1007  constexpr double spliteps = 1e-8;
1008  while (((find_nearest(z1, t1, y, k_min1, epsilonz, epsilont) && (k_min1 != k)) ||
1009  (find_nearest(z2, t2, y, k_min2, epsilonz, epsilont) && (k_min2 != k))) &&
1010  (std::abs(z2 - z1) > spliteps || std::abs(t2 - t1) > spliteps)) {
1011  z1 = 0.5 * (z1 + y.z_ptr[k]);
1012  t1 = 0.5 * (t1 + y.t_ptr[k]);
1013  z2 = 0.5 * (z2 + y.z_ptr[k]);
1014  t2 = 0.5 * (t2 + y.t_ptr[k]);
1015  }
1016 
1017 #ifdef DEBUG
1018  assert(k < nv);
1019  if (DEBUGLEVEL > 1) {
1020  if (std::fabs(y.z_ptr[k] - zdumpcenter_) < zdumpwidth_) {
1021  std::cout << " T= " << std::setw(10) << std::setprecision(1) << 1. / beta << " Tc= " << critical[ic].first
1022  << " direction =" << std::setprecision(4) << qsplit << " splitting (" << std::setw(8) << std::fixed
1023  << std::setprecision(4) << y.z_ptr[k] << "," << y.t_ptr[k] << ")"
1024  << " --> (" << z1 << ',' << t1 << "),(" << z2 << ',' << t2 << ") [" << p1 << "," << p2 << "]";
1025  if (std::fabs(z2 - z1) > epsilonz || std::fabs(t2 - t1) > epsilont) {
1026  std::cout << std::endl;
1027  } else {
1028  std::cout << " rejected " << std::endl;
1029  }
1030  }
1031  }
1032 #endif
1033 
1034  if (z1 > z2) {
1035  edm::LogInfo("DAClusterizerInZT") << "warning swapping z in split qsplit=" << qsplit << " cq=" << cq
1036  << " sq=" << sq << endl;
1037  auto ztemp = z1;
1038  auto ttemp = t1;
1039  auto ptemp = p1;
1040  z1 = z2;
1041  t1 = t2;
1042  p1 = p2;
1043  z2 = ztemp;
1044  t2 = ttemp;
1045  p2 = ptemp;
1046  }
1047 
1048  // split if the new subclusters are significantly separated
1049  if (std::fabs(z2 - z1) > epsilonz || std::fabs(t2 - t1) > epsilont) {
1050  split = true;
1051  double pk1 = p1 * y.pk_ptr[k] / (p1 + p2);
1052  double pk2 = p2 * y.pk_ptr[k] / (p1 + p2);
1053 
1054  // replace the original by (z2,t2)
1055  y.removeItem(k, tks);
1056  unsigned int k2 = y.insertOrdered(z2, t2, pk2, tks);
1057 
1058 #ifdef DEBUG
1059  if (k2 < k) {
1060  std::cout << "unexpected z-ordering in split" << std::endl;
1061  }
1062 #endif
1063 
1064  // adjust pointers if necessary
1065  if (!(k == k2)) {
1066  for (unsigned int jc = ic; jc < critical.size(); jc++) {
1067  if (critical[jc].second > k) {
1068  critical[jc].second--;
1069  }
1070  if (critical[jc].second >= k2) {
1071  critical[jc].second++;
1072  }
1073  }
1074  }
1075 
1076  // insert (z1,t1) where it belongs
1077  unsigned int k1 = y.insertOrdered(z1, t1, pk1, tks);
1078  nv++;
1079 
1080  // adjust remaining pointers
1081  for (unsigned int jc = ic; jc < critical.size(); jc++) {
1082  if (critical[jc].second >= k1) {
1083  critical[jc].second++;
1084  } // need to backport the ">-"?
1085  }
1086 
1087  } else {
1088 #ifdef DEBUG
1089  std::cout << "warning ! split rejected, too small." << endl;
1090 #endif
1091  }
1092  }
1093  return split;
1094 }

References funct::abs(), cms::cuda::assert(), zMuMuMuonUserData::beta, funct::cos(), gather_cfg::cout, hlt_jetmet_dqm_QT_fromfile_cfg::critical, DAClusterizerInZT_vect::track_t::dt2_ptr, DAClusterizerInZT_vect::track_t::dz2_ptr, MillePedeFileConverter_cfg::e, alignBH_cfg::fixed, DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, dqmdumpme::k, nt, AlCaHLTBitMon_ParallelJobs::p, p1, p2, DAClusterizerInZT_vect::track_t::pi_ptr, funct::pow(), q1, edm::second(), funct::sin(), cms::dd::split(), mathSSE::sqrt(), DAClusterizerInZT_vect::track_t::t, OrderedSet::t, RandomServiceHelper::t1, RandomServiceHelper::t2, DAClusterizerInZT_vect::track_t::t_ptr, remoteMonitoring_LED_IterMethod_cfg::threshold, testProducerWithPsetDescEmpty_cfi::z2, DAClusterizerInZT_vect::track_t::z_ptr, and DAClusterizerInZT_vect::track_t::Z_sum_ptr.

◆ thermalize()

unsigned int DAClusterizerInZT_vect::thermalize ( double  beta,
track_t &  gtracks,
vertex_t &  gvertices,
const double  delta_max,
const double  rho0 = 0. 
) const

Definition at line 632 of file DAClusterizerInZT_vect.cc.

633  {
634  unsigned int niter = 0;
635  double delta = 0;
636  double delta_max = delta_lowT_;
637 
638  if (convergence_mode_ == 0) {
639  delta_max = delta_max0;
640  } else if (convergence_mode_ == 1) {
641  delta_max = delta_lowT_ / sqrt(std::max(beta, 1.0));
642  }
643 
644  zorder(v);
645  set_vtx_range(beta, tks, v);
646 
647  double delta_sum_range = 0; // accumulate max(|delta-z|) as a lower bound
648  std::vector<double> z0 = v.z;
649 
650  while (niter++ < maxIterations_) {
651  delta = update(beta, tks, v, rho0);
652  delta_sum_range += delta;
653 
654  if (delta_sum_range > zrange_min_) {
655  for (unsigned int k = 0; k < v.getSize(); k++) {
656  if (std::abs(v.z[k] - z0[k]) > zrange_min_) {
657  zorder(v);
658  set_vtx_range(beta, tks, v);
659  delta_sum_range = 0;
660  z0 = v.z;
661  break;
662  }
663  }
664  }
665 
666  if (delta < delta_max) {
667  break;
668  }
669  }
670 
671 #ifdef DEBUG
672  if (DEBUGLEVEL > 0) {
673  std::cout << "DAClusterizerInZT_vect.thermalize niter = " << niter << " at T = " << 1 / beta
674  << " nv = " << v.getSize() << std::endl;
675  if (DEBUGLEVEL > 2)
676  dump(beta, v, tks, 0);
677  }
678 #endif
679 
680  return niter;
681 }

References funct::abs(), zMuMuMuonUserData::beta, gather_cfg::cout, dumpMFGeometry_cfg::delta, FrontierConditions_GlobalTag_cff::dump, dqmdumpme::k, SiStripPI::max, mathSSE::sqrt(), update, findQualityFiles::v, and HLTMuonOfflineAnalyzer_cfi::z0.

◆ update()

double DAClusterizerInZT_vect::update ( double  beta,
track_t &  gtracks,
vertex_t &  gvertices,
const double  rho0 = 0 
) const

Definition at line 351 of file DAClusterizerInZT_vect.cc.

351  {
352  //update weights and vertex positions
353  // mass constrained annealing without noise
354  // returns the maximum of changes of vertex positions
355 
356  const unsigned int nt = gtracks.getSize();
357  const unsigned int nv = gvertices.getSize();
358 
359  //initialize sums
360  double sumpi = 0.;
361 
362  // to return how much the prototype moved
363  double delta = 0.;
364 
365  // intial value of a sum
366  double Z_init = 0;
367  // independpent of loop
368  if (rho0 > 0) {
369  Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); // cut-off
370  }
371 
372  // define kernels
373 
374  auto kernel_calc_exp_arg_range = [beta](const unsigned int itrack,
375  track_t const& tracks,
376  vertex_t const& vertices,
377  const unsigned int kmin,
378  const unsigned int kmax) {
379  const auto track_z = tracks.z_ptr[itrack];
380  const auto track_t = tracks.t_ptr[itrack];
381  const auto botrack_dz2 = -beta * tracks.dz2_ptr[itrack];
382 #ifndef USEVTXDT2
383  const auto botrack_dt2 = -beta * tracks.dt2_ptr[itrack];
384 #endif
385  // auto-vectorized
386  for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
387  const auto mult_resz = track_z - vertices.z_ptr[ivertex];
388  const auto mult_rest = track_t - vertices.t_ptr[ivertex];
389 #ifdef USEVTXDT2
390  const auto botrack_dt2 = -beta * tracks.dt2_ptr[itrack] * vertices.dt2_ptr[ivertex] /
391  (tracks.dt2_ptr[itrack] + vertices.dt2_ptr[ivertex] + 1.e-10);
392 #endif
393  vertices.ei_cache_ptr[ivertex] = botrack_dz2 * (mult_resz * mult_resz) + botrack_dt2 * (mult_rest * mult_rest);
394  }
395  };
396 
397  auto kernel_add_Z_range = [Z_init](
398  vertex_t const& vertices, const unsigned int kmin, const unsigned int kmax) -> double {
399  double ZTemp = Z_init;
400  for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
401  ZTemp += vertices.pk_ptr[ivertex] * vertices.ei_ptr[ivertex];
402  }
403  return ZTemp;
404  };
405 
406  auto kernel_calc_normalization_range = [](const unsigned int track_num,
407  track_t& tks_vec,
408  vertex_t& y_vec,
409  const unsigned int kmin,
410  const unsigned int kmax) {
411  auto tmp_trk_pi = tks_vec.pi_ptr[track_num];
412  auto o_trk_Z_sum = 1. / tks_vec.Z_sum_ptr[track_num];
413  auto o_trk_err_z = tks_vec.dz2_ptr[track_num];
414  auto o_trk_err_t = tks_vec.dt2_ptr[track_num];
415  auto tmp_trk_z = tks_vec.z_ptr[track_num];
416  auto tmp_trk_t = tks_vec.t_ptr[track_num];
417 
418  // auto-vectorized
419  for (unsigned int k = kmin; k < kmax; ++k) {
420  // parens are important for numerical stability
421  y_vec.se_ptr[k] += tmp_trk_pi * (y_vec.ei_ptr[k] * o_trk_Z_sum);
422  const auto w = tmp_trk_pi * (y_vec.pk_ptr[k] * y_vec.ei_ptr[k] * o_trk_Z_sum); // p_{ik}
423  const auto wz = w * o_trk_err_z;
424  const auto wt = w * o_trk_err_t;
425 #ifdef USEVTXDT2
426  y_vec.sumw_ptr[k] += w; // for vtxdt2
427 #endif
428  y_vec.nuz_ptr[k] += wz;
429  y_vec.nut_ptr[k] += wt;
430  y_vec.swz_ptr[k] += wz * tmp_trk_z;
431  y_vec.swt_ptr[k] += wt * tmp_trk_t;
432  /* this is really only needed when we want to get Tc too, maybe better to do it elsewhere? */
433  const auto dsz = (tmp_trk_z - y_vec.z_ptr[k]) * o_trk_err_z;
434  const auto dst = (tmp_trk_t - y_vec.t_ptr[k]) * o_trk_err_t;
435  y_vec.szz_ptr[k] += w * dsz * dsz;
436  y_vec.stt_ptr[k] += w * dst * dst;
437  y_vec.szt_ptr[k] += w * dsz * dst;
438  }
439  };
440 
441  for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
442  gvertices.se_ptr[ivertex] = 0.0;
443  gvertices.nuz_ptr[ivertex] = 0.0;
444  gvertices.nut_ptr[ivertex] = 0.0;
445  gvertices.swz_ptr[ivertex] = 0.0;
446  gvertices.swt_ptr[ivertex] = 0.0;
447  gvertices.szz_ptr[ivertex] = 0.0;
448  gvertices.stt_ptr[ivertex] = 0.0;
449  gvertices.szt_ptr[ivertex] = 0.0;
450 #ifdef USEVTXDT2
451  gvertices.sumw_ptr[k] = 0.0;
452 #endif
453  }
454 
455  // loop over tracks
456  for (auto itrack = 0U; itrack < nt; ++itrack) {
457  unsigned int kmin = gtracks.kmin[itrack];
458  unsigned int kmax = gtracks.kmax[itrack];
459 
460 #ifdef DEBUG
461  assert((kmin < kmax) && (kmax <= nv));
462  assert(itrack < gtracks.Z_sum.size());
463 #endif
464 
465  kernel_calc_exp_arg_range(itrack, gtracks, gvertices, kmin, kmax);
466  local_exp_list_range(gvertices.ei_cache_ptr, gvertices.ei_ptr, kmin, kmax);
467 
468  gtracks.Z_sum_ptr[itrack] = kernel_add_Z_range(gvertices, kmin, kmax);
469  if (edm::isNotFinite(gtracks.Z_sum_ptr[itrack]))
470  gtracks.Z_sum_ptr[itrack] = 0.0;
471  // used in the next major loop to follow
472  sumpi += gtracks.pi_ptr[itrack];
473 
474  if (gtracks.Z_sum_ptr[itrack] > 1.e-100) {
475  kernel_calc_normalization_range(itrack, gtracks, gvertices, kmin, kmax);
476  }
477  }
478 
479  // now update z, t, and pk
480  auto kernel_calc_zt = [sumpi, nv](vertex_t& vertices) -> double {
481  double delta = 0;
482 
483  // does not vectorizes
484  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
485  if (vertices.nuz_ptr[ivertex] > 0.) {
486  auto znew = vertices.swz_ptr[ivertex] / vertices.nuz_ptr[ivertex];
487  delta = max(std::abs(vertices.z_ptr[ivertex] - znew), delta);
488  vertices.z_ptr[ivertex] = znew;
489  }
490 
491  if (vertices.nut_ptr[ivertex] > 0.) {
492  auto tnew = vertices.swt_ptr[ivertex] / vertices.nut_ptr[ivertex];
493  //delta = max(std::abs(vertices.t_ptr[ ivertex ] - tnew), delta); // FIXME
494  vertices.t_ptr[ivertex] = tnew;
495 #ifdef USEVTXDT2
496  vertices.dt2_ptr[ivertex] = vertices.nut_ptr[ivertex] / vertices.sumw_ptr[ivertex];
497 #endif
498  } else {
499  // FIXME
500  // apparently this cluster has not timing info attached
501  // this should be taken into account somehow, otherwise we might fail to attach
502  // new tracks that do have timing information
503  vertices.t_ptr[ivertex] = 0;
504 #ifdef USEVTXDT2
505  vertices.dt2_ptr[ivertex] = 0;
506 #endif
507  }
508 
509 #ifdef DEBUG
510  if ((vertices.nut_ptr[ivertex] <= 0.) && (vertices.nuz_ptr[ivertex] <= 0.)) {
511  edm::LogInfo("sumw") << "invalid sum of weights in fit: " << endl;
512  std::cout << " a cluster melted away ? "
513  << " zk=" << vertices.z_ptr[ivertex] << " pk=" << vertices.pk_ptr[ivertex]
514  << " sumw(z,t) =" << vertices.nuz_ptr[ivertex] << "," << vertices.nut_ptr[ivertex] << endl;
515  // FIXME: discard this cluster
516  }
517 #endif
518  }
519 
520  auto osumpi = 1. / sumpi;
521  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
522  vertices.pk_ptr[ivertex] = vertices.pk_ptr[ivertex] * vertices.se_ptr[ivertex] * osumpi;
523  }
524 
525  return delta;
526  };
527 
528  delta += kernel_calc_zt(gvertices);
529 
530  if (zorder(gvertices)) {
531  set_vtx_range(beta, gtracks, gvertices);
532  };
533 
534  // return how much the prototypes moved
535  return delta;
536 }

References funct::abs(), cms::cuda::assert(), zMuMuMuonUserData::beta, gather_cfg::cout, dumpMFGeometry_cfg::delta, math::cholesky::dst, DAClusterizerInZT_vect::vertex_t::ei_cache_ptr, DAClusterizerInZT_vect::vertex_t::ei_ptr, DAClusterizerInZT_vect::track_t::getSize(), DAClusterizerInZT_vect::vertex_t::getSize(), edm::isNotFinite(), dqmdumpme::k, DAClusterizerInZT_vect::track_t::kmax, DAClusterizerInZT_vect::track_t::kmin, SiStripPI::max, nt, DAClusterizerInZT_vect::vertex_t::nut_ptr, DAClusterizerInZT_vect::vertex_t::nuz_ptr, DAClusterizerInZT_vect::track_t::pi_ptr, DAClusterizerInZT_vect::vertex_t::se_ptr, DAClusterizerInZT_vect::vertex_t::stt_ptr, DAClusterizerInZT_vect::vertex_t::sumw_ptr, DAClusterizerInZT_vect::vertex_t::swt_ptr, DAClusterizerInZT_vect::vertex_t::swz_ptr, DAClusterizerInZT_vect::vertex_t::szt_ptr, DAClusterizerInZT_vect::vertex_t::szz_ptr, PDWG_EXOHSCP_cff::tracks, mitigatedMETSequence_cff::U, pwdgSkimBPark_cfi::vertices, w, DAClusterizerInZT_vect::track_t::Z_sum, and DAClusterizerInZT_vect::track_t::Z_sum_ptr.

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

◆ verify()

void DAClusterizerInZT_vect::verify ( const vertex_t &  v,
const track_t &  tks,
unsigned int  nv = 999999,
unsigned int  nt = 999999 
) const

Definition at line 121 of file DAClusterizerInZT_vect.cc.

121  {
122  if (nv == 999999) {
123  nv = v.getSize();
124  } else {
125  assert(nv == v.getSize());
126  }
127 
128  if (nt == 999999) {
129  nt = tks.getSize();
130  } else {
131  assert(nt == tks.getSize());
132  }
133 
134  // clusters
135  assert(v.z.size() == nv);
136  assert(v.t.size() == nv);
137 #ifdef USEVTXDT2
138  assert(v.dt2.size() == nv);
139 #endif
140  assert(v.sumw.size() == nv);
141  assert(v.pk.size() == nv);
142  assert(v.swz.size() == nv);
143  assert(v.swt.size() == nv);
144  assert(v.ei_cache.size() == nv);
145  assert(v.ei.size() == nv);
146  assert(v.se.size() == nv);
147  assert(v.nuz.size() == nv);
148  assert(v.nut.size() == nv);
149  assert(v.szz.size() == nv);
150  assert(v.stt.size() == nv);
151  assert(v.szt.size() == nv);
152 
153  assert(v.z_ptr == &v.z.front());
154  assert(v.t_ptr == &v.t.front());
155  assert(v.pk_ptr == &v.pk.front());
156  assert(v.ei_cache_ptr == &v.ei_cache.front());
157  assert(v.swz_ptr == &v.swz.front());
158  assert(v.swt_ptr == &v.swt.front());
159  assert(v.se_ptr == &v.se.front());
160  assert(v.nuz_ptr == &v.nuz.front());
161  assert(v.nut_ptr == &v.nut.front());
162  assert(v.szz_ptr == &v.szz.front());
163  assert(v.stt_ptr == &v.stt.front());
164  assert(v.szt_ptr == &v.szt.front());
165  assert(v.sumw_ptr == &v.sumw.front());
166 
167 #ifdef USEVTXDT2
168  assert(v.dt2_ptr == &v.dt2.front());
169 #endif
170 
171  for (unsigned int k = 0; k < nv - 1; k++) {
172  if (v.z[k] <= v.z[k + 1])
173  continue;
174  cout << " ZT, cluster z-ordering assertion failure z[" << k << "] =" << v.z[k] << " z[" << k + 1
175  << "] =" << v.z[k + 1] << endl;
176  }
177  //for(unsigned int k=0; k< nv-1; k++){
178  //assert( v.z[k] <= v.z[k+1]);
179  //}
180 
181  // tracks
182  assert(nt == tks.z.size());
183  assert(nt == tks.t.size());
184  assert(nt == tks.dz2.size());
185 #ifdef USEVTXDT2
186  assert(nt == tks.dt2.size());
187 #endif
188  assert(nt == tks.tt.size());
189  assert(nt == tks.pi.size());
190  assert(nt == tks.Z_sum.size());
191  assert(nt == tks.kmin.size());
192  assert(nt == tks.kmax.size());
193 
194  assert(tks.z_ptr == &tks.z.front());
195  assert(tks.t_ptr == &tks.t.front());
196  assert(tks.dz2_ptr == &tks.dz2.front());
197  assert(tks.dt2_ptr == &tks.dt2.front());
198  assert(tks.pi_ptr == &tks.pi.front());
199  assert(tks.Z_sum_ptr == &tks.Z_sum.front());
200 
201  for (unsigned int i = 0; i < nt; i++) {
202  if ((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv))
203  continue;
204  cout << "track vertex range assertion failure" << i << "/" << nt << " kmin,kmax=" << tks.kmin[i] << ", "
205  << tks.kmax[i] << " nv=" << nv << endl;
206  }
207 
208  for (unsigned int i = 0; i < nt; i++) {
209  assert((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv));
210  }
211 }

References cms::cuda::assert(), gather_cfg::cout, DAClusterizerInZT_vect::track_t::dt2, DAClusterizerInZT_vect::track_t::dt2_ptr, DAClusterizerInZT_vect::track_t::dz2, DAClusterizerInZT_vect::track_t::dz2_ptr, DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, dqmdumpme::k, DAClusterizerInZT_vect::track_t::kmax, DAClusterizerInZT_vect::track_t::kmin, nt, DAClusterizerInZT_vect::track_t::pi, DAClusterizerInZT_vect::track_t::pi_ptr, DAClusterizerInZT_vect::track_t::t, DAClusterizerInZT_vect::track_t::t_ptr, DAClusterizerInZT_vect::track_t::tt, findQualityFiles::v, DAClusterizerInZT_vect::track_t::z, DAClusterizerInZT_vect::track_t::z_ptr, DAClusterizerInZT_vect::track_t::Z_sum, and DAClusterizerInZT_vect::track_t::Z_sum_ptr.

◆ vertices()

vector< TransientVertex > DAClusterizerInZT_vect::vertices ( const std::vector< reco::TransientTrack > &  tracks,
const int  verbosity = 0 
) const

Definition at line 1096 of file DAClusterizerInZT_vect.cc.

1097  {
1098  track_t&& tks = fill(tracks);
1099  tks.extractRaw();
1100 
1101  unsigned int nt = tks.getSize();
1102  double rho0 = 0.0; // start with no outlier rejection
1103 
1104  vector<TransientVertex> clusters;
1105  if (tks.getSize() == 0)
1106  return clusters;
1107 
1108  vertex_t y; // the vertex prototypes
1109 
1110  // initialize:single vertex at infinite temperature
1111  y.addItem(0, 0, 1.0);
1112  clear_vtx_range(tks, y);
1113 
1114  // estimate first critical temperature
1115  double beta = beta0(betamax_, tks, y);
1116 #ifdef DEBUG
1117  if (DEBUGLEVEL > 0)
1118  std::cout << "Beta0 is " << beta << std::endl;
1119 #endif
1120 
1121  thermalize(beta, tks, y, delta_highT_, 0.);
1122 
1123  // annealing loop, stop when T<minT (i.e. beta>1/minT)
1124 
1125  double betafreeze = betamax_ * sqrt(coolingFactor_);
1126 
1127  while (beta < betafreeze) {
1128  update(beta, tks, y, rho0);
1129  while (merge(y, tks, beta)) {
1130  if (zorder(y))
1131  set_vtx_range(beta, tks, y);
1132  update(beta, tks, y, rho0);
1133  }
1134  split(beta, tks, y);
1135 
1136  beta = beta / coolingFactor_;
1137  thermalize(beta, tks, y, delta_highT_, 0.);
1138  }
1139 
1140 #ifdef DEBUG
1141  verify(y, tks);
1142 
1143  if (DEBUGLEVEL > 0) {
1144  std::cout << "DAClusterizerInZT_vect::vertices :"
1145  << "merging at T=" << 1 / beta << std::endl;
1146  }
1147 #endif
1148 
1149  zorder(y);
1150  set_vtx_range(beta, tks, y);
1151  update(beta, tks, y, rho0); // make sure Tc is up-to-date
1152 
1153  while (merge(y, tks, beta)) {
1154  zorder(y);
1155  set_vtx_range(beta, tks, y);
1156  update(beta, tks, y, rho0);
1157  }
1158 
1159 #ifdef DEBUG
1160  verify(y, tks);
1161  if (DEBUGLEVEL > 0) {
1162  std::cout << "DAClusterizerInZT_vect::vertices :"
1163  << "splitting/merging at T=" << 1 / beta << std::endl;
1164  }
1165 #endif
1166 
1167  unsigned int ntry = 0;
1168  double threshold = 1.0;
1169  while (split(beta, tks, y, threshold) && (ntry++ < 10)) {
1170  thermalize(beta, tks, y, delta_highT_, 0.);
1171 
1172  while (merge(y, tks, beta)) {
1173  if (zorder(y))
1174  set_vtx_range(beta, tks, y);
1175  update(beta, tks, y, rho0);
1176  }
1177 
1178 #ifdef DEBUG
1179  verify(y, tks);
1180 #endif
1181 
1182  // relax splitting a bit to reduce multiple split-merge cycles of the same cluster
1183  threshold *= 1.1;
1184  }
1185 
1186 #ifdef DEBUG
1187  verify(y, tks);
1188  if (DEBUGLEVEL > 0) {
1189  std::cout << "DAClusterizerInZT_vect::vertices :"
1190  << "turning on outlier rejection at T=" << 1 / beta << std::endl;
1191  }
1192 #endif
1193 
1194  // switch on outlier rejection at T=minT
1195  if (dzCutOff_ > 0) {
1196  rho0 = 1. / nt;
1197  for (unsigned int a = 0; a < 5; a++) {
1198  update(beta, tks, y, a * rho0 / 5); // adiabatic turn-on
1199  }
1200  }
1201 
1202  thermalize(beta, tks, y, delta_lowT_, rho0);
1203 
1204 #ifdef DEBUG
1205  verify(y, tks);
1206  if (DEBUGLEVEL > 0) {
1207  std::cout << "DAClusterizerInZT_vect::vertices :"
1208  << "merging with outlier rejection at T=" << 1 / beta << std::endl;
1209  }
1210  if (DEBUGLEVEL > 2)
1211  dump(beta, y, tks, 2);
1212 #endif
1213 
1214  // merge again (some cluster split by outliers collapse here)
1215  zorder(y);
1216  while (merge(y, tks, beta)) {
1217  set_vtx_range(beta, tks, y);
1218  update(beta, tks, y, rho0);
1219  }
1220 
1221 #ifdef DEBUG
1222  if (DEBUGLEVEL > 0) {
1223  std::cout << "DAClusterizerInZT_vect::vertices :"
1224  << "after merging with outlier rejection at T=" << 1 / beta << std::endl;
1225  }
1226  if (DEBUGLEVEL > 1)
1227  dump(beta, y, tks, 2);
1228 #endif
1229 
1230  // go down to the purging temperature (if it is lower than tmin)
1231  while (beta < betapurge_) {
1232  beta = min(beta / coolingFactor_, betapurge_);
1233  thermalize(beta, tks, y, delta_lowT_, rho0);
1234  }
1235 
1236 #ifdef DEBUG
1237  verify(y, tks);
1238  if (DEBUGLEVEL > 0) {
1239  std::cout << "DAClusterizerInZT_vect::vertices :"
1240  << "purging at T=" << 1 / beta << std::endl;
1241  }
1242 #endif
1243 
1244  // eliminate insigificant vertices, this is more restrictive at higher T
1245  while (purge(y, tks, rho0, beta)) {
1246  thermalize(beta, tks, y, delta_lowT_, rho0);
1247  if (zorder(y)) {
1248  set_vtx_range(beta, tks, y);
1249  };
1250  }
1251 
1252 #ifdef DEBUG
1253  verify(y, tks);
1254  if (DEBUGLEVEL > 0) {
1255  std::cout << "DAClusterizerInZT_vect::vertices :"
1256  << "last cooling T=" << 1 / beta << std::endl;
1257  }
1258 #endif
1259 
1260  // optionally cool some more without doing anything, to make the assignment harder
1261  while (beta < betastop_) {
1262  beta = min(beta / coolingFactor_, betastop_);
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  << "stop cooling at T=" << 1 / beta << std::endl;
1274  }
1275  if (DEBUGLEVEL > 2)
1276  dump(beta, y, tks, 2);
1277 #endif
1278 
1279  // new, merge here and not in "clusterize"
1280  // final merging step
1281  double betadummy = 1;
1282  while (merge(y, tks, betadummy))
1283  ;
1284 
1285  // select significant tracks and use a TransientVertex as a container
1286  GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01);
1287 
1288  // ensure correct normalization of probabilities, should makes double assignment reasonably impossible
1289  const unsigned int nv = y.getSize();
1290  for (unsigned int k = 0; k < nv; k++)
1291  if (edm::isNotFinite(y.pk_ptr[k]) || edm::isNotFinite(y.z_ptr[k])) {
1292  y.pk_ptr[k] = 0;
1293  y.z_ptr[k] = 0;
1294  }
1295 
1296  for (unsigned int i = 0; i < nt; i++) // initialize
1297  tks.Z_sum_ptr[i] = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
1298 
1299  // improve vectorization (does not require reduction ....)
1300  for (unsigned int k = 0; k < nv; k++) {
1301  for (unsigned int i = 0; i < nt; i++)
1302  tks.Z_sum_ptr[i] +=
1303  y.pk_ptr[k] *
1304  local_exp(-beta * Eik(tks.z_ptr[i], y.z_ptr[k], tks.dz2_ptr[i], tks.t_ptr[i], y.t_ptr[k], tks.dt2_ptr[i]));
1305  }
1306 
1307  for (unsigned int k = 0; k < nv; k++) {
1308  GlobalPoint pos(0, 0, y.z_ptr[k]);
1309 
1310  vector<reco::TransientTrack> vertexTracks;
1311  for (unsigned int i = 0; i < nt; i++) {
1312  if (tks.Z_sum_ptr[i] > 1e-100) {
1313  double p =
1314  y.pk_ptr[k] *
1315  local_exp(-beta * Eik(tks.z_ptr[i], y.z_ptr[k], tks.dz2_ptr[i], tks.t_ptr[i], y.t_ptr[k], tks.dt2_ptr[i])) /
1316  tks.Z_sum_ptr[i];
1317  if ((tks.pi_ptr[i] > 0) && (p > mintrkweight_)) {
1318  vertexTracks.push_back(*(tks.tt[i]));
1319  tks.Z_sum_ptr[i] = 0; // setting Z=0 excludes double assignment
1320  }
1321  }
1322  }
1323  TransientVertex v(pos, y.t_ptr[k], dummyError, vertexTracks, 0);
1324  clusters.push_back(v);
1325  }
1326 
1327  return clusters;
1328 }

References a, zMuMuMuonUserData::beta, bsc_activity_cfg::clusters, gather_cfg::cout, DAClusterizerInZT_vect::track_t::dt2_ptr, FrontierConditions_GlobalTag_cff::dump, DAClusterizerInZT_vect::track_t::dz2_ptr, MillePedeFileConverter_cfg::e, DAClusterizerInZT_vect::track_t::extractRaw(), ntuplemaker::fill, DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, edm::isNotFinite(), dqmdumpme::k, MatrixUtil::merge(), min(), nt, AlCaHLTBitMon_ParallelJobs::p, DAClusterizerInZT_vect::track_t::pi_ptr, cms::dd::split(), mathSSE::sqrt(), DAClusterizerInZT_vect::track_t::t_ptr, remoteMonitoring_LED_IterMethod_cfg::threshold, PDWG_EXOHSCP_cff::tracks, DAClusterizerInZT_vect::track_t::tt, update, findQualityFiles::v, DAClusterizerInZT_vect::track_t::z_ptr, and DAClusterizerInZT_vect::track_t::Z_sum_ptr.

◆ zorder()

bool DAClusterizerInZT_vect::zorder ( vertex_t &  y) const

Definition at line 538 of file DAClusterizerInZT_vect.cc.

538  {
539  const unsigned int nv = y.getSize();
540 
541 #ifdef DEBUG
542  assert(y.z.size() == nv);
543  assert(y.t.size() == nv);
544  // assert(y.dt2.size() == nv);
545  assert(y.pk.size() == nv);
546 #endif
547 
548  if (nv < 2)
549  return false;
550 
551  bool reordering = true;
552  bool modified = false;
553 
554  while (reordering) {
555  reordering = false;
556  for (unsigned int k = 0; (k + 1) < nv; k++) {
557  if (y.z[k + 1] < y.z[k]) {
558  auto ztemp = y.z[k];
559  y.z[k] = y.z[k + 1];
560  y.z[k + 1] = ztemp;
561  auto ttemp = y.t[k];
562  y.t[k] = y.t[k + 1];
563  y.t[k + 1] = ttemp;
564 #ifdef USEVTXDT2
565  auto dt2temp = y.dt2[k];
566  y.dt2[k] = y.dt2[k + 1];
567  y.dt2[k + 1] = dt2temp;
568 #endif
569  auto ptemp = y.pk[k];
570  y.pk[k] = y.pk[k + 1];
571  y.pk[k + 1] = ptemp;
572  reordering = true;
573  }
574  }
575  modified |= reordering;
576  }
577 
578  if (modified) {
579  y.extractRaw();
580  return true;
581  }
582 
583  return false;
584 }

References cms::cuda::assert(), and dqmdumpme::k.

Member Data Documentation

◆ betamax_

double DAClusterizerInZT_vect::betamax_
private

Definition at line 278 of file DAClusterizerInZT_vect.h.

◆ betapurge_

double DAClusterizerInZT_vect::betapurge_
private

Definition at line 289 of file DAClusterizerInZT_vect.h.

◆ betastop_

double DAClusterizerInZT_vect::betastop_
private

Definition at line 279 of file DAClusterizerInZT_vect.h.

◆ convergence_mode_

unsigned int DAClusterizerInZT_vect::convergence_mode_
private

Definition at line 291 of file DAClusterizerInZT_vect.h.

◆ coolingFactor_

double DAClusterizerInZT_vect::coolingFactor_
private

Definition at line 277 of file DAClusterizerInZT_vect.h.

◆ d0CutOff_

double DAClusterizerInZT_vect::d0CutOff_
private

Definition at line 281 of file DAClusterizerInZT_vect.h.

◆ delta_highT_

double DAClusterizerInZT_vect::delta_highT_
private

Definition at line 292 of file DAClusterizerInZT_vect.h.

◆ delta_lowT_

double DAClusterizerInZT_vect::delta_lowT_
private

Definition at line 293 of file DAClusterizerInZT_vect.h.

◆ dtCutOff_

double DAClusterizerInZT_vect::dtCutOff_
private

Definition at line 282 of file DAClusterizerInZT_vect.h.

◆ dzCutOff_

double DAClusterizerInZT_vect::dzCutOff_
private

Definition at line 280 of file DAClusterizerInZT_vect.h.

◆ maxIterations_

unsigned int DAClusterizerInZT_vect::maxIterations_
private

Definition at line 276 of file DAClusterizerInZT_vect.h.

◆ mintrkweight_

double DAClusterizerInZT_vect::mintrkweight_
private

Definition at line 285 of file DAClusterizerInZT_vect.h.

◆ sel_zrange_

double DAClusterizerInZT_vect::sel_zrange_
private

Definition at line 295 of file DAClusterizerInZT_vect.h.

◆ t0Max_

double DAClusterizerInZT_vect::t0Max_
private

Definition at line 283 of file DAClusterizerInZT_vect.h.

◆ tmerge_

double DAClusterizerInZT_vect::tmerge_
private

Definition at line 288 of file DAClusterizerInZT_vect.h.

◆ uniquetrkweight_

double DAClusterizerInZT_vect::uniquetrkweight_
private

Definition at line 286 of file DAClusterizerInZT_vect.h.

◆ verbose_

bool DAClusterizerInZT_vect::verbose_
private

Definition at line 270 of file DAClusterizerInZT_vect.h.

◆ vertexSize_

double DAClusterizerInZT_vect::vertexSize_
private

Definition at line 274 of file DAClusterizerInZT_vect.h.

◆ vertexSizeTime_

double DAClusterizerInZT_vect::vertexSizeTime_
private

Definition at line 275 of file DAClusterizerInZT_vect.h.

◆ zdumpcenter_

double DAClusterizerInZT_vect::zdumpcenter_
private

Definition at line 271 of file DAClusterizerInZT_vect.h.

◆ zdumpwidth_

double DAClusterizerInZT_vect::zdumpwidth_
private

Definition at line 272 of file DAClusterizerInZT_vect.h.

◆ zmerge_

double DAClusterizerInZT_vect::zmerge_
private

Definition at line 287 of file DAClusterizerInZT_vect.h.

◆ zrange_min_

const double DAClusterizerInZT_vect::zrange_min_ = 0.1
private

Definition at line 296 of file DAClusterizerInZT_vect.h.

DAClusterizerInZT_vect::coolingFactor_
double coolingFactor_
Definition: DAClusterizerInZT_vect.h:277
HIPAlignmentAlgorithm_cfi.verbosity
verbosity
Definition: HIPAlignmentAlgorithm_cfi.py:7
RandomServiceHelper.t2
t2
Definition: RandomServiceHelper.py:257
reco::HitPattern::MISSING_OUTER_HITS
Definition: HitPattern.h:155
DDAxes::y
alignBH_cfg.fixed
fixed
Definition: alignBH_cfg.py:54
PDWG_EXOHSCP_cff.tracks
tracks
Definition: PDWG_EXOHSCP_cff.py:28
Measurement1D
Definition: Measurement1D.h:11
mps_fire.i
i
Definition: mps_fire.py:355
DAClusterizerInZT_vect::uniquetrkweight_
double uniquetrkweight_
Definition: DAClusterizerInZT_vect.h:286
TkClusParameters_cff.zrange
zrange
Definition: TkClusParameters_cff.py:7
reco::ParticleMasses::k0
const double k0
Definition: ParticleMasses.h:7
edm::isNotFinite
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
nt
int nt
Definition: AMPTWrapper.h:42
DAClusterizerInZT_vect::convergence_mode_
unsigned int convergence_mode_
Definition: DAClusterizerInZT_vect.h:291
min
T min(T a, T b)
Definition: MathUtil.h:58
DAClusterizerInZT_vect::vertices
std::vector< TransientVertex > vertices(const std::vector< reco::TransientTrack > &tracks, const int verbosity=0) const
Definition: DAClusterizerInZT_vect.cc:1096
DAClusterizerInZT_vect::delta_highT_
double delta_highT_
Definition: DAClusterizerInZT_vect.h:292
zMuMuMuonUserData.beta
beta
Definition: zMuMuMuonUserData.py:10
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
gather_cfg.cout
cout
Definition: gather_cfg.py:144
DAClusterizerInZT_vect::betastop_
double betastop_
Definition: DAClusterizerInZT_vect.h:279
pos
Definition: PixelAliasList.h:18
AlignmentPI::t_z
Definition: AlignmentPayloadInspectorHelper.h:36
edm::LogInfo
Definition: MessageLogger.h:254
cms::cuda::assert
assert(be >=bs)
DAClusterizerInZT_vect::delta_lowT_
double delta_lowT_
Definition: DAClusterizerInZT_vect.h:293
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:215
DAClusterizerInZT_vect::zdumpcenter_
double zdumpcenter_
Definition: DAClusterizerInZT_vect.h:271
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
DAClusterizerInZT_vect::purge
bool purge(vertex_t &, track_t &, double &, const double) const
Definition: DAClusterizerInZT_vect.cc:742
DAClusterizerInZT_vect::verbose_
bool verbose_
Definition: DAClusterizerInZT_vect.h:270
findQualityFiles.v
v
Definition: findQualityFiles.py:179
SiStripMonitorCluster_cfi.zmin
zmin
Definition: SiStripMonitorCluster_cfi.py:200
DAClusterizerInZT_vect::sel_zrange_
double sel_zrange_
Definition: DAClusterizerInZT_vect.h:295
DAClusterizerInZT_vect::get_Tc
double get_Tc(const vertex_t &y, int k) const
Definition: DAClusterizerInZT_vect.cc:880
geometryDiff.epsilon
int epsilon
Definition: geometryDiff.py:26
F
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:163
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
testProducerWithPsetDescEmpty_cfi.z2
z2
Definition: testProducerWithPsetDescEmpty_cfi.py:41
DAClusterizerInZT_vect::vertexSize_
double vertexSize_
Definition: DAClusterizerInZT_vect.h:274
dt
float dt
Definition: AMPTWrapper.h:136
GeomDetEnumerators::PixelBarrel
Definition: GeomDetEnumerators.h:11
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
DAClusterizerInZT_vect::split
bool split(const double beta, track_t &t, vertex_t &y, double threshold=1.) const
Definition: DAClusterizerInZT_vect.cc:891
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
h
SiStripMonitorCluster_cfi.zmax
zmax
Definition: SiStripMonitorCluster_cfi.py:201
DAClusterizerInZT_vect::zdumpwidth_
double zdumpwidth_
Definition: DAClusterizerInZT_vect.h:272
w
const double w
Definition: UKUtility.cc:23
DAClusterizerInZT_vect::t0Max_
double t0Max_
Definition: DAClusterizerInZT_vect.h:283
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
reco::BeamSpot
Definition: BeamSpot.h:21
DDAxes::z
DAClusterizerInZT_vect::vertexSizeTime_
double vertexSizeTime_
Definition: DAClusterizerInZT_vect.h:275
DAClusterizerInZT_vect::betamax_
double betamax_
Definition: DAClusterizerInZT_vect.h:278
p2
double p2[4]
Definition: TauolaWrapper.h:90
beamspot
Definition: BeamSpotWrite2Txt.h:8
HLTMuonOfflineAnalyzer_cfi.z0
z0
Definition: HLTMuonOfflineAnalyzer_cfi.py:98
dqmdumpme.k
k
Definition: dqmdumpme.py:60
DAClusterizerInZT_vect::beta0
double beta0(const double betamax, track_t const &tks, vertex_t const &y) const
Definition: DAClusterizerInZT_vect.cc:813
Point3DBase< float, GlobalTag >
OrderedSet.t
t
Definition: OrderedSet.py:90
b
double b
Definition: hdecay.h:118
DAClusterizerInZT_vect::zorder
bool zorder(vertex_t &y) const
Definition: DAClusterizerInZT_vect.cc:538
mitigatedMETSequence_cff.U
U
Definition: mitigatedMETSequence_cff.py:36
DDAxes::rho
edm::LogWarning
Definition: MessageLogger.h:141
bsc_activity_cfg.clusters
clusters
Definition: bsc_activity_cfg.py:36
listHistos.IP
IP
Definition: listHistos.py:76
DAClusterizerInZT_vect::set_vtx_range
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const
Definition: DAClusterizerInZT_vect.cc:285
DAClusterizerInZT_vect::fill
track_t fill(const std::vector< reco::TransientTrack > &tracks) const
Definition: DAClusterizerInZT_vect.cc:214
q1
double q1[4]
Definition: TauolaWrapper.h:87
a
double a
Definition: hdecay.h:119
DAClusterizerInZT_vect::verify
void verify(const vertex_t &v, const track_t &tks, unsigned int nv=999999, unsigned int nt=999999) const
Definition: DAClusterizerInZT_vect.cc:121
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
DAClusterizerInZT_vect::zmerge_
double zmerge_
Definition: DAClusterizerInZT_vect.h:287
DAClusterizerInZT_vect::zrange_min_
const double zrange_min_
Definition: DAClusterizerInZT_vect.h:296
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
DAClusterizerInZT_vect::merge
bool merge(vertex_t &, track_t &, double &beta) const
Definition: DAClusterizerInZT_vect.cc:683
DAClusterizerInZT_vect::dzCutOff_
double dzCutOff_
Definition: DAClusterizerInZT_vect.h:280
createfilelist.int
int
Definition: createfilelist.py:10
MetAnalyzer.pv
def pv(vc)
Definition: MetAnalyzer.py:7
hlt_jetmet_dqm_QT_fromfile_cfg.critical
critical
Definition: hlt_jetmet_dqm_QT_fromfile_cfg.py:109
p1
double p1[4]
Definition: TauolaWrapper.h:89
GlobalErrorBase< double, ErrorMatrixTag >
TransientVertex
Definition: TransientVertex.h:18
DAClusterizerInZT_vect::mintrkweight_
double mintrkweight_
Definition: DAClusterizerInZT_vect.h:285
DAClusterizerInZT_vect::dump
void dump(const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0) const
Definition: DAClusterizerInZT_vect.cc:1359
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
funct::D
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
DAClusterizerInZT_vect::update
double update(double beta, track_t &gtracks, vertex_t &gvertices, const double rho0=0) const
Definition: DAClusterizerInZT_vect.cc:351
PVValHelper::dz
Definition: PVValidationHelpers.h:50
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
funct::arg
A arg
Definition: Factorize.h:36
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:30
DAClusterizerInZT_vect::clear_vtx_range
void clear_vtx_range(track_t &gtracks, vertex_t &gvertices) const
Definition: DAClusterizerInZT_vect.cc:342
math::cholesky::dst
M2 & dst
Definition: choleskyInversion.h:158
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
DAClusterizerInZT_vect::tmerge_
double tmerge_
Definition: DAClusterizerInZT_vect.h:288
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
remoteMonitoring_LED_IterMethod_cfg.threshold
threshold
Definition: remoteMonitoring_LED_IterMethod_cfg.py:426
DAClusterizerInZT_vect::thermalize
unsigned int thermalize(double beta, track_t &gtracks, vertex_t &gvertices, const double delta_max, const double rho0=0.) const
Definition: DAClusterizerInZT_vect.cc:632
HLT_2018_cff.atIP
atIP
Definition: HLT_2018_cff.py:44407
TauDecayModes.dec
dec
Definition: TauDecayModes.py:143
DAClusterizerInZT_vect::dtCutOff_
double dtCutOff_
Definition: DAClusterizerInZT_vect.h:282
DAClusterizerInZT_vect::find_nearest
bool find_nearest(double z, double t, vertex_t &y, unsigned int &k_min, double dz, double dt) const
Definition: DAClusterizerInZT_vect.cc:586
DAClusterizerInZT_vect::d0CutOff_
double d0CutOff_
Definition: DAClusterizerInZT_vect.h:281
DAClusterizerInZT_vect::maxIterations_
unsigned int maxIterations_
Definition: DAClusterizerInZT_vect.h:276
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
DAClusterizerInZT_vect::betapurge_
double betapurge_
Definition: DAClusterizerInZT_vect.h:289
reco::TrackBase::highPurity
Definition: TrackBase.h:154
Generated for CMSSW Reference Manual by doxygen 1.8.16