CMS 3D CMS Logo

List of all members | Classes | Public Member Functions | Static Public Member Functions | Private Attributes
DAClusterizerInZ_vect Class Referencefinal

#include <DAClusterizerInZ_vect.h>

Inheritance diagram for DAClusterizerInZ_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
 
 DAClusterizerInZ_vect (const edm::ParameterSet &conf)
 
void dump (const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0, const double rho0=0.) const
 
track_t fill (const std::vector< reco::TransientTrack > &tracks) const
 
bool merge (vertex_t &y, track_t &tks, 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
 
std::vector< TransientVertexvertices_in_blocks (const std::vector< reco::TransientTrack > &tracks) const
 
- Public Member Functions inherited from TrackClusterizerInZ
 TrackClusterizerInZ ()
 
 TrackClusterizerInZ (const edm::ParameterSet &conf)
 
virtual ~TrackClusterizerInZ ()
 

Static Public Member Functions

static void fillPSetDescription (edm::ParameterSetDescription &desc)
 

Private Attributes

double betamax_
 
double betapurge_
 
double betastop_
 
unsigned int block_size_
 
unsigned int convergence_mode_
 
double coolingFactor_
 
double d0CutOff_
 
double delta_highT_
 
double delta_lowT_
 
double dzCutOff_
 
unsigned int maxIterations_
 
double mintrkweight_
 
double overlap_frac_
 
bool runInBlocks_
 
double sel_zrange_
 
double uniquetrkminp_
 
double uniquetrkweight_
 
double vertexSize_
 
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 21 of file DAClusterizerInZ_vect.h.

Constructor & Destructor Documentation

◆ DAClusterizerInZ_vect()

DAClusterizerInZ_vect::DAClusterizerInZ_vect ( const edm::ParameterSet conf)

Definition at line 20 of file DAClusterizerInZ_vect.cc.

References gather_cfg::cout, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), SiStripPI::max, HLT_2023v12_cff::Tmin, HLT_2023v12_cff::Tpurge, and HLT_2023v12_cff::Tstop.

