Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-23 03:14:18

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