CMS 3D CMS Logo

List of all members | Classes | Public Member Functions | Static 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 bool updateTc=false) const
 
void verify (const vertex_t &v, const track_t &tks, unsigned int nv=999999, unsigned int nt=999999) const
 
std::vector< TransientVertexvertices (const std::vector< reco::TransientTrack > &tracks) const override
 
bool zorder (vertex_t &y) const
 
- Public Member Functions inherited from TrackClusterizerInZ
 TrackClusterizerInZ ()=default
 
 TrackClusterizerInZ (const edm::ParameterSet &conf)
 
virtual ~TrackClusterizerInZ ()=default
 

Static Public Member Functions

static void fillPSetDescription (edm::ParameterSetDescription &desc)
 

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 uniquetrkminp_
 
double uniquetrkweight_
 
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 26 of file DAClusterizerInZT_vect.h.

Constructor & Destructor Documentation

◆ DAClusterizerInZT_vect()

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

Definition at line 21 of file DAClusterizerInZT_vect.cc.

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

21  {
22  // hardcoded parameters
23  maxIterations_ = 1000;
24  mintrkweight_ = 0.5;
25 
26  // configurable debug output
27 #ifdef DEBUG
28  zdumpcenter_ = conf.getUntrackedParameter<double>("zdumpcenter", 0.);
29  zdumpwidth_ = conf.getUntrackedParameter<double>("zdumpwidth", 20.);
30 #endif
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  uniquetrkminp_ = conf.getParameter<double>("uniquetrkminp");
45  zmerge_ = conf.getParameter<double>("zmerge");
46  tmerge_ = conf.getParameter<double>("tmerge");
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 #ifdef DEBUG
53  std::cout << "DAClusterizerInZT_vect: mintrkweight = " << mintrkweight_ << std::endl;
54  std::cout << "DAClusterizerInZT_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl;
55  std::cout << "DAClusterizerInZT_vect: uniquetrkminp = " << uniquetrkminp_ << std::endl;
56  std::cout << "DAClusterizerInZT_vect: zmerge = " << zmerge_ << std::endl;
57  std::cout << "DAClusterizerInZT_vect: tmerge = " << tmerge_ << std::endl;
58  std::cout << "DAClusterizerInZT_vect: Tmin = " << minT << std::endl;
59  std::cout << "DAClusterizerInZT_vect: Tpurge = " << purgeT << std::endl;
60  std::cout << "DAClusterizerInZT_vect: Tstop = " << stopT << std::endl;
61  std::cout << "DAClusterizerInZT_vect: vertexSize = " << vertexSize_ << std::endl;
62  std::cout << "DAClusterizerInZT_vect: vertexSizeTime = " << vertexSizeTime_ << std::endl;
63  std::cout << "DAClusterizerInZT_vect: coolingFactor = " << coolingFactor_ << std::endl;
64  std::cout << "DAClusterizerInZT_vect: d0CutOff = " << d0CutOff_ << std::endl;
65  std::cout << "DAClusterizerInZT_vect: dzCutOff = " << dzCutOff_ << std::endl;
66  std::cout << "DAClusterizerInZT_vect: dtCutoff = " << dtCutOff_ << std::endl;
67  std::cout << "DAClusterizerInZT_vect: zrange = " << sel_zrange_ << std::endl;
68  std::cout << "DAClusterizerinZT_vect: convergence mode = " << convergence_mode_ << std::endl;
69  std::cout << "DAClusterizerinZT_vect: delta_highT = " << delta_highT_ << std::endl;
70  std::cout << "DAClusterizerinZT_vect: delta_lowT = " << delta_lowT_ << std::endl;
71  std::cout << "DAClusterizerinZT_vect: DEBUGLEVEL " << DEBUGLEVEL << std::endl;
72 #endif
73 
74  if (convergence_mode_ > 1) {
75  edm::LogWarning("DAClusterizerinZT_vect")
76  << "DAClusterizerInZT_vect: invalid convergence_mode " << convergence_mode_ << " reset to default " << 0;
78  }
79 
80  if (minT == 0) {
81  betamax_ = 1.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 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
T getUntrackedParameter(std::string const &, T const &) const
Log< level::Warning, false > LogWarning

Member Function Documentation

◆ beta0()

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

Definition at line 832 of file DAClusterizerInZT_vect.cc.

References gather_cfg::cout, dt, DAClusterizerInZT_vect::track_t::dt2, PVValHelper::dz, DAClusterizerInZT_vect::track_t::dz2, DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, createfilelist::int, dqmdumpme::k, dqm-mbProfile::log, nt, funct::pow(), mathSSE::sqrt(), DAClusterizerInZT_vect::track_t::tkwt, DAClusterizerInZT_vect::track_t::tpca, w(), and DAClusterizerInZT_vect::track_t::zpca.

832  {
833  double T0 = 0; // max Tc for beta=0
834  // estimate critical temperature from beta=0 (T=inf)
835  const unsigned int nt = tks.getSize();
836  const unsigned int nv = y.getSize();
837 
838  for (unsigned int k = 0; k < nv; k++) {
839  // vertex fit at T=inf
840  double sumwz = 0;
841  double sumwt = 0;
842  double sumw_z = 0;
843  double sumw_t = 0;
844  for (unsigned int i = 0; i < nt; i++) {
845  double w_z = tks.tkwt[i] * tks.dz2[i];
846  double w_t = tks.tkwt[i] * tks.dt2[i];
847  sumwz += w_z * tks.zpca[i];
848  sumwt += w_t * tks.tpca[i];
849  sumw_z += w_z;
850  sumw_t += w_t;
851  }
852  y.zvtx[k] = sumwz / sumw_z;
853  y.tvtx[k] = sumwt / sumw_t;
854 
855  // estimate Tc, eventually do this in the same loop
856  double szz = 0, stt = 0, szt = 0;
857  double nuz = 0, nut = 0;
858  for (unsigned int i = 0; i < nt; i++) {
859  double dz = (tks.zpca[i] - y.zvtx[k]) * tks.dz2[i];
860  double dt = (tks.tpca[i] - y.tvtx[k]) * tks.dt2[i];
861  double w = tks.tkwt[i];
862  szz += w * dz * dz;
863  stt += w * dt * dt;
864  szt += w * dz * dt;
865  nuz += w * tks.dz2[i];
866  nut += w * tks.dt2[i];
867  }
868  double Tz = szz / nuz;
869  double Tt = 0;
870  double Tc = 0;
871  if (nut > 0) {
872  Tt = stt / nut;
873  Tc = Tz + Tt + sqrt(pow(Tz - Tt, 2) + 4 * szt * szt / nuz / nut);
874  } else {
875  Tc = 2. * Tz;
876  }
877 
878  if (Tc > T0)
879  T0 = Tc;
880  } // vertex loop (normally there should be only one vertex at beta=0)
881 
882 #ifdef DEBUG
883  if (DEBUGLEVEL > 0) {
884  std::cout << "DAClusterizerInZT_vect.beta0: Tc = " << T0 << std::endl;
885  int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
886  std::cout << "DAClusterizerInZT_vect.beta0: nstep = " << coolingsteps << std::endl;
887  }
888 #endif
889 
890  if (T0 > 1. / betamax) {
891  int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
892  return betamax * std::pow(coolingFactor_, coolingsteps);
893  } else {
894  // ensure at least one annealing step
895  return betamax * coolingFactor_;
896  }
897 }
float dt
Definition: AMPTWrapper.h:136
T w() const
T sqrt(T t)
Definition: SSEVec.h:23
int nt
Definition: AMPTWrapper.h:42
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ clear_vtx_range()

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

Definition at line 344 of file DAClusterizerInZT_vect.cc.

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.

344  {
345  const unsigned int nt = gtracks.getSize();
346  const unsigned int nv = gvertices.getSize();
347  for (auto itrack = 0U; itrack < nt; ++itrack) {
348  gtracks.kmin[itrack] = 0;
349  gtracks.kmax[itrack] = nv;
350  }
351 }
int nt
Definition: AMPTWrapper.h:42

◆ clusterize()

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

Implements TrackClusterizerInZ.

Definition at line 1356 of file DAClusterizerInZT_vect.cc.

References bsc_activity_cfg::clusters, gather_cfg::cout, dqmdumpme::k, DiMuonV_cfg::tracks, and AlignmentTracksFromVertexSelector_cfi::vertices.

1357  {
1358  vector<vector<reco::TransientTrack> > clusters;
1359  vector<TransientVertex>&& pv = vertices(tracks);
1360 
1361 #ifdef DEBUG
1362  if (DEBUGLEVEL > 0) {
1363  std::cout << "###################################################" << endl;
1364  std::cout << "# vectorized DAClusterizerInZT_vect::clusterize nt=" << tracks.size() << endl;
1365  std::cout << "# DAClusterizerInZT_vect::clusterize pv.size=" << pv.size() << endl;
1366  std::cout << "###################################################" << endl;
1367  }
1368 #endif
1369 
1370  if (pv.empty()) {
1371  return clusters;
1372  }
1373 
1374  // fill into clusters, don't merge
1375  for (auto k = pv.begin(); k != pv.end(); k++) {
1376  vector<reco::TransientTrack> aCluster = k->originalTracks();
1377  if (aCluster.size() > 1) {
1378  clusters.push_back(aCluster);
1379  }
1380  }
1381 
1382  return clusters;
1383 }
std::vector< TransientVertex > vertices(const std::vector< reco::TransientTrack > &tracks) const override

◆ dump()

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

Definition at line 1385 of file DAClusterizerInZT_vect.cc.

References a, b, HLT_2024v14_cff::beta, gather_cfg::cout, TauDecayModes::dec, DAClusterizerInZT_vect::track_t::dt2, DAClusterizerInZT_vect::track_t::dz2, F(), alignBH_cfg::fixed, DAClusterizerInZT_vect::track_t::getSize(), h, 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, GeomDetEnumerators::PixelBarrel, jetUpdater_cfi::sort, mathSSE::sqrt(), DAClusterizerInZT_vect::track_t::sum_Z, DAClusterizerInZT_vect::track_t::tkwt, DAClusterizerInZT_vect::track_t::tpca, DAClusterizerInZT_vect::track_t::tt, verbosity, and DAClusterizerInZT_vect::track_t::zpca.

1385  {
1386 #ifdef DEBUG
1387  const unsigned int nv = y.getSize();
1388  const unsigned int nt = tks.getSize();
1389 
1390  std::vector<unsigned int> iz;
1391  for (unsigned int j = 0; j < nt; j++) {
1392  iz.push_back(j);
1393  }
1394  std::sort(iz.begin(), iz.end(), [tks](unsigned int a, unsigned int b) { return tks.zpca[a] < tks.zpca[b]; });
1395  std::cout << std::endl;
1396  std::cout << "-----DAClusterizerInZT::dump ----" << nv << " clusters " << std::endl;
1397  string h = " ";
1398  std::cout << h << " k= ";
1399  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1400  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1401  std::cout << setw(8) << fixed << ivertex;
1402  }
1403  }
1404  std::cout << endl;
1405 
1406  std::cout << h << " z= ";
1407  std::cout << setprecision(4);
1408  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1409  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1410  std::cout << setw(8) << fixed << y.zvtx[ivertex];
1411  }
1412  }
1413  std::cout << endl;
1414 
1415  std::cout << h << " t= ";
1416  std::cout << setprecision(4);
1417  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1418  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1419  std::cout << setw(8) << fixed << y.tvtx[ivertex];
1420  }
1421  }
1422  std::cout << endl;
1423 
1424  std::cout << "T=" << setw(15) << 1. / beta << " Tmin =" << setw(10) << 1. / betamax_
1425  << " Tc= ";
1426  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1427  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1428  double Tc = get_Tc(y, ivertex);
1429  std::cout << setw(8) << fixed << setprecision(1) << Tc;
1430  }
1431  }
1432  std::cout << endl;
1433 
1434  std::cout << h << "rho= ";
1435  double sumrho = 0;
1436  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1437  sumrho += y.rho[ivertex];
1438  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1439  continue;
1440  std::cout << setw(8) << setprecision(4) << fixed << y.rho[ivertex];
1441  }
1442  std::cout << endl;
1443 
1444  std::cout << h << "nt= ";
1445  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1446  sumrho += y.rho[ivertex];
1447  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1448  continue;
1449  std::cout << setw(8) << setprecision(1) << fixed << y.rho[ivertex] * nt;
1450  }
1451  std::cout << endl;
1452 
1453  if (verbosity > 0) {
1454  double E = 0, F = 0;
1455  std::cout << endl;
1456  std::cout << "---- z +/- dz t +/- dt ip +/-dip pt phi eta weights ----"
1457  << endl;
1458  std::cout << setprecision(4);
1459  for (unsigned int i0 = 0; i0 < nt; i0++) {
1460  unsigned int i = iz[i0];
1461  if (tks.sum_Z[i] > 0) {
1462  F -= std::log(tks.sum_Z[i]) / beta;
1463  }
1464  double tz = tks.zpca[i];
1465 
1466  if (std::fabs(tz - zdumpcenter_) > zdumpwidth_)
1467  continue;
1468  std::cout << setw(4) << i << ")" << setw(8) << fixed << setprecision(4) << tz << " +/-" << setw(6)
1469  << sqrt(1. / tks.dz2[i]);
1470 
1471  if (tks.dt2[i] > 0) {
1472  std::cout << setw(8) << fixed << setprecision(4) << tks.tpca[i] << " +/-" << setw(6) << sqrt(1. / tks.dt2[i]);
1473  } else {
1474  std::cout << " ";
1475  }
1476 
1477  if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) {
1478  std::cout << " *";
1479  } else {
1480  std::cout << " ";
1481  }
1482  if (tks.tt[i]->track().hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) {
1483  std::cout << "+";
1484  } else {
1485  std::cout << "-";
1486  }
1487  std::cout << setw(1)
1488  << tks.tt[i]
1489  ->track()
1490  .hitPattern()
1491  .pixelBarrelLayersWithMeasurement(); // see DataFormats/TrackReco/interface/HitPattern.h
1492  std::cout << setw(1) << tks.tt[i]->track().hitPattern().pixelEndcapLayersWithMeasurement();
1493  std::cout << setw(1) << hex
1494  << tks.tt[i]->track().hitPattern().trackerLayersWithMeasurement() -
1495  tks.tt[i]->track().hitPattern().pixelLayersWithMeasurement()
1496  << dec;
1497  std::cout << "=" << setw(1) << hex
1498  << tks.tt[i]->track().hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << dec;
1499 
1500  Measurement1D IP = tks.tt[i]->stateAtBeamLine().transverseImpactParameter();
1501  std::cout << setw(8) << IP.value() << "+/-" << setw(6) << IP.error();
1502  std::cout << " " << setw(6) << setprecision(2) << tks.tt[i]->track().pt() * tks.tt[i]->track().charge();
1503  std::cout << " " << setw(5) << setprecision(2) << tks.tt[i]->track().phi() << " " << setw(5) << setprecision(2)
1504  << tks.tt[i]->track().eta();
1505 
1506  double sump = 0.;
1507  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1508  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1509  continue;
1510 
1511  if ((tks.tkwt[i] > 0) && (tks.sum_Z[i] > 0)) {
1512  double p =
1513  y.rho[ivertex] *
1514  local_exp(-beta *
1515  Eik(tks.zpca[i], y.zvtx[ivertex], tks.dz2[i], tks.tpca[i], y.tvtx[ivertex], tks.dt2[i])) /
1516  tks.sum_Z[i];
1517 
1518  if ((ivertex >= tks.kmin[i]) && (ivertex < tks.kmax[i])) {
1519  if (p > 0.0001) {
1520  std::cout << setw(8) << setprecision(3) << p;
1521  } else {
1522  std::cout << " _ "; // tiny but in the cluster list
1523  }
1524  E += p * Eik(tks.zpca[i], y.zvtx[ivertex], tks.dz2[i], tks.tpca[i], y.tvtx[ivertex], tks.dt2[i]);
1525  sump += p;
1526  } else {
1527  if (p > 0.1) {
1528  std::cout << "XXXXXXXX"; // we have an inconsistency here
1529  } else if (p > 0.0001) {
1530  std::cout << "X" << setw(6) << setprecision(3) << p << "X";
1531  } else {
1532  std::cout << " . "; // not in the cluster list
1533  }
1534  }
1535  } else {
1536  std::cout << " ";
1537  }
1538  }
1539  std::cout << " ( " << std::setw(3) << tks.kmin[i] << "," << std::setw(3) << tks.kmax[i] - 1 << " ) ";
1540  std::cout << endl;
1541  }
1542  std::cout << " ";
1543  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1544  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1545  std::cout << " " << setw(3) << ivertex << " ";
1546  }
1547  }
1548  std::cout << endl;
1549  std::cout << " z= ";
1550  std::cout << setprecision(4);
1551  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1552  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1553  std::cout << setw(8) << fixed << y.zvtx[ivertex];
1554  }
1555  }
1556  std::cout << endl;
1557  std::cout << endl
1558  << "T=" << 1 / beta << " E=" << E << " n=" << y.getSize() << " F= " << F << endl
1559  << "----------" << endl;
1560  }
1561 #endif
1562 }
double get_Tc(const vertex_t &y, int k) const
T sqrt(T t)
Definition: SSEVec.h:23
int nt
Definition: AMPTWrapper.h:42
const int verbosity
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:163
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4