20  {
21  // hardcoded parameters
22  maxIterations_ = 1000;
23  mintrkweight_ = 0.5;
24 
25  // configurable debug output
26 #ifdef DEBUG
27  zdumpcenter_ = conf.getUntrackedParameter<double>("zdumpcenter", 0.);
28  zdumpwidth_ = conf.getUntrackedParameter<double>("zdumpwidth", 20.);
29 #endif
30 
31  // configurable parameters
32  double Tmin = conf.getParameter<double>("Tmin");
33  double Tpurge = conf.getParameter<double>("Tpurge");
34  double Tstop = conf.getParameter<double>("Tstop");
35  vertexSize_ = conf.getParameter<double>("vertexSize");
36  coolingFactor_ = conf.getParameter<double>("coolingFactor");
37  d0CutOff_ = conf.getParameter<double>("d0CutOff");
38  dzCutOff_ = conf.getParameter<double>("dzCutOff");
39  uniquetrkweight_ = conf.getParameter<double>("uniquetrkweight");
40  uniquetrkminp_ = conf.getParameter<double>("uniquetrkminp");
41  zmerge_ = conf.getParameter<double>("zmerge");
42  sel_zrange_ = conf.getParameter<double>("zrange");
43  convergence_mode_ = conf.getParameter<int>("convergence_mode");
44  delta_lowT_ = conf.getParameter<double>("delta_lowT");
45  delta_highT_ = conf.getParameter<double>("delta_highT");
46  runInBlocks_ = conf.getParameter<bool>("runInBlocks");
47  block_size_ = conf.getParameter<unsigned int>("block_size");
48  overlap_frac_ = conf.getParameter<double>("overlap_frac");
49 
50 #ifdef DEBUG
51  std::cout << "DAClusterizerinZ_vect: mintrkweight = " << mintrkweight << std::endl;
52  std::cout << "DAClusterizerinZ_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl;
53  std::cout << "DAClusterizerInZ_vect: uniquetrkminp = " << uniquetrkminp_ << std::endl;
54  std::cout << "DAClusterizerinZ_vect: zmerge = " << zmerge_ << std::endl;
55  std::cout << "DAClusterizerinZ_vect: Tmin = " << Tmin << std::endl;
56  std::cout << "DAClusterizerinZ_vect: Tpurge = " << Tpurge << std::endl;
57  std::cout << "DAClusterizerinZ_vect: Tstop = " << Tstop << std::endl;
58  std::cout << "DAClusterizerinZ_vect: vertexSize = " << vertexSize_ << std::endl;
59  std::cout << "DAClusterizerinZ_vect: coolingFactor = " << coolingFactor_ << std::endl;
60  std::cout << "DAClusterizerinZ_vect: d0CutOff = " << d0CutOff_ << std::endl;
61  std::cout << "DAClusterizerinZ_vect: dzCutOff = " << dzCutOff_ << std::endl;
62  std::cout << "DAClusterizerInZ_vect: zrange = " << sel_zrange_ << std::endl;
63  std::cout << "DAClusterizerinZ_vect: convergence mode = " << convergence_mode_ << std::endl;
64  std::cout << "DAClusterizerinZ_vect: delta_highT = " << delta_highT_ << std::endl;
65  std::cout << "DAClusterizerinZ_vect: delta_lowT = " << delta_lowT_ << std::endl;
66  std::cout << "DAClusterizerinZ_vect: DEBUGLEVEL " << DEBUGLEVEL << std::endl;
67 #endif
68 
69  if (convergence_mode_ > 1) {
70  edm::LogWarning("DAClusterizerinZ_vect")
71  << "DAClusterizerInZ_vect: invalid convergence_mode " << convergence_mode_ << " reset to default " << 0;
73  }
74 
75  if (Tmin == 0) {
76  betamax_ = 1.0;
77  edm::LogWarning("DAClusterizerinZ_vect")
78  << "DAClusterizerInZ_vect: invalid Tmin " << Tmin << " reset to default " << 1. / betamax_;
79  } else {
80  betamax_ = 1. / Tmin;
81  }
82 
83  if ((Tpurge > Tmin) || (Tpurge == 0)) {
84  edm::LogWarning("DAClusterizerinZ_vect")
85  << "DAClusterizerInZ_vect: invalid Tpurge " << Tpurge << " set to " << Tmin;
86  Tpurge = Tmin;
87  }
88  betapurge_ = 1. / Tpurge;
89 
90  if ((Tstop > Tpurge) || (Tstop == 0)) {
91  edm::LogWarning("DAClusterizerinZ_vect")
92  << "DAClusterizerInZ_vect: invalid Tstop " << Tstop << " set to " << max(1., Tpurge);
93  Tstop = max(1., Tpurge);
94  }
95  betastop_ = 1. / Tstop;
96 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
T getUntrackedParameter(std::string const &, T const &) const
Log< level::Warning, false > LogWarning

Member Function Documentation

◆ beta0()

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

Definition at line 593 of file DAClusterizerInZ_vect.cc.

References a, b, gather_cfg::cout, PVValHelper::dx, DAClusterizerInZ_vect::track_t::dz2, DAClusterizerInZ_vect::track_t::getSize(), mps_fire::i, createfilelist::int, dqmdumpme::k, dqm-mbProfile::log, nt, funct::pow(), DAClusterizerInZ_vect::track_t::tkwt, w(), and DAClusterizerInZ_vect::track_t::zpca.

593  {
594  double T0 = 0; // max Tc for beta=0
595  // estimate critical temperature from beta=0 (T=inf)
596  const unsigned int nt = tks.getSize();
597  const unsigned int nv = y.getSize();
598 
599  for (unsigned int k = 0; k < nv; k++) {
600  // vertex fit at T=inf
601  double sumwz = 0;
602  double sumw = 0;
603  for (unsigned int i = 0; i < nt; i++) {
604  double w = tks.tkwt[i] * tks.dz2[i];
605  sumwz += w * tks.zpca[i];
606  sumw += w;
607  }
608 
609  y.zvtx[k] = sumwz / sumw;
610 
611  // estimate Tcrit
612  double a = 0, b = 0;
613  for (unsigned int i = 0; i < nt; i++) {
614  double dx = tks.zpca[i] - y.zvtx[k];
615  double w = tks.tkwt[i] * tks.dz2[i];
616  a += w * std::pow(dx, 2) * tks.dz2[i];
617  b += w;
618  }
619  double Tc = 2. * a / b; // the critical temperature of this vertex
620 
621  if (Tc > T0)
622  T0 = Tc;
623 
624  } // vertex loop (normally there should be only one vertex at beta=0)
625 
626 #ifdef DEBUG
627  if (DEBUGLEVEL > 0) {
628  std::cout << "DAClusterizerInZ_vect.beta0: Tc = " << T0 << std::endl;
629  int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
630  std::cout << "DAClusterizerInZ_vect.beta0: nstep = " << coolingsteps << std::endl;
631  }
632 #endif
633 
634  if (T0 > 1. / betamax) {
635  int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
636 
637  return betamax * std::pow(coolingFactor_, coolingsteps);
638  } else {
639  // ensure at least one annealing step
640  return betamax * coolingFactor_;
641  }
642 }
T w() const
int nt
Definition: AMPTWrapper.h:42
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ clear_vtx_range()

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

Definition at line 272 of file DAClusterizerInZ_vect.cc.

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

272  {
273  const unsigned int nt = gtracks.getSize();
274  const unsigned int nv = gvertices.getSize();
275  for (auto itrack = 0U; itrack < nt; ++itrack) {
276  gtracks.kmin[itrack] = 0;
277  gtracks.kmax[itrack] = nv;
278  }
279 }
int nt
Definition: AMPTWrapper.h:42

◆ clusterize()

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

Implements TrackClusterizerInZ.

Definition at line 1275 of file DAClusterizerInZ_vect.cc.

References funct::abs(), bsc_activity_cfg::clusters, gather_cfg::cout, mps_fire::i, dqmdumpme::k, eostools::move(), position, MetAnalyzer::pv(), pwdgSkimBPark_cfi::tracks, and AlignmentTracksFromVertexSelector_cfi::vertices.

1276  {
1277  vector<vector<reco::TransientTrack>> clusters;
1278 
1279  vector<TransientVertex> pv;
1280  if (runInBlocks_ and (block_size_ < tracks.size())) //doesn't bother if low number of tracks
1282  else
1283  pv = vertices(tracks);
1284 
1285 #ifdef DEBUG
1286  if (DEBUGLEVEL > 0) {
1287  std::cout << "###################################################" << endl;
1288  std::cout << "# vectorized DAClusterizerInZ_vect::clusterize nt=" << tracks.size() << endl;
1289  std::cout << "# DAClusterizerInZ_vect::clusterize pv.size=" << pv.size() << endl;
1290  std::cout << "###################################################" << endl;
1291  }
1292 #endif
1293 
1294  if (pv.empty()) {
1295  return clusters;
1296  }
1297 
1298  // fill into clusters and merge
1299  vector<reco::TransientTrack> aCluster = pv.begin()->originalTracks();
1300 
1301  for (auto k = pv.begin() + 1; k != pv.end(); k++) {
1302  if (std::abs(k->position().z() - (k - 1)->position().z()) > (2 * vertexSize_)) {
1303  // close a cluster
1304  if (aCluster.size() > 1) {
1305  clusters.push_back(aCluster);
1306  }
1307 #ifdef DEBUG
1308  else {
1309  std::cout << " one track cluster at " << k->position().z() << " suppressed" << std::endl;
1310  }
1311 #endif
1312  aCluster.clear();
1313  }
1314  for (unsigned int i = 0; i < k->originalTracks().size(); i++) {
1315  aCluster.push_back(k->originalTracks()[i]);
1316  }
1317  }
1318  clusters.emplace_back(std::move(aCluster));
1319 
1320  return clusters;
1321 }
std::vector< TransientVertex > vertices(const std::vector< reco::TransientTrack > &tracks) const
def pv(vc)
Definition: MetAnalyzer.py:7
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< TransientVertex > vertices_in_blocks(const std::vector< reco::TransientTrack > &tracks) const
static int position[264][3]
Definition: ReadPGInfo.cc:289
def move(src, dest)
Definition: eostools.py:511

◆ dump()

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

Definition at line 1323 of file DAClusterizerInZ_vect.cc.

References a, b, HLT_2023v12_cff::beta, gather_cfg::cout, TauDecayModes::dec, DAClusterizerInZ_vect::track_t::dz2, F(), alignBH_cfg::fixed, DAClusterizerInZ_vect::track_t::getSize(), reco::TrackBase::highPurity, mps_fire::i, listHistos::IP, dqmiolumiharvest::j, DAClusterizerInZ_vect::track_t::kmax, DAClusterizerInZ_vect::track_t::kmin, dqm-mbProfile::log, reco::HitPattern::MISSING_OUTER_HITS, nt, AlCaHLTBitMon_ParallelJobs::p, GeomDetEnumerators::PixelBarrel, jetUpdater_cfi::sort, mathSSE::sqrt(), DAClusterizerInZ_vect::track_t::sum_Z, DAClusterizerInZ_vect::vertex_t::sw, DAClusterizerInZ_vect::vertex_t::swE, DAClusterizerInZ_vect::track_t::tkwt, DAClusterizerInZ_vect::track_t::tt, update, verbosity, and DAClusterizerInZ_vect::track_t::zpca.

1324  {
1325 #ifdef DEBUG
1326  const unsigned int nv = y.getSize();
1327  const unsigned int nt = tks.getSize();
1328  // make a local copies of clusters and tracks to update Tc [ = y_local.swE / y_local.sw ]
1329  vertex_t y_local = y;
1330  track_t tks_local = tks;
1331  update(beta, tks_local, y_local, rho0, true);
1332 
1333  std::vector<unsigned int> iz;
1334  for (unsigned int j = 0; j < nt; j++) {
1335  iz.push_back(j);
1336  }
1337  std::sort(iz.begin(), iz.end(), [tks](unsigned int a, unsigned int b) { return tks.zpca[a] < tks.zpca[b]; });
1338  std::cout << std::endl;
1339  std::cout << "-----DAClusterizerInZ::dump ----" << nv << " clusters " << std::endl;
1340  std::cout << " ";
1341  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1342  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1343  std::cout << " " << setw(3) << ivertex << " ";
1344  }
1345  }
1346  std::cout << endl;
1347  std::cout << " z= ";
1348  std::cout << setprecision(4);
1349  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1350  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1351  std::cout << setw(8) << fixed << y.zvtx[ivertex];
1352  }
1353  }
1354  std::cout << endl
1355  << "T=" << setw(15) << 1. / beta << " Tmin =" << setw(10) << 1. / betamax_
1356  << " Tc= ";
1357  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1358  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1359  double Tc = 2 * y_local.swE[ivertex] / y_local.sw[ivertex];
1360  std::cout << setw(8) << fixed << setprecision(1) << Tc;
1361  }
1362  }
1363  std::cout << endl;
1364 
1365  std::cout << " pk= ";
1366  double sumpk = 0;
1367  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1368  sumpk += y.rho[ivertex];
1369  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1370  continue;
1371  std::cout << setw(8) << setprecision(4) << fixed << y.rho[ivertex];
1372  }
1373  std::cout << endl;
1374 
1375  std::cout << " nt= ";
1376  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1377  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1378  continue;
1379  std::cout << setw(8) << setprecision(1) << fixed << y.rho[ivertex] * nt;
1380  }
1381  std::cout << endl;
1382 
1383  if (verbosity > 0) {
1384  double E = 0, F = 0;
1385  std::cout << endl;
1386  std::cout << "---- z +/- dz ip +/-dip pt phi eta weights ----" << endl;
1387  std::cout << setprecision(4);
1388  for (unsigned int i0 = 0; i0 < nt; i0++) {
1389  unsigned int i = iz[i0];
1390  if (tks.sum_Z[i] > 0) {
1391  F -= std::log(tks.sum_Z[i]) / beta;
1392  }
1393  double tz = tks.zpca[i];
1394 
1395  if (std::fabs(tz - zdumpcenter_) > zdumpwidth_)
1396  continue;
1397  std::cout << setw(4) << i << ")" << setw(8) << fixed << setprecision(4) << tz << " +/-" << setw(6)
1398  << sqrt(1. / tks.dz2[i]);
1399  if ((tks.tt[i] == nullptr)) {
1400  std::cout << " effective track ";
1401  } else {
1402  if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) {
1403  std::cout << " *";
1404  } else {
1405  std::cout << " ";
1406  }
1407  if (tks.tt[i]->track().hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) {
1408  std::cout << "+";
1409  } else {
1410  std::cout << "-";
1411  }
1412  std::cout << setw(1)
1413  << tks.tt[i]
1414  ->track()
1415  .hitPattern()
1416  .pixelBarrelLayersWithMeasurement(); // see DataFormats/TrackReco/interface/HitPattern.h
1417  std::cout << setw(1) << tks.tt[i]->track().hitPattern().pixelEndcapLayersWithMeasurement();
1418  std::cout << setw(1) << hex
1419  << tks.tt[i]->track().hitPattern().trackerLayersWithMeasurement() -
1420  tks.tt[i]->track().hitPattern().pixelLayersWithMeasurement()
1421  << dec;
1422  std::cout << "=" << setw(1) << hex
1423  << tks.tt[i]->track().hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << dec;
1424 
1425  Measurement1D IP = tks.tt[i]->stateAtBeamLine().transverseImpactParameter();
1426  std::cout << setw(8) << IP.value() << "+/-" << setw(6) << IP.error();
1427  std::cout << " " << setw(6) << setprecision(2) << tks.tt[i]->track().pt() * tks.tt[i]->track().charge();
1428  std::cout << " " << setw(5) << setprecision(2) << tks.tt[i]->track().phi() << " " << setw(5) << setprecision(2)
1429  << tks.tt[i]->track().eta();
1430  } // not a dummy track
1431 
1432  double sump = 0.;
1433  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1434  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1435  continue;
1436 
1437  if ((tks.tkwt[i] > 0) && (tks.sum_Z[i] > 0)) {
1438  //double p=pik(beta,tks[i],*k);
1439  double p = y.rho[ivertex] * local_exp(-beta * Eik(tks.zpca[i], y.zvtx[ivertex], tks.dz2[i])) / tks.sum_Z[i];
1440  if (p > 0.0001) {
1441  std::cout << setw(8) << setprecision(3) << p;
1442  } else {
1443  std::cout << " . ";
1444  }
1445  E += p * Eik(tks.zpca[i], y.zvtx[ivertex], tks.dz2[i]);
1446  sump += p;
1447  } else {
1448  std::cout << " ";
1449  }
1450  }
1451  std::cout << " ( " << std::setw(3) << tks.kmin[i] << "," << std::setw(3) << tks.kmax[i] - 1 << " ) ";
1452  std::cout << endl;
1453  }
1454  std::cout << " ";
1455  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1456  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1457  std::cout << " " << setw(3) << ivertex << " ";
1458  }
1459  }
1460  std::cout << endl;
1461  std::cout << " z= ";
1462  std::cout << setprecision(4);
1463  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1464  if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1465  std::cout << setw(8) << fixed << y.zvtx[ivertex];
1466  }
1467  }
1468  std::cout << endl;
1469  std::cout << endl
1470  << "T=" << 1 / beta << " E=" << E << " n=" << y.getSize() << " F= " << F << endl
1471  << "----------" << endl;
1472  }
1473 #endif
1474 }
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:19
int nt
Definition: AMPTWrapper.h:42
const int verbosity
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:163

