MuonTrackingRegionBuilder.cc

Go to the documentation of this file.

 
#include "RecoMuon/GlobalTrackingTools/interface/MuonTrackingRegionBuilder.h"
 
//-------------------------------
// Collaborating Class Headers --
//-------------------------------
 
#include "DataFormats/Common/interface/Handle.h"
 
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
 
#include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
 
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
#include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
 
//
// constructor
//
void MuonTrackingRegionBuilder::build(const edm::ParameterSet& par, edm::ConsumesCollector& iC) {
  // Adjust errors on Eta, Phi, Z
  theNsigmaEta = par.getParameter<double>("Rescale_eta");
  theNsigmaPhi = par.getParameter<double>("Rescale_phi");
  theNsigmaDz = par.getParameter<double>("Rescale_Dz");
 
  // Upper limits parameters
  theEtaRegionPar1 = par.getParameter<double>("EtaR_UpperLimit_Par1");
  theEtaRegionPar2 = par.getParameter<double>("EtaR_UpperLimit_Par2");
  thePhiRegionPar1 = par.getParameter<double>("PhiR_UpperLimit_Par1");
  thePhiRegionPar2 = par.getParameter<double>("PhiR_UpperLimit_Par2");
 
  // Flag to switch to use Vertices instead of BeamSpot
  useVertex = par.getParameter<bool>("UseVertex");
 
  // Flag to use fixed limits for Eta, Phi, Z, pT
  useFixedZ = par.getParameter<bool>("Z_fixed");
  useFixedPt = par.getParameter<bool>("Pt_fixed");
  useFixedPhi = par.getParameter<bool>("Phi_fixed");
  useFixedEta = par.getParameter<bool>("Eta_fixed");
 
  // Minimum value for pT
  thePtMin = par.getParameter<double>("Pt_min");
 
  // Minimum value for Phi
  thePhiMin = par.getParameter<double>("Phi_min");
 
  // Minimum value for Eta
  theEtaMin = par.getParameter<double>("Eta_min");
 
  // The static region size along the Z direction
  theHalfZ = par.getParameter<double>("DeltaZ");
 
  // The transverse distance of the region from the BS/PV
  theDeltaR = par.getParameter<double>("DeltaR");
 
  // The static region size in Eta
  theDeltaEta = par.getParameter<double>("DeltaEta");
 
  // The static region size in Phi
  theDeltaPhi = par.getParameter<double>("DeltaPhi");
 
  // Maximum number of regions to build when looping over Muons
  theMaxRegions = par.getParameter<int>("maxRegions");
 
  // Flag to use precise??
  thePrecise = par.getParameter<bool>("precise");
 
  // perigee reference point ToDo: Check this
  theOnDemand = RectangularEtaPhiTrackingRegion::intToUseMeasurementTracker(par.getParameter<int>("OnDemand"));
  if (theOnDemand != RectangularEtaPhiTrackingRegion::UseMeasurementTracker::kNever) {
    theMeasurementTrackerToken =
        iC.consumes<MeasurementTrackerEvent>(par.getParameter<edm::InputTag>("MeasurementTrackerName"));
  }
 
  // Vertex collection and Beam Spot
  beamSpotToken = iC.consumes<reco::BeamSpot>(par.getParameter<edm::InputTag>("beamSpot"));
  vertexCollectionToken = iC.consumes<reco::VertexCollection>(par.getParameter<edm::InputTag>("vertexCollection"));
 
  // Input muon collection
  inputCollectionToken = iC.consumes<reco::TrackCollection>(par.getParameter<edm::InputTag>("input"));
}
 
//
// Member function to be compatible with TrackingRegionProducerFactory: create many ROI for many tracks
//
std::vector<std::unique_ptr<TrackingRegion>> MuonTrackingRegionBuilder::regions(const edm::Event& ev,
                                                                                const edm::EventSetup& es) const {
  std::vector<std::unique_ptr<TrackingRegion>> result;
 
  edm::Handle<reco::TrackCollection> tracks;
  ev.getByToken(inputCollectionToken, tracks);
 
  int nRegions = 0;
  for (auto it = tracks->cbegin(), ed = tracks->cend(); it != ed && nRegions < theMaxRegions; ++it) {
    result.push_back(region(*it, ev));
    nRegions++;
  }
 
  return result;
}
 
//
// Call region on Track from TrackRef
//
std::unique_ptr<RectangularEtaPhiTrackingRegion> MuonTrackingRegionBuilder::region(const reco::TrackRef& track) const {
  return region(*track);
}
 
//
// ToDo: Not sure if this is needed?
//
void MuonTrackingRegionBuilder::setEvent(const edm::Event& event) { theEvent = &event; }
 
