Process the event and make histograms or summaries.
50 {
51 std::vector<double> bdt_features;
52
54 "SimParticles", sim_particles_pass_name_)};
55
56
57
58
59
60 std::vector<double> gamma_p(3);
61 std::vector<double> gamma_x0(3);
62
63 std::vector<double> recoil_p(3);
64 bool found_recoil_e = false;
65
66 if (recoil_from_tracking_) {
67 const auto &recoil_tracks{
68 event.getCollection<
ldmx::Track>(track_collection_, track_pass_name_)};
69
70 for (auto &track : recoil_tracks) {
71
72 if (track.q() == 1 && track.getNhits() == 5) {
73 gamma_x0 = track.getPosition();
74 gamma_p[0] = -1. * track.getMomentum()[0];
75 gamma_p[1] = -1. * track.getMomentum()[1];
76 gamma_p[2] = beam_energy_mev_ - track.getMomentum()[2];
77 }
78 }
79 } else {
80 if (event.
exists(sp_collection_, sp_pass_name_)) {
82 sp_collection_, sp_pass_name_);
83 bool found_rec = false;
84 for (auto const &it : particle_map) {
85 for (auto const &sphit : target_sp_hits) {
86 if (sphit.getPosition()[2] > 0) {
87 if (it.first == sphit.getTrackID()) {
88 if (it.second.getPdgID() == 622) {
89 std::vector<float> x0_float = sphit.getPosition();
90 std::vector<double> x0_double(x0_float.begin(), x0_float.end());
91 gamma_x0 = x0_double;
92 gamma_p = sphit.getMomentum();
93 found_rec = true;
94 }
95 if (it.second.getPdgID() == 11 &&
96 inList(it.second.getParents(), 0)) {
97 if (!found_rec) {
98 std::vector<float> x0_float = sphit.getPosition();
99 std::vector<double> x0_double(x0_float.begin(),
100 x0_float.end());
101 gamma_x0 = x0_double;
102 gamma_p[0] = -1. * sphit.getMomentum()[0];
103 gamma_p[1] = -1. * sphit.getMomentum()[1];
104 gamma_p[2] = beam_energy_mev_ - sphit.getMomentum()[2];
105 found_rec = true;
106 }
107 recoil_p = sphit.getMomentum();
108 found_recoil_e = true;
109 }
110 }
111 }
112 }
113 }
114 }
115 }
116
117 double p_mag = 0.;
118 if (found_recoil_e) {
119 p_mag = std::sqrt(recoil_p[0] * recoil_p[0] + recoil_p[1] * recoil_p[1] +
120 recoil_p[2] * recoil_p[2]);
121 }
122
123 const auto &ecal_rec_hits =
event.getCollection<
ldmx::EcalHit>(
124 ecal_rec_collection_, ecal_rec_pass_name_);
125 const auto &hcal_rec_hits =
event.getCollection<
ldmx::HcalHit>(
126 hcal_rec_collection_, hcal_rec_pass_name_);
127
128 double ecal_energy = 0.;
129 double hcal_energy = 0.;
130 bool hcal_containment = true;
131
133 if (hit.getEnergy() > 0.) {
134 ecal_energy += hit.getEnergy();
135 }
136 }
138 if (hit.getEnergy() > 0.) {
140 if (det_id.getSection() != 0) {
141 continue;
142 }
143 if (det_id.getLayerID() == 1 && hit.getPE() > 5) {
144 hcal_containment = false;
145 }
146 hcal_energy += 12. * hit.getEnergy();
147 }
148 }
149
150 if (ecal_energy < 3160 && hcal_energy > 4840 && hcal_containment &&
151 p_mag < 2400) {
152
153 int n_layers_hit = 0;
154 double x_std = 0.;
155 double y_std = 0.;
156 double z_std = 0.;
157 double x_mean = 0.;
158 double y_mean = 0.;
159 double r_mean = 0.;
160 int iso_hits = 0;
161 double iso_energy = 0.;
162 int n_readout_hits = 0;
163 double summed_det = 0.;
164 double r_mean_from_photon_track = 0.;
165
166 double z_mean = 0.;
167 std::vector<int> layers_hit;
168
170 if (hit.getEnergy() > 0.) {
172 if (det_id.getSection() != 0) {
173 continue;
174 }
175 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
176 continue;
177 }
178
179 n_readout_hits += 1;
180 double hit_x = hit.getXPos();
181 double hit_y = hit.getYPos();
182 double hit_z = hit.getZPos();
183 double hit_r = sqrt(hit_x * hit_x + hit_y * hit_y);
184
185 summed_det += hit.getEnergy();
186
187 x_mean += hit_x * hit.getEnergy();
188 y_mean += hit_y * hit.getEnergy();
189 z_mean += hit_z * hit.getEnergy();
190 r_mean += hit_r * hit.getEnergy();
191
192
193 if (!(std::find(layers_hit.begin(), layers_hit.end(),
194 det_id.getLayerID()) != layers_hit.end())) {
195 layers_hit.push_back(det_id.getLayerID());
196 }
197
198 double proj_x =
199 gamma_x0[0] + (hit_z - gamma_x0[2]) * gamma_p[0] / gamma_p[2];
200 double proj_y =
201 gamma_x0[1] + (hit_z - gamma_x0[2]) * gamma_p[1] / gamma_p[2];
202
203 r_mean_from_photon_track +=
204 hit.getEnergy() * sqrt((hit_x - proj_x) * (hit_x - proj_x) +
205 (hit_y - proj_y) * (hit_y - proj_y));
206
207
208 double closest_point = 9999.;
210 if (hit2.getEnergy() > 0.) {
212 if (fabs(hit2.getXPos()) > 1000 || fabs(hit2.getYPos()) > 1000) {
213 continue;
214 }
215 if (det_i_d2.getLayerID() == det_id.getLayerID()) {
216
217
218
219 if (hit.isOrientationX()) {
220 if (fabs(hit2.getYPos() - hit_y) > 0) {
221 if (fabs(hit2.getYPos() - hit_y) < closest_point) {
222 closest_point = fabs(hit2.getYPos() - hit_y);
223 }
224 }
225 }
226 if (hit.isOrientationY()) {
227 if (fabs(hit2.getXPos() - hit_x) > 0) {
228 if (fabs(hit2.getXPos() - hit_x) < closest_point) {
229 closest_point = fabs(hit2.getXPos() - hit_x);
230 }
231 }
232 }
233 }
234 }
235 }
236 if (closest_point > 50.) {
237 iso_hits += 1;
238 iso_energy += hit.getEnergy();
239 }
240 }
241 }
242
243 n_layers_hit = layers_hit.size();
244
245 if (summed_det > 0.) {
246 x_mean /= summed_det;
247 y_mean /= summed_det;
248 z_mean /= summed_det;
249 r_mean /= summed_det;
250
251 r_mean_from_photon_track /= summed_det;
252 }
253
255 if (hit.getEnergy() > 0.) {
256 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
257 continue;
258 }
260 if (det_id.getSection() == 0) {
261 x_std += hit.getEnergy() * (hit.getXPos() - x_mean) *
262 (hit.getXPos() - x_mean);
263 y_std += hit.getEnergy() * (hit.getYPos() - y_mean) *
264 (hit.getYPos() - y_mean);
265 z_std += hit.getEnergy() * (hit.getZPos() - z_mean) *
266 (hit.getZPos() - z_mean);
267 }
268 }
269 }
270
271 if (summed_det > 0.) {
272 x_std = sqrt(x_std / summed_det);
273 y_std = sqrt(y_std / summed_det);
274 z_std = sqrt(z_std / summed_det);
275 }
276
277
290
291 bdt_features.push_back(n_layers_hit);
292 bdt_features.push_back(x_std);
293 bdt_features.push_back(y_std);
294 bdt_features.push_back(z_std);
295 bdt_features.push_back(x_mean);
296 bdt_features.push_back(y_mean);
297 bdt_features.push_back(r_mean);
298 bdt_features.push_back(iso_hits);
299 bdt_features.push_back(iso_energy);
300 bdt_features.push_back(n_readout_hits);
301 bdt_features.push_back(hcal_energy);
302 bdt_features.push_back(r_mean_from_photon_track);
303
304 if (training_) {
305 std::ofstream file(training_file_, std::ios::app);
306 if (!file.is_open()) {
307 ldmx_log(fatal) << "Error: Could not open file " << training_file_;
308 return;
309 }
310 for (int i = 0; i < bdt_features.size(); ++i) {
311 file << bdt_features[i] << (i + 1 == bdt_features.size() ? "\n" : ", ");
312 }
313 }
314 }
315
316 return;
317}
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.