◆ fill()

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

Definition at line 173 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::track_t::addItemSorted(), HLT_2023v12_cff::atIP, pwdgSkimBPark_cfi::beamSpot, gather_cfg::cout, geometryDiff::epsilon, DAClusterizerInZ_vect::track_t::extractRaw(), DAClusterizerInZ_vect::track_t::getSize(), edm::isNotFinite(), SiStripPI::min, DAClusterizerInZ_vect::track_t::osumtkwt, position, funct::pow(), AlignmentPI::t_z, and pwdgSkimBPark_cfi::tracks.

173  {
174  // prepare track data for clustering
175  track_t tks;
176  double sumtkwt = 0.;
177  for (auto it = tracks.begin(); it != tracks.end(); it++) {
178  if (!(*it).isValid())
179  continue;
180  double t_tkwt = 1.;
181  double t_z = ((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
182  if (std::fabs(t_z) > 1000.)
183  continue;
184  auto const& t_mom = (*it).stateAtBeamLine().trackStateAtPCA().momentum();
185  // get the beam-spot
186  reco::BeamSpot beamspot = (it->stateAtBeamLine()).beamSpot();
187  double t_dz2 = std::pow((*it).track().dzError(), 2) // track errror
188  + (std::pow(beamspot.BeamWidthX() * t_mom.x(), 2) + std::pow(beamspot.BeamWidthY() * t_mom.y(), 2)) *
189  std::pow(t_mom.z(), 2) / std::pow(t_mom.perp2(), 2) // beam spot width
190  + std::pow(vertexSize_, 2); // intrinsic vertex size, safer for outliers and short lived decays
191  t_dz2 = 1. / t_dz2;
192  if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits<double>::min())
193  continue;
194  if (d0CutOff_ > 0) {
195  Measurement1D atIP = (*it).stateAtBeamLine().transverseImpactParameter(); // error contains beamspot
196  t_tkwt = 1. / (1. + local_exp(std::pow(atIP.value() / atIP.error(), 2) -
197  std::pow(d0CutOff_, 2))); // reduce weight for high ip tracks
199  edm::LogWarning("DAClusterizerinZ_vect") << "rejected track t_tkwt " << t_tkwt;
200  continue; // usually is > 0.99
201  }
202  }
203  tks.addItemSorted(t_z, t_dz2, &(*it), t_tkwt);
204  sumtkwt += t_tkwt;
205  }
206 
207  tks.extractRaw();
208  tks.osumtkwt = sumtkwt > 0 ? 1. / sumtkwt : 0.;
209 
210 #ifdef DEBUG
211  if (DEBUGLEVEL > 0) {
212  std::cout << "Track count (Z) " << tks.getSize() << std::endl;
213  }
214 #endif
215 
216  return tks;
217 }
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
static int position[264][3]
Definition: ReadPGInfo.cc:289
Log< level::Warning, false > LogWarning
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ fillPSetDescription()

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

Definition at line 1476 of file DAClusterizerInZ_vect.cc.

References submitPVResolutionJobs::desc.

Referenced by PrimaryVertexValidation::fillDescriptions(), and DAClusterizerInZT_vect::fillPSetDescription().

1476  {
1477  desc.addUntracked<double>("zdumpcenter", 0.);
1478  desc.addUntracked<double>("zdumpwidth", 20.);
1479  desc.add<double>("d0CutOff", 3.0);
1480  desc.add<double>("Tmin", 2.0);
1481  desc.add<double>("delta_lowT", 0.001);
1482  desc.add<double>("zmerge", 0.01);
1483  desc.add<double>("dzCutOff", 3.0);
1484  desc.add<double>("Tpurge", 2.0);
1485  desc.add<int>("convergence_mode", 0);
1486  desc.add<double>("delta_highT", 0.01);
1487  desc.add<double>("Tstop", 0.5);
1488  desc.add<double>("coolingFactor", 0.6);
1489  desc.add<double>("vertexSize", 0.006);
1490  desc.add<double>("uniquetrkweight", 0.8);
1491  desc.add<double>("uniquetrkminp", 0.0);
1492  desc.add<double>("zrange", 4.0);
1493  desc.add<bool>("runInBlocks", false);
1494  desc.add<unsigned int>("block_size", 512);
1495  desc.add<double>("overlap_frac", 0.5);
1496 }

◆ merge()

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

Definition at line 465 of file DAClusterizerInZ_vect.cc.

References cms::cuda::assert(), HLT_2023v12_cff::beta, gather_cfg::cout, hlt_jetmet_dqm_QT_fromfile_cfg::critical, alignBH_cfg::fixed, and dqmdumpme::k.

465  {
466  // merge clusters that collapsed or never separated,
467  // only merge if the estimated critical temperature of the merged vertex is below the current temperature
468  // return true if vertices were merged, false otherwise
469  const unsigned int nv = y.getSize();
470 
471  if (nv < 2)
472  return false;
473 
474  // merge the smallest distance clusters first
475  std::vector<std::pair<double, unsigned int>> critical;
476  for (unsigned int k = 0; (k + 1) < nv; k++) {
477  if (std::fabs(y.zvtx[k + 1] - y.zvtx[k]) < zmerge_) {
478  critical.push_back(make_pair(std::fabs(y.zvtx[k + 1] - y.zvtx[k]), k));
479  }
480  }
481  if (critical.empty())
482  return false;
483 
484  std::stable_sort(critical.begin(), critical.end(), std::less<std::pair<double, unsigned int>>());
485 
486  for (unsigned int ik = 0; ik < critical.size(); ik++) {
487  unsigned int k = critical[ik].second;
488  double rho = y.rho[k] + y.rho[k + 1];
489 
490 #ifdef DEBUG
491  assert((k + 1) < nv);
492  if (DEBUGLEVEL > 1) {
493  std::cout << "merging " << fixed << setprecision(4) << y.zvtx[k + 1] << " and " << y.zvtx[k]
494  << " sw = " << y.sw[k] + y.sw[k + 1] << std::endl;
495  }
496 #endif
497 
498  if (rho > 0) {
499  y.zvtx[k] = (y.rho[k] * y.zvtx[k] + y.rho[k + 1] * y.zvtx[k + 1]) / rho;
500  } else {
501  y.zvtx[k] = 0.5 * (y.zvtx[k] + y.zvtx[k + 1]);
502  }
503  y.rho[k] = rho;
504  y.sw[k] += y.sw[k + 1];
505  y.removeItem(k + 1, tks);
506  set_vtx_range(beta, tks, y);
507  y.extractRaw();
508  return true;
509  }
510 
511  return false;
512 }
assert(be >=bs)
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const

◆ purge()

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

Definition at line 514 of file DAClusterizerInZ_vect.cc.

References cms::cuda::assert(), HLT_2023v12_cff::beta, gather_cfg::cout, DAClusterizerInZ_vect::track_t::dz2, HLT_2023v12_cff::eps, cms::cuda::for(), DAClusterizerInZ_vect::track_t::getSize(), mps_fire::i, dqmdumpme::k, reco::ParticleMasses::k0, DAClusterizerInZ_vect::track_t::kmax, DAClusterizerInZ_vect::track_t::kmin, nt, AlCaHLTBitMon_ParallelJobs::p, DAClusterizerInZ_vect::track_t::sum_Z, DAClusterizerInZ_vect::track_t::tkwt, and DAClusterizerInZ_vect::track_t::zpca.

514  {
515  constexpr double eps = 1.e-100;
516  // eliminate clusters with only one significant/unique track
517  const unsigned int nv = y.getSize();
518  const unsigned int nt = tks.getSize();
519 
520  if (nv < 2)
521  return false;
522 
523  std::vector<double> sump_v(nv), arg_cache_v(nv), exp_cache_v(nv), pcut_cache_v(nv);
524  std::vector<int> nUnique_v(nv);
525  double* __restrict__ parg_cache;
526  double* __restrict__ pexp_cache;
527  double* __restrict__ ppcut_cache;
528  double* __restrict__ psump;
529  int* __restrict__ pnUnique;
530  int constexpr nunique_min_ = 2;
531 
532  set_vtx_range(beta, tks, y);
533 
534  parg_cache = arg_cache_v.data();
535  pexp_cache = exp_cache_v.data();
536  ppcut_cache = pcut_cache_v.data();
537  psump = sump_v.data();
538  pnUnique = nUnique_v.data();
539 
540  const auto rhoconst = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
541  for (unsigned int k = 0; k < nv; k++) {
542  const double pmax = y.rho[k] / (y.rho[k] + rhoconst);
543  ppcut_cache[k] = uniquetrkweight_ * pmax;
544  }
545 
546  for (unsigned int i = 0; i < nt; i++) {
547  const auto invZ = ((tks.sum_Z[i] > eps) && (tks.tkwt[i] > uniquetrkminp_)) ? 1. / tks.sum_Z[i] : 0.;
548  const auto track_z = tks.zpca[i];
549  const auto botrack_dz2 = -beta * tks.dz2[i];
550  const auto kmin = tks.kmin[i];
551  const auto kmax = tks.kmax[i];
552 
553  for (unsigned int k = kmin; k < kmax; k++) {
554  const auto mult_resz = track_z - y.zvtx[k];
555  parg_cache[k] = botrack_dz2 * (mult_resz * mult_resz);
556  }
557 
558  local_exp_list_range(parg_cache, pexp_cache, kmin, kmax);
559 
560  for (unsigned int k = kmin; k < kmax; k++) {
561  const double p = y.rho[k] * pexp_cache[k] * invZ;
562  psump[k] += p;
563  pnUnique[k] += (p > ppcut_cache[k]) ? 1 : 0;
564  }
565  }
566 
567  double sumpmin = nt;
568  unsigned int k0 = nv;
569  for (unsigned k = 0; k < nv; k++) {
570  if ((pnUnique[k] < nunique_min_) && (psump[k] < sumpmin)) {
571  sumpmin = psump[k];
572  k0 = k;
573  }
574  }
575 
576  if (k0 != nv) {
577 #ifdef DEBUG
578  assert(k0 < y.getSize());
579  if (DEBUGLEVEL > 1) {
580  std::cout << "eliminating prototype at " << std::setw(10) << std::setprecision(4) << y.zvtx[k0]
581  << " with sump=" << sumpmin << " rho*nt =" << y.rho[k0] * nt << endl;
582  }
583 #endif
584 
585  y.removeItem(k0, tks);
586  set_vtx_range(beta, tks, y);
587  return true;
588  } else {
589  return false;
590  }
591 }
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
assert(be >=bs)
int nt
Definition: AMPTWrapper.h:42
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const

