189 if (!hTracks->empty()) {
202 if (hOpticalFunctions->empty()) {
203 edm::LogInfo(
"CTPPSProtonProducer") <<
"No optical functions available, reconstruction disabled.";
215 std::ostringstream ssLog;
217 ssLog <<
"* input tracks:" << std::endl;
220 std::map<unsigned int, std::vector<unsigned int>> trackingSelection, timingSelection;
222 for (
unsigned int idx = 0; idx < hTracks->size(); ++idx) {
223 const auto &tr = hTracks->at(idx);
239 <<
"[" << idx <<
"] " << tr.
rpId() <<
" (" << (rpId.arm() * 100 + rpId.station() * 10 + rpId.rp())
241 <<
"x=" << tr.x() <<
" +- " << tr.xUnc() <<
" mm, "
242 <<
"y=" << tr.y() <<
" +- " << tr.yUnc() <<
" mm" << std::endl;
244 const bool trackerRP =
248 trackingSelection[rpId.arm()].push_back(idx);
250 timingSelection[rpId.arm()].push_back(idx);
254 for (
const auto &arm_it : trackingSelection) {
255 const auto &
indices = arm_it.second;
257 const auto &ac = ppsAssociationCuts->getAssociationCuts(arm_it.first);
260 std::map<unsigned int, reco::ForwardProton> singleRPResultsIndexed;
262 for (
const auto &idx :
indices) {
264 ssLog << std::endl <<
"* reconstruction from track " << idx << std::endl;
266 singleRPResultsIndexed[idx] =
274 std::set<unsigned int>
rpIds;
275 for (
const auto &idx : indices)
276 rpIds.insert(hTracks->at(idx).rpId());
281 std::vector<std::pair<unsigned int, unsigned int>> idx_pairs;
282 std::map<unsigned int, unsigned int> idx_pair_multiplicity;
283 for (
const auto &
i : indices) {
284 for (
const auto &
j : indices) {
285 const auto &tr_i = hTracks->at(
i);
286 const auto &tr_j = hTracks->at(
j);
288 const double z_i = hGeometry->rpTranslation(tr_i.rpId()).
z();
289 const double z_j = hGeometry->rpTranslation(tr_j.rpId()).
z();
291 const auto &pr_i = singleRPResultsIndexed[
i];
292 const auto &pr_j = singleRPResultsIndexed[
j];
294 if (tr_i.rpId() == tr_j.rpId())
300 bool matching =
true;
302 if (!ac.isSatisfied(ac.qX, tr_i.x(), tr_i.y(), hLHCInfo->crossingAngle(), tr_i.x() - tr_j.x()))
304 else if (!ac.isSatisfied(ac.qY, tr_i.x(), tr_i.y(), hLHCInfo->crossingAngle(), tr_i.y() - tr_j.y()))
306 else if (!ac.isSatisfied(ac.qXi, tr_i.x(), tr_i.y(), hLHCInfo->crossingAngle(), pr_i.xi() - pr_j.xi()))
308 else if (!ac.isSatisfied(
309 ac.qThetaY, tr_i.x(), tr_i.y(), hLHCInfo->crossingAngle(), pr_i.thetaY() - pr_j.thetaY()))
315 idx_pairs.emplace_back(
i,
j);
316 idx_pair_multiplicity[
i]++;
317 idx_pair_multiplicity[
j]++;
322 std::map<unsigned int, unsigned int> timing_RP_track_multiplicity;
323 for (
const auto &ti : timingSelection[arm_it.first]) {
324 const auto &tr = hTracks->at(ti);
325 timing_RP_track_multiplicity[tr.rpId()]++;
329 std::map<unsigned int, std::vector<unsigned int>> matched_timing_track_indices;
330 std::map<unsigned int, unsigned int> matched_timing_track_multiplicity;
331 for (
unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
332 const auto &
i = idx_pairs[pr_idx].first;
333 const auto &
j = idx_pairs[pr_idx].second;
336 if (idx_pair_multiplicity[
i] > 1 || idx_pair_multiplicity[
j] > 1)
339 const auto &tr_i = hTracks->at(
i);
340 const auto &tr_j = hTracks->at(
j);
342 const double z_i = hGeometry->rpTranslation(tr_i.rpId()).
z();
343 const double z_j = hGeometry->rpTranslation(tr_j.rpId()).
z();
345 for (
const auto &ti : timingSelection[arm_it.first]) {
346 const auto &tr_ti = hTracks->at(ti);
353 const double z_ti = hGeometry->rpTranslation(tr_ti.rpId()).
z();
354 const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j);
355 const double x_inter = f_i * tr_i.x() + f_j * tr_j.x();
356 const double x_inter_unc_sq =
357 f_i * f_i * tr_i.xUnc() * tr_i.xUnc() + f_j * f_j * tr_j.xUnc() * tr_j.xUnc();
359 const double de_x = tr_ti.x() - x_inter;
360 const double de_x_unc =
sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_sq);
361 const double r = (de_x_unc > 0.) ? de_x / de_x_unc : 1E100;
363 const bool matching = (ac.getTiTrMin() < r && r < ac.getTiTrMax());
366 ssLog <<
"ti=" << ti <<
", i=" <<
i <<
", j=" <<
j <<
" | z_ti=" << z_ti <<
", z_i=" << z_i
367 <<
", z_j=" << z_j <<
" | x_ti=" << tr_ti.x() <<
", x_inter=" << x_inter <<
", de_x=" << de_x
368 <<
", de_x_unc=" << de_x_unc <<
", matching=" << matching << std::endl;
373 matched_timing_track_indices[pr_idx].push_back(ti);
374 matched_timing_track_multiplicity[ti]++;
379 for (
unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
380 const auto &
i = idx_pairs[pr_idx].first;
381 const auto &
j = idx_pairs[pr_idx].second;
384 if (idx_pair_multiplicity[
i] > 1 || idx_pair_multiplicity[
j] > 1)
389 <<
"* reconstruction from tracking-RP tracks: " <<
i <<
", " <<
j <<
" and timing-RP tracks: ";
399 double sw = 0., swt = 0.;
400 for (
const auto &ti : matched_timing_track_indices[pr_idx]) {
402 if (matched_timing_track_multiplicity[ti] > 1)
410 const auto &tr = hTracks->at(ti);
411 const double t_unc = tr.timeUnc();
412 const double w = (t_unc > 0.) ? 1. / t_unc / t_unc : 1.;
414 swt += w * tr.time();
420 time_unc = 1. /
sqrt(sw);
424 ssLog << std::endl <<
" time = " << time <<
" +- " << time_unc << std::endl;
434 pOutMultiRP->emplace_back(proton);
439 for (
const auto &
p : singleRPResultsIndexed)
440 pOutSingleRP->emplace_back(
p.second);
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
reco::ForwardProton reconstructFromSingleRP(const CTPPSLocalTrackLiteRef &track, const LHCInfo &lhcInfo, std::ostream &os) const
run proton reconstruction using single-RP strategy
bool doMultiRPReconstruction_
bool pixelDiscardBXShiftedTracks_
edm::ESGetToken< CTPPSGeometry, VeryForwardRealGeometryRecord > geometryToken_
unsigned int max_n_timing_tracks_
void init(const LHCInterpolatedOpticalFunctionsSetCollection &opticalFunctions)
void setTimeError(float time_err)
reco::ForwardProton reconstructFromMultiRP(const CTPPSLocalTrackLiteRefVector &tracks, const LHCInfo &lhcInfo, std::ostream &os) const
run proton reconstruction using multiple-RP strategy
edm::ESWatcher< CTPPSInterpolatedOpticsRcd > opticsWatcher_
Abs< T >::type abs(const T &t)
edm::Ref< CTPPSLocalTrackLiteCollection > CTPPSLocalTrackLiteRef
Persistent reference to a CTPPSLocalTrackLite.
Log< level::Info, false > LogInfo
bool check(const edm::EventSetup &iSetup)
Base class for CTPPS detector IDs.
ProtonReconstructionAlgorithm algorithm_
edm::ESGetToken< PPSAssociationCuts, PPSAssociationCutsRcd > ppsAssociationCutsToken_
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
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
std::string singleRPReconstructionLabel_
std::string multiRPReconstructionLabel_
bool doSingleRPReconstruction_
edm::EDGetTokenT< CTPPSLocalTrackLiteCollection > tracksToken_
edm::ESGetToken< LHCInfo, LHCInfoRcd > lhcInfoToken_
edm::ESGetToken< LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd > opticalFunctionsToken_