LDMX Software
ecal::CLUE Class Reference

Classes

struct  Density
 

Public Member Functions

template<typename T >
dist (T x1, T y1, T x2, T y2)
 Euclidean distance between two points.
 
template<typename T >
dist (T x1, T y1, T z1, T x2, T y2, T z2)
 
std::vector< std::vector< const ldmx::EcalHit * > > createLayers (const std::vector< const ldmx::EcalHit * > &hits)
 
float roundToDecimal (float x, int num_decimal_precision_digits)
 
std::vector< std::shared_ptr< Density > > setup (const std::vector< const ldmx::EcalHit * > &hits)
 
std::vector< std::vector< const ldmx::EcalHit * > > clustering (std::vector< std::shared_ptr< Density > > &densities, bool connectingLayers, int layerTag=0)
 
std::vector< std::shared_ptr< Density > > setupForClue3D ()
 
void convertToIntermediateClusters (std::vector< std::vector< const ldmx::EcalHit * > > &clusters)
 
void cluster (const std::vector< ldmx::EcalHit > &hits, double dc, double rc, double deltac, double deltao, int nbrOfLayers, bool reclustering)
 
std::vector< double > getCentroidDistances () const
 
int getNLoops () const
 
int getInitialClusterNbr () const
 
std::vector< IntermediateClustergetClusters () const
 
std::vector< IntermediateClustergetFirstLayerCentroids () const
 

Private Member Functions

 enableLogging ("CLUE")
 

Private Attributes

int clustering_loops_
 
bool reclustering_
 
double dc_
 
double rhoc_
 
double deltac_
 
double deltao_
 
double dm_
 
int max_layers_ {32}
 
int nbr_of_layers_
 
std::vector< double > layer_rho_c_
 
std::vector< double > layer_delta_c_
 
std::vector< double > radius_
 
std::vector< double > centroid_distances_
 
IntermediateCluster event_centroid_
 
std::vector< IntermediateClusterfirst_layer_centroids_
 
int seed_index_ {0}
 
std::vector< std::vector< std::shared_ptr< Density > > > seeds_
 
int initial_cluster_nbr_ {-1}
 
std::vector< IntermediateClusterfinal_clusters_
 
std::vector< std::pair< double, double > > layer_centroid_separations_
 

Detailed Description

Definition at line 27 of file CLUE.h.

Member Function Documentation

◆ cluster()

void ecal::CLUE::cluster ( const std::vector< ldmx::EcalHit > & hits,
double dc,
double rc,
double deltac,
double deltao,
int nbrOfLayers,
bool reclustering )

Definition at line 558 of file CLUE.cxx.