◆ set_vtx_range()

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

Definition at line 223 of file DAClusterizerInZ_vect.cc.

References HLT_2023v12_cff::beta, DAClusterizerInZ_vect::track_t::dz2, DAClusterizerInZ_vect::track_t::getSize(), DAClusterizerInZ_vect::vertex_t::getSize(), DAClusterizerInZ_vect::track_t::kmax, DAClusterizerInZ_vect::track_t::kmin, SiStripPI::max, SiStripPI::min, nt, mathSSE::sqrt(), mitigatedMETSequence_cff::U, SiStripMonitorCluster_cfi::zmax, SiStripMonitorCluster_cfi::zmin, DAClusterizerInZ_vect::track_t::zpca, TkClusParameters_cff::zrange, and DAClusterizerInZ_vect::vertex_t::zvtx.

223  {
224  const unsigned int nv = gvertices.getSize();
225  const unsigned int nt = gtracks.getSize();
226 
227  if (nv == 0) {
228  edm::LogWarning("DAClusterizerinZ_vect") << "empty cluster list in set_vtx_range";
229  return;
230  }
231 
232  for (auto itrack = 0U; itrack < nt; ++itrack) {
233  double zrange = max(sel_zrange_ / sqrt(beta * gtracks.dz2[itrack]), zrange_min_);
234 
235  double zmin = gtracks.zpca[itrack] - zrange;
236  unsigned int kmin = min(nv - 1, gtracks.kmin[itrack]);
237  // find the smallest vertex_z that is larger than zmin
238  if (gvertices.zvtx[kmin] > zmin) {
239  while ((kmin > 0) && (gvertices.zvtx[kmin - 1] > zmin)) {
240  kmin--;
241  }
242  } else {
243  while ((kmin < (nv - 1)) && (gvertices.zvtx[kmin] < zmin)) {
244  kmin++;
245  }
246  }
247 
248  double zmax = gtracks.zpca[itrack] + zrange;
249  unsigned int kmax = min(nv - 1, gtracks.kmax[itrack] - 1);
250  // note: kmax points to the last vertex in the range, while gtracks.kmax points to the entry BEHIND the last vertex
251  // find the largest vertex_z that is smaller than zmax
252  if (gvertices.zvtx[kmax] < zmax) {
253  while ((kmax < (nv - 1)) && (gvertices.zvtx[kmax + 1] < zmax)) {
254  kmax++;
255  }
256  } else {
257  while ((kmax > 0) && (gvertices.zvtx[kmax] > zmax)) {
258  kmax--;
259  }
260  }
261 
262  if (kmin <= kmax) {
263  gtracks.kmin[itrack] = kmin;
264  gtracks.kmax[itrack] = kmax + 1;
265  } else {
266  gtracks.kmin[itrack] = max(0U, min(kmin, kmax));
267  gtracks.kmax[itrack] = min(nv, max(kmin, kmax) + 1);
268  }
269  }
270 }
T sqrt(T t)
Definition: SSEVec.h:19
int nt
Definition: AMPTWrapper.h:42
Log< level::Warning, false > LogWarning

◆ split()

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

Definition at line 644 of file DAClusterizerInZ_vect.cc.

References cms::cuda::assert(), HLT_2023v12_cff::beta, gather_cfg::cout, hlt_jetmet_dqm_QT_fromfile_cfg::critical, DAClusterizerInZ_vect::track_t::dz2, MillePedeFileConverter_cfg::e, geometryDiff::epsilon, alignBH_cfg::fixed, DAClusterizerInZ_vect::track_t::getSize(), mps_fire::i, dqmdumpme::k, nt, AlCaHLTBitMon_ParallelJobs::p, LaserDQM_cfg::p1, SiStripOfflineCRack_cfg::p2, edm::second(), submitPVValidationJobs::split(), mathSSE::sqrt(), DAClusterizerInZ_vect::track_t::sum_Z, submitPVValidationJobs::t, remoteMonitoring_LASER_era2018_cfg::threshold, DAClusterizerInZ_vect::track_t::tkwt, update, w(), w1, w2, testProducerWithPsetDescEmpty_cfi::z2, and DAClusterizerInZ_vect::track_t::zpca.

644  {
645  // split only critical vertices (Tc >~ T=1/beta <==> beta*Tc>~1)
646  // returns true if at least one cluster was split
647 
648  // update the sums needed for Tc, rho0 is never non-zero while splitting is still active
649  update(beta, tks, y, 0., true); // the "true" flag enables Tc evaluation
650 
651  constexpr double epsilon = 1e-3; // minimum split size
652  unsigned int nv = y.getSize();
653 
654  // avoid left-right biases by splitting highest Tc first
655 
656  std::vector<std::pair<double, unsigned int>> critical;
657  for (unsigned int k = 0; k < nv; k++) {
658  double Tc = 2 * y.swE[k] / y.sw[k];
659  if (beta * Tc > threshold) {
660  critical.push_back(make_pair(Tc, k));
661  }
662  }
663  if (critical.empty())
664  return false;
665 
666  std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int>>());
667 
668  bool split = false;
669  const unsigned int nt = tks.getSize();
670 
671  for (unsigned int ic = 0; ic < critical.size(); ic++) {
672  unsigned int k = critical[ic].second;
673 
674  // estimate subcluster positions and weight
675  double p1 = 0, z1 = 0, w1 = 0;
676  double p2 = 0, z2 = 0, w2 = 0;
677  for (unsigned int i = 0; i < nt; i++) {
678  if (tks.sum_Z[i] > 1.e-100) {
679  // winner-takes-all, usually overestimates splitting
680  double tl = tks.zpca[i] < y.zvtx[k] ? 1. : 0.;
681  double tr = 1. - tl;
682 
683  // soften it, especially at low T
684  double arg = (tks.zpca[i] - y.zvtx[k]) * sqrt(beta * tks.dz2[i]);
685  if (std::fabs(arg) < 20) {
686  double t = local_exp(-arg);
687  tl = t / (t + 1.);
688  tr = 1 / (t + 1.);
689  }
690 
691  double p = y.rho[k] * tks.tkwt[i] * local_exp(-beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i])) / tks.sum_Z[i];
692  double w = p * tks.dz2[i];
693  p1 += p * tl;
694  z1 += w * tl * tks.zpca[i];
695  w1 += w * tl;
696  p2 += p * tr;
697  z2 += w * tr * tks.zpca[i];
698  w2 += w * tr;
699  }
700  }
701 
702  if (w1 > 0) {
703  z1 = z1 / w1;
704  } else {
705  z1 = y.zvtx[k] - epsilon;
706  }
707  if (w2 > 0) {
708  z2 = z2 / w2;
709  } else {
710  z2 = y.zvtx[k] + epsilon;
711  }
712 
713  // reduce split size if there is not enough room
714  if ((k > 0) && (z1 < (0.6 * y.zvtx[k] + 0.4 * y.zvtx[k - 1]))) {
715  z1 = 0.6 * y.zvtx[k] + 0.4 * y.zvtx[k - 1];
716  }
717  if ((k + 1 < nv) && (z2 > (0.6 * y.zvtx[k] + 0.4 * y.zvtx[k + 1]))) {
718  z2 = 0.6 * y.zvtx[k] + 0.4 * y.zvtx[k + 1];
719  }
720 
721 #ifdef DEBUG
722  assert(k < nv);
723  if (DEBUGLEVEL > 1) {
724  if (std::fabs(y.zvtx[k] - zdumpcenter_) < zdumpwidth_) {
725  std::cout << " T= " << std::setw(8) << 1. / beta << " Tc= " << critical[ic].first << " splitting "
726  << std::fixed << std::setprecision(4) << y.zvtx[k] << " --> " << z1 << "," << z2 << " [" << p1
727  << "," << p2 << "]";
728  if (std::fabs(z2 - z1) > epsilon) {
729  std::cout << std::endl;
730  } else {
731  std::cout << " rejected " << std::endl;
732  }
733  }
734  }
735 #endif
736 
737  // split if the new subclusters are significantly separated
738  if ((z2 - z1) > epsilon) {
739  split = true;
740  double pk1 = p1 * y.rho[k] / (p1 + p2);
741  double pk2 = p2 * y.rho[k] / (p1 + p2);
742  y.zvtx[k] = z2;
743  y.rho[k] = pk2;
744  y.insertItem(k, z1, pk1, tks);
745  if (k == 0)
746  y.extractRaw();
747 
748  nv++;
749 
750  // adjust remaining pointers
751  for (unsigned int jc = ic; jc < critical.size(); jc++) {
752  if (critical[jc].second >= k) {
753  critical[jc].second++;
754  }
755  }
756  } else {
757 #ifdef DEBUG
758  std::cout << "warning ! split rejected, too small." << endl;
759 #endif
760  }
761  }
762 
763  return split;
764 }
common ppss p3p6s2 common epss epspn46 common const1 w2
Definition: inclppp.h:1
T w() const
assert(be >=bs)
A arg
Definition: Factorize.h:31
U second(std::pair< T, U > const &p)
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:19
int nt
Definition: AMPTWrapper.h:42
weight_default_t w1[2000]
Definition: w1.h:9
bool split(const double beta, track_t &t, vertex_t &y, double threshold=1.) const