◆ fill()

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

Definition at line 206 of file DAClusterizerInZT_vect.cc.

References funct::abs(), DAClusterizerInZT_vect::track_t::addItem(), HLT_2024v14_cff::atIP, gather_cfg::cout, TransientTrackBuilder::defaultInvalidTrackTimeReso, geometryDiff::epsilon, DAClusterizerInZT_vect::track_t::extractRaw(), DAClusterizerInZT_vect::track_t::getSize(), edm::isNotFinite(), SiStripPI::min, DAClusterizerInZT_vect::track_t::osumtkwt, funct::pow(), AlignmentPI::t_z, and DiMuonV_cfg::tracks.

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

◆ fillPSetDescription()

void DAClusterizerInZT_vect::fillPSetDescription ( edm::ParameterSetDescription desc)
static

Definition at line 1564 of file DAClusterizerInZT_vect.cc.

References submitPVResolutionJobs::desc, and DAClusterizerInZ_vect::fillPSetDescription().

Referenced by PrimaryVertexProducer::fillDescriptions().

1564  {
1566  desc.add<double>("tmerge", 0.1); // 4D only
1567  desc.add<double>("dtCutOff", 4.); // 4D only
1568  desc.add<double>("t0Max", 1.0); // 4D only
1569  desc.add<double>("vertexSizeTime", 0.008); // 4D only
1570 }
static void fillPSetDescription(edm::ParameterSetDescription &desc)

