Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:59:08 

0001 // -*- C++ -*-
0002 //
0003 // Package:    ElectronProducers
0004 // Class:      LowPtGsfElectronSeedProducer
0005 //
0006 /**\class LowPtGsfElectronSeedProducer RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedProducer.cc
0007  Description: EDProducer of ElectronSeed objects
0008  Implementation:
0009      <Notes on implementation>
0010 */
0011 // Original Author:  Robert Bainbridge
0012 
0013 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0014 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0015 #include "DataFormats/Common/interface/ValueMap.h"
0016 #include "DataFormats/Math/interface/deltaPhi.h"
0017 #include "DataFormats/Math/interface/deltaR.h"
0018 #include "DataFormats/Math/interface/LorentzVector.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFTrajectoryPoint.h"
0020 #include "DataFormats/TrackReco/interface/TrackBase.h"
0021 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0022 #include "CommonTools/BaseParticlePropagator/interface/RawParticle.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0025 #include "FWCore/Utilities/interface/InputTag.h"
0026 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0027 #include "RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h"
0028 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0029 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0030 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0031 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0032 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0033 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0034 #include "DataFormats/Common/interface/Handle.h"
0035 #include "DataFormats/Common/interface/RefProd.h"
0036 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0037 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0038 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
0039 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0040 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0041 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
0042 #include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
0043 #include "DataFormats/ParticleFlowReco/interface/PreId.h"
0044 #include "DataFormats/ParticleFlowReco/interface/PreIdFwd.h"
0045 #include "DataFormats/TrackReco/interface/Track.h"
0046 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0047 #include "FWCore/Framework/interface/ESHandle.h"
0048 #include "FWCore/Framework/interface/stream/EDProducer.h"
0049 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0050 #include "MagneticField/Engine/interface/MagneticField.h"
0051 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
0052 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
0053 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0054 
0055 #include "LowPtGsfElectronSeedHeavyObjectCache.h"
0056 
0057 class LowPtGsfElectronSeedProducer final
0058     : public edm::stream::EDProducer<edm::GlobalCache<lowptgsfeleseed::HeavyObjectCache> > {
0059 public:
0060   using TrackIndxMap = std::unordered_map<reco::TrackRef::key_type, size_t>;
0061   explicit LowPtGsfElectronSeedProducer(const edm::ParameterSet&, const lowptgsfeleseed::HeavyObjectCache*);
0062 
0063   static std::unique_ptr<lowptgsfeleseed::HeavyObjectCache> initializeGlobalCache(const edm::ParameterSet& conf) {
0064     return std::make_unique<lowptgsfeleseed::HeavyObjectCache>(lowptgsfeleseed::HeavyObjectCache(conf));
0065   }
0066 
0067   static void globalEndJob(lowptgsfeleseed::HeavyObjectCache const*) {}
0068 
0069   void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
0070 
0071   void produce(edm::Event&, const edm::EventSetup&) override;
0072 
0073   static void fillDescriptions(edm::ConfigurationDescriptions&);
0074 
0075 private:  // member functions
0076   template <typename T>
0077   void loop(const edm::Handle<std::vector<T> >& handle,
0078             edm::Handle<reco::PFClusterCollection>& hcalClusters,
0079             reco::ElectronSeedCollection& seeds,
0080             reco::PreIdCollection& ecalPreIds,
0081             reco::PreIdCollection& hcalPreIds,
0082             TrackIndxMap& trksToPreIdIndx,
0083             edm::Event&,
0084             const edm::EventSetup&);
0085 
0086   // Overloaded methods to retrieve reco::TrackRef
0087 
0088   reco::TrackRef getBaseRef(edm::Handle<std::vector<reco::Track> > handle, int idx) const;
0089   reco::TrackRef getBaseRef(edm::Handle<std::vector<reco::PFRecTrack> > handle, int idx) const;
0090 
0091   // Overloaded methods to populate PreIds (using PF or KF tracks)
0092 
0093   void propagateTrackToCalo(const reco::PFRecTrackRef& pfTrackRef,
0094                             const edm::Handle<reco::PFClusterCollection>& ecalClusters,
0095                             const edm::Handle<reco::PFClusterCollection>& hcalClusters,
0096                             std::vector<int>& matchedEcalClusters,
0097                             std::vector<int>& matchedHcalClusters,
0098                             reco::PreId& ecalPreId,
0099                             reco::PreId& hcalPreId);
0100 
0101   void propagateTrackToCalo(const reco::PFRecTrackRef& pfTrackRef,
0102                             const edm::Handle<reco::PFClusterCollection>& clusters,
0103                             std::vector<int>& matchedClusters,
0104                             reco::PreId& preId,
0105                             bool ecal);
0106 
0107   void propagateTrackToCalo(const reco::TrackRef& pfTrack,
0108                             const edm::Handle<reco::PFClusterCollection>& ecalClusters,
0109                             const edm::Handle<reco::PFClusterCollection>& hcalClusters,
0110                             std::vector<int>& matchedEcalClusters,
0111                             std::vector<int>& matchedHcalClusters,
0112                             reco::PreId& ecalPreId,
0113                             reco::PreId& hcalPreId);
0114   template <typename CollType>
0115   void fillPreIdRefValueMap(edm::Handle<CollType> tracksHandle,
0116                             const TrackIndxMap& trksToPreIdIndx,
0117                             const edm::OrphanHandle<reco::PreIdCollection>& preIdHandle,
0118                             edm::ValueMap<reco::PreIdRef>::Filler& filler);
0119 
0120   // Overloaded methods to evaluate BDTs (using PF or KF tracks)
0121 
0122   bool decision(const reco::PFRecTrackRef& pfTrackRef,
0123                 reco::PreId& ecal,
0124                 reco::PreId& hcal,
0125                 double rho,
0126                 const reco::BeamSpot& spot,
0127                 noZS::EcalClusterLazyTools& ecalTools);
0128 
0129   bool decision(const reco::TrackRef& kfTrackRef,
0130                 reco::PreId& ecal,
0131                 reco::PreId& hcal,
0132                 double rho,
0133                 const reco::BeamSpot& spot,
0134                 noZS::EcalClusterLazyTools& ecalTools);
0135 
0136   // Perform lightweight GSF tracking
0137   bool lightGsfTracking(reco::PreId&, const reco::TrackRef&, const reco::ElectronSeed&);
0138 
0139 private:  // member data
0140   edm::ESHandle<MagneticField> field_;
0141   std::unique_ptr<TrajectoryFitter> fitterPtr_;
0142   std::unique_ptr<TrajectorySmoother> smootherPtr_;
0143   edm::EDGetTokenT<reco::TrackCollection> kfTracks_;
0144   edm::EDGetTokenT<reco::PFRecTrackCollection> pfTracks_;
0145   const edm::EDGetTokenT<reco::PFClusterCollection> ecalClusters_;
0146   edm::EDGetTokenT<reco::PFClusterCollection> hcalClusters_;
0147   const edm::EDGetTokenT<EcalRecHitCollection> ebRecHits_;
0148   const edm::EDGetTokenT<EcalRecHitCollection> eeRecHits_;
0149   const edm::EDGetTokenT<double> rho_;
0150   const edm::EDGetTokenT<reco::BeamSpot> beamSpot_;
0151 
0152   const edm::ESGetToken<TrajectoryFitter, TrajectoryFitter::Record> trajectoryFitterToken_;
0153   const edm::ESGetToken<TrajectorySmoother, TrajectoryFitter::Record> trajectorySmootherToken_;
0154   const edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> builderToken_;
0155   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magToken_;
0156   const noZS::EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_;
0157 
0158   const bool passThrough_;
0159   const bool usePfTracks_;
0160   const double minPtThreshold_;
0161   const double maxPtThreshold_;
0162 
0163   // pow( sinh(1.65), 2. )
0164   static constexpr double boundary_ = 2.50746495928 * 2.50746495928;
0165   // pow( ele_mass, 2. )
0166   static constexpr double mass_ = 0.000511 * 0.000511;
0167 };
0168 
0169 //////////////////////////////////////////////////////////////////////////////////////////
0170 //
0171 LowPtGsfElectronSeedProducer::LowPtGsfElectronSeedProducer(const edm::ParameterSet& conf,
0172                                                            const lowptgsfeleseed::HeavyObjectCache*)
0173     : field_(),
0174       fitterPtr_(),
0175       smootherPtr_(),
0176       kfTracks_(),
0177       pfTracks_(),
0178       ecalClusters_{consumes(conf.getParameter<edm::InputTag>("ecalClusters"))},
0179       hcalClusters_(),
0180       ebRecHits_{consumes(conf.getParameter<edm::InputTag>("EBRecHits"))},
0181       eeRecHits_{consumes(conf.getParameter<edm::InputTag>("EERecHits"))},
0182       rho_(consumes(conf.getParameter<edm::InputTag>("rho"))),
0183       beamSpot_(consumes(conf.getParameter<edm::InputTag>("BeamSpot"))),
0184       trajectoryFitterToken_{esConsumes(conf.getParameter<edm::ESInputTag>("Fitter"))},
0185       trajectorySmootherToken_{esConsumes(conf.getParameter<edm::ESInputTag>("Smoother"))},
0186       builderToken_{esConsumes(conf.getParameter<edm::ESInputTag>("TTRHBuilder"))},
0187       magToken_{esConsumes<edm::Transition::BeginLuminosityBlock>()},
0188       ecalClusterToolsESGetTokens_{consumesCollector()},
0189       passThrough_(conf.getParameter<bool>("PassThrough")),
0190       usePfTracks_(conf.getParameter<bool>("UsePfTracks")),
0191       minPtThreshold_(conf.getParameter<double>("MinPtThreshold")),
0192       maxPtThreshold_(conf.getParameter<double>("MaxPtThreshold")) {
0193   if (usePfTracks_) {
0194     pfTracks_ = consumes(conf.getParameter<edm::InputTag>("pfTracks"));
0195     hcalClusters_ = consumes(conf.getParameter<edm::InputTag>("hcalClusters"));
0196   }
0197   kfTracks_ = consumes(conf.getParameter<edm::InputTag>("tracks"));
0198 
0199   produces<reco::ElectronSeedCollection>();
0200   produces<reco::PreIdCollection>();
0201   produces<reco::PreIdCollection>("HCAL");
0202   produces<edm::ValueMap<reco::PreIdRef> >();  // indexed by edm::Ref<ElectronSeed>.index()
0203 }
0204 
0205 //////////////////////////////////////////////////////////////////////////////////////////
0206 //
0207 void LowPtGsfElectronSeedProducer::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const& setup) {
0208   field_ = setup.getHandle(magToken_);
0209 }
0210 
0211 //////////////////////////////////////////////////////////////////////////////////////////
0212 //
0213 void LowPtGsfElectronSeedProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
0214   // Products
0215   auto seeds = std::make_unique<reco::ElectronSeedCollection>();
0216   auto ecalPreIds = std::make_unique<reco::PreIdCollection>();
0217   auto hcalPreIds = std::make_unique<reco::PreIdCollection>();
0218 
0219   const edm::RefProd<reco::PreIdCollection> preIdsRefProd = event.getRefBeforePut<reco::PreIdCollection>();
0220 
0221   // HCAL clusters (only used with PF tracks)
0222   edm::Handle<reco::PFClusterCollection> hcalClusters;
0223 
0224   //we always need kftracks as we link the preids to them
0225   edm::Handle<reco::TrackCollection> kfTracks;
0226   event.getByToken(kfTracks_, kfTracks);
0227 
0228   TrackIndxMap trksToPreIdIndx;
0229   if (usePfTracks_) {
0230     edm::Handle<reco::PFRecTrackCollection> pfTracks;
0231     event.getByToken(pfTracks_, pfTracks);
0232     event.getByToken(hcalClusters_, hcalClusters);
0233 
0234     //check consistency between kfTracks and pfTracks collection
0235     for (auto& trk : *pfTracks) {
0236       if (trk.trackRef().isNonnull()) {
0237         if (trk.trackRef().id() != kfTracks.id()) {
0238           throw cms::Exception("ConfigError")
0239               << "kfTracks is not the collection that pfTracks was built from, please fix this";
0240         }
0241         break;  //we only need one valid track ref to check this so end the loop
0242       }
0243     }
0244 
0245     loop(pfTracks,  // PF tracks
0246          hcalClusters,
0247          *seeds,
0248          *ecalPreIds,
0249          *hcalPreIds,
0250          trksToPreIdIndx,
0251          event,
0252          setup);
0253 
0254   } else {
0255     loop(kfTracks,  // KF tracks
0256          hcalClusters,
0257          *seeds,
0258          *ecalPreIds,
0259          *hcalPreIds,
0260          trksToPreIdIndx,
0261          event,
0262          setup);
0263   }
0264 
0265   auto ecalPreIdsHandle = event.put(std::move(ecalPreIds));
0266   event.put(std::move(hcalPreIds), "HCAL");
0267   event.put(std::move(seeds));
0268 
0269   auto preIdVMOut = std::make_unique<edm::ValueMap<reco::PreIdRef> >();
0270   edm::ValueMap<reco::PreIdRef>::Filler mapFiller(*preIdVMOut);
0271   fillPreIdRefValueMap(kfTracks, trksToPreIdIndx, ecalPreIdsHandle, mapFiller);
0272   mapFiller.fill();
0273   event.put(std::move(preIdVMOut));
0274 }
0275 
0276 //////////////////////////////////////////////////////////////////////////////////////////
0277 // Return reco::Track from edm::Ref<T>
0278 
0279 reco::TrackRef LowPtGsfElectronSeedProducer::getBaseRef(edm::Handle<std::vector<reco::Track> > handle, int idx) const {
0280   return reco::TrackRef(handle, idx);
0281 }
0282 
0283 reco::TrackRef LowPtGsfElectronSeedProducer::getBaseRef(edm::Handle<std::vector<reco::PFRecTrack> > handle,
0284                                                         int idx) const {
0285   return reco::PFRecTrackRef(handle, idx)->trackRef();
0286 }
0287 
0288 //////////////////////////////////////////////////////////////////////////////////////////
0289 // Template function, instantiated for both reco::Tracks and reco::PFRecTracks
0290 template <typename T>
0291 void LowPtGsfElectronSeedProducer::loop(const edm::Handle<std::vector<T> >& handle,  // PF or KF tracks
0292                                         edm::Handle<reco::PFClusterCollection>& hcalClusters,
0293                                         reco::ElectronSeedCollection& seeds,
0294                                         reco::PreIdCollection& ecalPreIds,
0295                                         reco::PreIdCollection& hcalPreIds,
0296                                         TrackIndxMap& trksToPreIdIndx,
0297                                         edm::Event& event,
0298                                         const edm::EventSetup& setup) {
0299   // Pileup
0300   auto const& rho = event.get(rho_);
0301 
0302   // Beam spot
0303   auto const& spot = event.get(beamSpot_);
0304 
0305   // Track fitter
0306   fitterPtr_ = setup.getData(trajectoryFitterToken_).clone();
0307 
0308   // Track smoother
0309   smootherPtr_.reset(setup.getData(trajectorySmootherToken_).clone());
0310 
0311   // RecHit cloner
0312   TkClonerImpl hitCloner = static_cast<TkTransientTrackingRecHitBuilder const&>(setup.getData(builderToken_)).cloner();
0313   fitterPtr_->setHitCloner(&hitCloner);
0314   smootherPtr_->setHitCloner(&hitCloner);
0315 
0316   // ECAL clusters
0317   auto ecalClusters = event.getHandle(ecalClusters_);
0318 
0319   // Utility to access to shower shape vars
0320   noZS::EcalClusterLazyTools ecalTools(event, ecalClusterToolsESGetTokens_.get(setup), ebRecHits_, eeRecHits_);
0321 
0322   // Ensure each cluster is only matched once to a track
0323   std::vector<int> matchedEcalClusters;
0324   std::vector<int> matchedHcalClusters;
0325 
0326   // Reserve
0327   seeds.reserve(handle->size());
0328   ecalPreIds.reserve(handle->size());
0329   hcalPreIds.reserve(handle->size());
0330 
0331   // Iterate through (PF or KF) tracks
0332   for (unsigned int itrk = 0; itrk < handle.product()->size(); itrk++) {
0333     edm::Ref<std::vector<T> > templatedRef(handle, itrk);  // TrackRef or PFRecTrackRef
0334     reco::TrackRef trackRef = getBaseRef(handle, itrk);
0335 
0336     if (!(trackRef->quality(reco::TrackBase::qualityByName("highPurity")))) {
0337       continue;
0338     }
0339     if (!passThrough_ && (trackRef->pt() < minPtThreshold_)) {
0340       continue;
0341     }
0342 
0343     // Create ElectronSeed
0344     reco::ElectronSeed seed(*(trackRef->seedRef()));
0345     if (seed.nHits() == 0) {  //if DeepCore is used in jetCore iteration the seed are hitless, in case skip
0346       continue;
0347     }
0348     seed.setCtfTrack(trackRef);
0349 
0350     // Create PreIds
0351     unsigned int nModels = globalCache()->modelNames().size();
0352     reco::PreId ecalPreId(nModels);
0353     reco::PreId hcalPreId(nModels);
0354 
0355     // Add track ref to PreId
0356     ecalPreId.setTrack(trackRef);
0357     hcalPreId.setTrack(trackRef);
0358 
0359     // Add Track-Calo matching variables to PreIds
0360     propagateTrackToCalo(
0361         templatedRef, ecalClusters, hcalClusters, matchedEcalClusters, matchedHcalClusters, ecalPreId, hcalPreId);
0362 
0363     // Add variables related to GSF tracks to PreId
0364     lightGsfTracking(ecalPreId, trackRef, seed);
0365 
0366     // Decision based on BDT
0367     bool result = decision(templatedRef, ecalPreId, hcalPreId, rho, spot, ecalTools);
0368 
0369     // If fails BDT, do not store seed
0370     if (!result) {
0371       continue;
0372     }
0373 
0374     // Store PreId
0375     ecalPreIds.push_back(ecalPreId);
0376     hcalPreIds.push_back(hcalPreId);
0377     trksToPreIdIndx[trackRef.index()] = ecalPreIds.size() - 1;
0378 
0379     // Store ElectronSeed and corresponding edm::Ref<PreId>.index()
0380     seeds.push_back(seed);
0381   }
0382 }
0383 
0384 //////////////////////////////////////////////////////////////////////////////////////////
0385 // Template instantiation for reco::Tracks
0386 template void LowPtGsfElectronSeedProducer::loop<reco::Track>(const edm::Handle<std::vector<reco::Track> >&,
0387                                                               edm::Handle<reco::PFClusterCollection>& hcalClusters,
0388                                                               reco::ElectronSeedCollection& seeds,
0389                                                               reco::PreIdCollection& ecalPreIds,
0390                                                               reco::PreIdCollection& hcalPreIds,
0391                                                               TrackIndxMap& trksToPreIdIndx,
0392                                                               edm::Event&,
0393                                                               const edm::EventSetup&);
0394 
0395 //////////////////////////////////////////////////////////////////////////////////////////
0396 // Template instantiation for reco::PFRecTracks
0397 template void LowPtGsfElectronSeedProducer::loop<reco::PFRecTrack>(const edm::Handle<std::vector<reco::PFRecTrack> >&,
0398                                                                    edm::Handle<reco::PFClusterCollection>& hcalClusters,
0399                                                                    reco::ElectronSeedCollection& seeds,
0400                                                                    reco::PreIdCollection& ecalPreIds,
0401                                                                    reco::PreIdCollection& hcalPreIds,
0402                                                                    TrackIndxMap& trksToPreIdIndx,
0403                                                                    edm::Event&,
0404                                                                    const edm::EventSetup&);
0405 
0406 //////////////////////////////////////////////////////////////////////////////////////////
0407 // Loops through both ECAL and HCAL clusters
0408 void LowPtGsfElectronSeedProducer::propagateTrackToCalo(const reco::PFRecTrackRef& pfTrackRef,
0409                                                         const edm::Handle<reco::PFClusterCollection>& ecalClusters,
0410                                                         const edm::Handle<reco::PFClusterCollection>& hcalClusters,
0411                                                         std::vector<int>& matchedEcalClusters,
0412                                                         std::vector<int>& matchedHcalClusters,
0413                                                         reco::PreId& ecalPreId,
0414                                                         reco::PreId& hcalPreId) {
0415   propagateTrackToCalo(pfTrackRef, ecalClusters, matchedEcalClusters, ecalPreId, true);
0416   propagateTrackToCalo(pfTrackRef, hcalClusters, matchedHcalClusters, hcalPreId, false);
0417 }
0418 
0419 //////////////////////////////////////////////////////////////////////////////////////////
0420 // Loops through ECAL or HCAL clusters (called twice)
0421 void LowPtGsfElectronSeedProducer::propagateTrackToCalo(const reco::PFRecTrackRef& pfTrackRef,
0422                                                         const edm::Handle<reco::PFClusterCollection>& clusters,
0423                                                         std::vector<int>& matched,
0424                                                         reco::PreId& preId,
0425                                                         bool ecal) {
0426   // Store info for PreId
0427   struct Info {
0428     reco::PFClusterRef cluRef = reco::PFClusterRef();
0429     float dr2min = 1.e6;
0430     float deta = 1.e6;
0431     float dphi = 1.e6;
0432     math::XYZPoint showerPos = math::XYZPoint(0., 0., 0.);
0433   } info;
0434 
0435   // Find closest "seed cluster" to KF track extrapolated to ECAL (or HCAL)
0436   reco::PFTrajectoryPoint point;
0437   if (ecal) {
0438     point = pfTrackRef->extrapolatedPoint(reco::PFTrajectoryPoint::LayerType::ECALShowerMax);
0439   } else {
0440     point = pfTrackRef->extrapolatedPoint(reco::PFTrajectoryPoint::LayerType::HCALEntrance);
0441   }
0442 
0443   if (point.isValid()) {
0444     Info info;
0445     for (unsigned int iclu = 0; iclu < clusters.product()->size(); iclu++) {
0446       if (std::find(matched.begin(), matched.end(), iclu) == matched.end()) {
0447         reco::PFClusterRef cluRef(clusters, iclu);
0448 
0449         // Determine dR squared
0450         float dr2 = reco::deltaR2(cluRef->positionREP(), point.positionREP());
0451 
0452         if (dr2 < info.dr2min) {
0453           info.dr2min = dr2;
0454           info.cluRef = cluRef;
0455           info.deta = cluRef->positionREP().eta() - point.positionREP().eta();
0456           info.dphi =
0457               reco::deltaPhi(cluRef->positionREP().phi(), point.positionREP().phi()) * pfTrackRef->trackRef()->charge();
0458           info.showerPos = point.position();
0459         }
0460       }
0461     }
0462 
0463     // Set PreId content if match found
0464     if (info.dr2min < 1.e5) {
0465       float ep = info.cluRef->correctedEnergy() / std::sqrt(pfTrackRef->trackRef()->innerMomentum().mag2());
0466       preId.setECALMatchingProperties(info.cluRef,
0467                                       point.position(),  // ECAL or HCAL surface
0468                                       info.showerPos,    //
0469                                       info.deta,
0470                                       info.dphi,
0471                                       0.f,                                       // chieta
0472                                       0.f,                                       // chiphi
0473                                       pfTrackRef->trackRef()->normalizedChi2(),  // chi2
0474                                       ep);
0475     }
0476 
0477   }  // clusters
0478 }
0479 
0480 //////////////////////////////////////////////////////////////////////////////////////////
0481 // Original implementation in GoodSeedProducer, loops over ECAL clusters only
0482 void LowPtGsfElectronSeedProducer::propagateTrackToCalo(
0483     const reco::TrackRef& kfTrackRef,
0484     const edm::Handle<reco::PFClusterCollection>& ecalClusters,
0485     const edm::Handle<reco::PFClusterCollection>& hcalClusters,  // not used
0486     std::vector<int>& matchedEcalClusters,
0487     std::vector<int>& matchedHcalClusters,  // not used
0488     reco::PreId& ecalPreId,
0489     reco::PreId& hcalPreId /* not used */) {
0490   // Store info for PreId
0491   struct Info {
0492     reco::PFClusterRef cluRef = reco::PFClusterRef();
0493     float dr2min = 1.e6;
0494     float deta = 1.e6;
0495     float dphi = 1.e6;
0496     math::XYZPoint showerPos = math::XYZPoint(0., 0., 0.);
0497   } info;
0498 
0499   // Propagate 'electron' to ECAL surface
0500   float energy = sqrt(mass_ + kfTrackRef->outerMomentum().Mag2());
0501   XYZTLorentzVector mom = XYZTLorentzVector(
0502       kfTrackRef->outerMomentum().x(), kfTrackRef->outerMomentum().y(), kfTrackRef->outerMomentum().z(), energy);
0503   XYZTLorentzVector pos = XYZTLorentzVector(
0504       kfTrackRef->outerPosition().x(), kfTrackRef->outerPosition().y(), kfTrackRef->outerPosition().z(), 0.);
0505   math::XYZVector field(field_->inTesla(GlobalPoint(0, 0, 0)));
0506   BaseParticlePropagator particle(RawParticle(mom, pos, kfTrackRef->charge()), 0, 0, field.z());
0507   particle.propagateToEcalEntrance(false);
0508   if (particle.getSuccess() == 0) {
0509     return;
0510   }
0511 
0512   // ECAL entry point for track
0513   GlobalPoint ecal_pos(
0514       particle.particle().vertex().x(), particle.particle().vertex().y(), particle.particle().vertex().z());
0515 
0516   // Preshower limit
0517   bool below_ps = pow(ecal_pos.z(), 2.) > boundary_ * ecal_pos.perp2();
0518 
0519   // Iterate through ECAL clusters
0520   for (unsigned int iclu = 0; iclu < ecalClusters.product()->size(); iclu++) {
0521     reco::PFClusterRef cluRef(ecalClusters, iclu);
0522 
0523     // Correct ecal_pos for shower depth
0524     double shower_depth = reco::PFCluster::getDepthCorrection(cluRef->correctedEnergy(), below_ps, false);
0525     GlobalPoint showerPos = ecal_pos + GlobalVector(particle.particle().momentum().x(),
0526                                                     particle.particle().momentum().y(),
0527                                                     particle.particle().momentum().z())
0528                                                .unit() *
0529                                            shower_depth;
0530 
0531     // Determine dR squared
0532     float dr2 = reco::deltaR2(cluRef->positionREP(), showerPos);
0533 
0534     // Find nearest ECAL cluster
0535     if (dr2 < info.dr2min) {
0536       info.dr2min = dr2;
0537       info.cluRef = cluRef;
0538       info.deta = std::abs(cluRef->positionREP().eta() - showerPos.eta());
0539       info.dphi = std::abs(reco::deltaPhi(cluRef->positionREP().phi(), showerPos.phi())) * kfTrackRef->charge();
0540       info.showerPos = showerPos;
0541     }
0542   }
0543 
0544   // Populate PreId object
0545   math::XYZPoint point(ecal_pos.x(), ecal_pos.y(), ecal_pos.z());
0546 
0547   // Set PreId content
0548   ecalPreId.setECALMatchingProperties(
0549       info.cluRef,
0550       point,
0551       info.showerPos,
0552       info.deta,
0553       info.dphi,
0554       0.f,                                                                              // chieta
0555       0.f,                                                                              // chiphi
0556       kfTrackRef->normalizedChi2(),                                                     // chi2
0557       info.cluRef->correctedEnergy() / std::sqrt(kfTrackRef->innerMomentum().mag2()));  // E/p
0558 }
0559 
0560 //////////////////////////////////////////////////////////////////////////////////////////
0561 // Original implementation for "lightweight" GSF tracking
0562 bool LowPtGsfElectronSeedProducer::lightGsfTracking(reco::PreId& preId,
0563                                                     const reco::TrackRef& trackRef,
0564                                                     const reco::ElectronSeed& seed) {
0565   Trajectory::ConstRecHitContainer hits;
0566   for (unsigned int ihit = 0; ihit < trackRef->recHitsSize(); ++ihit) {
0567     hits.push_back(trackRef->recHit(ihit)->cloneSH());
0568   }
0569 
0570   GlobalVector gv(trackRef->innerMomentum().x(), trackRef->innerMomentum().y(), trackRef->innerMomentum().z());
0571   GlobalPoint gp(trackRef->innerPosition().x(), trackRef->innerPosition().y(), trackRef->innerPosition().z());
0572 
0573   GlobalTrajectoryParameters gtps(gp, gv, trackRef->charge(), &*field_);
0574 
0575   TrajectoryStateOnSurface tsos(gtps, trackRef->innerStateCovariance(), *hits[0]->surface());
0576 
0577   // Track fitted and smoothed under electron hypothesis
0578   Trajectory traj1 = fitterPtr_->fitOne(seed, hits, tsos);
0579   if (!traj1.isValid()) {
0580     return false;
0581   }
0582   Trajectory traj2 = smootherPtr_->trajectory(traj1);
0583   if (!traj2.isValid()) {
0584     return false;
0585   }
0586 
0587   // Set PreId content
0588   float chi2Ratio = trackRef->chi2() > 0. ? traj2.chiSquared() / trackRef->chi2() : -1.;
0589   float gsfReducedChi2 = chi2Ratio > -1. ? chi2Ratio * trackRef->normalizedChi2() : -1.;
0590   float ptOut = traj2.firstMeasurement().updatedState().globalMomentum().perp();
0591   float ptIn = traj2.lastMeasurement().updatedState().globalMomentum().perp();
0592   float gsfDpt = (ptIn > 0) ? fabs(ptOut - ptIn) / ptIn : 0.;
0593   preId.setTrackProperties(gsfReducedChi2, chi2Ratio, gsfDpt);
0594 
0595   return true;
0596 }
0597 
0598 //////////////////////////////////////////////////////////////////////////////////////////
0599 // Decision based on OR of outputs from list of models
0600 bool LowPtGsfElectronSeedProducer::decision(const reco::PFRecTrackRef& pfTrackRef,
0601                                             reco::PreId& ecalPreId,
0602                                             reco::PreId& hcalPreId,
0603                                             double rho,
0604                                             const reco::BeamSpot& spot,
0605                                             noZS::EcalClusterLazyTools& ecalTools) {
0606   bool result = false;
0607   for (auto& name : globalCache()->modelNames()) {
0608     result |= globalCache()->eval(name, ecalPreId, hcalPreId, rho, spot, ecalTools);
0609   }
0610   return passThrough_ || (pfTrackRef->trackRef()->pt() > maxPtThreshold_) || result;
0611 }
0612 
0613 //////////////////////////////////////////////////////////////////////////////////////////
0614 //
0615 bool LowPtGsfElectronSeedProducer::decision(const reco::TrackRef& kfTrackRef,
0616                                             reco::PreId& ecalPreId,
0617                                             reco::PreId& hcalPreId,
0618                                             double rho,
0619                                             const reco::BeamSpot& spot,
0620                                             noZS::EcalClusterLazyTools& ecalTools) {
0621   // No implementation currently
0622   return passThrough_;
0623 }
0624 
0625 template <typename CollType>
0626 void LowPtGsfElectronSeedProducer::fillPreIdRefValueMap(edm::Handle<CollType> tracksHandle,
0627                                                         const TrackIndxMap& trksToPreIdIndx,
0628                                                         const edm::OrphanHandle<reco::PreIdCollection>& preIdHandle,
0629                                                         edm::ValueMap<reco::PreIdRef>::Filler& filler) {
0630   std::vector<reco::PreIdRef> values;
0631 
0632   unsigned ntracks = tracksHandle->size();
0633   for (unsigned itrack = 0; itrack < ntracks; ++itrack) {
0634     edm::Ref<CollType> trackRef(tracksHandle, itrack);
0635     auto preIdRefIt = trksToPreIdIndx.find(trackRef.index());
0636     if (preIdRefIt == trksToPreIdIndx.end()) {
0637       values.push_back(reco::PreIdRef());
0638     } else {
0639       edm::Ref<reco::PreIdCollection> preIdRef(preIdHandle, preIdRefIt->second);
0640       values.push_back(preIdRef);
0641     }
0642   }
0643   filler.insert(tracksHandle, values.begin(), values.end());
0644 }
0645 
0646 //////////////////////////////////////////////////////////////////////////////////////////
0647 //
0648 void LowPtGsfElectronSeedProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0649   edm::ParameterSetDescription desc;
0650   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0651   desc.add<edm::InputTag>("pfTracks", edm::InputTag("lowPtGsfElePfTracks"));
0652   desc.add<edm::InputTag>("ecalClusters", edm::InputTag("particleFlowClusterECAL"));
0653   desc.add<edm::InputTag>("hcalClusters", edm::InputTag("particleFlowClusterHCAL"));
0654   desc.add<edm::InputTag>("EBRecHits", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0655   desc.add<edm::InputTag>("EERecHits", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0656   desc.add<edm::InputTag>("rho", edm::InputTag("fixedGridRhoFastjetAllTmp"));
0657   desc.add<edm::InputTag>("BeamSpot", edm::InputTag("offlineBeamSpot"));
0658   desc.add<edm::ESInputTag>("Fitter", edm::ESInputTag{"", "GsfTrajectoryFitter_forPreId"});
0659   desc.add<edm::ESInputTag>("Smoother", edm::ESInputTag{"", "GsfTrajectorySmoother_forPreId"});
0660   desc.add<edm::ESInputTag>("TTRHBuilder", edm::ESInputTag{"", "WithAngleAndTemplate"});
0661   desc.add<std::vector<std::string> >("ModelNames", {});
0662   desc.add<std::vector<std::string> >("ModelWeights", {});
0663   desc.add<std::vector<double> >("ModelThresholds", {});
0664   desc.add<bool>("PassThrough", false);
0665   desc.add<bool>("UsePfTracks", true);
0666   desc.add<double>("MinPtThreshold", 1.0);
0667   desc.add<double>("MaxPtThreshold", 15.);
0668   descriptions.add("lowPtGsfElectronSeeds", desc);
0669 }
0670 
0671 //////////////////////////////////////////////////////////////////////////////////////////
0672 //
0673 #include "FWCore/Framework/interface/MakerMacros.h"
0674 DEFINE_FWK_MODULE(LowPtGsfElectronSeedProducer);
[ 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!