Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:31 

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       muons_token_(consumes<std::vector<reco::Muon>>(ps.getParameter<edm::InputTag>("muons"))),
0122       useMTDTiming_(ps.getParameter<bool>("useMTDTiming")),
0123       useTimingAverage_(ps.getParameter<bool>("useTimingAverage")),
0124       timingQualityThreshold_(ps.getParameter<double>("timingQualityThreshold")),
0125       geometry_token_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0126       bfield_token_(esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>()),
0127       detector_(ps.getParameter<std::string>("detector")),
0128       propName_(ps.getParameter<std::string>("propagator")),
0129       propagator_token_(
0130           esConsumes<Propagator, TrackingComponentsRecord, edm::Transition::BeginRun>(edm::ESInputTag("", propName_))),
0131       trackingGeometry_token_(
0132           esConsumes<GlobalTrackingGeometry, GlobalTrackingGeometryRecord, edm::Transition::BeginRun>()),
0133       cutTk_(ps.getParameter<std::string>("cutTk")) {
0134   // These are the CLUE3DEM Tracksters put in the event by the TracksterLinksProducer with the superclustering algorithm
0135   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("egamma_tracksters_collections")) {
0136     egamma_tracksters_tokens_.emplace_back(consumes<std::vector<Trackster>>(tag));
0137   }
0138 
0139   // These are the links put in the event by the TracksterLinksProducer with the superclustering algorithm
0140   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("egamma_tracksterlinks_collections")) {
0141     egamma_tracksterlinks_tokens_.emplace_back(consumes<std::vector<std::vector<unsigned int>>>(tag));
0142   }
0143 
0144   //make sure that the number of tracksters collections and tracksterlinks collections is the same
0145   assert(egamma_tracksters_tokens_.size() == egamma_tracksterlinks_tokens_.size());
0146 
0147   // Loop over the edm::VInputTag and append the token to general_tracksters_tokens_
0148   // These, instead, are the tracksters already merged by the TrackstersLinksProducer
0149   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("general_tracksters_collections")) {
0150     general_tracksters_tokens_.emplace_back(consumes<std::vector<Trackster>>(tag));
0151   }
0152 
0153   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("general_tracksterlinks_collections")) {
0154     general_tracksterlinks_tokens_.emplace_back(consumes<std::vector<std::vector<unsigned int>>>(tag));
0155   }
0156 
0157   //make sure that the number of tracksters collections and tracksterlinks collections is the same
0158   assert(general_tracksters_tokens_.size() == general_tracksterlinks_tokens_.size());
0159 
0160   //Loop over the edm::VInputTag of masks and append the token to original_masks_tokens_
0161   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("original_masks")) {
0162     original_masks_tokens_.emplace_back(consumes<std::vector<float>>(tag));
0163   }
0164 
0165   std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive";
0166   hdc_token_ =
0167       esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag("", detectorName_));
0168   if (useMTDTiming_) {
0169     inputTimingToken_ = consumes<MtdHostCollection>(ps.getParameter<edm::InputTag>("timingSoA"));
0170   }
0171 
0172   produces<std::vector<TICLCandidate>>();
0173 
0174   // New trackster collection after linking
0175   produces<std::vector<Trackster>>();
0176 
0177   auto interpretationPSet = ps.getParameter<edm::ParameterSet>("interpretationDescPSet");
0178   auto algoType = interpretationPSet.getParameter<std::string>("type");
0179   generalInterpretationAlgo_ =
0180       TICLGeneralInterpretationPluginFactory::get()->create(algoType, interpretationPSet, consumesCollector());
0181 }
0182 
0183 void TICLCandidateProducer::beginRun(edm::Run const &iEvent, edm::EventSetup const &es) {
0184   edm::ESHandle<HGCalDDDConstants> hdc = es.getHandle(hdc_token_);
0185   hgcons_ = hdc.product();
0186 
0187   edm::ESHandle<CaloGeometry> geom = es.getHandle(geometry_token_);
0188   rhtools_.setGeometry(*geom);
0189 
0190   bfield_ = es.getHandle(bfield_token_);
0191   propagator_ = es.getHandle(propagator_token_);
0192   generalInterpretationAlgo_->initialize(hgcons_, rhtools_, bfield_, propagator_);
0193 
0194   trackingGeometry_ = es.getHandle(trackingGeometry_token_);
0195 };
0196 
0197 void filterTracks(edm::Handle<std::vector<reco::Track>> tkH,
0198                   const edm::Handle<std::vector<reco::Muon>> &muons_h,
0199                   const StringCutObjectSelector<reco::Track> cutTk_,
0200                   const float tkEnergyCut_,
0201                   std::vector<bool> &maskTracks) {
0202   auto const &tracks = *tkH;
0203   for (unsigned i = 0; i < tracks.size(); ++i) {
0204     const auto &tk = tracks[i];
0205     reco::TrackRef trackref = reco::TrackRef(tkH, i);
0206 
0207     // veto tracks associated to muons
0208     int muId = PFMuonAlgo::muAssocToTrack(trackref, *muons_h);
0209     const reco::MuonRef muonref = reco::MuonRef(muons_h, muId);
0210 
0211     if (!cutTk_((tk)) or (muId != -1 and PFMuonAlgo::isMuon(muonref) and not(*muons_h)[muId].isTrackerMuon())) {
0212       maskTracks[i] = false;
0213       continue;
0214     }
0215 
0216     // don't consider tracks below 2 GeV for linking
0217     if (std::sqrt(tk.p() * tk.p() + ticl::mpion2) < tkEnergyCut_) {
0218       maskTracks[i] = false;
0219       continue;
0220     }
0221 
0222     // record tracks that can be used to make a ticlcandidate
0223     maskTracks[i] = true;
0224   }
0225 }
0226 
0227 void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) {
0228   auto resultTracksters = std::make_unique<std::vector<Trackster>>();
0229   auto resultTrackstersMerged = std::make_unique<std::vector<Trackster>>();
0230   auto linkedResultTracksters = std::make_unique<std::vector<std::vector<unsigned int>>>();
0231 
0232   const auto &layerClusters = evt.get(clusters_token_);
0233   const auto &layerClustersTimes = evt.get(clustersTime_token_);
0234   edm::Handle<reco::MuonCollection> muons_h;
0235   evt.getByToken(muons_token_, muons_h);
0236 
0237   edm::Handle<std::vector<reco::Track>> tracks_h;
0238   evt.getByToken(tracks_token_, tracks_h);
0239   const auto &tracks = *tracks_h;
0240 
0241   edm::Handle<MtdHostCollection> inputTiming_h;
0242   MtdHostCollection::ConstView inputTimingView;
0243   if (useMTDTiming_) {
0244     evt.getByToken(inputTimingToken_, inputTiming_h);
0245     inputTimingView = (*inputTiming_h).const_view();
0246   }
0247 
0248   auto const bFieldProd = bfield_.product();
0249   const Propagator *propagator = propagator_.product();
0250 
0251   // loop over the original_masks_tokens_ and get the original masks collections and multiply them
0252   // to get the global mask
0253   std::vector<float> original_global_mask(layerClusters.size(), 1.f);
0254   for (unsigned int i = 0; i < original_masks_tokens_.size(); ++i) {
0255     const auto &tmp_mask = evt.get(original_masks_tokens_[i]);
0256     for (unsigned int j = 0; j < tmp_mask.size(); ++j) {
0257       original_global_mask[j] *= tmp_mask[j];
0258     }
0259   }
0260 
0261   auto resultMask = std::make_unique<std::vector<float>>(original_global_mask);
0262 
0263   std::vector<edm::Handle<std::vector<Trackster>>> general_tracksters_h(general_tracksters_tokens_.size());
0264   MultiVectorManager<Trackster> generalTrackstersManager;
0265   for (unsigned int i = 0; i < general_tracksters_tokens_.size(); ++i) {
0266     evt.getByToken(general_tracksters_tokens_[i], general_tracksters_h[i]);
0267     //Fill MultiVectorManager
0268     generalTrackstersManager.addVector(*general_tracksters_h[i]);
0269   }
0270   //now get the general_tracksterlinks_tokens_
0271   std::vector<edm::Handle<std::vector<std::vector<unsigned>>>> general_tracksterlinks_h(
0272       general_tracksterlinks_tokens_.size());
0273   std::vector<std::vector<unsigned>> generalTracksterLinksGlobalId;
0274   for (unsigned int i = 0; i < general_tracksterlinks_tokens_.size(); ++i) {
0275     evt.getByToken(general_tracksterlinks_tokens_[i], general_tracksterlinks_h[i]);
0276     for (unsigned int j = 0; j < general_tracksterlinks_h[i]->size(); ++j) {
0277       generalTracksterLinksGlobalId.emplace_back();
0278       auto &links_vector = generalTracksterLinksGlobalId.back();
0279       links_vector.resize((*general_tracksterlinks_h[i])[j].size());
0280       for (unsigned int k = 0; k < links_vector.size(); ++k) {
0281         links_vector[k] = generalTrackstersManager.getGlobalIndex(i, (*general_tracksterlinks_h[i])[j][k]);
0282       }
0283     }
0284   }
0285 
0286   std::vector<bool> maskTracks;
0287   maskTracks.resize(tracks.size());
0288   filterTracks(tracks_h, muons_h, cutTk_, tkEnergyCut_, maskTracks);
0289 
0290   const typename TICLInterpretationAlgoBase<reco::Track>::Inputs input(evt,
0291                                                                        es,
0292                                                                        layerClusters,
0293                                                                        layerClustersTimes,
0294                                                                        generalTrackstersManager,
0295                                                                        generalTracksterLinksGlobalId,
0296                                                                        tracks_h,
0297                                                                        maskTracks);
0298 
0299   auto resultCandidates = std::make_unique<std::vector<TICLCandidate>>();
0300   std::vector<int> trackstersInTrackIndices(tracks.size(), -1);
0301 
0302   //TODO
0303   //egammaInterpretationAlg_->makecandidates(inputGSF, inputTiming, *resultTrackstersMerged, trackstersInGSFTrackIndices)
0304   // mask generalTracks associated to GSFTrack linked in egammaInterpretationAlgo_
0305 
0306   generalInterpretationAlgo_->makeCandidates(input, inputTiming_h, *resultTracksters, trackstersInTrackIndices);
0307 
0308   assignPCAtoTracksters(*resultTracksters,
0309                         layerClusters,
0310                         layerClustersTimes,
0311                         rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(),
0312                         rhtools_,
0313                         true);
0314 
0315   std::vector<bool> maskTracksters(resultTracksters->size(), true);
0316   edm::OrphanHandle<std::vector<Trackster>> resultTracksters_h = evt.put(std::move(resultTracksters));
0317   //create ChargedCandidates
0318   for (size_t iTrack = 0; iTrack < tracks.size(); iTrack++) {
0319     if (maskTracks[iTrack]) {
0320       auto const tracksterId = trackstersInTrackIndices[iTrack];
0321       auto trackPtr = edm::Ptr<reco::Track>(tracks_h, iTrack);
0322       if (tracksterId != -1 and !maskTracksters.empty()) {
0323         auto tracksterPtr = edm::Ptr<Trackster>(resultTracksters_h, tracksterId);
0324         TICLCandidate chargedCandidate(trackPtr, tracksterPtr);
0325         resultCandidates->push_back(chargedCandidate);
0326         maskTracksters[tracksterId] = false;
0327       } else {
0328         //charged candidates track only
0329         edm::Ptr<Trackster> tracksterPtr;
0330         TICLCandidate chargedCandidate(trackPtr, tracksterPtr);
0331         auto trackRef = edm::Ref<reco::TrackCollection>(tracks_h, iTrack);
0332         const int muId = PFMuonAlgo::muAssocToTrack(trackRef, *muons_h);
0333         const reco::MuonRef muonRef = reco::MuonRef(muons_h, muId);
0334         if (muonRef.isNonnull() and muonRef->isGlobalMuon()) {
0335           // create muon candidate
0336           chargedCandidate.setPdgId(13 * trackPtr.get()->charge());
0337         }
0338         resultCandidates->push_back(chargedCandidate);
0339       }
0340     }
0341   }
0342 
0343   //Neutral Candidate
0344   for (size_t iTrackster = 0; iTrackster < maskTracksters.size(); iTrackster++) {
0345     if (maskTracksters[iTrackster]) {
0346       edm::Ptr<Trackster> tracksterPtr(resultTracksters_h, iTrackster);
0347       edm::Ptr<reco::Track> trackPtr;
0348       TICLCandidate neutralCandidate(trackPtr, tracksterPtr);
0349       resultCandidates->push_back(neutralCandidate);
0350     }
0351   }
0352 
0353   auto getPathLength =
0354       [&](const reco::Track &track, float zVal) {
0355         const auto &fts_inn = trajectoryStateTransform::innerFreeState(track, bFieldProd);
0356         const auto &fts_out = trajectoryStateTransform::outerFreeState(track, bFieldProd);
0357         const auto &surf_inn = trajectoryStateTransform::innerStateOnSurface(track, *trackingGeometry_, bFieldProd);
0358         const auto &surf_out = trajectoryStateTransform::outerStateOnSurface(track, *trackingGeometry_, bFieldProd);
0359 
0360         Basic3DVector<float> pos(track.referencePoint());
0361         Basic3DVector<float> mom(track.momentum());
0362         FreeTrajectoryState stateAtBeamspot{GlobalPoint(pos), GlobalVector(mom), track.charge(), bFieldProd};
0363 
0364         float pathlength = propagator->propagateWithPath(stateAtBeamspot, surf_inn.surface()).second;
0365 
0366         if (pathlength) {
0367           const auto &t_inn_out = propagator->propagateWithPath(fts_inn, surf_out.surface());
0368 
0369           if (t_inn_out.first.isValid()) {
0370             pathlength += t_inn_out.second;
0371 
0372             std::pair<float, float> rMinMax = hgcons_->rangeR(zVal, true);
0373 
0374             int iSide = int(track.eta() > 0);
0375             float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
0376             const auto &disk = std::make_unique<GeomDet>(
0377                 Disk::build(Disk::PositionType(0, 0, zSide),
0378                             Disk::RotationType(),
0379                             SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5))
0380                     .get());
0381             const auto &tsos = propagator->propagateWithPath(fts_out, disk->surface());
0382 
0383             if (tsos.first.isValid()) {
0384               pathlength += tsos.second;
0385               return pathlength;
0386             }
0387           }
0388         }
0389 #ifdef EDM_ML_DEBUG
0390         LogDebug("TICLCandidateProducer")
0391             << "Not able to use the track to compute the path length. A straight line will be used instead.";
0392 #endif
0393         return 0.f;
0394       };
0395 
0396   assignTimeToCandidates(*resultCandidates, tracks_h, inputTimingView, getPathLength);
0397 
0398   evt.put(std::move(resultCandidates));
0399 }
0400 
0401 template <typename F>
0402 void TICLCandidateProducer::assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates,
0403                                                    edm::Handle<std::vector<reco::Track>> track_h,
0404                                                    MtdHostCollection::ConstView &inputTimingView,
0405                                                    F func) const {
0406   for (auto &cand : resultCandidates) {
0407     float beta = 1;
0408     float time = 0.f;
0409     float invTimeErr = 0.f;
0410     float timeErr = -1.f;
0411 
0412     for (const auto &tr : cand.tracksters()) {
0413       if (tr->timeError() > 0) {
0414         const auto invTimeESq = pow(tr->timeError(), -2);
0415         const auto x = tr->barycenter().X();
0416         const auto y = tr->barycenter().Y();
0417         const auto z = tr->barycenter().Z();
0418         auto path = std::sqrt(x * x + y * y + z * z);
0419         if (cand.trackPtr().get() != nullptr) {
0420           const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
0421           if (useMTDTiming_ and inputTimingView.timeErr()[trackIndex] > 0) {
0422             const auto xMtd = inputTimingView.posInMTD_x()[trackIndex];
0423             const auto yMtd = inputTimingView.posInMTD_y()[trackIndex];
0424             const auto zMtd = inputTimingView.posInMTD_z()[trackIndex];
0425 
0426             beta = inputTimingView.beta()[trackIndex];
0427             path = std::sqrt((x - xMtd) * (x - xMtd) + (y - yMtd) * (y - yMtd) + (z - zMtd) * (z - zMtd)) +
0428                    inputTimingView.pathLength()[trackIndex];
0429           } else {
0430             float pathLength = func(*(cand.trackPtr().get()), z);
0431             if (pathLength) {
0432               path = pathLength;
0433             }
0434           }
0435         }
0436         time += (tr->time() - path / (beta * c_light_)) * invTimeESq;
0437         invTimeErr += invTimeESq;
0438       }
0439     }
0440     if (invTimeErr > 0) {
0441       time = time / invTimeErr;
0442       // FIXME_ set a liminf of 0.02 ns on the ts error (based on residuals)
0443       timeErr = sqrt(1.f / invTimeErr);
0444       if (timeErr < timeRes)
0445         timeErr = timeRes;
0446       cand.setTime(time, timeErr);
0447     }
0448 
0449     if (useMTDTiming_ and cand.charge()) {
0450       // Check MTD timing availability
0451       const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
0452       const bool assocQuality = inputTimingView.MVAquality()[trackIndex] > timingQualityThreshold_;
0453       if (assocQuality) {
0454         const auto timeHGC = cand.time();
0455         const auto timeEHGC = cand.timeError();
0456         const auto timeMTD = inputTimingView.time0()[trackIndex];
0457         const auto timeEMTD = inputTimingView.time0Err()[trackIndex];
0458 
0459         if (useTimingAverage_ && (timeEMTD > 0 && timeEHGC > 0)) {
0460           // Compute weighted average between HGCAL and MTD timing
0461           const auto invTimeESqHGC = pow(timeEHGC, -2);
0462           const auto invTimeESqMTD = pow(timeEMTD, -2);
0463           timeErr = 1.f / (invTimeESqHGC + invTimeESqMTD);
0464           time = (timeHGC * invTimeESqHGC + timeMTD * invTimeESqMTD) * timeErr;
0465           timeErr = sqrt(timeErr);
0466         } else if (timeEMTD > 0) {
0467           time = timeMTD;
0468           timeErr = timeEMTD;
0469         }
0470       }
0471       cand.setTime(time, timeErr);
0472       cand.setMTDTime(inputTimingView.time()[trackIndex], inputTimingView.timeErr()[trackIndex]);
0473     }
0474   }
0475 }
0476 
0477 void TICLCandidateProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0478   edm::ParameterSetDescription desc;
0479   edm::ParameterSetDescription interpretationDesc;
0480   interpretationDesc.addNode(edm::PluginDescription<TICLGeneralInterpretationPluginFactory>("type", "General", true));
0481   desc.add<edm::ParameterSetDescription>("interpretationDescPSet", interpretationDesc);
0482   desc.add<std::vector<edm::InputTag>>("egamma_tracksters_collections", {edm::InputTag("ticlTracksterLinks")});
0483   desc.add<std::vector<edm::InputTag>>("egamma_tracksterlinks_collections", {edm::InputTag("ticlTracksterLinks")});
0484   desc.add<std::vector<edm::InputTag>>("general_tracksters_collections", {edm::InputTag("ticlTracksterLinks")});
0485   desc.add<std::vector<edm::InputTag>>("general_tracksterlinks_collections", {edm::InputTag("ticlTracksterLinks")});
0486   desc.add<std::vector<edm::InputTag>>("original_masks",
0487                                        {edm::InputTag("hgcalMergeLayerClusters", "InitialLayerClustersMask")});
0488   desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalMergeLayerClusters"));
0489   desc.add<edm::InputTag>("layer_clustersTime", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster"));
0490   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0491   desc.add<edm::InputTag>("timingSoA", edm::InputTag("mtdSoA"));
0492   desc.add<edm::InputTag>("muons", edm::InputTag("muons1stStep"));
0493   desc.add<std::string>("detector", "HGCAL");
0494   desc.add<std::string>("propagator", "PropagatorWithMaterial");
0495   desc.add<bool>("useMTDTiming", true);
0496   desc.add<bool>("useTimingAverage", true);
0497   desc.add<double>("timingQualityThreshold", 0.5);
0498   desc.add<std::string>("cutTk",
0499                         "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
0500                         "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
0501   descriptions.add("ticlCandidateProducer", desc);
0502 }
0503 
0504 DEFINE_FWK_MODULE(TICLCandidateProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!