560 {
561 ldmx_log(info) << "Starting CLUE clustering with parameters:" << "dc " << dc
562 << ", rc " << rc << ", delta_c " << delta_c << ", delta_o "
563 << delta_o << ", nbr_of_layers " << nbr_of_layers
564 << ", reclustering " << reclustering;
565 // cutoff distance for local density
566 dc_ = dc;
567 // min density to promote as seed/max density to demote as outlier
568 rhoc_ = rc;
569 // min separation distance for seeds
570 deltac_ = delta_c;
571 // min separation distance for outliers
572 deltao_ = delta_o;
573 // max distance to parent for both seeds and followers
574 dm_ = std::max(delta_c, delta_o);
575
576 // Recluster merged clusters or not
577 reclustering_ = reclustering;
578 nbr_of_layers_ = nbr_of_layers;
579
580 if (nbr_of_layers_ < 1) {
581 // anything below 1 => include all layers
582 nbr_of_layers_ = max_layers_;
583 } else if (nbr_of_layers_ > max_layers_) {
584 ldmx_log(warn) << "nbr_of_layers_ " << nbr_of_layers_
585 << " exceeds max layers " << max_layers_
586 << ", setting to max layers";
587 nbr_of_layers_ = max_layers_;
588 }
589
590 // first copy *addresses* so we are only ever passing around pointers
591 std::vector<const ldmx::EcalHit*> hits;
592 hits.reserve(unsorted_hits.size());
593 ldmx_log(debug) << "Clustering " << unsorted_hits.size() << " hits";
594 for (const auto& unsorted_hit : unsorted_hits) {
595 hits.push_back(&unsorted_hit);
596 }
597 // sort hits by Z position
598 ldmx_log(debug) << "Sorting hits by Z position";
599 std::sort(hits.begin(), hits.end(),
600 [](const ldmx::EcalHit* a, const ldmx::EcalHit* b) {
601 return a->getZPos() < b->getZPos();
602 });
603
604 seeds_.resize(nbr_of_layers);
605
606 if (nbr_of_layers_ > 1) {
607 ldmx_log(debug) << "Creating layers";
608 // returns a vector of layers, with each layer having a vector of hits
609 const auto layers = createLayers(hits);
610 ldmx_log(debug) << "Doing layer-wise clustering on " << layers.size()
611 << " layers";
612 for (int i = 0; i < layers.size(); i++) {
613 ldmx_log(trace) << "--- LAYER " << i + 1 << " ---";
614 auto densities = setup(layers[i]);
615 auto clusters = clustering(densities, false, i);
616 convertToIntermediateClusters(clusters);
617 // clustering without 3D
618 }
619 // Below for CLUE3D, comment for just layer clustering
620 // This does not work properly yet
621 // auto densities = setupForClue3D();
622 // auto clusters = clustering(densities, true);
623 // convertToIntermediateClusters(clusters);
624 } else {
625 ldmx_log(debug) << "Only one layer, doing 2D clustering";
626 auto densities = setup(hits);
627 auto clusters = clustering(densities, false);
628 convertToIntermediateClusters(clusters);
629 }
630}
Stores reconstructed hit information from the ECAL.
Definition EcalHit.h:19

◆ clustering()

std::vector< std::vector< const ldmx::EcalHit * > > ecal::CLUE::clustering ( std::vector< std::shared_ptr< Density > > & densities,
bool connectingLayers,
int layerTag = 0 )

Definition at line 247 of file CLUE.cxx.

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 // if layerwise clustering override rhoc_ and deltac_ with per-layer values
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 // if doing 3D clustering, override deltac_ and rhoc_ with other values
266 // NOT IMPLEMENTED YET
267 // if currently doing 3D clustering
268 // IF 3D clustering
269 // deltao_ = 200.;
270 // deltac_ = 100.;
271 // rhoc_ = 1000.;
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 // stores seeds of this layer
281 std::vector<std::shared_ptr<Density>>& layer_seeds = seeds_[layer_index];
282
283 // stores hits in cluster
284 std::vector<std::vector<const ldmx::EcalHit*>> clusters;
285 // keeps track of which densities have merged; only used if reclustering
286 std::vector<bool> merged_densities; // index_= cluster id
287 merged_densities.resize(densities.size());
288 // keeps track of cluster energies
289 std::vector<double> cluster_energies;
290 do {
291 // while no cluster has merged
292 if (energy_overload) {
293 // makes delta_c smaller if clusters have merged
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 // cluster index
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 // stores followers of densities at corr index_
314 std::vector<std::vector<int>> followers;
315 followers.resize(densities.size());
316
317 // Mark as seed, follower, or outlier
318 for (auto& density : densities) {
319 // funky line to generalize this function for both 2D and 3D case
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 // if energy has been overloaded and this density belongs to cluster
331 // that was overloaded and this density is close enough to event
332 // centroid use modded delta c
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 // get the index of the seed density
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 // Go through all seeds and add followers, then follower's followers, etc.
371 while (cluster_stack.size() > 0) {
372 auto& d = densities[cluster_stack.top()];
373 cluster_stack.pop();
374 auto& cid = d->cluster_id_;
375 // for indices of followers of dp
376 for (const auto follower_index : followers[d->index_]) {
377 auto& follower = densities[follower_index];
378 // Set cluster index of the follower to the cluster index of `d`
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 // If reclustering is on and cluster energy is too high,
385 // delta_c_mod is not too low, and we haven't tried for too long
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; // Don't break on the first loop to save the initial
393 // cluster number
394 }
395 }
396
397 clusters[cid].insert(std::end(clusters[cid]),
398 std::begin(follower->hits_),
399 std::end(follower->hits_));
400 // Add follower to the stack so its followers can also get the correct
401 // cluster index
402 cluster_stack.push(follower_index);
403 }
404 }
405 // for first clusteringloop, we want to save number of clusters before
406 // reclustering
407 if (clustering_loops_ == 1 && energy_overload)
408 initial_cluster_nbr_ = clusters.size();
409 endwhile:;
410 } while (energy_overload);
411 // if we have more than one layer and we are not currently doing CLUE3D
412 if (!connecting_layers && nbr_of_layers_ > 1) {
413 // Overwrite seed densities' properties with cluster properties
414 // Might be cleaner to just create new densities for cluster seeds
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 // Sort seeds in layer based on energy
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} // end of clustering
T dist(T x1, T y1, T x2, T y2)
Euclidean distance between two points.
Definition CLUE.cxx:15