◆ 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 595 of file DAClusterizerInZT_vect.cc.

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

596  {
597  // find the cluster nearest to (z,t)
598  // distance measure is delta = (delta_z / dz )**2 + (delta_t/ d_t)**2
599  // assumes that clusters are ordered n z
600  // return value is false if no neighbour with distance < 1 is found
601 
602  unsigned int nv = y.getSize();
603  if (nv < 2) {
604  k_min = 0;
605  return false;
606  }
607 
608  // find nearest in z, binary search later
609  unsigned int k = 0;
610  for (unsigned int k0 = 1; k0 < nv; k0++) {
611  if (std::abs(y.zvtx[k0] - z) < std::abs(y.zvtx[k] - z)) {
612  k = k0;
613  }
614  }
615 
616  double delta_min = 1.;
617 
618  //search left
619  unsigned int k1 = k;
620  while ((k1 > 0) && ((y.zvtx[k] - y.zvtx[--k1]) < dz)) {
621  auto delta = std::pow((y.zvtx[k] - y.zvtx[k1]) / dz, 2) + std::pow((y.tvtx[k] - y.tvtx[k1]) / dt, 2);
622  if (delta < delta_min) {
623  k_min = k1;
624  delta_min = delta;
625  }
626  }
627 
628  //search right
629  k1 = k;
630  while (((++k1) < nv) && ((y.zvtx[k1] - y.zvtx[k]) < dz)) {
631  auto delta = std::pow((y.zvtx[k1] - y.zvtx[k]) / dz, 2) + std::pow((y.tvtx[k1] - y.tvtx[k]) / dt, 2);
632  if (delta < delta_min) {
633  k_min = k1;
634  delta_min = delta;
635  }
636  }
637 
638  return (delta_min < 1.);
639 }
float dt
Definition: AMPTWrapper.h:136
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ get_Tc()

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

Definition at line 899 of file DAClusterizerInZT_vect.cc.

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

899  {
900  double Tz = y.szz[k] / y.nuz[k]; // actually 0.5*Tc(z)
901  double Tt = 0.;
902  if (y.nut[k] > 0) {
903  Tt = y.stt[k] / y.nut[k];
904  double mx = y.szt[k] / y.nuz[k] * y.szt[k] / y.nut[k];
905  return Tz + Tt + sqrt(pow(Tz - Tt, 2) + 4 * mx);
906  }
907  return 2 * Tz;
908 }
T sqrt(T t)
Definition: SSEVec.h:23
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ merge()

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

Definition at line 692 of file DAClusterizerInZT_vect.cc.

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

692  {
693  // merge clusters that collapsed or never separated,
694  // return true if vertices were merged, false otherwise
695  const unsigned int nv = y.getSize();
696 
697  if (nv < 2)
698  return false;
699 
700  // merge the smallest distance clusters first
701  unsigned int k1_min = 0, k2_min = 0;
702  double delta_min = 0;
703 
704  for (unsigned int k1 = 0; (k1 + 1) < nv; k1++) {
705  unsigned int k2 = k1;
706  while ((++k2 < nv) && (std::fabs(y.zvtx[k2] - y.zvtx[k1]) < zmerge_)) {
707  auto delta = std::pow((y.zvtx[k2] - y.zvtx[k1]) / zmerge_, 2) + std::pow((y.tvtx[k2] - y.tvtx[k1]) / tmerge_, 2);
708  if ((delta < delta_min) || (k1_min == k2_min)) {
709  k1_min = k1;
710  k2_min = k2;
711  delta_min = delta;
712  }
713  }
714  }
715 
716  if ((k1_min == k2_min) || (delta_min > 1)) {
717  return false;
718  }
719 
720  double rho = y.rho[k1_min] + y.rho[k2_min];
721 
722 #ifdef DEBUG
723  assert((k1_min < nv) && (k2_min < nv));
724  if (DEBUGLEVEL > 1) {
725  std::cout << "merging (" << setw(8) << fixed << setprecision(4) << y.zvtx[k1_min] << ',' << y.tvtx[k1_min]
726  << ") and (" << y.zvtx[k2_min] << ',' << y.tvtx[k2_min] << ")"
727  << " idx=" << k1_min << "," << k2_min << std::endl;
728  }
729 #endif
730 
731  if (rho > 0) {
732  y.zvtx[k1_min] = (y.rho[k1_min] * y.zvtx[k1_min] + y.rho[k2_min] * y.zvtx[k2_min]) / rho;
733  y.tvtx[k1_min] = (y.rho[k1_min] * y.tvtx[k1_min] + y.rho[k2_min] * y.tvtx[k2_min]) / rho;
734 #ifdef USEVTXDT2
735  y.dt2[k1_min] = (y.rho[k1_min] * y.dt2[k1_min] + y.rho[k2_min] * y.dt2[k2_min]) / rho;
736 #endif
737  } else {
738  y.zvtx[k1_min] = 0.5 * (y.zvtx[k1_min] + y.zvtx[k2_min]);
739  y.tvtx[k1_min] = 0.5 * (y.tvtx[k1_min] + y.tvtx[k2_min]);
740  }
741  y.rho[k1_min] = rho;
742  y.removeItem(k2_min, tks);
743 
744  zorder(y);
745  set_vtx_range(beta, tks, y);
746  y.extractRaw();
747  return true;
748 }
assert(be >=bs)
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const
bool zorder(vertex_t &y) const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ purge()

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

Definition at line 750 of file DAClusterizerInZT_vect.cc.

References cms::cuda::assert(), HLT_2024v14_cff::beta, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), gather_cfg::cout, DAClusterizerInZT_vect::track_t::dt2, DAClusterizerInZT_vect::track_t::dz2, HLT_2024v14_cff::eps, cms::cuda::for(), DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, dqmdumpme::k, reco::ParticleMasses::k0, DAClusterizerInZT_vect::track_t::kmax, DAClusterizerInZT_vect::track_t::kmin, nt, AlCaHLTBitMon_ParallelJobs::p, DAClusterizerInZT_vect::track_t::sum_Z, DAClusterizerInZT_vect::track_t::tkwt, DAClusterizerInZT_vect::track_t::tpca, and DAClusterizerInZT_vect::track_t::zpca.

750  {
751  constexpr double eps = 1.e-100;
752  // eliminate clusters with only one significant/unique track
753  const unsigned int nv = y.getSize();
754  const unsigned int nt = tks.getSize();
755 
756  if (nv < 2)
757  return false;
758 
759  std::vector<double> sump_v(nv), arg_cache_v(nv), exp_cache_v(nv), pcut_cache_v(nv);
760  std::vector<int> nUnique_v(nv);
761  double* __restrict__ parg_cache;
762  double* __restrict__ pexp_cache;
763  double* __restrict__ ppcut_cache;
764  double* __restrict__ psump;
765  int* __restrict__ pnUnique;
766  int constexpr nunique_min_ = 2;
767 
768  zorder(y);
769  set_vtx_range(beta, tks, y);
770 
771  parg_cache = arg_cache_v.data();
772  pexp_cache = exp_cache_v.data();
773  ppcut_cache = pcut_cache_v.data();
774  psump = sump_v.data();
775  pnUnique = nUnique_v.data();
776 
777  const auto rhoconst = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
778  for (unsigned int k = 0; k < nv; ++k) {
779  const double pmax = y.rho[k] / (y.rho[k] + rhoconst);
780  ppcut_cache[k] = uniquetrkweight_ * pmax;
781  }
782 
783  for (unsigned int i = 0; i < nt; i++) {
784  const auto invZ = ((tks.sum_Z[i] > eps) && (tks.tkwt[i] > uniquetrkminp_)) ? 1. / tks.sum_Z[i] : 0.;
785  const auto track_z = tks.zpca[i];
786  const auto track_t = tks.tpca[i];
787  const auto botrack_dz2 = -beta * tks.dz2[i];
788  const auto botrack_dt2 = -beta * tks.dt2[i];
789  const auto kmin = tks.kmin[i];
790  const auto kmax = tks.kmax[i];
791 
792  for (unsigned int k = kmin; k < kmax; k++) {
793  const auto mult_resz = track_z - y.zvtx[k];
794  const auto mult_rest = track_t - y.tvtx[k];
795  parg_cache[k] = botrack_dz2 * (mult_resz * mult_resz) + botrack_dt2 * (mult_rest * mult_rest);
796  }
797 
798  local_exp_list_range(parg_cache, pexp_cache, kmin, kmax);
799 
800  for (unsigned int k = kmin; k < kmax; k++) {
801  const double p = y.rho[k] * pexp_cache[k] * invZ;
802  psump[k] += p;
803  pnUnique[k] += (p > ppcut_cache[k]) ? 1 : 0;
804  }
805  }
806 
807  double sumpmin = nt;
808  unsigned int k0 = nv;
809  for (unsigned k = 0; k < nv; k++) {
810  if ((pnUnique[k] < nunique_min_) && (psump[k] < sumpmin)) {
811  sumpmin = psump[k];
812  k0 = k;
813  }
814  }
815 
816  if (k0 != nv) {
817 #ifdef DEBUG
818  assert(k0 < y.getSize());
819  if (DEBUGLEVEL > 1) {
820  std::cout << "eliminating prototype at " << std::setw(10) << std::setprecision(4) << y.zvtx[k0] << ","
821  << y.tvtx[k0] << " with sump=" << sumpmin << " rho*nt =" << y.rho[k0] * nt << endl;
822  }
823 #endif
824  y.removeItem(k0, tks);
825  set_vtx_range(beta, tks, y);
826  return true;
827  } else {
828  return false;
829  }
830 }
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
assert(be >=bs)
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const
bool zorder(vertex_t &y) const
int nt
Definition: AMPTWrapper.h:42

