File indexing completed on 2023-03-17 11:18:04
0001
0002
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 eidInputName_(conf.getParameter<std::string>("eid_input_name")),
0044 eidOutputNameEnergy_(conf.getParameter<std::string>("eid_output_name_energy")),
0045 eidOutputNameId_(conf.getParameter<std::string>("eid_output_name_id")),
0046 eidMinClusterEnergy_(conf.getParameter<double>("eid_min_cluster_energy")),
0047 eidNLayers_(conf.getParameter<int>("eid_n_layers")),
0048 eidNClusters_(conf.getParameter<int>("eid_n_clusters")),
0049 siblings_maxRSquared_(conf.getParameter<std::vector<double>>("siblings_maxRSquared")){};
0050
0051 template <typename TILES>
0052 PatternRecognitionbyCA<TILES>::~PatternRecognitionbyCA(){};
0053
0054 template <typename TILES>
0055 void PatternRecognitionbyCA<TILES>::makeTracksters(
0056 const typename PatternRecognitionAlgoBaseT<TILES>::Inputs &input,
0057 std::vector<Trackster> &result,
0058 std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) {
0059
0060 if (input.regions.empty())
0061 return;
0062
0063 edm::EventSetup const &es = input.es;
0064 const CaloGeometry &geom = es.getData(caloGeomToken_);
0065 rhtools_.setGeometry(geom);
0066
0067 theGraph_->setVerbosity(PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_);
0068 theGraph_->clear();
0069 if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::None) {
0070 LogDebug("HGCPatternRecoByCA") << "Making Tracksters with CA" << std::endl;
0071 }
0072
0073 constexpr auto isHFnose = std::is_same<TILES, TICLLayerTilesHFNose>::value;
0074 constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0075 constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0076
0077 bool isRegionalIter = (input.regions[0].index != -1);
0078 std::vector<HGCDoublet::HGCntuplet> foundNtuplets;
0079 std::vector<int> seedIndices;
0080 std::vector<uint8_t> layer_cluster_usage(input.layerClusters.size(), 0);
0081 theGraph_->makeAndConnectDoublets(input.tiles,
0082 input.regions,
0083 nEtaBin,
0084 nPhiBin,
0085 input.layerClusters,
0086 input.mask,
0087 input.layerClustersTime,
0088 1,
0089 1,
0090 min_cos_theta_,
0091 min_cos_pointing_,
0092 root_doublet_max_distance_from_seed_squared_,
0093 etaLimitIncreaseWindow_,
0094 skip_layers_,
0095 rhtools_.lastLayer(isHFnose),
0096 max_delta_time_,
0097 rhtools_.lastLayerEE(isHFnose),
0098 rhtools_.lastLayerFH(),
0099 siblings_maxRSquared_);
0100
0101 theGraph_->findNtuplets(foundNtuplets, seedIndices, min_clusters_per_ntuplet_, out_in_dfs_, max_out_in_hops_);
0102
0103 const auto &doublets = theGraph_->getAllDoublets();
0104 int tracksterId = -1;
0105
0106
0107 std::vector<Trackster> tmpTracksters;
0108 tmpTracksters.reserve(foundNtuplets.size());
0109
0110 for (auto const &ntuplet : foundNtuplets) {
0111 tracksterId++;
0112
0113 std::set<unsigned int> effective_cluster_idx;
0114
0115 for (auto const &doublet : ntuplet) {
0116 auto innerCluster = doublets[doublet].innerClusterId();
0117 auto outerCluster = doublets[doublet].outerClusterId();
0118
0119 effective_cluster_idx.insert(innerCluster);
0120 effective_cluster_idx.insert(outerCluster);
0121
0122 if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0123 LogDebug("HGCPatternRecoByCA") << " New doublet " << doublet << " for trackster: " << result.size()
0124 << " InnerCl " << innerCluster << " " << input.layerClusters[innerCluster].x()
0125 << " " << input.layerClusters[innerCluster].y() << " "
0126 << input.layerClusters[innerCluster].z() << " OuterCl " << outerCluster << " "
0127 << input.layerClusters[outerCluster].x() << " "
0128 << input.layerClusters[outerCluster].y() << " "
0129 << input.layerClusters[outerCluster].z() << " " << tracksterId << std::endl;
0130 }
0131 }
0132 unsigned showerMinLayerId = 99999;
0133 std::vector<unsigned int> uniqueLayerIds;
0134 uniqueLayerIds.reserve(effective_cluster_idx.size());
0135 std::vector<std::pair<unsigned int, unsigned int>> lcIdAndLayer;
0136 lcIdAndLayer.reserve(effective_cluster_idx.size());
0137 for (auto const i : effective_cluster_idx) {
0138 auto const &haf = input.layerClusters[i].hitsAndFractions();
0139 auto layerId = rhtools_.getLayerWithOffset(haf[0].first);
0140 showerMinLayerId = std::min(layerId, showerMinLayerId);
0141 uniqueLayerIds.push_back(layerId);
0142 lcIdAndLayer.emplace_back(i, layerId);
0143 }
0144 std::sort(uniqueLayerIds.begin(), uniqueLayerIds.end());
0145 uniqueLayerIds.erase(std::unique(uniqueLayerIds.begin(), uniqueLayerIds.end()), uniqueLayerIds.end());
0146 unsigned int numberOfLayersInTrackster = uniqueLayerIds.size();
0147 if (check_missing_layers_) {
0148 int numberOfMissingLayers = 0;
0149 unsigned int j = showerMinLayerId;
0150 unsigned int indexInVec = 0;
0151 for (const auto &layer : uniqueLayerIds) {
0152 if (layer != j) {
0153 numberOfMissingLayers++;
0154 j++;
0155 if (numberOfMissingLayers > max_missing_layers_in_trackster_) {
0156 numberOfLayersInTrackster = indexInVec;
0157 for (auto &llpair : lcIdAndLayer) {
0158 if (llpair.second >= layer) {
0159 effective_cluster_idx.erase(llpair.first);
0160 }
0161 }
0162 break;
0163 }
0164 }
0165 indexInVec++;
0166 j++;
0167 }
0168 }
0169 if ((numberOfLayersInTrackster >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_)) {
0170
0171 Trackster tmp;
0172 tmp.vertices().reserve(effective_cluster_idx.size());
0173 tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 1);
0174
0175
0176 tmp.setSeed(input.regions[0].collectionID, seedIndices[tracksterId]);
0177
0178 std::copy(std::begin(effective_cluster_idx), std::end(effective_cluster_idx), std::back_inserter(tmp.vertices()));
0179 tmpTracksters.push_back(tmp);
0180 }
0181 }
0182 ticl::assignPCAtoTracksters(tmpTracksters,
0183 input.layerClusters,
0184 input.layerClustersTime,
0185 rhtools_.getPositionLayer(rhtools_.lastLayerEE(isHFnose), isHFnose).z());
0186
0187
0188 energyRegressionAndID(input.layerClusters, input.tfSession, tmpTracksters);
0189
0190
0191
0192
0193
0194
0195 auto filter_on_pids = [&](Trackster &t) -> bool {
0196 auto cumulative_prob = 0.;
0197 for (auto index : filter_on_categories_) {
0198 cumulative_prob += t.id_probabilities(index);
0199 }
0200 return (cumulative_prob <= pid_threshold_) &&
0201 (t.raw_em_energy() < energy_em_over_total_threshold_ * t.raw_energy());
0202 };
0203
0204 std::vector<unsigned int> selectedTrackstersIds;
0205 for (unsigned i = 0; i < tmpTracksters.size(); ++i) {
0206 if (!filter_on_pids(tmpTracksters[i]) and tmpTracksters[i].sigmasPCA()[0] < max_longitudinal_sigmaPCA_) {
0207 selectedTrackstersIds.push_back(i);
0208 }
0209 }
0210
0211 result.reserve(selectedTrackstersIds.size());
0212
0213 for (unsigned i = 0; i < selectedTrackstersIds.size(); ++i) {
0214 const auto &t = tmpTracksters[selectedTrackstersIds[i]];
0215 for (auto const lcId : t.vertices()) {
0216 layer_cluster_usage[lcId]++;
0217 if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Basic)
0218 LogDebug("HGCPatternRecoByCA") << "LayerID: " << lcId << " count: " << (int)layer_cluster_usage[lcId]
0219 << std::endl;
0220 }
0221 if (isRegionalIter) {
0222 seedToTracksterAssociation[t.seedIndex()].push_back(i);
0223 }
0224 result.push_back(t);
0225 }
0226
0227 for (auto &trackster : result) {
0228 assert(trackster.vertices().size() <= trackster.vertex_multiplicity().size());
0229 for (size_t i = 0; i < trackster.vertices().size(); ++i) {
0230 trackster.vertex_multiplicity()[i] = layer_cluster_usage[trackster.vertices(i)];
0231 if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Basic)
0232 LogDebug("HGCPatternRecoByCA") << "LayerID: " << trackster.vertices(i)
0233 << " count: " << (int)trackster.vertex_multiplicity(i) << std::endl;
0234 }
0235 }
0236
0237 if (oneTracksterPerTrackSeed_) {
0238 std::vector<Trackster> tmp;
0239 mergeTrackstersTRK(result, input.layerClusters, tmp, seedToTracksterAssociation);
0240 tmp.swap(result);
0241 }
0242
0243 ticl::assignPCAtoTracksters(result,
0244 input.layerClusters,
0245 input.layerClustersTime,
0246 rhtools_.getPositionLayer(rhtools_.lastLayerEE(isHFnose), isHFnose).z());
0247
0248
0249 energyRegressionAndID(input.layerClusters, input.tfSession, result);
0250
0251
0252
0253 if (promoteEmptyRegionToTrackster_) {
0254 emptyTrackstersFromSeedsTRK(result, seedToTracksterAssociation, input.regions[0].collectionID);
0255 }
0256
0257 if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0258 for (auto &trackster : result) {
0259 LogDebug("HGCPatternRecoByCA") << "Trackster characteristics: " << std::endl;
0260 LogDebug("HGCPatternRecoByCA") << "Size: " << trackster.vertices().size() << std::endl;
0261 auto counter = 0;
0262 for (auto const &p : trackster.id_probabilities()) {
0263 LogDebug("HGCPatternRecoByCA") << counter++ << ": " << p << std::endl;
0264 }
0265 }
0266 }
0267 theGraph_->clear();
0268 }
0269
0270 template <typename TILES>
0271 void PatternRecognitionbyCA<TILES>::mergeTrackstersTRK(
0272 const std::vector<Trackster> &input,
0273 const std::vector<reco::CaloCluster> &layerClusters,
0274 std::vector<Trackster> &output,
0275 std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) const {
0276 output.reserve(input.size());
0277 for (auto &thisSeed : seedToTracksterAssociation) {
0278 auto &tracksters = thisSeed.second;
0279 if (!tracksters.empty()) {
0280 auto numberOfTrackstersInSeed = tracksters.size();
0281 output.emplace_back(input[tracksters[0]]);
0282 auto &outTrackster = output.back();
0283 tracksters[0] = output.size() - 1;
0284 auto updated_size = outTrackster.vertices().size();
0285 for (unsigned int j = 1; j < numberOfTrackstersInSeed; ++j) {
0286 auto &thisTrackster = input[tracksters[j]];
0287 updated_size += thisTrackster.vertices().size();
0288 if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Basic) {
0289 LogDebug("HGCPatternRecoByCA") << "Updated size: " << updated_size << std::endl;
0290 }
0291 outTrackster.vertices().reserve(updated_size);
0292 outTrackster.vertex_multiplicity().reserve(updated_size);
0293 std::copy(std::begin(thisTrackster.vertices()),
0294 std::end(thisTrackster.vertices()),
0295 std::back_inserter(outTrackster.vertices()));
0296 std::copy(std::begin(thisTrackster.vertex_multiplicity()),
0297 std::end(thisTrackster.vertex_multiplicity()),
0298 std::back_inserter(outTrackster.vertex_multiplicity()));
0299 }
0300 tracksters.resize(1);
0301
0302
0303 auto &orig_vtx = outTrackster.vertices();
0304 auto vtx_sorted{orig_vtx};
0305 std::sort(std::begin(vtx_sorted), std::end(vtx_sorted));
0306 for (unsigned int iLC = 1; iLC < vtx_sorted.size(); ++iLC) {
0307 if (vtx_sorted[iLC] == vtx_sorted[iLC - 1]) {
0308
0309 const auto lcIdx = vtx_sorted[iLC];
0310 const auto firstEl = std::find(orig_vtx.begin(), orig_vtx.end(), lcIdx);
0311 const auto firstPos = std::distance(std::begin(orig_vtx), firstEl);
0312 auto iDup = std::find(std::next(firstEl), orig_vtx.end(), lcIdx);
0313 while (iDup != orig_vtx.end()) {
0314 orig_vtx.erase(iDup);
0315 outTrackster.vertex_multiplicity().erase(outTrackster.vertex_multiplicity().begin() +
0316 std::distance(std::begin(orig_vtx), iDup));
0317 outTrackster.vertex_multiplicity()[firstPos] -= 1;
0318 iDup = std::find(std::next(firstEl), orig_vtx.end(), lcIdx);
0319 };
0320 }
0321 }
0322 }
0323 }
0324 output.shrink_to_fit();
0325 }
0326
0327 template <typename TILES>
0328 void PatternRecognitionbyCA<TILES>::emptyTrackstersFromSeedsTRK(
0329 std::vector<Trackster> &tracksters,
0330 std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation,
0331 const edm::ProductID &collectionID) const {
0332 for (auto &thisSeed : seedToTracksterAssociation) {
0333 if (thisSeed.second.empty()) {
0334 Trackster t;
0335 t.setRegressedEnergy(0.f);
0336 t.zeroProbabilities();
0337 t.setIdProbability(ticl::Trackster::ParticleType::charged_hadron, 1.f);
0338 t.setSeed(collectionID, thisSeed.first);
0339 tracksters.emplace_back(t);
0340 thisSeed.second.emplace_back(tracksters.size() - 1);
0341 }
0342 }
0343 }
0344
0345 template <typename TILES>
0346 void PatternRecognitionbyCA<TILES>::energyRegressionAndID(const std::vector<reco::CaloCluster> &layerClusters,
0347 const tensorflow::Session *eidSession,
0348 std::vector<Trackster> &tracksters) {
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372 std::vector<int> tracksterIndices;
0373 for (int i = 0; i < (int)tracksters.size(); i++) {
0374
0375
0376
0377 float sumClusterEnergy = 0.;
0378 for (const unsigned int &vertex : tracksters[i].vertices()) {
0379 sumClusterEnergy += (float)layerClusters[vertex].energy();
0380
0381 if (sumClusterEnergy >= eidMinClusterEnergy_) {
0382
0383 tracksters[i].setRegressedEnergy(0.f);
0384 tracksters[i].zeroProbabilities();
0385 tracksterIndices.push_back(i);
0386 break;
0387 }
0388 }
0389 }
0390
0391
0392 int batchSize = (int)tracksterIndices.size();
0393 if (batchSize == 0) {
0394 return;
0395 }
0396
0397
0398 tensorflow::TensorShape shape({batchSize, eidNLayers_, eidNClusters_, eidNFeatures_});
0399 tensorflow::Tensor input(tensorflow::DT_FLOAT, shape);
0400 tensorflow::NamedTensorList inputList = {{eidInputName_, input}};
0401
0402 std::vector<tensorflow::Tensor> outputs;
0403 std::vector<std::string> outputNames;
0404 if (!eidOutputNameEnergy_.empty()) {
0405 outputNames.push_back(eidOutputNameEnergy_);
0406 }
0407 if (!eidOutputNameId_.empty()) {
0408 outputNames.push_back(eidOutputNameId_);
0409 }
0410
0411
0412 for (int i = 0; i < batchSize; i++) {
0413 const Trackster &trackster = tracksters[tracksterIndices[i]];
0414
0415
0416
0417
0418 std::vector<int> clusterIndices(trackster.vertices().size());
0419 for (int k = 0; k < (int)trackster.vertices().size(); k++) {
0420 clusterIndices[k] = k;
0421 }
0422 sort(clusterIndices.begin(), clusterIndices.end(), [&layerClusters, &trackster](const int &a, const int &b) {
0423 return layerClusters[trackster.vertices(a)].energy() > layerClusters[trackster.vertices(b)].energy();
0424 });
0425
0426
0427 std::vector<int> seenClusters(eidNLayers_);
0428
0429
0430 for (const int &k : clusterIndices) {
0431
0432 const reco::CaloCluster &cluster = layerClusters[trackster.vertices(k)];
0433 int j = rhtools_.getLayerWithOffset(cluster.hitsAndFractions()[0].first) - 1;
0434 if (j < eidNLayers_ && seenClusters[j] < eidNClusters_) {
0435
0436 float *features = &input.tensor<float, 4>()(i, j, seenClusters[j], 0);
0437
0438
0439 *(features++) = float(cluster.energy() / float(trackster.vertex_multiplicity(k)));
0440 *(features++) = float(std::abs(cluster.eta()));
0441 *(features) = float(cluster.phi());
0442
0443
0444 seenClusters[j]++;
0445 }
0446 }
0447
0448
0449 for (int j = 0; j < eidNLayers_; j++) {
0450 for (int k = seenClusters[j]; k < eidNClusters_; k++) {
0451 float *features = &input.tensor<float, 4>()(i, j, k, 0);
0452 for (int l = 0; l < eidNFeatures_; l++) {
0453 *(features++) = 0.f;
0454 }
0455 }
0456 }
0457 }
0458
0459
0460 tensorflow::run(eidSession, inputList, outputNames, &outputs);
0461
0462
0463 if (!eidOutputNameEnergy_.empty()) {
0464
0465 float *energy = outputs[0].flat<float>().data();
0466
0467 for (const int &i : tracksterIndices) {
0468 tracksters[i].setRegressedEnergy(*(energy++));
0469 }
0470 }
0471
0472
0473 if (!eidOutputNameId_.empty()) {
0474
0475 int probsIdx = eidOutputNameEnergy_.empty() ? 0 : 1;
0476 float *probs = outputs[probsIdx].flat<float>().data();
0477
0478 for (const int &i : tracksterIndices) {
0479 tracksters[i].setProbabilities(probs);
0480 probs += tracksters[i].id_probabilities().size();
0481 }
0482 }
0483 }
0484
0485 template <typename TILES>
0486 void PatternRecognitionbyCA<TILES>::fillPSetDescription(edm::ParameterSetDescription &iDesc) {
0487 iDesc.add<int>("algo_verbosity", 0);
0488 iDesc.add<bool>("oneTracksterPerTrackSeed", false);
0489 iDesc.add<bool>("promoteEmptyRegionToTrackster", false);
0490 iDesc.add<bool>("out_in_dfs", true);
0491 iDesc.add<int>("max_out_in_hops", 10);
0492 iDesc.add<double>("min_cos_theta", 0.915);
0493 iDesc.add<double>("min_cos_pointing", -1.);
0494 iDesc.add<double>("root_doublet_max_distance_from_seed_squared", 9999);
0495 iDesc.add<double>("etaLimitIncreaseWindow", 2.1);
0496 iDesc.add<int>("skip_layers", 0);
0497 iDesc.add<int>("max_missing_layers_in_trackster", 9999);
0498 iDesc.add<int>("shower_start_max_layer", 9999)->setComment("make default such that no filtering is applied");
0499 iDesc.add<int>("min_layers_per_trackster", 10);
0500 iDesc.add<std::vector<int>>("filter_on_categories", {0});
0501 iDesc.add<double>("pid_threshold", 0.)->setComment("make default such that no filtering is applied");
0502 iDesc.add<double>("energy_em_over_total_threshold", -1.)
0503 ->setComment("make default such that no filtering is applied");
0504 iDesc.add<double>("max_longitudinal_sigmaPCA", 9999);
0505 iDesc.add<double>("max_delta_time", 3.)->setComment("nsigma");
0506 iDesc.add<std::string>("eid_input_name", "input");
0507 iDesc.add<std::string>("eid_output_name_energy", "output/regressed_energy");
0508 iDesc.add<std::string>("eid_output_name_id", "output/id_probabilities");
0509 iDesc.add<double>("eid_min_cluster_energy", 1.);
0510 iDesc.add<int>("eid_n_layers", 50);
0511 iDesc.add<int>("eid_n_clusters", 10);
0512 iDesc.add<std::vector<double>>("siblings_maxRSquared", {6e-4, 6e-4, 6e-4});
0513 }
0514
0515 template class ticl::PatternRecognitionbyCA<TICLLayerTiles>;
0516 template class ticl::PatternRecognitionbyCA<TICLLayerTilesHFNose>;