249 {
250 ldmx_log(trace) << "--- CLUSTERING ---";
251 ldmx_log(trace) << "Number of densities: " << densities.size()
252 << "; connecting_layers: " << connecting_layers
253 << "; layer_index: " << layer_index;
254 if (!connecting_layers && nbr_of_layers_ > 1) {
255
256 rhoc_ = layer_rho_c_[layer_index];
257 ldmx_log(trace) << "Setting rho_c on layer " << layer_index << " to "
258 << rhoc_;
259 if (layer_index * 2 - 1 < radius_.size()) {
260 deltac_ = radius_[layer_index * 2 - 1];
261 ldmx_log(trace) << "Setting delta_c on layer " << layer_index << " to "
262 << deltac_;
263 }
264 } else if (connecting_layers) {
265
266
267
268
269
270
271
272 }
273
274 bool energy_overload = false;
275 double max_energy = 10000.;
276 clustering_loops_ = 0;
277 double delta_c_mod = deltac_;
278 double centroid_radius = 10.;
279
280
281 std::vector<std::shared_ptr<Density>>& layer_seeds = seeds_[layer_index];
282
283
284 std::vector<std::vector<const ldmx::EcalHit*>> clusters;
285
286 std::vector<bool> merged_densities;
287 merged_densities.resize(densities.size());
288
289 std::vector<double> cluster_energies;
290 do {
291
292 if (energy_overload) {
293
294 delta_c_mod = delta_c_mod / 1.1;
295 ldmx_log(trace) << "Energy overload, new delta_cmod: " << delta_c_mod;
296 energy_overload = false;
297 }
298
299 clustering_loops_++;
300 ldmx_log(trace) << "Clustering loop " << clustering_loops_;
301
302
303 int k = 0;
304
305 layer_seeds.clear();
306 layer_seeds.reserve(densities.size());
307 clusters.clear();
308 clusters.reserve(densities.size());
309 cluster_energies.clear();
310 cluster_energies.reserve(densities.size());
311
312 std::stack<int> cluster_stack;
313
314 std::vector<std::vector<int>> followers;
315 followers.resize(densities.size());
316
317
318 for (auto& density : densities) {
319
320 ldmx_log(trace) << " Index: " << density->index_
321 << "; x: " << density->x_ << "; y: " << density->y_
322 << "; Energy: " << density->total_energy_
323 << " Parent ID: " << density->follower_of_
324 << "; Delta: " << density->delta_;
325
326 bool is_seed;
327 if (delta_c_mod != deltac_ && merged_densities[density->cluster_id_] &&
328 dist(density->x_, density->y_, event_centroid_.centroid().x(),
329 event_centroid_.centroid().y()) < centroid_radius) {
330
331
332
333 is_seed =
334 density->total_energy_ > rhoc_ && density->delta_ > delta_c_mod;
335 } else {
336 is_seed = density->total_energy_ > rhoc_ && density->delta_ > deltac_;
337 }
338
339 bool is_outlier =
340 (density->total_energy_ < rhoc_) && (density->delta_ > deltao_);
341 density->cluster_id_ = -1;
342 if (is_seed) {
343 ldmx_log(trace) << " This is a Seed";
344 ldmx_log(trace) << " Distance to centroid: "
345 <<
dist(density->x_, density->y_,
346 event_centroid_.centroid().x(),
347 event_centroid_.centroid().y())
348 << "; with delta " << density->delta_;
349 ldmx_log(trace) << " Setting cluster ID to " << k;
350 density->cluster_id_ = k;
351 k++;
352
353 cluster_stack.push(density->index_);
354 clusters.push_back(density->hits_);
355 cluster_energies.push_back(density->total_energy_);
356 layer_seeds.push_back(density);
357 } else if (!is_outlier) {
358 ldmx_log(trace) << " This is a Follower";
359 int& parent_index = density->follower_of_;
360 if (parent_index != -1)
361 followers[parent_index].push_back(density->index_);
362 } else {
363 ldmx_log(trace) << " This is an Outlier";
364 }
365 }
366
367 merged_densities.clear();
368 merged_densities.resize(densities.size());
369
370
371 while (cluster_stack.size() > 0) {
372 auto& d = densities[cluster_stack.top()];
373 cluster_stack.pop();
374 auto& cid = d->cluster_id_;
375
376 for (const auto follower_index : followers[d->index_]) {
377 auto& follower = densities[follower_index];
378
379 follower->cluster_id_ = cid;
380 cluster_energies[cid] += follower->total_energy_;
381
382 if (reclustering_ && cluster_energies[cid] > max_energy &&
383 delta_c_mod > 0.5 && clustering_loops_ < 100) {
384
385
386 merged_densities[cid] = true;
387 if (!energy_overload && clustering_loops_ == 99) {
388 ldmx_log(warn) << "Merging clusters, max cluster loops hit";
389 }
390 energy_overload = true;
391 if (clustering_loops_ != 1) {
392 goto endwhile;
393
394 }
395 }
396
397 clusters[cid].insert(std::end(clusters[cid]),
398 std::begin(follower->hits_),
399 std::end(follower->hits_));
400
401
402 cluster_stack.push(follower_index);
403 }
404 }
405
406
407 if (clustering_loops_ == 1 && energy_overload)
408 initial_cluster_nbr_ = clusters.size();
409 endwhile:;
410 } while (energy_overload);
411
412 if (!connecting_layers && nbr_of_layers_ > 1) {
413
414
415 for (auto& seed : layer_seeds) {
416 seed->delta_ = std::numeric_limits<float>::max();
417 seed->hits_ = clusters[seed->cluster_id_];
418 seed->total_energy_ = cluster_energies[seed->cluster_id_];
419 seed->index_ = seed_index_;
420 seed_index_++;
421 }
422
423 std::sort(layer_seeds.begin(), layer_seeds.end(),
424 [](const std::shared_ptr<Density>& a,
425 const std::shared_ptr<Density>& b) {
426 return a->total_energy_ > b->total_energy_;
427 });
428 }
429 return clusters;
430}
T dist(T x1, T y1, T x2, T y2)
Euclidean distance between two points.