◆ set_vtx_range()

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

Definition at line 287 of file DAClusterizerInZT_vect.cc.

References HLT_2024v14_cff::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, SiStripPI::min, nt, mathSSE::sqrt(), mitigatedMETSequence_cff::U, SiStripMonitorCluster_cfi::zmax, SiStripMonitorCluster_cfi::zmin, DAClusterizerInZT_vect::track_t::zpca, OfflinePixel3DPrimaryVertices_cfi::zrange, and DAClusterizerInZT_vect::vertex_t::zvtx.

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

◆ split()

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

Definition at line 910 of file DAClusterizerInZT_vect.cc.

References funct::abs(), cms::cuda::assert(), HLT_2024v14_cff::beta, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), funct::cos(), gather_cfg::cout, hlt_jetmet_dqm_QT_fromfile_cfg::critical, DAClusterizerInZT_vect::track_t::dt2, DAClusterizerInZT_vect::track_t::dz2, MillePedeFileConverter_cfg::e, alignBH_cfg::fixed, DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, dqmdumpme::k, nt, AlCaHLTBitMon_ParallelJobs::p, LaserDQM_cfg::p1, SiStripOfflineCRack_cfg::p2, funct::pow(), edm::second(), funct::sin(), submitPVValidationJobs::split(), mathSSE::sqrt(), DAClusterizerInZT_vect::track_t::sum_Z, submitPVValidationJobs::t, RandomServiceHelper::t1, RandomServiceHelper::t2, DiMuonV_cfg::threshold, DAClusterizerInZT_vect::track_t::tkwt, DAClusterizerInZT_vect::track_t::tpca, testProducerWithPsetDescEmpty_cfi::z2, and DAClusterizerInZT_vect::track_t::zpca.

910  {
911  // split only critical vertices (Tc >~ T=1/beta <==> beta*Tc>~1)
912  // an update must have been made just before doing this (same beta, no merging)
913  // returns true if at least one cluster was split
914 
915  constexpr double epsilonz = 1e-3; // minimum split size z
916  constexpr double epsilont = 1e-2; // minimum split size t
917  unsigned int nv = y.getSize();
918  const double twoBeta = 2.0 * beta;
919 
920  // avoid left-right biases by splitting highest Tc first
921 
922  std::vector<std::pair<double, unsigned int> > critical;
923  for (unsigned int k = 0; k < nv; k++) {
924  double Tc = get_Tc(y, k);
925  if (beta * Tc > threshold) {
926  critical.push_back(make_pair(Tc, k));
927  }
928  }
929  if (critical.empty())
930  return false;
931 
932  std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int> >());
933 
934  bool split = false;
935  const unsigned int nt = tks.getSize();
936 
937  for (unsigned int ic = 0; ic < critical.size(); ic++) {
938  unsigned int k = critical[ic].second;
939 
940  // split direction in the (z,t)-plane
941 
942  double Mzz = y.nuz[k] - twoBeta * y.szz[k];
943  double Mtt = y.nut[k] - twoBeta * y.stt[k];
944  double Mzt = -twoBeta * y.szt[k];
945  const double twoMzt = 2.0 * Mzt;
946  double D = sqrt(pow(Mtt - Mzz, 2) + twoMzt * twoMzt);
947  double q1 = atan2(-Mtt + Mzz + D, -twoMzt);
948  double l1 = 0.5 * (-Mzz - Mtt + D);
949  double l2 = 0.5 * (-Mzz - Mtt - D);
950  if ((std::abs(l1) < 1e-4) && (std::abs(l2) < 1e-4)) {
951  edm::LogWarning("DAClusterizerInZT_vect") << "warning, bad eigenvalues! idx=" << k << " z= " << y.zvtx[k]
952  << " Mzz=" << Mzz << " Mtt=" << Mtt << " Mzt=" << Mzt << endl;
953  }
954 
955  double qsplit = q1;
956  double cq = cos(qsplit);
957  double sq = sin(qsplit);
958  if (cq < 0) {
959  cq = -cq;
960  sq = -sq;
961  } // prefer cq>0 to keep z-ordering
962 
963  // estimate subcluster positions and weight
964  double p1 = 0, z1 = 0, t1 = 0, wz1 = 0, wt1 = 0;
965  double p2 = 0, z2 = 0, t2 = 0, wz2 = 0, wt2 = 0;
966  for (unsigned int i = 0; i < nt; ++i) {
967  if (tks.sum_Z[i] > 1.e-100) {
968  double lr = (tks.zpca[i] - y.zvtx[k]) * cq + (tks.tpca[i] - y.tvtx[k]) * sq;
969  // winner-takes-all, usually overestimates splitting
970  double tl = lr < 0 ? 1. : 0.;
971  double tr = 1. - tl;
972 
973  // soften it, especially at low T
974  double arg = lr * std::sqrt(beta * (cq * cq * tks.dz2[i] + sq * sq * tks.dt2[i]));
975  if (std::abs(arg) < 20) {
976  double t = local_exp(-arg);
977  tl = t / (t + 1.);
978  tr = 1 / (t + 1.);
979  }
980 
981  double p = y.rho[k] * tks.tkwt[i] *
982  local_exp(-beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i], tks.tpca[i], y.tvtx[k], tks.dt2[i])) /
983  tks.sum_Z[i];
984  double wz = p * tks.dz2[i];
985  double wt = p * tks.dt2[i];
986  p1 += p * tl;
987  z1 += wz * tl * tks.zpca[i];
988  t1 += wt * tl * tks.tpca[i];
989  wz1 += wz * tl;
990  wt1 += wt * tl;
991  p2 += p * tr;
992  z2 += wz * tr * tks.zpca[i];
993  t2 += wt * tr * tks.tpca[i];
994  wz2 += wz * tr;
995  wt2 += wt * tr;
996  }
997  }
998 
999  if (wz1 > 0) {
1000  z1 /= wz1;
1001  } else {
1002  z1 = y.zvtx[k] - epsilonz * cq;
1003  edm::LogWarning("DAClusterizerInZT_vect") << "warning, wz1 = " << scientific << wz1 << endl;
1004  }
1005  if (wt1 > 0) {
1006  t1 /= wt1;
1007  } else {
1008  t1 = y.tvtx[k] - epsilont * sq;
1009  edm::LogWarning("DAClusterizerInZT_vect") << "warning, wt1 = " << scientific << wt1 << endl;
1010  }
1011  if (wz2 > 0) {
1012  z2 /= wz2;
1013  } else {
1014  z2 = y.zvtx[k] + epsilonz * cq;
1015  edm::LogWarning("DAClusterizerInZT_vect") << "warning, wz2 = " << scientific << wz2 << endl;
1016  }
1017  if (wt2 > 0) {
1018  t2 /= wt2;
1019  } else {
1020  t2 = y.tvtx[k] + epsilont * sq;
1021  edm::LogWarning("DAClusterizerInZT_vect") << "warning, wt2 = " << scientific << wt2 << endl;
1022  }
1023 
1024  unsigned int k_min1 = k, k_min2 = k;
1025  constexpr double spliteps = 1e-8;
1026  while (((find_nearest(z1, t1, y, k_min1, epsilonz, epsilont) && (k_min1 != k)) ||
1027  (find_nearest(z2, t2, y, k_min2, epsilonz, epsilont) && (k_min2 != k))) &&
1028  (std::abs(z2 - z1) > spliteps || std::abs(t2 - t1) > spliteps)) {
1029  z1 = 0.5 * (z1 + y.zvtx[k]);
1030  t1 = 0.5 * (t1 + y.tvtx[k]);
1031  z2 = 0.5 * (z2 + y.zvtx[k]);
1032  t2 = 0.5 * (t2 + y.tvtx[k]);
1033  }
1034 
1035 #ifdef DEBUG
1036  assert(k < nv);
1037  if (DEBUGLEVEL > 1) {
1038  if (std::fabs(y.zvtx[k] - zdumpcenter_) < zdumpwidth_) {
1039  std::cout << " T= " << std::setw(10) << std::setprecision(1) << 1. / beta << " Tc= " << critical[ic].first
1040  << " direction =" << std::setprecision(4) << qsplit << " splitting (" << std::setw(8) << std::fixed
1041  << std::setprecision(4) << y.zvtx[k] << "," << y.tvtx[k] << ")"
1042  << " --> (" << z1 << ',' << t1 << "),(" << z2 << ',' << t2 << ") [" << p1 << "," << p2 << "]";
1043  if (std::fabs(z2 - z1) > epsilonz || std::fabs(t2 - t1) > epsilont) {
1044  std::cout << std::endl;
1045  } else {
1046  std::cout << " rejected " << std::endl;
1047  }
1048  }
1049  }
1050 #endif
1051 
1052  if (z1 > z2) {
1053  edm::LogInfo("DAClusterizerInZT") << "warning swapping z in split qsplit=" << qsplit << " cq=" << cq
1054  << " sq=" << sq << endl;
1055  auto ztemp = z1;
1056  auto ttemp = t1;
1057  auto ptemp = p1;
1058  z1 = z2;
1059  t1 = t2;
1060  p1 = p2;
1061  z2 = ztemp;
1062  t2 = ttemp;
1063  p2 = ptemp;
1064  }
1065 
1066  // split if the new subclusters are significantly separated
1067  if (std::fabs(z2 - z1) > epsilonz || std::fabs(t2 - t1) > epsilont) {
1068  split = true;
1069  double rho1 = p1 * y.rho[k] / (p1 + p2);
1070  double rho2 = p2 * y.rho[k] / (p1 + p2);
1071 
1072  // replace the original by (z2,t2)
1073  y.removeItem(k, tks);
1074  unsigned int k2 = y.insertOrdered(z2, t2, rho2, tks);
1075 
1076 #ifdef DEBUG
1077  if (k2 < k) {
1078  std::cout << "unexpected z-ordering in split" << std::endl;
1079  }
1080 #endif
1081 
1082  // adjust pointers if necessary
1083  if (!(k == k2)) {
1084  for (unsigned int jc = ic; jc < critical.size(); jc++) {
1085  if (critical[jc].second > k) {
1086  critical[jc].second--;
1087  }
1088  if (critical[jc].second >= k2) {
1089  critical[jc].second++;
1090  }
1091  }
1092  }
1093 
1094  // insert (z1,t1) where it belongs
1095  unsigned int k1 = y.insertOrdered(z1, t1, rho1, tks);
1096  nv++;
1097 
1098  // adjust remaining pointers
1099  for (unsigned int jc = ic; jc < critical.size(); jc++) {
1100  if (critical[jc].second >= k1) {
1101  critical[jc].second++;
1102  }
1103  }
1104  } else {
1105 #ifdef DEBUG
1106  std::cout << "warning ! split rejected, too small." << endl;
1107 #endif
1108  }
1109  }
1110 
1111  return split;
1112 }
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
double get_Tc(const vertex_t &y, int k) const
assert(be >=bs)
A arg
Definition: Factorize.h:31
U second(std::pair< T, U > const &p)
bool split(const double beta, track_t &t, vertex_t &y, double threshold=1.) const
T sqrt(T t)
Definition: SSEVec.h:23
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int nt
Definition: AMPTWrapper.h:42
Log< level::Info, false > LogInfo
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
Log< level::Warning, false > LogWarning
bool find_nearest(double z, double t, vertex_t &y, unsigned int &k_min, double dz, double dt) const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ 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 641 of file DAClusterizerInZT_vect.cc.

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

