Merge pull request #22005 from JanFSchulte/candidateSeededRecovery10X

Candidate Seeded option for automated doublet recovery in HLT tracking

RecoTracker/TkTrackingRegions/plugins/AreaSeededTrackingRegionsBuilder.cc

@@ -4,6 +4,7 @@
 #include "DataFormats/Common/interface/Handle.h"
 #include "DataFormats/Math/interface/PtEtaPhiMass.h"
 #include "DataFormats/Math/interface/deltaPhi.h"
+#include "DataFormats/Math/interface/normalizedPhi.h"
 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
 
 #include <array>
@@ -15,7 +16,9 @@ namespace {
   }
 }
 
-AreaSeededTrackingRegionsBuilder::AreaSeededTrackingRegionsBuilder(const edm::ParameterSet& regPSet, edm::ConsumesCollector& iC) {
+AreaSeededTrackingRegionsBuilder::AreaSeededTrackingRegionsBuilder(const edm::ParameterSet& regPSet, edm::ConsumesCollector& iC):
+ candidates_(regPSet, iC)
+ {
   m_extraPhi = regPSet.getParameter<double>("extraPhi");
   m_extraEta = regPSet.getParameter<double>("extraEta");
 
@@ -40,7 +43,7 @@ void AreaSeededTrackingRegionsBuilder::fillDescriptions(edm::ParameterSetDescrip
 
   desc.add<std::string>("whereToUseMeasurementTracker", "Never");
   desc.add<edm::InputTag>("measurementTrackerName", edm::InputTag(""));
-
+  TrackingSeedCandidates::fillDescriptions(desc);
   desc.add<bool>("searchOpt", false);
 }
 
@@ -52,7 +55,7 @@ AreaSeededTrackingRegionsBuilder::Builder AreaSeededTrackingRegionsBuilder::begi
     e.getByToken(token_measurementTracker, hmte);
     builder.setMeasurementTracker(hmte.product());
   }
-
+  builder.setCandidates( ( candidates_.objects(e)));
   return builder;
 }
 
@@ -63,7 +66,9 @@ std::vector<std::unique_ptr<TrackingRegion> > AreaSeededTrackingRegionsBuilder::
   // create tracking regions in directions of the points of interest
   int n_regions = 0;
   for(const auto& origin: origins) {
-    result.push_back(region(origin, areas));
+    auto reg = region(origin, areas);
+    if (!reg) continue;
+    result.push_back(std::move(reg));
     ++n_regions;
   }
   LogDebug("AreaSeededTrackingRegionsBuilder") << "produced "<<n_regions<<" regions";
@@ -228,29 +233,93 @@ std::unique_ptr<TrackingRegion> AreaSeededTrackingRegionsBuilder::Builder::regio
 
   const auto meanEta = (minEta+maxEta)/2.f;
   const auto meanPhi = (minPhi+maxPhi)/2.f;
