43 auto start = std::chrono::high_resolution_clock::now();
50 ldmx::EcalGeometry::CONDITIONS_OBJECT_NAME);
54 ecal_collection_name_, ecal_pass_name_);
56 mip_collection_name_, mip_pass_name_);
57 std::vector<XYCoords> ele_trajectory, photon_trajectory,
58 ele_trajectory_at_target;
59 std::vector<ldmx::HitData> tracking_hit_list;
60 ldmx_log(trace) <<
"EcalMipTrackingProcessor::produce() called";
61 ele_trajectory = ecal_trajectory_info.getEleTrajectory();
62 photon_trajectory = ecal_trajectory_info.getPhotonTrajectory();
63 tracking_hit_list = ecal_trajectory_info.getTrackingHitList();
64 n_readout_hits_ = ecal_veto_result.getNReadoutHits();
76 ROOT::Math::XYZVector e_traj_start;
77 ROOT::Math::XYZVector e_traj_end;
78 ROOT::Math::XYZVector p_traj_start;
79 ROOT::Math::XYZVector p_traj_end;
80 if (!ele_trajectory.empty() && !photon_trajectory.empty()) {
83 e_traj_start.SetXYZ(ele_trajectory[0].first, ele_trajectory[0].second,
85 e_traj_end.SetXYZ(ele_trajectory[(n_ecal_layers_ - 1)].first,
86 ele_trajectory[(n_ecal_layers_ - 1)].second,
88 p_traj_start.SetXYZ(photon_trajectory[0].first, photon_trajectory[0].second,
90 p_traj_end.SetXYZ(photon_trajectory[(n_ecal_layers_ - 1)].first,
91 photon_trajectory[(n_ecal_layers_ - 1)].second,
97 e_traj_end = ROOT::Math::XYZVector(
101 p_traj_end = ROOT::Math::XYZVector(
111 ldmx_log(trace) <<
"Finding first near photon layer_";
112 if (!photon_trajectory.empty()) {
113 for (std::vector<ldmx::HitData>::iterator it = tracking_hit_list.begin();
114 it != tracking_hit_list.end(); ++it) {
116 sqrt(pow((*it).pos_.X() - photon_trajectory[(*it).layer_].first, 2) +
117 pow((*it).pos_.Y() - photon_trajectory[(*it).layer_].second, 2));
130 ROOT::Math::XYZVector g_toe = (e_traj_start - p_traj_start).Unit();
132 ROOT::Math::XYZVector origin = p_traj_start + 0.5 * 8.7 * g_toe;
133 ldmx_log(trace) <<
"Origin of photon territory: " << origin.X() <<
", "
134 << origin.Y() <<
", " << origin.Z();
135 if (!ele_trajectory.empty()) {
136 for (
auto& hit_data : tracking_hit_list) {
137 ROOT::Math::XYZVector hit_pos = hit_data.pos_;
138 ROOT::Math::XYZVector hit_prime = hit_pos - origin;
139 if (hit_prime.Dot(g_toe) <= 0) {
152 tracking_hit_list.begin(), tracking_hit_list.end(),
158 std::vector<std::vector<ldmx::HitData>> track_list;
162 ldmx_log(trace) <<
"====== Tracking hit list (original) length "
163 << tracking_hit_list.size() <<
" ======";
164 for (
int i = 0; i < tracking_hit_list.size(); i++) {
165 ldmx_log(trace) <<
"[" << tracking_hit_list[i].pos_.X() <<
", "
166 << tracking_hit_list[i].pos_.Y() <<
", "
167 << tracking_hit_list[i].layer_ <<
"], ";
169 ldmx_log(trace) <<
"====== END OF Tracking hit list ======";
174 for (
int i_hit = 0; i_hit < tracking_hit_list.size(); i_hit++) {
177 int current_hit{i_hit};
187 while (j_hit < tracking_hit_list.size()) {
188 if ((tracking_hit_list[j_hit].layer_ ==
189 tracking_hit_list[current_hit].layer_ - 1 ||
190 tracking_hit_list[j_hit].layer_ ==
191 tracking_hit_list[current_hit].layer_ - 2) &&
192 std::abs(tracking_hit_list[j_hit].pos_.X() -
193 tracking_hit_list[current_hit].pos_.X()) <=
195 std::abs(tracking_hit_list[j_hit].pos_.Y() -
196 tracking_hit_list[current_hit].pos_.Y()) <=
198 track[track_len] = j_hit;
206 if (track_len < 2)
continue;
208 tracking_hit_list[track[0]].pos_,
209 tracking_hit_list[track[track_len - 1]].pos_, e_traj_start, e_traj_end);
211 tracking_hit_list[track[0]].pos_,
212 tracking_hit_list[track[track_len - 1]].pos_, p_traj_start, p_traj_end);
215 if (closest_p > cell_width and closest_e < 2 * cell_width)
continue;
216 if (track_len < 4 and closest_e > closest_p)
continue;
218 ldmx_log(debug) <<
"====== After rejection for MIP tracking ======";
219 ldmx_log(debug) <<
"current hit: [" << tracking_hit_list[i_hit].pos_.X()
220 <<
", " << tracking_hit_list[i_hit].pos_.Y() <<
", "
221 << tracking_hit_list[i_hit].layer_ <<
"]";
223 for (
int k = 0; k < track_len; k++) {
224 ldmx_log(debug) <<
"track[" << k <<
"] position = ["
225 << tracking_hit_list[track[k]].pos_.X() <<
", "
226 << tracking_hit_list[track[k]].pos_.Y() <<
", "
227 << tracking_hit_list[track[k]].layer_ <<
"]";
232 if (track_len >= 2) {
233 std::vector<ldmx::HitData> temp_track_list;
235 for (
int k_hit = 0; k_hit < track_len; k_hit++) {
236 temp_track_list.push_back(tracking_hit_list[track[k_hit] - n_remove]);
237 tracking_hit_list.erase(tracking_hit_list.begin() + track[k_hit] -
243 ldmx_log(trace) <<
"====== Tracking hit list (after erase) length "
244 << tracking_hit_list.size() <<
" ======";
245 for (
int i = 0; i < tracking_hit_list.size(); i++) {
246 ldmx_log(trace) <<
"[" << tracking_hit_list[i].pos_.X() <<
", "
247 << tracking_hit_list[i].pos_.Y() <<
", "
248 << tracking_hit_list[i].layer_ <<
"] ";
250 ldmx_log(trace) <<
"====== END OF Tracking hit list ======";
252 track_list.push_back(temp_track_list);
260 ldmx_log(debug) <<
"Straight tracks found (before merge): "
261 << track_list.size();
263 for (
int i_track = 0; i_track < track_list.size(); i_track++) {
264 ldmx_log(trace) <<
"Track " << i_track <<
":";
265 for (
int i_hit = 0; i_hit < track_list[i_track].size(); i_hit++) {
266 ldmx_log(trace) <<
" Hit " << i_hit <<
": ["
267 << track_list[i_track][i_hit].pos_.X() <<
", "
268 << track_list[i_track][i_hit].pos_.Y() <<
", "
269 << track_list[i_track][i_hit].layer_ <<
"]" << std::endl;
276 ldmx_log(debug) <<
"Beginning track merging using " << track_list.size()
279 for (
int track_i = 0; track_i < track_list.size(); track_i++) {
282 std::vector<ldmx::HitData> base_track = track_list[track_i];
285 ldmx_log(trace) <<
" Considering track " << track_i;
286 for (
int track_j = track_i + 1; track_j < track_list.size(); track_j++) {
287 std::vector<ldmx::HitData> checking_track = track_list[track_j];
288 if (checking_track.empty()) {
289 ldmx_log(error) <<
"Broken logic: a straight ecal track had no hits in "
295 if ((head_hitdata.layer_ == tail_hitdata.layer_ + 1 ||
296 head_hitdata.layer_ == tail_hitdata.layer_ + 2) &&
297 pow(pow(head_hitdata.pos_.X() - tail_hitdata.pos_.X(), 2) +
298 pow(head_hitdata.pos_.Y() - tail_hitdata.pos_.Y(), 2),
299 0.5) <= cell_width) {
303 ldmx_log(trace) <<
" ** Compatible track found at index_ "
305 ldmx_log(trace) <<
" Tail xylayer: " << head_hitdata.pos_.X() <<
","
306 << head_hitdata.pos_.Y() <<
"," << head_hitdata.layer_;
307 ldmx_log(trace) <<
" Head xylayer: " << tail_hitdata.pos_.X() <<
","
308 << tail_hitdata.pos_.Y() <<
"," << tail_hitdata.layer_;
309 for (
int hit_k = 0; hit_k < checking_track.size(); hit_k++) {
310 base_track.push_back(track_list[track_j][hit_k]);
312 track_list[track_i] = base_track;
313 track_list.erase(track_list.begin() + track_j);
320 ldmx_log(debug) <<
"Straight tracks found (after merge): "
322 for (
int track_i = 0; track_i < track_list.size(); track_i++) {
323 ldmx_log(debug) <<
"Track " << track_i <<
":";
324 for (
int hit_i = 0; hit_i < track_list[track_i].size(); hit_i++) {
325 ldmx_log(debug) <<
" Hit " << hit_i <<
": ["
326 << track_list[track_i][hit_i].pos_.X() <<
", "
327 << track_list[track_i][hit_i].pos_.Y() <<
", "
328 << track_list[track_i][hit_i].layer_ <<
"]";
332 auto straight_tracks = std::chrono::high_resolution_clock::now();
333 profiling_map_[
"straight_tracks"] +=
334 std::chrono::duration<double, std::milli>(straight_tracks - start)
338 ldmx_log(info) <<
"Finding linreg tracks from a total of "
339 << tracking_hit_list.size() <<
" hits using a radius of "
340 << linreg_radius_ <<
" mm";
342 for (
int i_hit = 0; i_hit < 0; i_hit++) {
345 std::vector<int> hits_in_region;
348 ROOT::Math::XYZVector slope_vec;
349 ROOT::Math::XYZVector h_mean;
350 ROOT::Math::XYZVector h_point;
351 float r_corr_best{0.0};
355 int best_hit_nums[3];
357 hits_in_region.push_back(i_hit);
359 for (
int j_hit = 0; j_hit < tracking_hit_list.size(); j_hit++) {
361 if (tracking_hit_list[i_hit].pos_.Z() ==
362 tracking_hit_list[j_hit].pos_.Z()) {
366 (tracking_hit_list[i_hit].pos_ - tracking_hit_list[j_hit].pos_).R();
370 if (dist_to_hit <= 2 * linreg_radius_) {
371 hits_in_region.push_back(j_hit);
375 bool best_lin_reg_found{
false};
377 ldmx_log(debug) <<
"There are " << hits_in_region.size()
378 <<
" hits within a radius of " << linreg_radius_ <<
" mm";
382 for (
int j_hit_in_reg = 1; j_hit_in_reg < hits_in_region.size() - 1;
385 if (hits_in_region.size() < 3)
break;
386 hit_nums[1] = hits_in_region[j_hit_in_reg];
387 for (
int k_hit_reg = j_hit_in_reg + 1; k_hit_reg < hits_in_region.size();
389 hit_nums[2] = hits_in_region[k_hit_reg];
390 const auto& p0 = tracking_hit_list[hit_nums[0]].pos_;
391 const auto& p1 = tracking_hit_list[hit_nums[1]].pos_;
392 const auto& p2 = tracking_hit_list[hit_nums[2]].pos_;
394 h_mean = (p0 + p1 + p2) / 3.0;
396 double p_arr[3][3] = {{p0.X(), p0.Y(), p0.Z()},
397 {p1.X(), p1.Y(), p1.Z()},
398 {p2.X(), p2.Y(), p2.Z()}};
401 double mean_arr[3] = {(p0.X() + p1.X() + p2.X()) / 3.0,
402 (p0.Y() + p1.Y() + p2.Y()) / 3.0,
403 (p0.Z() + p1.Z() + p2.Z()) / 3.0};
405 for (
int h_ind = 0; h_ind < 3; ++h_ind) {
406 for (
int l_ind = 0; l_ind < 3; ++l_ind) {
407 hdt(h_ind, l_ind) = p_arr[h_ind][l_ind] - mean_arr[l_ind];
414 hdt(0, 0) * (hdt(1, 1) * hdt(2, 2) - hdt(1, 2) * hdt(2, 1)) -
415 hdt(0, 1) * (hdt(1, 0) * hdt(2, 2) - hdt(1, 2) * hdt(2, 0)) +
416 hdt(0, 2) * (hdt(1, 0) * hdt(2, 1) - hdt(1, 1) * hdt(2, 0));
418 if (determinant == 0)
continue;
420 TDecompSVD svd_obj(hdt);
421 bool decomposed = svd_obj.Decompose();
422 if (!decomposed)
continue;
426 slope_vec.SetX(vm[0][0]);
427 slope_vec.SetY(vm[0][1]);
428 slope_vec.SetZ(vm[0][2]);
430 h_point = slope_vec + h_mean;
440 if (closest_p > cell_width or closest_e < 1.5 * cell_width)
continue;
443 float vrnc = (tracking_hit_list[hit_nums[0]].pos_ - h_mean).R() +
444 (tracking_hit_list[hit_nums[1]].pos_ - h_mean).R() +
445 (tracking_hit_list[hit_nums[2]].pos_ - h_mean).R();
453 float r_corr = 1 - sumerr / vrnc;
457 if (r_corr > r_corr_best and r_corr > .6) {
458 r_corr_best = r_corr;
460 best_lin_reg_found =
true;
461 for (
int k = 0; k < 3; k++) {
462 best_hit_nums[k] = hit_nums[k];
469 if (!best_lin_reg_found)
continue;
473 for (
int final_hit_index = 0; final_hit_index < 3; final_hit_index++) {
475 <<
" Hit " << final_hit_index <<
" ["
476 << tracking_hit_list[best_hit_nums[final_hit_index]].pos_.X() <<
", "
477 << tracking_hit_list[best_hit_nums[final_hit_index]].pos_.Y() <<
", "
478 << tracking_hit_list[best_hit_nums[final_hit_index]].pos_.Z() <<
"] ";
482 for (
int l_hit = 0; l_hit < 3; l_hit++) {
483 tracking_hit_list.erase(tracking_hit_list.begin() + best_hit_nums[l_hit]);
489 <<
" lin-reg tracks";
491 auto linreg_tracks = std::chrono::high_resolution_clock::now();
492 profiling_map_[
"linreg_tracks"] +=
493 std::chrono::duration<double, std::milli>(linreg_tracks - straight_tracks)
500 event.add(mip_result_name_, mip_result);
502 auto end = std::chrono::high_resolution_clock::now();
503 auto time_diff = end - start;
505 std::chrono::duration<double, std::milli>(time_diff).count();