642  {
643  unsigned int niter = 0;
644  double delta = 0;
645  double delta_max = delta_lowT_;
646 
647  if (convergence_mode_ == 0) {
648  delta_max = delta_max0;
649  } else if (convergence_mode_ == 1) {
650  delta_max = delta_lowT_ / sqrt(std::max(beta, 1.0));
651  }
652 
653  zorder(v);
654  set_vtx_range(beta, tks, v);
655 
656  double delta_sum_range = 0; // accumulate max(|delta-z|) as a lower bound
657  std::vector<double> z0 = v.zvtx_vec;
658 
659  while (niter++ < maxIterations_) {
660  delta = update(beta, tks, v, rho0);
661  delta_sum_range += delta;
662 
663  if (delta_sum_range > zrange_min_) {
664  for (unsigned int k = 0; k < v.getSize(); k++) {
665  if (std::abs(v.zvtx_vec[k] - z0[k]) > zrange_min_) {
666  zorder(v);
667  set_vtx_range(beta, tks, v);
668  delta_sum_range = 0;
669  z0 = v.zvtx_vec;
670  break;
671  }
672  }
673  }
674 
675  if (delta < delta_max) {
676  break;
677  }
678  }
679 
680 #ifdef DEBUG
681  if (DEBUGLEVEL > 0) {
682  std::cout << "DAClusterizerInZT_vect.thermalize niter = " << niter << " at T = " << 1 / beta
683  << " nv = " << v.getSize() << std::endl;
684  if (DEBUGLEVEL > 2)
685  dump(beta, v, tks, 0);
686  }
687 #endif
688 
689  return niter;
690 }
double update(double beta, track_t &gtracks, vertex_t &gvertices, const double rho0=0, const bool updateTc=false) const
T sqrt(T t)
Definition: SSEVec.h:23
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool zorder(vertex_t &y) const
void dump(const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0) const

◆ update()

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

Definition at line 353 of file DAClusterizerInZT_vect.cc.

References funct::abs(), cms::cuda::assert(), HLT_2024v14_cff::beta, gather_cfg::cout, dumpMFGeometry_cfg::delta, math::cholesky::dst, run3scouting_cff::dsz, DAClusterizerInZT_vect::vertex_t::exp, DAClusterizerInZT_vect::vertex_t::exp_arg, 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, DAClusterizerInZT_vect::vertex_t::nuz, DAClusterizerInZT_vect::track_t::osumtkwt, DAClusterizerInZT_vect::vertex_t::se, DAClusterizerInZT_vect::vertex_t::stt, DAClusterizerInZT_vect::track_t::sum_Z, DAClusterizerInZT_vect::track_t::sum_Z_vec, DAClusterizerInZT_vect::vertex_t::swt, DAClusterizerInZT_vect::vertex_t::swz, DAClusterizerInZT_vect::vertex_t::szt, DAClusterizerInZT_vect::vertex_t::szz, DiMuonV_cfg::tracks, mitigatedMETSequence_cff::U, AlignmentTracksFromVertexSelector_cfi::vertices, and w().

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