-  const auto deltaEta = maxEta-meanEta + m_conf->m_extraEta;
-  const auto deltaPhi = maxPhi-meanPhi + m_conf->m_extraPhi;
+  const auto dEta = maxEta-meanEta + m_conf->m_extraEta;
+  const auto dPhi = maxPhi-meanPhi + m_conf->m_extraPhi;
+
+  int n_objects = 0;
+  if (candidates.first) n_objects = candidates.first->size();
+
+  if (n_objects > 0){
+	// If we have objected used for seeding, loop over objects and find overlap with the found region. Return overlaps as tracking regions to use
+	for(const auto& object : *candidates.first) {
+		float dEta_Cand = candidates.second.first;
+		float dPhi_Cand = candidates.second.second;
+	       	float eta_Cand = object.eta();
+		float phi_Cand = object.phi();
+          	float dEta_Cand_Point = std::abs(eta_Cand-meanEta);
+          	float dPhi_Cand_Point = std::abs(deltaPhi(phi_Cand,meanPhi));
+
+		if(dEta_Cand_Point > (dEta_Cand + dEta) || dPhi_Cand_Point > (dPhi_Cand + dPhi)) continue;
+
+		float  etaMin_RoI = std::max(eta_Cand-dEta_Cand,meanEta-dEta);
+		float  etaMax_RoI = std::min(eta_Cand+dEta_Cand,meanEta+dEta);
+
+		float  phi_Cand_minus  = normalizedPhi(phi_Cand-dPhi_Cand);
+		float  phi_Point_minus = normalizedPhi(meanPhi-dPhi);
+		float  phi_Cand_plus  = normalizedPhi(phi_Cand+dPhi_Cand);
+		float  phi_Point_plus = normalizedPhi(meanPhi+dPhi);
+
+		float phiMin_RoI = deltaPhi(phi_Cand_minus,phi_Point_minus)>0. ? phi_Cand_minus : phi_Point_minus ;	
+		float phiMax_RoI = deltaPhi(phi_Cand_plus,phi_Point_plus)<0. ? phi_Cand_plus : phi_Point_plus;
+
+
+		const auto meanEtaTemp = (etaMin_RoI+etaMax_RoI)/2.f;
+		auto meanPhiTemp = (phiMin_RoI+phiMax_RoI)/2.f;
+		if( phiMax_RoI < phiMin_RoI ) meanPhiTemp-=M_PI;
+	        meanPhiTemp = normalizedPhi(meanPhiTemp);
+
+		const auto dPhiTemp = deltaPhi(phiMax_RoI,meanPhiTemp);
+		const auto dEtaTemp = etaMax_RoI-meanEtaTemp;
+
+		const auto x = std::cos(meanPhiTemp);
+  	  	const auto y = std::sin(meanPhiTemp);
+  	  	const auto z = 1./std::tan(2.f*std::atan(std::exp(-meanEtaTemp)));
+
+  		LogTrace("AreaSeededTrackingRegionsBuilder") << "Direction x,y,z " << x << "," << y << "," << z
+                                               << " eta,phi " << meanEtaTemp << "," << meanPhiTemp
+                                               << " window eta " << (meanEtaTemp-dEtaTemp) << "," << (meanEtaTemp+dEtaTemp)
+                                               << " phi " << (meanPhiTemp-dPhiTemp) << "," << (meanPhiTemp+dPhiTemp);
+
+  		return std::make_unique<RectangularEtaPhiTrackingRegion>(
+      			GlobalVector(x,y,z),
+      			origin.first, // GlobalPoint
+      			m_conf->m_ptMin,
+      			m_conf->m_originRadius,
+      			origin.second,
+      			dEtaTemp,
+      			dPhiTemp,
+      			m_conf->m_whereToUseMeasurementTracker,
+      			m_conf->m_precise,
+      			m_measurementTracker,
+      			m_conf->m_searchOpt
+  		);  
+	}
+	// Have to retun nullptr here to ensure that we always return something
+	return nullptr;
 
-  const auto x = std::cos(meanPhi);
-  const auto y = std::sin(meanPhi);
-  const auto z = (x*x+y*y)/std::tan(2.f*std::atan(std::exp(-meanEta))); // simplify?
+  }
+  else{
+  	const auto x = std::cos(meanPhi);
+  	const auto y = std::sin(meanPhi);
+  	const auto z = 1./std::tan(2.f*std::atan(std::exp(-meanEta))); 
 
-  LogTrace("AreaSeededTrackingRegionsBuilder") << "Direction x,y,z " << x << "," << y << "," << z
+  	LogTrace("AreaSeededTrackingRegionsBuilder") << "Direction x,y,z " << x << "," << y << "," << z
                                                << " eta,phi " << meanEta << "," << meanPhi
-                                               << " window eta " << (meanEta-deltaEta) << "," << (meanEta+deltaEta)
-                                               << " phi " << (meanPhi-deltaPhi) << "," << (meanPhi+deltaPhi);
-
-  return std::make_unique<RectangularEtaPhiTrackingRegion>(
-      GlobalVector(x,y,z),
-      origin.first, // GlobalPoint
-      m_conf->m_ptMin,
-      m_conf->m_originRadius,
-      origin.second,
-      deltaEta,
-      deltaPhi,
-      m_conf->m_whereToUseMeasurementTracker,
-      m_conf->m_precise,
-      m_measurementTracker,
-      m_conf->m_searchOpt
-  );  
+                                               << " window eta " << (meanEta-dEta) << "," << (meanEta+dEta)
+                                               << " phi " << (meanPhi-dPhi) << "," << (meanPhi+dPhi);
+
+  	return std::make_unique<RectangularEtaPhiTrackingRegion>(
+      		GlobalVector(x,y,z),
+      		origin.first, // GlobalPoint
+      		m_conf->m_ptMin,
+      		m_conf->m_originRadius,
+      		origin.second,
+      		dEta,
+      		dPhi,
+      		m_conf->m_whereToUseMeasurementTracker,
+      		m_conf->m_precise,
+      		m_measurementTracker,
+      		m_conf->m_searchOpt
+  	);  
+  }
 }

