Process the event and make histograms or summaries.
61 {
62 std::vector<float> bdt_features;
63
64 double decay_z;
65
67 sim_particles_coll_name_, sim_particles_pass_name_)};
68
69 for (auto const &it : particle_map) {
70 std::vector<int> parents = it.second.getParents();
71 for (int track_id : parents) {
72 if (track_id == 0 && it.second.getPdgID() == -11) {
73 decay_z = it.second.getVertex()[2];
74 }
75 }
76 }
77
79
80
81
82
83
84 std::vector<double> gamma_p(3);
85 double gamma_energy = 90000.;
86 std::vector<double> gamma_x0(3);
87
88 std::vector<double> recoil_p(3);
89 bool found_recoil_e = false;
90
91 if (recoil_from_tracking_) {
92 const auto &recoil_tracks{
93 event.getCollection<
ldmx::Track>(track_collection_, track_pass_name_)};
94 for (auto &track : recoil_tracks) {
95
96 if (track.q() == 1 && track.getNhits() == 5) {
97 gamma_x0 = track.getPosition();
98 gamma_p[0] = -1. * track.getMomentum()[0];
99 gamma_p[1] = -1. * track.getMomentum()[1];
100 gamma_p[2] = 8000. - track.getMomentum()[2];
101 }
102 }
103 } else {
104 if (event.
exists(sp_collection_, sp_pass_name_)) {
106 sp_collection_, sp_pass_name_);
107 bool found_rec = false;
108 for (auto const &it : particle_map) {
109 for (auto const &sphit : target_sp_hits) {
110 if (sphit.getPosition()[2] > 0) {
111 if (it.first == sphit.getTrackID()) {
112 if (it.second.getPdgID() == 622) {
113 std::vector<float> x0_float = sphit.getPosition();
114 std::vector<double> x0_double(x0_float.begin(), x0_float.end());
115 gamma_x0 = x0_double;
116 gamma_p = sphit.getMomentum();
117 gamma_energy = sphit.getEnergy();
118 found_rec = true;
119 }
120 if (it.second.getPdgID() == 11 &&
121 inList(it.second.getParents(), 0)) {
122 if (!found_rec) {
123 std::vector<float> x0_float = sphit.getPosition();
124 std::vector<double> x0_double(x0_float.begin(),
125 x0_float.end());
126 gamma_x0 = x0_double;
127 gamma_p[0] = -1. * sphit.getMomentum()[0];
128 gamma_p[1] = -1. * sphit.getMomentum()[1];
129 gamma_p[2] = beam_energy_mev_ - sphit.getMomentum()[2];
130 gamma_energy = beam_energy_mev_ - sphit.getEnergy();
131 found_rec = true;
132 }
133 recoil_p = sphit.getMomentum();
134 found_recoil_e = true;
135 }
136 }
137 }
138 }
139 }
140 }
141 }
144 std::acos(gamma_p[2] / std::sqrt(gamma_p[0] * gamma_p[0] +
145 gamma_p[1] * gamma_p[1] +
146 gamma_p[2] * gamma_p[2])));
147
148 double p_mag = 0.;
149 if (found_recoil_e) {
150 p_mag = std::sqrt(recoil_p[0] * recoil_p[0] + recoil_p[1] * recoil_p[1] +
151 recoil_p[2] * recoil_p[2]);
152 }
153
154 if (p_mag < 50. && all_cuts_) {
155 return;
156 }
158
159
160 const auto &ecal_rec_hits =
event.getCollection<
ldmx::EcalHit>(
161 ecal_rec_collection_, ecal_rec_pass_name_);
162 const auto &hcal_rec_hits =
event.getCollection<
ldmx::HcalHit>(
163 hcal_rec_collection_, hcal_rec_pass_name_);
164
165 double trigger = 0.;
166 double ecal_energy = 0.;
167 double hcal_energy = 0.;
168 bool hcal_containment = true;
169
171 if (hit.getEnergy() > 0.) {
172 ecal_energy += hit.getEnergy();
173 if (hit.getZPos() <= 541.722) {
174 trigger += hit.getEnergy();
175 }
176 }
177 }
179 if (hit.getEnergy() > 0.) {
181 if (det_id.getSection() != 0) {
182 continue;
183 }
184 if (det_id.getLayerID() == 1 && hit.getPE() > 5) {
185 hcal_containment = false;
186 }
187 hcal_energy += 12. * hit.getEnergy();
188 }
189 }
190
191 if (trigger > 3160. && all_cuts_) {
192 return;
193 }
195
196 if (ecal_energy > 3160. && all_cuts_) {
197 return;
198 }
200
201 if (p_mag > 2400. && all_cuts_) {
202 return;
203 }
205
207 ecal_veto_collection_, ecal_veto_pass_)};
208 if (ecal_veto.getDisc() < ecal_bdt_cut_val_ && all_cuts_) {
209 return;
210 }
212
213 if (hcal_energy < 4840 && all_cuts_) {
214 return;
215 }
217
218 if (!hcal_containment && all_cuts_) {
219 return;
220 }
222
223
224 int n_layers_hit = 0;
225 double x_std = 0.;
226 double y_std = 0.;
227 double z_std = 0.;
228 double x_mean = 0.;
229 double y_mean = 0.;
230 double r_mean = 0.;
231 int iso_hits = 0;
232 double iso_energy = 0.;
233 int n_readout_hits = 0;
234 double summed_det = 0.;
235 double r_mean_from_photon_track = 0.;
236
237 double z_mean = 0.;
238 std::vector<int> layers_hit;
240 if (hit.getEnergy() > 0.) {
242 if (det_id.getSection() != 0) {
243 continue;
244 }
245 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
246 continue;
247 }
248 n_readout_hits += 1;
249 double hit_x = hit.getXPos();
250 double hit_y = hit.getYPos();
251 double hit_z = hit.getZPos();
252 double hit_r = sqrt(hit_x * hit_x + hit_y * hit_y);
253
254 summed_det += hit.getEnergy();
255
256 x_mean += hit_x * hit.getEnergy();
257 y_mean += hit_y * hit.getEnergy();
258 z_mean += hit_z * hit.getEnergy();
259 r_mean += hit_r * hit.getEnergy();
260
261
262 if (!(std::find(layers_hit.begin(), layers_hit.end(),
263 det_id.getLayerID()) != layers_hit.end())) {
264 layers_hit.push_back(det_id.getLayerID());
265 }
266
267 double proj_x =
268 gamma_x0[0] + (hit_z - gamma_x0[2]) * gamma_p[0] / gamma_p[2];
269 double proj_y =
270 gamma_x0[1] + (hit_z - gamma_x0[2]) * gamma_p[1] / gamma_p[2];
271
272 r_mean_from_photon_track +=
273 hit.getEnergy() * sqrt((hit_x - proj_x) * (hit_x - proj_x) +
274 (hit_y - proj_y) * (hit_y - proj_y));
275
276
277 double closest_point = 9999.;
279 if (hit2.getEnergy() > 0.) {
281 if (fabs(hit2.getXPos()) > 1000 || fabs(hit2.getYPos()) > 1000) {
282 continue;
283 }
284 if (det_i_d2.getLayerID() == det_id.getLayerID()) {
285
286
287
288 if (hit.isOrientationX()) {
289 if (fabs(hit2.getYPos() - hit_y) > 0) {
290 if (fabs(hit2.getYPos() - hit_y) < closest_point) {
291 closest_point = abs(hit2.getYPos() - hit_y);
292 }
293 }
294 }
295 if (hit.isOrientationY()) {
296 if (fabs(hit2.getXPos() - hit_x) > 0) {
297 if (fabs(hit2.getXPos() - hit_x) < closest_point) {
298 closest_point = fabs(hit2.getXPos() - hit_x);
299 }
300 }
301 }
302 }
303 }
304 }
305 if (closest_point > 50.) {
306 iso_hits += 1;
307 iso_energy += hit.getEnergy();
308 }
309 }
310 }
311
312 n_layers_hit = layers_hit.size();
313
314 if (summed_det > 0.) {
315 x_mean /= summed_det;
316 y_mean /= summed_det;
317 z_mean /= summed_det;
318 r_mean /= summed_det;
319
320 r_mean_from_photon_track /= summed_det;
321 }
322
324 if (hit.getEnergy() > 0.) {
325 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
326 continue;
327 }
329 if (det_id.getSection() == 0) {
330 x_std += hit.getEnergy() * (hit.getXPos() - x_mean) *
331 (hit.getXPos() - x_mean);
332 y_std += hit.getEnergy() * (hit.getYPos() - y_mean) *
333 (hit.getYPos() - y_mean);
334 z_std += hit.getEnergy() * (hit.getZPos() - z_mean) *
335 (hit.getZPos() - z_mean);
336 }
337 }
338 }
339
340 if (summed_det > 0.) {
341 x_std = sqrt(x_std / summed_det);
342 y_std = sqrt(y_std / summed_det);
343 z_std = sqrt(z_std / summed_det);
344 }
345
346
359
360 bdt_features.push_back(n_layers_hit);
361 bdt_features.push_back(x_std);
362 bdt_features.push_back(y_std);
363 bdt_features.push_back(z_std);
364 bdt_features.push_back(x_mean);
365 bdt_features.push_back(y_mean);
366 bdt_features.push_back(r_mean);
367 bdt_features.push_back(iso_hits);
368 bdt_features.push_back(iso_energy);
369 bdt_features.push_back(n_readout_hits);
370 bdt_features.push_back(hcal_energy);
371 bdt_features.push_back(r_mean_from_photon_track);
372
373 ldmx::ort::FloatArrays inputs({bdt_features});
374 float pred =
375 rt_->run({feature_list_name_}, inputs, {"probabilities"})[0].at(1);
376
378
379 for (int i = 0; i < 10000; i++) {
380 double disc = 0.999 + (double(i) / 10000.) * (1. - 0.999);
382 0.999 + ((double(i) + 0.5) / 10000.) * (1. - 0.999));
383 double disc_roc = 0.99 + (double(i) / 10000.) * (1. - 0.99);
384 if (pred >= disc) {
386 0.999 + ((double(i) + 0.5) / 10000.) * (1. - 0.999));
387 }
388 if (pred >= disc_roc) {
390 }
391 }
392
394
395 if (pred < bdt_cut_val_) {
396 return;
397 }
399
400 return;
401}
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.