Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-07-02 00:54:00

0001 // Author: Felice Pantaleo, Wahid Redjeb, Aurora Perego (CERN) - felice.pantaleo@cern.ch, wahid.redjeb@cern.ch, aurora.perego@cern.ch
0002 // Date: 12/2023
0003 #include <memory>  // unique_ptr
0004 #include "CommonTools/RecoAlgos/interface/MultiVectorManager.h"
0005 #include "FWCore/Framework/interface/stream/EDProducer.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0008 #include "FWCore/ParameterSet/interface/PluginDescription.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Utilities/interface/ESGetToken.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0015 #include "FWCore/Framework/interface/ConsumesCollector.h"
0016 #include "DataFormats/Common/interface/OrphanHandle.h"
0017 
0018 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0019 #include "DataFormats/HGCalReco/interface/Common.h"
0020 #include "DataFormats/HGCalReco/interface/MtdHostCollection.h"
0021 #include "DataFormats/HGCalReco/interface/TICLLayerTile.h"
0022 #include "DataFormats/HGCalReco/interface/Trackster.h"
0023 #include "DataFormats/TrackReco/interface/Track.h"
0024 #include "DataFormats/MuonReco/interface/Muon.h"
0025 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0026 #include "DataFormats/HGCalReco/interface/TICLCandidate.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0029 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0030 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0031 
0032 #include "RecoHGCal/TICL/interface/TICLInterpretationAlgoBase.h"
0033 #include "RecoHGCal/TICL/plugins/TICLInterpretationPluginFactory.h"
0034 #include "RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.h"
0035 
0036 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
0037 
0038 #include "RecoHGCal/TICL/interface/GlobalCache.h"
0039 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0040 
0041 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0042 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0043 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0044 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToBeamLine.h"
0045 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0046 
0047 #include "MagneticField/Engine/interface/MagneticField.h"
0048 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0049 
0050 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0051 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0052 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0053 
0054 #include "TrackstersPCA.h"
0055 
0056 using namespace ticl;
0057 
0058 class TICLCandidateProducer : public edm::stream::EDProducer<> {
0059 public:
0060   explicit TICLCandidateProducer(const edm::ParameterSet &ps);
0061   ~TICLCandidateProducer() override{};
0062   void produce(edm::Event &, const edm::EventSetup &) override;
0063   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0064 
0065   void beginRun(edm::Run const &iEvent, edm::EventSetup const &es) override;
0066 
0067 private:
0068   void dumpCandidate(const TICLCandidate &) const;
0069 
0070   template <typename F>
0071   void assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates,
0072                               edm::Handle<std::vector<reco::Track>> track_h,
0073                               MtdHostCollection::ConstView &inputTimingView,
0074                               F func) const;
0075 
0076   std::unique_ptr<TICLInterpretationAlgoBase<reco::Track>> generalInterpretationAlgo_;
0077   std::vector<edm::EDGetTokenT<std::vector<Trackster>>> egamma_tracksters_tokens_;
0078   std::vector<edm::EDGetTokenT<std::vector<std::vector<unsigned>>>> egamma_tracksterlinks_tokens_;
0079 
0080   std::vector<edm::EDGetTokenT<std::vector<Trackster>>> general_tracksters_tokens_;
0081   std::vector<edm::EDGetTokenT<std::vector<std::vector<unsigned>>>> general_tracksterlinks_tokens_;
0082 
0083   const edm::EDGetTokenT<std::vector<reco::CaloCluster>> clusters_token_;
0084   const edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> clustersTime_token_;
0085 
0086   std::vector<edm::EDGetTokenT<std::vector<float>>> original_masks_tokens_;
0087 
0088   const edm::EDGetTokenT<std::vector<reco::Track>> tracks_token_;
0089   edm::EDGetTokenT<MtdHostCollection> inputTimingToken_;
0090 
0091   const edm::EDGetTokenT<std::vector<reco::Muon>> muons_token_;
0092   const bool useMTDTiming_;
0093   const bool useTimingAverage_;
0094   const float timingQualityThreshold_;
0095   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometry_token_;
0096 
0097   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bfield_token_;
0098   const std::string detector_;
0099   const std::string propName_;
0100   const edm::ESGetToken<Propagator, TrackingComponentsRecord> propagator_token_;
0101   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> trackingGeometry_token_;
0102 
0103   std::once_flag initializeGeometry_;
0104   const HGCalDDDConstants *hgcons_;
0105   hgcal::RecHitTools rhtools_;
0106   const float tkEnergyCut_ = 2.0f;
0107   const StringCutObjectSelector<reco::Track> cutTk_;
0108   edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> hdc_token_;
0109   edm::ESHandle<MagneticField> bfield_;
0110   edm::ESHandle<Propagator> propagator_;
0111   edm::ESHandle<GlobalTrackingGeometry> trackingGeometry_;
0112   static constexpr float c_light_ = CLHEP::c_light * CLHEP::ns / CLHEP::cm;
0113   static constexpr float timeRes = 0.02f;
0114 };
0115 
0116 TICLCandidateProducer::TICLCandidateProducer(const edm::ParameterSet &ps)
0117     : clusters_token_(consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("layer_clusters"))),
0118       clustersTime_token_(
0119           consumes<edm::ValueMap<std::pair<float, float>>>(ps.getParameter<edm::InputTag>("layer_clustersTime"))),
0120       tracks_token_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("tracks"))),
0121       inputTimingToken_(consumes<MtdHostCollection>(ps.getParameter<edm::InputTag>("timingSoA"))),
0122       muons_token_(consumes<std::vector<reco::Muon>>(ps.getParameter<edm::InputTag>("muons"))),
0123       useMTDTiming_(ps.getParameter<bool>("useMTDTiming")),
0124       useTimingAverage_(ps.getParameter<bool>("useTimingAverage")),
0125       timingQualityThreshold_(ps.getParameter<double>("timingQualityThreshold")),
0126       geometry_token_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0127       bfield_token_(esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>()),
0128       detector_(ps.getParameter<std::string>("detector")),
0129       propName_(ps.getParameter<std::string>("propagator")),
0130       propagator_token_(
0131           esConsumes<Propagator, TrackingComponentsRecord, edm::Transition::BeginRun>(edm::ESInputTag("", propName_))),
0132       trackingGeometry_token_(
0133           esConsumes<GlobalTrackingGeometry, GlobalTrackingGeometryRecord, edm::Transition::BeginRun>()),
0134       cutTk_(ps.getParameter<std::string>("cutTk")) {
0135   // These are the CLUE3DEM Tracksters put in the event by the TracksterLinksProducer with the superclustering algorithm
0136   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("egamma_tracksters_collections")) {
0137     egamma_tracksters_tokens_.emplace_back(consumes<std::vector<Trackster>>(tag));
0138   }
0139 
0140   // These are the links put in the event by the TracksterLinksProducer with the superclustering algorithm
0141   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("egamma_tracksterlinks_collections")) {
0142     egamma_tracksterlinks_tokens_.emplace_back(consumes<std::vector<std::vector<unsigned int>>>(tag));
0143   }
0144 
0145   //make sure that the number of tracksters collections and tracksterlinks collections is the same
0146   assert(egamma_tracksters_tokens_.size() == egamma_tracksterlinks_tokens_.size());
0147 
0148   // Loop over the edm::VInputTag and append the token to general_tracksters_tokens_
0149   // These, instead, are the tracksters already merged by the TrackstersLinksProducer
0150   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("general_tracksters_collections")) {
0151     general_tracksters_tokens_.emplace_back(consumes<std::vector<Trackster>>(tag));
0152   }
0153 
0154   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("general_tracksterlinks_collections")) {
0155     general_tracksterlinks_tokens_.emplace_back(consumes<std::vector<std::vector<unsigned int>>>(tag));
0156   }
0157 
0158   //make sure that the number of tracksters collections and tracksterlinks collections is the same
0159   assert(general_tracksters_tokens_.size() == general_tracksterlinks_tokens_.size());
0160 
0161   //Loop over the edm::VInputTag of masks and append the token to original_masks_tokens_
0162   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("original_masks")) {
0163     original_masks_tokens_.emplace_back(consumes<std::vector<float>>(tag));
0164   }
0165 
0166   std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive";
0167   hdc_token_ =
0168       esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag("", detectorName_));
0169   if (useMTDTiming_) {
0170     inputTimingToken_ = consumes<MtdHostCollection>(ps.getParameter<edm::InputTag>("timingSoA"));
0171   }
0172 
0173   produces<std::vector<TICLCandidate>>();
0174 
0175   // New trackster collection after linking
0176   produces<std::vector<Trackster>>();
0177 
0178   auto interpretationPSet = ps.getParameter<edm::ParameterSet>("interpretationDescPSet");
0179   auto algoType = interpretationPSet.getParameter<std::string>("type");
0180   generalInterpretationAlgo_ =
0181       TICLGeneralInterpretationPluginFactory::get()->create(algoType, interpretationPSet, consumesCollector());
0182 }
0183 
0184 void TICLCandidateProducer::beginRun(edm::Run const &iEvent, edm::EventSetup const &es) {
0185   edm::ESHandle<HGCalDDDConstants> hdc = es.getHandle(hdc_token_);
0186   hgcons_ = hdc.product();
0187 
0188   edm::ESHandle<CaloGeometry> geom = es.getHandle(geometry_token_);
0189   rhtools_.setGeometry(*geom);
0190 
0191   bfield_ = es.getHandle(bfield_token_);
0192   propagator_ = es.getHandle(propagator_token_);
0193   generalInterpretationAlgo_->initialize(hgcons_, rhtools_, bfield_, propagator_);
0194 
0195   trackingGeometry_ = es.getHandle(trackingGeometry_token_);
0196 };
0197 
0198 void filterTracks(edm::Handle<std::vector<reco::Track>> tkH,
0199                   const edm::Handle<std::vector<reco::Muon>> &muons_h,
0200                   const StringCutObjectSelector<reco::Track> cutTk_,
0201                   const float tkEnergyCut_,
0202                   std::vector<bool> &maskTracks) {
0203   auto const &tracks = *tkH;
0204   for (unsigned i = 0; i < tracks.size(); ++i) {
0205     const auto &tk = tracks[i];
0206     reco::TrackRef trackref = reco::TrackRef(tkH, i);
0207 
0208     // veto tracks associated to muons
0209     int muId = PFMuonAlgo::muAssocToTrack(trackref, *muons_h);
0210     const reco::MuonRef muonref = reco::MuonRef(muons_h, muId);
0211 
0212     if (!cutTk_((tk)) or (muId != -1 and PFMuonAlgo::isMuon(muonref) and not(*muons_h)[muId].isTrackerMuon())) {
0213       maskTracks[i] = false;
0214       continue;
0215     }
0216 
0217     // don't consider tracks below 2 GeV for linking
0218     if (std::sqrt(tk.p() * tk.p() + ticl::mpion2) < tkEnergyCut_) {
0219       maskTracks[i] = false;
0220       continue;
0221     }
0222 
0223     // record tracks that can be used to make a ticlcandidate
0224     maskTracks[i] = true;
0225   }
0226 }
0227 
0228 void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) {
0229   auto resultTracksters = std::make_unique<std::vector<Trackster>>();
0230   auto resultTrackstersMerged = std::make_unique<std::vector<Trackster>>();
0231   auto linkedResultTracksters = std::make_unique<std::vector<std::vector<unsigned int>>>();
0232 
0233   const auto &layerClusters = evt.get(clusters_token_);
0234   const auto &layerClustersTimes = evt.get(clustersTime_token_);
0235   edm::Handle<reco::MuonCollection> muons_h;
0236   evt.getByToken(muons_token_, muons_h);
0237 
0238   edm::Handle<std::vector<reco::Track>> tracks_h;
0239   evt.getByToken(tracks_token_, tracks_h);
0240   const auto &tracks = *tracks_h;
0241 
0242   edm::Handle<MtdHostCollection> inputTiming_h;
0243   MtdHostCollection::ConstView inputTimingView;
0244   if (useMTDTiming_) {
0245     evt.getByToken(inputTimingToken_, inputTiming_h);
0246     inputTimingView = (*inputTiming_h).const_view();
0247   }
0248 
0249   auto const bFieldProd = bfield_.product();
0250   const Propagator *propagator = propagator_.product();
0251 
0252   // loop over the original_masks_tokens_ and get the original masks collections and multiply them
0253   // to get the global mask
0254   std::vector<float> original_global_mask(layerClusters.size(), 1.f);
0255   for (unsigned int i = 0; i < original_masks_tokens_.size(); ++i) {
0256     const auto &tmp_mask = evt.get(original_masks_tokens_[i]);
0257     for (unsigned int j = 0; j < tmp_mask.size(); ++j) {
0258       original_global_mask[j] *= tmp_mask[j];
0259     }
0260   }
0261 
0262   auto resultMask = std::make_unique<std::vector<float>>(original_global_mask);
0263 
0264   std::vector<edm::Handle<std::vector<Trackster>>> general_tracksters_h(general_tracksters_tokens_.size());
0265   MultiVectorManager<Trackster> generalTrackstersManager;
0266   for (unsigned int i = 0; i < general_tracksters_tokens_.size(); ++i) {
0267     evt.getByToken(general_tracksters_tokens_[i], general_tracksters_h[i]);
0268     //Fill MultiVectorManager
0269     generalTrackstersManager.addVector(*general_tracksters_h[i]);
0270   }
0271   //now get the general_tracksterlinks_tokens_
0272   std::vector<edm::Handle<std::vector<std::vector<unsigned>>>> general_tracksterlinks_h(
0273       general_tracksterlinks_tokens_.size());
0274   std::vector<std::vector<unsigned>> generalTracksterLinksGlobalId;
0275   for (unsigned int i = 0; i < general_tracksterlinks_tokens_.size(); ++i) {
0276     evt.getByToken(general_tracksterlinks_tokens_[i], general_tracksterlinks_h[i]);
0277     for (unsigned int j = 0; j < general_tracksterlinks_h[i]->size(); ++j) {
0278       generalTracksterLinksGlobalId.emplace_back();
0279       auto &links_vector = generalTracksterLinksGlobalId.back();
0280       links_vector.resize((*general_tracksterlinks_h[i])[j].size());
0281       for (unsigned int k = 0; k < links_vector.size(); ++k) {
0282         links_vector[k] = generalTrackstersManager.getGlobalIndex(i, (*general_tracksterlinks_h[i])[j][k]);
0283       }
0284     }
0285   }
0286 
0287   std::vector<bool> maskTracks;
0288   maskTracks.resize(tracks.size());
0289   filterTracks(tracks_h, muons_h, cutTk_, tkEnergyCut_, maskTracks);
0290 
0291   const typename TICLInterpretationAlgoBase<reco::Track>::Inputs input(evt,
0292                                                                        es,
0293                                                                        layerClusters,
0294                                                                        layerClustersTimes,
0295                                                                        generalTrackstersManager,
0296                                                                        generalTracksterLinksGlobalId,
0297                                                                        tracks_h,
0298                                                                        maskTracks);
0299 
0300   auto resultCandidates = std::make_unique<std::vector<TICLCandidate>>();
0301   std::vector<int> trackstersInTrackIndices(tracks.size(), -1);
0302 
0303   //TODO
0304   //egammaInterpretationAlg_->makecandidates(inputGSF, inputTiming, *resultTrackstersMerged, trackstersInGSFTrackIndices)
0305   // mask generalTracks associated to GSFTrack linked in egammaInterpretationAlgo_
0306 
0307   generalInterpretationAlgo_->makeCandidates(input, inputTiming_h, *resultTracksters, trackstersInTrackIndices);
0308 
0309   assignPCAtoTracksters(
0310       *resultTracksters, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(), true);
0311 
0312   std::vector<bool> maskTracksters(resultTracksters->size(), true);
0313   edm::OrphanHandle<std::vector<Trackster>> resultTracksters_h = evt.put(std::move(resultTracksters));
0314   //create ChargedCandidates
0315   for (size_t iTrack = 0; iTrack < tracks.size(); iTrack++) {
0316     if (maskTracks[iTrack]) {
0317       auto const tracksterId = trackstersInTrackIndices[iTrack];
0318       auto trackPtr = edm::Ptr<reco::Track>(tracks_h, iTrack);
0319       if (tracksterId != -1 and !maskTracksters.empty()) {
0320         auto tracksterPtr = edm::Ptr<Trackster>(resultTracksters_h, tracksterId);
0321         TICLCandidate chargedCandidate(trackPtr, tracksterPtr);
0322         resultCandidates->push_back(chargedCandidate);
0323         maskTracksters[tracksterId] = false;
0324       } else {
0325         //charged candidates track only
0326         edm::Ptr<Trackster> tracksterPtr;
0327         TICLCandidate chargedCandidate(trackPtr, tracksterPtr);
0328         auto trackRef = edm::Ref<reco::TrackCollection>(tracks_h, iTrack);
0329         const int muId = PFMuonAlgo::muAssocToTrack(trackRef, *muons_h);
0330         const reco::MuonRef muonRef = reco::MuonRef(muons_h, muId);
0331         if (muonRef.isNonnull() and muonRef->isGlobalMuon()) {
0332           // create muon candidate
0333           chargedCandidate.setPdgId(13 * trackPtr.get()->charge());
0334         }
0335         resultCandidates->push_back(chargedCandidate);
0336       }
0337     }
0338   }
0339 
0340   //Neutral Candidate
0341   for (size_t iTrackster = 0; iTrackster < maskTracksters.size(); iTrackster++) {
0342     if (maskTracksters[iTrackster]) {
0343       edm::Ptr<Trackster> tracksterPtr(resultTracksters_h, iTrackster);
0344       edm::Ptr<reco::Track> trackPtr;
0345       TICLCandidate neutralCandidate(trackPtr, tracksterPtr);
0346       resultCandidates->push_back(neutralCandidate);
0347     }
0348   }
0349 
0350   auto getPathLength = [&](const reco::Track &track, float zVal) {
0351     const auto &fts_inn = trajectoryStateTransform::innerFreeState(track, bFieldProd);
0352     const auto &fts_out = trajectoryStateTransform::outerFreeState(track, bFieldProd);
0353     const auto &surf_inn = trajectoryStateTransform::innerStateOnSurface(track, *trackingGeometry_, bFieldProd);
0354     const auto &surf_out = trajectoryStateTransform::outerStateOnSurface(track, *trackingGeometry_, bFieldProd);
0355 
0356     Basic3DVector<float> pos(track.referencePoint());
0357     Basic3DVector<float> mom(track.momentum());
0358     FreeTrajectoryState stateAtBeamspot{GlobalPoint(pos), GlobalVector(mom), track.charge(), bFieldProd};
0359 
0360     float pathlength = propagator->propagateWithPath(stateAtBeamspot, surf_inn.surface()).second;
0361 
0362     if (pathlength) {
0363       const auto &t_inn_out = propagator->propagateWithPath(fts_inn, surf_out.surface());
0364 
0365       if (t_inn_out.first.isValid()) {
0366         pathlength += t_inn_out.second;
0367 
0368         std::pair<float, float> rMinMax = hgcons_->rangeR(zVal, true);
0369 
0370         int iSide = int(track.eta() > 0);
0371         float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
0372         const auto &disk = std::make_unique<GeomDet>(
0373             Disk::build(Disk::PositionType(0, 0, zSide),
0374                         Disk::RotationType(),
0375                         SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5))
0376                 .get());
0377         const auto &tsos = propagator->propagateWithPath(fts_out, disk->surface());
0378 
0379         if (tsos.first.isValid()) {
0380           pathlength += tsos.second;
0381           return pathlength;
0382         }
0383       }
0384     }
0385     edm::LogWarning("TICLCandidateProducer")
0386         << "Not able to use the track to compute the path length. A straight line will be used instead.";
0387     return 0.f;
0388   };
0389 
0390   assignTimeToCandidates(*resultCandidates, tracks_h, inputTimingView, getPathLength);
0391 
0392   evt.put(std::move(resultCandidates));
0393 }
0394 
0395 template <typename F>
0396 void TICLCandidateProducer::assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates,
0397                                                    edm::Handle<std::vector<reco::Track>> track_h,
0398                                                    MtdHostCollection::ConstView &inputTimingView,
0399                                                    F func) const {
0400   for (auto &cand : resultCandidates) {
0401     float beta = 1;
0402     float time = 0.f;
0403     float invTimeErr = 0.f;
0404     float timeErr = -1.f;
0405 
0406     for (const auto &tr : cand.tracksters()) {
0407       if (tr->timeError() > 0) {
0408         const auto invTimeESq = pow(tr->timeError(), -2);
0409         const auto x = tr->barycenter().X();
0410         const auto y = tr->barycenter().Y();
0411         const auto z = tr->barycenter().Z();
0412         auto path = std::sqrt(x * x + y * y + z * z);
0413         if (cand.trackPtr().get() != nullptr) {
0414           const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
0415           if (useMTDTiming_ and inputTimingView.timeErr()[trackIndex] > 0) {
0416             const auto xMtd = inputTimingView.posInMTD_x()[trackIndex];
0417             const auto yMtd = inputTimingView.posInMTD_y()[trackIndex];
0418             const auto zMtd = inputTimingView.posInMTD_z()[trackIndex];
0419 
0420             beta = inputTimingView.beta()[trackIndex];
0421             path = std::sqrt((x - xMtd) * (x - xMtd) + (y - yMtd) * (y - yMtd) + (z - zMtd) * (z - zMtd)) +
0422                    inputTimingView.pathLength()[trackIndex];
0423           } else {
0424             float pathLength = func(*(cand.trackPtr().get()), z);
0425             if (pathLength) {
0426               path = pathLength;
0427             }
0428           }
0429         }
0430         time += (tr->time() - path / (beta * c_light_)) * invTimeESq;
0431         invTimeErr += invTimeESq;
0432       }
0433     }
0434     if (invTimeErr > 0) {
0435       time = time / invTimeErr;
0436       // FIXME_ set a liminf of 0.02 ns on the ts error (based on residuals)
0437       timeErr = sqrt(1.f / invTimeErr);
0438       if (timeErr < timeRes)
0439         timeErr = timeRes;
0440       cand.setTime(time, timeErr);
0441     }
0442 
0443     if (useMTDTiming_ and cand.charge()) {
0444       // Check MTD timing availability
0445       const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
0446       const bool assocQuality = inputTimingView.MVAquality()[trackIndex] > timingQualityThreshold_;
0447       if (assocQuality) {
0448         const auto timeHGC = cand.time();
0449         const auto timeEHGC = cand.timeError();
0450         const auto timeMTD = inputTimingView.time0()[trackIndex];
0451         const auto timeEMTD = inputTimingView.time0Err()[trackIndex];
0452 
0453         if (useTimingAverage_ && (timeEMTD > 0 && timeEHGC > 0)) {
0454           // Compute weighted average between HGCAL and MTD timing
0455           const auto invTimeESqHGC = pow(timeEHGC, -2);
0456           const auto invTimeESqMTD = pow(timeEMTD, -2);
0457           timeErr = 1.f / (invTimeESqHGC + invTimeESqMTD);
0458           time = (timeHGC * invTimeESqHGC + timeMTD * invTimeESqMTD) * timeErr;
0459           timeErr = sqrt(timeErr);
0460         } else if (timeEMTD > 0) {
0461           time = timeMTD;
0462           timeErr = timeEMTD;
0463         }
0464       }
0465       cand.setTime(time, timeErr);
0466       cand.setMTDTime(inputTimingView.time()[trackIndex], inputTimingView.timeErr()[trackIndex]);
0467     }
0468   }
0469 }
0470 
0471 void TICLCandidateProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0472   edm::ParameterSetDescription desc;
0473   edm::ParameterSetDescription interpretationDesc;
0474   interpretationDesc.addNode(edm::PluginDescription<TICLGeneralInterpretationPluginFactory>("type", "General", true));
0475   desc.add<edm::ParameterSetDescription>("interpretationDescPSet", interpretationDesc);
0476   desc.add<std::vector<edm::InputTag>>("egamma_tracksters_collections", {edm::InputTag("ticlTracksterLinks")});
0477   desc.add<std::vector<edm::InputTag>>("egamma_tracksterlinks_collections", {edm::InputTag("ticlTracksterLinks")});
0478   desc.add<std::vector<edm::InputTag>>("general_tracksters_collections", {edm::InputTag("ticlTracksterLinks")});
0479   desc.add<std::vector<edm::InputTag>>("general_tracksterlinks_collections", {edm::InputTag("ticlTracksterLinks")});
0480   desc.add<std::vector<edm::InputTag>>("original_masks",
0481                                        {edm::InputTag("hgcalMergeLayerClusters", "InitialLayerClustersMask")});
0482   desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalMergeLayerClusters"));
0483   desc.add<edm::InputTag>("layer_clustersTime", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster"));
0484   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0485   desc.add<edm::InputTag>("timingSoA", edm::InputTag("mtdSoA"));
0486   desc.add<edm::InputTag>("muons", edm::InputTag("muons1stStep"));
0487   desc.add<std::string>("detector", "HGCAL");
0488   desc.add<std::string>("propagator", "PropagatorWithMaterial");
0489   desc.add<bool>("useMTDTiming", true);
0490   desc.add<bool>("useTimingAverage", true);
0491   desc.add<double>("timingQualityThreshold", 0.5);
0492   desc.add<std::string>("cutTk",
0493                         "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
0494                         "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
0495   descriptions.add("ticlCandidateProducer", desc);
0496 }
0497 
0498 DEFINE_FWK_MODULE(TICLCandidateProducer);