◆ convertToIntermediateClusters()

void ecal::CLUE::convertToIntermediateClusters ( std::vector< std::vector< const ldmx::EcalHit * > > & clusters)

Definition at line 528 of file CLUE.cxx.

529 {
530 // Convert to workingecalclusters to ensure compatibility with
531 // EcalClusterProducer
532 for (const auto& cluster : clusters) {
533 IntermediateCluster intermediate_cluster{},
534 intermediate_cluster_first_layer{};
535
536 for (const auto& hit : cluster) {
537 intermediate_cluster.add(hit);
538 // if hit is in first layer, add to first layer cluster
539 ldmx::EcalID ecal_id(hit->getID());
540 auto layer = ecal_id.layer();
541 intermediate_cluster.setLayer(layer);
542 if (layer == 0) {
543 intermediate_cluster_first_layer.add(hit);
544 intermediate_cluster_first_layer.setLayer(layer);
545 }
546 }
547 final_clusters_.push_back(intermediate_cluster);
548 first_layer_centroids_.push_back(intermediate_cluster_first_layer);
549 auto cent_x = intermediate_cluster.centroid().x();
550 auto cent_y = intermediate_cluster.centroid().y();
551 auto event_cent_x = event_centroid_.centroid().x();
552 auto event_cent_y = event_centroid_.centroid().y();
553 const auto& distance = dist(cent_x, cent_y, event_cent_x, event_cent_y);
554 centroid_distances_.push_back(distance);
555 }
556}
Extension of DetectorID providing access to ECal layers and cell numbers in a hex grid.
Definition EcalID.h:20

◆ createLayers()

std::vector< std::vector< const ldmx::EcalHit * > > ecal::CLUE::createLayers ( const std::vector< const ldmx::EcalHit * > & hits)

Definition at line 97 of file CLUE.cxx.

98 {
99 ldmx_log(trace) << "--- LAYER CREATION ---";
100 ldmx_log(trace) << "Number of layers: " << nbr_of_layers_;
101
102 // vector of layers, each layer having a vector of hits
103 // initialize with nbr_of_layers_ empty vectors
104 std::vector<std::vector<const ldmx::EcalHit*>> hits_per_layer(nbr_of_layers_);
105
106 // Track highest energy per layer for proper rho_c calculation
107 std::vector<double> layer_max_energies(nbr_of_layers_, 0.0);
108
109 // Clear any existing layer_rho_c_ values
110 layer_rho_c_.clear();
111
112 // This is rather ad-hoc, but this way it's still tunable via rhoc_ parameter
113 double rhoc_factor = rhoc_ / 250.;
114
115 for (const auto& hit : hits) {
116 ldmx::EcalID ecal_id(hit->getID());
117 // layer number from EcalID, starting at 0
118 int layer = ecal_id.layer();
119
120 if (layer >= nbr_of_layers_) {
121 ldmx_log(trace) << "Skipping hit in layer " << layer
122 << " (beyond nbr_of_layers_ = " << nbr_of_layers_ << ")";
123 continue;
124 }
125
126 // Track the highest energy in this specific layer
127 if (hit->getEnergy() > layer_max_energies[layer]) {
128 layer_max_energies[layer] = hit->getEnergy();
129 }
130
131 ldmx_log(trace) << " Adding hit with energy " << hit->getEnergy()
132 << " to layer " << layer;
133 hits_per_layer[layer].push_back(hit);
134 } // end of loop over hits
135
136 // Calculate rhoc for each layer based on that layer's maximum energy
137 for (int i = 0; i < nbr_of_layers_; i++) {
138 double rho_c = layer_max_energies[i] / rhoc_factor;
139 layer_rho_c_.push_back(rho_c);
140 ldmx_log(trace) << "Layer " << i << " has " << hits_per_layer[i].size()
141 << " hits, max energy " << layer_max_energies[i]
142 << ", rho_c " << rho_c;
143 }
144
145 ldmx_log(trace) << "Created " << hits_per_layer.size() << " layers";
146 return hits_per_layer;
147} // end of createLayers

