TICL/plugins/HGCDoublet.cc

0001 #include "HGCDoublet.h"
0002
0003 bool HGCDoublet::checkCompatibilityAndTag(std::vector<HGCDoublet> &allDoublets,
0004                                           const std::vector<int> &innerDoublets,
0005                                           const GlobalVector &refDir,
0006                                           float minCosTheta,
0007                                           float minCosPointing,
0008                                           bool debug) {
0009   int nDoublets = innerDoublets.size();
0010   int constexpr VSIZE = 4;
0011   int ok[VSIZE];
0012   double xi[VSIZE];
0013   double yi[VSIZE];
0014   double zi[VSIZE];
0015   int seedi[VSIZE];
0016   bool siblingsClusters[VSIZE];
0017   auto xo = outerX();
0018   auto yo = outerY();
0019   auto zo = outerZ();
0020   unsigned int doubletId = theDoubletId_;
0021
0022   auto loop = [&](int i, int vs) {
0023     for (int j = 0; j < vs; ++j) {
0024       auto otherDoubletId = innerDoublets[i + j];
0025       auto &otherDoublet = allDoublets[otherDoubletId];
0026       xi[j] = otherDoublet.innerX();
0027       yi[j] = otherDoublet.innerY();
0028       zi[j] = otherDoublet.innerZ();
0029       seedi[j] = otherDoublet.seedIndex();
0030       siblingsClusters[j] = otherDoublet.areSiblingClusters();
0031       if (debug) {
0032         LogDebug("HGCDoublet") << i + j << " is doublet " << otherDoubletId << std::endl;
0033       }
0034     }
0035     for (int j = 0; j < vs; ++j) {
0036       if (seedi[j] != seedIndex_) {
0037         ok[j] = 0;
0038         continue;
0039       }
0040       if (areSiblingClusters_ and siblingsClusters[j]) {
0041         ok[j] = 1;
0042         continue;
0043       }
0044       ok[j] = areAligned(xi[j], yi[j], zi[j], xo, yo, zo, minCosTheta, minCosPointing, refDir, debug);
0045       if (debug) {
0046         LogDebug("HGCDoublet") << "Are aligned for InnerDoubletId: " << i + j << " is " << ok[j] << std::endl;
0047       }
0048     }
0049     for (int j = 0; j < vs; ++j) {
0050       auto otherDoubletId = innerDoublets[i + j];
0051       auto &otherDoublet = allDoublets[otherDoubletId];
0052       if (ok[j]) {
0053         otherDoublet.tagAsOuterNeighbor(doubletId);
0054         allDoublets[doubletId].tagAsInnerNeighbor(otherDoubletId);
0055       }
0056     }
0057   };
0058   auto lim = VSIZE * (nDoublets / VSIZE);
0059   for (int i = 0; i < lim; i += VSIZE)
0060     loop(i, VSIZE);
0061   loop(lim, nDoublets - lim);
0062
0063   if (debug) {
0064     LogDebug("HGCDoublet") << "Found " << innerNeighbors_.size() << " compatible doublets out of " << nDoublets
0065                            << " considered" << std::endl;
0066   }
0067   return innerNeighbors_.empty();
0068 }
0069
0070 int HGCDoublet::areAligned(double xi,
0071                            double yi,
0072                            double zi,
0073                            double xo,
0074                            double yo,
0075                            double zo,
0076                            float minCosTheta,
0077                            float minCosPointing,
0078                            const GlobalVector &refDir,
0079                            bool debug) const {
0080   auto dx1 = xo - xi;
0081   auto dy1 = yo - yi;
0082   auto dz1 = zo - zi;
0083
0084   auto dx2 = innerX() - xi;
0085   auto dy2 = innerY() - yi;
0086   auto dz2 = innerZ() - zi;
0087
0088   // inner product
0089   auto dot = dx1 * dx2 + dy1 * dy2 + dz1 * dz2;
0090   auto dotsq = dot * dot;
0091   // magnitudes
0092   auto mag1sq = dx1 * dx1 + dy1 * dy1 + dz1 * dz1;
0093   auto mag2sq = dx2 * dx2 + dy2 * dy2 + dz2 * dz2;
0094
0095   auto minCosTheta_sq = minCosTheta * minCosTheta;
0096   bool isWithinLimits = (dotsq > minCosTheta_sq * mag1sq * mag2sq);
0097
0098   if (debug) {
0099     LogDebug("HGCDoublet") << "-- Are Aligned -- dotsq: " << dotsq << " mag1sq: " << mag1sq << " mag2sq: " << mag2sq
0100                            << "minCosTheta_sq:" << minCosTheta_sq << " isWithinLimits: " << isWithinLimits << std::endl;
0101   }
0102
0103   // Now check the compatibility with the pointing origin.
0104   // Pointing origin is a vector obtained by the seed (track extrapolation i.e.)
0105   // or the direction wrt (0,0,0)
0106   // The compatibility is checked only for the innermost doublets: the
0107   // one with the outer doublets comes in by the alignment requirement of
0108   // the doublets themeselves
0109
0110   const GlobalVector firstDoublet(dx2, dy2, dz2);
0111   const GlobalVector pointingDir = (seedIndex_ == -1) ? GlobalVector(innerX(), innerY(), innerZ()) : refDir;
0112
0113   auto dot_pointing = pointingDir.dot(firstDoublet);
0114   auto dot_pointing_sq = dot_pointing * dot_pointing;
0115   auto mag_pointing_sq = pointingDir.mag2();
0116   auto minCosPointing_sq = minCosPointing * minCosPointing;
0117   bool isWithinLimitsPointing = (dot_pointing_sq > minCosPointing_sq * mag_pointing_sq * mag2sq);
0118   if (debug) {
0119     LogDebug("HGCDoublet") << "Pointing direction: " << pointingDir << std::endl;
0120     LogDebug("HGCDoublet") << "-- Are Aligned -- dot_pointing_sq: " << dot_pointing_sq
0121                            << " mag_pointing_sq: " << mag_pointing_sq << " mag2sq: " << mag2sq
0122                            << " isWithinLimitsPointing: " << isWithinLimitsPointing << std::endl;
0123   }
0124   // by squaring cosTheta and multiplying by the squares of the magnitudes
0125   // an equivalent comparison is made without the division and square root which are costly FP ops.
0126   return isWithinLimits && isWithinLimitsPointing;
0127 }
0128
0129 void HGCDoublet::findNtuplets(std::vector<HGCDoublet> &allDoublets,
0130                               HGCntuplet &tmpNtuplet,
0131                               int seedIndex,
0132                               const bool outInDFS,
0133                               const unsigned int outInHops,
0134                               const unsigned int maxOutInHops,
0135                               std::vector<std::pair<unsigned int, unsigned int> > &outInToVisit) {
0136   if (!alreadyVisited_ && seedIndex == seedIndex_) {
0137     alreadyVisited_ = true;
0138     tmpNtuplet.push_back(theDoubletId_);
0139     unsigned int numberOfOuterNeighbors = outerNeighbors_.size();
0140     for (unsigned int i = 0; i < numberOfOuterNeighbors; ++i) {
0141       allDoublets[outerNeighbors_[i]].findNtuplets(
0142           allDoublets, tmpNtuplet, seedIndex, outInDFS, outInHops, maxOutInHops, outInToVisit);
0143     }
0144     if (outInDFS && outInHops < maxOutInHops) {
0145       for (auto inN : innerNeighbors_) {
0146         outInToVisit.emplace_back(inN, outInHops + 1);
0147       }
0148     }
0149   }
0150 }