TICL/plugins/PatternRecognitionbyCA.cc

0001 // Author: Felice Pantaleo, Marco Rovere - felice.pantaleo@cern.ch, marco.rovere@cern.ch
0002 // Date: 11/2018
0003 #include <algorithm>
0004 #include <set>
0005 #include <vector>
0006
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "PatternRecognitionbyCA.h"
0010
0011 #include "TrackstersPCA.h"
0012 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0013 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015
0016 using namespace ticl;
0017
0018 template <typename TILES>
0019 PatternRecognitionbyCA<TILES>::PatternRecognitionbyCA(const edm::ParameterSet &conf, edm::ConsumesCollector iC)
0020     : PatternRecognitionAlgoBaseT<TILES>(conf, iC),
0021       caloGeomToken_(iC.esConsumes<CaloGeometry, CaloGeometryRecord>()),
0022       theGraph_(std::make_unique<HGCGraphT<TILES>>()),
0023       oneTracksterPerTrackSeed_(conf.getParameter<bool>("oneTracksterPerTrackSeed")),
0024       promoteEmptyRegionToTrackster_(conf.getParameter<bool>("promoteEmptyRegionToTrackster")),
0025       out_in_dfs_(conf.getParameter<bool>("out_in_dfs")),
0026       max_out_in_hops_(conf.getParameter<int>("max_out_in_hops")),
0027       min_cos_theta_(conf.getParameter<double>("min_cos_theta")),
0028       min_cos_pointing_(conf.getParameter<double>("min_cos_pointing")),
0029       root_doublet_max_distance_from_seed_squared_(
0030           conf.getParameter<double>("root_doublet_max_distance_from_seed_squared")),
0031       etaLimitIncreaseWindow_(conf.getParameter<double>("etaLimitIncreaseWindow")),
0032       skip_layers_(conf.getParameter<int>("skip_layers")),
0033       max_missing_layers_in_trackster_(conf.getParameter<int>("max_missing_layers_in_trackster")),
0034       check_missing_layers_(max_missing_layers_in_trackster_ < 100),
0035       shower_start_max_layer_(conf.getParameter<int>("shower_start_max_layer")),
0036       min_layers_per_trackster_(conf.getParameter<int>("min_layers_per_trackster")),
0037       filter_on_categories_(conf.getParameter<std::vector<int>>("filter_on_categories")),
0038       pid_threshold_(conf.getParameter<double>("pid_threshold")),
0039       energy_em_over_total_threshold_(conf.getParameter<double>("energy_em_over_total_threshold")),
0040       max_longitudinal_sigmaPCA_(conf.getParameter<double>("max_longitudinal_sigmaPCA")),
0041       min_clusters_per_ntuplet_(min_layers_per_trackster_),
0042       max_delta_time_(conf.getParameter<double>("max_delta_time")),
0043       computeLocalTime_(conf.getParameter<bool>("computeLocalTime")),
0044       siblings_maxRSquared_(conf.getParameter<std::vector<double>>("siblings_maxRSquared")){};
0045
0046 template <typename TILES>
0047 PatternRecognitionbyCA<TILES>::~PatternRecognitionbyCA(){};
0048
0049 template <typename TILES>
0050 void PatternRecognitionbyCA<TILES>::makeTracksters(
0051     const typename PatternRecognitionAlgoBaseT<TILES>::Inputs &input,
0052     std::vector<Trackster> &result,
0053     std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) {
0054   // Protect from events with no seeding regions
0055   if (input.regions.empty())
0056     return;
0057
0058   edm::EventSetup const &es = input.es;
0059   const CaloGeometry &geom = es.getData(caloGeomToken_);
0060   rhtools_.setGeometry(geom);
0061
0062   theGraph_->setVerbosity(PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_);
0063   theGraph_->clear();
0064   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::None) {
0065     LogDebug("HGCPatternRecoByCA") << "Making Tracksters with CA" << std::endl;
0066   }
0067
0068   constexpr auto isHFnose = std::is_same<TILES, TICLLayerTilesHFNose>::value;
0069   constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0070   constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0071
0072   std::vector<HGCDoublet::HGCntuplet> foundNtuplets;
0073   std::vector<int> seedIndices;
0074   std::vector<uint8_t> layer_cluster_usage(input.layerClusters.size(), 0);
0075   theGraph_->makeAndConnectDoublets(input.tiles,
0076                                     input.regions,
0077                                     nEtaBin,
0078                                     nPhiBin,
0079                                     input.layerClusters,
0080                                     input.mask,
0081                                     input.layerClustersTime,
0082                                     1,
0083                                     1,
0084                                     min_cos_theta_,
0085                                     min_cos_pointing_,
0086                                     root_doublet_max_distance_from_seed_squared_,
0087                                     etaLimitIncreaseWindow_,
0088                                     skip_layers_,
0089                                     rhtools_.lastLayer(isHFnose),
0090                                     max_delta_time_,
0091                                     rhtools_.lastLayerEE(isHFnose),
0092                                     rhtools_.lastLayerFH(),
0093                                     siblings_maxRSquared_);
0094
0095   theGraph_->findNtuplets(foundNtuplets, seedIndices, min_clusters_per_ntuplet_, out_in_dfs_, max_out_in_hops_);
0096   //#ifdef FP_DEBUG
0097   const auto &doublets = theGraph_->getAllDoublets();
0098   int tracksterId = -1;
0099
0100   // container for holding tracksters before selection
0101   result.reserve(foundNtuplets.size());
0102
0103   for (auto const &ntuplet : foundNtuplets) {
0104     tracksterId++;
0105
0106     std::set<unsigned int> effective_cluster_idx;
0107
0108     for (auto const &doublet : ntuplet) {
0109       auto innerCluster = doublets[doublet].innerClusterId();
0110       auto outerCluster = doublets[doublet].outerClusterId();
0111
0112       effective_cluster_idx.insert(innerCluster);
0113       effective_cluster_idx.insert(outerCluster);
0114
0115       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0116         LogDebug("HGCPatternRecoByCA") << " New doublet " << doublet << " for trackster: " << result.size()
0117                                        << " InnerCl " << innerCluster << " " << input.layerClusters[innerCluster].x()
0118                                        << " " << input.layerClusters[innerCluster].y() << " "
0119                                        << input.layerClusters[innerCluster].z() << " OuterCl " << outerCluster << " "
0120                                        << input.layerClusters[outerCluster].x() << " "
0121                                        << input.layerClusters[outerCluster].y() << " "
0122                                        << input.layerClusters[outerCluster].z() << " " << tracksterId << std::endl;
0123       }
0124     }
0125     unsigned showerMinLayerId = 99999;
0126     std::vector<unsigned int> uniqueLayerIds;
0127     uniqueLayerIds.reserve(effective_cluster_idx.size());
0128     std::vector<std::pair<unsigned int, unsigned int>> lcIdAndLayer;
0129     lcIdAndLayer.reserve(effective_cluster_idx.size());
0130     for (auto const i : effective_cluster_idx) {
0131       auto const &haf = input.layerClusters[i].hitsAndFractions();
0132       auto layerId = rhtools_.getLayerWithOffset(haf[0].first);
0133       showerMinLayerId = std::min(layerId, showerMinLayerId);
0134       uniqueLayerIds.push_back(layerId);
0135       lcIdAndLayer.emplace_back(i, layerId);
0136     }
0137     std::sort(uniqueLayerIds.begin(), uniqueLayerIds.end());
0138     uniqueLayerIds.erase(std::unique(uniqueLayerIds.begin(), uniqueLayerIds.end()), uniqueLayerIds.end());
0139     unsigned int numberOfLayersInTrackster = uniqueLayerIds.size();
0140     if (check_missing_layers_) {
0141       int numberOfMissingLayers = 0;
0142       unsigned int j = showerMinLayerId;
0143       unsigned int indexInVec = 0;
0144       for (const auto &layer : uniqueLayerIds) {
0145         if (layer != j) {
0146           numberOfMissingLayers++;
0147           j++;
0148           if (numberOfMissingLayers > max_missing_layers_in_trackster_) {
0149             numberOfLayersInTrackster = indexInVec;
0150             for (auto &llpair : lcIdAndLayer) {
0151               if (llpair.second >= layer) {
0152                 effective_cluster_idx.erase(llpair.first);
0153               }
0154             }
0155             break;
0156           }
0157         }
0158         indexInVec++;
0159         j++;
0160       }
0161     }
0162     if ((numberOfLayersInTrackster >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_)) {
0163       // Put back indices, in the form of a Trackster, into the results vector
0164       Trackster tmp;
0165       tmp.vertices().reserve(effective_cluster_idx.size());
0166       tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 1);
0167       //regions and seedIndices can have different size
0168       //if a seeding region does not lead to any trackster
0169       tmp.setSeed(input.regions[0].collectionID, seedIndices[tracksterId]);
0170
0171       std::copy(std::begin(effective_cluster_idx), std::end(effective_cluster_idx), std::back_inserter(tmp.vertices()));
0172       result.push_back(tmp);
0173     }
0174   }
0175   ticl::assignPCAtoTracksters(result,
0176                               input.layerClusters,
0177                               input.layerClustersTime,
0178                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(isHFnose), isHFnose).z(),
0179                               rhtools_,
0180                               computeLocalTime_);
0181
0182   theGraph_->clear();
0183 }
0184
0185 template <typename TILES>
0186 void PatternRecognitionbyCA<TILES>::filter(std::vector<Trackster> &output,
0187                                            const std::vector<Trackster> &inTracksters,
0188                                            const typename PatternRecognitionAlgoBaseT<TILES>::Inputs &input,
0189                                            std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) {
0190   auto filter_on_pids = [&](const ticl::Trackster &t) -> bool {
0191     auto cumulative_prob = 0.;
0192     for (auto index : filter_on_categories_) {
0193       cumulative_prob += t.id_probabilities(index);
0194     }
0195     return (cumulative_prob <= pid_threshold_) &&
0196            (t.raw_em_energy() < energy_em_over_total_threshold_ * t.raw_energy());
0197   };
0198
0199   std::vector<unsigned int> selectedTrackstersIds;
0200   for (unsigned i = 0; i < inTracksters.size(); ++i) {
0201     auto &t = inTracksters[i];
0202     if (!filter_on_pids(t) and t.sigmasPCA()[0] < max_longitudinal_sigmaPCA_) {
0203       selectedTrackstersIds.push_back(i);
0204     }
0205   }
0206   output.reserve(selectedTrackstersIds.size());
0207   bool isRegionalIter = !input.regions.empty() && (input.regions[0].index != -1);
0208   for (unsigned i = 0; i < selectedTrackstersIds.size(); ++i) {
0209     const auto &t = inTracksters[selectedTrackstersIds[i]];
0210     if (isRegionalIter) {
0211       seedToTracksterAssociation[t.seedIndex()].push_back(i);
0212     }
0213     output.push_back(t);
0214   }
0215
0216   // Now decide if the tracksters from the track-based iterations have to be merged
0217   if (oneTracksterPerTrackSeed_) {
0218     std::vector<Trackster> tmp;
0219     mergeTrackstersTRK(output, input.layerClusters, tmp, seedToTracksterAssociation);
0220     tmp.swap(output);
0221   }
0222
0223   // now adding dummy tracksters from seeds not connected to any shower in the result collection
0224   // these are marked as charged hadrons with probability 1.
0225   if (promoteEmptyRegionToTrackster_) {
0226     emptyTrackstersFromSeedsTRK(output, seedToTracksterAssociation, input.regions[0].collectionID);
0227   }
0228
0229   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0230     for (auto &trackster : output) {
0231       LogDebug("HGCPatternRecoByCA") << "Trackster characteristics: " << std::endl;
0232       LogDebug("HGCPatternRecoByCA") << "Size: " << trackster.vertices().size() << std::endl;
0233       auto counter = 0;
0234       for (auto const &p : trackster.id_probabilities()) {
0235         LogDebug("HGCPatternRecoByCA") << counter++ << ": " << p << std::endl;
0236       }
0237     }
0238   }
0239 }
0240 template <typename TILES>
0241 void PatternRecognitionbyCA<TILES>::mergeTrackstersTRK(
0242     const std::vector<Trackster> &input,
0243     const std::vector<reco::CaloCluster> &layerClusters,
0244     std::vector<Trackster> &output,
0245     std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) const {
0246   output.reserve(input.size());
0247   for (auto &thisSeed : seedToTracksterAssociation) {
0248     auto &tracksters = thisSeed.second;
0249     if (!tracksters.empty()) {
0250       auto numberOfTrackstersInSeed = tracksters.size();
0251       output.emplace_back(input[tracksters[0]]);
0252       auto &outTrackster = output.back();
0253       tracksters[0] = output.size() - 1;
0254       auto updated_size = outTrackster.vertices().size();
0255       for (unsigned int j = 1; j < numberOfTrackstersInSeed; ++j) {
0256         auto &thisTrackster = input[tracksters[j]];
0257         updated_size += thisTrackster.vertices().size();
0258         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Basic) {
0259           LogDebug("HGCPatternRecoByCA") << "Updated size: " << updated_size << std::endl;
0260         }
0261         outTrackster.vertices().reserve(updated_size);
0262         outTrackster.vertex_multiplicity().reserve(updated_size);
0263         std::copy(std::begin(thisTrackster.vertices()),
0264                   std::end(thisTrackster.vertices()),
0265                   std::back_inserter(outTrackster.vertices()));
0266         std::copy(std::begin(thisTrackster.vertex_multiplicity()),
0267                   std::end(thisTrackster.vertex_multiplicity()),
0268                   std::back_inserter(outTrackster.vertex_multiplicity()));
0269       }
0270       tracksters.resize(1);
0271
0272       // Find duplicate LCs
0273       auto &orig_vtx = outTrackster.vertices();
0274       auto vtx_sorted{orig_vtx};
0275       std::sort(std::begin(vtx_sorted), std::end(vtx_sorted));
0276       for (unsigned int iLC = 1; iLC < vtx_sorted.size(); ++iLC) {
0277         if (vtx_sorted[iLC] == vtx_sorted[iLC - 1]) {
0278           // Clean up duplicate LCs
0279           const auto lcIdx = vtx_sorted[iLC];
0280           const auto firstEl = std::find(orig_vtx.begin(), orig_vtx.end(), lcIdx);
0281           const auto firstPos = std::distance(std::begin(orig_vtx), firstEl);
0282           auto iDup = std::find(std::next(firstEl), orig_vtx.end(), lcIdx);
0283           while (iDup != orig_vtx.end()) {
0284             orig_vtx.erase(iDup);
0285             outTrackster.vertex_multiplicity().erase(outTrackster.vertex_multiplicity().begin() +
0286                                                      std::distance(std::begin(orig_vtx), iDup));
0287             outTrackster.vertex_multiplicity()[firstPos] -= 1;
0288             iDup = std::find(std::next(firstEl), orig_vtx.end(), lcIdx);
0289           };
0290         }
0291       }
0292     }
0293   }
0294   output.shrink_to_fit();
0295 }
0296
0297 template <typename TILES>
0298 void PatternRecognitionbyCA<TILES>::emptyTrackstersFromSeedsTRK(
0299     std::vector<Trackster> &tracksters,
0300     std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation,
0301     const edm::ProductID &collectionID) const {
0302   for (auto &thisSeed : seedToTracksterAssociation) {
0303     if (thisSeed.second.empty()) {
0304       Trackster t;
0305       t.setRegressedEnergy(0.f);
0306       t.zeroProbabilities();
0307       t.setIdProbability(ticl::Trackster::ParticleType::charged_hadron, 1.f);
0308       t.setSeed(collectionID, thisSeed.first);
0309       tracksters.emplace_back(t);
0310       thisSeed.second.emplace_back(tracksters.size() - 1);
0311     }
0312   }
0313 }
0314
0315 template <typename TILES>
0316 void PatternRecognitionbyCA<TILES>::fillPSetDescription(edm::ParameterSetDescription &iDesc) {
0317   iDesc.add<int>("algo_verbosity", 0);
0318   iDesc.add<bool>("oneTracksterPerTrackSeed", false);
0319   iDesc.add<bool>("promoteEmptyRegionToTrackster", false);
0320   iDesc.add<bool>("out_in_dfs", true);
0321   iDesc.add<int>("max_out_in_hops", 10);
0322   iDesc.add<double>("min_cos_theta", 0.915);
0323   iDesc.add<double>("min_cos_pointing", -1.);
0324   iDesc.add<double>("root_doublet_max_distance_from_seed_squared", 9999);
0325   iDesc.add<double>("etaLimitIncreaseWindow", 2.1);
0326   iDesc.add<int>("skip_layers", 0);
0327   iDesc.add<int>("max_missing_layers_in_trackster", 9999);
0328   iDesc.add<int>("shower_start_max_layer", 9999)->setComment("make default such that no filtering is applied");
0329   iDesc.add<int>("min_layers_per_trackster", 10);
0330   iDesc.add<std::vector<int>>("filter_on_categories", {0});
0331   iDesc.add<double>("pid_threshold", 0.)->setComment("make default such that no filtering is applied");
0332   iDesc.add<double>("energy_em_over_total_threshold", -1.)
0333       ->setComment("make default such that no filtering is applied");
0334   iDesc.add<double>("max_longitudinal_sigmaPCA", 9999);
0335   iDesc.add<double>("max_delta_time", 3.)->setComment("nsigma");
0336   iDesc.add<bool>("computeLocalTime", false);
0337   iDesc.add<std::vector<double>>("siblings_maxRSquared", {6e-4, 6e-4, 6e-4});
0338 }
0339
0340 template class ticl::PatternRecognitionbyCA<TICLLayerTiles>;
0341 template class ticl::PatternRecognitionbyCA<TICLLayerTilesHFNose>;