◆ thermalize()

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

Definition at line 417 of file DAClusterizerInZ_vect.cc.

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

418  {
419  unsigned int niter = 0;
420  double delta = 0;
421  double delta_max = delta_lowT_;
422 
423  if (convergence_mode_ == 0) {
424  delta_max = delta_max0;
425  } else if (convergence_mode_ == 1) {
426  delta_max = delta_lowT_ / sqrt(std::max(beta, 1.0));
427  }
428 
429  set_vtx_range(beta, tks, v);
430  double delta_sum_range = 0; // accumulate max(|delta-z|) as a lower bound
431  std::vector<double> z0 = v.zvtx_vec;
432 
433  while (niter++ < maxIterations_) {
434  delta = update(beta, tks, v, rho0);
435  delta_sum_range += delta;
436 
437  if (delta_sum_range > zrange_min_) {
438  for (unsigned int k = 0; k < v.getSize(); k++) {
439  if (std::abs(v.zvtx_vec[k] - z0[k]) > zrange_min_) {
440  set_vtx_range(beta, tks, v);
441  delta_sum_range = 0;
442  z0 = v.zvtx_vec;
443  break;
444  }
445  }
446  }
447 
448  if (delta < delta_max) {
449  break;
450  }
451  }
452 
453 #ifdef DEBUG
454  if (DEBUGLEVEL > 0) {
455  std::cout << "DAClusterizerInZ_vect.thermalize niter = " << niter << " at T = " << 1 / beta
456  << " nv = " << v.getSize() << std::endl;
457  if (DEBUGLEVEL > 2)
458  dump(beta, v, tks, 0, rho0);
459  }
460 #endif
461 
462  return niter;
463 }
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:19
void dump(const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0, const double rho0=0.) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const

◆ update()

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

Definition at line 281 of file DAClusterizerInZ_vect.cc.

References funct::abs(), HLT_2023v12_cff::beta, dumpMFGeometry_cfg::delta, DAClusterizerInZ_vect::vertex_t::exp, DAClusterizerInZ_vect::vertex_t::exp_arg, DAClusterizerInZ_vect::track_t::getSize(), DAClusterizerInZ_vect::vertex_t::getSize(), edm::isNotFinite(), dqmdumpme::k, DAClusterizerInZ_vect::track_t::kmax, DAClusterizerInZ_vect::track_t::kmin, SiStripPI::max, nt, DAClusterizerInZ_vect::track_t::osumtkwt, DAClusterizerInZ_vect::vertex_t::se, DAClusterizerInZ_vect::track_t::sum_Z, DAClusterizerInZ_vect::vertex_t::sw, DAClusterizerInZ_vect::vertex_t::swE, DAClusterizerInZ_vect::vertex_t::swz, pwdgSkimBPark_cfi::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().

282  {
283  // update weights and vertex positions
284  // returns the maximum of changes of vertex positions
285  // sums needed for Tc are only updated if updateTC == true
286 
287  const unsigned int nt = gtracks.getSize();
288  const unsigned int nv = gvertices.getSize();
289  auto osumtkwt = gtracks.osumtkwt;
290 
291  double Z_init = 0;
292  if (rho0 > 0) {
293  Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); // cut-off
294  }
295 
296  // define kernels
297  auto kernel_calc_exp_arg_range = [beta](const unsigned int itrack,
298  track_t const& tracks,
299  vertex_t const& vertices,
300  const unsigned int kmin,
301  const unsigned int kmax) {
302  const double track_z = tracks.zpca[itrack];
303  const double botrack_dz2 = -beta * tracks.dz2[itrack];
304 
305  // auto-vectorized
306  for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
307  auto mult_res = track_z - vertices.zvtx[ivertex];
308  vertices.exp_arg[ivertex] = botrack_dz2 * (mult_res * mult_res);
309  }
310  };
311 
312  auto kernel_add_Z_range = [Z_init](
313  vertex_t const& vertices, const unsigned int kmin, const unsigned int kmax) -> double {
314  double ZTemp = Z_init;
315  for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
316  ZTemp += vertices.rho[ivertex] * vertices.exp[ivertex];
317  }
318  return ZTemp;
319  };
320 
321  auto kernel_calc_normalization_range = [updateTc](const unsigned int track_num,
322  track_t& tracks,
323  vertex_t& vertices,
324  const unsigned int kmin,
325  const unsigned int kmax) {
326  auto o_trk_sum_Z = tracks.tkwt[track_num] / tracks.sum_Z[track_num];
327  auto o_trk_dz2 = tracks.dz2[track_num];
328  auto tmp_trk_z = tracks.zpca[track_num];
329 
330  // auto-vectorized
331  if (updateTc) {
332 #pragma GCC ivdep
333  for (unsigned int k = kmin; k < kmax; ++k) {
334  vertices.se[k] += vertices.exp[k] * o_trk_sum_Z;
335  auto w = vertices.rho[k] * vertices.exp[k] * (o_trk_sum_Z * o_trk_dz2);
336  vertices.sw[k] += w;
337  vertices.swz[k] += w * tmp_trk_z;
338  vertices.swE[k] += w * vertices.exp_arg[k];
339  }
340  } else {
341  // same loop but without updating sWE
342 #pragma GCC ivdep
343  for (unsigned int k = kmin; k < kmax; ++k) {
344  vertices.se[k] += vertices.exp[k] * o_trk_sum_Z;
345  auto w = vertices.rho[k] * vertices.exp[k] * (o_trk_sum_Z * o_trk_dz2);
346  vertices.sw[k] += w;
347  vertices.swz[k] += w * tmp_trk_z;
348  }
349  }
350  };
351 
352  if (updateTc) {
353  for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
354  gvertices.se[ivertex] = 0.0;
355  gvertices.sw[ivertex] = 0.0;
356  gvertices.swz[ivertex] = 0.0;
357  gvertices.swE[ivertex] = 0.0;
358  }
359  } else {
360  for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
361  gvertices.se[ivertex] = 0.0;
362  gvertices.sw[ivertex] = 0.0;
363  gvertices.swz[ivertex] = 0.0;
364  }
365  }
366 
367  // loop over tracks
368  for (auto itrack = 0U; itrack < nt; ++itrack) {
369  const unsigned int kmin = gtracks.kmin[itrack];
370  const unsigned int kmax = gtracks.kmax[itrack];
371 
372  kernel_calc_exp_arg_range(itrack, gtracks, gvertices, kmin, kmax);
373  local_exp_list_range(gvertices.exp_arg, gvertices.exp, kmin, kmax);
374  gtracks.sum_Z[itrack] = kernel_add_Z_range(gvertices, kmin, kmax);
375 
376  if (edm::isNotFinite(gtracks.sum_Z[itrack]))
377  gtracks.sum_Z[itrack] = 0.0;
378 
379  if (gtracks.sum_Z[itrack] > 1.e-100) {
380  kernel_calc_normalization_range(itrack, gtracks, gvertices, kmin, kmax);
381  }
382  }
383 
384  // (un-)apply the factor -beta which is needed in exp_arg, but not in swE
385  if (updateTc) {
386  auto obeta = -1. / beta;
387  for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
388  gvertices.swE[ivertex] *= obeta;
389  }
390  }
391 
392  // now update z and rho
393  auto kernel_calc_z = [osumtkwt, nv](vertex_t& vertices) -> double {
394  double delta = 0;
395  // does not vectorize
396  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
397  if (vertices.sw[ivertex] > 0) {
398  auto znew = vertices.swz[ivertex] / vertices.sw[ivertex];
399  // prevents from vectorizing if
400  delta = max(std::abs(vertices.zvtx[ivertex] - znew), delta);
401  vertices.zvtx[ivertex] = znew;
402  }
403  }
404 
405  for (unsigned int ivertex = 0; ivertex < nv; ++ivertex)
406  vertices.rho[ivertex] = vertices.rho[ivertex] * vertices.se[ivertex] * osumtkwt;
407 
408  return delta;
409  };
410 
411  double delta = kernel_calc_z(gvertices);
412 
413  // return how much the prototypes moved
414  return delta;
415 }
T w() const
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
std::vector< TransientVertex > vertices(const std::vector< reco::TransientTrack > &tracks) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int nt
Definition: AMPTWrapper.h:42

◆ verify()

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

Definition at line 111 of file DAClusterizerInZ_vect.cc.

References cms::cuda::assert(), gather_cfg::cout, DAClusterizerInZ_vect::track_t::dz2, DAClusterizerInZ_vect::track_t::dz2_vec, DAClusterizerInZ_vect::track_t::getSize(), mps_fire::i, dqmdumpme::k, DAClusterizerInZ_vect::track_t::kmax, DAClusterizerInZ_vect::track_t::kmin, nt, DAClusterizerInZ_vect::track_t::sum_Z, DAClusterizerInZ_vect::track_t::sum_Z_vec, DAClusterizerInZ_vect::track_t::tkwt, DAClusterizerInZ_vect::track_t::tkwt_vec, DAClusterizerInZ_vect::track_t::tt, findQualityFiles::v, DAClusterizerInZ_vect::track_t::zpca, and DAClusterizerInZ_vect::track_t::zpca_vec.

