252 if (!hTracks->empty()) {
263 if (hOpticalFunctions->empty()) {
264 edm::LogInfo(
"CTPPSProtonProducer") <<
"No optical functions available, reconstruction disabled.";
276 std::ostringstream ssLog;
278 ssLog <<
"* input tracks:" << std::endl;
281 std::map<unsigned int, std::vector<unsigned int>> trackingSelection, timingSelection;
283 for (
unsigned int idx = 0;
idx < hTracks->size(); ++
idx) {
284 const auto &tr = hTracks->at(
idx);
300 <<
"[" <<
idx <<
"] " << tr.rpId() <<
" (" << (
rpId.arm() * 100 +
rpId.station() * 10 +
rpId.rp())
302 <<
"x=" << tr.x() <<
" +- " << tr.xUnc() <<
" mm, "
303 <<
"y=" << tr.y() <<
" +- " << tr.yUnc() <<
" mm" << std::endl;
305 const bool trackerRP =
309 trackingSelection[
rpId.arm()].push_back(
idx);
311 timingSelection[
rpId.arm()].push_back(
idx);
315 for (
const auto &arm_it : trackingSelection) {
316 const auto &
indices = arm_it.second;
321 std::map<unsigned int, reco::ForwardProton> singleRPResultsIndexed;
325 ssLog << std::endl <<
"* reconstruction from track " <<
idx << std::endl;
327 singleRPResultsIndexed[
idx] =
335 std::set<unsigned int>
rpIds;
337 rpIds.insert(hTracks->at(
idx).rpId());
342 std::vector<std::pair<unsigned int, unsigned int>> idx_pairs;
343 std::map<unsigned int, unsigned int> idx_pair_multiplicity;
346 const auto &tr_i = hTracks->at(
i);
347 const auto &tr_j = hTracks->at(
j);
352 const auto &pr_i = singleRPResultsIndexed[
i];
353 const auto &pr_j = singleRPResultsIndexed[
j];
355 if (tr_i.rpId() == tr_j.rpId())
363 if (ac.x_cut_apply &&
std::abs(tr_i.x() - tr_j.x() - ac.x_cut_mean) > ac.x_cut_value)
365 else if (ac.y_cut_apply &&
std::abs(tr_i.y() - tr_j.y() - ac.y_cut_mean) > ac.y_cut_value)
367 else if (ac.xi_cut_apply &&
std::abs(pr_i.xi() - pr_j.xi() - ac.xi_cut_mean) > ac.xi_cut_value)
369 else if (ac.th_y_cut_apply &&
370 std::abs(pr_i.thetaY() - pr_j.thetaY() - ac.th_y_cut_mean) > ac.th_y_cut_value)
376 idx_pairs.emplace_back(
i,
j);
377 idx_pair_multiplicity[
i]++;
378 idx_pair_multiplicity[
j]++;
383 std::map<unsigned int, unsigned int> timing_RP_track_multiplicity;
384 for (
const auto &ti : timingSelection[arm_it.first]) {
385 const auto &tr = hTracks->at(ti);
386 timing_RP_track_multiplicity[tr.rpId()]++;
390 std::map<unsigned int, std::vector<unsigned int>> matched_timing_track_indices;
391 std::map<unsigned int, unsigned int> matched_timing_track_multiplicity;
392 for (
unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
393 const auto &
i = idx_pairs[pr_idx].first;
394 const auto &
j = idx_pairs[pr_idx].second;
397 if (idx_pair_multiplicity[
i] > 1 || idx_pair_multiplicity[
j] > 1)
400 const auto &tr_i = hTracks->at(
i);
401 const auto &tr_j = hTracks->at(
j);
406 for (
const auto &ti : timingSelection[arm_it.first]) {
407 const auto &tr_ti = hTracks->at(ti);
415 const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j);
416 const double x_inter = f_i * tr_i.x() + f_j * tr_j.x();
417 const double x_inter_unc_sq =
418 f_i * f_i * tr_i.xUnc() * tr_i.xUnc() + f_j * f_j * tr_j.xUnc() * tr_j.xUnc();
420 const double de_x = tr_ti.x() - x_inter;
421 const double de_x_unc =
sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_sq);
422 const double r = (de_x_unc > 0.) ?
de_x / de_x_unc : 1E100;
424 const bool matching = (ac.ti_tr_min <
r &&
r < ac.ti_tr_max);
427 ssLog <<
"ti=" << ti <<
", i=" <<
i <<
", j=" <<
j <<
" | z_ti=" << z_ti <<
", z_i=" << z_i
428 <<
", z_j=" << z_j <<
" | x_ti=" << tr_ti.x() <<
", x_inter=" << x_inter <<
", de_x=" <<
de_x
429 <<
", de_x_unc=" << de_x_unc <<
", matching=" <<
matching << std::endl;
434 matched_timing_track_indices[pr_idx].push_back(ti);
435 matched_timing_track_multiplicity[ti]++;
440 for (
unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
441 const auto &
i = idx_pairs[pr_idx].first;
442 const auto &
j = idx_pairs[pr_idx].second;
445 if (idx_pair_multiplicity[
i] > 1 || idx_pair_multiplicity[
j] > 1)
450 <<
"* reconstruction from tracking-RP tracks: " <<
i <<
", " <<
j <<
" and timing-RP tracks: ";
460 double sw = 0., swt = 0.;
461 for (
const auto &ti : matched_timing_track_indices[pr_idx]) {
463 if (matched_timing_track_multiplicity[ti] > 1)
471 const auto &tr = hTracks->at(ti);
472 const double t_unc = tr.timeUnc();
473 const double w = (t_unc > 0.) ? 1. / t_unc / t_unc : 1.;
475 swt +=
w * tr.time();
481 time_unc = 1. /
sqrt(sw);
485 ssLog << std::endl <<
" time = " <<
time <<
" +- " << time_unc << std::endl;
495 pOutMultiRP->emplace_back(proton);
500 for (
const auto &
p : singleRPResultsIndexed)
501 pOutSingleRP->emplace_back(
p.second);