354  {
355  //update weights and vertex positions
356  // returns the maximum of changes of vertex positions
357  // sums needed for Tc are only updated if updateTC == true
358 
359  const unsigned int nt = gtracks.getSize();
360  const unsigned int nv = gvertices.getSize();
361  auto osumtkwt = gtracks.osumtkwt;
362 
363  double Z_init = 0;
364  if (rho0 > 0) {
365  Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); // cut-off
366  }
367 
368  // define kernels
369 
370  auto kernel_calc_exp_arg_range = [beta](const unsigned int itrack,
371  track_t const& tracks,
372  vertex_t const& vertices,
373  const unsigned int kmin,
374  const unsigned int kmax) {
375  const auto track_z = tracks.zpca[itrack];
376  const auto track_t = tracks.tpca[itrack];
377  const auto botrack_dz2 = -beta * tracks.dz2[itrack];
378 #ifndef USEVTXDT2
379  const auto botrack_dt2 = -beta * tracks.dt2[itrack];
380 #endif
381  // auto-vectorized
382  for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
383  const auto mult_resz = track_z - vertices.zvtx[ivertex];
384  const auto mult_rest = track_t - vertices.tvtx[ivertex];
385 #ifdef USEVTXDT2
386  const auto botrack_dt2 =
387  -beta * tracks.dt2[itrack] * vertices.dt2[ivertex] / (tracks.dt2[itrack] + vertices.dt2[ivertex] + 1.e-10);
388 #endif
389  vertices.exp_arg[ivertex] = botrack_dz2 * (mult_resz * mult_resz) + botrack_dt2 * (mult_rest * mult_rest);
390  }
391  };
392 
393  auto kernel_add_Z_range = [Z_init](
394  vertex_t const& vertices, const unsigned int kmin, const unsigned int kmax) -> double {
395  double ZTemp = Z_init;
396  for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
397  ZTemp += vertices.rho[ivertex] * vertices.exp[ivertex];
398  }
399  return ZTemp;
400  };
401 
402  auto kernel_calc_normalization_range = [updateTc](const unsigned int track_num,
403  track_t& tracks,
404  vertex_t& vertices,
405  const unsigned int kmin,
406  const unsigned int kmax) {
407  auto tmp_trk_tkwt = tracks.tkwt[track_num];
408  auto o_trk_sum_Z = 1. / tracks.sum_Z[track_num];
409  auto o_trk_err_z = tracks.dz2[track_num];
410  auto o_trk_err_t = tracks.dt2[track_num];
411  auto tmp_trk_z = tracks.zpca[track_num];
412  auto tmp_trk_t = tracks.tpca[track_num];
413  // auto-vectorized
414  if (updateTc) {
415 #pragma GCC ivdep
416  for (unsigned int k = kmin; k < kmax; ++k) {
417  // parens are important for numerical stability
418  vertices.se[k] += tmp_trk_tkwt * (vertices.exp[k] * o_trk_sum_Z);
419  const auto w = tmp_trk_tkwt * (vertices.rho[k] * vertices.exp[k] * o_trk_sum_Z); // p_{ik}
420  const auto wz = w * o_trk_err_z;
421  const auto wt = w * o_trk_err_t;
422 #ifdef USEVTXDT2
423  vertices.sumw[k] += w; // for vtxdt2
424 #endif
425  vertices.nuz[k] += wz;
426  vertices.nut[k] += wt;
427  vertices.swz[k] += wz * tmp_trk_z;
428  vertices.swt[k] += wt * tmp_trk_t;
429  // the following lines are only needed for Tc updates
430  const auto dsz = (tmp_trk_z - vertices.zvtx[k]) * o_trk_err_z;
431  const auto dst = (tmp_trk_t - vertices.tvtx[k]) * o_trk_err_t;
432  vertices.szz[k] += w * dsz * dsz;
433  vertices.stt[k] += w * dst * dst;
434  vertices.szt[k] += w * dsz * dst;
435  }
436  } else {
437 #pragma GCC ivdep
438  for (unsigned int k = kmin; k < kmax; ++k) {
439  vertices.se[k] += tmp_trk_tkwt * (vertices.exp[k] * o_trk_sum_Z);
440  const auto w = tmp_trk_tkwt * (vertices.rho[k] * vertices.exp[k] * o_trk_sum_Z); // p_{ik}
441  const auto wz = w * o_trk_err_z;
442  const auto wt = w * o_trk_err_t;
443 #ifdef USEVTXDT2
444  vertices.sumw[k] += w; // for vtxdt2
445 #endif
446  vertices.nuz[k] += wz;
447  vertices.nut[k] += wt;
448  vertices.swz[k] += wz * tmp_trk_z;
449  vertices.swt[k] += wt * tmp_trk_t;
450  }
451  }
452  };
453 
454  for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
455  gvertices.se[ivertex] = 0.0;
456  gvertices.nuz[ivertex] = 0.0;
457  gvertices.nut[ivertex] = 0.0;
458  gvertices.swz[ivertex] = 0.0;
459  gvertices.swt[ivertex] = 0.0;
460  gvertices.szz[ivertex] = 0.0;
461  gvertices.stt[ivertex] = 0.0;
462  gvertices.szt[ivertex] = 0.0;
463 #ifdef USEVTXDT2
464  gvertices.sumw[ivertex] = 0.0;
465 #endif
466  }
467 
468  // loop over tracks
469  for (auto itrack = 0U; itrack < nt; ++itrack) {
470  const unsigned int kmin = gtracks.kmin[itrack];
471  const unsigned int kmax = gtracks.kmax[itrack];
472 
473 #ifdef DEBUG
474  assert((kmin < kmax) && (kmax <= nv));
475  assert(itrack < gtracks.sum_Z_vec.size());
476 #endif
477 
478  kernel_calc_exp_arg_range(itrack, gtracks, gvertices, kmin, kmax);
479  local_exp_list_range(gvertices.exp_arg, gvertices.exp, kmin, kmax);
480  gtracks.sum_Z[itrack] = kernel_add_Z_range(gvertices, kmin, kmax);
481 
482  if (edm::isNotFinite(gtracks.sum_Z[itrack]))
483  gtracks.sum_Z[itrack] = 0.0;
484 
485  if (gtracks.sum_Z[itrack] > 1.e-100) {
486  kernel_calc_normalization_range(itrack, gtracks, gvertices, kmin, kmax);
487  }
488  }
489 
490  // now update z, t, and rho
491  auto kernel_calc_zt = [osumtkwt, nv](vertex_t& vertices) -> double {
492  double delta = 0;
493 
494  // does not vectorize
495  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
496  if (vertices.nuz[ivertex] > 0.) {
497  auto znew = vertices.swz[ivertex] / vertices.nuz[ivertex];
498  delta = max(std::abs(vertices.zvtx[ivertex] - znew), delta);
499  vertices.zvtx[ivertex] = znew;
500  }
501 
502  if (vertices.nut[ivertex] > 0.) {
503  auto tnew = vertices.swt[ivertex] / vertices.nut[ivertex];
504  //delta = max(std::abs(vertices.t[ ivertex ] - tnew), delta); // FIXME
505  vertices.tvtx[ivertex] = tnew;
506 #ifdef USEVTXDT2
507  vertices.dt2[ivertex] = vertices.nut[ivertex] / vertices.sumw[ivertex];
508 #endif
509  } else {
510  // FIXME
511  // apparently this cluster has not timing info attached
512  // this should be taken into account somehow, otherwise we might fail to attach
513  // new tracks that do have timing information
514  vertices.tvtx[ivertex] = 0;
515 #ifdef USEVTXDT2
516  vertices.dt2[ivertex] = 0;
517 #endif
518  }
519 
520 #ifdef DEBUG
521  if ((vertices.nut[ivertex] <= 0.) && (vertices.nuz[ivertex] <= 0.)) {
522  edm::LogInfo("sumw") << "invalid sum of weights in fit: " << endl;
523  std::cout << " a cluster melted away ? "
524  << " zk=" << vertices.zvtx[ivertex] << " rho=" << vertices.rho[ivertex]
525  << " sumw(z,t) =" << vertices.nuz[ivertex] << "," << vertices.nut[ivertex] << endl;
526  // FIXME: discard this cluster
527  }
528 #endif
529  }
530 
531  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
532  vertices.rho[ivertex] = vertices.rho[ivertex] * vertices.se[ivertex] * osumtkwt;
533  }
534 
535  return delta;
536  };
537 
538  double delta = kernel_calc_zt(gvertices);
539 
540  if (zorder(gvertices)) {
541  set_vtx_range(beta, gtracks, gvertices);
542  };
543 
544  // return how much the prototypes moved
545  return delta;
546 }
T w() const
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
assert(be >=bs)
std::vector< TransientVertex > vertices(const std::vector< reco::TransientTrack > &tracks) const override
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool zorder(vertex_t &y) const
int nt
Definition: AMPTWrapper.h:42
Log< level::Info, false > LogInfo
constexpr void M2 & dst

◆ 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 116 of file DAClusterizerInZT_vect.cc.

References cms::cuda::assert(), gather_cfg::cout, DAClusterizerInZT_vect::track_t::dt2, DAClusterizerInZT_vect::track_t::dt2_vec, DAClusterizerInZT_vect::track_t::dz2, DAClusterizerInZT_vect::track_t::dz2_vec, 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::sum_Z, DAClusterizerInZT_vect::track_t::sum_Z_vec, DAClusterizerInZT_vect::track_t::tkwt, DAClusterizerInZT_vect::track_t::tkwt_vec, DAClusterizerInZT_vect::track_t::tpca, DAClusterizerInZT_vect::track_t::tpca_vec, DAClusterizerInZT_vect::track_t::tt, findQualityFiles::v, DAClusterizerInZT_vect::track_t::zpca, and DAClusterizerInZT_vect::track_t::zpca_vec.

116  {
117  if (nv == 999999) {
118  nv = v.getSize();
119  } else {
120  assert(nv == v.getSize());
121  }
122 
123  if (nt == 999999) {
124  nt = tks.getSize();
125  } else {
126  assert(nt == tks.getSize());
127  }
128 
129  // clusters
130  assert(v.zvtx_vec.size() == nv);
131  assert(v.tvtx_vec.size() == nv);
132 #ifdef USEVTXDT2
133  assert(v.dt2_vec.size() == nv);
134  assert(v.sumw_vec.size() == nv);
135 #endif
136  assert(v.rho_vec.size() == nv);
137  assert(v.swz_vec.size() == nv);
138  assert(v.swt_vec.size() == nv);
139  assert(v.exp_arg_vec.size() == nv);
140  assert(v.exp_vec.size() == nv);
141  assert(v.se_vec.size() == nv);
142  assert(v.nuz_vec.size() == nv);
143  assert(v.nut_vec.size() == nv);
144  assert(v.szz_vec.size() == nv);
145  assert(v.stt_vec.size() == nv);
146  assert(v.szt_vec.size() == nv);
147 
148  assert(v.zvtx == &v.zvtx_vec.front());
149  assert(v.tvtx == &v.tvtx_vec.front());
150  assert(v.rho == &v.rho_vec.front());
151  assert(v.exp_arg == &v.exp_arg_vec.front());
152  assert(v.swz == &v.swz_vec.front());
153  assert(v.swt == &v.swt_vec.front());
154  assert(v.se == &v.se_vec.front());
155  assert(v.nuz == &v.nuz_vec.front());
156  assert(v.nut == &v.nut_vec.front());
157  assert(v.szz == &v.szz_vec.front());
158  assert(v.stt == &v.stt_vec.front());
159  assert(v.szt == &v.szt_vec.front());
160 #ifdef USEVTXDT2
161  assert(v.sumw == &v.sumw_vec.front());
162  assert(v.dt2 == &v.dt2_vec.front());
163 #endif
164 
165  for (unsigned int k = 0; k < nv - 1; k++) {
166  if (v.zvtx[k] <= v.zvtx[k + 1])
167  continue;
168  cout << " ZT, cluster z-ordering assertion failure z[" << k << "] =" << v.zvtx[k] << " z[" << k + 1
169  << "] =" << v.zvtx[k + 1] << endl;
170  }
171  //for(unsigned int k=0; k< nv-1; k++){
172  //assert( v.z[k] <= v.z[k+1]);
173  //}
174 
175  // tracks
176  assert(nt == tks.zpca_vec.size());
177  assert(nt == tks.tpca_vec.size());
178  assert(nt == tks.dz2_vec.size());
179  assert(nt == tks.dt2_vec.size());
180  assert(nt == tks.tt.size());
181  assert(nt == tks.tkwt_vec.size());
182  assert(nt == tks.sum_Z_vec.size());
183  assert(nt == tks.kmin.size());
184  assert(nt == tks.kmax.size());
185 
186  assert(tks.zpca == &tks.zpca_vec.front());
187  assert(tks.tpca == &tks.tpca_vec.front());
188  assert(tks.dz2 == &tks.dz2_vec.front());
189  assert(tks.dt2 == &tks.dt2_vec.front());
190  assert(tks.tkwt == &tks.tkwt_vec.front());
191  assert(tks.sum_Z == &tks.sum_Z_vec.front());
192 
193  for (unsigned int i = 0; i < nt; i++) {
194  if ((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv))
195  continue;
196  cout << "track vertex range assertion failure" << i << "/" << nt << " kmin,kmax=" << tks.kmin[i] << ", "
197  << tks.kmax[i] << " nv=" << nv << endl;
198  }
199 
200  for (unsigned int i = 0; i < nt; i++) {
201  assert((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv));
202  }
203 }
assert(be >=bs)
int nt
Definition: AMPTWrapper.h:42