◆ dist() [1/2]

template<typename T >
T ecal::CLUE::dist ( T x1,
T y1,
T x2,
T y2 )

Euclidean distance between two points.

Euclidean distance between two points Using x*x and std::sqrt since they are more specialized and faster than std::pow.

Template Parameters
Tfloating point type

Definition at line 15 of file CLUE.cxx.

15 {
16 auto delta_x = x1 - x2;
17 auto delta_y = y1 - y2;
18 auto r_square = delta_x * delta_x + delta_y * delta_y;
19 return std::sqrt(r_square);
20}

◆ dist() [2/2]

template<typename T >
T ecal::CLUE::dist ( T x1,
T y1,
T z1,
T x2,
T y2,
T z2 )

Definition at line 24 of file CLUE.cxx.

24 {
25 auto delta_x = x1 - x2;
26 auto delta_y = y1 - y2;
27 auto delta_z = z1 - z2;
28 auto r_square = delta_x * delta_x + delta_y * delta_y + delta_z * delta_z;
29 return std::sqrt(r_square);
30}

◆ getCentroidDistances()

std::vector< double > ecal::CLUE::getCentroidDistances ( ) const
inline

Definition at line 98 of file CLUE.h.

98 {
99 return centroid_distances_;
100 }

◆ getClusters()

std::vector< IntermediateCluster > ecal::CLUE::getClusters ( ) const
inline

Definition at line 106 of file CLUE.h.

106 {
107 return final_clusters_;
108 }

◆ getFirstLayerCentroids()

std::vector< IntermediateCluster > ecal::CLUE::getFirstLayerCentroids ( ) const
inline

Definition at line 112 of file CLUE.h.

112 {
113 return first_layer_centroids_;
114 }

◆ getInitialClusterNbr()

int ecal::CLUE::getInitialClusterNbr ( ) const
inline

Definition at line 104 of file CLUE.h.

104{ return initial_cluster_nbr_; }

◆ getNLoops()

int ecal::CLUE::getNLoops ( ) const
inline

Definition at line 102 of file CLUE.h.

102{ return clustering_loops_; }

◆ roundToDecimal()

float ecal::CLUE::roundToDecimal ( float x,
int num_decimal_precision_digits )

Definition at line 149 of file CLUE.cxx.

149 {
150 float power_of_10 = std::pow(10, num_decimal_precision_digits);
151 return std::round(x_ * power_of_10) / power_of_10;
152}

◆ setup()

std::vector< std::shared_ptr< CLUE::Density > > ecal::CLUE::setup ( const std::vector< const ldmx::EcalHit * > & hits)

Definition at line 154 of file CLUE.cxx.

