Process the event and put new data products into it.
42 {
43 const auto& ecal_hits{
event.getCollection<
ldmx::EcalHit>(rec_hit_coll_name_,
44 rec_hit_pass_name_)};
45 if (ecal_hits.size() == 0) {
46
47 ldmx_log(fatal) << "No ECal hits found... exiting";
48 return;
49 }
50
51 if (clue_) {
52 ldmx_log(info) << "Using CLUE clustering algorithm";
53 CLUE cf;
54 cf.cluster(ecal_hits, dc_, rhoc_, deltac_, deltao_, nbr_of_layers_,
55 reclustering_);
56 ldmx_log(debug) << "CLUE algorithm finished clustering";
57 std::vector<IntermediateCluster> interm_cluster = cf.getClusters();
58 std::vector<IntermediateCluster> f_interm_cluster =
59 cf.getFirstLayerCentroids();
60 ldmx_log(debug) << "Got " << interm_cluster.size() << " clusters with "
61 << f_interm_cluster.size() << " first layer centroids";
62
63 auto n_loops = cf.getNLoops();
64 ldmx_log(debug) << "Number of clustering loops: " << n_loops;
65
66 if (reclustering_) {
67 ldmx_log(debug) << "Reclustererd initial number of clusters: "
68 << cf.getInitialClusterNbr()
69 << ", final number of clusters: "
70 << interm_cluster.size();
71 }
72
73 std::vector<ldmx::EcalCluster> ecal_clusters;
74 ldmx_log(debug) << "Filling " << interm_cluster.size()
75 << " clusters into ecal_clusters";
76 for (size_t cluster_indx = 0; cluster_indx < interm_cluster.size();
77 cluster_indx++) {
79
80 cluster.
setEnergy(interm_cluster[cluster_indx].energy());
82 interm_cluster[cluster_indx].centroidY(),
83 interm_cluster[cluster_indx].centroidZ());
84 cluster.setFirstLayerCentroidXYZ(
85 f_interm_cluster[cluster_indx].centroidX(),
86 f_interm_cluster[cluster_indx].centroidY(),
87 f_interm_cluster[cluster_indx].centroidZ());
88 cluster.
setNHits(interm_cluster[cluster_indx].hits().size());
89 cluster.
addHits(interm_cluster[cluster_indx].hits());
90 cluster.addFirstLayerHits(f_interm_cluster[cluster_indx].hits());
91
92 float cl_x(0), cl_y(0), cl_z(0), cl_xx(0), cl_yy(0), cl_zz(0);
93 float cl_w = 1;
94 float sumw = 0;
95
96 for (auto hit : interm_cluster[cluster_indx].hits()) {
97 if (hit->getEnergy() < min_hit_energy_) continue;
98 cl_w = log(hit->getEnergy()) - log(min_hit_energy_);
99 cl_x += cl_w * hit->getXPos();
100 cl_y += cl_w * hit->getYPos();
101 cl_z += cl_w * hit->getZPos();
102 cl_xx += cl_w * hit->getXPos() * hit->getXPos();
103 cl_yy += cl_w * hit->getYPos() * hit->getYPos();
104 cl_zz += cl_w * hit->getZPos() * hit->getZPos();
105 sumw += cl_w;
106 }
107
108 cl_x /= sumw;
109 cl_y /= sumw;
110 cl_z /= sumw;
111 cl_xx /= sumw;
112 cl_yy /= sumw;
113 cl_zz /= sumw;
114 cl_xx = sqrt(cl_xx - cl_x * cl_x);
115 cl_yy = sqrt(cl_yy - cl_y * cl_y);
116 cl_zz = sqrt(cl_zz - cl_z * cl_z);
117
118 cluster.setRMSXYZ(cl_xx, cl_yy, cl_zz);
119 cluster.setLayer(interm_cluster[cluster_indx].layer());
120 ldmx_log(trace) << "Cluster " << cluster_indx
121 << " energy: " << cluster.getEnergy()
122 << ", nHits: " << cluster.getNHits() << ", centroid: ("
127 << cluster.
getRMSZ() <<
")" <<
", first layer centroid: ("
128 << cluster.getFirstLayerCentroidX() << ", "
129 << cluster.getFirstLayerCentroidY() << ", "
130 << cluster.getFirstLayerCentroidZ() << ")";
131
132 ecal_clusters.push_back(cluster);
133 }
134 ldmx_log(debug) << "Filled " << ecal_clusters.size()
135 << " clusters into ecal_clusters";
136 event.add(cluster_coll_name_, ecal_clusters);
137 } else {
138 ldmx_log(info) <<
"Using simple clustering algorithm " <<
algo_name_;
140
142
143 if (hit.getEnergy() == 0) {
144 continue;
145 }
147 }
148
149 cf.
cluster(seed_threshold_, cutoff_);
151 std::map<int, double> c_weights = cf.
getWeights();
152
158
159 std::map<int, double>::iterator it = c_weights.begin();
160 for (it = c_weights.begin(); it != c_weights.end(); it++) {
161 algo_result.
setWeight(it->first, it->second / 100);
162 }
163
164 std::vector<ldmx::EcalCluster> ecal_clusters;
165 for (size_t cluster_indx = 0; cluster_indx < interm_cluster.size();
166 cluster_indx++) {
168
169 cluster.
setEnergy(interm_cluster[cluster_indx].energy());
171 interm_cluster[cluster_indx].centroidY(),
172 interm_cluster[cluster_indx].centroidZ());
173 cluster.
setNHits(interm_cluster[cluster_indx].hits().size());
174 cluster.
addHits(interm_cluster[cluster_indx].hits());
175 ecal_clusters.push_back(cluster);
176 }
177
178 event.add(cluster_coll_name_, ecal_clusters);
179 event.add(algo_coll_name_, algo_result);
180 }
181}
double getRMSX() const
rms in x
void setNHits(int nHits)
Sets total number of hits in the cluster.
double getCentroidZ() const
centroid z-location
double getCentroidX() const
centroid x-location
void setCentroidXYZ(double centroid_x, double centroid_y, double centroid_z)
Sets the three coordinates of the cluster centroid.
double getRMSZ() const
rms in z
double getRMSY() const
rms in y
void setEnergy(double energy)
Sets total energy for the cluster.
double getCentroidY() const
centroid y-location
Contains details about the clustering algorithm.
void setAlgoVar(int element, double value)
Set an algorithm variable.
void set(const TString &name, int nvar)
Set name and number of variables of cluster algo.
void setWeight(int nClusters, double weight)
Set a weight when number of clusters reached.
Stores cluster information from the ECal.
void addHits(const std::vector< const ldmx::EcalHit * > &hits)
Take in the hits that make up the cluster.
Stores reconstructed hit information from the ECAL.
A templated agglomerative clustering algorithm.
void cluster(double seed_threshold, double cutoff)
Run the clustering algorithm.
std::vector< ClusterType > getClusters() const
Get the final clusters after filtering by seed threshold.
std::map< int, double > getWeights() const
Get the transition weights (cluster count -> minimum weight at that step).
int getNSeeds() const
Get the number of seed clusters found.
void add(const HitType &hit)
Add a hit to be clustered using its stored position.