111  {
112  if (!(nv == 999999)) {
113  assert(nv == v.getSize());
114  } else {
115  nv = v.getSize();
116  }
117 
118  if (!(nt == 999999)) {
119  assert(nt == tks.getSize());
120  } else {
121  nt = tks.getSize();
122  }
123 
124  assert(v.zvtx_vec.size() == nv);
125  assert(v.rho_vec.size() == nv);
126  assert(v.swz_vec.size() == nv);
127  assert(v.exp_arg_vec.size() == nv);
128  assert(v.exp_vec.size() == nv);
129  assert(v.se_vec.size() == nv);
130  assert(v.swz_vec.size() == nv);
131  assert(v.swE_vec.size() == nv);
132 
133  assert(v.zvtx == &v.zvtx_vec.front());
134  assert(v.rho == &v.rho_vec.front());
135  assert(v.exp_arg == &v.exp_arg_vec.front());
136  assert(v.sw == &v.sw_vec.front());
137  assert(v.swz == &v.swz_vec.front());
138  assert(v.se == &v.se_vec.front());
139  assert(v.swE == &v.swE_vec.front());
140 
141  for (unsigned int k = 0; k < nv - 1; k++) {
142  if (v.zvtx_vec[k] <= v.zvtx_vec[k + 1])
143  continue;
144  cout << " Z, cluster z-ordering assertion failure z[" << k << "] =" << v.zvtx_vec[k] << " z[" << k + 1
145  << "] =" << v.zvtx_vec[k + 1] << endl;
146  }
147 
148  assert(nt == tks.zpca_vec.size());
149  assert(nt == tks.dz2_vec.size());
150  assert(nt == tks.tt.size());
151  assert(nt == tks.tkwt_vec.size());
152  assert(nt == tks.sum_Z_vec.size());
153  assert(nt == tks.kmin.size());
154  assert(nt == tks.kmax.size());
155 
156  assert(tks.zpca == &tks.zpca_vec.front());
157  assert(tks.dz2 == &tks.dz2_vec.front());
158  assert(tks.tkwt == &tks.tkwt_vec.front());
159  assert(tks.sum_Z == &tks.sum_Z_vec.front());
160 
161  for (unsigned int i = 0; i < nt; i++) {
162  if ((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv))
163  continue;
164  cout << "track vertex range assertion failure" << i << "/" << nt << " kmin,kmax=" << tks.kmin[i] << ", "
165  << tks.kmax[i] << " nv=" << nv << endl;
166  }
167 
168  for (unsigned int i = 0; i < nt; i++) {
169  assert((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv));
170  }
171 }
assert(be >=bs)
int nt
Definition: AMPTWrapper.h:42

◆ vertices()

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

Definition at line 766 of file DAClusterizerInZ_vect.cc.

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

766  {
767  track_t&& tks = fill(tracks);
768  tks.extractRaw();
769 
770  unsigned int nt = tks.getSize();
771  double rho0 = 0.0; // start with no outlier rejection
772 
773  vector<TransientVertex> clusters;
774  if (tks.getSize() == 0)
775  return clusters;
776 
777  vertex_t y; // the vertex prototypes
778 
779  // initialize:single vertex at infinite temperature
780  y.addItem(0, 1.0);
781  clear_vtx_range(tks, y);
782 
783  // estimate first critical temperature
784  double beta = beta0(betamax_, tks, y);
785 #ifdef DEBUG
786  if (DEBUGLEVEL > 0)
787  std::cout << "Beta0 is " << beta << std::endl;
788 #endif
789 
790  thermalize(beta, tks, y, delta_highT_);
791 
792  // annealing loop, stop when T<Tmin (i.e. beta>1/Tmin)
793 
794  double betafreeze = betamax_ * sqrt(coolingFactor_);
795 
796  while (beta < betafreeze) {
797  while (merge(y, tks, beta)) {
798  update(beta, tks, y, rho0, false);
799  }
800  split(beta, tks, y);
801 
803  thermalize(beta, tks, y, delta_highT_);
804  }
805 
806 #ifdef DEBUG
807  verify(y, tks);
808 
809  if (DEBUGLEVEL > 0) {
810  std::cout << "DAClusterizerInZ_vect::vertices :"
811  << "last round of splitting" << std::endl;
812  }
813 #endif
814 
815  set_vtx_range(beta, tks, y);
816  update(beta, tks, y, rho0, false);
817 
818  while (merge(y, tks, beta)) {
819  set_vtx_range(beta, tks, y);
820  update(beta, tks, y, rho0, false);
821  }
822 
823  unsigned int ntry = 0;
824  double threshold = 1.0;
825  while (split(beta, tks, y, threshold) && (ntry++ < 10)) {
826  thermalize(beta, tks, y, delta_highT_, rho0); // rho0 = 0. here
827  while (merge(y, tks, beta)) {
828  update(beta, tks, y, rho0, false);
829  }
830 
831  // relax splitting a bit to reduce multiple split-merge cycles of the same cluster
832  threshold *= 1.1;
833  }
834 
835 #ifdef DEBUG
836  verify(y, tks);
837  if (DEBUGLEVEL > 0) {
838  std::cout << "DAClusterizerInZ_vect::vertices :"
839  << "turning on outlier rejection at T=" << 1 / beta << std::endl;
840  }
841 #endif
842 
843  // switch on outlier rejection at T=Tmin
844  if (dzCutOff_ > 0) {
845  rho0 = y.getSize() > 1 ? 1. / y.getSize() : 1.;
846  for (unsigned int a = 0; a < 5; a++) {
847  update(beta, tks, y, a * rho0 / 5.); // adiabatic turn-on
848  }
849  }
850 
851  thermalize(beta, tks, y, delta_lowT_, rho0);
852 
853 #ifdef DEBUG
854  verify(y, tks);
855  if (DEBUGLEVEL > 0) {
856  std::cout << "DAClusterizerInZ_vect::vertices :"
857  << "merging with outlier rejection at T=" << 1 / beta << std::endl;
858  }
859  if (DEBUGLEVEL > 2)
860  dump(beta, y, tks, 2, rho0);
861 #endif
862 
863  // merge again (some cluster split by outliers collapse here)
864  while (merge(y, tks, beta)) {
865  set_vtx_range(beta, tks, y);
866  update(beta, tks, y, rho0, false);
867  }
868 
869 #ifdef DEBUG
870  verify(y, tks);
871  if (DEBUGLEVEL > 0) {
872  std::cout << "DAClusterizerInZ_vect::vertices :"
873  << "after merging with outlier rejection at T=" << 1 / beta << std::endl;
874  }
875  if (DEBUGLEVEL > 2)
876  dump(beta, y, tks, 2, rho0);
877 #endif
878 
879  // go down to the purging temperature (if it is lower than tmin)
880  while (beta < betapurge_) {
882  thermalize(beta, tks, y, delta_lowT_, rho0);
883  }
884 
885 #ifdef DEBUG
886  verify(y, tks);
887  if (DEBUGLEVEL > 0) {
888  std::cout << "DAClusterizerInZ_vect::vertices :"
889  << "purging at T=" << 1 / beta << std::endl;
890  }
891 #endif
892 
893  // eliminate insignificant vertices, this is more restrictive at higher T
894  while (purge(y, tks, rho0, beta)) {
895  thermalize(beta, tks, y, delta_lowT_, rho0);
896  }
897 
898 #ifdef DEBUG
899  verify(y, tks);
900  if (DEBUGLEVEL > 0) {
901  std::cout << "DAClusterizerInZ_vect::vertices :"
902  << "last cooling T=" << 1 / beta << std::endl;
903  }
904 #endif
905 
906  // optionally cool some more without doing anything, to make the assignment harder
907  while (beta < betastop_) {
909  thermalize(beta, tks, y, delta_lowT_, rho0);
910  }
911 
912 #ifdef DEBUG
913  verify(y, tks);
914  if (DEBUGLEVEL > 0) {
915  std::cout << "DAClusterizerInZ_vect::vertices :"
916  << "stop cooling at T=" << 1 / beta << std::endl;
917  }
918  if (DEBUGLEVEL > 2)
919  dump(beta, y, tks, 2, rho0);
920 #endif
921 
922  // select significant tracks and use a TransientVertex as a container
923 
924  set_vtx_range(beta, tks, y);
925  const unsigned int nv = y.getSize();
926  for (unsigned int k = 0; k < nv; k++) {
927  if (edm::isNotFinite(y.rho[k]) || edm::isNotFinite(y.zvtx[k])) {
928  y.rho[k] = 0;
929  y.zvtx[k] = 0;
930  }
931  }
932 
933  const auto z_sum_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
934  std::vector<std::vector<unsigned int>> vtx_track_indices(nv);
935  for (unsigned int i = 0; i < nt; i++) {
936  const auto kmin = tks.kmin[i];
937  const auto kmax = tks.kmax[i];
938  for (auto k = kmin; k < kmax; k++) {
939  y.exp_arg[k] = -beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i]);
940  }
941 
942  local_exp_list_range(y.exp_arg, y.exp, kmin, kmax);
943 
944  tks.sum_Z[i] = z_sum_init;
945  for (auto k = kmin; k < kmax; k++) {
946  tks.sum_Z[i] += y.rho[k] * y.exp[k];
947  }
948  const double invZ = tks.sum_Z[i] > 1e-100 ? 1. / tks.sum_Z[i] : 0.0;
949 
950  for (auto k = kmin; k < kmax; k++) {
951  double p = y.rho[k] * y.exp[k] * invZ;
952  if (p > mintrkweight_) {
953  // assign track i -> vertex k (hard, mintrkweight should be >= 0.5 here
954  vtx_track_indices[k].push_back(i);
955  break;
956  }
957  }
958 
959  } // track loop
960 
961  GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01);
962  for (unsigned int k = 0; k < nv; k++) {
963  if (!vtx_track_indices[k].empty()) {
964  GlobalPoint pos(0, 0, y.zvtx[k]);
965  vector<reco::TransientTrack> vertexTracks;
966  for (auto i : vtx_track_indices[k]) {
967  vertexTracks.push_back(*(tks.tt[i]));
968  }
969  TransientVertex v(pos, dummyError, vertexTracks, 0);
970  clusters.push_back(v);
971  }
972  }
973 
974  return clusters;
975 }
double beta0(const double betamax, track_t const &tks, vertex_t const &y) const
bool merge(vertex_t &y, track_t &tks, double &beta) const
void clear_vtx_range(track_t &gtracks, vertex_t &gvertices) const
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
T sqrt(T t)
Definition: SSEVec.h:19
void dump(const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0, const double rho0=0.) const
track_t fill(const std::vector< reco::TransientTrack > &tracks) const
int nt
Definition: AMPTWrapper.h:42
unsigned int thermalize(double beta, track_t &gtracks, vertex_t &gvertices, const double delta_max, const double rho0=0.) const
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const
double a
Definition: hdecay.h:119
void verify(const vertex_t &v, const track_t &tks, unsigned int nv=999999, unsigned int nt=999999) const
bool split(const double beta, track_t &t, vertex_t &y, double threshold=1.) const

