242 if (!hTracks->empty()) {
255 if (hOpticalFunctions->empty()) {
256 edm::LogInfo(
"CTPPSProtonProducer") <<
"No optical functions available, reconstruction disabled.";
268 std::ostringstream ssLog;
270 ssLog <<
"* input tracks:" << std::endl;
273 std::map<unsigned int, std::vector<unsigned int>> trackingSelection, timingSelection;
275 for (
unsigned int idx = 0;
idx < hTracks->size(); ++
idx) {
276 const auto &tr = hTracks->at(
idx);
292 <<
"[" <<
idx <<
"] " << tr.rpId() <<
" (" << (
rpId.arm() * 100 +
rpId.station() * 10 +
rpId.rp())
294 <<
"x=" << tr.x() <<
" +- " << tr.xUnc() <<
" mm, "
295 <<
"y=" << tr.y() <<
" +- " << tr.yUnc() <<
" mm" << std::endl;
297 const bool trackerRP =
301 trackingSelection[
rpId.arm()].push_back(
idx);
303 timingSelection[
rpId.arm()].push_back(
idx);
307 for (
const auto &arm_it : trackingSelection) {
308 const auto &
indices = arm_it.second;
313 std::map<unsigned int, reco::ForwardProton> singleRPResultsIndexed;
317 ssLog << std::endl <<
"* reconstruction from track " <<
idx << std::endl;
319 singleRPResultsIndexed[
idx] =
327 std::set<unsigned int>
rpIds;
329 rpIds.insert(hTracks->at(
idx).rpId());
334 std::vector<std::pair<unsigned int, unsigned int>> idx_pairs;
335 std::map<unsigned int, unsigned int> idx_pair_multiplicity;
338 const auto &tr_i = hTracks->at(
i);
339 const auto &tr_j = hTracks->at(
j);
344 const auto &pr_i = singleRPResultsIndexed[
i];
345 const auto &pr_j = singleRPResultsIndexed[
j];
347 if (tr_i.rpId() == tr_j.rpId())
355 if (ac.x_cut_apply &&
std::abs(tr_i.x() - tr_j.x() - ac.x_cut_mean) > ac.x_cut_value)
357 else if (ac.y_cut_apply &&
std::abs(tr_i.y() - tr_j.y() - ac.y_cut_mean) > ac.y_cut_value)
359 else if (ac.xi_cut_apply &&
std::abs(pr_i.xi() - pr_j.xi() - ac.xi_cut_mean) > ac.xi_cut_value)
361 else if (ac.th_y_cut_apply &&
362 std::abs(pr_i.thetaY() - pr_j.thetaY() - ac.th_y_cut_mean) > ac.th_y_cut_value)
368 idx_pairs.emplace_back(
i,
j);
369 idx_pair_multiplicity[
i]++;
370 idx_pair_multiplicity[
j]++;
375 std::map<unsigned int, unsigned int> timing_RP_track_multiplicity;
376 for (
const auto &ti : timingSelection[arm_it.first]) {
377 const auto &tr = hTracks->at(ti);
378 timing_RP_track_multiplicity[tr.rpId()]++;
382 std::map<unsigned int, std::vector<unsigned int>> matched_timing_track_indices;
383 std::map<unsigned int, unsigned int> matched_timing_track_multiplicity;
384 for (
unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
385 const auto &
i = idx_pairs[pr_idx].first;
386 const auto &
j = idx_pairs[pr_idx].second;
389 if (idx_pair_multiplicity[
i] > 1 || idx_pair_multiplicity[
j] > 1)
392 const auto &tr_i = hTracks->at(
i);
393 const auto &tr_j = hTracks->at(
j);
398 for (
const auto &ti : timingSelection[arm_it.first]) {
399 const auto &tr_ti = hTracks->at(ti);
408 const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j);
409 const double x_inter = f_i * tr_i.x() + f_j * tr_j.x();
410 const double x_inter_unc_sq =
411 f_i * f_i * tr_i.xUnc() * tr_i.xUnc() + f_j * f_j * tr_j.xUnc() * tr_j.xUnc();
413 const double de_x = tr_ti.x() - x_inter;
414 const double de_x_unc =
sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_sq);
415 const double r = (de_x_unc > 0.) ?
de_x / de_x_unc : 1E100;
417 const bool matching = (ac.ti_tr_min <
r &&
r < ac.ti_tr_max);
420 ssLog <<
"ti=" << ti <<
", i=" <<
i <<
", j=" <<
j <<
" | z_ti=" << z_ti <<
", z_i=" << z_i
421 <<
", z_j=" << z_j <<
" | x_ti=" << tr_ti.x() <<
", x_inter=" << x_inter <<
", de_x=" <<
de_x
422 <<
", de_x_unc=" << de_x_unc <<
", matching=" <<
matching << std::endl;
427 matched_timing_track_indices[pr_idx].push_back(ti);
428 matched_timing_track_multiplicity[ti]++;
433 for (
unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
434 const auto &
i = idx_pairs[pr_idx].first;
435 const auto &
j = idx_pairs[pr_idx].second;
438 if (idx_pair_multiplicity[
i] > 1 || idx_pair_multiplicity[
j] > 1)
443 <<
"* reconstruction from tracking-RP tracks: " <<
i <<
", " <<
j <<
" and timing-RP tracks: ";
453 double sw = 0., swt = 0.;
454 for (
const auto &ti : matched_timing_track_indices[pr_idx]) {
456 if (matched_timing_track_multiplicity[ti] > 1)
464 const auto &tr = hTracks->at(ti);
465 const double t_unc = tr.timeUnc();
466 const double w = (t_unc > 0.) ? 1. / t_unc / t_unc : 1.;
468 swt +=
w * tr.time();
474 time_unc = 1. /
sqrt(sw);
478 ssLog << std::endl <<
" time = " <<
time <<
" +- " << time_unc << std::endl;
488 pOutMultiRP->emplace_back(proton);
493 for (
const auto &
p : singleRPResultsIndexed)
494 pOutSingleRP->emplace_back(
std::move(
p.second));