Process the event and make histograms or summaries.
53 {
54 std::vector<double> bdt_features;
55
57 sim_particles_coll_name_, sim_particles_pass_name_)};
58
59
60
61
62
63 std::vector<double> gamma_p(3);
64 std::vector<double> gamma_x0(3);
65
66 std::vector<double> recoil_p(3);
67 bool found_recoil_e = false;
68
69 if (recoil_from_tracking_) {
70 const auto &recoil_tracks{
71 event.getCollection<
ldmx::Track>(track_collection_, track_pass_name_)};
72
73 for (auto &track : recoil_tracks) {
74
75 auto trk_pos = track.getPositionAtTarget();
76 auto trk_mom = track.getMomentumAtTarget();
77 if (track.getCharge() == 1 && track.getNhits() == 5 &&
78 trk_pos.size() == 3 && trk_mom.size() == 3) {
79 gamma_x0 = trk_pos;
80 gamma_p[0] = -1. * trk_mom[0];
81 gamma_p[1] = -1. * trk_mom[1];
82 gamma_p[2] = beam_energy_mev_ - trk_mom[2];
83 }
84 }
85 } else {
86 if (event.
exists(sp_collection_, sp_pass_name_)) {
88 sp_collection_, sp_pass_name_);
89 bool found_rec = false;
90 for (auto const &it : particle_map) {
91 for (auto const &sphit : target_sp_hits) {
92 if (sphit.getPosition()[2] > 0) {
93 if (it.first == sphit.getTrackID()) {
94 if (it.second.getPdgID() == 622) {
95 std::vector<float> x0_float = sphit.getPosition();
96 std::vector<double> x0_double(x0_float.begin(), x0_float.end());
97 gamma_x0 = x0_double;
98 gamma_p = sphit.getMomentum();
99 found_rec = true;
100 }
101 if (it.second.getPdgID() == 11 &&
102 inList(it.second.getParents(), 0)) {
103 if (!found_rec) {
104 std::vector<float> x0_float = sphit.getPosition();
105 std::vector<double> x0_double(x0_float.begin(),
106 x0_float.end());
107 gamma_x0 = x0_double;
108 gamma_p[0] = -1. * sphit.getMomentum()[0];
109 gamma_p[1] = -1. * sphit.getMomentum()[1];
110 gamma_p[2] = beam_energy_mev_ - sphit.getMomentum()[2];
111 found_rec = true;
112 }
113 recoil_p = sphit.getMomentum();
114 found_recoil_e = true;
115 }
116 }
117 }
118 }
119 }
120 }
121 }
122
123 double p_mag = 0.;
124 if (found_recoil_e) {
125 p_mag = std::sqrt(recoil_p[0] * recoil_p[0] + recoil_p[1] * recoil_p[1] +
126 recoil_p[2] * recoil_p[2]);
127 }
128
129 const auto &ecal_rec_hits =
event.getCollection<
ldmx::EcalHit>(
130 ecal_rec_collection_, ecal_rec_pass_name_);
131 const auto &hcal_rec_hits =
event.getCollection<
ldmx::HcalHit>(
132 hcal_rec_collection_, hcal_rec_pass_name_);
133
134 double ecal_energy = 0.;
135 double hcal_energy = 0.;
136 bool hcal_containment = true;
137
139 if (hit.getEnergy() > 0.) {
140 ecal_energy += hit.getEnergy();
141 }
142 }
144 if (hit.getEnergy() > 0.) {
146 if (det_id.getSection() != 0) {
147 continue;
148 }
149 if (det_id.getLayerID() == 1 && hit.getPE() > 5) {
150 hcal_containment = false;
151 }
152 hcal_energy += 12. * hit.getEnergy();
153 }
154 }
155
156 if (ecal_energy < 3160 && hcal_energy > 4840 && hcal_containment &&
157 p_mag < 2400) {
158
159 int n_layers_hit = 0;
160 double x_std = 0.;
161 double y_std = 0.;
162 double z_std = 0.;
163 double x_mean = 0.;
164 double y_mean = 0.;
165 double r_mean = 0.;
166 int iso_hits = 0;
167 double iso_energy = 0.;
168 int n_readout_hits = 0;
169 double summed_det = 0.;
170 double r_mean_from_photon_track = 0.;
171
172 double z_mean = 0.;
173 std::vector<int> layers_hit;
174
176 if (hit.getEnergy() > 0.) {
178 if (det_id.getSection() != 0) {
179 continue;
180 }
181 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
182 continue;
183 }
184
185 n_readout_hits += 1;
186 double hit_x = hit.getXPos();
187 double hit_y = hit.getYPos();
188 double hit_z = hit.getZPos();
189 double hit_r = sqrt(hit_x * hit_x + hit_y * hit_y);
190
191 summed_det += hit.getEnergy();
192
193 x_mean += hit_x * hit.getEnergy();
194 y_mean += hit_y * hit.getEnergy();
195 z_mean += hit_z * hit.getEnergy();
196 r_mean += hit_r * hit.getEnergy();
197
198
199 if (!(std::find(layers_hit.begin(), layers_hit.end(),
200 det_id.getLayerID()) != layers_hit.end())) {
201 layers_hit.push_back(det_id.getLayerID());
202 }
203
204 double proj_x =
205 gamma_x0[0] + (hit_z - gamma_x0[2]) * gamma_p[0] / gamma_p[2];
206 double proj_y =
207 gamma_x0[1] + (hit_z - gamma_x0[2]) * gamma_p[1] / gamma_p[2];
208
209 r_mean_from_photon_track +=
210 hit.getEnergy() * sqrt((hit_x - proj_x) * (hit_x - proj_x) +
211 (hit_y - proj_y) * (hit_y - proj_y));
212
213
214 double closest_point = 9999.;
216 if (hit2.getEnergy() > 0.) {
218 if (fabs(hit2.getXPos()) > 1000 || fabs(hit2.getYPos()) > 1000) {
219 continue;
220 }
221 if (det_i_d2.getLayerID() == det_id.getLayerID()) {
222
223
224
225 if (hit.isOrientationX()) {
226 if (fabs(hit2.getYPos() - hit_y) > 0) {
227 if (fabs(hit2.getYPos() - hit_y) < closest_point) {
228 closest_point = fabs(hit2.getYPos() - hit_y);
229 }
230 }
231 }
232 if (hit.isOrientationY()) {
233 if (fabs(hit2.getXPos() - hit_x) > 0) {
234 if (fabs(hit2.getXPos() - hit_x) < closest_point) {
235 closest_point = fabs(hit2.getXPos() - hit_x);
236 }
237 }
238 }
239 }
240 }
241 }
242 if (closest_point > 50.) {
243 iso_hits += 1;
244 iso_energy += hit.getEnergy();
245 }
246 }
247 }
248
249 n_layers_hit = layers_hit.size();
250
251 if (summed_det > 0.) {
252 x_mean /= summed_det;
253 y_mean /= summed_det;
254 z_mean /= summed_det;
255 r_mean /= summed_det;
256
257 r_mean_from_photon_track /= summed_det;
258 }
259
261 if (hit.getEnergy() > 0.) {
262 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
263 continue;
264 }
266 if (det_id.getSection() == 0) {
267 x_std += hit.getEnergy() * (hit.getXPos() - x_mean) *
268 (hit.getXPos() - x_mean);
269 y_std += hit.getEnergy() * (hit.getYPos() - y_mean) *
270 (hit.getYPos() - y_mean);
271 z_std += hit.getEnergy() * (hit.getZPos() - z_mean) *
272 (hit.getZPos() - z_mean);
273 }
274 }
275 }
276
277 if (summed_det > 0.) {
278 x_std = sqrt(x_std / summed_det);
279 y_std = sqrt(y_std / summed_det);
280 z_std = sqrt(z_std / summed_det);
281 }
282
283
296
297 bdt_features.push_back(n_layers_hit);
298 bdt_features.push_back(x_std);
299 bdt_features.push_back(y_std);
300 bdt_features.push_back(z_std);
301 bdt_features.push_back(x_mean);
302 bdt_features.push_back(y_mean);
303 bdt_features.push_back(r_mean);
304 bdt_features.push_back(iso_hits);
305 bdt_features.push_back(iso_energy);
306 bdt_features.push_back(n_readout_hits);
307 bdt_features.push_back(hcal_energy);
308 bdt_features.push_back(r_mean_from_photon_track);
309
310 if (training_) {
311 std::ofstream file(training_file_, std::ios::app);
312 if (!file.is_open()) {
313 ldmx_log(fatal) << "Error: Could not open file " << training_file_;
314 return;
315 }
316 for (int i = 0; i < bdt_features.size(); ++i) {
317 file << bdt_features[i] << (i + 1 == bdt_features.size() ? "\n" : ", ");
318 }
319 }
320 }
321
322 return;
323}
HistogramPool histograms_
helper object for making and filling histograms
bool exists(const std::string &name, const std::string &passName, bool unique=true) const
Check for the existence of an object or collection with the given name and pass name in the event.
void fill(const std::string &name, const T &val)
Fill a 1D histogram.
Stores reconstructed hit information from the ECAL.
Stores reconstructed hit information from the HCAL.
Implements detector ids for HCal subdetector.
Class representing a simulated particle.
Represents a simulated tracker hit in the simulation.
Implementation of a track object.