//
//  Main member function called to create the ROI
//
std::unique_ptr<RectangularEtaPhiTrackingRegion> MuonTrackingRegionBuilder::region(const reco::Track& staTrack,
                                                                                   const edm::Event& ev) const {
  // get track momentum/direction at vertex
  const math::XYZVector& mom = staTrack.momentum();
  GlobalVector dirVector(mom.x(), mom.y(), mom.z());
  double pt = staTrack.pt();
 
  // Fix for StandAlone tracks with low momentum
  const math::XYZVector& innerMomentum = staTrack.innerMomentum();
  GlobalVector forSmallMomentum(innerMomentum.x(), innerMomentum.y(), innerMomentum.z());
  if (staTrack.p() <= 1.5) {
    pt = std::abs(forSmallMomentum.perp());
  }
 
  // initial vertex position - in the following it is replaced with beamspot/vertexing
  GlobalPoint vertexPos(0.0, 0.0, 0.0);
  // standard 15.9, if useVertex than use error from  vertex
  double deltaZ = theHalfZ;
 
  // retrieve beam spot information
  edm::Handle<reco::BeamSpot> bs;
  bool bsHandleFlag = ev.getByToken(beamSpotToken, bs);
 
  // check the validity, otherwise vertexing
  if (bsHandleFlag && bs.isValid() && !useVertex) {
    vertexPos = GlobalPoint(bs->x0(), bs->y0(), bs->z0());
    deltaZ = useFixedZ ? theHalfZ : bs->sigmaZ() * theNsigmaDz;
  } else {
    // get originZPos from list of reconstructed vertices (first or all)
    edm::Handle<reco::VertexCollection> vertexCollection;
    bool vtxHandleFlag = ev.getByToken(vertexCollectionToken, vertexCollection);
    // check if there exists at least one reconstructed vertex
    if (vtxHandleFlag && !vertexCollection->empty()) {
      // use the first vertex in the collection and assume it is the primary event vertex
      reco::VertexCollection::const_iterator vtx = vertexCollection->begin();
      if (!vtx->isFake() && vtx->isValid()) {
        vertexPos = GlobalPoint(vtx->x(), vtx->y(), vtx->z());
        deltaZ = useFixedZ ? theHalfZ : vtx->zError() * theNsigmaDz;
      }
    }
  }
 
  // inizialize to the maximum possible value
  double deta = 0.4;
  double dphi = 0.6;
 
  // evaluate the dynamical region if possible
  deta = theNsigmaEta * (staTrack.etaError());
  dphi = theNsigmaPhi * (staTrack.phiError());
 
  // Region_Parametrizations to take into account possible L2 error matrix inconsistencies
  double region_dEta = 0;
  double region_dPhi = 0;
  double eta = 0;
  double phi = 0;
 
  // eta, pt parametrization from MC study (circa 2009?)
  if (pt <= 10.) {
    // angular coefficients
    float acoeff_Phi = (thePhiRegionPar2 - thePhiRegionPar1) / 5;
    float acoeff_Eta = (theEtaRegionPar2 - theEtaRegionPar1) / 5;
 
    eta = theEtaRegionPar1 + (acoeff_Eta) * (pt - 5.);
    phi = thePhiRegionPar1 + (acoeff_Phi) * (pt - 5.);
  }
  // parametrization 2nd bin in pt from MC study
  if (pt > 10. && pt < 100.) {
    eta = theEtaRegionPar2;
    phi = thePhiRegionPar2;
  }
  // parametrization 3rd bin in pt from MC study
  if (pt >= 100.) {
    // angular coefficients
    float acoeff_Phi = (thePhiRegionPar1 - thePhiRegionPar2) / 900;
    float acoeff_Eta = (theEtaRegionPar1 - theEtaRegionPar2) / 900;
 
    eta = theEtaRegionPar2 + (acoeff_Eta) * (pt - 100.);
    phi = thePhiRegionPar2 + (acoeff_Phi) * (pt - 100.);
  }
 
  double region_dPhi1 = std::min(phi, dphi);
  double region_dEta1 = std::min(eta, deta);
 
  // decide to use either a parametrization or a dynamical region
  region_dPhi = useFixedPhi ? theDeltaPhi : std::max(thePhiMin, region_dPhi1);
  region_dEta = useFixedEta ? theDeltaEta : std::max(theEtaMin, region_dEta1);
 
  float deltaR = theDeltaR;
  double minPt = useFixedPt ? thePtMin : std::max(thePtMin, pt * 0.6);
 
  const MeasurementTrackerEvent* measurementTracker = nullptr;
  if (!theMeasurementTrackerToken.isUninitialized()) {
    edm::Handle<MeasurementTrackerEvent> hmte;
    ev.getByToken(theMeasurementTrackerToken, hmte);
    measurementTracker = hmte.product();
  }
 
  auto region = std::make_unique<RectangularEtaPhiTrackingRegion>(
      dirVector, vertexPos, minPt, deltaR, deltaZ, region_dEta, region_dPhi, theOnDemand, thePrecise, measurementTracker);
 
  LogDebug("MuonTrackingRegionBuilder") << "the region parameters are:\n"
                                        << "\n dirVector: " << dirVector << "\n vertexPos: " << vertexPos
                                        << "\n minPt: " << minPt << "\n deltaR:" << deltaR << "\n deltaZ:" << deltaZ
                                        << "\n region_dEta:" << region_dEta << "\n region_dPhi:" << region_dPhi
                                        << "\n on demand parameter: " << static_cast<int>(theOnDemand);
 
  return region;
}
 