◆ vertices()

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

Implements TrackClusterizerInZ.

Definition at line 1114 of file DAClusterizerInZT_vect.cc.

References a, HLT_2024v14_cff::beta, bsc_activity_cfg::clusters, gather_cfg::cout, DAClusterizerInZT_vect::track_t::dt2, GCP_Ntuples_cfg::dump, DAClusterizerInZT_vect::track_t::dz2, MillePedeFileConverter_cfg::e, relativeConstraints::empty, DAClusterizerInZT_vect::track_t::extractRaw(), ntuplemaker::fill, DAClusterizerInZT_vect::track_t::getSize(), mps_fire::i, edm::isNotFinite(), dqmdumpme::k, DAClusterizerInZT_vect::track_t::kmax, DAClusterizerInZT_vect::track_t::kmin, merge(), SiStripPI::min, nt, AlCaHLTBitMon_ParallelJobs::p, submitPVValidationJobs::split(), mathSSE::sqrt(), DAClusterizerInZT_vect::track_t::sum_Z, DiMuonV_cfg::threshold, DAClusterizerInZT_vect::track_t::tpca, DiMuonV_cfg::tracks, DAClusterizerInZT_vect::track_t::tt, update, findQualityFiles::v, and DAClusterizerInZT_vect::track_t::zpca.

1114  {
1115  track_t&& tks = fill(tracks);
1116  vector<TransientVertex> clusters;
1117  if (tks.getSize() == 0)
1118  return clusters;
1119  tks.extractRaw();
1120 
1121  unsigned int nt = tks.getSize();
1122  double rho0 = 0.0; // start with no outlier rejection
1123 
1124  vertex_t y; // the vertex prototypes
1125 
1126  // initialize:single vertex at infinite temperature
1127  y.addItem(0, 0, 1.0);
1128  clear_vtx_range(tks, y);
1129 
1130  // estimate first critical temperature
1131  double beta = beta0(betamax_, tks, y);
1132 #ifdef DEBUG
1133  if (DEBUGLEVEL > 0)
1134  std::cout << "Beta0 is " << beta << std::endl;
1135 #endif
1136 
1137  thermalize(beta, tks, y, delta_highT_, 0.);
1138 
1139  // annealing loop, stop when T<minT (i.e. beta>1/minT)
1140 
1141  double betafreeze = betamax_ * sqrt(coolingFactor_);
1142 
1143  while (beta < betafreeze) {
1144  update(beta, tks, y, rho0, true);
1145  while (merge(y, tks, beta)) {
1146  if (zorder(y))
1147  set_vtx_range(beta, tks, y);
1148  update(beta, tks, y, rho0, true);
1149  }
1150  split(beta, tks, y);
1151 
1152  beta = beta / coolingFactor_;
1153  thermalize(beta, tks, y, delta_highT_, 0.);
1154  }
1155 
1156 #ifdef DEBUG
1157  verify(y, tks);
1158 
1159  if (DEBUGLEVEL > 0) {
1160  std::cout << "DAClusterizerInZT_vect::vertices :"
1161  << "merging at T=" << 1 / beta << std::endl;
1162  }
1163 #endif
1164 
1165  zorder(y);
1166  set_vtx_range(beta, tks, y);
1167  // ZT_vect version of merge() does not use Tc, but split() does
1168  update(beta, tks, y, rho0, true);
1169 
1170  while (merge(y, tks, beta)) {
1171  zorder(y);
1172  set_vtx_range(beta, tks, y);
1173  update(beta, tks, y, rho0, true);
1174  }
1175 
1176 #ifdef DEBUG
1177  verify(y, tks);
1178  if (DEBUGLEVEL > 0) {
1179  std::cout << "DAClusterizerInZT_vect::vertices :"
1180  << "splitting/merging at T=" << 1 / beta << std::endl;
1181  }
1182 #endif
1183 
1184  unsigned int ntry = 0;
1185  double threshold = 1.0;
1186  while (split(beta, tks, y, threshold) && (ntry++ < 10)) {
1187  thermalize(beta, tks, y, delta_highT_, 0.);
1188 
1189  while (merge(y, tks, beta)) {
1190  if (zorder(y))
1191  set_vtx_range(beta, tks, y);
1192  update(beta, tks, y, rho0, true); // updateTc in case we go back to split
1193  }
1194 
1195 #ifdef DEBUG
1196  verify(y, tks);
1197 #endif
1198 
1199  // relax splitting a bit to reduce multiple split-merge cycles of the same cluster
1200  threshold *= 1.1;
1201  }
1202 
1203 #ifdef DEBUG
1204  verify(y, tks);
1205  if (DEBUGLEVEL > 0) {
1206  std::cout << "DAClusterizerInZT_vect::vertices :"
1207  << "turning on outlier rejection at T=" << 1 / beta << std::endl;
1208  }
1209 #endif
1210 
1211  // switch on outlier rejection at T=minT
1212  if (dzCutOff_ > 0) {
1213  rho0 = 1. / nt;
1214  for (unsigned int a = 0; a < 5; a++) {
1215  update(beta, tks, y, a * rho0 / 5); // adiabatic turn-on
1216  }
1217  }
1218 
1219  thermalize(beta, tks, y, delta_lowT_, rho0);
1220 
1221 #ifdef DEBUG
1222  verify(y, tks);
1223  if (DEBUGLEVEL > 0) {
1224  std::cout << "DAClusterizerInZT_vect::vertices :"
1225  << "merging with outlier rejection at T=" << 1 / beta << std::endl;
1226  }
1227  if (DEBUGLEVEL > 2)
1228  dump(beta, y, tks, 2);
1229 #endif
1230 
1231  // merge again (some cluster split by outliers collapse here)
1232  zorder(y);
1233  while (merge(y, tks, beta)) {
1234  set_vtx_range(beta, tks, y);
1235  update(beta, tks, y, rho0);
1236  }
1237 
1238 #ifdef DEBUG
1239  if (DEBUGLEVEL > 0) {
1240  std::cout << "DAClusterizerInZT_vect::vertices :"
1241  << "after merging with outlier rejection at T=" << 1 / beta << std::endl;
1242  }
1243  if (DEBUGLEVEL > 1)
1244  dump(beta, y, tks, 2);
1245 #endif
1246 
1247  // go down to the purging temperature (if it is lower than tmin)
1248  while (beta < betapurge_) {
1250  thermalize(beta, tks, y, delta_lowT_, rho0);
1251  }
1252 
1253 #ifdef DEBUG
1254  verify(y, tks);
1255  if (DEBUGLEVEL > 0) {
1256  std::cout << "DAClusterizerInZT_vect::vertices :"
1257  << "purging at T=" << 1 / beta << std::endl;
1258  }
1259 #endif
1260 
1261  // eliminate insignificant vertices, this is more restrictive at higher T
1262  while (purge(y, tks, rho0, beta)) {
1263  thermalize(beta, tks, y, delta_lowT_, rho0);
1264  if (zorder(y)) {
1265  set_vtx_range(beta, tks, y);
1266  };
1267  }
1268 
1269 #ifdef DEBUG
1270  verify(y, tks);
1271  if (DEBUGLEVEL > 0) {
1272  std::cout << "DAClusterizerInZT_vect::vertices :"
1273  << "last cooling T=" << 1 / beta << std::endl;
1274  }
1275 #endif
1276 
1277  // optionally cool some more without doing anything, to make the assignment harder
1278  while (beta < betastop_) {
1280  thermalize(beta, tks, y, delta_lowT_, rho0);
1281  if (zorder(y)) {
1282  set_vtx_range(beta, tks, y);
1283  };
1284  }
1285 
1286 #ifdef DEBUG
1287  verify(y, tks);
1288  if (DEBUGLEVEL > 0) {
1289  std::cout << "DAClusterizerInZT_vect::vertices :"
1290  << "stop cooling at T=" << 1 / beta << std::endl;
1291  }
1292  if (DEBUGLEVEL > 2)
1293  dump(beta, y, tks, 2);
1294 #endif
1295 
1296  // new, merge here and not in "clusterize"
1297  // final merging step
1298  double betadummy = 1;
1299  while (merge(y, tks, betadummy))
1300  ;
1301 
1302  // select significant tracks and use a TransientVertex as a container
1303  zorder(y);
1304  set_vtx_range(beta, tks, y);
1305  const unsigned int nv = y.getSize();
1306  for (unsigned int k = 0; k < nv; k++)
1307  if (edm::isNotFinite(y.rho[k]) || edm::isNotFinite(y.zvtx[k])) {
1308  y.rho[k] = 0;
1309  y.zvtx[k] = 0;
1310  }
1311 
1312  const auto z_sum_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
1313  std::vector<std::vector<unsigned int> > vtx_track_indices(nv);
1314  for (unsigned int i = 0; i < nt; i++) {
1315  const auto kmin = tks.kmin[i];
1316  const auto kmax = tks.kmax[i];
1317  for (auto k = kmin; k < kmax; k++) {
1318  y.exp_arg[k] = -beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i], tks.tpca[i], y.tvtx[k], tks.dt2[i]);
1319  }
1320 
1321  local_exp_list_range(y.exp_arg, y.exp, kmin, kmax);
1322 
1323  tks.sum_Z[i] = z_sum_init;
1324  for (auto k = kmin; k < kmax; k++) {
1325  tks.sum_Z[i] += y.rho[k] * y.exp[k];
1326  }
1327  const double invZ = tks.sum_Z[i] > 1e-100 ? 1. / tks.sum_Z[i] : 0.0;
1328 
1329  for (auto k = kmin; k < kmax; k++) {
1330  double p = y.rho[k] * y.exp[k] * invZ;
1331  if (p > mintrkweight_) {
1332  // assign track i -> vertex k (hard, mintrkweight_ should be >= 0.5 here
1333  vtx_track_indices[k].push_back(i);
1334  break;
1335  }
1336  }
1337 
1338  } // track loop
1339 
1340  GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01);
1341  for (unsigned int k = 0; k < nv; k++) {
1342  if (!vtx_track_indices[k].empty()) {
1343  GlobalPoint pos(0, 0, y.zvtx[k]);
1344  vector<reco::TransientTrack> vertexTracks;
1345  for (auto i : vtx_track_indices[k]) {
1346  vertexTracks.push_back(*(tks.tt[i]));
1347  }
1348  TransientVertex v(pos, dummyError, vertexTracks, 0);
1349  clusters.push_back(v);
1350  }
1351  }
1352 
1353  return clusters;
1354 }
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
bool purge(vertex_t &, track_t &, double &, const double) const
double update(double beta, track_t &gtracks, vertex_t &gvertices, const double rho0=0, const bool updateTc=false) const
bool split(const double beta, track_t &t, vertex_t &y, double threshold=1.) const
double beta0(const double betamax, track_t const &tks, vertex_t const &y) const
T sqrt(T t)
Definition: SSEVec.h:23
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const
bool zorder(vertex_t &y) const
bool merge(vertex_t &, track_t &, double &beta) const
track_t fill(const std::vector< reco::TransientTrack > &tracks) const
void verify(const vertex_t &v, const track_t &tks, unsigned int nv=999999, unsigned int nt=999999) const
int nt
Definition: AMPTWrapper.h:42
void dump(const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0) const
double a
Definition: hdecay.h:121
unsigned int thermalize(double beta, track_t &gtracks, vertex_t &gvertices, const double delta_max, const double rho0=0.) const
void clear_vtx_range(track_t &gtracks, vertex_t &gvertices) const

