Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-01 22:41:04

0001 #ifndef CandidatePointSeededTrackingRegionsProducer_h
0002 #define CandidatePointSeededTrackingRegionsProducer_h
0003 
0004 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionProducer.h"
0005 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0006 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0007 
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Framework/interface/ConsumesCollector.h"
0013 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0014 #include "DataFormats/Common/interface/Handle.h"
0015 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0016 #include "DataFormats/VertexReco/interface/Vertex.h"
0017 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0018 #include "DataFormats/Candidate/interface/Candidate.h"
0019 #include "DataFormats/Math/interface/deltaPhi.h"
0020 #include "DataFormats/Math/interface/normalizedPhi.h"
0021 
0022 #include "VertexBeamspotOrigins.h"
0023 
0024 /** class CandidatePointSeededTrackingRegionsProducer
0025  *
0026  * eta-phi TrackingRegions producer in directions defined by Candidate-based objects of interest
0027  * from a collection defined by the "input" parameter.
0028  *
0029  * Four operational modes are supported ("operationMode" parameter):
0030  *
0031  *   BeamSpotFixed:
0032  *     origin is defined by the beam spot
0033  *     z-half-length is defined by a fixed zErrorBeamSpot parameter
0034  *   BeamSpotSigma:
0035  *     origin is defined by the beam spot
0036  *     z-half-length is defined by nSigmaZBeamSpot * beamSpot.sigmaZ
0037  *   VerticesFixed:
0038  *     origins are defined by vertices from VertexCollection (use maximum MaxNVertices of them)
0039  *     z-half-length is defined by a fixed zErrorVetex parameter
0040  *   VerticesSigma:
0041  *     origins are defined by vertices from VertexCollection (use maximum MaxNVertices of them)
0042  *     z-half-length is defined by nSigmaZVertex * vetex.zError
0043  *
0044  *   If, while using one of the "Vertices" modes, there's no vertices in an event, we fall back into
0045  *   either BeamSpotSigma or BeamSpotFixed mode, depending on the positiveness of nSigmaZBeamSpot.
0046  *
0047  * Three seeding modes are supported ("seedingMode" parameter):
0048  *
0049  *   Candidate-seeded:
0050  *      defines regions around candidates from the "input" collection
0051  *   Point-seeded:
0052  *      defines regions around fixed points in the detector (previously in PointSeededTrackingRegionsProducer)
0053  *   Candidate+point-seeded:
0054  *      defines regions as intersections of regions around candidates from the "input" collections and around fixed points in the detector
0055  *
0056  *   \authors Vadim Khotilovich + Thomas Strebler
0057  */
0058 class CandidatePointSeededTrackingRegionsProducer : public TrackingRegionProducer {
0059 public:
0060   enum class SeedingMode { CANDIDATE_SEEDED, POINT_SEEDED, CANDIDATE_POINT_SEEDED };
0061 
0062   explicit CandidatePointSeededTrackingRegionsProducer(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0063       : m_origins(conf.getParameter<edm::ParameterSet>("RegionPSet"), iC), token_field(iC.esConsumes()) {
0064     edm::ParameterSet regPSet = conf.getParameter<edm::ParameterSet>("RegionPSet");
0065 
0066     // seeding mode
0067     std::string seedingModeString = regPSet.getParameter<std::string>("seedingMode");
0068     if (seedingModeString == "Candidate")
0069       m_seedingMode = SeedingMode::CANDIDATE_SEEDED;
0070     else if (seedingModeString == "Point")
0071       m_seedingMode = SeedingMode::POINT_SEEDED;
0072     else if (seedingModeString == "CandidatePoint")
0073       m_seedingMode = SeedingMode::CANDIDATE_POINT_SEEDED;
0074     else
0075       throw edm::Exception(edm::errors::Configuration) << "Unknown seeding mode string: " << seedingModeString;
0076 
0077     // basic inputs
0078     if (m_seedingMode == SeedingMode::CANDIDATE_SEEDED || m_seedingMode == SeedingMode::CANDIDATE_POINT_SEEDED)
0079       m_token_input = iC.consumes<reco::CandidateView>(regPSet.getParameter<edm::InputTag>("input"));
0080 
0081     // Specific points in the detector
0082     if (m_seedingMode == SeedingMode::POINT_SEEDED || m_seedingMode == SeedingMode::CANDIDATE_POINT_SEEDED) {
0083       edm::ParameterSet points = regPSet.getParameter<edm::ParameterSet>("points");
0084       std::vector<double> etaPoints = points.getParameter<std::vector<double>>("eta");
0085       std::vector<double> phiPoints = points.getParameter<std::vector<double>>("phi");
0086 
0087       if (!(etaPoints.size() == phiPoints.size()))
0088         throw edm::Exception(edm::errors::Configuration) << "The parameters 'eta' and 'phi' must have the same size";
0089       if (etaPoints.empty())
0090         throw edm::Exception(edm::errors::Configuration)
0091             << "At least one point should be defined for point or candidate+point seeding modes";
0092 
0093       for (size_t i = 0; i < etaPoints.size(); ++i) {
0094         m_etaPhiPoints.push_back(std::make_pair(etaPoints[i], phiPoints[i]));
0095 
0096         double x = std::cos(phiPoints[i]);
0097         double y = std::sin(phiPoints[i]);
0098         double theta = 2 * std::atan(std::exp(-etaPoints[i]));
0099         double z = 1. / std::tan(theta);
0100         GlobalVector direction(x, y, z);
0101         m_directionPoints.push_back(direction);
0102       }
0103     }
0104 
0105     m_maxNRegions = regPSet.getParameter<unsigned int>("maxNRegions");
0106     if (m_maxNRegions == 0)
0107       throw edm::Exception(edm::errors::Configuration) << "maxNRegions should be greater than or equal to 1";
0108 
0109     // RectangularEtaPhiTrackingRegion parameters:
0110     m_ptMin = regPSet.getParameter<double>("ptMin");
0111     m_originRadius = regPSet.getParameter<double>("originRadius");
0112 
0113     if (m_seedingMode == SeedingMode::CANDIDATE_SEEDED) {
0114       m_deltaEta_Cand = regPSet.getParameter<double>("deltaEta_Cand");
0115       m_deltaPhi_Cand = regPSet.getParameter<double>("deltaPhi_Cand");
0116       if (m_deltaEta_Cand < 0 || m_deltaPhi_Cand < 0)
0117         throw edm::Exception(edm::errors::Configuration)
0118             << "Delta eta and phi parameters must be set for candidates in candidate seeding mode";
0119     } else if (m_seedingMode == SeedingMode::POINT_SEEDED) {
0120       m_deltaEta_Point = regPSet.getParameter<double>("deltaEta_Point");
0121       m_deltaPhi_Point = regPSet.getParameter<double>("deltaPhi_Point");
0122       if (m_deltaEta_Point < 0 || m_deltaPhi_Point < 0)
0123         throw edm::Exception(edm::errors::Configuration)
0124             << "Delta eta and phi parameters must be set for points in point seeding mode";
0125     } else if (m_seedingMode == SeedingMode::CANDIDATE_POINT_SEEDED) {
0126       m_deltaEta_Cand = regPSet.getParameter<double>("deltaEta_Cand");
0127       m_deltaPhi_Cand = regPSet.getParameter<double>("deltaPhi_Cand");
0128       m_deltaEta_Point = regPSet.getParameter<double>("deltaEta_Point");
0129       m_deltaPhi_Point = regPSet.getParameter<double>("deltaPhi_Point");
0130       if (m_deltaEta_Cand < 0 || m_deltaPhi_Cand < 0 || m_deltaEta_Point < 0 || m_deltaPhi_Point < 0)
0131         throw edm::Exception(edm::errors::Configuration) << "Delta eta and phi parameters must be set separately for "
0132                                                             "candidates and points in candidate+point seeding mode";
0133     }
0134 
0135     m_precise = regPSet.getParameter<bool>("precise");
0136     m_whereToUseMeasurementTracker = RectangularEtaPhiTrackingRegion::stringToUseMeasurementTracker(
0137         regPSet.getParameter<std::string>("whereToUseMeasurementTracker"));
0138     if (m_whereToUseMeasurementTracker != RectangularEtaPhiTrackingRegion::UseMeasurementTracker::kNever) {
0139       m_token_measurementTracker =
0140           iC.consumes<MeasurementTrackerEvent>(regPSet.getParameter<edm::InputTag>("measurementTrackerName"));
0141     }
0142     m_searchOpt = false;
0143     if (regPSet.exists("searchOpt"))
0144       m_searchOpt = regPSet.getParameter<bool>("searchOpt");
0145 
0146     if (m_precise) {
0147       token_msmaker = iC.esConsumes();
0148     }
0149   }
0150 
0151   ~CandidatePointSeededTrackingRegionsProducer() override = default;
0152 
0153   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0154     edm::ParameterSetDescription desc;
0155 
0156     desc.add<std::string>("seedingMode", "Candidate");
0157 
0158     desc.add<edm::InputTag>("input", edm::InputTag(""));
0159     edm::ParameterSetDescription descPoints;
0160     descPoints.add<std::vector<double>>("eta", {});
0161     descPoints.add<std::vector<double>>("phi", {});
0162     desc.add<edm::ParameterSetDescription>("points", descPoints);
0163 
0164     desc.add<unsigned int>("maxNRegions", 10);
0165 
0166     VertexBeamspotOrigins::fillDescriptions(desc, "hltOnlineBeamSpot", "hltPixelVertices", 1);
0167 
0168     desc.add<double>("ptMin", 0.9);
0169     desc.add<double>("originRadius", 0.2);
0170     desc.add<double>("deltaEta_Cand", -1.);
0171     desc.add<double>("deltaPhi_Cand", -1.);
0172     desc.add<double>("deltaEta_Point", -1.);
0173     desc.add<double>("deltaPhi_Point", -1.);
0174     desc.add<bool>("precise", true);
0175 
0176     desc.add<std::string>("whereToUseMeasurementTracker", "ForSiStrips");
0177     desc.add<edm::InputTag>("measurementTrackerName", edm::InputTag(""));
0178 
0179     desc.add<bool>("searchOpt", false);
0180 
0181     // Only for backwards-compatibility
0182     edm::ParameterSetDescription descRegion;
0183     descRegion.add<edm::ParameterSetDescription>("RegionPSet", desc);
0184 
0185     descriptions.add("candidatePointSeededTrackingRegionsFromBeamSpot", descRegion);
0186   }
0187 
0188   std::vector<std::unique_ptr<TrackingRegion>> regions(const edm::Event& e, const edm::EventSetup& es) const override {
0189     std::vector<std::unique_ptr<TrackingRegion>> result;
0190 
0191     // pick up the candidate objects of interest
0192     edm::Handle<reco::CandidateView> objects;
0193     size_t n_objects = 0;
0194 
0195     if (m_seedingMode == SeedingMode::CANDIDATE_SEEDED || m_seedingMode == SeedingMode::CANDIDATE_POINT_SEEDED) {
0196       e.getByToken(m_token_input, objects);
0197       n_objects = objects->size();
0198       if (n_objects == 0)
0199         return result;
0200     }
0201 
0202     const auto& objs = *objects;
0203 
0204     const auto& origins = m_origins.origins(e);
0205     if (origins.empty()) {
0206       return result;
0207     }
0208 
0209     const MeasurementTrackerEvent* measurementTracker = nullptr;
0210     if (!m_token_measurementTracker.isUninitialized()) {
0211       edm::Handle<MeasurementTrackerEvent> hmte;
0212       e.getByToken(m_token_measurementTracker, hmte);
0213       measurementTracker = hmte.product();
0214     }
0215 
0216     const auto& field = es.getData(token_field);
0217     const MultipleScatteringParametrisationMaker* msmaker = nullptr;
0218     if (m_precise) {
0219       msmaker = &es.getData(token_msmaker);
0220     }
0221 
0222     // create tracking regions (maximum MaxNRegions of them) in directions of the
0223     // objects of interest (we expect that the collection was sorted in decreasing pt order)
0224     int n_regions = 0;
0225 
0226     if (m_seedingMode == SeedingMode::CANDIDATE_SEEDED) {
0227       for (const auto& object : objs) {
0228         GlobalVector direction(object.momentum().x(), object.momentum().y(), object.momentum().z());
0229 
0230         for (const auto& origin : origins) {
0231           result.push_back(std::make_unique<RectangularEtaPhiTrackingRegion>(direction,
0232                                                                              origin.first,
0233                                                                              m_ptMin,
0234                                                                              m_originRadius,
0235                                                                              origin.second,
0236                                                                              m_deltaEta_Cand,
0237                                                                              m_deltaPhi_Cand,
0238                                                                              field,
0239                                                                              msmaker,
0240                                                                              m_precise,
0241                                                                              m_whereToUseMeasurementTracker,
0242                                                                              measurementTracker,
0243                                                                              m_searchOpt));
0244           ++n_regions;
0245           if (n_regions >= m_maxNRegions)
0246             break;
0247         }
0248 
0249         if (n_regions >= m_maxNRegions)
0250           break;
0251       }
0252 
0253     }
0254 
0255     else if (m_seedingMode == SeedingMode::POINT_SEEDED) {
0256       for (const auto& direction : m_directionPoints) {
0257         for (const auto& origin : origins) {
0258           result.push_back(std::make_unique<RectangularEtaPhiTrackingRegion>(direction,     // GlobalVector
0259                                                                              origin.first,  // GlobalPoint
0260                                                                              m_ptMin,
0261                                                                              m_originRadius,
0262                                                                              origin.second,
0263                                                                              m_deltaEta_Point,
0264                                                                              m_deltaPhi_Point,
0265                                                                              field,
0266                                                                              msmaker,
0267                                                                              m_precise,
0268                                                                              m_whereToUseMeasurementTracker,
0269                                                                              measurementTracker,
0270                                                                              m_searchOpt));
0271           ++n_regions;
0272           if (n_regions >= m_maxNRegions)
0273             break;
0274         }
0275 
0276         if (n_regions >= m_maxNRegions)
0277           break;
0278       }
0279 
0280     }
0281 
0282     else if (m_seedingMode == SeedingMode::CANDIDATE_POINT_SEEDED) {
0283       for (const auto& object : objs) {
0284         double eta_Cand = object.eta();
0285         double phi_Cand = object.phi();
0286 
0287         for (const auto& etaPhiPoint : m_etaPhiPoints) {
0288           double eta_Point = etaPhiPoint.first;
0289           double phi_Point = etaPhiPoint.second;
0290           double dEta_Cand_Point = std::abs(eta_Cand - eta_Point);
0291           double dPhi_Cand_Point = std::abs(deltaPhi(phi_Cand, phi_Point));
0292 
0293           //Check if there is an overlap between Candidate- and Point-based regions of interest
0294           if (dEta_Cand_Point > (m_deltaEta_Cand + m_deltaEta_Point) ||
0295               dPhi_Cand_Point > (m_deltaPhi_Cand + m_deltaPhi_Point))
0296             continue;
0297 
0298           //Determines boundaries of intersection of RoIs
0299           double etaMin_RoI = std::max(eta_Cand - m_deltaEta_Cand, eta_Point - m_deltaEta_Point);
0300           double etaMax_RoI = std::min(eta_Cand + m_deltaEta_Cand, eta_Point + m_deltaEta_Point);
0301 
0302           double phi_Cand_minus = normalizedPhi(phi_Cand - m_deltaPhi_Cand);
0303           double phi_Point_minus = normalizedPhi(phi_Point - m_deltaPhi_Point);
0304           double phi_Cand_plus = normalizedPhi(phi_Cand + m_deltaPhi_Cand);
0305           double phi_Point_plus = normalizedPhi(phi_Point + m_deltaPhi_Point);
0306 
0307           double phiMin_RoI = deltaPhi(phi_Cand_minus, phi_Point_minus) > 0. ? phi_Cand_minus : phi_Point_minus;
0308           double phiMax_RoI = deltaPhi(phi_Cand_plus, phi_Point_plus) < 0. ? phi_Cand_plus : phi_Point_plus;
0309 
0310           //Determines position and width of new RoI
0311           double eta_RoI = 0.5 * (etaMax_RoI + etaMin_RoI);
0312           double deltaEta_RoI = etaMax_RoI - eta_RoI;
0313 
0314           double phi_RoI = 0.5 * (phiMax_RoI + phiMin_RoI);
0315           if (phiMax_RoI < phiMin_RoI)
0316             phi_RoI -= M_PI;
0317           phi_RoI = normalizedPhi(phi_RoI);
0318           double deltaPhi_RoI = deltaPhi(phiMax_RoI, phi_RoI);
0319 
0320           double x = std::cos(phi_RoI);
0321           double y = std::sin(phi_RoI);
0322           double theta = 2 * std::atan(std::exp(-eta_RoI));
0323           double z = 1. / std::tan(theta);
0324 
0325           GlobalVector direction(x, y, z);
0326 
0327           for (const auto& origin : origins) {
0328             result.push_back(std::make_unique<RectangularEtaPhiTrackingRegion>(direction,
0329                                                                                origin.first,
0330                                                                                m_ptMin,
0331                                                                                m_originRadius,
0332                                                                                origin.second,
0333                                                                                deltaEta_RoI,
0334                                                                                deltaPhi_RoI,
0335                                                                                field,
0336                                                                                msmaker,
0337                                                                                m_precise,
0338                                                                                m_whereToUseMeasurementTracker,
0339                                                                                measurementTracker,
0340                                                                                m_searchOpt));
0341             ++n_regions;
0342             if (n_regions >= m_maxNRegions)
0343               break;
0344           }
0345 
0346           if (n_regions >= m_maxNRegions)
0347             break;
0348         }
0349 
0350         if (n_regions >= m_maxNRegions)
0351           break;
0352       }
0353     }
0354 
0355     edm::LogInfo("CandidatePointSeededTrackingRegionsProducer") << "produced " << n_regions << " regions";
0356 
0357     return result;
0358   }
0359 
0360 private:
0361   VertexBeamspotOrigins m_origins;
0362   SeedingMode m_seedingMode;
0363 
0364   int m_maxNRegions;
0365   edm::EDGetTokenT<reco::CandidateView> m_token_input;
0366 
0367   std::vector<std::pair<double, double>> m_etaPhiPoints;
0368   std::vector<GlobalVector> m_directionPoints;
0369 
0370   float m_ptMin;
0371   float m_originRadius;
0372   float m_deltaEta_Cand;
0373   float m_deltaPhi_Cand;
0374   float m_deltaEta_Point;
0375   float m_deltaPhi_Point;
0376   bool m_precise;
0377   edm::EDGetTokenT<MeasurementTrackerEvent> m_token_measurementTracker;
0378   RectangularEtaPhiTrackingRegion::UseMeasurementTracker m_whereToUseMeasurementTracker;
0379   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> token_field;
0380   edm::ESGetToken<MultipleScatteringParametrisationMaker, TrackerMultipleScatteringRecord> token_msmaker;
0381   bool m_searchOpt;
0382 };
0383 
0384 #endif