From 1faba747f6a3e07eb9fca9d7febf94d644b00fa3 Mon Sep 17 00:00:00 2001 From: Carlo Battilana Date: Wed, 27 Sep 2017 21:56:03 +0200 Subject: [PATCH 1/6] make TSGForOI phase2 aware --- RecoMuon/TrackerSeedGenerator/plugins/TSGForOI.cc | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/RecoMuon/TrackerSeedGenerator/plugins/TSGForOI.cc b/RecoMuon/TrackerSeedGenerator/plugins/TSGForOI.cc index 02012ec2213fa..338adc259176e 100644 --- a/RecoMuon/TrackerSeedGenerator/plugins/TSGForOI.cc +++ b/RecoMuon/TrackerSeedGenerator/plugins/TSGForOI.cc @@ -58,16 +58,23 @@ void TSGForOI::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { iSetup.get().get(geometry_); iSetup.get().get(estimatorName_,estimator_); iEvent.getByToken(measurementTrackerTag_, measurementTracker_); + + edm::ESHandle tmpTkGeometry; + iSetup.get().get(tmpTkGeometry); + edm::Handle l2TrackCol; iEvent.getByToken(src_, l2TrackCol); - // The product: + //The product: std::unique_ptr > result(new std::vector()); - // Get vector of Detector layers once: std::vector const& tob = measurementTracker_->geometricSearchTracker()->tobLayers(); - std::vector const& tecPositive = measurementTracker_->geometricSearchTracker()->posTecLayers(); - std::vector const& tecNegative = measurementTracker_->geometricSearchTracker()->negTecLayers(); + std::vector const& tecPositive = tmpTkGeometry->isThere(GeomDetEnumerators::P2OTEC) ? + measurementTracker_->geometricSearchTracker()->posTidLayers() : + measurementTracker_->geometricSearchTracker()->posTecLayers(); + std::vector const& tecNegative = tmpTkGeometry->isThere(GeomDetEnumerators::P2OTEC) ? + measurementTracker_->geometricSearchTracker()->negTidLayers() : + measurementTracker_->geometricSearchTracker()->negTecLayers(); // Get the suitable propagators: std::unique_ptr propagatorAlong = SetPropagationDirection(*propagatorAlong_,alongMomentum); From b3658e976cd9e1d8dcf416348128d22bdd7d124b Mon Sep 17 00:00:00 2001 From: Carlo Battilana Date: Wed, 27 Sep 2017 21:59:37 +0200 Subject: [PATCH 2/6] Minor fix: initialize bools --- .../plugins/MaskedMeasurementTrackerEventProducer.cc | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/RecoTracker/MeasurementDet/plugins/MaskedMeasurementTrackerEventProducer.cc b/RecoTracker/MeasurementDet/plugins/MaskedMeasurementTrackerEventProducer.cc index 59d5bb1675378..a5e3ca9227d7d 100644 --- a/RecoTracker/MeasurementDet/plugins/MaskedMeasurementTrackerEventProducer.cc +++ b/RecoTracker/MeasurementDet/plugins/MaskedMeasurementTrackerEventProducer.cc @@ -18,17 +18,18 @@ class dso_hidden MaskedMeasurementTrackerEventProducer final : public edm::strea edm::EDGetTokenT src_; + bool skipClusters_; + bool phase2skipClusters_; + edm::EDGetTokenT maskStrips_; edm::EDGetTokenT maskPixels_; edm::EDGetTokenT maskPhase2OTs_; - - bool skipClusters_; - bool phase2skipClusters_; }; MaskedMeasurementTrackerEventProducer::MaskedMeasurementTrackerEventProducer(const edm::ParameterSet &iConfig) : - src_(consumes(iConfig.getParameter("src"))) + src_(consumes(iConfig.getParameter("src"))), + skipClusters_(false), phase2skipClusters_(false) { //FIXME:temporary solution in order to use this class for both phase0/1 and phase2 if (iConfig.existsAs("clustersToSkip")) { From 92e0e187e37f640188ddd4ca89e1fe01131e814e Mon Sep 17 00:00:00 2001 From: Carlo Battilana Date: Sun, 1 Oct 2017 19:10:03 +0200 Subject: [PATCH 3/6] TrackerHitAssociator from validationFixes_90X --- .../interface/TrackerHitAssociator.h | 11 +- .../src/TrackerHitAssociator.cc | 179 ++++++++++++------ 2 files changed, 124 insertions(+), 66 deletions(-) diff --git a/SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h b/SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h index fb59d31e8177d..be3977b4da1ea 100644 --- a/SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h +++ b/SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h @@ -41,6 +41,7 @@ #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" +#include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h" #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h" #include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h" #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h" @@ -57,9 +58,9 @@ class TrackerHitAssociator { struct Config { Config(const edm::ParameterSet& conf, edm::ConsumesCollector && iC); Config(edm::ConsumesCollector && iC); - bool doPixel_, doStrip_, doTrackAssoc_, assocHitbySimTrack_; + bool doPixel_, doStrip_, useOTph2_, doTrackAssoc_, assocHitbySimTrack_; edm::EDGetTokenT > stripToken_; - edm::EDGetTokenT > pixelToken_; + edm::EDGetTokenT > pixelToken_, ph2OTrToken_; std::vector > > cfTokens_; std::vector > > simHitTokens_; }; @@ -92,7 +93,8 @@ class TrackerHitAssociator { std::vector associateMatchedRecHit(const SiStripMatchedRecHit2D * matchedrechit, std::vector* simhitCFPos=0) const; std::vector associateProjectedRecHit(const ProjectedSiStripRecHit2D * projectedrechit, std::vector* simhitCFPos=0) const; - void associatePixelRecHit(const SiPixelRecHit * pixelrechit, std::vector & simhitid, std::vector* simhitCFPos=0) const; + void associatePhase2TrackerRecHit(const Phase2TrackerRecHit1D* rechit, std::vector & simtrackid, std::vector* simhitCFPos=0) const; + void associatePixelRecHit(const SiPixelRecHit * pixelrechit, std::vector & simtrackid, std::vector* simhitCFPos=0) const; std::vector associateFastRecHit(const FastTrackerRecHit * rechit) const; std::vector associateMultiRecHitId(const SiTrackerMultiRecHit * multirechit, std::vector* simhitCFPos=0) const; std::vector associateMultiRecHit(const SiTrackerMultiRecHit * multirechit) const; @@ -109,7 +111,8 @@ class TrackerHitAssociator { void makeMaps(const edm::Event& theEvent, const Config& config); edm::Handle< edm::DetSetVector > stripdigisimlink; edm::Handle< edm::DetSetVector > pixeldigisimlink; - bool doPixel_, doStrip_, doTrackAssoc_, assocHitbySimTrack_; + edm::Handle< edm::DetSetVector > ph2trackerdigisimlink; + bool doPixel_, doStrip_, useOTph2_, doTrackAssoc_, assocHitbySimTrack_; }; #endif diff --git a/SimTracker/TrackerHitAssociation/src/TrackerHitAssociator.cc b/SimTracker/TrackerHitAssociation/src/TrackerHitAssociator.cc index 6366e4f04f428..64d91a5e4fd30 100644 --- a/SimTracker/TrackerHitAssociation/src/TrackerHitAssociator.cc +++ b/SimTracker/TrackerHitAssociation/src/TrackerHitAssociator.cc @@ -26,11 +26,18 @@ using namespace edm; TrackerHitAssociator::Config::Config(edm::ConsumesCollector && iC) : doPixel_(true), doStrip_(true), + useOTph2_(false), doTrackAssoc_(false), assocHitbySimTrack_(false) { - if(doStrip_) stripToken_ = iC.consumes >(edm::InputTag("simSiStripDigis")); - if(doPixel_) pixelToken_ = iC.consumes >(edm::InputTag("simSiPixelDigis")); + if(doStrip_) { + if (useOTph2_) ph2OTrToken_ = iC.consumes >(edm::InputTag("simSiPixelDigis","Tracker")); + else stripToken_ = iC.consumes >(edm::InputTag("simSiStripDigis")); + } + if(doPixel_) { + if (useOTph2_) pixelToken_ = iC.consumes >(edm::InputTag("simSiPixelDigis","Pixel")); + else pixelToken_ = iC.consumes >(edm::InputTag("simSiPixelDigis")); + } if(!doTrackAssoc_) { std::vector trackerContainers; trackerContainers.reserve(12); @@ -61,10 +68,15 @@ TrackerHitAssociator::Config::Config(edm::ConsumesCollector && iC) : TrackerHitAssociator::Config::Config(const edm::ParameterSet& conf, edm::ConsumesCollector && iC) : doPixel_( conf.getParameter("associatePixel") ), doStrip_( conf.getParameter("associateStrip") ), + useOTph2_( conf.existsAs("usePhase2Tracker") ? conf.getParameter("usePhase2Tracker") : false), + // doTrackAssoc_( conf.getParameter("associateRecoTracks") ), assocHitbySimTrack_(conf.existsAs("associateHitbySimTrack") ? conf.getParameter("associateHitbySimTrack") : false) { - if(doStrip_) stripToken_ = iC.consumes >(conf.getParameter("stripSimLinkSrc")); + if(doStrip_) { + if (useOTph2_) ph2OTrToken_ = iC.consumes >(conf.getParameter("phase2TrackerSimLinkSrc")); + else stripToken_ = iC.consumes >(conf.getParameter("stripSimLinkSrc")); + } if(doPixel_) pixelToken_ = iC.consumes >(conf.getParameter("pixelSimLinkSrc")); if(!doTrackAssoc_) { std::vector trackerContainers(conf.getParameter >("ROUList")); @@ -83,6 +95,7 @@ TrackerHitAssociator::Config::Config(const edm::ParameterSet& conf, edm::Consume TrackerHitAssociator::TrackerHitAssociator(const edm::Event& e, const TrackerHitAssociator::Config& config) : doPixel_(config.doPixel_), doStrip_(config.doStrip_), + useOTph2_(config.useOTph2_), doTrackAssoc_(config.doTrackAssoc_), assocHitbySimTrack_(config.assocHitbySimTrack_) { //if track association there is no need to access the input collections @@ -90,7 +103,10 @@ TrackerHitAssociator::TrackerHitAssociator(const edm::Event& e, const TrackerHit makeMaps(e, config); } - if(doStrip_) e.getByToken(config.stripToken_, stripdigisimlink); + if(doStrip_) { + if (useOTph2_) e.getByToken(config.ph2OTrToken_, ph2trackerdigisimlink); + else e.getByToken(config.stripToken_, stripdigisimlink); + } if(doPixel_) e.getByToken(config.pixelToken_, pixeldigisimlink); } @@ -325,65 +341,38 @@ std::vector< SimHitIdpr > TrackerHitAssociator::associateHitId(const TrackingRec void TrackerHitAssociator::associateHitId(const TrackingRecHit & thit, std::vector< SimHitIdpr > & simtkid, std::vector* simhitCFPos) const { - - simtkid.clear(); + simtkid.clear(); - //get the Detector type of the rechit - DetId detid= thit.geographicalId(); - if (const SiTrackerMultiRecHit * rechit = dynamic_cast(&thit)){ - simtkid=associateMultiRecHitId(rechit, simhitCFPos); - } + if (const SiTrackerMultiRecHit * rechit = dynamic_cast(&thit)) + simtkid = associateMultiRecHitId(rechit, simhitCFPos); - // cout << "Associator ---> get Detid " << detID << endl; - //check we are in the strip tracker - if(detid.subdetId() == StripSubdetector::TIB || - detid.subdetId() == StripSubdetector::TOB || - detid.subdetId() == StripSubdetector::TID || - detid.subdetId() == StripSubdetector::TEC) - { - //check if it is a simple SiStripRecHit2D - if(const SiStripRecHit2D * rechit = - dynamic_cast(&thit)) - { - associateSiStripRecHit(rechit, simtkid, simhitCFPos); - } - //check if it is a SiStripRecHit1D - else if(const SiStripRecHit1D * rechit = - dynamic_cast(&thit)) - { - associateSiStripRecHit(rechit, simtkid, simhitCFPos); - } - //check if it is a SiStripMatchedRecHit2D - else if(const SiStripMatchedRecHit2D * rechit = - dynamic_cast(&thit)) - { - simtkid = associateMatchedRecHit(rechit, simhitCFPos); - } - //check if it is a ProjectedSiStripRecHit2D - else if(const ProjectedSiStripRecHit2D * rechit = - dynamic_cast(&thit)) - { - simtkid = associateProjectedRecHit(rechit, simhitCFPos); - } - else{ - //std::cout << "associate to invalid" << std::endl; - //throw cms::Exception("Unknown RecHit Type") << "TrackerHitAssociator failed second casting of " << typeid(thit).name() << " type "; - } - } - //check we are in the pixel tracker - else if( (unsigned int)(detid.subdetId()) == PixelSubdetector::PixelBarrel || - (unsigned int)(detid.subdetId()) == PixelSubdetector::PixelEndcap) - { - if(const SiPixelRecHit * rechit = dynamic_cast(&thit)) - { - associatePixelRecHit(rechit, simtkid, simhitCFPos); - } - } - //check if these are GSRecHits (from FastSim) - if(trackerHitRTTI::isFast(thit)) - { - simtkid = associateFastRecHit(static_cast(&thit)); - } + //check if it is a simple SiStripRecHit2D + if (const SiStripRecHit2D * rechit = dynamic_cast(&thit)) + associateSiStripRecHit(rechit, simtkid, simhitCFPos); + + //check if it is a SiStripRecHit1D + else if(const SiStripRecHit1D * rechit = dynamic_cast(&thit)) + associateSiStripRecHit(rechit, simtkid, simhitCFPos); + + //check if it is a SiStripMatchedRecHit2D + else if(const SiStripMatchedRecHit2D * rechit = dynamic_cast(&thit)) + simtkid = associateMatchedRecHit(rechit, simhitCFPos); + + //check if it is a ProjectedSiStripRecHit2D + else if(const ProjectedSiStripRecHit2D * rechit = dynamic_cast(&thit)) + simtkid = associateProjectedRecHit(rechit, simhitCFPos); + + //check if it is a Phase2TrackerRecHit1D + else if(const Phase2TrackerRecHit1D * rechit = dynamic_cast(&thit)) + associatePhase2TrackerRecHit(rechit, simtkid, simhitCFPos); + + //check if it is a SiPixelRecHit + else if(const SiPixelRecHit * rechit = dynamic_cast(&thit)) + associatePixelRecHit(rechit, simtkid, simhitCFPos); + + //check if these are GSRecHits (from FastSim) + if(trackerHitRTTI::isFast(thit)) + simtkid = associateFastRecHit(static_cast(&thit)); } template @@ -533,7 +522,6 @@ std::vector TrackerHitAssociator::associateMatchedRecHit(const SiSt return simtrackid; } - std::vector TrackerHitAssociator::associateProjectedRecHit(const ProjectedSiStripRecHit2D * projectedrechit, std::vector* simhitCFPos) const { @@ -546,6 +534,73 @@ std::vector TrackerHitAssociator::associateProjectedRecHit(const Pr return matched_mono; } +void TrackerHitAssociator::associatePhase2TrackerRecHit(const Phase2TrackerRecHit1D* rechit, + std::vector & simtrackid, + std::vector* simhitCFPos) const +{ + // + // Phase 2 outer tracker associator + // + DetId detid= rechit->geographicalId(); + uint32_t detID = detid.rawId(); + + edm::DetSetVector::const_iterator isearch = ph2trackerdigisimlink->find(detID); + if(isearch != ph2trackerdigisimlink->end()) { //if it is not empty + edm::DetSet link_detset = (*isearch); + Phase2TrackerRecHit1D::CluRef const& cluster = rechit->cluster(); + + //check the reference is valid + + if(!(cluster.isNull())){//if the cluster is valid + int minRow = (*cluster).firstStrip(); + int maxRow = (*cluster).firstStrip() + (*cluster).size(); + int Col = (*cluster).column(); + // std::cout << " Cluster minRow " << minRow << " maxRow " << maxRow << " column " << Col << std::endl; + edm::DetSet::const_iterator linkiter = link_detset.data.begin(), linkEnd = link_detset.data.end(); + int dsl = 0; + std::vector idcachev; + std::vector CFposcachev; + for( ; linkiter != linkEnd; ++linkiter) { + ++dsl; + std::pair coord = Phase2TrackerDigi::channelToPixel(linkiter->channel()); + // std::cout << " " << dsl << ") Digi link: row " << pixel_coord.first << " col " << pixel_coord.second << std::endl; + if( coord.first <= maxRow && + coord.first >= minRow && + coord.second == Col ) { + // std::cout << " !-> trackid " << linkiter->SimTrackId() << endl; + // std::cout << " fraction " << linkiter->fraction() << endl; + SimHitIdpr currentId(linkiter->SimTrackId(), linkiter->eventId()); + if(find(idcachev.begin(),idcachev.end(),currentId) == idcachev.end()){ + simtrackid.push_back(currentId); + idcachev.push_back(currentId); + } + + if (simhitCFPos != 0) { + //create a vector that contains all the position (in the MixCollection) of the SimHits that contributed to the RecHit + //write position only once + unsigned int currentCFPos = linkiter->CFposition(); + unsigned int tofBin = linkiter->TofBin(); + subDetTofBin theSubDetTofBin = std::make_pair(detid.subdetId(), tofBin); + simhit_collectionMap::const_iterator it = SimHitCollMap.find(theSubDetTofBin); + if (it!= SimHitCollMap.end()) { + simhitAddr currentAddr = std::make_pair(it->second, currentCFPos); + if(find(CFposcachev.begin(), CFposcachev.end(), currentAddr) == CFposcachev.end()) { + CFposcachev.push_back(currentAddr); + simhitCFPos->push_back(currentAddr); + } + } + } + + } + } // end of simlink loop + } + else{ + edm::LogError("TrackerHitAssociator")<<"no Pixel cluster reference attached"; + + } + } +} + void TrackerHitAssociator::associatePixelRecHit(const SiPixelRecHit * pixelrechit, std::vector & simtrackid, std::vector* simhitCFPos) const From 9fafe5e136b8fbc9a2a004e673917c6ecda50655 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Fri, 11 Aug 2017 11:55:41 +0100 Subject: [PATCH 4/6] backporting 92X pixel matching to 903 --- .../EgammaReco/interface/ElectronSeed.h | 203 +++++--- .../EgammaReco/interface/ElectronSeedFwd.h | 2 +- DataFormats/EgammaReco/src/ElectronSeed.cc | 261 ++++++---- DataFormats/EgammaReco/src/classes.h | 4 +- DataFormats/EgammaReco/src/classes_def.xml | 13 +- .../interface/TrajSeedMatcher.h | 301 +++++++++++ .../src/ElectronSeedGenerator.cc | 3 +- .../src/TrajSeedMatcher.cc | 486 ++++++++++++++++++ .../plugins/ElectronNHitSeedProducer.cc | 179 +++++++ ...ckingRegionsFromSuperClustersEDProducer.cc | 9 + ...TrackingRegionsFromSuperClustersProducer.h | 314 +++++++++++ .../EgammaHLTFilteredSuperClusterProducer.cc | 168 ++++++ .../src/EgammaHLTPixelMatchVarProducer.cc | 164 +++++- 13 files changed, 1926 insertions(+), 181 deletions(-) create mode 100644 RecoEgamma/EgammaElectronAlgos/interface/TrajSeedMatcher.h create mode 100644 RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc create mode 100644 RecoEgamma/EgammaElectronProducers/plugins/ElectronNHitSeedProducer.cc create mode 100644 RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersEDProducer.cc create mode 100644 RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersProducer.h create mode 100644 RecoEgamma/EgammaHLTProducers/src/EgammaHLTFilteredSuperClusterProducer.cc diff --git a/DataFormats/EgammaReco/interface/ElectronSeed.h b/DataFormats/EgammaReco/interface/ElectronSeed.h index 173150a215332..2c6098d1ec2d1 100644 --- a/DataFormats/EgammaReco/interface/ElectronSeed.h +++ b/DataFormats/EgammaReco/interface/ElectronSeed.h @@ -1,14 +1,38 @@ -#ifndef ElectronSeed_h -#define ElectronSeed_h 1 - -/** \class reco::ElectronSeed - * - * ElectronSeed is a seed for gsf tracking, constructed from - * either a supercluster or a ctf track. - * - * \author D.Chamont, U.Berthon, C.Charlot, LLR Palaiseau - * - ************************************************************/ +#ifndef DataFormats_EgammaReco_ElectronSeed_h +#define DataFormats_EgammaReco_ElectronSeed_h + +//******************************************************************** +// +// A verson of reco::ElectronSeed which can have N hits as part of the +// 2017 upgrade of E/gamma pixel matching for the phaseI pixels +// +// author: S. Harper (RAL), 2017 +// +//notes: +// While it is technically named ElectronSeed, it is effectively a new class +// However to simplify things, the name ElectronSeed was kept +// (trust me it was simplier...) +// +// Noticed that h/e values never seem to used anywhere and they are a +// mild pain to propagate in the new framework so they were removed +// +// infinities are used to mark invalid unset values to maintain +// compatibilty with the orginal ElectronSeed class +// +//description: +// An ElectronSeed is a TrajectorySeed with E/gamma specific information +// A TrajectorySeed has a series of hits associated with it +// (accessed by TrajectorySeed::nHits(), TrajectorySeed::recHits()) +// and ElectronSeed stores which of those hits match well to a supercluster +// together with the matching parameters (this is known as EcalDriven). +// ElectronSeeds can be TrackerDriven in which case the matching is not done. +// It used to be fixed to two matched hits, now this is an arbitary number +// Its designed with pixel matching with mind but tries to be generally +// applicable to strips as well. +// It is worth noting that due to different ways ElectronSeeds can be created +// they do not always have all parameters filled +// +//******************************************************************** #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h" #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h" @@ -23,93 +47,124 @@ #include namespace reco - { - -class ElectronSeed : public TrajectorySeed - { +{ + + class ElectronSeed : public TrajectorySeed { public : - + struct PMVars { + float dRZPos; + float dRZNeg; + float dPhiPos; + float dPhiNeg; + int detId; //this is already stored as the hit is stored in traj seed but a useful sanity check + int layerOrDiskNr;//redundant as stored in detId but its a huge pain to hence why its saved here + + PMVars(); + + void setDPhi(float pos,float neg); + void setDRZ(float pos,float neg); + void setDet(int iDetId,int iLayerOrDiskNr); + + }; + typedef edm::OwnVector RecHitContainer ; typedef edm::RefToBase CaloClusterRef ; typedef edm::Ref CtfTrackRef ; - + static std::string const & name() - { + { static std::string const name_("ElectronSeed") ; return name_; - } - + } + //! Construction of base attributes ElectronSeed() ; ElectronSeed( const TrajectorySeed & ) ; ElectronSeed( PTrajectoryStateOnDet & pts, RecHitContainer & rh, PropagationDirection & dir ) ; ElectronSeed * clone() const { return new ElectronSeed(*this) ; } - virtual ~ElectronSeed() ; + virtual ~ElectronSeed(); //! Set additional info void setCtfTrack( const CtfTrackRef & ) ; - void setCaloCluster - ( const CaloClusterRef &, - unsigned char hitsMask =0, - int subDet2 =0, int subDet1 =0, - float hoe1 =std::numeric_limits::infinity(), - float hoe2 =std::numeric_limits::infinity() ) ; - void setNegAttributes - ( float dRz2 =std::numeric_limits::infinity(), - float dPhi2 =std::numeric_limits::infinity(), - float dRz1 =std::numeric_limits::infinity(), - float dPhi1 =std::numeric_limits::infinity() ) ; - void setPosAttributes - ( float dRz2 =std::numeric_limits::infinity(), - float dPhi2 =std::numeric_limits::infinity(), - float dRz1 =std::numeric_limits::infinity(), - float dPhi1 =std::numeric_limits::infinity() ) ; - + void setCaloCluster( const CaloClusterRef& clus){caloCluster_=clus;isEcalDriven_=true;} + void addHitInfo(const PMVars& hitVars){hitInfo_.push_back(hitVars);} + void setNrLayersAlongTraj(int val){nrLayersAlongTraj_=val;} //! Accessors const CtfTrackRef& ctfTrack() const { return ctfTrack_ ; } const CaloClusterRef& caloCluster() const { return caloCluster_ ; } - unsigned char hitsMask() const { return hitsMask_ ; } - int subDet2() const { return subDet2_ ; } - float dRz2() const { return dRz2_ ; } - float dPhi2() const { return dPhi2_ ; } - float dRz2Pos() const { return dRz2Pos_ ; } - float dPhi2Pos() const { return dPhi2Pos_ ; } - int subDet1() const { return subDet1_ ; } - float dRz1() const { return dRz1_ ; } - float dPhi1() const { return dPhi1_ ; } - float dRz1Pos() const { return dRz1Pos_ ; } - float dPhi1Pos() const { return dPhi1Pos_ ; } - float hoe1() const { return hcalDepth1OverEcal_ ; } - float hoe2() const { return hcalDepth2OverEcal_ ; } - + //! Utility TrackCharge getCharge() const { return startingState().parameters().charge() ; } bool isEcalDriven() const { return isEcalDriven_ ; } bool isTrackerDriven() const { return isTrackerDriven_ ; } + const std::vector& hitInfo()const{return hitInfo_;} + float dPhiNeg(size_t hitNr)const{return getVal(hitNr,&PMVars::dPhiNeg);} + float dPhiPos(size_t hitNr)const{return getVal(hitNr,&PMVars::dPhiPos);} + float dPhiBest(size_t hitNr)const{return bestVal(dPhiNeg(hitNr),dPhiPos(hitNr));} + float dRZPos(size_t hitNr)const{return getVal(hitNr,&PMVars::dRZPos);} + float dRZNeg(size_t hitNr)const{return getVal(hitNr,&PMVars::dRZNeg);} + float dRZBest(size_t hitNr)const{return bestVal(dRZNeg(hitNr),dRZPos(hitNr));} + int detId(size_t hitNr)const{return hitNr::infinity(), + const float dPhi2=std::numeric_limits::infinity(), + const float dRZ1=std::numeric_limits::infinity(), + const float dPhi1=std::numeric_limits::infinity()); + void setPosAttributes(const float dRZ2=std::numeric_limits::infinity(), + const float dPhi2=std::numeric_limits::infinity(), + const float dRZ1=std::numeric_limits::infinity(), + const float dPhi1=std::numeric_limits::infinity()); + + //this is a backwards compatible function designed to + //convert old format ElectronSeeds to the new format + //only public due to root io rules, not intended for any other use + //also in theory not necessary to part of this class + static std::vector createHitInfo(const float dPhi1Pos,const float dPhi1Neg, + const float dRZ1Pos,const float dRZ1Neg, + const float dPhi2Pos,const float dPhi2Neg, + const float dRZ2Pos,const float dRZ2Neg, + const char hitMask,const TrajectorySeed::range recHits); private: - - CtfTrackRef ctfTrack_ ; - CaloClusterRef caloCluster_ ; - unsigned char hitsMask_ ; - int subDet2_ ; - float dRz2_ ; - float dPhi2_ ; - float dRz2Pos_ ; - float dPhi2Pos_ ; - int subDet1_ ; - float dRz1_ ; - float dPhi1_ ; - float dRz1Pos_ ; - float dPhi1Pos_ ; - float hcalDepth1OverEcal_ ; // hcal over ecal seed cluster energy using first hcal depth - float hcalDepth2OverEcal_ ; // hcal over ecal seed cluster energy using 2nd hcal depth - bool isEcalDriven_ ; - bool isTrackerDriven_ ; - - } ; - - } + static float bestVal(float val1,float val2){return std::abs(val1) + T getVal(unsigned int hitNr,T PMVars::*val)const{ + return hitNr::infinity(); + } + static std::vector hitNrsFromMask(unsigned int hitMask); + + private: + CtfTrackRef ctfTrack_; + CaloClusterRef caloCluster_; + std::vector hitInfo_; + int nrLayersAlongTraj_; + + bool isEcalDriven_; + bool isTrackerDriven_; + + }; +} #endif diff --git a/DataFormats/EgammaReco/interface/ElectronSeedFwd.h b/DataFormats/EgammaReco/interface/ElectronSeedFwd.h index 9f18dbe8a6509..8562974255b66 100644 --- a/DataFormats/EgammaReco/interface/ElectronSeedFwd.h +++ b/DataFormats/EgammaReco/interface/ElectronSeedFwd.h @@ -17,7 +17,7 @@ namespace reco { /// vector of objects in the same collection of ElectronSeed objects typedef edm::RefVector ElectronSeedRefVector; /// iterator over a vector of reference to ElectronSeed objects - typedef ElectronSeedRefVector::iterator electronephltseed_iterator; + typedef ElectronSeedRefVector::iterator electronseed_iterator; } #endif diff --git a/DataFormats/EgammaReco/src/ElectronSeed.cc b/DataFormats/EgammaReco/src/ElectronSeed.cc index d4b758ced108a..176a4d368a204 100644 --- a/DataFormats/EgammaReco/src/ElectronSeed.cc +++ b/DataFormats/EgammaReco/src/ElectronSeed.cc @@ -1,105 +1,190 @@ #include "DataFormats/EgammaReco/interface/ElectronSeed.h" -#include "DataFormats/TrackReco/interface/Track.h" + +#include using namespace reco ; ElectronSeed::ElectronSeed() - : TrajectorySeed(), ctfTrack_(), caloCluster_(), hitsMask_(0), - subDet2_(0), - dRz2_(std::numeric_limits::infinity()), - dPhi2_(std::numeric_limits::infinity()), - dRz2Pos_(std::numeric_limits::infinity()), - dPhi2Pos_(std::numeric_limits::infinity()), - subDet1_(0), - dRz1_(std::numeric_limits::infinity()), - dPhi1_(std::numeric_limits::infinity()), - dRz1Pos_(std::numeric_limits::infinity()), - dPhi1Pos_(std::numeric_limits::infinity()), - hcalDepth1OverEcal_(std::numeric_limits::infinity()), - hcalDepth2OverEcal_(std::numeric_limits::infinity()), - isEcalDriven_(false), isTrackerDriven_(false) - {} + : TrajectorySeed(), ctfTrack_(), caloCluster_(), hitInfo_(), + nrLayersAlongTraj_(0), + isEcalDriven_(false), isTrackerDriven_(false) + +{} ElectronSeed::ElectronSeed - ( const TrajectorySeed & seed ) - : TrajectorySeed(seed), - ctfTrack_(), caloCluster_(), hitsMask_(0), - subDet2_(0), - dRz2_(std::numeric_limits::infinity()), - dPhi2_(std::numeric_limits::infinity()), - dRz2Pos_(std::numeric_limits::infinity()), - dPhi2Pos_(std::numeric_limits::infinity()), - subDet1_(0), - dRz1_(std::numeric_limits::infinity()), - dPhi1_(std::numeric_limits::infinity()), - dRz1Pos_(std::numeric_limits::infinity()), - dPhi1Pos_(std::numeric_limits::infinity()), - hcalDepth1OverEcal_(std::numeric_limits::infinity()), - hcalDepth2OverEcal_(std::numeric_limits::infinity()), - isEcalDriven_(false), isTrackerDriven_(false) - {} +( const TrajectorySeed & seed ) + : TrajectorySeed(seed), + ctfTrack_(), caloCluster_(), hitInfo_(), + nrLayersAlongTraj_(0), + isEcalDriven_(false), isTrackerDriven_(false) +{} ElectronSeed::ElectronSeed ( PTrajectoryStateOnDet & pts, recHitContainer & rh, PropagationDirection & dir ) - : TrajectorySeed(pts,rh,dir), - ctfTrack_(), caloCluster_(), hitsMask_(0), - subDet2_(0), - dRz2_(std::numeric_limits::infinity()), - dPhi2_(std::numeric_limits::infinity()), - dRz2Pos_(std::numeric_limits::infinity()), - dPhi2Pos_(std::numeric_limits::infinity()), - subDet1_(0), - dRz1_(std::numeric_limits::infinity()), - dPhi1_(std::numeric_limits::infinity()), - dRz1Pos_(std::numeric_limits::infinity()), - dPhi1Pos_(std::numeric_limits::infinity()), - hcalDepth1OverEcal_(std::numeric_limits::infinity()), - hcalDepth2OverEcal_(std::numeric_limits::infinity()), - isEcalDriven_(false), isTrackerDriven_(false) - {} + : TrajectorySeed(pts,rh,dir), + ctfTrack_(), caloCluster_(), hitInfo_(), + nrLayersAlongTraj_(0), + isEcalDriven_(false), isTrackerDriven_(false) +{} + +ElectronSeed::~ElectronSeed()=default; void ElectronSeed::setCtfTrack - ( const CtfTrackRef & ctfTrack ) - { +( const CtfTrackRef & ctfTrack ) +{ ctfTrack_ = ctfTrack ; isTrackerDriven_ = true ; - } - -void ElectronSeed::setCaloCluster - ( const CaloClusterRef & scl, - unsigned char hitsMask, - int subDet2, int subDet1, - float hoe1, float hoe2 ) - { - caloCluster_ = scl ; - hitsMask_ = hitsMask ; - isEcalDriven_ = true ; - subDet2_ = subDet2 ; - subDet1_ = subDet1 ; - hcalDepth1OverEcal_ = hoe1 ; - hcalDepth2OverEcal_ = hoe2 ; - } - -void ElectronSeed::setNegAttributes - ( float dRz2, float dPhi2, float dRz1, float dPhi1 ) - { - dRz2_ = dRz2 ; - dPhi2_ = dPhi2 ; - dRz1_ = dRz1 ; - dPhi1_ = dPhi1 ; - } - -void ElectronSeed::setPosAttributes - ( float dRz2, float dPhi2, float dRz1, float dPhi1 ) - { - dRz2Pos_ = dRz2 ; - dPhi2Pos_ = dPhi2 ; - dRz1Pos_ = dRz1 ; - dPhi1Pos_ = dPhi1 ; - } - -ElectronSeed::~ElectronSeed() - {} +} + +//the hit mask tells us which hits were used in the seed +//typically all are used at the HLT but this could change in the future +//RECO only uses some of them +unsigned int ElectronSeed::hitsMask()const +{ + int mask=0; + for(size_t hitNr=0;hitNrgeographicalId().rawId(); + auto detIdMatcher = [hitDetId](const ElectronSeed::PMVars& var){return hitDetId==var.detId;}; + if(std::find_if(hitInfo_.begin(),hitInfo_.end(),detIdMatcher)!=hitInfo_.end()){ + mask|=bitNr; + } + } + return mask; +} + +void ElectronSeed::initTwoHitSeed(const unsigned char hitMask) +{ + hitInfo_.resize(2); + + std::vector hitNrs = hitNrsFromMask(hitMask); + if(hitNrs.size()!=2){ + throw cms::Exception("LogicError") + <<"in ElectronSeed::"<<__FUNCTION__<<","<<__LINE__ + <<": number of hits in hit mask is "<(hitMask)<=nHits() || hitNrs[1]>=nHits()){ + throw cms::Exception("LogicError") + <<"in ElectronSeed::"<<__FUNCTION__<<","<<__LINE__ + <<": hits are "<(hitMask)<::infinity(),std::numeric_limits::infinity()); + info.setDRZ(std::numeric_limits::infinity(),std::numeric_limits::infinity()); + info.setDet((recHits().first+hitNrs[hitNr])->geographicalId(),-1); + } + +} + +void ElectronSeed::setNegAttributes(float dRZ2,float dPhi2,float dRZ1,float dPhi1) +{ + if(hitInfo_.size()!=2){ + throw cms::Exception("LogicError") <<"ElectronSeed::setNegAttributes should only operate on seeds with exactly two hits. This is because it is a legacy function to preverse backwards compatiblity and should not be used on new code which matches variable number of hits"; + } + hitInfo_[0].dRZNeg = dRZ1; + hitInfo_[1].dRZNeg = dRZ2; + hitInfo_[0].dPhiNeg = dPhi1; + hitInfo_[1].dPhiNeg = dPhi2; + +} + +void ElectronSeed::setPosAttributes(float dRZ2,float dPhi2,float dRZ1,float dPhi1) +{ + if(hitInfo_.size()!=2){ + throw cms::Exception("LogicError") <<"ElectronSeed::setPosAttributes should only operate on seeds with exactly two hits. This is because it is a legacy function to preverse backwards compatiblity and should not be used on new code which matches variable number of hits"; + } + hitInfo_[0].dRZPos = dRZ1; + hitInfo_[1].dRZPos = dRZ2; + hitInfo_[0].dPhiPos = dPhi1; + hitInfo_[1].dPhiPos = dPhi2; +} + +std::vector +ElectronSeed::hitNrsFromMask(unsigned int hitMask) +{ + std::vector hitNrs; + for(size_t bitNr=0; bitNr +ElectronSeed::createHitInfo(const float dPhi1Pos,const float dPhi1Neg, + const float dRZ1Pos,const float dRZ1Neg, + const float dPhi2Pos,const float dPhi2Neg, + const float dRZ2Pos,const float dRZ2Neg, + const char hitMask,const TrajectorySeed::range recHits) +{ + if(hitMask==0) return std::vector(); //was trackerDriven so no matched hits + + size_t nrRecHits = std::distance(recHits.first,recHits.second); + std::vector hitNrs = hitNrsFromMask(hitMask); + + if(hitNrs.size()!=2){ + throw cms::Exception("LogicError") + <<"in ElectronSeed::"<<__FUNCTION__<<","<<__LINE__ + <<": number of hits in hit mask is "<(hitMask)<=nrRecHits || hitNrs[1]>=nrRecHits){ + throw cms::Exception("LogicError") + <<"in ElectronSeed::"<<__FUNCTION__<<","<<__LINE__ + <<": hits are "<(hitMask)< hitInfo(2); + hitInfo[0].setDPhi(dPhi1Pos,dPhi1Neg); + hitInfo[0].setDRZ(dRZ1Pos,dRZ1Neg); + hitInfo[0].setDet((recHits.first+hitNrs[0])->geographicalId(),-1); //getting the layer information needs tracker topo, hence why its stored in the first as its a pain to access + hitInfo[1].setDPhi(dPhi2Pos,dPhi2Neg); + hitInfo[1].setDRZ(dRZ2Pos,dRZ2Neg); + hitInfo[1].setDet((recHits.first+hitNrs[1])->geographicalId(),-1); //getting the layer information needs tracker topo, hence why its stored in the first as its a pain to access + return hitInfo; + +} + + +ElectronSeed::PMVars::PMVars(): + dRZPos(std::numeric_limits::infinity()), + dRZNeg(std::numeric_limits::infinity()), + dPhiPos(std::numeric_limits::infinity()), + dPhiNeg(std::numeric_limits::infinity()), + detId(0), + layerOrDiskNr(-1) +{} + +void ElectronSeed::PMVars::setDPhi(float pos,float neg) +{ + dPhiPos = pos; + dPhiNeg = neg; +} + +void ElectronSeed::PMVars::setDRZ(float pos,float neg) +{ + dRZPos = pos; + dRZNeg = neg; +} + +void ElectronSeed::PMVars::setDet(int iDetId,int iLayerOrDiskNr) +{ + detId = iDetId; + layerOrDiskNr = iLayerOrDiskNr; +} + + + diff --git a/DataFormats/EgammaReco/src/classes.h b/DataFormats/EgammaReco/src/classes.h index e52a650d1f337..d4122b8050142 100644 --- a/DataFormats/EgammaReco/src/classes.h +++ b/DataFormats/EgammaReco/src/classes.h @@ -43,8 +43,8 @@ namespace DataFormats_EgammaReco { edm::Ref r3; edm::RefProd rp3; edm::Wrapper > wrv3; - - + std::vector,reco::SuperCluster,edm::refhelper::FindUsingAdvance,reco::SuperCluster> > >vr3; + edm::Wrapper,reco::SuperCluster,edm::refhelper::FindUsingAdvance,reco::SuperCluster> > > > wvr3; reco::EgammaTriggerCollection v4; edm::Wrapper w4; edm::Ref r4; diff --git a/DataFormats/EgammaReco/src/classes_def.xml b/DataFormats/EgammaReco/src/classes_def.xml index 5dcd4ac4c0173..a89006b5deb8a 100644 --- a/DataFormats/EgammaReco/src/classes_def.xml +++ b/DataFormats/EgammaReco/src/classes_def.xml @@ -34,6 +34,8 @@ + + @@ -86,11 +88,19 @@ - + + + + recHits());]]> + + + + + @@ -103,7 +113,6 @@ - diff --git a/RecoEgamma/EgammaElectronAlgos/interface/TrajSeedMatcher.h b/RecoEgamma/EgammaElectronAlgos/interface/TrajSeedMatcher.h new file mode 100644 index 0000000000000..57798b54e1743 --- /dev/null +++ b/RecoEgamma/EgammaElectronAlgos/interface/TrajSeedMatcher.h @@ -0,0 +1,301 @@ +#ifndef RecoEgamma_EgammaElectronAlgos_TrajSeedMatcher_h +#define RecoEgamma_EgammaElectronAlgos_TrajSeedMatcher_h + + +//****************************************************************************** +// +// Part of the refactorisation of of the E/gamma pixel matching for 2017 pixels +// This refactorisation converts the monolithic approach to a series of +// independent producer modules, with each modules performing a specific +// job as recommended by the 2017 tracker framework +// +// +// The module is based of PixelHitMatcher (the seed based functions) but +// extended to match on an arbitary number of hits rather than just doublets. +// It is also aware of how many layers the supercluster trajectory passed through +// and uses that information to determine how many hits to require +// Other than that, its a direct port and follows what PixelHitMatcher did +// +// +// Author : Sam Harper (RAL), 2017 +// +//******************************************************************************* + + + + +#include "FWCore/Framework/interface/ESHandle.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" + +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h" +#include "TrackingTools/DetLayers/interface/NavigationSchool.h" +#include "TrackingTools/RecoGeometry/interface/GlobalDetLayerGeometry.h" +#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" + +#include + +namespace edm{ + class EventSetup; + class ConfigurationDescriptions; + class ParameterSet; + class ParameterSetDescription; +} + +class FreeTrajectoryState; +class TrackingRecHit; + + +//stolen from PixelHitMatcher +//decide if its evil or not later +//actually I think the answer is, yes, yes its evil +//maybe replace with less evil? +namespace std{ + template<> + struct hash > { + std::size_t operator()(const std::pair& g) const { + auto h1 = std::hash()((unsigned long long)g.first); + unsigned long long k; memcpy(&k, &g.second,sizeof(k)); + auto h2 = std::hash()(k); + return h1 ^ (h2 << 1); + } + }; + template<> + struct equal_to > : public std::binary_function,std::pair,bool> { + bool operator()(const std::pair& a, + const std::pair& b) const { + return (a.first == b.first) & (a.second == b.second); + } + }; +} + +class TrajSeedMatcher { +public: + class SCHitMatch { + public: + SCHitMatch():detId_(0), + dRZ_(std::numeric_limits::max()), + dPhi_(std::numeric_limits::max()), + hit_(nullptr), + et_(0),eta_(0),phi_(0),charge_(0),nrClus_(0){} + + //does not set charge,et,nrclus + SCHitMatch(const GlobalPoint& vtxPos, + const TrajectoryStateOnSurface& trajState, + const TrackingRecHit& hit + ); + ~SCHitMatch()=default; + + void setExtra(float et, float eta, float phi, int charge, int nrClus){ + et_=et; + eta_=eta; + phi_=phi; + charge_=charge; + nrClus_=nrClus; + } + + int subdetId()const{return detId_.subdetId();} + DetId detId()const{return detId_;} + float dRZ()const{return dRZ_;} + float dPhi()const{return dPhi_;} + const GlobalPoint& hitPos()const{return hitPos_;} + float et()const{return et_;} + float eta()const{return eta_;} + float phi()const{return phi_;} + int charge()const{return charge_;} + int nrClus()const{return nrClus_;} + const TrackingRecHit* hit()const{return hit_;} + private: + DetId detId_; + GlobalPoint hitPos_; + float dRZ_; + float dPhi_; + const TrackingRecHit* hit_; //we do not own this + //extra quanities which are set later + float et_; + float eta_; + float phi_; + int charge_; + int nrClus_; + }; + + + + struct MatchInfo { + public: + DetId detId; + float dRZPos,dRZNeg; + float dPhiPos,dPhiNeg; + + MatchInfo(const DetId& iDetId, + float iDRZPos, float iDRZNeg, + float iDPhiPos, float iDPhiNeg): + detId(iDetId),dRZPos(iDRZPos),dRZNeg(iDRZNeg), + dPhiPos(iDPhiPos),dPhiNeg(iDPhiNeg){} + }; + + class SeedWithInfo { + public: + SeedWithInfo(const TrajectorySeed& seed, + const std::vector& posCharge, + const std::vector& negCharge, + int nrValidLayers); + ~SeedWithInfo()=default; + + const TrajectorySeed& seed()const{return seed_;} + float dRZPos(size_t hitNr)const{return getVal(hitNr,&MatchInfo::dRZPos);} + float dRZNeg(size_t hitNr)const{return getVal(hitNr,&MatchInfo::dRZNeg);} + float dPhiPos(size_t hitNr)const{return getVal(hitNr,&MatchInfo::dPhiPos);} + float dPhiNeg(size_t hitNr)const{return getVal(hitNr,&MatchInfo::dPhiNeg);} + DetId detId(size_t hitNr)const{return hitNr& matches()const{return matchInfo_;} + int nrValidLayers()const{return nrValidLayers_;} + private: + float getVal(size_t hitNr,float MatchInfo::*val)const{ + return hitNr::max(); + } + + private: + const TrajectorySeed& seed_; + std::vector matchInfo_; + int nrValidLayers_; + }; + + class MatchingCuts { + public: + MatchingCuts(){} + virtual ~MatchingCuts(){} + virtual bool operator()(const SCHitMatch& scHitMatch)const=0; + }; + + class MatchingCutsV1 : public MatchingCuts { + public: + explicit MatchingCutsV1(const edm::ParameterSet& pset); + bool operator()(const SCHitMatch& scHitMatch)const; + private: + float getDRZCutValue(const float scEt, const float scEta)const; + private: + const double dPhiMax_; + const double dRZMax_; + const double dRZMaxLowEtThres_; + const std::vector dRZMaxLowEtEtaBins_; + const std::vector dRZMaxLowEt_; + }; + + class MatchingCutsV2 : public MatchingCuts { + public: + explicit MatchingCutsV2(const edm::ParameterSet& pset); + bool operator()(const SCHitMatch& scHitMatch)const; + private: + size_t getBinNr(float eta)const; + float getCutValue(float et, float highEt, float highEtThres, float lowEtGrad)const{ + return highEt + std::min(0.f,et-highEtThres)*lowEtGrad; + } + private: + std::vector dPhiHighEt_,dPhiHighEtThres_,dPhiLowEtGrad_; + std::vector dRZHighEt_,dRZHighEtThres_,dRZLowEtGrad_; + std::vector etaBins_; + }; + + +public: + explicit TrajSeedMatcher(const edm::ParameterSet& pset); + ~TrajSeedMatcher()=default; + + static edm::ParameterSetDescription makePSetDescription(); + + void doEventSetup(const edm::EventSetup& iSetup); + + std::vector + compatibleSeeds(const TrajectorySeedCollection& seeds, const GlobalPoint& candPos, + const GlobalPoint & vprim, const float energy); + + void setMeasTkEvtHandle(edm::Handle handle){measTkEvt_=std::move(handle);} + +private: + + std::vector processSeed(const TrajectorySeed& seed, const GlobalPoint& candPos, + const GlobalPoint & vprim, const float energy, const int charge ); + + static float getZVtxFromExtrapolation(const GlobalPoint& primeVtxPos, const GlobalPoint& hitPos, + const GlobalPoint& candPos); + + bool passTrajPreSel(const GlobalPoint& hitPos,const GlobalPoint& candPos)const; + + TrajSeedMatcher::SCHitMatch matchFirstHit(const TrajectorySeed& seed, + const TrajectoryStateOnSurface& trajState, + const GlobalPoint& vtxPos, + const PropagatorWithMaterial& propagator); + + TrajSeedMatcher::SCHitMatch match2ndToNthHit(const TrajectorySeed& seed, + const FreeTrajectoryState& trajState, + const size_t hitNr, + const GlobalPoint& prevHitPos, + const GlobalPoint& vtxPos, + const PropagatorWithMaterial& propagator); + + const TrajectoryStateOnSurface& getTrajStateFromVtx(const TrackingRecHit& hit, + const TrajectoryStateOnSurface& initialState, + const PropagatorWithMaterial& propagator); + + const TrajectoryStateOnSurface& getTrajStateFromPoint(const TrackingRecHit& hit, + const FreeTrajectoryState& initialState, + const GlobalPoint& point, + const PropagatorWithMaterial& propagator); + + void clearCache(); + + bool passesMatchSel(const SCHitMatch& hit, const size_t hitNr)const; + + int getNrValidLayersAlongTraj(const SCHitMatch& hit1, const SCHitMatch& hit2, + const GlobalPoint& candPos, + const GlobalPoint & vprim, + const float energy, const int charge); + + int getNrValidLayersAlongTraj(const DetId& hitId, + const TrajectoryStateOnSurface& hitTrajState)const; + + bool layerHasValidHits(const DetLayer& layer,const TrajectoryStateOnSurface& hitSurState, + const Propagator& propToLayerFromState)const; + + size_t getNrHitsRequired(const int nrValidLayers)const; + +private: + static constexpr float kElectronMass_ = 0.000511; + static constexpr float kPhiCut_ = -0.801144;//cos(2.5) + std::unique_ptr forwardPropagator_; + std::unique_ptr backwardPropagator_; + unsigned long long cacheIDMagField_; + edm::ESHandle magField_; + edm::Handle measTkEvt_; + edm::ESHandle navSchool_; + edm::ESHandle detLayerGeom_; + std::string navSchoolLabel_; + std::string detLayerGeomLabel_; + + bool useRecoVertex_; + std::vector > matchingCuts_; + + //these two varibles determine how hits we require + //based on how many valid layers we had + //right now we always need atleast two hits + //also highly dependent on the seeds you pass in + //which also require a given number of hits + const std::vector minNrHits_; + const std::vector minNrHitsValidLayerBins_; + + std::unordered_map trajStateFromVtxPosChargeCache_; + std::unordered_map trajStateFromVtxNegChargeCache_; + + std::unordered_map,TrajectoryStateOnSurface> trajStateFromPointPosChargeCache_; + std::unordered_map,TrajectoryStateOnSurface> trajStateFromPointNegChargeCache_; + +}; + +#endif diff --git a/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc b/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc index 81809107c5600..6f8e78412f6dd 100644 --- a/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc +++ b/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc @@ -475,7 +475,8 @@ void ElectronSeedGenerator::seedsFromTrajectorySeeds for ( s = pixelSeeds.begin() ; s != pixelSeeds.end() ; s++ ) { reco::ElectronSeed seed(s->seed()) ; - seed.setCaloCluster(cluster,s->hitsMask(),s->subDet2(),s->subDet1(),hoe1,hoe2) ; + seed.setCaloCluster(cluster); + seed.initTwoHitSeed(s->hitsMask()); addSeed(seed,&*s,positron,out) ; } } diff --git a/RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc b/RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc new file mode 100644 index 0000000000000..506895a96816f --- /dev/null +++ b/RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc @@ -0,0 +1,486 @@ +#include "RecoEgamma/EgammaElectronAlgos/interface/TrajSeedMatcher.h" + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" + +#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/PerpendicularBoundPlaneBuilder.h" +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" + +#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" +#include "TrackingTools/RecoGeometry/interface/RecoGeometryRecord.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" + +#include "RecoTracker/Record/interface/NavigationSchoolRecord.h" + +#include "RecoEgamma/EgammaElectronAlgos/interface/FTSFromVertexToPointFactory.h" +#include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h" + +constexpr float TrajSeedMatcher::kElectronMass_; + +TrajSeedMatcher::TrajSeedMatcher(const edm::ParameterSet& pset): + cacheIDMagField_(0), + minNrHits_(pset.getParameter >("minNrHits")), + minNrHitsValidLayerBins_(pset.getParameter >("minNrHitsValidLayerBins")) +{ + useRecoVertex_ = pset.getParameter("useRecoVertex"); + navSchoolLabel_ = pset.getParameter("navSchool"); + detLayerGeomLabel_ = pset.getParameter("detLayerGeom"); + const auto cutsPSets=pset.getParameter >("matchingCuts"); + for(const auto & cutPSet : cutsPSets){ + int version=cutPSet.getParameter("version"); + switch(version){ + case 1: + matchingCuts_.emplace_back(std::make_unique(cutPSet)); + break; + case 2: + matchingCuts_.emplace_back(std::make_unique(cutPSet)); + break; + default: + throw cms::Exception("InvalidConfig") <<" Error TrajSeedMatcher::TrajSeedMatcher pixel match cuts version "<("useRecoVertex",false); + desc.add("navSchool","SimpleNavigationSchool"); + desc.add("detLayerGeom","hltESPGlobalDetLayerGeometry"); + desc.add >("minNrHitsValidLayerBins",{4}); + desc.add >("minNrHits",{2,3}); + + edm::ParameterSetDescription cutsDesc; + auto cutDescCases = + 1 >> + (edm::ParameterDescription("dPhiMax",0.04,true) and + edm::ParameterDescription("dRZMax",0.09,true) and + edm::ParameterDescription("dRZMaxLowEtThres",20.,true) and + edm::ParameterDescription >("dRZMaxLowEtEtaBins",std::vector{1.,1.5},true) and + edm::ParameterDescription >("dRZMaxLowEt",std::vector{0.09,0.15,0.09},true)) or + 2 >> + (edm::ParameterDescription >("dPhiMaxHighEt",{0.003},true) and + edm::ParameterDescription >("dPhiMaxHighEtThres",{0.0},true) and + edm::ParameterDescription >("dPhiMaxLowEtGrad",{0.0},true) and + edm::ParameterDescription >("dRZMaxHighEt",{0.005},true) and + edm::ParameterDescription >("dRZMaxHighEtThres",{30},true) and + edm::ParameterDescription >("dRZMaxLowEtGrad",{-0.002},true) and + edm::ParameterDescription >("etaBins",{},true)); + cutsDesc.ifValue(edm::ParameterDescription("version",1,true), std::move(cutDescCases)); + + edm::ParameterSet defaults; + defaults.addParameter("dPhiMax",0.04); + defaults.addParameter("dRZMax",0.09); + defaults.addParameter("dRZMaxLowEtThres",0.09); + defaults.addParameter >("dRZMaxLowEtEtaBins",std::vector{1.,1.5}); + defaults.addParameter >("dRZMaxLowEt",std::vector{0.09,0.09,0.09}); + defaults.addParameter("version",1); + desc.addVPSet("matchingCuts",cutsDesc,std::vector{defaults,defaults,defaults}); + return desc; +} + +void TrajSeedMatcher::doEventSetup(const edm::EventSetup& iSetup) +{ + if (cacheIDMagField_!=iSetup.get().cacheIdentifier()) { + iSetup.get().get(magField_); + cacheIDMagField_=iSetup.get().cacheIdentifier(); + forwardPropagator_=std::make_unique(alongMomentum,kElectronMass_,&*(magField_)); + backwardPropagator_=std::make_unique(oppositeToMomentum,kElectronMass_,&*(magField_)); + } + iSetup.get().get(navSchoolLabel_,navSchool_); + iSetup.get().get(detLayerGeomLabel_,detLayerGeom_); +} + + +std::vector +TrajSeedMatcher::compatibleSeeds(const TrajectorySeedCollection& seeds, const GlobalPoint& candPos, + const GlobalPoint & vprim, const float energy) +{ + if(!forwardPropagator_ || !backwardPropagator_ || !magField_.isValid()){ + throw cms::Exception("LogicError") <<__FUNCTION__<<" can not make pixel seeds as event setup has not properly been called"; + } + + clearCache(); + + std::vector matchedSeeds; + for(const auto& seed : seeds) { + std::vector matchedHitsNeg = processSeed(seed,candPos,vprim,energy,-1); + std::vector matchedHitsPos = processSeed(seed,candPos,vprim,energy,+1); + int nrValidLayersPos = 0; + int nrValidLayersNeg = 0; + if(matchedHitsNeg.size()>=2){ + nrValidLayersNeg = getNrValidLayersAlongTraj(matchedHitsNeg[0], + matchedHitsNeg[1], + candPos,vprim,energy,-1); + } + if(matchedHitsPos.size()>=2){ + nrValidLayersPos = getNrValidLayersAlongTraj(matchedHitsPos[0], + matchedHitsPos[1], + candPos,vprim,energy,+1); + } + + int nrValidLayers = std::max(nrValidLayersNeg,nrValidLayersPos); + size_t nrHitsRequired = getNrHitsRequired(nrValidLayers); + //so we require the number of hits to exactly match, this is because we currently do not + //do any duplicate cleaning for the input seeds + //this means is a hit pair with a 3rd hit will appear twice (or however many hits it has) + //so if you did >=nrHitsRequired, you would get the same seed multiple times + //ideally we should fix this and clean our input seed collection so each seed is only in once + //also it should be studied what impact having a 3rd hit has on a GsfTrack + //ie will we get a significantly different result seeding with a hit pair + //and the same the hit pair with a 3rd hit added + //in that case, perhaps it should be >= + if(matchedHitsNeg.size()==nrHitsRequired || + matchedHitsPos.size()==nrHitsRequired){ + matchedSeeds.push_back({seed,matchedHitsPos,matchedHitsNeg,nrValidLayers}); + } + + + } + return matchedSeeds; +} + +//checks if the hits of the seed match within requested selection +//matched hits are required to be consecutive, as soon as hit isnt matched, +//the function returns, it doesnt allow skipping hits +std::vector +TrajSeedMatcher::processSeed(const TrajectorySeed& seed, const GlobalPoint& candPos, + const GlobalPoint & vprim, const float energy, const int charge ) +{ + const float candEta = candPos.eta(); + const float candEt = energy*std::sin(candPos.theta()); + + FreeTrajectoryState trajStateFromVtx = FTSFromVertexToPointFactory::get(*magField_, candPos, vprim, energy, charge); + PerpendicularBoundPlaneBuilder bpb; + TrajectoryStateOnSurface initialTrajState(trajStateFromVtx,*bpb(trajStateFromVtx.position(), + trajStateFromVtx.momentum())); + + std::vector matchedHits; + SCHitMatch firstHit = matchFirstHit(seed,initialTrajState,vprim,*backwardPropagator_); + firstHit.setExtra(candEt,candEta,candPos.phi(),charge,1); + if(passesMatchSel(firstHit,0)){ + matchedHits.push_back(firstHit); + + //now we can figure out the z vertex + double zVertex = useRecoVertex_ ? vprim.z() : getZVtxFromExtrapolation(vprim,firstHit.hitPos(),candPos); + GlobalPoint vertex(vprim.x(),vprim.y(),zVertex); + + FreeTrajectoryState firstHitFreeTraj = FTSFromVertexToPointFactory::get(*magField_, firstHit.hitPos(), + vertex, energy, charge) ; + + GlobalPoint prevHitPos = firstHit.hitPos(); + for(size_t hitNr=1;hitNrgdetIndex(); + auto res = trajStateFromVtxCache.find(key); + if(res!=trajStateFromVtxCache.end()) return res->second; + else{ //doesnt exist, need to make it + //FIXME: check for efficiency + auto val = trajStateFromVtxCache.emplace(key,propagator.propagate(initialState,hit.det()->surface())); + return val.first->second; + } +} + +const TrajectoryStateOnSurface& TrajSeedMatcher::getTrajStateFromPoint(const TrackingRecHit& hit, + const FreeTrajectoryState& initialState, + const GlobalPoint& point, + const PropagatorWithMaterial& propagator) +{ + + auto& trajStateFromPointCache = initialState.charge()==1 ? trajStateFromPointPosChargeCache_ : + trajStateFromPointNegChargeCache_; + + auto key = std::make_pair(hit.det()->gdetIndex(),point); + auto res = trajStateFromPointCache.find(key); + if(res!=trajStateFromPointCache.end()) return res->second; + else{ //doesnt exist, need to make it + //FIXME: check for efficiency + auto val = trajStateFromPointCache.emplace(key,propagator.propagate(initialState,hit.det()->surface())); + return val.first->second; + } +} + +TrajSeedMatcher::SCHitMatch TrajSeedMatcher::matchFirstHit(const TrajectorySeed& seed, + const TrajectoryStateOnSurface& initialState, + const GlobalPoint& vtxPos, + const PropagatorWithMaterial& propagator) +{ + const TrajectorySeed::range& hits = seed.recHits(); + auto hitIt = hits.first; + + if(hitIt->isValid()){ + const TrajectoryStateOnSurface& trajStateFromVtx = getTrajStateFromVtx(*hitIt,initialState,propagator); + if(trajStateFromVtx.isValid()) return SCHitMatch(vtxPos,trajStateFromVtx,*hitIt); + } + return SCHitMatch(); +} + +TrajSeedMatcher::SCHitMatch TrajSeedMatcher::match2ndToNthHit(const TrajectorySeed& seed, + const FreeTrajectoryState& initialState, + const size_t hitNr, + const GlobalPoint& prevHitPos, + const GlobalPoint& vtxPos, + const PropagatorWithMaterial& propagator) +{ + const TrajectorySeed::range& hits = seed.recHits(); + auto hitIt = hits.first+hitNr; + + if(hitIt->isValid()){ + const TrajectoryStateOnSurface& trajState = getTrajStateFromPoint(*hitIt,initialState,prevHitPos,propagator); + + if(trajState.isValid()){ + return SCHitMatch(vtxPos,trajState,*hitIt); + } + } + return SCHitMatch(); + +} + +void TrajSeedMatcher::clearCache() +{ + trajStateFromVtxPosChargeCache_.clear(); + trajStateFromVtxNegChargeCache_.clear(); + trajStateFromPointPosChargeCache_.clear(); + trajStateFromPointNegChargeCache_.clear(); +} + +bool TrajSeedMatcher::passesMatchSel(const TrajSeedMatcher::SCHitMatch& hit, const size_t hitNr)const +{ + if(hitNrgeographicalId(),secondHitTraj); +} + +int TrajSeedMatcher::getNrValidLayersAlongTraj(const DetId& hitId, const TrajectoryStateOnSurface& hitTrajState)const +{ + + const DetLayer* detLayer = detLayerGeom_->idToLayer(hitId); + if(detLayer==nullptr) return 0; + + const FreeTrajectoryState& hitFreeState = *hitTrajState.freeState(); + const std::vector inLayers = navSchool_->compatibleLayers(*detLayer,hitFreeState,oppositeToMomentum); + const std::vector outLayers = navSchool_->compatibleLayers(*detLayer,hitFreeState,alongMomentum); + + int nrValidLayers=1; //because our current hit is also valid and wont be included in the count otherwise + int nrPixInLayers=0; + int nrPixOutLayers=0; + for(auto layer : inLayers){ + if(GeomDetEnumerators::isTrackerPixel(layer->subDetector())){ + nrPixInLayers++; + if(layerHasValidHits(*layer,hitTrajState,*backwardPropagator_)) nrValidLayers++; + } + } + for(auto layer : outLayers){ + if(GeomDetEnumerators::isTrackerPixel(layer->subDetector())){ + nrPixOutLayers++; + if(layerHasValidHits(*layer,hitTrajState,*forwardPropagator_)) nrValidLayers++; + } + } + return nrValidLayers; +} + +bool TrajSeedMatcher::layerHasValidHits(const DetLayer& layer, const TrajectoryStateOnSurface& hitSurState, + const Propagator& propToLayerFromState)const +{ + //FIXME: do not hardcode with werid magic numbers stolen from ancient tracking code + //its taken from https://cmssdt.cern.ch/dxr/CMSSW/source/RecoTracker/TrackProducer/interface/TrackProducerBase.icc#165 + //which inspires this code + Chi2MeasurementEstimator estimator(30.,-3.0,0.5,2.0,0.5,1.e12); // same as defauts.... + + const std::vector& detWithState = layer.compatibleDets(hitSurState,propToLayerFromState,estimator); + if(detWithState.empty()) return false; + else{ + DetId id = detWithState.front().first->geographicalId(); + MeasurementDetWithData measDet = measTkEvt_->idToDet(id); + if(measDet.isActive()) return true; + else return false; + } +} + + +size_t TrajSeedMatcher::getNrHitsRequired(const int nrValidLayers)const +{ + for(size_t binNr=0;binNr& posCharge, + const std::vector& negCharge, + int nrValidLayers): + seed_(seed),nrValidLayers_(nrValidLayers) +{ + size_t nrHitsMax = std::max(posCharge.size(),negCharge.size()); + for(size_t hitNr=0;hitNr::max(); + float dPhiPos = hitNr::max(); + + DetId detIdNeg = hitNr::max(); + float dPhiNeg = hitNr::max(); + + if(detIdPos!=detIdNeg && (detIdPos.rawId()!=0 && detIdNeg.rawId()!=0)){ + cms::Exception("LogicError")<<" error in "<<__FILE__<<", "<<__LINE__<<" hits to be combined have different detIDs, this should not be possible and nothing good will come of it"; + } + DetId detId = detIdPos.rawId()!=0 ? detIdPos : detIdNeg; + matchInfo_.push_back(MatchInfo(detId,dRZPos,dRZNeg,dPhiPos,dPhiNeg)); + } +} + +TrajSeedMatcher::MatchingCutsV1::MatchingCutsV1(const edm::ParameterSet& pset): + dPhiMax_(pset.getParameter("dPhiMax")), + dRZMax_(pset.getParameter("dRZMax")), + dRZMaxLowEtThres_(pset.getParameter("dRZMaxLowEtThres")), + dRZMaxLowEtEtaBins_(pset.getParameter >("dRZMaxLowEtEtaBins")), + dRZMaxLowEt_(pset.getParameter >("dRZMaxLowEt")) +{ + if(dRZMaxLowEtEtaBins_.size()+1!=dRZMaxLowEt_.size()){ + throw cms::Exception("InvalidConfig")<<" dRZMaxLowEtEtaBins should be 1 less than dRZMaxLowEt when its "<=0 && std::abs(scHitMatch.dPhi()) > dPhiMax_) return false; + + const float dRZMax = getDRZCutValue(scHitMatch.et(),scHitMatch.eta()); + if(dRZMax_>=0 && std::abs(scHitMatch.dRZ()) > dRZMax) return false; + + return true; +} + +float TrajSeedMatcher::MatchingCutsV1::getDRZCutValue(const float scEt, const float scEta)const +{ + if(scEt>=dRZMaxLowEtThres_) return dRZMax_; + else{ + const float absEta = std::abs(scEta); + for(size_t etaNr=0;etaNr >("dPhiMaxHighEt")), + dPhiHighEtThres_(pset.getParameter >("dPhiMaxHighEtThres")), + dPhiLowEtGrad_(pset.getParameter >("dPhiMaxLowEtGrad")), + dRZHighEt_(pset.getParameter >("dRZMaxHighEt")), + dRZHighEtThres_(pset.getParameter >("dRZMaxHighEtThres")), + dRZLowEtGrad_(pset.getParameter >("dRZMaxLowEtGrad")), + etaBins_(pset.getParameter >("etaBins")) +{ + auto binSizeCheck = [](size_t sizeEtaBins,const std::vector& vec,const std::string& name){ + if(vec.size()!=sizeEtaBins+1){ + throw cms::Exception("InvalidConfig")<<" when constructing TrajSeedMatcher::MatchingCutsV2 "<< name<<" has "<=0 && std::abs(scHitMatch.dPhi()) > dPhiMax) return false; + float dRZMax = getCutValue(scHitMatch.et(),dRZHighEt_[binNr],dRZHighEtThres_[binNr],dRZLowEtGrad_[binNr]); + if(dRZMax>=0 && std::abs(scHitMatch.dRZ()) > dRZMax) return false; + + return true; +} + +//eta bins is exactly 1 smaller than the vectors which will be accessed by this bin nr +size_t TrajSeedMatcher::MatchingCutsV2::getBinNr(float eta)const +{ + const float absEta = std::abs(eta); + for(size_t etaNr=0;etaNr { +public: + + + explicit ElectronNHitSeedProducer( const edm::ParameterSet & ) ; + virtual ~ElectronNHitSeedProducer()=default; + + virtual void produce( edm::Event &, const edm::EventSetup & ) override final; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + + TrajSeedMatcher matcher_; + + std::vector> > superClustersTokens_; + edm::EDGetTokenT initialSeedsToken_ ; + edm::EDGetTokenT > verticesToken_; + edm::EDGetTokenT beamSpotToken_ ; + edm::EDGetTokenT measTkEvtToken_; + +}; + +namespace { + template + edm::Handle getHandle(const edm::Event& event,const edm::EDGetTokenT& token) + { + edm::Handle handle; + event.getByToken(token,handle); + return handle; + } + + template + GlobalPoint convertToGP(const T& orgPoint){ + return GlobalPoint(orgPoint.x(),orgPoint.y(),orgPoint.z()); + } + + int getLayerOrDiskNr(DetId detId,const TrackerTopology& trackerTopo) + { + if(detId.subdetId()==PixelSubdetector::PixelBarrel){ + return trackerTopo.pxbLayer(detId); + }else if(detId.subdetId()==PixelSubdetector::PixelEndcap){ + return trackerTopo.pxfDisk(detId); + }else return -1; + } + + reco::ElectronSeed::PMVars + makeSeedPixelVar(const TrajSeedMatcher::MatchInfo& matchInfo, + const TrackerTopology& trackerTopo) + { + + int layerOrDisk = getLayerOrDiskNr(matchInfo.detId,trackerTopo); + reco::ElectronSeed::PMVars pmVars; + pmVars.setDet(matchInfo.detId,layerOrDisk); + pmVars.setDPhi(matchInfo.dPhiPos,matchInfo.dPhiNeg); + pmVars.setDRZ(matchInfo.dRZPos,matchInfo.dRZNeg); + + return pmVars; + } + +} + +ElectronNHitSeedProducer::ElectronNHitSeedProducer( const edm::ParameterSet& pset): + matcher_(pset.getParameter("matcherConfig")), + initialSeedsToken_(consumes(pset.getParameter("initialSeeds"))), + verticesToken_(consumes >(pset.getParameter("vertices"))), + beamSpotToken_(consumes(pset.getParameter("beamSpot"))), + measTkEvtToken_(consumes(pset.getParameter("measTkEvt"))) +{ + const auto superClusTags = pset.getParameter >("superClusters"); + for(const auto& scTag : superClusTags){ + superClustersTokens_.emplace_back(consumes>(scTag)); + } + produces() ; +} + +void ElectronNHitSeedProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) +{ + edm::ParameterSetDescription desc; + desc.add("initialSeeds",edm::InputTag("hltElePixelSeedsCombined")); + desc.add("vertices",edm::InputTag()); + desc.add("beamSpot",edm::InputTag("hltOnlineBeamSpot")); + desc.add("measTkEvt",edm::InputTag("hltSiStripClusters")); + desc.add >("superClusters",std::vector{edm::InputTag{"hltEgammaSuperClustersToPixelMatch"}}); + desc.add("matcherConfig",TrajSeedMatcher::makePSetDescription()); + + descriptions.add("electronNHitSeedProducer",desc); +} + +void ElectronNHitSeedProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + edm::ESHandle trackerTopoHandle; + iSetup.get().get(trackerTopoHandle); + + matcher_.doEventSetup(iSetup); + matcher_.setMeasTkEvtHandle(getHandle(iEvent,measTkEvtToken_)); + + auto eleSeeds = std::make_unique(); + + auto initialSeedsHandle = getHandle(iEvent,initialSeedsToken_); + + auto beamSpotHandle = getHandle(iEvent,beamSpotToken_); + GlobalPoint primVtxPos = convertToGP(beamSpotHandle->position()); + + for(const auto& superClustersToken : superClustersTokens_){ + auto superClustersHandle = getHandle(iEvent,superClustersToken); + for(auto& superClusRef : *superClustersHandle){ + + //the eta of the supercluster when mustache clustered is slightly biased due to bending in magnetic field + //the eta of its seed cluster is a better estimate of the orginal position + GlobalPoint caloPosition(GlobalPoint::Polar(superClusRef->seed()->position().theta(), //seed theta + superClusRef->position().phi(), //supercluster phi + superClusRef->position().r())); //supercluster r + + + const std::vector matchedSeeds = + matcher_.compatibleSeeds(*initialSeedsHandle,caloPosition, + primVtxPos,superClusRef->energy()); + + for(auto& matchedSeed : matchedSeeds){ + reco::ElectronSeed eleSeed(matchedSeed.seed()); + reco::ElectronSeed::CaloClusterRef caloClusRef(superClusRef); + eleSeed.setCaloCluster(caloClusRef); + eleSeed.setNrLayersAlongTraj(matchedSeed.nrValidLayers()); + for(auto& matchInfo : matchedSeed.matches()){ + eleSeed.addHitInfo(makeSeedPixelVar(matchInfo,*trackerTopoHandle)); + } + eleSeeds->emplace_back(eleSeed); + } + } + + } + iEvent.put(std::move(eleSeeds)); +} + + + +DEFINE_FWK_MODULE(ElectronNHitSeedProducer); diff --git a/RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersEDProducer.cc b/RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersEDProducer.cc new file mode 100644 index 0000000000000..3c1e0f9743b54 --- /dev/null +++ b/RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersEDProducer.cc @@ -0,0 +1,9 @@ +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "RecoTracker/TkTrackingRegions/interface/TrackingRegionProducerFactory.h" +#include "RecoTracker/TkTrackingRegions/interface/TrackingRegionEDProducerT.h" +#include "RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersProducer.h" +DEFINE_EDM_PLUGIN(TrackingRegionProducerFactory, TrackingRegionsFromSuperClustersProducer, "TrackingRegionsFromSuperClustersProducer"); +using TrackingRegionsFromSuperClustersEDProducer = TrackingRegionEDProducerT; +DEFINE_FWK_MODULE(TrackingRegionsFromSuperClustersEDProducer); diff --git a/RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersProducer.h b/RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersProducer.h new file mode 100644 index 0000000000000..a0426aae87871 --- /dev/null +++ b/RecoEgamma/EgammaElectronProducers/plugins/TrackingRegionsFromSuperClustersProducer.h @@ -0,0 +1,314 @@ +#ifndef RecoTracker_TkTrackingRegions_TrackingRegionsFromSuperClustersProducer_H +#define RecoTracker_TkTrackingRegions_TrackingRegionsFromSuperClustersProducer_H + +//****************************************************************************** +// +// Part of the refactorisation of of the E/gamma pixel matching pre-2017 +// This refactorisation converts the monolithic approach to a series of +// independent producer modules, with each modules performing a specific +// job as recommended by the 2017 tracker framework +// +// This module is called a Producer even though its not an ED producer +// This was done to be consistant with other TrackingRegion producers +// in RecoTracker/TkTrackingRegions +// +// The module closely follows the other TrackingRegion producers +// in RecoTracker/TkTrackingRegions and is intended to become an EDProducer +// by TrackingRegionEDProducerT + +// This module c tracking regions from the superclusters. It mostly +// replicates the functionality of the SeedFilter class +// although unlike that class, it does not actually create seeds +// +// Author : Sam Harper (RAL), 2017 +// +//******************************************************************************* + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" + +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/GeometryVector/interface/GlobalVector.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" + + +#include "RecoTracker/TkTrackingRegions/interface/TrackingRegionProducer.h" +#include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h" +#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" + +#include "RecoEgamma/EgammaElectronAlgos/interface/FTSFromVertexToPointFactory.h" + +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "CommonTools/Utils/interface/StringToEnumValue.h" + +//stick this in common tools +#include "TEnum.h" +#include "TEnumConstant.h" +namespace{ + template MyEnum strToEnum(std::string const& enumConstName){ + TEnum* en = TEnum::GetEnum(typeid(MyEnum)); + if (en != nullptr){ + if (TEnumConstant const* enc = en->GetConstant(enumConstName.c_str())){ + return static_cast(enc->GetValue()); + }else{ + throw cms::Exception("Configuration") < RectangularEtaPhiTrackingRegion::UseMeasurementTracker strToEnum(std::string const& enumConstName){ + using MyEnum = RectangularEtaPhiTrackingRegion::UseMeasurementTracker; + if(enumConstName=="kNever") return MyEnum::kNever; + else if(enumConstName=="kForSiStrips") return MyEnum::kForSiStrips; + else if(enumConstName=="kAlways") return MyEnum::kAlways; + else{ + throw cms::Exception("InvalidConfiguration") < > + regions (const edm::Event& iEvent, const edm::EventSetup& iSetup)const override; + +private: + GlobalPoint getVtxPos(const edm::Event& iEvent,double& deltaZVertex)const; + + std::unique_ptr + createTrackingRegion(const reco::SuperCluster& superCluster,const GlobalPoint& vtxPos, + const double deltaZVertex,const Charge charge, + const MeasurementTrackerEvent* measTrackerEvent, + const MagneticField& magField)const; + +private: + void validateConfigSettings()const; + +private: + //there are 3 modes in which to define the Z area of the tracking region + //1) from first vertex in the passed vertices collection +/- originHalfLength in z (useZInVertex=true) + //2) the beamspot +/- nrSigmaForBSDeltaZ* zSigma of beamspot (useZInBeamspot=true) + // the zSigma of the beamspot can have a minimum value specified + // we do this because a common error is that beamspot has too small of a value + //3) defaultZ_ +/- originHalfLength (if useZInVertex and useZInBeamspot are both false) + double ptMin_; + double originRadius_; + double originHalfLength_; + double deltaEtaRegion_; + double deltaPhiRegion_; + bool useZInVertex_; + bool useZInBeamspot_; + double nrSigmaForBSDeltaZ_; + double defaultZ_; + double minBSDeltaZ_; + bool precise_; + RectangularEtaPhiTrackingRegion::UseMeasurementTracker whereToUseMeasTracker_; + + edm::EDGetTokenT verticesToken_; + edm::EDGetTokenT beamSpotToken_; + edm::EDGetTokenT measTrackerEventToken_; + std::vector> > superClustersTokens_; + +}; + + + +namespace { + template + edm::Handle getHandle(const edm::Event& event,const edm::EDGetTokenT & token){ + edm::Handle handle; + event.getByToken(token,handle); + return handle; + } +} + +TrackingRegionsFromSuperClustersProducer:: +TrackingRegionsFromSuperClustersProducer(const edm::ParameterSet& cfg, + edm::ConsumesCollector && iC) +{ + edm::ParameterSet regionPSet = cfg.getParameter("RegionPSet"); + + ptMin_ = regionPSet.getParameter("ptMin"); + originRadius_ = regionPSet.getParameter("originRadius"); + originHalfLength_ = regionPSet.getParameter("originHalfLength"); + deltaPhiRegion_ = regionPSet.getParameter("deltaPhiRegion"); + deltaEtaRegion_ = regionPSet.getParameter("deltaEtaRegion"); + useZInVertex_ = regionPSet.getParameter("useZInVertex"); + useZInBeamspot_ = regionPSet.getParameter("useZInBeamspot"); + nrSigmaForBSDeltaZ_ = regionPSet.getParameter("nrSigmaForBSDeltaZ"); + defaultZ_ = regionPSet.getParameter("defaultZ"); + minBSDeltaZ_ = regionPSet.getParameter("minBSDeltaZ"); + precise_ = regionPSet.getParameter("precise"); + whereToUseMeasTracker_ = strToEnum(regionPSet.getParameter("whereToUseMeasTracker")); + + validateConfigSettings(); + + auto verticesTag = regionPSet.getParameter("vertices"); + auto beamSpotTag = regionPSet.getParameter("beamSpot"); + auto superClustersTags = regionPSet.getParameter >("superClusters"); + auto measTrackerEventTag = regionPSet.getParameter("measurementTrackerEvent"); + + if(useZInVertex_){ + verticesToken_ = iC.consumes(verticesTag); + }else{ + beamSpotToken_ = iC.consumes(beamSpotTag); + } + if(whereToUseMeasTracker_ != RectangularEtaPhiTrackingRegion::UseMeasurementTracker::kNever){ + measTrackerEventToken_ = iC.consumes(measTrackerEventTag); + } + for(const auto& tag : superClustersTags){ + superClustersTokens_.emplace_back(iC.consumes>(tag)); + } +} + + + +void TrackingRegionsFromSuperClustersProducer:: +fillDescriptions(edm::ConfigurationDescriptions& descriptions) +{ + edm::ParameterSetDescription desc; + + desc.add("ptMin", 1.5); + desc.add("originRadius", 0.2); + desc.add("originHalfLength", 15.0)->setComment("z range is +/- this value except when using the beamspot (useZInBeamspot=true)"); + desc.add("deltaPhiRegion",0.4); + desc.add("deltaEtaRegion",0.1); + desc.add("useZInVertex", false)->setComment("use the leading vertex position +/-orginHalfLength, mutually exclusive with useZInBeamspot"); + desc.add("useZInBeamspot", true)->setComment("use the beamspot position +/- nrSigmaForBSDeltaZ* sigmaZ_{bs}, mutually exclusive with useZInVertex"); + desc.add("nrSigmaForBSDeltaZ",3.0)->setComment("# of sigma to extend the z region when using the beamspot, only active if useZInBeamspot=true"); + desc.add("minBSDeltaZ",0.0)->setComment("a minimum value of the beamspot sigma z to use, only active if useZInBeamspot=true"); + desc.add("defaultZ",0.)->setComment("the default z position, only used if useZInVertex and useZInBeamspot are both false"); + desc.add("precise", true); + desc.add("whereToUseMeasTracker","kNever"); + desc.add("beamSpot", edm::InputTag("hltOnlineBeamSpot"))->setComment("only used if useZInBeamspot is true"); + desc.add("vertices", edm::InputTag())->setComment("only used if useZInVertex is true"); + desc.add >("superClusters", std::vector{edm::InputTag{"hltEgammaSuperClustersToPixelMatch"}}); + desc.add("measurementTrackerEvent",edm::InputTag()); + + edm::ParameterSetDescription descRegion; + descRegion.add("RegionPSet", desc); + + descriptions.add("trackingRegionsFromSuperClusters", descRegion); +} + + + +std::vector > +TrackingRegionsFromSuperClustersProducer:: +regions(const edm::Event& iEvent, const edm::EventSetup& iSetup)const +{ + std::vector > trackingRegions; + + double deltaZVertex=0; + GlobalPoint vtxPos = getVtxPos(iEvent,deltaZVertex); + + const MeasurementTrackerEvent *measTrackerEvent = nullptr; + if(!measTrackerEventToken_.isUninitialized()){ + measTrackerEvent = getHandle(iEvent,measTrackerEventToken_).product(); + } + edm::ESHandle magFieldHandle; + iSetup.get().get(magFieldHandle); + + for(auto& superClustersToken : superClustersTokens_){ + auto superClustersHandle = getHandle(iEvent,superClustersToken); + for(auto& superClusterRef : *superClustersHandle){ + //do both charge hypothesises + trackingRegions.emplace_back(createTrackingRegion(*superClusterRef,vtxPos,deltaZVertex,Charge::POS,measTrackerEvent,*magFieldHandle)); + trackingRegions.emplace_back(createTrackingRegion(*superClusterRef,vtxPos,deltaZVertex,Charge::NEG,measTrackerEvent,*magFieldHandle)); + } + } + return trackingRegions; +} + +GlobalPoint TrackingRegionsFromSuperClustersProducer:: +getVtxPos(const edm::Event& iEvent,double& deltaZVertex)const +{ + if(useZInVertex_){ + auto verticesHandle = getHandle(iEvent,verticesToken_); + //we throw if the vertices are not there but if no vertex is + //recoed in the event, we default to 0,0,defaultZ as the vertex + if(!verticesHandle->empty()){ + deltaZVertex = originHalfLength_; + const auto& pv = verticesHandle->front(); + return GlobalPoint(pv.x(),pv.y(),pv.z()); + }else{ + deltaZVertex = originHalfLength_; + return GlobalPoint(0,0,defaultZ_); + } + }else{ + + auto beamSpotHandle = getHandle(iEvent,beamSpotToken_); + const reco::BeamSpot::Point& bsPos = beamSpotHandle->position(); + + if(useZInBeamspot_){ + //as this is what has been done traditionally for e/gamma, others just use sigmaZ + const double bsSigmaZ = std::sqrt(beamSpotHandle->sigmaZ()*beamSpotHandle->sigmaZ() + + beamSpotHandle->sigmaZ0Error()*beamSpotHandle->sigmaZ0Error()); + const double sigmaZ = std::max(bsSigmaZ,minBSDeltaZ_); + deltaZVertex = nrSigmaForBSDeltaZ_*sigmaZ; + + return GlobalPoint(bsPos.x(),bsPos.y(),bsPos.z()); + }else{ + deltaZVertex = originHalfLength_; + return GlobalPoint(bsPos.x(),bsPos.y(),defaultZ_); + } + } + +} + +std::unique_ptr +TrackingRegionsFromSuperClustersProducer:: +createTrackingRegion(const reco::SuperCluster& superCluster,const GlobalPoint& vtxPos, + const double deltaZVertex,const Charge charge, + const MeasurementTrackerEvent* measTrackerEvent, + const MagneticField& magField)const +{ + const GlobalPoint clusterPos(superCluster.position().x(), superCluster.position().y(), superCluster.position().z()); + const double energy = superCluster.energy(); + + FreeTrajectoryState freeTrajState = FTSFromVertexToPointFactory::get(magField, clusterPos, vtxPos, energy, static_cast(charge)); + return std::make_unique(freeTrajState.momentum(), + vtxPos, + ptMin_, + originRadius_, + deltaZVertex, + deltaEtaRegion_, + deltaPhiRegion_, + whereToUseMeasTracker_, + precise_, + measTrackerEvent); +} + +void TrackingRegionsFromSuperClustersProducer::validateConfigSettings()const +{ + if(useZInVertex_ && useZInBeamspot_){ + throw cms::Exception("InvalidConfiguration") <<" when constructing TrackingRegionsFromSuperClustersProducer both useZInVertex ("< edm::Handle getHandle(const edm::Event& event,const edm::EDGetTokenT& token) + { + edm::Handle handle; + event.getByToken(token,handle); + return handle; + } + +} + +class EgammaHLTFilteredSuperClusterProducer : public edm::stream::EDProducer<>{ +public: + class SelectionCut { + private: + struct CutValues { + float cut; + float cutOverE; + float cutOverE2; + bool useEt; + std::function compFunc; + + CutValues(const edm::ParameterSet& pset): + cut(pset.getParameter("cut")), + cutOverE(pset.getParameter("cutOverE")), + cutOverE2(pset.getParameter("cutOverE2")), + useEt(pset.getParameter("useEt")), + compFunc(std::less()) {} + + bool operator()(const reco::RecoEcalCandidate& cand,float value)const{ + float energy = useEt ? cand.et() : cand.energy(); + return compFunc(value,cut) || compFunc(value/energy,cutOverE) || + compFunc(value/energy/energy,cutOverE2); + } + }; + public: + SelectionCut(const edm::ParameterSet& pset,edm::ConsumesCollector && iC): + ebCut_(pset.getParameter("barrelCut")), + eeCut_(pset.getParameter("endcapCut")), + varToken_(iC.consumes(pset.getParameter("var"))) + {} + + ~SelectionCut()=default; + + bool operator()(const reco::RecoEcalCandidateRef& cand)const{ + CutValues cut = std::abs(cand->eta())<1.479 ? ebCut_ : eeCut_; + return cut(*cand,getVar(cand)); + } + + float getVar(const reco::RecoEcalCandidateRef& cand)const{ + auto res = varHandle_->find(cand); + if(res!=varHandle_->end()) return res->val; + else{ + //FIX ME: add some provenance info to this + throw cms::Exception("LogicError") <<" candidate not found in collection "; + } + } + + void getHandles(const edm::Event& event){ + event.getByToken(varToken_,varHandle_); + } + private: + CutValues ebCut_; + CutValues eeCut_; + edm::EDGetTokenT varToken_; + edm::Handle varHandle_; + }; + + explicit EgammaHLTFilteredSuperClusterProducer(const edm::ParameterSet& pset); + ~EgammaHLTFilteredSuperClusterProducer()=default; + static void fillDescriptions(edm::ConfigurationDescriptions & descriptions); + + void produce(edm::Event&, const edm::EventSetup&) override; + + +private: + edm::EDGetTokenT candsToken_; + std::vector cuts_; +}; + +EgammaHLTFilteredSuperClusterProducer:: +EgammaHLTFilteredSuperClusterProducer(const edm::ParameterSet& pset): + candsToken_(consumes(pset.getParameter("cands"))) +{ + const auto& cutPsets = pset.getParameter >("cuts"); + for(auto& cutPset : cutPsets){ + cuts_.push_back(SelectionCut(cutPset,consumesCollector())); + } + + produces>(); + +} + +void EgammaHLTFilteredSuperClusterProducer:: +fillDescriptions(edm::ConfigurationDescriptions & descriptions) +{ + edm::ParameterSetDescription desc; + desc.add("cands",edm::InputTag("hltEgammaCandidates")); + + edm::ParameterSetDescription cutsDesc; + edm::ParameterSetDescription regionCutsDesc; + regionCutsDesc.add("cut",-1); + regionCutsDesc.add("cutOverE",-1); + regionCutsDesc.add("cutOverE2",-1); + regionCutsDesc.add("useEt",false); + edm::ParameterSet cutDefaults; + cutDefaults.addParameter("cutOverE",0.2); + cutDefaults.addParameter("useEt",false); + + cutsDesc.add("barrelCut",regionCutsDesc); + cutsDesc.add("endcapCut",regionCutsDesc); + cutsDesc.add("var",edm::InputTag("hltEgammaHoverE")); + + edm::ParameterSet defaults; + defaults.addParameter("var",edm::InputTag("hltEgammaHoverE")); + defaults.addParameter("barrelCut",cutDefaults); + defaults.addParameter("endcapCut",cutDefaults); + desc.addVPSet("cuts",cutsDesc,std::vector{defaults}); + + descriptions.add("egammaHLTFilteredSuperClusterProducer",desc); + +} + +void EgammaHLTFilteredSuperClusterProducer:: +produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + for(auto& cut : cuts_) cut.getHandles(iEvent); + auto candsHandle = getHandle(iEvent,candsToken_); + + auto outputSCs = std::make_unique>(); + + for(size_t candNr=0;candNrsize();candNr++){ + reco::RecoEcalCandidateRef candRef(candsHandle,candNr); + bool passAllCuts=true; + for(const auto& cut: cuts_){ + if(!cut(candRef)){ + passAllCuts=false; + break; + } + } + if(passAllCuts) outputSCs->push_back(candRef->superCluster()); + } + + iEvent.put(std::move(outputSCs)); +} + +DEFINE_FWK_MODULE(EgammaHLTFilteredSuperClusterProducer); + +#endif diff --git a/RecoEgamma/EgammaHLTProducers/src/EgammaHLTPixelMatchVarProducer.cc b/RecoEgamma/EgammaHLTProducers/src/EgammaHLTPixelMatchVarProducer.cc index 2a90e39891b9a..a2aa595f9d581 100644 --- a/RecoEgamma/EgammaHLTProducers/src/EgammaHLTPixelMatchVarProducer.cc +++ b/RecoEgamma/EgammaHLTProducers/src/EgammaHLTPixelMatchVarProducer.cc @@ -1,7 +1,7 @@ #include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -24,14 +24,88 @@ #include "RecoEgamma/EgammaHLTProducers/interface/EgammaHLTPixelMatchParamObjects.h" -class EgammaHLTPixelMatchVarProducer : public edm::global::EDProducer<> { +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "FWCore/Framework/interface/ESHandle.h" + +namespace { + //first 4 bits are sub detect of each hit (0=barrel, 1 = endcap) + //next 8 bits are layer information (0=no hit, 1 = hit), first 4 are barrel, next 4 are endcap (theres an empty bit here + //next 4 bits are nr of layers info + int makeSeedInfo(const reco::ElectronSeed& seed){ + int info = 0; + for(size_t hitNr=0;hitNr& candHandle): + name_(std::move(name)),hitNr_(hitNr),func_(func), + val_(std::numeric_limits::max()), + valInfo_(0) + { + valMap_=std::make_unique(candHandle); + valInfoMap_=std::make_unique(candHandle); + } + PixelData(PixelData&& rhs)=default; + + void resetVal(){val_=std::numeric_limits::max();valInfo_=0;} + void fill(const reco::ElectronSeed& seed){ + if(hitNr_insert(candRef,val_); + valInfoMap_->insert(candRef,valInfo_); + val_ = std::numeric_limits::max(); + valInfo_ = 0; + } + + void putInto(edm::Event& event){ + event.put(std::move(valMap_),name_+std::to_string(hitNr_+1)); + event.put(std::move(valInfoMap_),name_+std::to_string(hitNr_+1)+"Info"); + } + +private: + std::unique_ptr valMap_; + std::unique_ptr valInfoMap_; + std::string name_; + size_t hitNr_; + float (reco::ElectronSeed::*func_)(size_t)const; + float val_; + float valInfo_; + +}; + + +class EgammaHLTPixelMatchVarProducer : public edm::stream::EDProducer<> { public: explicit EgammaHLTPixelMatchVarProducer(const edm::ParameterSet&); ~EgammaHLTPixelMatchVarProducer(); static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - void produce(edm::StreamID sid, edm::Event&, const edm::EventSetup&) const override; + void produce(edm::Event&, const edm::EventSetup&) override; std::array calS2(const reco::ElectronSeed& seed,int charge)const; private: @@ -45,7 +119,7 @@ class EgammaHLTPixelMatchVarProducer : public edm::global::EDProducer<> { egPM::Param dRZ2Para_; int productsToWrite_; - + size_t nrHits_; }; EgammaHLTPixelMatchVarProducer::EgammaHLTPixelMatchVarProducer(const edm::ParameterSet& config) : @@ -54,7 +128,8 @@ EgammaHLTPixelMatchVarProducer::EgammaHLTPixelMatchVarProducer(const edm::Parame dPhi1Para_(config.getParameter("dPhi1SParams")), dPhi2Para_(config.getParameter("dPhi2SParams")), dRZ2Para_(config.getParameter("dRZ2SParams")), - productsToWrite_(config.getParameter("productsToWrite")) + productsToWrite_(config.getParameter("productsToWrite")), + nrHits_(4) { //register your products @@ -64,7 +139,19 @@ EgammaHLTPixelMatchVarProducer::EgammaHLTPixelMatchVarProducer(const edm::Parame produces < reco::RecoEcalCandidateIsolationMap >("dPhi2BestS2"); produces < reco::RecoEcalCandidateIsolationMap >("dzBestS2"); } - + if(productsToWrite_>=2){ + //note for product names we start from index 1 + for(size_t hitNr=1;hitNr<=nrHits_;hitNr++){ + produces < reco::RecoEcalCandidateIsolationMap >("dPhi"+std::to_string(hitNr)); + produces < reco::RecoEcalCandidateIsolationMap >("dPhi"+std::to_string(hitNr)+"Info"); + produces < reco::RecoEcalCandidateIsolationMap >("dRZ"+std::to_string(hitNr)); + produces < reco::RecoEcalCandidateIsolationMap >("dRZ"+std::to_string(hitNr)+"Info"); + } + produces < reco::RecoEcalCandidateIsolationMap >("nrClus"); + produces < reco::RecoEcalCandidateIsolationMap >("seedClusEFrac"); + produces < reco::RecoEcalCandidateIsolationMap >("phiWidth"); + produces < reco::RecoEcalCandidateIsolationMap >("etaWidth"); + } } @@ -114,7 +201,7 @@ void EgammaHLTPixelMatchVarProducer::fillDescriptions(edm::ConfigurationDescript descriptions.add(("hltEgammaHLTPixelMatchVarProducer"), desc); } -void EgammaHLTPixelMatchVarProducer::produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const { +void EgammaHLTPixelMatchVarProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup){ // Get the HLT filtered objects edm::Handle recoEcalCandHandle; @@ -124,15 +211,38 @@ void EgammaHLTPixelMatchVarProducer::produce(edm::StreamID sid, edm::Event& iEve edm::Handle pixelSeedsHandle; iEvent.getByToken(pixelSeedsToken_,pixelSeedsHandle); - if(!recoEcalCandHandle.isValid() || !pixelSeedsHandle.isValid()) return; + if(!recoEcalCandHandle.isValid()) return; + else if(!pixelSeedsHandle.isValid()){ + auto s2Map = std::make_unique(recoEcalCandHandle); + for(unsigned int candNr = 0; candNrsize(); candNr++) { + reco::RecoEcalCandidateRef candRef(recoEcalCandHandle,candNr); + s2Map->insert(candRef,0); + } + iEvent.put(std::move(s2Map),"s2"); + return; + } + edm::ESHandle trackerTopoHandle; + iSetup.get().get(trackerTopoHandle); + auto dPhi1BestS2Map = std::make_unique(recoEcalCandHandle); auto dPhi2BestS2Map = std::make_unique(recoEcalCandHandle); auto dzBestS2Map = std::make_unique(recoEcalCandHandle); auto s2Map = std::make_unique(recoEcalCandHandle); + auto nrClusMap = std::make_unique(recoEcalCandHandle); + auto seedClusEFracMap = std::make_unique(recoEcalCandHandle); + auto phiWidthMap = std::make_unique(recoEcalCandHandle); + auto etaWidthMap = std::make_unique(recoEcalCandHandle); + + std::vector pixelData; + for(size_t hitNr=0;hitNrsize(); candNr++) { - + reco::RecoEcalCandidateRef candRef(recoEcalCandHandle,candNr); reco::SuperClusterRef candSCRef = candRef->superCluster(); @@ -147,6 +257,11 @@ void EgammaHLTPixelMatchVarProducer::produce(edm::StreamID sid, edm::Event& iEve if(s2Data[0]=2){ + for(auto& pixelDatum : pixelData){ + pixelDatum.fill(seed); + } + } } } @@ -157,7 +272,21 @@ void EgammaHLTPixelMatchVarProducer::produce(edm::StreamID sid, edm::Event& iEve dPhi2BestS2Map->insert(candRef,bestS2[2]); dzBestS2Map->insert(candRef,bestS2[3]); } - + if(productsToWrite_>=2){ + nrClusMap->insert(candRef,candSCRef->clustersSize()); + float seedClusEFrac = candSCRef->rawEnergy()>0 ? candSCRef->seed()->energy() / candSCRef->rawEnergy() : 0.; + // std::cout <<"cand "<energy()<<" E Corr "<correctedEnergyUncertainty()<<" "<correctedEnergy()<<" width "<phiWidth()<phiWidth(); + seedClusEFracMap->insert(candRef,seedClusEFrac); + float phiWidth = candSCRef->phiWidth(); + float etaWidth = candSCRef->etaWidth(); + phiWidthMap->insert(candRef,phiWidth); + etaWidthMap->insert(candRef,etaWidth); + + for(auto& pixelDatum : pixelData){ + pixelDatum.fill(candRef); + } + } } iEvent.put(std::move(s2Map),"s2"); @@ -166,6 +295,15 @@ void EgammaHLTPixelMatchVarProducer::produce(edm::StreamID sid, edm::Event& iEve iEvent.put(std::move(dPhi2BestS2Map),"dPhi2BestS2"); iEvent.put(std::move(dzBestS2Map),"dzBestS2"); } + if(productsToWrite_>=2){ + for(auto& pixelDatum : pixelData){ + pixelDatum.putInto(iEvent); + } + iEvent.put(std::move(nrClusMap),"nrClus"); + iEvent.put(std::move(seedClusEFracMap),"seedClusEFrac"); + iEvent.put(std::move(phiWidthMap),"phiWidth"); + iEvent.put(std::move(etaWidthMap),"etaWidth"); + } } std::array EgammaHLTPixelMatchVarProducer::calS2(const reco::ElectronSeed& seed,int charge)const @@ -174,9 +312,9 @@ std::array EgammaHLTPixelMatchVarProducer::calS2(const reco::ElectronSe const float dPhi2Const = dPhi2Para_(seed); const float dRZ2Const = dRZ2Para_(seed); - float dPhi1 = (charge <0 ? seed.dPhi1() : seed.dPhi1Pos())/dPhi1Const; - float dPhi2 = (charge <0 ? seed.dPhi2() : seed.dPhi2Pos())/dPhi2Const; - float dRz2 = (charge <0 ? seed.dRz2() : seed.dRz2Pos())/dRZ2Const; + float dPhi1 = (charge <0 ? seed.dPhiNeg(0) : seed.dPhiPos(0))/dPhi1Const; + float dPhi2 = (charge <0 ? seed.dPhiNeg(1) : seed.dPhiPos(1))/dPhi2Const; + float dRz2 = (charge <0 ? seed.dRZNeg(1) : seed.dRZPos(1))/dRZ2Const; float s2 = dPhi1*dPhi1+dPhi2*dPhi2+dRz2*dRz2; return std::array{{s2,dPhi1,dPhi2,dRz2}}; From f99b4fc44eeb0d64fdc8efebe536833f9f974804 Mon Sep 17 00:00:00 2001 From: Carlo Battilana Date: Sun, 1 Oct 2017 21:54:04 +0200 Subject: [PATCH 5/6] muon associators using phase2 tracker --- SimMuon/MCTruth/python/MuonAssociatorByHits_cfi.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/SimMuon/MCTruth/python/MuonAssociatorByHits_cfi.py b/SimMuon/MCTruth/python/MuonAssociatorByHits_cfi.py index 2990c6e086317..6869cae8b632e 100644 --- a/SimMuon/MCTruth/python/MuonAssociatorByHits_cfi.py +++ b/SimMuon/MCTruth/python/MuonAssociatorByHits_cfi.py @@ -70,8 +70,10 @@ # associatePixel = cms.bool(True), associateStrip = cms.bool(True), + usePhase2Tracker = cms.bool(False), pixelSimLinkSrc = cms.InputTag("simSiPixelDigis"), stripSimLinkSrc = cms.InputTag("simSiStripDigis"), + phase2TrackerSimLinkSrc = cms.InputTag("simSiPixelDigis","Tracker"), associateRecoTracks = cms.bool(True), # ROUList = cms.vstring('TrackerHitsTIBLowTof', @@ -142,6 +144,5 @@ from Configuration.Eras.Modifier_run3_GEM_cff import run3_GEM run3_GEM.toModify( muonAssociatorByHits, useGEMs = cms.bool(True) ) from Configuration.Eras.Modifier_phase2_tracker_cff import phase2_tracker +phase2_tracker.toModify( muonAssociatorByHits, usePhase2Tracker = cms.bool(True) ) phase2_tracker.toModify( muonAssociatorByHits, pixelSimLinkSrc = "simSiPixelDigis:Pixel" ) -phase2_tracker.toModify( muonAssociatorByHits, stripSimLinkSrc = "simSiPixelDigis:Tracker" ) - From 6e272dbda5804fa90db6f3eaffb7cdc8569a55bf Mon Sep 17 00:00:00 2001 From: Carlo Battilana Date: Mon, 2 Oct 2017 22:02:05 +0200 Subject: [PATCH 6/6] Backport fix from PR 18801 --- .../TkDetLayers/src/Phase2OTtiltedBarrelLayer.cc | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/RecoTracker/TkDetLayers/src/Phase2OTtiltedBarrelLayer.cc b/RecoTracker/TkDetLayers/src/Phase2OTtiltedBarrelLayer.cc index 1cbde6eceff75..cc288fc5f0921 100644 --- a/RecoTracker/TkDetLayers/src/Phase2OTtiltedBarrelLayer.cc +++ b/RecoTracker/TkDetLayers/src/Phase2OTtiltedBarrelLayer.cc @@ -97,11 +97,14 @@ Phase2OTtiltedBarrelLayer::groupedCompatibleDetsV( const TrajectoryStateOnSurfac vector closestResultNeg; vector closestResultPos; Phase2OTBarrelLayer::groupedCompatibleDetsV(tsos, prop, est, closestResultRods); - for(auto ring : theNegativeRingsComps){ - ring->groupedCompatibleDetsV(tsos, prop, est, closestResultNeg); - } - for(auto ring : thePositiveRingsComps){ - ring->groupedCompatibleDetsV(tsos, prop, est, closestResultPos); + if(tsos.globalPosition().z()<0){ + for(auto& ring : theNegativeRingsComps){ + ring->groupedCompatibleDetsV(tsos, prop, est, closestResultNeg); + } + } else { + for(auto& ring : thePositiveRingsComps){ + ring->groupedCompatibleDetsV(tsos, prop, est, closestResultPos); + } } result.assign(closestResultRods.begin(),closestResultRods.end());