Process the event and put new data products into it.
86 {
87 const auto& ecal_hits{
event.getCollection<
ldmx::EcalHit>(rec_hit_coll_name_,
88 rec_hit_pass_name_)};
89 if (ecal_hits.size() == 0) {
90
91 ldmx_log(debug) << "No ECal hits found";
92 std::vector<ldmx::EcalCluster> empty;
93 event.add(cluster_coll_name_, empty);
94 return;
95 }
96
97 if (algorithm_ == "clue3d") {
98 CLUE3D algo(clue3d_);
99 auto clusters = algo.cluster(ecal_hits);
100 std::vector<ldmx::EcalCluster> ecal_clusters;
101 for (const auto& c : clusters) {
105 cluster.setDXDZ(c.dxdz_);
106 cluster.setDYDZ(c.dydz_);
107 cluster.setRMSXYZ(c.rms_x_, c.rms_y_, c.rms_z_);
108 cluster.setFirstLayerCentroidXYZ(c.x0_, c.y0_, c.z0_);
109
110 cluster.setLayer(0);
111 std::vector<unsigned int> ids;
112 std::vector<float> hx, hy, hz, he;
113 for (std::size_t i = 0; i < c.hits_.size(); i++) {
114 const auto& h = ecal_hits[c.hits_[i]];
115 ids.push_back(h.getID());
116 hx.push_back(h.getXPos());
117 hy.push_back(h.getYPos());
118 hz.push_back(h.getZPos());
119
120 he.push_back(c.frac_[i] * h.getEnergy());
121 }
124 cluster.setHitValsX(hx);
125 cluster.setHitValsY(hy);
126 cluster.setHitValsZ(hz);
127 cluster.setHitValsE(he);
128 ecal_clusters.push_back(cluster);
129 }
130 event.add(cluster_coll_name_, ecal_clusters);
131 } else if (clue_) {
132 ldmx_log(info) << "Using CLUE clustering algorithm";
133 CLUE cf;
134 cf.cluster(ecal_hits, dc_, rhoc_, deltac_, deltao_, nbr_of_layers_,
135 reclustering_);
136 ldmx_log(debug) << "CLUE algorithm finished clustering";
137 std::vector<IntermediateCluster> interm_cluster = cf.getClusters();
138 std::vector<IntermediateCluster> f_interm_cluster =
139 cf.getFirstLayerCentroids();
140 ldmx_log(debug) << "Got " << interm_cluster.size() << " clusters with "
141 << f_interm_cluster.size() << " first layer centroids";
142
143 auto n_loops = cf.getNLoops();
144 ldmx_log(debug) << "Number of clustering loops: " << n_loops;
145
146 if (reclustering_) {
147 ldmx_log(debug) << "Reclustererd initial number of clusters: "
148 << cf.getInitialClusterNbr()
149 << ", final number of clusters: "
150 << interm_cluster.size();
151 }
152
153 std::vector<ldmx::EcalCluster> ecal_clusters;
154 ldmx_log(debug) << "Filling " << interm_cluster.size()
155 << " clusters into ecal_clusters";
156 for (size_t cluster_indx = 0; cluster_indx < interm_cluster.size();
157 cluster_indx++) {
159
160 cluster.
setEnergy(interm_cluster[cluster_indx].energy());
162 interm_cluster[cluster_indx].centroidY(),
163 interm_cluster[cluster_indx].centroidZ());
164 cluster.setFirstLayerCentroidXYZ(
165 f_interm_cluster[cluster_indx].centroidX(),
166 f_interm_cluster[cluster_indx].centroidY(),
167 f_interm_cluster[cluster_indx].centroidZ());
168 cluster.
setNHits(interm_cluster[cluster_indx].hits().size());
169 cluster.
addHits(interm_cluster[cluster_indx].hits());
170 cluster.addFirstLayerHits(f_interm_cluster[cluster_indx].hits());
171
172 float cl_x(0), cl_y(0), cl_z(0), cl_xx(0), cl_yy(0), cl_zz(0);
173 float cl_w = 1;
174 float sumw = 0;
175
176 for (auto hit : interm_cluster[cluster_indx].hits()) {
177 if (hit->getEnergy() < min_hit_energy_) continue;
178 cl_w = log(hit->getEnergy()) - log(min_hit_energy_);
179 cl_x += cl_w * hit->getXPos();
180 cl_y += cl_w * hit->getYPos();
181 cl_z += cl_w * hit->getZPos();
182 cl_xx += cl_w * hit->getXPos() * hit->getXPos();
183 cl_yy += cl_w * hit->getYPos() * hit->getYPos();
184 cl_zz += cl_w * hit->getZPos() * hit->getZPos();
185 sumw += cl_w;
186 }
187
188 cl_x /= sumw;
189 cl_y /= sumw;
190 cl_z /= sumw;
191 cl_xx /= sumw;
192 cl_yy /= sumw;
193 cl_zz /= sumw;
194 cl_xx = sqrt(cl_xx - cl_x * cl_x);
195 cl_yy = sqrt(cl_yy - cl_y * cl_y);
196 cl_zz = sqrt(cl_zz - cl_z * cl_z);
197
198 cluster.setRMSXYZ(cl_xx, cl_yy, cl_zz);
199 cluster.setLayer(interm_cluster[cluster_indx].layer());
200 ldmx_log(trace) << "Cluster " << cluster_indx
201 << " energy: " << cluster.getEnergy()
202 << ", nHits: " << cluster.getNHits() << ", centroid: ("
207 << cluster.
getRMSZ() <<
")" <<
", first layer centroid: ("
208 << cluster.getFirstLayerCentroidX() << ", "
209 << cluster.getFirstLayerCentroidY() << ", "
210 << cluster.getFirstLayerCentroidZ() << ")";
211
212 ecal_clusters.push_back(cluster);
213 }
214 ldmx_log(debug) << "Filled " << ecal_clusters.size()
215 << " clusters into ecal_clusters";
216 event.add(cluster_coll_name_, ecal_clusters);
217 } else {
218 ldmx_log(info) <<
"Using simple clustering algorithm " <<
algo_name_;
220
222
223 if (hit.getEnergy() == 0) {
224 continue;
225 }
227 }
228
229 cf.
cluster(seed_threshold_, cutoff_);
231 std::map<int, double> c_weights = cf.
getWeights();
232
238
239 std::map<int, double>::iterator it = c_weights.begin();
240 for (it = c_weights.begin(); it != c_weights.end(); it++) {
241 algo_result.
setWeight(it->first, it->second / 100);
242 }
243
244 std::vector<ldmx::EcalCluster> ecal_clusters;
245 for (size_t cluster_indx = 0; cluster_indx < interm_cluster.size();
246 cluster_indx++) {
248
249 cluster.
setEnergy(interm_cluster[cluster_indx].energy());
251 interm_cluster[cluster_indx].centroidY(),
252 interm_cluster[cluster_indx].centroidZ());
253 cluster.
setNHits(interm_cluster[cluster_indx].hits().size());
254 cluster.
addHits(interm_cluster[cluster_indx].hits());
255 ecal_clusters.push_back(cluster);
256 }
257
258 event.add(cluster_coll_name_, ecal_clusters);
259 event.add(algo_coll_name_, algo_result);
260 }
261}
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.
void setIDs(std::vector< unsigned int > &hitIDs)
Sets a sorted vector for the IDs of the hits_ that make up 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.