PixelSeeding/src/CellularAutomaton.cc

0001 #include <queue>
0002
0003 #include "CellularAutomaton.h"
0004
0005 void CellularAutomaton::createAndConnectCells(const std::vector<const HitDoublets *> &hitDoublets,
0006                                               const TrackingRegion &region,
0007                                               const CACut &thetaCut,
0008                                               const CACut &phiCut,
0009                                               const float hardPtCut) {
0010   int tsize = 0;
0011   for (auto hd : hitDoublets) {
0012     tsize += hd->size();
0013   }
0014   allCells.reserve(tsize);
0015   unsigned int cellId = 0;
0016   float ptmin = region.ptMin();
0017   float region_origin_x = region.origin().x();
0018   float region_origin_y = region.origin().y();
0019   float region_origin_radius = region.originRBound();
0020
0021   std::vector<bool> alreadyVisitedLayerPairs;
0022   alreadyVisitedLayerPairs.resize(theLayerGraph.theLayerPairs.size());
0023   for (auto visited : alreadyVisitedLayerPairs) {
0024     visited = false;
0025   }
0026   for (int rootVertex : theLayerGraph.theRootLayers) {
0027     std::queue<int> LayerPairsToVisit;
0028
0029     for (int LayerPair : theLayerGraph.theLayers[rootVertex].theOuterLayerPairs) {
0030       LayerPairsToVisit.push(LayerPair);
0031     }
0032
0033     unsigned int numberOfLayerPairsToVisitAtThisDepth = LayerPairsToVisit.size();
0034
0035     while (not LayerPairsToVisit.empty()) {
0036       auto currentLayerPair = LayerPairsToVisit.front();
0037       auto &currentLayerPairRef = theLayerGraph.theLayerPairs[currentLayerPair];
0038       auto &currentInnerLayerRef = theLayerGraph.theLayers[currentLayerPairRef.theLayers[0]];
0039       auto &currentOuterLayerRef = theLayerGraph.theLayers[currentLayerPairRef.theLayers[1]];
0040       bool allInnerLayerPairsAlreadyVisited{true};
0041
0042       CACut::CAValuesByInnerLayerIds caThetaCut =
0043           thetaCut.getCutsByInnerLayer(currentInnerLayerRef.seqNum(), currentOuterLayerRef.seqNum());
0044       CACut::CAValuesByInnerLayerIds caPhiCut =
0045           phiCut.getCutsByInnerLayer(currentInnerLayerRef.seqNum(), currentOuterLayerRef.seqNum());
0046
0047       for (auto innerLayerPair : currentInnerLayerRef.theInnerLayerPairs) {
0048         allInnerLayerPairsAlreadyVisited &= alreadyVisitedLayerPairs[innerLayerPair];
0049       }
0050
0051       if (alreadyVisitedLayerPairs[currentLayerPair] == false and allInnerLayerPairsAlreadyVisited) {
0052         const HitDoublets *doubletLayerPairId = hitDoublets[currentLayerPair];
0053         auto numberOfDoublets = doubletLayerPairId->size();
0054         currentLayerPairRef.theFoundCells[0] = cellId;
0055         currentLayerPairRef.theFoundCells[1] = cellId + numberOfDoublets;
0056         for (unsigned int i = 0; i < numberOfDoublets; ++i) {
0057           allCells.emplace_back(
0058               doubletLayerPairId, i, doubletLayerPairId->innerHitId(i), doubletLayerPairId->outerHitId(i));
0059
0060           currentOuterLayerRef.isOuterHitOfCell[doubletLayerPairId->outerHitId(i)].push_back(cellId);
0061
0062           cellId++;
0063
0064           auto &neigCells = currentInnerLayerRef.isOuterHitOfCell[doubletLayerPairId->innerHitId(i)];
0065           allCells.back().checkAlignmentAndTag(allCells,
0066                                                neigCells,
0067                                                ptmin,
0068                                                region_origin_x,
0069                                                region_origin_y,
0070                                                region_origin_radius,
0071                                                caThetaCut,
0072                                                caPhiCut,
0073                                                hardPtCut);
0074         }
0075         assert(cellId == currentLayerPairRef.theFoundCells[1]);
0076         for (auto outerLayerPair : currentOuterLayerRef.theOuterLayerPairs) {
0077           LayerPairsToVisit.push(outerLayerPair);
0078         }
0079
0080         alreadyVisitedLayerPairs[currentLayerPair] = true;
0081       }
0082       LayerPairsToVisit.pop();
0083       numberOfLayerPairsToVisitAtThisDepth--;
0084       if (numberOfLayerPairsToVisitAtThisDepth == 0) {
0085         numberOfLayerPairsToVisitAtThisDepth = LayerPairsToVisit.size();
0086       }
0087     }
0088   }
0089 }
0090
0091 void CellularAutomaton::evolve(const unsigned int minHitsPerNtuplet) {
0092   allStatus.resize(allCells.size());
0093
0094   unsigned int numberOfIterations = minHitsPerNtuplet - 2;
0095   // keeping the last iteration for later
0096   for (unsigned int iteration = 0; iteration < numberOfIterations - 1; ++iteration) {
0097     for (auto &layerPair : theLayerGraph.theLayerPairs) {
0098       for (auto i = layerPair.theFoundCells[0]; i < layerPair.theFoundCells[1]; ++i) {
0099         allCells[i].evolve(i, allStatus);
0100       }
0101     }
0102
0103     for (auto &layerPair : theLayerGraph.theLayerPairs) {
0104       for (auto i = layerPair.theFoundCells[0]; i < layerPair.theFoundCells[1]; ++i) {
0105         allStatus[i].updateState();
0106       }
0107     }
0108   }
0109
0110   // last iteration
0111
0112   for (int rootLayerId : theLayerGraph.theRootLayers) {
0113     for (int rootLayerPair : theLayerGraph.theLayers[rootLayerId].theOuterLayerPairs) {
0114       auto foundCells = theLayerGraph.theLayerPairs[rootLayerPair].theFoundCells;
0115       for (auto i = foundCells[0]; i < foundCells[1]; ++i) {
0116         auto &cell = allStatus[i];
0117         allCells[i].evolve(i, allStatus);
0118         cell.updateState();
0119         if (cell.isRootCell(minHitsPerNtuplet - 2)) {
0120           theRootCells.push_back(i);
0121         }
0122       }
0123     }
0124   }
0125 }
0126
0127 void CellularAutomaton::findNtuplets(std::vector<CACell::CAntuplet> &foundNtuplets,
0128                                      const unsigned int minHitsPerNtuplet) {
0129   CACell::CAntuple tmpNtuplet;
0130   tmpNtuplet.reserve(minHitsPerNtuplet);
0131
0132   for (auto root_cell : theRootCells) {
0133     tmpNtuplet.clear();
0134     tmpNtuplet.push_back(root_cell);
0135     allCells[root_cell].findNtuplets(allCells, foundNtuplets, tmpNtuplet, minHitsPerNtuplet);
0136   }
0137 }
0138
0139 void CellularAutomaton::findTriplets(std::vector<const HitDoublets *> const &hitDoublets,
0140                                      std::vector<CACell::CAntuplet> &foundTriplets,
0141                                      TrackingRegion const &region,
0142                                      const CACut &thetaCut,
0143                                      const CACut &phiCut,
0144                                      const float hardPtCut) {
0145   int tsize = 0;
0146   for (auto hd : hitDoublets) {
0147     tsize += hd->size();
0148   }
0149   allCells.reserve(tsize);
0150
0151   unsigned int cellId = 0;
0152   float ptmin = region.ptMin();
0153   float region_origin_x = region.origin().x();
0154   float region_origin_y = region.origin().y();
0155   float region_origin_radius = region.originRBound();
0156
0157   std::vector<bool> alreadyVisitedLayerPairs;
0158   alreadyVisitedLayerPairs.resize(theLayerGraph.theLayerPairs.size());
0159   for (auto visited : alreadyVisitedLayerPairs) {
0160     visited = false;
0161   }
0162   for (int rootVertex : theLayerGraph.theRootLayers) {
0163     std::queue<int> LayerPairsToVisit;
0164
0165     for (int LayerPair : theLayerGraph.theLayers[rootVertex].theOuterLayerPairs) {
0166       LayerPairsToVisit.push(LayerPair);
0167     }
0168
0169     unsigned int numberOfLayerPairsToVisitAtThisDepth = LayerPairsToVisit.size();
0170
0171     while (not LayerPairsToVisit.empty()) {
0172       auto currentLayerPair = LayerPairsToVisit.front();
0173       auto &currentLayerPairRef = theLayerGraph.theLayerPairs[currentLayerPair];
0174       auto &currentInnerLayerRef = theLayerGraph.theLayers[currentLayerPairRef.theLayers[0]];
0175       auto &currentOuterLayerRef = theLayerGraph.theLayers[currentLayerPairRef.theLayers[1]];
0176       bool allInnerLayerPairsAlreadyVisited{true};
0177
0178       CACut::CAValuesByInnerLayerIds caThetaCut =
0179           thetaCut.getCutsByInnerLayer(currentInnerLayerRef.seqNum(), currentOuterLayerRef.seqNum());
0180       CACut::CAValuesByInnerLayerIds caPhiCut =
0181           phiCut.getCutsByInnerLayer(currentInnerLayerRef.seqNum(), currentOuterLayerRef.seqNum());
0182
0183       for (auto innerLayerPair : currentInnerLayerRef.theInnerLayerPairs) {
0184         allInnerLayerPairsAlreadyVisited &= alreadyVisitedLayerPairs[innerLayerPair];
0185       }
0186
0187       if (alreadyVisitedLayerPairs[currentLayerPair] == false and allInnerLayerPairsAlreadyVisited) {
0188         const HitDoublets *doubletLayerPairId = hitDoublets[currentLayerPair];
0189         auto numberOfDoublets = doubletLayerPairId->size();
0190         currentLayerPairRef.theFoundCells[0] = cellId;
0191         currentLayerPairRef.theFoundCells[1] = cellId + numberOfDoublets;
0192         for (unsigned int i = 0; i < numberOfDoublets; ++i) {
0193           allCells.emplace_back(
0194               doubletLayerPairId, i, doubletLayerPairId->innerHitId(i), doubletLayerPairId->outerHitId(i));
0195
0196           currentOuterLayerRef.isOuterHitOfCell[doubletLayerPairId->outerHitId(i)].push_back(cellId);
0197
0198           cellId++;
0199
0200           auto &neigCells = currentInnerLayerRef.isOuterHitOfCell[doubletLayerPairId->innerHitId(i)];
0201           allCells.back().checkAlignmentAndPushTriplet(allCells,
0202                                                        neigCells,
0203                                                        foundTriplets,
0204                                                        ptmin,
0205                                                        region_origin_x,
0206                                                        region_origin_y,
0207                                                        region_origin_radius,
0208                                                        caThetaCut,
0209                                                        caPhiCut,
0210                                                        hardPtCut);
0211         }
0212         assert(cellId == currentLayerPairRef.theFoundCells[1]);
0213         for (auto outerLayerPair : currentOuterLayerRef.theOuterLayerPairs) {
0214           LayerPairsToVisit.push(outerLayerPair);
0215         }
0216
0217         alreadyVisitedLayerPairs[currentLayerPair] = true;
0218       }
0219       LayerPairsToVisit.pop();
0220       numberOfLayerPairsToVisitAtThisDepth--;
0221       if (numberOfLayerPairsToVisitAtThisDepth == 0) {
0222         numberOfLayerPairsToVisitAtThisDepth = LayerPairsToVisit.size();
0223       }
0224     }
0225   }
0226 }