208 if (!hTracks->empty()) {
221 if (hOpticalFunctions->empty()) {
222 edm::LogInfo(
"CTPPSProtonProducer") <<
"No optical functions available, reconstruction disabled.";
234 std::ostringstream ssLog;
236 ssLog <<
"* input tracks:" << std::endl;
239 std::map<unsigned int, std::vector<unsigned int>> trackingSelection, timingSelection;
241 for (
unsigned int idx = 0;
idx < hTracks->size(); ++
idx) {
242 const auto &tr = hTracks->at(
idx);
258 <<
"[" <<
idx <<
"] " << tr.rpId() <<
" (" << (rpId.arm() * 100 + rpId.station() * 10 + rpId.rp())
260 <<
"x=" << tr.x() <<
" +- " << tr.xUnc() <<
" mm, " 261 <<
"y=" << tr.y() <<
" +- " << tr.yUnc() <<
" mm" << std::endl;
263 const bool trackerRP =
267 trackingSelection[rpId.arm()].push_back(
idx);
269 timingSelection[rpId.arm()].push_back(
idx);
273 for (
const auto &arm_it : trackingSelection) {
274 const auto &
indices = arm_it.second;
279 std::map<unsigned int, reco::ForwardProton> singleRPResultsIndexed;
283 ssLog << std::endl <<
"* reconstruction from track " <<
idx << std::endl;
285 singleRPResultsIndexed[
idx] =
293 std::set<unsigned int>
rpIds;
295 rpIds.insert(hTracks->at(
idx).rpId());
300 std::vector<std::pair<unsigned int, unsigned int>> idx_pairs;
301 std::map<unsigned int, unsigned int> idx_pair_multiplicity;
304 const auto &tr_i = hTracks->at(
i);
305 const auto &tr_j = hTracks->at(
j);
310 const auto &pr_i = singleRPResultsIndexed[
i];
311 const auto &pr_j = singleRPResultsIndexed[
j];
313 if (tr_i.rpId() == tr_j.rpId())
319 bool matching =
true;
321 if (!ac.isSatisfied(ac.qX, tr_i.x(), tr_i.y(), lhcInfoCombined.crossingAngle(), tr_i.x() - tr_j.x()))
323 else if (!ac.isSatisfied(ac.qY, tr_i.x(), tr_i.y(), lhcInfoCombined.crossingAngle(), tr_i.y() - tr_j.y()))
325 else if (!ac.isSatisfied(
326 ac.qXi, tr_i.x(), tr_i.y(), lhcInfoCombined.crossingAngle(), pr_i.xi() - pr_j.xi()))
328 else if (!ac.isSatisfied(ac.qThetaY,
331 lhcInfoCombined.crossingAngle(),
332 pr_i.thetaY() - pr_j.thetaY()))
338 idx_pairs.emplace_back(
i,
j);
339 idx_pair_multiplicity[
i]++;
340 idx_pair_multiplicity[
j]++;
345 std::map<unsigned int, unsigned int> timing_RP_track_multiplicity;
346 for (
const auto &ti : timingSelection[arm_it.first]) {
347 const auto &tr = hTracks->at(ti);
348 timing_RP_track_multiplicity[tr.rpId()]++;
352 std::map<unsigned int, std::vector<unsigned int>> matched_timing_track_indices;
353 std::map<unsigned int, unsigned int> matched_timing_track_multiplicity;
354 for (
unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
355 const auto &
i = idx_pairs[pr_idx].first;
356 const auto &
j = idx_pairs[pr_idx].second;
359 if (idx_pair_multiplicity[
i] > 1 || idx_pair_multiplicity[
j] > 1)
362 const auto &tr_i = hTracks->at(
i);
363 const auto &tr_j = hTracks->at(
j);
368 for (
const auto &ti : timingSelection[arm_it.first]) {
369 const auto &tr_ti = hTracks->at(ti);
377 const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j);
378 const double x_inter = f_i * tr_i.x() + f_j * tr_j.x();
379 const double x_inter_unc_sq =
380 f_i * f_i * tr_i.xUnc() * tr_i.xUnc() + f_j * f_j * tr_j.xUnc() * tr_j.xUnc();
382 const double de_x = tr_ti.x() - x_inter;
383 const double de_x_unc =
sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_sq);
384 const double r = (de_x_unc > 0.) ?
de_x / de_x_unc : 1E100;
386 const bool matching = (ac.getTiTrMin() <
r &&
r < ac.getTiTrMax());
389 ssLog <<
"ti=" << ti <<
", i=" <<
i <<
", j=" <<
j <<
" | z_ti=" << z_ti <<
", z_i=" << z_i
390 <<
", z_j=" << z_j <<
" | x_ti=" << tr_ti.x() <<
", x_inter=" << x_inter <<
", de_x=" <<
de_x 391 <<
", de_x_unc=" << de_x_unc <<
", matching=" << matching << std::endl;
396 matched_timing_track_indices[pr_idx].push_back(ti);
397 matched_timing_track_multiplicity[ti]++;
402 for (
unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
403 const auto &
i = idx_pairs[pr_idx].first;
404 const auto &
j = idx_pairs[pr_idx].second;
407 if (idx_pair_multiplicity[
i] > 1 || idx_pair_multiplicity[
j] > 1)
412 <<
"* reconstruction from tracking-RP tracks: " <<
i <<
", " <<
j <<
" and timing-RP tracks: ";
422 double sw = 0., swt = 0.;
423 for (
const auto &ti : matched_timing_track_indices[pr_idx]) {
425 if (matched_timing_track_multiplicity[ti] > 1)
433 const auto &tr = hTracks->at(ti);
434 const double t_unc = tr.timeUnc();
435 const double w = (t_unc > 0.) ? 1. / t_unc / t_unc : 1.;
437 swt +=
w * tr.time();
443 time_unc = 1. /
sqrt(sw);
447 ssLog << std::endl <<
" time = " <<
time <<
" +- " << time_unc << std::endl;
458 pOutMultiRP->emplace_back(proton);
463 for (
const auto &
p : singleRPResultsIndexed)
464 pOutSingleRP->emplace_back(
p.second);
const bool useNewLHCInfo_
const edm::ESGetToken< LHCInfoPerFill, LHCInfoPerFillRcd > lhcInfoPerFillToken_
reco::ForwardProton reconstructFromSingleRP(const CTPPSLocalTrackLiteRef &track, const float energy, std::ostream &os) const
run proton reconstruction using single-RP strategy
bool doMultiRPReconstruction_
bool pixelDiscardBXShiftedTracks_
const edm::ESGetToken< PPSAssociationCuts, PPSAssociationCutsRcd > ppsAssociationCutsToken_
unsigned int max_n_timing_tracks_
const edm::ESGetToken< LHCInfoPerLS, LHCInfoPerLSRcd > lhcInfoPerLSToken_
void init(const LHCInterpolatedOpticalFunctionsSetCollection &opticalFunctions)
const CutsPerArm & getAssociationCuts(const int sector) const
void setTimeError(float time_err)
const edm::ESGetToken< LHCInfo, LHCInfoRcd > lhcInfoToken_
edm::ESWatcher< CTPPSInterpolatedOpticsRcd > opticsWatcher_
Abs< T >::type abs(const T &t)
edm::Ref< CTPPSLocalTrackLiteCollection > CTPPSLocalTrackLiteRef
Persistent reference to a CTPPSLocalTrackLite.
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
reco::ForwardProton reconstructFromMultiRP(const CTPPSLocalTrackLiteRefVector &tracks, const float energy, std::ostream &os) const
run proton reconstruction using multiple-RP strategy
Log< level::Info, false > LogInfo
Vector rpTranslation(unsigned int id) const
bool check(const edm::EventSetup &iSetup)
Base class for CTPPS detector IDs.
ProtonReconstructionAlgorithm algorithm_
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
std::vector< ForwardProton > ForwardProtonCollection
Collection of ForwardProton objects.
void setContributingLocalTracks(const CTPPSLocalTrackLiteRefVector &v)
store the list of RP tracks that contributed to this global track
const edm::ESGetToken< CTPPSGeometry, VeryForwardRealGeometryRecord > geometryToken_
std::string singleRPReconstructionLabel_
const edm::ESGetToken< LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd > opticalFunctionsToken_
std::string multiRPReconstructionLabel_
bool doSingleRPReconstruction_
edm::EDGetTokenT< CTPPSLocalTrackLiteCollection > tracksToken_