155 {
156 std::vector<std::shared_ptr<Density>> densities;
157 std::map<std::pair<float, float>, std::shared_ptr<Density>> density_map;
158 event_centroid_ = IntermediateCluster();
159 ldmx_log(trace) << "--- SETUP ---";
160 ldmx_log(trace) << "Building densities";
161 for (const auto& hit : hits) {
162 // collapse z dimension
163 float x = roundToDecimal(hit->getXPos(), 4);
164 float y = roundToDecimal(hit->getYPos(), 4);
165 ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}";
166 std::pair<float, float> coords;
167 if (dc_ != 0 && nbr_of_layers_ > 1) {
168 // if more than one layer, divide hit into densities with side dc
169 double i = std::ceil(std::abs(x) / dc_);
170 double j = std::ceil(std::abs(y) / dc_);
171 if (x < 0) {
172 i = -i;
173 x = (i + 0.5) * dc_;
174 } else {
175 x = (i - 0.5) * dc_;
176 }
177 if (y < 0) {
178 j = -j;
179 y = (j + 0.5) * dc_;
180 } else {
181 y = (j - 0.5) * dc_;
182 }
183 coords = {i, j};
184 ldmx_log(trace) << " Index " << i << ", " << j << "; x: " << x
185 << " y: " << y;
186 } else {
187 // if just one layer, have all densities with the same x,y be in same
188 // density
189 coords = {x, y};
190 }
191 if (density_map.find(coords) == density_map.end()) {
192 density_map.emplace(coords, std::make_shared<CLUE::Density>(x, y));
193 ldmx_log(trace) << " New density created";
194 } else {
195 ldmx_log(trace) << " Found density with x: " << density_map[coords]->x_
196 << " y: " << density_map[coords]->y_;
197 }
198 density_map[coords]->hits_.push_back(hit);
199 density_map[coords]->total_energy_ += hit->getEnergy();
200 density_map[coords]->z_ += hit->getZPos();
201
202 event_centroid_.add(hit);
203 } // end of loop over hits
204
205 densities.reserve(density_map.size());
206 for (const auto& entry : density_map) {
207 densities.push_back(std::move(entry.second));
208 }
209 // sort according to energy
210 std::sort(densities.begin(), densities.end(),
211 [](const std::shared_ptr<CLUE::Density>& a,
212 const std::shared_ptr<CLUE::Density>& b) {
213 return a->total_energy_ > b->total_energy_;
214 });
215
216 ldmx_log(trace) << "Decide parents";
217
218 // decide delta_ and follower_of_
219 for (int i = 0; i < densities.size(); i++) {
220 densities[i]->index_ = i;
221 // avg z position
222 densities[i]->z_ = densities[i]->z_ / densities[i]->hits_.size();
223 ldmx_log(trace) << " Index: " << i << "; x: " << densities[i]->x_
224 << "; y: " << densities[i]->y_
225 << "; Energy: " << densities[i]->total_energy_;
226 // loop through all higher energy densities
227 for (int j = 0; j < i; j++) {
228 float distance_2d = dist(densities[i]->x_, densities[i]->y_,
229 densities[j]->x_, densities[j]->y_);
230 // condition energyJ > energyI but this should be baked in as we sorted
231 // according to energy
232 if ((distance_2d < dm_) && (distance_2d < densities[i]->delta_)) {
233 densities[i]->delta_ = distance_2d;
234 densities[i]->follower_of_ = j;
235 ldmx_log(trace) << " New parent, index " << j
236 << "; delta_2d: " << std::setprecision(4)
237 << distance_2d;
238 }
239 }
240 }
241 return densities;
242} // end of setup

◆ setupForClue3D()

std::vector< std::shared_ptr< CLUE::Density > > ecal::CLUE::setupForClue3D ( )

Definition at line 434 of file CLUE.cxx.