◆ zorder()

bool DAClusterizerInZT_vect::zorder ( vertex_t y) const

Definition at line 548 of file DAClusterizerInZT_vect.cc.

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

548  {
549  const unsigned int nv = y.getSize();
550 
551 #ifdef DEBUG
552  assert(y.zvtx_vec.size() == nv);
553  assert(y.tvtx_vec.size() == nv);
554  assert(y.rho_vec.size() == nv);
555 #endif
556 
557  if (nv < 2)
558  return false;
559 
560  bool reordering = true;
561  bool modified = false;
562 
563  while (reordering) {
564  reordering = false;
565  for (unsigned int k = 0; (k + 1) < nv; k++) {
566  if (y.zvtx[k + 1] < y.zvtx[k]) {
567  auto ztemp = y.zvtx[k];
568  y.zvtx[k] = y.zvtx[k + 1];
569  y.zvtx[k + 1] = ztemp;
570  auto ttemp = y.tvtx[k];
571  y.tvtx[k] = y.tvtx[k + 1];
572  y.tvtx[k + 1] = ttemp;
573 #ifdef USEVTXDT2
574  auto dt2temp = y.dt2[k];
575  y.dt2[k] = y.dt2[k + 1];
576  y.dt2[k + 1] = dt2temp;
577 #endif
578  auto ptemp = y.rho[k];
579  y.rho[k] = y.rho[k + 1];
580  y.rho[k + 1] = ptemp;
581  reordering = true;
582  }
583  }
584  modified |= reordering;
585  }
586 
587  if (modified) {
588  y.extractRaw();
589  return true;
590  }
591 
592  return false;
593 }
assert(be >=bs)

Member Data Documentation

◆ betamax_

double DAClusterizerInZT_vect::betamax_
private

Definition at line 302 of file DAClusterizerInZT_vect.h.

◆ betapurge_

double DAClusterizerInZT_vect::betapurge_
private

Definition at line 314 of file DAClusterizerInZT_vect.h.

◆ betastop_

double DAClusterizerInZT_vect::betastop_
private

Definition at line 303 of file DAClusterizerInZT_vect.h.

◆ convergence_mode_

unsigned int DAClusterizerInZT_vect::convergence_mode_
private

Definition at line 316 of file DAClusterizerInZT_vect.h.

◆ coolingFactor_

double DAClusterizerInZT_vect::coolingFactor_
private

Definition at line 301 of file DAClusterizerInZT_vect.h.

◆ d0CutOff_

double DAClusterizerInZT_vect::d0CutOff_
private

Definition at line 305 of file DAClusterizerInZT_vect.h.

◆ delta_highT_

double DAClusterizerInZT_vect::delta_highT_
private

Definition at line 317 of file DAClusterizerInZT_vect.h.

◆ delta_lowT_

double DAClusterizerInZT_vect::delta_lowT_
private

Definition at line 318 of file DAClusterizerInZT_vect.h.

◆ dtCutOff_

double DAClusterizerInZT_vect::dtCutOff_
private

Definition at line 306 of file DAClusterizerInZT_vect.h.

◆ dzCutOff_

double DAClusterizerInZT_vect::dzCutOff_
private

Definition at line 304 of file DAClusterizerInZT_vect.h.

◆ maxIterations_

unsigned int DAClusterizerInZT_vect::maxIterations_
private

Definition at line 300 of file DAClusterizerInZT_vect.h.

◆ mintrkweight_

double DAClusterizerInZT_vect::mintrkweight_
private

Definition at line 309 of file DAClusterizerInZT_vect.h.

◆ sel_zrange_

double DAClusterizerInZT_vect::sel_zrange_
private

Definition at line 320 of file DAClusterizerInZT_vect.h.

◆ t0Max_

double DAClusterizerInZT_vect::t0Max_
private

Definition at line 307 of file DAClusterizerInZT_vect.h.

◆ tmerge_

double DAClusterizerInZT_vect::tmerge_
private

Definition at line 313 of file DAClusterizerInZT_vect.h.

◆ uniquetrkminp_

double DAClusterizerInZT_vect::uniquetrkminp_
private

Definition at line 311 of file DAClusterizerInZT_vect.h.

◆ uniquetrkweight_

double DAClusterizerInZT_vect::uniquetrkweight_
private

Definition at line 310 of file DAClusterizerInZT_vect.h.

◆ vertexSize_

double DAClusterizerInZT_vect::vertexSize_
private

Definition at line 298 of file DAClusterizerInZT_vect.h.

◆ vertexSizeTime_

double DAClusterizerInZT_vect::vertexSizeTime_
private

Definition at line 299 of file DAClusterizerInZT_vect.h.

◆ zdumpcenter_

double DAClusterizerInZT_vect::zdumpcenter_
private

Definition at line 295 of file DAClusterizerInZT_vect.h.

◆ zdumpwidth_

double DAClusterizerInZT_vect::zdumpwidth_
private

Definition at line 296 of file DAClusterizerInZT_vect.h.

◆ zmerge_

double DAClusterizerInZT_vect::zmerge_
private

Definition at line 312 of file DAClusterizerInZT_vect.h.

◆ zrange_min_

const double DAClusterizerInZT_vect::zrange_min_ = 0.1
private

Definition at line 321 of file DAClusterizerInZT_vect.h.