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/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
#include "RecoTracker/Record/interface/TrackerMultipleScatteringRecord.h"
#include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.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"));

  bfieldToken = iC.esConsumes();
  if (thePrecise) {
    msmakerToken = iC.esConsumes();
  }
}

//
// 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, es));
    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, const edm::EventSetup& es) {
  theEvent = &event;
  theEventSetup = &es;
}

//
//  Main member function called to create the ROI
//
std::unique_ptr<RectangularEtaPhiTrackingRegion> MuonTrackingRegionBuilder::region(const reco::Track& staTrack,
                                                                                   const edm::Event& ev,
                                                                                   const edm::EventSetup& es) 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();
  }

  const auto& bfield = es.getData(bfieldToken);
  const MultipleScatteringParametrisationMaker* msmaker = nullptr;
  if (thePrecise) {
    msmaker = &es.getData(msmakerToken);
  }

  auto region = std::make_unique<RectangularEtaPhiTrackingRegion>(dirVector,
                                                                  vertexPos,
                                                                  minPt,
                                                                  deltaR,
                                                                  deltaZ,
                                                                  region_dEta,
                                                                  region_dPhi,
                                                                  bfield,
                                                                  msmaker,
                                                                  thePrecise,
                                                                  theOnDemand,
                                                                  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(""));
}