◆ vertices_in_blocks()

vector< TransientVertex > DAClusterizerInZ_vect::vertices_in_blocks ( const std::vector< reco::TransientTrack > &  tracks) const

Definition at line 977 of file DAClusterizerInZ_vect.cc.

References a, b, HLT_2023v12_cff::beta, groupFilesInBlocks::block, bsc_activity_cfg::clusters, gather_cfg::cout, GCP_Ntuples_cfg::dump, MillePedeFileConverter_cfg::e, relativeConstraints::empty, DAClusterizerInZ_vect::track_t::extractRaw(), ntuplemaker::fill, dqmdumpme::first, mps_fire::i, createfilelist::int, dqmdumpme::k, SiStripPI::max, merge(), SiStripPI::min, nt, AlCaHLTBitMon_ParallelJobs::p, position, push_back(), jetUpdater_cfi::sort, submitPVValidationJobs::split(), mathSSE::sqrt(), remoteMonitoring_LASER_era2018_cfg::threshold, pwdgSkimBPark_cfi::tracks, mitigatedMETSequence_cff::U, update, findQualityFiles::v, SiStripMonitorCluster_cfi::zmax, SiStripMonitorCluster_cfi::zmin, and TkClusParameters_cff::zrange.

977  {
978  vector<reco::TransientTrack> sorted_tracks;
979  vector<pair<float, float>> vertices_tot; // z, rho for each vertex
980  for (unsigned int i = 0; i < tracks.size(); i++) {
981  sorted_tracks.push_back(tracks[i]);
982  }
983  double rho0, beta;
984  std::sort(sorted_tracks.begin(),
985  sorted_tracks.end(),
986  [](const reco::TransientTrack& a, const reco::TransientTrack& b) -> bool {
987  return (a.stateAtBeamLine().trackStateAtPCA()).position().z() <
988  (b.stateAtBeamLine().trackStateAtPCA()).position().z();
989  });
990 
991  unsigned int nBlocks = (unsigned int)std::floor(sorted_tracks.size() / (block_size_ * (1 - overlap_frac_)));
992  if (nBlocks < 1) {
993  nBlocks = 1;
994  edm::LogWarning("DAClusterizerinZ_vect")
995  << "Warning nBlocks was 0 with ntracks = " << sorted_tracks.size() << " block_size = " << block_size_
996  << " and overlap fraction = " << overlap_frac_ << ". Setting nBlocks = 1";
997  }
998  for (unsigned int block = 0; block < nBlocks; block++) {
999  vector<reco::TransientTrack> block_tracks;
1000  unsigned int begin = (unsigned int)(block * block_size_ * (1 - overlap_frac_));
1001  unsigned int end = (unsigned int)std::min(begin + block_size_, (unsigned int)sorted_tracks.size());
1002  for (unsigned int i = begin; i < end; i++) {
1003  block_tracks.push_back(sorted_tracks[i]);
1004  }
1005  if (block_tracks.empty()) {
1006  continue;
1007  }
1008 
1009 #ifdef DEBUG
1010  std::cout << "Running vertices_in_blocks on" << std::endl;
1011  std::cout << "- block no." << block << " on " << nBlocks << " blocks " << std::endl;
1012  std::cout << "- block track size: " << sorted_tracks.size() << " - block size: " << block_size_ << std::endl;
1013 #endif
1014  track_t&& tks = fill(block_tracks);
1015  tks.extractRaw();
1016 
1017  rho0 = 0.0; // start with no outlier rejection
1018 
1019  vertex_t y; // the vertex prototypes
1020 
1021  // initialize:single vertex at infinite temperature
1022  y.addItem(0, 1.0);
1023  clear_vtx_range(tks, y);
1024 
1025  // estimate first critical temperature
1026  beta = beta0(betamax_, tks, y);
1027 #ifdef DEBUG
1028  if (DEBUGLEVEL > 0)
1029  std::cout << "Beta0 is " << beta << std::endl;
1030 #endif
1031 
1032  thermalize(beta, tks, y, delta_highT_);
1033 
1034  // annealing loop, stop when T<Tmin (i.e. beta>1/Tmin)
1035 
1036  double betafreeze = betamax_ * sqrt(coolingFactor_);
1037  while (beta < betafreeze) {
1038  while (merge(y, tks, beta)) {
1039  update(beta, tks, y, rho0, false);
1040  }
1041  split(beta, tks, y);
1042 
1043  beta = beta / coolingFactor_;
1044  thermalize(beta, tks, y, delta_highT_);
1045  }
1046 
1047 #ifdef DEBUG
1048  verify(y, tks);
1049 
1050  if (DEBUGLEVEL > 0) {
1051  std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1052  << "last round of splitting" << std::endl;
1053  }
1054 #endif
1055 
1056  set_vtx_range(beta, tks, y);
1057  update(beta, tks, y, rho0, false);
1058 
1059  while (merge(y, tks, beta)) {
1060  set_vtx_range(beta, tks, y);
1061  update(beta, tks, y, rho0, false);
1062  }
1063 
1064  unsigned int ntry = 0;
1065  double threshold = 1.0;
1066  while (split(beta, tks, y, threshold) && (ntry++ < 10)) {
1067  thermalize(beta, tks, y, delta_highT_, rho0); // rho0 = 0. here
1068  while (merge(y, tks, beta)) {
1069  update(beta, tks, y, rho0, false);
1070  }
1071 
1072  // relax splitting a bit to reduce multiple split-merge cycles of the same cluster
1073  threshold *= 1.1;
1074  }
1075 
1076 #ifdef DEBUG
1077  verify(y, tks);
1078  if (DEBUGLEVEL > 0) {
1079  std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1080  << "turning on outlier rejection at T=" << 1 / beta << std::endl;
1081  }
1082 #endif
1083 
1084  // switch on outlier rejection at T=Tmin
1085  if (dzCutOff_ > 0) {
1086  rho0 = y.getSize() > 1 ? 1. / y.getSize() : 1.;
1087  for (unsigned int a = 0; a < 5; a++) {
1088  update(beta, tks, y, a * rho0 / 5.); // adiabatic turn-on
1089  }
1090  }
1091 
1092  thermalize(beta, tks, y, delta_lowT_, rho0);
1093 
1094 #ifdef DEBUG
1095  verify(y, tks);
1096  if (DEBUGLEVEL > 0) {
1097  std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1098  << "merging with outlier rejection at T=" << 1 / beta << std::endl;
1099  }
1100  if (DEBUGLEVEL > 2)
1101  dump(beta, y, tks, 2, rho0);
1102 #endif
1103 
1104  // merge again (some cluster split by outliers collapse here)
1105  while (merge(y, tks, beta)) {
1106  set_vtx_range(beta, tks, y);
1107  update(beta, tks, y, rho0, false);
1108  }
1109 
1110 #ifdef DEBUG
1111  verify(y, tks);
1112  if (DEBUGLEVEL > 0) {
1113  std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1114  << "after merging with outlier rejection at T=" << 1 / beta << std::endl;
1115  }
1116  if (DEBUGLEVEL > 2)
1117  dump(beta, y, tks, 2, rho0);
1118 #endif
1119 
1120  // go down to the purging temperature (if it is lower than tmin)
1121  while (beta < betapurge_) {
1123  thermalize(beta, tks, y, delta_lowT_, rho0);
1124  }
1125 
1126 #ifdef DEBUG
1127  verify(y, tks);
1128  if (DEBUGLEVEL > 0) {
1129  std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1130  << "purging at T=" << 1 / beta << std::endl;
1131  }
1132 #endif
1133 
1134  // eliminate insignificant vertices, this is more restrictive at higher T
1135  while (purge(y, tks, rho0, beta)) {
1136  thermalize(beta, tks, y, delta_lowT_, rho0);
1137  }
1138 
1139 #ifdef DEBUG
1140  verify(y, tks);
1141  if (DEBUGLEVEL > 0) {
1142  std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1143  << "last cooling T=" << 1 / beta << std::endl;
1144  }
1145 #endif
1146 
1147  // optionally cool some more without doing anything, to make the assignment harder
1148  while (beta < betastop_) {
1150  thermalize(beta, tks, y, delta_lowT_, rho0);
1151  }
1152 
1153 #ifdef DEBUG
1154  verify(y, tks);
1155  if (DEBUGLEVEL > 0) {
1156  std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1157  << "stop cooling at T=" << 1 / beta << std::endl;
1158  }
1159  if (DEBUGLEVEL > 2)
1160  dump(beta, y, tks, 2, rho0);
1161 #endif
1162 
1163  for (unsigned int ivertex = 0; ivertex < y.getSize(); ivertex++) {
1164  if (y.zvtx_vec[ivertex] != 0 && y.rho_vec[ivertex] != 0) {
1165  vertices_tot.push_back(pair(y.zvtx_vec[ivertex], y.rho_vec[ivertex]));
1166 #ifdef DEBUG
1167  std::cout << "Found new vertex " << y.zvtx_vec[ivertex] << " , " << y.rho_vec[ivertex] << std::endl;
1168 #endif
1169  }
1170  }
1171  }
1172 
1173  std::sort(vertices_tot.begin(),
1174  vertices_tot.end(),
1175  [](const pair<float, float>& a, const pair<float, float>& b) -> bool { return a.first < b.first; });
1176 
1177  // reassign tracks to vertices
1178  track_t&& tracks_tot = fill(tracks);
1179  const unsigned int nv = vertices_tot.size();
1180  const unsigned int nt = tracks_tot.getSize();
1181 
1182  for (auto itrack = 0U; itrack < nt; ++itrack) {
1183  double zrange = max(sel_zrange_ / sqrt(beta * tracks_tot.dz2[itrack]), zrange_min_);
1184 
1185  double zmin = tracks_tot.zpca[itrack] - zrange;
1186  unsigned int kmin = min(nv - 1, tracks_tot.kmin[itrack]);
1187  // find the smallest vertex_z that is larger than zmin
1188  if (vertices_tot[kmin].first > zmin) {
1189  while ((kmin > 0) && (vertices_tot[kmin - 1].first > zmin)) {
1190  kmin--;
1191  }
1192  } else {
1193  while ((kmin < (nv - 1)) && (vertices_tot[kmin].first < zmin)) {
1194  kmin++;
1195  }
1196  }
1197 
1198  double zmax = tracks_tot.zpca[itrack] + zrange;
1199  unsigned int kmax = min(nv - 1, tracks_tot.kmax[itrack] - 1);
1200  // note: kmax points to the last vertex in the range, while gtracks.kmax points to the entry BEHIND the last vertex
1201  // find the largest vertex_z that is smaller than zmax
1202  if (vertices_tot[kmax].first < zmax) {
1203  while ((kmax < (nv - 1)) && (vertices_tot[kmax + 1].first < zmax)) {
1204  kmax++;
1205  }
1206  } else {
1207  while ((kmax > 0) && (vertices_tot[kmax].first > zmax)) {
1208  kmax--;
1209  }
1210  }
1211 
1212  if (kmin <= kmax) {
1213  tracks_tot.kmin[itrack] = kmin;
1214  tracks_tot.kmax[itrack] = kmax + 1;
1215  } else {
1216  tracks_tot.kmin[itrack] = max(0U, min(kmin, kmax));
1217  tracks_tot.kmax[itrack] = min(nv, max(kmin, kmax) + 1);
1218  }
1219  }
1220 
1221  rho0 = nv > 1 ? 1. / nv : 1.;
1222  const auto z_sum_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
1223 
1224  std::vector<std::vector<unsigned int>> vtx_track_indices(nv);
1225  for (unsigned int i = 0; i < nt; i++) {
1226  const auto kmin = tracks_tot.kmin[i];
1227  const auto kmax = tracks_tot.kmax[i];
1228  double p_max = -1;
1229  unsigned int iMax = 10000; //just a "big" number w.r.t. number of vertices
1230  float sum_Z = z_sum_init;
1231  for (auto k = kmin; k < kmax; k++) {
1232  float v_exp = local_exp(-beta * Eik(tracks_tot.zpca[i], vertices_tot[k].first, tracks_tot.dz2[i]));
1233  sum_Z += vertices_tot[k].second * v_exp;
1234  }
1235  double invZ = sum_Z > 1e-100 ? 1. / sum_Z : 0.0;
1236  for (auto k = kmin; k < kmax && invZ != 0.0; k++) {
1237  float v_exp = local_exp(-beta * Eik(tracks_tot.zpca[i], vertices_tot[k].first, tracks_tot.dz2[i]));
1238  double p = vertices_tot[k].second * v_exp * invZ;
1239  if (p > p_max && p > mintrkweight_) {
1240  p_max = p;
1241  iMax = k;
1242  }
1243  }
1244  if (iMax < vtx_track_indices.size())
1245  vtx_track_indices[iMax].push_back(i);
1246  }
1247 #ifdef DEBUG
1248  for (auto itrack = 0U; itrack < nt; ++itrack) {
1249  std::cout << "itrack " << itrack << " , " << tracks_tot.kmin[itrack] << " , " << tracks_tot.kmax[itrack]
1250  << std::endl;
1251  }
1252 #endif
1253 
1254  vector<TransientVertex> clusters;
1255 
1256  GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01);
1257  for (unsigned int k = 0; k < nv; k++) {
1258  if (!vtx_track_indices[k].empty()) {
1259  GlobalPoint pos(0, 0, vertices_tot[k].first);
1260  vector<reco::TransientTrack> vertexTracks;
1261  for (auto i : vtx_track_indices[k]) {
1262  vertexTracks.push_back(*(tracks_tot.tt[i]));
1263 #ifdef DEBUG
1264  std::cout << y.zvtx[k] << "," << (*tks.tt[i]).stateAtBeamLine().trackStateAtPCA().position().z() << std::endl;
1265 #endif
1266  }
1267  TransientVertex v(pos, dummyError, vertexTracks, 0);
1268  clusters.push_back(v);
1269  }
1270  }
1271 
1272  return clusters;
1273 }
double beta0(const double betamax, track_t const &tks, vertex_t const &y) const
bool merge(vertex_t &y, track_t &tks, double &beta) const
void clear_vtx_range(track_t &gtracks, vertex_t &gvertices) const
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
T sqrt(T t)
Definition: SSEVec.h:19
void dump(const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0, const double rho0=0.) const
track_t fill(const std::vector< reco::TransientTrack > &tracks) const
int nt
Definition: AMPTWrapper.h:42
unsigned int thermalize(double beta, track_t &gtracks, vertex_t &gvertices, const double delta_max, const double rho0=0.) const
void set_vtx_range(double beta, track_t &gtracks, vertex_t &gvertices) const
double b
Definition: hdecay.h:118
deadvectors [0] push_back({0.0175431, 0.538005, 6.80997, 13.29})
double a
Definition: hdecay.h:119
static int position[264][3]
Definition: ReadPGInfo.cc:289
Log< level::Warning, false > LogWarning
void verify(const vertex_t &v, const track_t &tks, unsigned int nv=999999, unsigned int nt=999999) const
bool split(const double beta, track_t &t, vertex_t &y, double threshold=1.) const

