Process the event and put new data products into it.
36 {
37 const auto& ecal_hits{
event.getCollection<
ldmx::EcalHit>(rec_hit_coll_name_,
38 rec_hit_pass_name_)};
39 if (ecal_hits.size() == 0) {
40
41 ldmx_log(fatal) << "No ECal hits found... exiting";
42 return;
43 }
44
45 if (clue_) {
46 ldmx_log(info) << "Using CLUE clustering algorithm";
47 CLUE cf;
48 cf.cluster(ecal_hits, dc_, rhoc_, deltac_, deltao_, nbr_of_layers_,
49 reclustering_);
50 ldmx_log(debug) << "CLUE algorithm finished clustering";
51 std::vector<IntermediateCluster> interm_cluster = cf.getClusters();
52 std::vector<IntermediateCluster> f_interm_cluster =
53 cf.getFirstLayerCentroids();
54 ldmx_log(debug) << "Got " << interm_cluster.size() << " clusters with "
55 << f_interm_cluster.size() << " first layer centroids";
56
57 auto n_loops = cf.getNLoops();
58 ldmx_log(debug) << "Number of clustering loops: " << n_loops;
59
60 if (reclustering_) {
61 ldmx_log(debug) << "Reclustererd initial number of clusters: "
62 << cf.getInitialClusterNbr()
63 << ", final number of clusters: "
64 << interm_cluster.size();
65 }
66
67 std::vector<ldmx::EcalCluster> ecal_clusters;
68 ldmx_log(debug) << "Filling " << interm_cluster.size()
69 << " clusters into ecal_clusters";
70 for (int cluster_indx = 0; cluster_indx < interm_cluster.size();
71 cluster_indx++) {
73
74 cluster.
setEnergy(interm_cluster[cluster_indx].centroid().E());
75 cluster.
setCentroidXYZ(interm_cluster[cluster_indx].centroid().x(),
76 interm_cluster[cluster_indx].centroid().y(),
77 interm_cluster[cluster_indx].centroid().z());
78 cluster.setFirstLayerCentroidXYZ(
79 f_interm_cluster[cluster_indx].centroid().x(),
80 f_interm_cluster[cluster_indx].centroid().y(),
81 f_interm_cluster[cluster_indx].centroid().z());
82 cluster.
setNHits(interm_cluster[cluster_indx].getHits().size());
83 cluster.
addHits(interm_cluster[cluster_indx].getHits());
84 cluster.addFirstLayerHits(f_interm_cluster[cluster_indx].getHits());
85
86 float cl_x(0), cl_y(0), cl_z(0), cl_xx(0), cl_yy(0), cl_zz(0);
87 float cl_w = 1;
88 float sumw = 0;
89
90 for (auto hit : interm_cluster[cluster_indx].getHits()) {
91 if (hit->getEnergy() < min_hit_energy_) continue;
92 cl_w = log(hit->getEnergy()) - log(min_hit_energy_);
93 cl_x += cl_w * hit->getXPos();
94 cl_y += cl_w * hit->getYPos();
95 cl_z += cl_w * hit->getZPos();
96 cl_xx += cl_w * hit->getXPos() * hit->getXPos();
97 cl_yy += cl_w * hit->getYPos() * hit->getYPos();
98 cl_zz += cl_w * hit->getZPos() * hit->getZPos();
99 sumw += cl_w;
100 }
101
102 cl_x /= sumw;
103 cl_y /= sumw;
104 cl_z /= sumw;
105 cl_xx /= sumw;
106 cl_yy /= sumw;
107 cl_zz /= sumw;
108 cl_xx = sqrt(cl_xx - cl_x * cl_x);
109 cl_yy = sqrt(cl_yy - cl_y * cl_y);
110 cl_zz = sqrt(cl_zz - cl_z * cl_z);
111
112 cluster.setRMSXYZ(cl_xx, cl_yy, cl_zz);
113 cluster.setLayer(interm_cluster[cluster_indx].getLayer());
114 ldmx_log(trace) << "Cluster " << cluster_indx
115 << " energy: " << cluster.getEnergy()
116 << ", nHits: " << cluster.getNHits() << ", centroid: ("
121 << cluster.
getRMSZ() <<
")" <<
", first layer centroid: ("
122 << cluster.getFirstLayerCentroidX() << ", "
123 << cluster.getFirstLayerCentroidY() << ", "
124 << cluster.getFirstLayerCentroidZ() << ")";
125
126 ecal_clusters.push_back(cluster);
127 }
128 ldmx_log(debug) << "Filled " << ecal_clusters.size()
129 << " clusters into ecal_clusters";
130 event.add(cluster_coll_name_, ecal_clusters);
131 } else {
132 ldmx_log(info) <<
"Using simple clustering algorithm " <<
algo_name_;
133 TemplatedClusterFinder<MyClusterWeight> cf;
134
136
137 if (hit.getEnergy() == 0) {
138 continue;
139 }
140 cf.add(hit);
141 }
142
143 cf.cluster(seed_threshold_, cutoff_);
144 std::vector<IntermediateCluster> interm_cluster = cf.getClusters();
145 std::map<int, double> c_weights = cf.getWeights();
146
152
153 std::map<int, double>::iterator it = c_weights.begin();
154 for (it = c_weights.begin(); it != c_weights.end(); it++) {
155 algo_result.
setWeight(it->first, it->second / 100);
156 }
157
158 std::vector<ldmx::EcalCluster> ecal_clusters;
159 for (int cluster_indx = 0; cluster_indx < interm_cluster.size();
160 cluster_indx++) {
162
163 cluster.
setEnergy(interm_cluster[cluster_indx].centroid().E());
164 cluster.
setCentroidXYZ(interm_cluster[cluster_indx].centroid().x(),
165 interm_cluster[cluster_indx].centroid().y(),
166 interm_cluster[cluster_indx].centroid().z());
167 cluster.setLayer(interm_cluster[cluster_indx].getLayer());
168 std::cout <<
"Cluster layer: " << cluster.
getLayer() << std::endl;
169 cluster.
setNHits(interm_cluster[cluster_indx].getHits().size());
170 cluster.
addHits(interm_cluster[cluster_indx].getHits());
171 ecal_clusters.push_back(cluster);
172 }
173
174 event.add(cluster_coll_name_, ecal_clusters);
175 event.add(algo_coll_name_, algo_result);
176 }
177}
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.
int getLayer() const
layer 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.