RecoTracker/TkTrackingRegions/plugins/AreaSeededTrackingRegionsBuilder.h

@@ -12,6 +12,7 @@
 #include "DataFormats/VertexReco/interface/Vertex.h"
 #include "DataFormats/VertexReco/interface/VertexFwd.h"
 
+#include "TrackingSeedCandidates.h"
 class AreaSeededTrackingRegionsBuilder {
 public:
   using Origin = std::pair<GlobalPoint, float>; // (origin, half-length in z)
@@ -78,6 +79,7 @@ class AreaSeededTrackingRegionsBuilder {
     ~Builder() = default;
 
     void setMeasurementTracker(const MeasurementTrackerEvent *mte) { m_measurementTracker = mte; }
+    void setCandidates(const TrackingSeedCandidates::Objects cands) { candidates = cands; }
 
     std::vector<std::unique_ptr<TrackingRegion> > regions(const Origins& origins, const std::vector<Area>& areas) const;
     std::unique_ptr<TrackingRegion> region(const Origin& origin, const std::vector<Area>& areas) const;
@@ -89,6 +91,7 @@ class AreaSeededTrackingRegionsBuilder {
 
     const AreaSeededTrackingRegionsBuilder *m_conf = nullptr;
     const MeasurementTrackerEvent *m_measurementTracker = nullptr;
+    TrackingSeedCandidates::Objects candidates;
   };
 
   AreaSeededTrackingRegionsBuilder(const edm::ParameterSet& regPSet, edm::ConsumesCollector&& iC): AreaSeededTrackingRegionsBuilder(regPSet, iC) {}
@@ -101,7 +104,7 @@ class AreaSeededTrackingRegionsBuilder {
 
 private:
   std::vector<Area> m_areas;
-
+  TrackingSeedCandidates candidates_;
   float m_extraPhi;
   float m_extraEta;
   float m_ptMin;

RecoTracker/TkTrackingRegions/plugins/PixelInactiveAreaTrackingRegionsSeedingLayersProducer.cc

@@ -75,6 +75,7 @@ void PixelInactiveAreaTrackingRegionsSeedingLayersProducer::produce(edm::Event&
     LogTrace("PixelInactiveAreaTrackingRegionsSeedingLayersProducer") << "Origin " << origin.first.x() << "," << origin.first.y() << "," << origin.first.z() << " z half lengh " << origin.second;
     for(auto& areasLayerSet: areasLayerSets) {
       auto region = builder.region(origin, areasLayerSet.first);
+      if(!region) continue;
 #ifdef EDM_ML_DEBUG
       auto etaPhiRegion = dynamic_cast<const RectangularEtaPhiTrackingRegion *>(region.get());
       std::stringstream ss;

RecoTracker/TkTrackingRegions/plugins/TrackingSeedCandidates.cc

@@ -0,0 +1,40 @@
+#include "TrackingSeedCandidates.h"
+
+TrackingSeedCandidates::TrackingSeedCandidates(const edm::ParameterSet& regPSet, edm::ConsumesCollector& iC) {
+  // operation mode
+  //
+    std::string seedingModeString = regPSet.getParameter<std::string>("seedingMode");
+    if      (seedingModeString == "Candidate")      m_seedingMode = SeedingMode::CANDIDATE_SEEDED;
+    else if (seedingModeString == "Global")          m_seedingMode = SeedingMode::GLOBAL;
+    else throw edm::Exception(edm::errors::Configuration) << "Unknown seeding mode string: "<<seedingModeString;
+
+    m_deltaEta_Cand = regPSet.getParameter<double>("deltaEta_Cand"); 
+    m_deltaPhi_Cand = regPSet.getParameter<double>("deltaPhi_Cand");
+
+    // basic inputs
+    if(m_seedingMode == SeedingMode::CANDIDATE_SEEDED){
+      m_token_input        = iC.consumes<reco::CandidateView>(regPSet.getParameter<edm::InputTag>("input"));
+      if (m_deltaEta_Cand<0 || m_deltaPhi_Cand<0) throw edm::Exception(edm::errors::Configuration) << "Delta eta and phi parameters must be set for candidates in candidate seeding mode"; 
+    }
+}
+
+void TrackingSeedCandidates::fillDescriptions(edm::ParameterSetDescription& desc) {
+  desc.add<std::string>("seedingMode", "Global");
+  desc.add<edm::InputTag>("input",edm::InputTag());
+  desc.add<double>("deltaEta_Cand", -1.);
+  desc.add<double>("deltaPhi_Cand", -1.);
+}
+
+TrackingSeedCandidates::Objects TrackingSeedCandidates::objects(const edm::Event& iEvent) const {
+
+ Objects result;
+ std::pair <float,float> dimensions = std::make_pair(m_deltaEta_Cand,m_deltaPhi_Cand);
+ edm::Handle< reco::CandidateView > objects;
+
+ if(m_seedingMode == SeedingMode::CANDIDATE_SEEDED){
+      iEvent.getByToken( m_token_input, objects ); 
+      result = std::make_pair(objects.product(),dimensions); 
+  }
+  else result = std::make_pair(nullptr,dimensions); 
+  return result;
+}

RecoTracker/TkTrackingRegions/plugins/TrackingSeedCandidates.h

@@ -0,0 +1,38 @@
+#ifndef RecoTracker_TkTrackingRegions_TrackingSeedCandidates_h
+#define RecoTracker_TkTrackingRegions_TrackingSeedCandidates_h
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "DataFormats/Candidate/interface/Candidate.h"
+#include "DataFormats/Common/interface/Handle.h"
+
+#include <vector>
+#include <utility>
+
+class TrackingSeedCandidates {
+public:
+
+  enum class SeedingMode {CANDIDATE_SEEDED, GLOBAL};
+  //using Objects = std::pair< edm::Handle< reco::CandidateView > , std::pair < float, float > > ; // (origin, half-length in z)
+  using Objects = std::pair< const reco::CandidateView* , std::pair < float, float > > ; // (origin, half-length in z)
+  TrackingSeedCandidates(const edm::ParameterSet& regPSet, edm::ConsumesCollector&& iC): TrackingSeedCandidates(regPSet, iC) {}
+  TrackingSeedCandidates(const edm::ParameterSet& regPSet, edm::ConsumesCollector& iC);
+  ~TrackingSeedCandidates() = default;
+
+  static void fillDescriptions(edm::ParameterSetDescription& desc);
+
+  Objects objects(const edm::Event& iEvent) const;
+
+private:
+  SeedingMode m_seedingMode;
+  float m_deltaEta_Cand;
+  float m_deltaPhi_Cand;
+
+  edm::EDGetTokenT<reco::CandidateView> m_token_input;
+
+};
+
+
+#endif