Member Data Documentation

◆ betamax_

double DAClusterizerInZ_vect::betamax_
private

Definition at line 209 of file DAClusterizerInZ_vect.h.

◆ betapurge_

double DAClusterizerInZ_vect::betapurge_
private

Definition at line 218 of file DAClusterizerInZ_vect.h.

◆ betastop_

double DAClusterizerInZ_vect::betastop_
private

Definition at line 210 of file DAClusterizerInZ_vect.h.

◆ block_size_

unsigned int DAClusterizerInZ_vect::block_size_
private

Definition at line 228 of file DAClusterizerInZ_vect.h.

◆ convergence_mode_

unsigned int DAClusterizerInZ_vect::convergence_mode_
private

Definition at line 220 of file DAClusterizerInZ_vect.h.

◆ coolingFactor_

double DAClusterizerInZ_vect::coolingFactor_
private

Definition at line 208 of file DAClusterizerInZ_vect.h.

◆ d0CutOff_

double DAClusterizerInZ_vect::d0CutOff_
private

Definition at line 212 of file DAClusterizerInZ_vect.h.

◆ delta_highT_

double DAClusterizerInZ_vect::delta_highT_
private

Definition at line 221 of file DAClusterizerInZ_vect.h.

◆ delta_lowT_

double DAClusterizerInZ_vect::delta_lowT_
private

Definition at line 222 of file DAClusterizerInZ_vect.h.

◆ dzCutOff_

double DAClusterizerInZ_vect::dzCutOff_
private

Definition at line 211 of file DAClusterizerInZ_vect.h.

◆ maxIterations_

unsigned int DAClusterizerInZ_vect::maxIterations_
private

Definition at line 207 of file DAClusterizerInZ_vect.h.

◆ mintrkweight_

double DAClusterizerInZ_vect::mintrkweight_
private

Definition at line 214 of file DAClusterizerInZ_vect.h.

◆ overlap_frac_

double DAClusterizerInZ_vect::overlap_frac_
private

Definition at line 229 of file DAClusterizerInZ_vect.h.

◆ runInBlocks_

bool DAClusterizerInZ_vect::runInBlocks_
private

Definition at line 227 of file DAClusterizerInZ_vect.h.

◆ sel_zrange_

double DAClusterizerInZ_vect::sel_zrange_
private

Definition at line 224 of file DAClusterizerInZ_vect.h.

◆ uniquetrkminp_

double DAClusterizerInZ_vect::uniquetrkminp_
private

Definition at line 216 of file DAClusterizerInZ_vect.h.

◆ uniquetrkweight_

double DAClusterizerInZ_vect::uniquetrkweight_
private

Definition at line 215 of file DAClusterizerInZ_vect.h.

◆ vertexSize_

double DAClusterizerInZ_vect::vertexSize_
private

Definition at line 206 of file DAClusterizerInZ_vect.h.

◆ zdumpcenter_

double DAClusterizerInZ_vect::zdumpcenter_
private

Definition at line 203 of file DAClusterizerInZ_vect.h.

◆ zdumpwidth_

double DAClusterizerInZ_vect::zdumpwidth_
private

Definition at line 204 of file DAClusterizerInZ_vect.h.

◆ zmerge_

double DAClusterizerInZ_vect::zmerge_
private

Definition at line 217 of file DAClusterizerInZ_vect.h.

◆ zrange_min_

const double DAClusterizerInZ_vect::zrange_min_ = 0.1
private

Definition at line 225 of file DAClusterizerInZ_vect.h.