void MuonTrackingRegionBuilder::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  {
    edm::ParameterSetDescription desc;
    fillDescriptionsOffline(desc);
    descriptions.add("MuonTrackingRegionBuilder", desc);
  }
  {
    edm::ParameterSetDescription desc;
    fillDescriptionsHLT(desc);
    descriptions.add("MuonTrackingRegionBuilderHLT", desc);
  }
  descriptions.setComment(
      "Build a TrackingRegion around a standalone muon. Options to define region around beamspot or primary vertex and "
      "dynamic regions are included.");
}
void MuonTrackingRegionBuilder::fillDescriptionsHLT(edm::ParameterSetDescription& desc) {
  desc.add<double>("EtaR_UpperLimit_Par1", 0.25);
  desc.add<double>("DeltaR", 0.2);
  desc.add<edm::InputTag>("beamSpot", edm::InputTag("hltOnlineBeamSpot"));
  desc.add<int>("OnDemand", -1);
  desc.add<edm::InputTag>("vertexCollection", edm::InputTag("pixelVertices"));
  desc.add<double>("Rescale_phi", 3.0);
  desc.add<bool>("Eta_fixed", false);
  desc.add<double>("Rescale_eta", 3.0);
  desc.add<double>("PhiR_UpperLimit_Par2", 0.2);
  desc.add<double>("Eta_min", 0.05);
  desc.add<bool>("Phi_fixed", false);
  desc.add<double>("Phi_min", 0.05);
  desc.add<double>("PhiR_UpperLimit_Par1", 0.6);
  desc.add<double>("EtaR_UpperLimit_Par2", 0.15);
  desc.add<edm::InputTag>("MeasurementTrackerName", edm::InputTag("hltESPMeasurementTracker"));
  desc.add<bool>("UseVertex", false);
  desc.add<double>("Rescale_Dz", 3.0);
  desc.add<bool>("Pt_fixed", false);
  desc.add<bool>("Z_fixed", true);
  desc.add<double>("Pt_min", 1.5);
  desc.add<double>("DeltaZ", 15.9);
  desc.add<double>("DeltaEta", 0.2);
  desc.add<double>("DeltaPhi", 0.2);
  desc.add<int>("maxRegions", 1);
  desc.add<bool>("precise", true);
  desc.add<edm::InputTag>("input", edm::InputTag("hltL2Muons", "UpdatedAtVtx"));
}
 
void MuonTrackingRegionBuilder::fillDescriptionsOffline(edm::ParameterSetDescription& desc) {
  desc.add<double>("EtaR_UpperLimit_Par1", 0.25);
  desc.add<double>("DeltaR", 0.2);
  desc.add<edm::InputTag>("beamSpot", edm::InputTag(""));
  desc.add<int>("OnDemand", -1);
  desc.add<edm::InputTag>("vertexCollection", edm::InputTag(""));
  desc.add<double>("Rescale_phi", 3.0);
  desc.add<bool>("Eta_fixed", false);
  desc.add<double>("Rescale_eta", 3.0);
  desc.add<double>("PhiR_UpperLimit_Par2", 0.2);
  desc.add<double>("Eta_min", 0.05);
  desc.add<bool>("Phi_fixed", false);
  desc.add<double>("Phi_min", 0.05);
  desc.add<double>("PhiR_UpperLimit_Par1", 0.6);
  desc.add<double>("EtaR_UpperLimit_Par2", 0.15);
  desc.add<edm::InputTag>("MeasurementTrackerName", edm::InputTag(""));
  desc.add<bool>("UseVertex", false);
  desc.add<double>("Rescale_Dz", 3.0);
  desc.add<bool>("Pt_fixed", false);
  desc.add<bool>("Z_fixed", true);
  desc.add<double>("Pt_min", 1.5);
  desc.add<double>("DeltaZ", 15.9);
  desc.add<double>("DeltaEta", 0.2);
  desc.add<double>("DeltaPhi", 0.2);
  desc.add<int>("maxRegions", 1);
  desc.add<bool>("precise", true);
  desc.add<edm::InputTag>("input", edm::InputTag(""));
}