434 {
435 ldmx_log(trace) << "--- LAYER SETUP ---";
436 std::vector<std::shared_ptr<CLUE::Density>> densities;
437 layer_rho_c_.clear();
438 for (int layer = 0; layer < nbr_of_layers_; layer++) {
439 ldmx_log(trace) << " LAYER " << layer << " with " << seeds_[layer].size()
440 << " seeds";
441 auto& seeds_in_current_layer = seeds_[layer];
442 double highest_energy = 0.;
443 for (const auto& current_seed : seeds_in_current_layer) {
444 // for each seed in layer
445 current_seed->layer_ = layer;
446 if (current_seed->total_energy_ > highest_energy)
447 highest_energy = current_seed->total_energy_;
448 ldmx_log(trace) << " Density with index " << current_seed->index_
449 << ", energy: " << current_seed->total_energy_;
450 int depth = 1;
451 // decide delta_ and followerof from seeds in previous and next layer_
452 // do {
453 // depth++;
454 if ((layer - depth >= 0) && (layer - depth < seeds_.size())) {
455 ldmx_log(trace) << " Looking at pre-layer: " << layer - depth;
456 // look at previous layer
457 auto& previous_layer = seeds_[layer - depth];
458 for (const auto& prev_seed : previous_layer) {
459 // for each seed in previous layer
460 auto distance_2d_prev = dist(current_seed->x_, current_seed->y_,
461 prev_seed->x_, prev_seed->y_);
462 auto dz_prev = std::abs(current_seed->z_ - prev_seed->z_);
463 ldmx_log(trace) << " DeltaXY to index " << prev_seed->index_
464 << ": " << std::setprecision(4) << distance_2d_prev;
465 ldmx_log(trace) << " DeltaZ to index " << prev_seed->index_
466 << ": " << std::setprecision(4) << dz_prev;
467 if (prev_seed->total_energy_ > current_seed->total_energy_ &&
468 distance_2d_prev < current_seed->delta_ &&
469 dz_prev < current_seed->z_delta_) {
470 ldmx_log(trace) << " New parent: index " << prev_seed->index_
471 << " on layer " << layer - depth << "; energy "
472 << prev_seed->total_energy_;
473 ldmx_log(trace) << " New 2D distance to prev-layer: "
474 << std::setprecision(4) << distance_2d_prev;
475 ldmx_log(trace)
476 << " New delta Z: " << std::setprecision(4) << dz_prev;
477 current_seed->delta_ = distance_2d_prev;
478 current_seed->z_delta_ = dz_prev;
479 current_seed->follower_of_ = prev_seed->index_;
480 } else if (prev_seed->total_energy_ < current_seed->total_energy_) {
481 ldmx_log(trace) << " Breaking on index " << prev_seed->index_
482 << " with energy " << prev_seed->total_energy_;
483 break;
484 }
485 }
486 }
487 if (layer + depth < nbr_of_layers_ && layer + depth < seeds_.size()) {
488 ldmx_log(trace) << " Looking at post-layer: " << layer + depth;
489 auto& next_layer = seeds_[layer + depth];
490 for (const auto& next_seed : next_layer) {
491 auto distance_2d_next = dist(current_seed->x_, current_seed->y_,
492 next_seed->x_, next_seed->y_);
493 auto dz_next = std::abs(current_seed->z_ - next_seed->z_);
494 ldmx_log(trace) << " DeltaXY to index " << next_seed->index_
495 << ": " << std::setprecision(4) << distance_2d_next;
496 ldmx_log(trace) << " DeltaZ to index_" << next_seed->index_
497 << ": " << std::setprecision(4) << dz_next;
498 if (next_seed->total_energy_ > current_seed->total_energy_ &&
499 distance_2d_next < current_seed->delta_ &&
500 dz_next < current_seed->z_delta_) {
501 ldmx_log(trace) << " New parent: index_" << next_seed->index_
502 << " on layer " << layer + depth << "; energy "
503 << next_seed->total_energy_;
504 ldmx_log(trace) << " New 2D distance to next-layer: "
505 << std::setprecision(4) << distance_2d_next;
506 ldmx_log(trace)
507 << " New delta_Z: " << std::setprecision(4) << dz_next;
508 current_seed->delta_ = distance_2d_next;
509 current_seed->z_delta_ = dz_next;
510 current_seed->follower_of_ = next_seed->index_;
511 } else if (next_seed->total_energy_ < current_seed->total_energy_) {
512 ldmx_log(trace) << " Breaking on index_" << next_seed->index_
513 << " with energy " << next_seed->total_energy_;
514 break;
515 }
516 } // end loop on seeds in next layer
517 } // end of looking at next layer
518 // } while (currentLayer[i]->layerFollowerOf == -1 && (layer - depth >=
519 // 0 || layer + depth < nbr_of_layers_));
520 densities.push_back(current_seed);
521 }
522 // TODO: This 2 needs to be configurable
523 layer_rho_c_.push_back(highest_energy / 2);
524 } // end loop over layers
525 return densities;
526} // end of setupForClue3D

Member Data Documentation

◆ centroid_distances_

