ProtonReconstruction/plugins/CTPPSProtonProducer.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  *   Laurent Forthomme
0005  ****************************************************************************/
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/ESWatcher.h"
0012 #include "FWCore/Utilities/interface/ESGetToken.h"
0013 #include "FWCore/Utilities/interface/ESInputTag.h"
0014
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016
0017 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0018 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0019 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h"
0020
0021 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0022 #include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
0023
0024 #include "RecoPPS/ProtonReconstruction/interface/ProtonReconstructionAlgorithm.h"
0025
0026 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0027 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0028 #include "CondTools/RunInfo/interface/LHCInfoCombined.h"
0029 #include "CondFormats/RunInfo/interface/LHCInfoPerLS.h"
0030 #include "CondFormats/DataRecord/interface/LHCInfoPerLSRcd.h"
0031 #include "CondFormats/RunInfo/interface/LHCInfoPerFill.h"
0032 #include "CondFormats/DataRecord/interface/LHCInfoPerFillRcd.h"
0033
0034 #include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
0035 #include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"
0036
0037 #include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
0038 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0039
0040 #include "CondFormats/DataRecord/interface/PPSAssociationCutsRcd.h"
0041 #include "CondFormats/PPSObjects/interface/PPSAssociationCuts.h"
0042
0043 //----------------------------------------------------------------------------------------------------
0044
0045 class CTPPSProtonProducer : public edm::stream::EDProducer<> {
0046 public:
0047   explicit CTPPSProtonProducer(const edm::ParameterSet &);
0048   ~CTPPSProtonProducer() override = default;
0049
0050   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0051
0052 private:
0053   void produce(edm::Event &, const edm::EventSetup &) override;
0054
0055   edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tracksToken_;
0056
0057   bool pixelDiscardBXShiftedTracks_;
0058
0059   const std::string lhcInfoPerLSLabel_;
0060   const std::string lhcInfoPerFillLabel_;
0061   const std::string lhcInfoLabel_;
0062   const std::string opticsLabel_;
0063   const std::string ppsAssociationCutsLabel_;
0064   const bool useNewLHCInfo_;
0065
0066   unsigned int verbosity_;
0067
0068   bool doSingleRPReconstruction_;
0069   bool doMultiRPReconstruction_;
0070
0071   std::string singleRPReconstructionLabel_;
0072   std::string multiRPReconstructionLabel_;
0073
0074   double localAngleXMin_, localAngleXMax_, localAngleYMin_, localAngleYMax_;
0075
0076   unsigned int max_n_timing_tracks_;
0077   double default_time_;
0078
0079   ProtonReconstructionAlgorithm algorithm_;
0080
0081   bool opticsValid_;
0082   edm::ESWatcher<CTPPSInterpolatedOpticsRcd> opticsWatcher_;
0083
0084   const edm::ESGetToken<LHCInfoPerLS, LHCInfoPerLSRcd> lhcInfoPerLSToken_;
0085   const edm::ESGetToken<LHCInfoPerFill, LHCInfoPerFillRcd> lhcInfoPerFillToken_;
0086   const edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoToken_;
0087   const edm::ESGetToken<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd> opticalFunctionsToken_;
0088   const edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> geometryToken_;
0089   const edm::ESGetToken<PPSAssociationCuts, PPSAssociationCutsRcd> ppsAssociationCutsToken_;
0090 };
0091
0092 //----------------------------------------------------------------------------------------------------
0093
0094 CTPPSProtonProducer::CTPPSProtonProducer(const edm::ParameterSet &iConfig)
0095     : tracksToken_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagLocalTrackLite"))),
0096
0097       pixelDiscardBXShiftedTracks_(iConfig.getParameter<bool>("pixelDiscardBXShiftedTracks")),
0098
0099       lhcInfoPerLSLabel_(iConfig.getParameter<std::string>("lhcInfoPerLSLabel")),
0100       lhcInfoPerFillLabel_(iConfig.getParameter<std::string>("lhcInfoPerFillLabel")),
0101       lhcInfoLabel_(iConfig.getParameter<std::string>("lhcInfoLabel")),
0102       opticsLabel_(iConfig.getParameter<std::string>("opticsLabel")),
0103       ppsAssociationCutsLabel_(iConfig.getParameter<std::string>("ppsAssociationCutsLabel")),
0104       useNewLHCInfo_(iConfig.getParameter<bool>("useNewLHCInfo")),
0105       verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity", 0)),
0106       doSingleRPReconstruction_(iConfig.getParameter<bool>("doSingleRPReconstruction")),
0107       doMultiRPReconstruction_(iConfig.getParameter<bool>("doMultiRPReconstruction")),
0108       singleRPReconstructionLabel_(iConfig.getParameter<std::string>("singleRPReconstructionLabel")),
0109       multiRPReconstructionLabel_(iConfig.getParameter<std::string>("multiRPReconstructionLabel")),
0110
0111       localAngleXMin_(iConfig.getParameter<double>("localAngleXMin")),
0112       localAngleXMax_(iConfig.getParameter<double>("localAngleXMax")),
0113       localAngleYMin_(iConfig.getParameter<double>("localAngleYMin")),
0114       localAngleYMax_(iConfig.getParameter<double>("localAngleYMax")),
0115
0116       max_n_timing_tracks_(iConfig.getParameter<unsigned int>("max_n_timing_tracks")),
0117       default_time_(iConfig.getParameter<double>("default_time")),
0118
0119       algorithm_(iConfig.getParameter<bool>("fitVtxY"),
0120                  iConfig.getParameter<bool>("useImprovedInitialEstimate"),
0121                  iConfig.getParameter<std::string>("multiRPAlgorithm"),
0122                  verbosity_),
0123       opticsValid_(false),
0124       lhcInfoPerLSToken_(esConsumes<LHCInfoPerLS, LHCInfoPerLSRcd>(edm::ESInputTag("", lhcInfoPerLSLabel_))),
0125       lhcInfoPerFillToken_(esConsumes<LHCInfoPerFill, LHCInfoPerFillRcd>(edm::ESInputTag("", lhcInfoPerFillLabel_))),
0126       lhcInfoToken_(esConsumes<LHCInfo, LHCInfoRcd>(edm::ESInputTag("", lhcInfoLabel_))),
0127       opticalFunctionsToken_(esConsumes<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd>(
0128           edm::ESInputTag("", opticsLabel_))),
0129       geometryToken_(esConsumes<CTPPSGeometry, VeryForwardRealGeometryRecord>()),
0130       ppsAssociationCutsToken_(
0131           esConsumes<PPSAssociationCuts, PPSAssociationCutsRcd>(edm::ESInputTag("", ppsAssociationCutsLabel_))) {
0132   if (doSingleRPReconstruction_)
0133     produces<reco::ForwardProtonCollection>(singleRPReconstructionLabel_);
0134
0135   if (doMultiRPReconstruction_)
0136     produces<reco::ForwardProtonCollection>(multiRPReconstructionLabel_);
0137 }
0138
0139 //----------------------------------------------------------------------------------------------------
0140
0141 void CTPPSProtonProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0142   edm::ParameterSetDescription desc;
0143
0144   desc.add<edm::InputTag>("tagLocalTrackLite", edm::InputTag("ctppsLocalTrackLiteProducer"))
0145       ->setComment("specification of the input lite-track collection");
0146
0147   desc.add<bool>("pixelDiscardBXShiftedTracks", false)
0148       ->setComment("whether to discard pixel tracks built from BX-shifted planes");
0149
0150   desc.add<std::string>("lhcInfoPerFillLabel", "")->setComment("label of the LHCInfoPerFill record");
0151   desc.add<std::string>("lhcInfoPerLSLabel", "")->setComment("label of the LHCInfoPerLS record");
0152   desc.add<std::string>("lhcInfoLabel", "")->setComment("label of the LHCInfo record");
0153   desc.add<bool>("useNewLHCInfo", false)->setComment("flag whether to use new LHCInfoPer* records or old LHCInfo");
0154
0155   desc.add<std::string>("opticsLabel", "")->setComment("label of the optics record");
0156   desc.add<std::string>("ppsAssociationCutsLabel", "")->setComment("label of the association cuts record");
0157
0158   desc.addUntracked<unsigned int>("verbosity", 0)->setComment("verbosity level");
0159
0160   desc.add<bool>("doSingleRPReconstruction", true)
0161       ->setComment("flag whether to apply single-RP reconstruction strategy");
0162
0163   desc.add<bool>("doMultiRPReconstruction", true)->setComment("flag whether to apply multi-RP reconstruction strategy");
0164
0165   desc.add<std::string>("singleRPReconstructionLabel", "singleRP")
0166       ->setComment("output label for single-RP reconstruction products");
0167
0168   desc.add<std::string>("multiRPReconstructionLabel", "multiRP")
0169       ->setComment("output label for multi-RP reconstruction products");
0170
0171   desc.add<double>("localAngleXMin", -0.03)->setComment("minimal accepted value of local horizontal angle (rad)");
0172   desc.add<double>("localAngleXMax", +0.03)->setComment("maximal accepted value of local horizontal angle (rad)");
0173   desc.add<double>("localAngleYMin", -0.04)->setComment("minimal accepted value of local vertical angle (rad)");
0174   desc.add<double>("localAngleYMax", +0.04)->setComment("maximal accepted value of local vertical angle (rad)");
0175
0176   std::vector<edm::ParameterSet> config;
0177
0178   desc.add<unsigned int>("max_n_timing_tracks", 5)->setComment("maximum number of timing tracks per RP");
0179
0180   desc.add<double>("default_time", 0.)->setComment("proton time to be used when no timing information available");
0181
0182   desc.add<bool>("fitVtxY", true)
0183       ->setComment("for multi-RP reconstruction, flag whether y* should be free fit parameter");
0184
0185   desc.add<bool>("useImprovedInitialEstimate", true)
0186       ->setComment(
0187           "for multi-RP reconstruction, flag whether a quadratic estimate of the initial point should be used");
0188
0189   desc.add<std::string>("multiRPAlgorithm", "chi2")
0190       ->setComment("algorithm for multi-RP reco, options include chi2, newton, anal-iter");
0191
0192   descriptions.add("ctppsProtonsDefault", desc);
0193 }
0194
0195 //----------------------------------------------------------------------------------------------------
0196
0197 void CTPPSProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0198   // get input
0199   edm::Handle<CTPPSLocalTrackLiteCollection> hTracks;
0200   iEvent.getByToken(tracksToken_, hTracks);
0201
0202   // book output
0203   std::unique_ptr<reco::ForwardProtonCollection> pOutSingleRP(new reco::ForwardProtonCollection);
0204   std::unique_ptr<reco::ForwardProtonCollection> pOutMultiRP(new reco::ForwardProtonCollection);
0205
0206   // continue only if there is something to process
0207   // NB: this avoids loading (possibly non-existing) conditions in workflows without proton data
0208   if (!hTracks->empty()) {
0209     // get conditions
0210     LHCInfoCombined lhcInfoCombined(iSetup, lhcInfoPerLSToken_, lhcInfoPerFillToken_, lhcInfoToken_, useNewLHCInfo_);
0211
0212     edm::ESHandle<LHCInterpolatedOpticalFunctionsSetCollection> hOpticalFunctions =
0213         iSetup.getHandle(opticalFunctionsToken_);
0214
0215     edm::ESHandle<CTPPSGeometry> hGeometry = iSetup.getHandle(geometryToken_);
0216
0217     edm::ESHandle<PPSAssociationCuts> ppsAssociationCuts = iSetup.getHandle(ppsAssociationCutsToken_);
0218
0219     // re-initialise algorithm upon crossing-angle change
0220     if (opticsWatcher_.check(iSetup)) {
0221       if (hOpticalFunctions->empty()) {
0222         edm::LogInfo("CTPPSProtonProducer") << "No optical functions available, reconstruction disabled.";
0223         algorithm_.release();
0224         opticsValid_ = false;
0225       } else {
0226         algorithm_.init(*hOpticalFunctions);
0227         opticsValid_ = true;
0228       }
0229     }
0230
0231     // do reconstruction only if optics is valid
0232     if (opticsValid_) {
0233       // prepare log
0234       std::ostringstream ssLog;
0235       if (verbosity_)
0236         ssLog << "* input tracks:" << std::endl;
0237
0238       // select tracks with small local angles, split them by LHC sector and tracker/timing RPs
0239       std::map<unsigned int, std::vector<unsigned int>> trackingSelection, timingSelection;
0240
0241       for (unsigned int idx = 0; idx < hTracks->size(); ++idx) {
0242         const auto &tr = hTracks->at(idx);
0243
0244         if (tr.tx() < localAngleXMin_ || tr.tx() > localAngleXMax_ || tr.ty() < localAngleYMin_ ||
0245             tr.ty() > localAngleYMax_)
0246           continue;
0247
0248         if (pixelDiscardBXShiftedTracks_) {
0249           if (tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes ||
0250               tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes)
0251             continue;
0252         }
0253
0254         const CTPPSDetId rpId(tr.rpId());
0255
0256         if (verbosity_)
0257           ssLog << "\t"
0258                 << "[" << idx << "] " << tr.rpId() << " (" << (rpId.arm() * 100 + rpId.station() * 10 + rpId.rp())
0259                 << "): "
0260                 << "x=" << tr.x() << " +- " << tr.xUnc() << " mm, "
0261                 << "y=" << tr.y() << " +- " << tr.yUnc() << " mm" << std::endl;
0262
0263         const bool trackerRP =
0264             (rpId.subdetId() == CTPPSDetId::sdTrackingStrip || rpId.subdetId() == CTPPSDetId::sdTrackingPixel);
0265
0266         if (trackerRP)
0267           trackingSelection[rpId.arm()].push_back(idx);
0268         else
0269           timingSelection[rpId.arm()].push_back(idx);
0270       }
0271
0272       // process each arm
0273       for (const auto &arm_it : trackingSelection) {
0274         const auto &indices = arm_it.second;
0275
0276         const auto &ac = ppsAssociationCuts->getAssociationCuts(arm_it.first);
0277
0278         // do single-RP reco if needed
0279         std::map<unsigned int, reco::ForwardProton> singleRPResultsIndexed;
0280         if (doSingleRPReconstruction_ || ac.isApplied(ac.qXi) || ac.isApplied(ac.qThetaY)) {
0281           for (const auto &idx : indices) {
0282             if (verbosity_)
0283               ssLog << std::endl << "* reconstruction from track " << idx << std::endl;
0284
0285             singleRPResultsIndexed[idx] =
0286                 algorithm_.reconstructFromSingleRP(CTPPSLocalTrackLiteRef(hTracks, idx), lhcInfoCombined.energy, ssLog);
0287           }
0288         }
0289
0290         // check that exactly two tracking RPs are involved
0291         //    - 1 is insufficient for multi-RP reconstruction
0292         //    - PPS did not use more than 2 tracking RPs per arm -> algorithms are tuned to this
0293         std::set<unsigned int> rpIds;
0294         for (const auto &idx : indices)
0295           rpIds.insert(hTracks->at(idx).rpId());
0296
0297         // do multi-RP reco if chosen
0298         if (doMultiRPReconstruction_ && rpIds.size() == 2) {
0299           // find matching track pairs from different tracking RPs, ordered: i=near, j=far RP
0300           std::vector<std::pair<unsigned int, unsigned int>> idx_pairs;
0301           std::map<unsigned int, unsigned int> idx_pair_multiplicity;
0302           for (const auto &i : indices) {
0303             for (const auto &j : indices) {
0304               const auto &tr_i = hTracks->at(i);
0305               const auto &tr_j = hTracks->at(j);
0306
0307               const double z_i = hGeometry->rpTranslation(tr_i.rpId()).z();
0308               const double z_j = hGeometry->rpTranslation(tr_j.rpId()).z();
0309
0310               const auto &pr_i = singleRPResultsIndexed[i];
0311               const auto &pr_j = singleRPResultsIndexed[j];
0312
0313               if (tr_i.rpId() == tr_j.rpId())
0314                 continue;
0315
0316               if (std::abs(z_i) >= std::abs(z_j))
0317                 continue;
0318
0319               bool matching = true;
0320
0321               if (!ac.isSatisfied(ac.qX, tr_i.x(), tr_i.y(), lhcInfoCombined.crossingAngle(), tr_i.x() - tr_j.x()))
0322                 matching = false;
0323               else if (!ac.isSatisfied(ac.qY, tr_i.x(), tr_i.y(), lhcInfoCombined.crossingAngle(), tr_i.y() - tr_j.y()))
0324                 matching = false;
0325               else if (!ac.isSatisfied(
0326                            ac.qXi, tr_i.x(), tr_i.y(), lhcInfoCombined.crossingAngle(), pr_i.xi() - pr_j.xi()))
0327                 matching = false;
0328               else if (!ac.isSatisfied(ac.qThetaY,
0329                                        tr_i.x(),
0330                                        tr_i.y(),
0331                                        lhcInfoCombined.crossingAngle(),
0332                                        pr_i.thetaY() - pr_j.thetaY()))
0333                 matching = false;
0334
0335               if (!matching)
0336                 continue;
0337
0338               idx_pairs.emplace_back(i, j);
0339               idx_pair_multiplicity[i]++;
0340               idx_pair_multiplicity[j]++;
0341             }
0342           }
0343
0344           // evaluate track multiplicity in each timing RP
0345           std::map<unsigned int, unsigned int> timing_RP_track_multiplicity;
0346           for (const auto &ti : timingSelection[arm_it.first]) {
0347             const auto &tr = hTracks->at(ti);
0348             timing_RP_track_multiplicity[tr.rpId()]++;
0349           }
0350
0351           // associate tracking-RP pairs with timing-RP tracks
0352           std::map<unsigned int, std::vector<unsigned int>> matched_timing_track_indices;
0353           std::map<unsigned int, unsigned int> matched_timing_track_multiplicity;
0354           for (unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
0355             const auto &i = idx_pairs[pr_idx].first;
0356             const auto &j = idx_pairs[pr_idx].second;
0357
0358             // skip non-unique associations
0359             if (idx_pair_multiplicity[i] > 1 || idx_pair_multiplicity[j] > 1)
0360               continue;
0361
0362             const auto &tr_i = hTracks->at(i);
0363             const auto &tr_j = hTracks->at(j);
0364
0365             const double z_i = hGeometry->rpTranslation(tr_i.rpId()).z();
0366             const double z_j = hGeometry->rpTranslation(tr_j.rpId()).z();
0367
0368             for (const auto &ti : timingSelection[arm_it.first]) {
0369               const auto &tr_ti = hTracks->at(ti);
0370
0371               // skip if timing RP saturated (high track multiplicity)
0372               if (timing_RP_track_multiplicity[tr_ti.rpId()] > max_n_timing_tracks_)
0373                 continue;
0374
0375               // interpolation from tracking RPs
0376               const double z_ti = hGeometry->rpTranslation(tr_ti.rpId()).z();
0377               const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j);
0378               const double x_inter = f_i * tr_i.x() + f_j * tr_j.x();
0379               const double x_inter_unc_sq =
0380                   f_i * f_i * tr_i.xUnc() * tr_i.xUnc() + f_j * f_j * tr_j.xUnc() * tr_j.xUnc();
0381
0382               const double de_x = tr_ti.x() - x_inter;
0383               const double de_x_unc = sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_sq);
0384               const double r = (de_x_unc > 0.) ? de_x / de_x_unc : 1E100;
0385
0386               const bool matching = (ac.getTiTrMin() < r && r < ac.getTiTrMax());
0387
0388               if (verbosity_)
0389                 ssLog << "ti=" << ti << ", i=" << i << ", j=" << j << " | z_ti=" << z_ti << ", z_i=" << z_i
0390                       << ", z_j=" << z_j << " | x_ti=" << tr_ti.x() << ", x_inter=" << x_inter << ", de_x=" << de_x
0391                       << ", de_x_unc=" << de_x_unc << ", matching=" << matching << std::endl;
0392
0393               if (!matching)
0394                 continue;
0395
0396               matched_timing_track_indices[pr_idx].push_back(ti);
0397               matched_timing_track_multiplicity[ti]++;
0398             }
0399           }
0400
0401           // process associated tracks
0402           for (unsigned int pr_idx = 0; pr_idx < idx_pairs.size(); ++pr_idx) {
0403             const auto &i = idx_pairs[pr_idx].first;
0404             const auto &j = idx_pairs[pr_idx].second;
0405
0406             // skip non-unique associations of tracking-RP tracks
0407             if (idx_pair_multiplicity[i] > 1 || idx_pair_multiplicity[j] > 1)
0408               continue;
0409
0410             if (verbosity_)
0411               ssLog << std::endl
0412                     << "* reconstruction from tracking-RP tracks: " << i << ", " << j << " and timing-RP tracks: ";
0413
0414             // buffer contributing tracks
0415             CTPPSLocalTrackLiteRefVector sel_tracks;
0416             sel_tracks.push_back(CTPPSLocalTrackLiteRef(hTracks, i));
0417             sel_tracks.push_back(CTPPSLocalTrackLiteRef(hTracks, j));
0418
0419             CTPPSLocalTrackLiteRefVector sel_track_for_kin_reco = sel_tracks;
0420
0421             // process timing-RP data
0422             double sw = 0., swt = 0.;
0423             for (const auto &ti : matched_timing_track_indices[pr_idx]) {
0424               // skip non-unique associations of timing-RP tracks
0425               if (matched_timing_track_multiplicity[ti] > 1)
0426                 continue;
0427
0428               sel_tracks.push_back(CTPPSLocalTrackLiteRef(hTracks, ti));
0429
0430               if (verbosity_)
0431                 ssLog << ti << ", ";
0432
0433               const auto &tr = hTracks->at(ti);
0434               const double t_unc = tr.timeUnc();
0435               const double w = (t_unc > 0.) ? 1. / t_unc / t_unc : 1.;
0436               sw += w;
0437               swt += w * tr.time();
0438             }
0439
0440             float time = default_time_, time_unc = 0.;
0441             if (sw > 0.) {
0442               time = swt / sw;
0443               time_unc = 1. / sqrt(sw);
0444             }
0445
0446             if (verbosity_)
0447               ssLog << std::endl << "    time = " << time << " +- " << time_unc << std::endl;
0448
0449             // process tracking-RP data
0450             reco::ForwardProton proton =
0451                 algorithm_.reconstructFromMultiRP(sel_track_for_kin_reco, lhcInfoCombined.energy, ssLog);
0452
0453             // save combined output
0454             proton.setContributingLocalTracks(sel_tracks);
0455             proton.setTime(time);
0456             proton.setTimeError(time_unc);
0457
0458             pOutMultiRP->emplace_back(proton);
0459           }
0460         }
0461
0462         // save single-RP results (un-indexed)
0463         for (const auto &p : singleRPResultsIndexed)
0464           pOutSingleRP->emplace_back(p.second);
0465       }
0466
0467       // dump log
0468       if (verbosity_)
0469         edm::LogInfo("CTPPSProtonProducer") << ssLog.str();
0470     }
0471   }
0472
0473   // save output
0474   if (doSingleRPReconstruction_)
0475     iEvent.put(std::move(pOutSingleRP), singleRPReconstructionLabel_);
0476
0477   if (doMultiRPReconstruction_)
0478     iEvent.put(std::move(pOutMultiRP), multiRPReconstructionLabel_);
0479 }
0480
0481 //----------------------------------------------------------------------------------------------------
0482
0483 DEFINE_FWK_MODULE(CTPPSProtonProducer);