18 template <
typename TILES>
23 oneTracksterPerTrackSeed_(conf.getParameter<
bool>(
"oneTracksterPerTrackSeed")),
24 promoteEmptyRegionToTrackster_(conf.getParameter<
bool>(
"promoteEmptyRegionToTrackster")),
25 out_in_dfs_(conf.getParameter<
bool>(
"out_in_dfs")),
26 max_out_in_hops_(conf.getParameter<
int>(
"max_out_in_hops")),
27 min_cos_theta_(conf.getParameter<double>(
"min_cos_theta")),
28 min_cos_pointing_(conf.getParameter<double>(
"min_cos_pointing")),
29 root_doublet_max_distance_from_seed_squared_(
30 conf.getParameter<double>(
"root_doublet_max_distance_from_seed_squared")),
31 etaLimitIncreaseWindow_(conf.getParameter<double>(
"etaLimitIncreaseWindow")),
32 skip_layers_(conf.getParameter<
int>(
"skip_layers")),
33 max_missing_layers_in_trackster_(conf.getParameter<
int>(
"max_missing_layers_in_trackster")),
34 check_missing_layers_(max_missing_layers_in_trackster_ < 100),
35 shower_start_max_layer_(conf.getParameter<
int>(
"shower_start_max_layer")),
36 min_layers_per_trackster_(conf.getParameter<
int>(
"min_layers_per_trackster")),
37 filter_on_categories_(conf.getParameter<
std::
vector<
int>>(
"filter_on_categories")),
38 pid_threshold_(conf.getParameter<double>(
"pid_threshold")),
39 energy_em_over_total_threshold_(conf.getParameter<double>(
"energy_em_over_total_threshold")),
40 max_longitudinal_sigmaPCA_(conf.getParameter<double>(
"max_longitudinal_sigmaPCA")),
41 min_clusters_per_ntuplet_(min_layers_per_trackster_),
42 max_delta_time_(conf.getParameter<double>(
"max_delta_time")),
43 eidInputName_(conf.getParameter<
std::
string>(
"eid_input_name")),
44 eidOutputNameEnergy_(conf.getParameter<
std::
string>(
"eid_output_name_energy")),
45 eidOutputNameId_(conf.getParameter<
std::
string>(
"eid_output_name_id")),
46 eidMinClusterEnergy_(conf.getParameter<double>(
"eid_min_cluster_energy")),
47 eidNLayers_(conf.getParameter<
int>(
"eid_n_layers")),
48 eidNClusters_(conf.getParameter<
int>(
"eid_n_clusters")),
49 computeLocalTime_(conf.getParameter<
bool>(
"computeLocalTime")),
50 siblings_maxRSquared_(conf.getParameter<
std::
vector<double>>(
"siblings_maxRSquared")){};
52 template <
typename TILES>
55 template <
typename TILES>
58 std::vector<Trackster> &
result,
59 std::unordered_map<
int, std::vector<int>> &seedToTracksterAssociation) {
61 if (
input.regions.empty())
66 rhtools_.setGeometry(
geom);
71 LogDebug(
"HGCPatternRecoByCA") <<
"Making Tracksters with CA" << std::endl;
78 bool isRegionalIter = (
input.regions[0].index != -1);
80 std::vector<int> seedIndices;
81 std::vector<uint8_t> layer_cluster_usage(
input.layerClusters.size(), 0);
82 theGraph_->makeAndConnectDoublets(
input.tiles,
88 input.layerClustersTime,
93 root_doublet_max_distance_from_seed_squared_,
94 etaLimitIncreaseWindow_,
96 rhtools_.lastLayer(isHFnose),
98 rhtools_.lastLayerEE(isHFnose),
99 rhtools_.lastLayerFH(),
100 siblings_maxRSquared_);
102 theGraph_->findNtuplets(
foundNtuplets, seedIndices, min_clusters_per_ntuplet_, out_in_dfs_, max_out_in_hops_);
104 const auto &
doublets = theGraph_->getAllDoublets();
105 int tracksterId = -1;
108 std::vector<Trackster> tmpTracksters;
114 std::set<unsigned int> effective_cluster_idx;
116 for (
auto const &doublet :
ntuplet) {
117 auto innerCluster =
doublets[doublet].innerClusterId();
118 auto outerCluster =
doublets[doublet].outerClusterId();
120 effective_cluster_idx.insert(innerCluster);
121 effective_cluster_idx.insert(outerCluster);
124 LogDebug(
"HGCPatternRecoByCA") <<
" New doublet " << doublet <<
" for trackster: " <<
result.size()
125 <<
" InnerCl " << innerCluster <<
" " <<
input.layerClusters[innerCluster].x()
126 <<
" " <<
input.layerClusters[innerCluster].y() <<
" " 127 <<
input.layerClusters[innerCluster].z() <<
" OuterCl " << outerCluster <<
" " 128 <<
input.layerClusters[outerCluster].x() <<
" " 129 <<
input.layerClusters[outerCluster].y() <<
" " 130 <<
input.layerClusters[outerCluster].z() <<
" " << tracksterId << std::endl;
133 unsigned showerMinLayerId = 99999;
134 std::vector<unsigned int> uniqueLayerIds;
135 uniqueLayerIds.reserve(effective_cluster_idx.size());
136 std::vector<std::pair<unsigned int, unsigned int>> lcIdAndLayer;
137 lcIdAndLayer.reserve(effective_cluster_idx.size());
138 for (
auto const i : effective_cluster_idx) {
139 auto const &haf =
input.layerClusters[
i].hitsAndFractions();
140 auto layerId = rhtools_.getLayerWithOffset(haf[0].
first);
141 showerMinLayerId =
std::min(layerId, showerMinLayerId);
142 uniqueLayerIds.push_back(layerId);
143 lcIdAndLayer.emplace_back(
i, layerId);
145 std::sort(uniqueLayerIds.begin(), uniqueLayerIds.end());
146 uniqueLayerIds.erase(
std::unique(uniqueLayerIds.begin(), uniqueLayerIds.end()), uniqueLayerIds.end());
147 unsigned int numberOfLayersInTrackster = uniqueLayerIds.size();
148 if (check_missing_layers_) {
149 int numberOfMissingLayers = 0;
150 unsigned int j = showerMinLayerId;
151 unsigned int indexInVec = 0;
152 for (
const auto &layer : uniqueLayerIds) {
154 numberOfMissingLayers++;
156 if (numberOfMissingLayers > max_missing_layers_in_trackster_) {
157 numberOfLayersInTrackster = indexInVec;
158 for (
auto &llpair : lcIdAndLayer) {
159 if (llpair.second >= layer) {
160 effective_cluster_idx.erase(llpair.first);
170 if ((numberOfLayersInTrackster >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_)) {
173 tmp.vertices().reserve(effective_cluster_idx.size());
174 tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 1);
177 tmp.setSeed(
input.regions[0].collectionID, seedIndices[tracksterId]);
179 std::copy(std::begin(effective_cluster_idx), std::end(effective_cluster_idx), std::back_inserter(
tmp.vertices()));
180 tmpTracksters.push_back(
tmp);
185 input.layerClustersTime,
186 rhtools_.getPositionLayer(rhtools_.lastLayerEE(isHFnose), isHFnose).z(),
190 energyRegressionAndID(
input.layerClusters,
input.tfSession, tmpTracksters);
197 auto filter_on_pids = [&](
Trackster &
t) ->
bool {
198 auto cumulative_prob = 0.;
199 for (
auto index : filter_on_categories_) {
200 cumulative_prob +=
t.id_probabilities(
index);
202 return (cumulative_prob <= pid_threshold_) &&
203 (
t.raw_em_energy() < energy_em_over_total_threshold_ *
t.raw_energy());
206 std::vector<unsigned int> selectedTrackstersIds;
207 for (
unsigned i = 0;
i < tmpTracksters.size(); ++
i) {
208 if (!filter_on_pids(tmpTracksters[
i]) and tmpTracksters[
i].sigmasPCA()[0] < max_longitudinal_sigmaPCA_) {
209 selectedTrackstersIds.push_back(
i);
213 result.reserve(selectedTrackstersIds.size());
215 for (
unsigned i = 0;
i < selectedTrackstersIds.size(); ++
i) {
216 const auto &
t = tmpTracksters[selectedTrackstersIds[
i]];
217 for (
auto const lcId :
t.vertices()) {
218 layer_cluster_usage[lcId]++;
220 LogDebug(
"HGCPatternRecoByCA") <<
"LayerID: " << lcId <<
" count: " << (
int)layer_cluster_usage[lcId]
223 if (isRegionalIter) {
224 seedToTracksterAssociation[
t.seedIndex()].push_back(
i);
229 for (
auto &trackster :
result) {
230 assert(trackster.vertices().size() <= trackster.vertex_multiplicity().size());
231 for (
size_t i = 0;
i < trackster.vertices().size(); ++
i) {
232 trackster.vertex_multiplicity()[
i] = layer_cluster_usage[trackster.vertices(
i)];
234 LogDebug(
"HGCPatternRecoByCA") <<
"LayerID: " << trackster.vertices(
i)
235 <<
" count: " << (
int)trackster.vertex_multiplicity(
i) << std::endl;
239 if (oneTracksterPerTrackSeed_) {
240 std::vector<Trackster>
tmp;
241 mergeTrackstersTRK(
result,
input.layerClusters,
tmp, seedToTracksterAssociation);
247 input.layerClustersTime,
248 rhtools_.getPositionLayer(rhtools_.lastLayerEE(isHFnose), isHFnose).z(),
256 if (promoteEmptyRegionToTrackster_) {
257 emptyTrackstersFromSeedsTRK(
result, seedToTracksterAssociation,
input.regions[0].collectionID);
261 for (
auto &trackster :
result) {
262 LogDebug(
"HGCPatternRecoByCA") <<
"Trackster characteristics: " << std::endl;
263 LogDebug(
"HGCPatternRecoByCA") <<
"Size: " << trackster.vertices().size() << std::endl;
265 for (
auto const &
p : trackster.id_probabilities()) {
273 template <
typename TILES>
275 const std::vector<Trackster> &
input,
277 std::vector<Trackster> &
output,
278 std::unordered_map<
int, std::vector<int>> &seedToTracksterAssociation)
const {
280 for (
auto &thisSeed : seedToTracksterAssociation) {
281 auto &tracksters = thisSeed.second;
282 if (!tracksters.empty()) {
283 auto numberOfTrackstersInSeed = tracksters.size();
285 auto &outTrackster =
output.back();
286 tracksters[0] =
output.size() - 1;
287 auto updated_size = outTrackster.vertices().size();
288 for (
unsigned int j = 1;
j < numberOfTrackstersInSeed; ++
j) {
289 auto &thisTrackster =
input[tracksters[
j]];
290 updated_size += thisTrackster.vertices().size();
292 LogDebug(
"HGCPatternRecoByCA") <<
"Updated size: " << updated_size << std::endl;
294 outTrackster.vertices().reserve(updated_size);
295 outTrackster.vertex_multiplicity().reserve(updated_size);
296 std::copy(std::begin(thisTrackster.vertices()),
297 std::end(thisTrackster.vertices()),
298 std::back_inserter(outTrackster.vertices()));
299 std::copy(std::begin(thisTrackster.vertex_multiplicity()),
300 std::end(thisTrackster.vertex_multiplicity()),
301 std::back_inserter(outTrackster.vertex_multiplicity()));
303 tracksters.resize(1);
306 auto &orig_vtx = outTrackster.vertices();
307 auto vtx_sorted{orig_vtx};
308 std::sort(std::begin(vtx_sorted), std::end(vtx_sorted));
309 for (
unsigned int iLC = 1; iLC < vtx_sorted.size(); ++iLC) {
310 if (vtx_sorted[iLC] == vtx_sorted[iLC - 1]) {
312 const auto lcIdx = vtx_sorted[iLC];
313 const auto firstEl =
std::find(orig_vtx.begin(), orig_vtx.end(), lcIdx);
314 const auto firstPos =
std::distance(std::begin(orig_vtx), firstEl);
316 while (iDup != orig_vtx.end()) {
317 orig_vtx.erase(iDup);
318 outTrackster.vertex_multiplicity().erase(outTrackster.vertex_multiplicity().begin() +
320 outTrackster.vertex_multiplicity()[firstPos] -= 1;
330 template <
typename TILES>
332 std::vector<Trackster> &tracksters,
333 std::unordered_map<
int, std::vector<int>> &seedToTracksterAssociation,
335 for (
auto &thisSeed : seedToTracksterAssociation) {
336 if (thisSeed.second.empty()) {
338 t.setRegressedEnergy(0.
f);
339 t.zeroProbabilities();
341 t.setSeed(collectionID, thisSeed.first);
342 tracksters.emplace_back(
t);
343 thisSeed.second.emplace_back(tracksters.size() - 1);
348 template <
typename TILES>
350 const tensorflow::Session *eidSession,
351 std::vector<Trackster> &tracksters) {
375 std::vector<int> tracksterIndices;
376 for (
int i = 0;
i < (
int)tracksters.size();
i++) {
380 float sumClusterEnergy = 0.;
384 if (sumClusterEnergy >= eidMinClusterEnergy_) {
386 tracksters[
i].setRegressedEnergy(0.
f);
387 tracksters[
i].zeroProbabilities();
388 tracksterIndices.push_back(
i);
401 tensorflow::TensorShape
shape({
batchSize, eidNLayers_, eidNClusters_, eidNFeatures_});
402 tensorflow::Tensor
input(tensorflow::DT_FLOAT,
shape);
405 std::vector<tensorflow::Tensor>
outputs;
407 if (!eidOutputNameEnergy_.empty()) {
410 if (!eidOutputNameId_.empty()) {
416 const Trackster &trackster = tracksters[tracksterIndices[
i]];
421 std::vector<int> clusterIndices(trackster.
vertices().size());
423 clusterIndices[
k] =
k;
425 sort(clusterIndices.begin(), clusterIndices.end(), [&
layerClusters, &trackster](
const int &
a,
const int &
b) {
430 std::vector<int> seenClusters(eidNLayers_);
433 for (
const int &
k : clusterIndices) {
437 if (
j < eidNLayers_ && seenClusters[
j] < eidNClusters_) {
452 for (
int j = 0;
j < eidNLayers_;
j++) {
453 for (
int k = seenClusters[
j];
k < eidNClusters_;
k++) {
455 for (
int l = 0;
l < eidNFeatures_;
l++) {
466 if (!eidOutputNameEnergy_.empty()) {
470 for (
const int &
i : tracksterIndices) {
471 tracksters[
i].setRegressedEnergy(*(
energy++));
476 if (!eidOutputNameId_.empty()) {
478 int probsIdx = eidOutputNameEnergy_.empty() ? 0 : 1;
479 float *probs =
outputs[probsIdx].flat<
float>().
data();
481 for (
const int &
i : tracksterIndices) {
482 tracksters[
i].setProbabilities(probs);
483 probs += tracksters[
i].id_probabilities().size();
488 template <
typename TILES>
490 iDesc.
add<
int>(
"algo_verbosity", 0);
491 iDesc.
add<
bool>(
"oneTracksterPerTrackSeed",
false);
492 iDesc.
add<
bool>(
"promoteEmptyRegionToTrackster",
false);
493 iDesc.
add<
bool>(
"out_in_dfs",
true);
494 iDesc.
add<
int>(
"max_out_in_hops", 10);
495 iDesc.
add<
double>(
"min_cos_theta", 0.915);
496 iDesc.
add<
double>(
"min_cos_pointing", -1.);
497 iDesc.
add<
double>(
"root_doublet_max_distance_from_seed_squared", 9999);
498 iDesc.
add<
double>(
"etaLimitIncreaseWindow", 2.1);
499 iDesc.
add<
int>(
"skip_layers", 0);
500 iDesc.
add<
int>(
"max_missing_layers_in_trackster", 9999);
501 iDesc.
add<
int>(
"shower_start_max_layer", 9999)->setComment(
"make default such that no filtering is applied");
502 iDesc.
add<
int>(
"min_layers_per_trackster", 10);
503 iDesc.
add<std::vector<int>>(
"filter_on_categories", {0});
504 iDesc.
add<
double>(
"pid_threshold", 0.)->
setComment(
"make default such that no filtering is applied");
505 iDesc.
add<
double>(
"energy_em_over_total_threshold", -1.)
506 ->
setComment(
"make default such that no filtering is applied");
507 iDesc.
add<
double>(
"max_longitudinal_sigmaPCA", 9999);
508 iDesc.
add<
double>(
"max_delta_time", 3.)->
setComment(
"nsigma");
510 iDesc.
add<
std::string>(
"eid_output_name_energy",
"output/regressed_energy");
511 iDesc.
add<
std::string>(
"eid_output_name_id",
"output/id_probabilities");
512 iDesc.
add<
double>(
"eid_min_cluster_energy", 1.);
513 iDesc.
add<
int>(
"eid_n_layers", 50);
514 iDesc.
add<
int>(
"eid_n_clusters", 10);
515 iDesc.
add<
bool>(
"computeLocalTime",
false);
516 iDesc.
add<std::vector<double>>(
"siblings_maxRSquared", {6
e-4, 6
e-4, 6
e-4});
void setComment(std::string const &value)
std::vector< NamedTensor > NamedTensorList
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
void mergeTrackstersTRK(const std::vector< Trackster > &, const std::vector< reco::CaloCluster > &, std::vector< Trackster > &, std::unordered_map< int, std::vector< int >> &seedToTracksterAssociation) const
void emptyTrackstersFromSeedsTRK(std::vector< Trackster > &tracksters, std::unordered_map< int, std::vector< int >> &seedToTracksterAssociation, const edm::ProductID &collectionID) const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
~PatternRecognitionbyCA() override
double phi() const
azimuthal angle of cluster centroid
void makeTracksters(const typename PatternRecognitionAlgoBaseT< TILES >::Inputs &input, std::vector< Trackster > &result, std::unordered_map< int, std::vector< int >> &seedToTracksterAssociation) override
static std::string const input
std::vector< float > features(const reco::PreId &ecal, const reco::PreId &hcal, double rho, const reco::BeamSpot &spot, noZS::EcalClusterLazyTools &ecalTools)
def unique(seq, keepstr=True)
void run(Session *session, const NamedTensorList &inputs, const std::vector< std::string > &outputNames, std::vector< Tensor > *outputs, const thread::ThreadPoolOptions &threadPoolOptions)
Abs< T >::type abs(const T &t)
ParameterDescriptionBase * add(U const &iLabel, T const &value)
double energy() const
cluster energy
std::vector< unsigned int > & vertices()
void energyRegressionAndID(const std::vector< reco::CaloCluster > &layerClusters, const tensorflow::Session *, std::vector< Trackster > &result)
PatternRecognitionbyCA(const edm::ParameterSet &conf, edm::ConsumesCollector iC)
std::vector< float > & vertex_multiplicity()
auto const & foundNtuplets
char data[epos_bytes_allocation]
static void fillPSetDescription(edm::ParameterSetDescription &iDesc)
static constexpr int nPhiBins
void assignPCAtoTracksters(std::vector< Trackster > &, const std::vector< reco::CaloCluster > &, const edm::ValueMap< std::pair< float, float >> &, double, bool computeLocalTime=false, bool energyWeight=true)
double eta() const
pseudorapidity of cluster centroid