std::vector<double> ecal::CLUE::centroid_distances_
private

Definition at line 146 of file CLUE.h.

◆ clustering_loops_

int ecal::CLUE::clustering_loops_
private

Definition at line 117 of file CLUE.h.

◆ dc_

double ecal::CLUE::dc_
private

Definition at line 121 of file CLUE.h.

◆ deltac_

double ecal::CLUE::deltac_
private

Definition at line 123 of file CLUE.h.

◆ deltao_

double ecal::CLUE::deltao_
private

Definition at line 124 of file CLUE.h.

◆ dm_

double ecal::CLUE::dm_
private

Definition at line 125 of file CLUE.h.

◆ event_centroid_

IntermediateCluster ecal::CLUE::event_centroid_
private

Definition at line 147 of file CLUE.h.

◆ final_clusters_

std::vector<IntermediateCluster> ecal::CLUE::final_clusters_
private

Definition at line 155 of file CLUE.h.

◆ first_layer_centroids_

std::vector<IntermediateCluster> ecal::CLUE::first_layer_centroids_
private

Definition at line 149 of file CLUE.h.

◆ initial_cluster_nbr_

int ecal::CLUE::initial_cluster_nbr_ {-1}
private

Definition at line 154 of file CLUE.h.

154{-1};

◆ layer_centroid_separations_

std::vector<std::pair<double, double> > ecal::CLUE::layer_centroid_separations_
private

Definition at line 156 of file CLUE.h.

◆ layer_delta_c_

std::vector<double> ecal::CLUE::layer_delta_c_
private

Definition at line 132 of file CLUE.h.

◆ layer_rho_c_

std::vector<double> ecal::CLUE::layer_rho_c_
private

Definition at line 131 of file CLUE.h.

◆ max_layers_

int ecal::CLUE::max_layers_ {32}
private

Definition at line 128 of file CLUE.h.

128{32};

◆ nbr_of_layers_

int ecal::CLUE::nbr_of_layers_
private

Definition at line 129 of file CLUE.h.

◆ radius_

std::vector<double> ecal::CLUE::radius_
private
Initial value:
{
5.723387467629167, 5.190678018534044, 5.927290663506518,
6.182560329200212, 7.907549398117859, 8.606100542857211,
10.93381822596916, 12.043201938160239, 14.784548371508041,
16.102403056546482, 18.986402399412817, 20.224453740305716,
23.048820910305643, 24.11202594672678, 26.765135236851666,
27.78700483852502, 30.291794353801293, 31.409870873194464,
33.91006482486666, 35.173073672355926, 38.172422630271,
40.880288341493205, 44.696485719120005, 49.23802839743545,
53.789910813378675, 60.87843355562641, 66.32931132415688,
75.78117972604727, 86.04697356716805, 96.90360704034346}

Definition at line 134 of file CLUE.h.

134 {
135 5.723387467629167, 5.190678018534044, 5.927290663506518,
136 6.182560329200212, 7.907549398117859, 8.606100542857211,
137 10.93381822596916, 12.043201938160239, 14.784548371508041,
138 16.102403056546482, 18.986402399412817, 20.224453740305716,
139 23.048820910305643, 24.11202594672678, 26.765135236851666,
140 27.78700483852502, 30.291794353801293, 31.409870873194464,
141 33.91006482486666, 35.173073672355926, 38.172422630271,
142 40.880288341493205, 44.696485719120005, 49.23802839743545,
143 53.789910813378675, 60.87843355562641, 66.32931132415688,
144 75.78117972604727, 86.04697356716805, 96.90360704034346};

◆ reclustering_

bool ecal::CLUE::reclustering_
private

Definition at line 119 of file CLUE.h.

◆ rhoc_

double ecal::CLUE::rhoc_
private

Definition at line 122 of file CLUE.h.

◆ seed_index_

int ecal::CLUE::seed_index_ {0}
private

Definition at line 151 of file CLUE.h.

151{0};

◆ seeds_

std::vector<std::vector<std::shared_ptr<Density> > > ecal::CLUE::seeds_
private

Definition at line 152 of file CLUE.h.


The documentation for this class was generated from the following files: