GlobalTrackingTools/src/MuonTrackingRegionBuilder.cc

0001 /**

0002  *  Class: MuonTrackingRegionBuilder

0003  *

0004  *  Build a TrackingRegion around a standalone muon

0005  *

0006  *  \author N. Neumeister   Purdue University

0007  *  \author A. Everett      Purdue University

0008  */
0009
0010 #include "RecoMuon/GlobalTrackingTools/interface/MuonTrackingRegionBuilder.h"
0011
0012 //-------------------------------

0013 // Collaborating Class Headers --

0014 //-------------------------------

0015
0016 #include "DataFormats/Common/interface/Handle.h"
0017
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022
0023 #include "MagneticField/Engine/interface/MagneticField.h"
0024 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0025 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0026 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0027 #include "RecoTracker/Record/interface/TrackerMultipleScatteringRecord.h"
0028 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.h"
0029
0030 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0031 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0032 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0033
0034 //

0035 // constructor

0036 //

0037 void MuonTrackingRegionBuilder::build(const edm::ParameterSet& par, edm::ConsumesCollector& iC) {
0038   // Adjust errors on Eta, Phi, Z

0039   theNsigmaEta = par.getParameter<double>("Rescale_eta");
0040   theNsigmaPhi = par.getParameter<double>("Rescale_phi");
0041   theNsigmaDz = par.getParameter<double>("Rescale_Dz");
0042
0043   // Upper limits parameters

0044   theEtaRegionPar1 = par.getParameter<double>("EtaR_UpperLimit_Par1");
0045   theEtaRegionPar2 = par.getParameter<double>("EtaR_UpperLimit_Par2");
0046   thePhiRegionPar1 = par.getParameter<double>("PhiR_UpperLimit_Par1");
0047   thePhiRegionPar2 = par.getParameter<double>("PhiR_UpperLimit_Par2");
0048
0049   // Flag to switch to use Vertices instead of BeamSpot

0050   useVertex = par.getParameter<bool>("UseVertex");
0051
0052   // Flag to use fixed limits for Eta, Phi, Z, pT

0053   useFixedZ = par.getParameter<bool>("Z_fixed");
0054   useFixedPt = par.getParameter<bool>("Pt_fixed");
0055   useFixedPhi = par.getParameter<bool>("Phi_fixed");
0056   useFixedEta = par.getParameter<bool>("Eta_fixed");
0057
0058   // Minimum value for pT

0059   thePtMin = par.getParameter<double>("Pt_min");
0060
0061   // Minimum value for Phi

0062   thePhiMin = par.getParameter<double>("Phi_min");
0063
0064   // Minimum value for Eta

0065   theEtaMin = par.getParameter<double>("Eta_min");
0066
0067   // The static region size along the Z direction

0068   theHalfZ = par.getParameter<double>("DeltaZ");
0069
0070   // The transverse distance of the region from the BS/PV

0071   theDeltaR = par.getParameter<double>("DeltaR");
0072
0073   // The static region size in Eta

0074   theDeltaEta = par.getParameter<double>("DeltaEta");
0075
0076   // The static region size in Phi

0077   theDeltaPhi = par.getParameter<double>("DeltaPhi");
0078
0079   // Maximum number of regions to build when looping over Muons

0080   theMaxRegions = par.getParameter<int>("maxRegions");
0081
0082   // Flag to use precise??

0083   thePrecise = par.getParameter<bool>("precise");
0084
0085   // perigee reference point ToDo: Check this

0086   theOnDemand = RectangularEtaPhiTrackingRegion::intToUseMeasurementTracker(par.getParameter<int>("OnDemand"));
0087   if (theOnDemand != RectangularEtaPhiTrackingRegion::UseMeasurementTracker::kNever) {
0088     theMeasurementTrackerToken =
0089         iC.consumes<MeasurementTrackerEvent>(par.getParameter<edm::InputTag>("MeasurementTrackerName"));
0090   }
0091
0092   // Vertex collection and Beam Spot

0093   beamSpotToken = iC.consumes<reco::BeamSpot>(par.getParameter<edm::InputTag>("beamSpot"));
0094   vertexCollectionToken = iC.consumes<reco::VertexCollection>(par.getParameter<edm::InputTag>("vertexCollection"));
0095
0096   // Input muon collection

0097   inputCollectionToken = iC.consumes<reco::TrackCollection>(par.getParameter<edm::InputTag>("input"));
0098
0099   bfieldToken = iC.esConsumes();
0100   if (thePrecise) {
0101     msmakerToken = iC.esConsumes();
0102   }
0103 }
0104
0105 //

0106 // Member function to be compatible with TrackingRegionProducerFactory: create many ROI for many tracks

0107 //

0108 std::vector<std::unique_ptr<TrackingRegion>> MuonTrackingRegionBuilder::regions(const edm::Event& ev,
0109                                                                                 const edm::EventSetup& es) const {
0110   std::vector<std::unique_ptr<TrackingRegion>> result;
0111
0112   edm::Handle<reco::TrackCollection> tracks;
0113   ev.getByToken(inputCollectionToken, tracks);
0114
0115   int nRegions = 0;
0116   for (auto it = tracks->cbegin(), ed = tracks->cend(); it != ed && nRegions < theMaxRegions; ++it) {
0117     result.push_back(region(*it, ev, es));
0118     nRegions++;
0119   }
0120
0121   return result;
0122 }
0123
0124 //

0125 // Call region on Track from TrackRef

0126 //

0127 std::unique_ptr<RectangularEtaPhiTrackingRegion> MuonTrackingRegionBuilder::region(const reco::TrackRef& track) const {
0128   return region(*track);
0129 }
0130
0131 //

0132 // ToDo: Not sure if this is needed?

0133 //

0134 void MuonTrackingRegionBuilder::setEvent(const edm::Event& event, const edm::EventSetup& es) {
0135   theEvent = &event;
0136   theEventSetup = &es;
0137 }
0138
0139 //

0140 //  Main member function called to create the ROI

0141 //

0142 std::unique_ptr<RectangularEtaPhiTrackingRegion> MuonTrackingRegionBuilder::region(const reco::Track& staTrack,
0143                                                                                    const edm::Event& ev,
0144                                                                                    const edm::EventSetup& es) const {
0145   // get track momentum/direction at vertex

0146   const math::XYZVector& mom = staTrack.momentum();
0147   GlobalVector dirVector(mom.x(), mom.y(), mom.z());
0148   double pt = staTrack.pt();
0149
0150   // Fix for StandAlone tracks with low momentum

0151   const math::XYZVector& innerMomentum = staTrack.innerMomentum();
0152   GlobalVector forSmallMomentum(innerMomentum.x(), innerMomentum.y(), innerMomentum.z());
0153   if (staTrack.p() <= 1.5) {
0154     pt = std::abs(forSmallMomentum.perp());
0155   }
0156
0157   // initial vertex position - in the following it is replaced with beamspot/vertexing

0158   GlobalPoint vertexPos(0.0, 0.0, 0.0);
0159   // standard 15.9, if useVertex than use error from  vertex

0160   double deltaZ = theHalfZ;
0161
0162   // retrieve beam spot information

0163   edm::Handle<reco::BeamSpot> bs;
0164   bool bsHandleFlag = ev.getByToken(beamSpotToken, bs);
0165
0166   // check the validity, otherwise vertexing

0167   if (bsHandleFlag && bs.isValid() && !useVertex) {
0168     vertexPos = GlobalPoint(bs->x0(), bs->y0(), bs->z0());
0169     deltaZ = useFixedZ ? theHalfZ : bs->sigmaZ() * theNsigmaDz;
0170   } else {
0171     // get originZPos from list of reconstructed vertices (first or all)

0172     edm::Handle<reco::VertexCollection> vertexCollection;
0173     bool vtxHandleFlag = ev.getByToken(vertexCollectionToken, vertexCollection);
0174     // check if there exists at least one reconstructed vertex

0175     if (vtxHandleFlag && !vertexCollection->empty()) {
0176       // use the first vertex in the collection and assume it is the primary event vertex

0177       reco::VertexCollection::const_iterator vtx = vertexCollection->begin();
0178       if (!vtx->isFake() && vtx->isValid()) {
0179         vertexPos = GlobalPoint(vtx->x(), vtx->y(), vtx->z());
0180         deltaZ = useFixedZ ? theHalfZ : vtx->zError() * theNsigmaDz;
0181       }
0182     }
0183   }
0184
0185   // inizialize to the maximum possible value

0186   double deta = 0.4;
0187   double dphi = 0.6;
0188
0189   // evaluate the dynamical region if possible

0190   deta = theNsigmaEta * (staTrack.etaError());
0191   dphi = theNsigmaPhi * (staTrack.phiError());
0192
0193   // Region_Parametrizations to take into account possible L2 error matrix inconsistencies

0194   double region_dEta = 0;
0195   double region_dPhi = 0;
0196   double eta = 0;
0197   double phi = 0;
0198
0199   // eta, pt parametrization from MC study (circa 2009?)

0200   if (pt <= 10.) {
0201     // angular coefficients

0202     float acoeff_Phi = (thePhiRegionPar2 - thePhiRegionPar1) / 5;
0203     float acoeff_Eta = (theEtaRegionPar2 - theEtaRegionPar1) / 5;
0204
0205     eta = theEtaRegionPar1 + (acoeff_Eta) * (pt - 5.);
0206     phi = thePhiRegionPar1 + (acoeff_Phi) * (pt - 5.);
0207   }
0208   // parametrization 2nd bin in pt from MC study

0209   if (pt > 10. && pt < 100.) {
0210     eta = theEtaRegionPar2;
0211     phi = thePhiRegionPar2;
0212   }
0213   // parametrization 3rd bin in pt from MC study

0214   if (pt >= 100.) {
0215     // angular coefficients

0216     float acoeff_Phi = (thePhiRegionPar1 - thePhiRegionPar2) / 900;
0217     float acoeff_Eta = (theEtaRegionPar1 - theEtaRegionPar2) / 900;
0218
0219     eta = theEtaRegionPar2 + (acoeff_Eta) * (pt - 100.);
0220     phi = thePhiRegionPar2 + (acoeff_Phi) * (pt - 100.);
0221   }
0222
0223   double region_dPhi1 = std::min(phi, dphi);
0224   double region_dEta1 = std::min(eta, deta);
0225
0226   // decide to use either a parametrization or a dynamical region

0227   region_dPhi = useFixedPhi ? theDeltaPhi : std::max(thePhiMin, region_dPhi1);
0228   region_dEta = useFixedEta ? theDeltaEta : std::max(theEtaMin, region_dEta1);
0229
0230   float deltaR = theDeltaR;
0231   double minPt = useFixedPt ? thePtMin : std::max(thePtMin, pt * 0.6);
0232
0233   const MeasurementTrackerEvent* measurementTracker = nullptr;
0234   if (!theMeasurementTrackerToken.isUninitialized()) {
0235     edm::Handle<MeasurementTrackerEvent> hmte;
0236     ev.getByToken(theMeasurementTrackerToken, hmte);
0237     measurementTracker = hmte.product();
0238   }
0239
0240   const auto& bfield = es.getData(bfieldToken);
0241   const MultipleScatteringParametrisationMaker* msmaker = nullptr;
0242   if (thePrecise) {
0243     msmaker = &es.getData(msmakerToken);
0244   }
0245
0246   auto region = std::make_unique<RectangularEtaPhiTrackingRegion>(dirVector,
0247                                                                   vertexPos,
0248                                                                   minPt,
0249                                                                   deltaR,
0250                                                                   deltaZ,
0251                                                                   region_dEta,
0252                                                                   region_dPhi,
0253                                                                   bfield,
0254                                                                   msmaker,
0255                                                                   thePrecise,
0256                                                                   theOnDemand,
0257                                                                   measurementTracker);
0258
0259   LogDebug("MuonTrackingRegionBuilder") << "the region parameters are:\n"
0260                                         << "\n dirVector: " << dirVector << "\n vertexPos: " << vertexPos
0261                                         << "\n minPt: " << minPt << "\n deltaR:" << deltaR << "\n deltaZ:" << deltaZ
0262                                         << "\n region_dEta:" << region_dEta << "\n region_dPhi:" << region_dPhi
0263                                         << "\n on demand parameter: " << static_cast<int>(theOnDemand);
0264
0265   return region;
0266 }
0267
0268 void MuonTrackingRegionBuilder::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0269   {
0270     edm::ParameterSetDescription desc;
0271     fillDescriptionsOffline(desc);
0272     descriptions.add("MuonTrackingRegionBuilder", desc);
0273   }
0274   {
0275     edm::ParameterSetDescription desc;
0276     fillDescriptionsHLT(desc);
0277     descriptions.add("MuonTrackingRegionBuilderHLT", desc);
0278   }
0279   descriptions.setComment(
0280       "Build a TrackingRegion around a standalone muon. Options to define region around beamspot or primary vertex and "
0281       "dynamic regions are included.");
0282 }
0283 void MuonTrackingRegionBuilder::fillDescriptionsHLT(edm::ParameterSetDescription& desc) {
0284   desc.add<double>("EtaR_UpperLimit_Par1", 0.25);
0285   desc.add<double>("DeltaR", 0.2);
0286   desc.add<edm::InputTag>("beamSpot", edm::InputTag("hltOnlineBeamSpot"));
0287   desc.add<int>("OnDemand", -1);
0288   desc.add<edm::InputTag>("vertexCollection", edm::InputTag("pixelVertices"));
0289   desc.add<double>("Rescale_phi", 3.0);
0290   desc.add<bool>("Eta_fixed", false);
0291   desc.add<double>("Rescale_eta", 3.0);
0292   desc.add<double>("PhiR_UpperLimit_Par2", 0.2);
0293   desc.add<double>("Eta_min", 0.05);
0294   desc.add<bool>("Phi_fixed", false);
0295   desc.add<double>("Phi_min", 0.05);
0296   desc.add<double>("PhiR_UpperLimit_Par1", 0.6);
0297   desc.add<double>("EtaR_UpperLimit_Par2", 0.15);
0298   desc.add<edm::InputTag>("MeasurementTrackerName", edm::InputTag("hltESPMeasurementTracker"));
0299   desc.add<bool>("UseVertex", false);
0300   desc.add<double>("Rescale_Dz", 3.0);
0301   desc.add<bool>("Pt_fixed", false);
0302   desc.add<bool>("Z_fixed", true);
0303   desc.add<double>("Pt_min", 1.5);
0304   desc.add<double>("DeltaZ", 15.9);
0305   desc.add<double>("DeltaEta", 0.2);
0306   desc.add<double>("DeltaPhi", 0.2);
0307   desc.add<int>("maxRegions", 1);
0308   desc.add<bool>("precise", true);
0309   desc.add<edm::InputTag>("input", edm::InputTag("hltL2Muons", "UpdatedAtVtx"));
0310 }
0311
0312 void MuonTrackingRegionBuilder::fillDescriptionsOffline(edm::ParameterSetDescription& desc) {
0313   desc.add<double>("EtaR_UpperLimit_Par1", 0.25);
0314   desc.add<double>("DeltaR", 0.2);
0315   desc.add<edm::InputTag>("beamSpot", edm::InputTag(""));
0316   desc.add<int>("OnDemand", -1);
0317   desc.add<edm::InputTag>("vertexCollection", edm::InputTag(""));
0318   desc.add<double>("Rescale_phi", 3.0);
0319   desc.add<bool>("Eta_fixed", false);
0320   desc.add<double>("Rescale_eta", 3.0);
0321   desc.add<double>("PhiR_UpperLimit_Par2", 0.2);
0322   desc.add<double>("Eta_min", 0.05);
0323   desc.add<bool>("Phi_fixed", false);
0324   desc.add<double>("Phi_min", 0.05);
0325   desc.add<double>("PhiR_UpperLimit_Par1", 0.6);
0326   desc.add<double>("EtaR_UpperLimit_Par2", 0.15);
0327   desc.add<edm::InputTag>("MeasurementTrackerName", edm::InputTag(""));
0328   desc.add<bool>("UseVertex", false);
0329   desc.add<double>("Rescale_Dz", 3.0);
0330   desc.add<bool>("Pt_fixed", false);
0331   desc.add<bool>("Z_fixed", true);
0332   desc.add<double>("Pt_min", 1.5);
0333   desc.add<double>("DeltaZ", 15.9);
0334   desc.add<double>("DeltaEta", 0.2);
0335   desc.add<double>("DeltaPhi", 0.2);
0336   desc.add<int>("maxRegions", 1);
0337   desc.add<bool>("precise", true);
0338   desc.add<edm::InputTag>("input", edm::InputTag(""));
0339 }