Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:11 

0001 // Author: Felice Pantaleo - felice.pantaleo@cern.ch
0002 // Date: 11/2018
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/HGCalReco/interface/Common.h"
0005 #include "PatternRecognitionbyCA.h"
0006 #include "HGCDoublet.h"
0007 #include "HGCGraph.h"
0008 #include "DataFormats/Common/interface/ValueMap.h"
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 
0011 template <typename TILES>
0012 void HGCGraphT<TILES>::makeAndConnectDoublets(const TILES &histo,
0013                                               const std::vector<TICLSeedingRegion> &regions,
0014                                               int nEtaBins,
0015                                               int nPhiBins,
0016                                               const std::vector<reco::CaloCluster> &layerClusters,
0017                                               const std::vector<float> &mask,
0018                                               const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
0019                                               int deltaIEta,
0020                                               int deltaIPhi,
0021                                               float minCosTheta,
0022                                               float minCosPointing,
0023                                               float root_doublet_max_distance_from_seed_squared,
0024                                               float etaLimitIncreaseWindow,
0025                                               int skip_layers,
0026                                               int maxNumberOfLayers,
0027                                               float maxDeltaTime,
0028                                               int lastLayerEE,
0029                                               int lastLayerFH,
0030                                               const std::vector<double> &siblings_maxRSquared) {
0031   isOuterClusterOfDoublets_.clear();
0032   isOuterClusterOfDoublets_.resize(layerClusters.size());
0033   allDoublets_.clear();
0034   theRootDoublets_.clear();
0035   bool checkDistanceRootDoubletVsSeed = root_doublet_max_distance_from_seed_squared < 9999;
0036   float origin_eta;
0037   float origin_phi;
0038   float maxRSquared;
0039   for (const auto &r : regions) {
0040     bool isGlobal = (r.index == -1);
0041     auto zSide = r.zSide;
0042     int startEtaBin, endEtaBin, startPhiBin, endPhiBin;
0043 
0044     if (isGlobal) {
0045       startEtaBin = 0;
0046       startPhiBin = 0;
0047       endEtaBin = nEtaBins;
0048       endPhiBin = nPhiBins;
0049       origin_eta = 0;
0050       origin_phi = 0;
0051     } else {
0052       auto firstLayerOnZSide = maxNumberOfLayers * zSide;
0053       const auto &firstLayerHisto = histo[firstLayerOnZSide];
0054       origin_eta = r.origin.eta();
0055       origin_phi = r.origin.phi();
0056       int entryEtaBin = firstLayerHisto.etaBin(origin_eta);
0057       int entryPhiBin = firstLayerHisto.phiBin(origin_phi);
0058       // For track-seeded iterations, if the impact point is below a certain
0059       // eta-threshold, i.e., it has higher eta, make the initial search
0060       // window bigger in both eta and phi by one bin, to contain better low
0061       // energy showers.
0062       auto etaWindow = deltaIEta;
0063       auto phiWindow = deltaIPhi;
0064       if (std::abs(origin_eta) > etaLimitIncreaseWindow) {
0065         etaWindow++;
0066         phiWindow++;
0067         LogDebug("HGCGraph") << "Limit of Eta for increase: " << etaLimitIncreaseWindow
0068                              << " reached! Increasing inner search window" << std::endl;
0069       }
0070       startEtaBin = std::max(entryEtaBin - etaWindow, 0);
0071       endEtaBin = std::min(entryEtaBin + etaWindow + 1, nEtaBins);
0072       startPhiBin = entryPhiBin - phiWindow;
0073       endPhiBin = entryPhiBin + phiWindow + 1;
0074       if (verbosity_ > ticl::VerbosityLevel::Guru) {
0075         LogDebug("HGCGraph") << " Entrance eta, phi: " << origin_eta << ", " << origin_phi
0076                              << " entryEtaBin: " << entryEtaBin << " entryPhiBin: " << entryPhiBin
0077                              << " globalBin: " << firstLayerHisto.globalBin(origin_eta, origin_phi)
0078                              << " on layer: " << firstLayerOnZSide << " startEtaBin: " << startEtaBin
0079                              << " endEtaBin: " << endEtaBin << " startPhiBin: " << startPhiBin
0080                              << " endPhiBin: " << endPhiBin << " phiBin(0): " << firstLayerHisto.phiBin(0.)
0081                              << " phiBin(" << M_PI / 2. << "): " << firstLayerHisto.phiBin(M_PI / 2.) << " phiBin("
0082                              << M_PI << "): " << firstLayerHisto.phiBin(M_PI) << " phiBin(" << -M_PI / 2.
0083                              << "): " << firstLayerHisto.phiBin(-M_PI / 2.) << " phiBin(" << -M_PI
0084                              << "): " << firstLayerHisto.phiBin(-M_PI) << " phiBin(" << 2. * M_PI
0085                              << "): " << firstLayerHisto.phiBin(2. * M_PI) << std::endl;
0086       }
0087     }
0088 
0089     for (int il = 0; il < maxNumberOfLayers - 1; ++il) {
0090       for (int outer_layer = 0; outer_layer < std::min(1 + skip_layers, maxNumberOfLayers - 1 - il); ++outer_layer) {
0091         int currentInnerLayerId = il + maxNumberOfLayers * zSide;
0092         int currentOuterLayerId = currentInnerLayerId + 1 + outer_layer;
0093         auto const &outerLayerHisto = histo[currentOuterLayerId];
0094         auto const &innerLayerHisto = histo[currentInnerLayerId];
0095         maxRSquared = (il <= lastLayerEE)   ? siblings_maxRSquared[0]
0096                       : (il <= lastLayerFH) ? siblings_maxRSquared[1]
0097                                             : siblings_maxRSquared[2];
0098         const int etaLimitIncreaseWindowBin = innerLayerHisto.etaBin(etaLimitIncreaseWindow);
0099         if (verbosity_ > ticl::VerbosityLevel::Advanced) {
0100           LogDebug("HGCGraph") << "Limit of Eta for increase: " << etaLimitIncreaseWindow
0101                                << " at etaBin: " << etaLimitIncreaseWindowBin << std::endl;
0102         }
0103 
0104         for (int ieta = startEtaBin; ieta < endEtaBin; ++ieta) {
0105           auto offset = ieta * nPhiBins;
0106           for (int iphi_it = startPhiBin; iphi_it < endPhiBin; ++iphi_it) {
0107             int iphi = ((iphi_it % nPhiBins + nPhiBins) % nPhiBins);
0108             if (verbosity_ > ticl::VerbosityLevel::Advanced) {
0109               LogDebug("HGCGraph") << "Inner Global Bin: " << (offset + iphi)
0110                                    << " on layers I/O: " << currentInnerLayerId << "/" << currentOuterLayerId
0111                                    << " with clusters: " << innerLayerHisto[offset + iphi].size() << std::endl;
0112             }
0113             for (auto innerClusterId : innerLayerHisto[offset + iphi]) {
0114               // Skip masked clusters
0115               if (mask[innerClusterId] == 0.) {
0116                 if (verbosity_ > ticl::VerbosityLevel::Advanced)
0117                   LogDebug("HGCGraph") << "Skipping inner masked cluster " << innerClusterId << std::endl;
0118                 continue;
0119               }
0120 
0121               // For global-seeded iterations, if the inner cluster is above a certain
0122               // eta-threshold, i.e., it has higher eta, make the outer search
0123               // window bigger in both eta and phi by one bin, to contain better low
0124               // energy showers. Track-Seeded iterations are excluded since, in
0125               // that case, the inner search window has already been enlarged.
0126               auto etaWindow = deltaIEta;
0127               auto phiWindow = deltaIPhi;
0128               if (isGlobal && ieta > etaLimitIncreaseWindowBin) {
0129                 etaWindow++;
0130                 phiWindow++;
0131                 if (verbosity_ > ticl::VerbosityLevel::Advanced) {
0132                   LogDebug("HGCGraph") << "Eta and Phi window increased by one" << std::endl;
0133                 }
0134               }
0135               const auto etaRangeMin = std::max(0, ieta - etaWindow);
0136               const auto etaRangeMax = std::min(ieta + etaWindow + 1, nEtaBins);
0137 
0138               for (int oeta = etaRangeMin; oeta < etaRangeMax; ++oeta) {
0139                 // wrap phi bin
0140                 for (int phiRange = 0; phiRange < 2 * phiWindow + 1; ++phiRange) {
0141                   // The first wrapping is to take into account the
0142                   // cases in which we would have to seach in
0143                   // negative bins. The second wrap is mandatory to
0144                   // account for all other cases, since we add in
0145                   // between a full nPhiBins slot.
0146                   auto ophi = ((iphi + phiRange - phiWindow) % nPhiBins + nPhiBins) % nPhiBins;
0147                   if (verbosity_ > ticl::VerbosityLevel::Guru) {
0148                     LogDebug("HGCGraph") << "Outer Global Bin: " << (oeta * nPhiBins + ophi)
0149                                          << " on layers I/O: " << currentInnerLayerId << "/" << currentOuterLayerId
0150                                          << " with clusters: " << innerLayerHisto[oeta * nPhiBins + ophi].size()
0151                                          << std::endl;
0152                   }
0153                   for (auto outerClusterId : outerLayerHisto[oeta * nPhiBins + ophi]) {
0154                     // Skip masked clusters
0155                     if (mask[outerClusterId] == 0.) {
0156                       if (verbosity_ > ticl::VerbosityLevel::Advanced)
0157                         LogDebug("HGCGraph") << "Skipping outer masked cluster " << outerClusterId << std::endl;
0158                       continue;
0159                     }
0160                     auto doubletId = allDoublets_.size();
0161                     if (maxDeltaTime != -1 &&
0162                         !areTimeCompatible(innerClusterId, outerClusterId, layerClustersTime, maxDeltaTime)) {
0163                       if (verbosity_ > ticl::VerbosityLevel::Advanced)
0164                         LogDebug("HGCGraph") << "Rejecting doublets due to timing!" << std::endl;
0165                       continue;
0166                     }
0167                     if (currentOuterLayerId - currentInnerLayerId == 1) {
0168                       if (areOverlappingOnSiblingLayers(innerClusterId, outerClusterId, layerClusters, maxRSquared)) {
0169                         allDoublets_.emplace_back(
0170                             innerClusterId, outerClusterId, doubletId, &layerClusters, r.index, true);
0171                       } else {
0172                         continue;
0173                       }
0174                     } else {
0175                       allDoublets_.emplace_back(
0176                           innerClusterId, outerClusterId, doubletId, &layerClusters, r.index, false);
0177                     }
0178                     if (verbosity_ > ticl::VerbosityLevel::Advanced) {
0179                       LogDebug("HGCGraph")
0180                           << "Creating doubletsId: " << doubletId << " layerLink in-out: [" << currentInnerLayerId
0181                           << ", " << currentOuterLayerId << "] clusterLink in-out: [" << innerClusterId << ", "
0182                           << outerClusterId << "]" << std::endl;
0183                     }
0184                     isOuterClusterOfDoublets_[outerClusterId].push_back(doubletId);
0185                     auto &neigDoublets = isOuterClusterOfDoublets_[innerClusterId];
0186                     auto &thisDoublet = allDoublets_[doubletId];
0187                     if (verbosity_ > ticl::VerbosityLevel::Expert) {
0188                       LogDebug("HGCGraph")
0189                           << "Checking compatibility of doubletId: " << doubletId
0190                           << " with all possible inners doublets link by the innerClusterId: " << innerClusterId
0191                           << std::endl;
0192                     }
0193                     bool isRootDoublet =
0194                         thisDoublet.checkCompatibilityAndTag(allDoublets_,
0195                                                              neigDoublets,
0196                                                              r.directionAtOrigin,
0197                                                              minCosTheta,
0198                                                              minCosPointing,
0199                                                              verbosity_ > ticl::VerbosityLevel::Advanced);
0200                     if (isRootDoublet and checkDistanceRootDoubletVsSeed) {
0201                       if (reco::deltaR2(layerClusters[innerClusterId].eta(),
0202                                         layerClusters[innerClusterId].phi(),
0203                                         origin_eta,
0204                                         origin_phi) > root_doublet_max_distance_from_seed_squared) {
0205                         isRootDoublet = false;
0206                       }
0207                     }
0208                     if (isRootDoublet) {
0209                       theRootDoublets_.push_back(doubletId);
0210                     }
0211                   }
0212                 }
0213               }
0214             }
0215           }
0216         }
0217       }
0218     }
0219   }
0220   // #ifdef FP_DEBUG
0221   if (verbosity_ > ticl::VerbosityLevel::None) {
0222     LogDebug("HGCGraph") << "number of Root doublets " << theRootDoublets_.size() << " over a total number of doublets "
0223                          << allDoublets_.size() << std::endl;
0224   }
0225   // #endif
0226 }
0227 
0228 template <typename TILES>
0229 bool HGCGraphT<TILES>::areTimeCompatible(int innerIdx,
0230                                          int outerIdx,
0231                                          const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
0232                                          float maxDeltaTime) {
0233   float timeIn = layerClustersTime.get(innerIdx).first;
0234   float timeInE = layerClustersTime.get(innerIdx).second;
0235   float timeOut = layerClustersTime.get(outerIdx).first;
0236   float timeOutE = layerClustersTime.get(outerIdx).second;
0237 
0238   return (timeIn == -99. || timeOut == -99. ||
0239           std::abs(timeIn - timeOut) < maxDeltaTime * sqrt(timeInE * timeInE + timeOutE * timeOutE));
0240 }
0241 
0242 template <typename TILES>
0243 bool HGCGraphT<TILES>::areOverlappingOnSiblingLayers(int innerIdx,
0244                                                      int outerIdx,
0245                                                      const std::vector<reco::CaloCluster> &layerClusters,
0246                                                      float maxRSquared) {
0247   return reco::deltaR2(layerClusters[outerIdx].eta(),
0248                        layerClusters[outerIdx].phi(),
0249                        layerClusters[innerIdx].eta(),
0250                        layerClusters[innerIdx].phi()) < maxRSquared;
0251 }
0252 
0253 //also return a vector of seedIndex for the reconstructed tracksters
0254 template <typename TILES>
0255 void HGCGraphT<TILES>::findNtuplets(std::vector<HGCDoublet::HGCntuplet> &foundNtuplets,
0256                                     std::vector<int> &seedIndices,
0257                                     const unsigned int minClustersPerNtuplet,
0258                                     const bool outInDFS,
0259                                     unsigned int maxOutInHops) {
0260   HGCDoublet::HGCntuplet tmpNtuplet;
0261   tmpNtuplet.reserve(minClustersPerNtuplet);
0262   std::vector<std::pair<unsigned int, unsigned int>> outInToVisit;
0263   for (auto rootDoublet : theRootDoublets_) {
0264     tmpNtuplet.clear();
0265     outInToVisit.clear();
0266     int seedIndex = allDoublets_[rootDoublet].seedIndex();
0267     int outInHops = 0;
0268     allDoublets_[rootDoublet].findNtuplets(
0269         allDoublets_, tmpNtuplet, seedIndex, outInDFS, outInHops, maxOutInHops, outInToVisit);
0270     while (!outInToVisit.empty()) {
0271       allDoublets_[outInToVisit.back().first].findNtuplets(
0272           allDoublets_, tmpNtuplet, seedIndex, outInDFS, outInToVisit.back().second, maxOutInHops, outInToVisit);
0273       outInToVisit.pop_back();
0274     }
0275 
0276     if (tmpNtuplet.size() > minClustersPerNtuplet) {
0277       foundNtuplets.push_back(tmpNtuplet);
0278       seedIndices.push_back(seedIndex);
0279     }
0280   }
0281 }
0282 
0283 template class HGCGraphT<TICLLayerTiles>;
0284 template class HGCGraphT<TICLLayerTilesHFNose>;
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!