TSGForOI.cc

Go to the documentation of this file.

#include "RecoMuon/TrackerSeedGenerator/plugins/TSGForOI.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"

#include <memory>

using namespace edm;
using namespace std;

TSGForOI::TSGForOI(const edm::ParameterSet& iConfig)
    : src_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("src"))),
      numOfMaxSeedsParam_(iConfig.getParameter<uint32_t>("maxSeeds")),
      numOfLayersToTry_(iConfig.getParameter<int32_t>("layersToTry")),
      numOfHitsToTry_(iConfig.getParameter<int32_t>("hitsToTry")),
      fixedErrorRescalingForHits_(iConfig.getParameter<double>("fixedErrorRescaleFactorForHits")),
      fixedErrorRescalingForHitless_(iConfig.getParameter<double>("fixedErrorRescaleFactorForHitless")),
      adjustErrorsDynamicallyForHits_(iConfig.getParameter<bool>("adjustErrorsDynamicallyForHits")),
      adjustErrorsDynamicallyForHitless_(iConfig.getParameter<bool>("adjustErrorsDynamicallyForHitless")),
      estimatorName_(iConfig.getParameter<std::string>("estimator")),
      minEtaForTEC_(iConfig.getParameter<double>("minEtaForTEC")),
      maxEtaForTOB_(iConfig.getParameter<double>("maxEtaForTOB")),
      useHitLessSeeds_(iConfig.getParameter<bool>("UseHitLessSeeds")),
      useStereoLayersInTEC_(iConfig.getParameter<bool>("UseStereoLayersInTEC")),
      updator_(new KFUpdator()),
      measurementTrackerTag_(
          consumes<MeasurementTrackerEvent>(iConfig.getParameter<edm::InputTag>("MeasurementTrackerEvent"))),
      pT1_(iConfig.getParameter<double>("pT1")),
      pT2_(iConfig.getParameter<double>("pT2")),
      pT3_(iConfig.getParameter<double>("pT3")),
      eta1_(iConfig.getParameter<double>("eta1")),
      eta2_(iConfig.getParameter<double>("eta2")),
      SF1_(iConfig.getParameter<double>("SF1")),
      SF2_(iConfig.getParameter<double>("SF2")),
      SF3_(iConfig.getParameter<double>("SF3")),
      SF4_(iConfig.getParameter<double>("SF4")),
      SF5_(iConfig.getParameter<double>("SF5")),
      tsosDiff_(iConfig.getParameter<double>("tsosDiff")),
      propagatorName_(iConfig.getParameter<std::string>("propagatorName")),
      theCategory(string("Muon|RecoMuon|TSGForOI")) {
  produces<std::vector<TrajectorySeed> >();
}

TSGForOI::~TSGForOI() {}

void TSGForOI::produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
  unsigned int numOfMaxSeeds = numOfMaxSeedsParam_;
  unsigned int numSeedsMade = 0;
  bool analysedL2 = false;
  unsigned int layerCount = 0;

  Plane::PlanePointer dummyPlane = Plane::build(Plane::PositionType(), Plane::RotationType());

  edm::Handle<MeasurementTrackerEvent> measurementTrackerH;
  edm::ESHandle<Chi2MeasurementEstimatorBase> estimatorH;
  edm::ESHandle<MagneticField> magfieldH;
  edm::ESHandle<Propagator> propagatorAlongH;
  edm::ESHandle<Propagator> propagatorOppositeH;
  edm::ESHandle<TrackerGeometry> tmpTkGeometryH;
  edm::ESHandle<GlobalTrackingGeometry> geometryH;

  iSetup.get<IdealMagneticFieldRecord>().get(magfieldH);
  iSetup.get<TrackingComponentsRecord>().get(propagatorName_, propagatorOppositeH);
  iSetup.get<TrackingComponentsRecord>().get(propagatorName_, propagatorAlongH);
  iSetup.get<GlobalTrackingGeometryRecord>().get(geometryH);
  iSetup.get<TrackerDigiGeometryRecord>().get(tmpTkGeometryH);
  iSetup.get<TrackingComponentsRecord>().get(estimatorName_, estimatorH);
  iEvent.getByToken(measurementTrackerTag_, measurementTrackerH);

  edm::Handle<reco::TrackCollection> l2TrackCol;
  iEvent.getByToken(src_, l2TrackCol);

  //    The product:
  std::unique_ptr<std::vector<TrajectorySeed> > result(new std::vector<TrajectorySeed>());

  //    Get vector of Detector layers once:
  std::vector<BarrelDetLayer const*> const& tob = measurementTrackerH->geometricSearchTracker()->tobLayers();
  std::vector<ForwardDetLayer const*> const& tecPositive =
      tmpTkGeometryH->isThere(GeomDetEnumerators::P2OTEC)
          ? measurementTrackerH->geometricSearchTracker()->posTidLayers()
          : measurementTrackerH->geometricSearchTracker()->posTecLayers();
  std::vector<ForwardDetLayer const*> const& tecNegative =
      tmpTkGeometryH->isThere(GeomDetEnumerators::P2OTEC)
          ? measurementTrackerH->geometricSearchTracker()->negTidLayers()
          : measurementTrackerH->geometricSearchTracker()->negTecLayers();
  edm::ESHandle<TrackerTopology> tTopo_handle;
  iSetup.get<TrackerTopologyRcd>().get(tTopo_handle);
  const TrackerTopology* tTopo = tTopo_handle.product();

  //    Get the suitable propagators:
  std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(*propagatorAlongH, alongMomentum);
  std::unique_ptr<Propagator> propagatorOpposite = SetPropagationDirection(*propagatorOppositeH, oppositeToMomentum);

  edm::ESHandle<Propagator> SmartOpposite;
  edm::ESHandle<Propagator> SHPOpposite;
  iSetup.get<TrackingComponentsRecord>().get("hltESPSmartPropagatorAnyOpposite", SmartOpposite);
  iSetup.get<TrackingComponentsRecord>().get("hltESPSteppingHelixPropagatorOpposite", SHPOpposite);

  //    Loop over the L2's and make seeds for all of them:
  LogTrace(theCategory) << "TSGForOI::produce: Number of L2's: " << l2TrackCol->size();
  for (unsigned int l2TrackColIndex(0); l2TrackColIndex != l2TrackCol->size(); ++l2TrackColIndex) {
    const reco::TrackRef l2(l2TrackCol, l2TrackColIndex);
    std::unique_ptr<std::vector<TrajectorySeed> > out(new std::vector<TrajectorySeed>());
    LogTrace("TSGForOI") << "TSGForOI::produce: L2 muon pT, eta, phi --> " << l2->pt() << " , " << l2->eta() << " , "
                         << l2->phi() << endl;

    FreeTrajectoryState fts = trajectoryStateTransform::initialFreeState(*l2, magfieldH.product());
    dummyPlane->move(fts.position() - dummyPlane->position());
    TrajectoryStateOnSurface tsosAtIP = TrajectoryStateOnSurface(fts, *dummyPlane);
    LogTrace("TSGForOI") << "TSGForOI::produce: Created TSOSatIP: " << tsosAtIP << std::endl;

    // get the TSOS on the innermost layer of the L2.
    TrajectoryStateOnSurface tsosAtMuonSystem =
        trajectoryStateTransform::innerStateOnSurface(*l2, *geometryH, magfieldH.product());
    LogTrace("TSGForOI") << "TSGForOI::produce: Created TSOSatMuonSystem: " << tsosAtMuonSystem << endl;

    LogTrace("TSGForOI") << "TSGForOI::produce: Check the error of the L2 parameter and use hit seeds if big errors"
                         << endl;
    StateOnTrackerBound fromInside(propagatorAlong.get());
    TrajectoryStateOnSurface outerTkStateInside = fromInside(fts);

    StateOnTrackerBound fromOutside(&*SmartOpposite);
    TrajectoryStateOnSurface outerTkStateOutside = fromOutside(tsosAtMuonSystem);

    // for now only checking if the two positions (using updated and not-updated) agree withing certain extent,
    // will probably have to design something fancier for the future.
    auto dist = 0.0;
    bool useBoth = false;
    if (outerTkStateInside.isValid() && outerTkStateOutside.isValid()) {
      dist = match_Chi2(outerTkStateInside, outerTkStateOutside);
    }
    if (dist > tsosDiff_) {
      useBoth = true;
    }

    numSeedsMade = 0;
    analysedL2 = false;

    // if both TSOSs agree, use only the one at vertex, as it uses more information. If they do not agree, search for seed based on both

    //      BARREL
    if (std::abs(l2->eta()) < maxEtaForTOB_) {
      layerCount = 0;
      for (auto it = tob.rbegin(); it != tob.rend(); ++it) {  //This goes from outermost to innermost layer
        LogTrace("TSGForOI") << "TSGForOI::produce: looping in TOB layer " << layerCount << endl;
        findSeedsOnLayer(tTopo,
                         **it,
                         tsosAtIP,
                         *(propagatorAlong.get()),
                         *(propagatorOpposite.get()),
                         l2,
                         estimatorH,
                         measurementTrackerH,
                         numSeedsMade,
                         numOfMaxSeeds,
                         layerCount,
                         analysedL2,
                         out);
      }
      if (useBoth) {
        numSeedsMade = 0;
        for (auto it = tob.rbegin(); it != tob.rend(); ++it) {  //This goes from outermost to innermost layer
          LogTrace("TSGForOI") << "TSGForOI::produce: looping in TOB layer " << layerCount << endl;
          findSeedsOnLayer(tTopo,
                           **it,
                           tsosAtMuonSystem,
                           *(propagatorOpposite.get()),
                           *(propagatorOpposite.get()),
                           l2,
                           estimatorH,
                           measurementTrackerH,
                           numSeedsMade,
                           numOfMaxSeeds,
                           layerCount,
                           analysedL2,
                           out);
        }
      }
    }
    //      Reset Number of seeds if in overlap region:

    if (std::abs(l2->eta()) > minEtaForTEC_ && std::abs(l2->eta()) < maxEtaForTOB_) {
      numSeedsMade = 0;
    }

    //      ENDCAP+
    if (l2->eta() > minEtaForTEC_) {
      layerCount = 0;
      for (auto it = tecPositive.rbegin(); it != tecPositive.rend(); ++it) {
        LogTrace("TSGForOI") << "TSGForOI::produce: looping in TEC+ layer " << layerCount << endl;
        findSeedsOnLayer(tTopo,
                         **it,
                         tsosAtIP,
                         *(propagatorAlong.get()),
                         *(propagatorOpposite.get()),
                         l2,
                         estimatorH,
                         measurementTrackerH,
                         numSeedsMade,
                         numOfMaxSeeds,
                         layerCount,
                         analysedL2,
                         out);
      }
      if (useBoth) {
        numSeedsMade = 0;
        for (auto it = tecPositive.rbegin(); it != tecPositive.rend(); ++it) {
          LogTrace("TSGForOI") << "TSGForOI::produce: looping in TEC+ layer " << layerCount << endl;
          findSeedsOnLayer(tTopo,
                           **it,
                           tsosAtMuonSystem,
                           *(propagatorOpposite.get()),
                           *(propagatorOpposite.get()),
                           l2,
                           estimatorH,
                           measurementTrackerH,
                           numSeedsMade,
                           numOfMaxSeeds,
                           layerCount,
                           analysedL2,
                           out);
        }
      }
    }

    //      ENDCAP-
    if (l2->eta() < -minEtaForTEC_) {
      layerCount = 0;
      for (auto it = tecNegative.rbegin(); it != tecNegative.rend(); ++it) {
        LogTrace("TSGForOI") << "TSGForOI::produce: looping in TEC- layer " << layerCount << endl;
        findSeedsOnLayer(tTopo,
                         **it,
                         tsosAtIP,
                         *(propagatorAlong.get()),
                         *(propagatorOpposite.get()),
                         l2,
                         estimatorH,
                         measurementTrackerH,
                         numSeedsMade,
                         numOfMaxSeeds,
                         layerCount,
                         analysedL2,
                         out);
      }
      if (useBoth) {
        numSeedsMade = 0;
        for (auto it = tecNegative.rbegin(); it != tecNegative.rend(); ++it) {
          LogTrace("TSGForOI") << "TSGForOI::produce: looping in TEC- layer " << layerCount << endl;
          findSeedsOnLayer(tTopo,
                           **it,
                           tsosAtMuonSystem,
                           *(propagatorOpposite.get()),
                           *(propagatorOpposite.get()),
                           l2,
                           estimatorH,
                           measurementTrackerH,
                           numSeedsMade,
                           numOfMaxSeeds,
                           layerCount,
                           analysedL2,
                           out);
        }
      }
    }

    for (std::vector<TrajectorySeed>::iterator it = out->begin(); it != out->end(); ++it) {
      result->push_back(*it);
    }
  }  //L2Collection
  edm::LogInfo(theCategory) << "TSGForOI::produce: number of seeds made: " << result->size();

  iEvent.put(std::move(result));
}

void TSGForOI::findSeedsOnLayer(const TrackerTopology* tTopo,
                                const GeometricSearchDet& layer,
                                const TrajectoryStateOnSurface& tsosAtIP,
                                const Propagator& propagatorAlong,
                                const Propagator& propagatorOpposite,
                                const reco::TrackRef l2,
                                edm::ESHandle<Chi2MeasurementEstimatorBase>& estimatorH,
                                edm::Handle<MeasurementTrackerEvent>& measurementTrackerH,
                                unsigned int& numSeedsMade,
                                unsigned int& numOfMaxSeeds,
                                unsigned int& layerCount,
                                bool& analysedL2,
                                std::unique_ptr<std::vector<TrajectorySeed> >& out) const {
  if (numSeedsMade > numOfMaxSeeds)
    return;
  LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: numSeedsMade = " << numSeedsMade
                       << " , layerCount = " << layerCount << endl;

  double errorSFHits = 1.0;
  double errorSFHitless = 1.0;
  if (!adjustErrorsDynamicallyForHits_)
    errorSFHits = fixedErrorRescalingForHits_;
  else
    errorSFHits = calculateSFFromL2(l2);
  if (!adjustErrorsDynamicallyForHitless_)
    errorSFHitless = fixedErrorRescalingForHitless_;

  // Hitless:  TO Be discarded from here at some point.
  if (useHitLessSeeds_) {
    LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: Start hitless" << endl;
    std::vector<GeometricSearchDet::DetWithState> dets;
    layer.compatibleDetsV(tsosAtIP, propagatorAlong, *estimatorH, dets);
    if (!dets.empty()) {
      auto const& detOnLayer = dets.front().first;
      auto const& tsosOnLayer = dets.front().second;
      LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: tsosOnLayer " << tsosOnLayer << endl;
      if (!tsosOnLayer.isValid()) {
        edm::LogInfo(theCategory) << "ERROR!: Hitless TSOS is not valid!";
      } else {
        // calculate SF from L2 (only once -- if needed)
        if (!analysedL2 && adjustErrorsDynamicallyForHitless_) {
          errorSFHitless = calculateSFFromL2(l2);
          analysedL2 = true;
        }

        dets.front().second.rescaleError(errorSFHitless);
        PTrajectoryStateOnDet const& ptsod =
            trajectoryStateTransform::persistentState(tsosOnLayer, detOnLayer->geographicalId().rawId());
        TrajectorySeed::recHitContainer rHC;
        out->push_back(TrajectorySeed(ptsod, rHC, oppositeToMomentum));
        LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: TSOD (Hitless) done " << endl;
        numSeedsMade++;
      }
    }
  }

  // Hits:
  if (layerCount > numOfLayersToTry_)
    return;
  LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: Start Hits" << endl;
  if (makeSeedsFromHits(tTopo,
                        layer,
                        tsosAtIP,
                        *out,
                        propagatorAlong,
                        *measurementTrackerH,
                        estimatorH,
                        numSeedsMade,
                        errorSFHits,
                        l2->eta()))
    ++layerCount;
}

double TSGForOI::calculateSFFromL2(const reco::TrackRef track) const {
  double theSF = 1.0;
  //    L2 direction vs pT blowup - as was previously done:
  //    Split into 4 pT ranges: <pT1_, pT1_<pT2_, pT2_<pT3_, <pT4_: 13,30,70
  //    Split into 2 eta ranges for the middle two pT ranges: 1.0,1.4
  double abseta = std::abs(track->eta());
  if (track->pt() <= pT1_)
    theSF = SF1_;
  if (track->pt() > pT1_ && track->pt() <= pT2_) {
    if (abseta <= eta1_)
      theSF = SF3_;
    if (abseta > eta1_ && abseta <= eta2_)
      theSF = SF2_;
    if (abseta > eta2_)
      theSF = SF3_;
  }
  if (track->pt() > pT2_ && track->pt() <= pT3_) {
    if (abseta <= eta1_)
      theSF = SF5_;
    if (abseta > eta1_ && abseta <= eta2_)
      theSF = SF4_;
    if (abseta > eta2_)
      theSF = SF5_;
  }
  if (track->pt() > pT3_)
    theSF = SF5_;

  LogTrace(theCategory) << "TSGForOI::calculateSFFromL2: SF has been calculated as: " << theSF;
  return theSF;
}

int TSGForOI::makeSeedsFromHits(const TrackerTopology* tTopo,
                                const GeometricSearchDet& layer,
                                const TrajectoryStateOnSurface& tsosAtIP,
                                std::vector<TrajectorySeed>& out,
                                const Propagator& propagatorAlong,
                                const MeasurementTrackerEvent& measurementTracker,
                                edm::ESHandle<Chi2MeasurementEstimatorBase>& estimatorH,
                                unsigned int& numSeedsMade,
                                const double errorSF,
                                const double l2Eta) const {
  //        Error Rescaling:
  TrajectoryStateOnSurface onLayer(tsosAtIP);
  onLayer.rescaleError(errorSF);

  std::vector<GeometricSearchDet::DetWithState> dets;
  layer.compatibleDetsV(onLayer, propagatorAlong, *estimatorH, dets);

  //    Find Measurements on each DetWithState:
  LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: find measurements on each detWithState  " << dets.size() << endl;
  std::vector<TrajectoryMeasurement> meas;
  for (std::vector<GeometricSearchDet::DetWithState>::iterator it = dets.begin(); it != dets.end(); ++it) {
    MeasurementDetWithData det = measurementTracker.idToDet(it->first->geographicalId());
    if (det.isNull()) {
      continue;
    }
    if (!it->second.isValid())
      continue;  //Skip if TSOS is not valid

    std::vector<TrajectoryMeasurement> mymeas =
        det.fastMeasurements(it->second, onLayer, propagatorAlong, *estimatorH);  //Second TSOS is not used
    for (std::vector<TrajectoryMeasurement>::const_iterator it2 = mymeas.begin(), ed2 = mymeas.end(); it2 != ed2;
         ++it2) {
      if (it2->recHit()->isValid())
        meas.push_back(*it2);  //Only save those which are valid
    }
  }

  //    Update TSOS using TMs after sorting, then create Trajectory Seed and put into vector:
  LogTrace("TSGForOI")
      << "TSGForOI::findSeedsOnLayer: Update TSOS using TMs after sorting, then create Trajectory Seed, number of TM = "
      << meas.size() << endl;
  unsigned int found = 0;
  std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
  for (std::vector<TrajectoryMeasurement>::const_iterator it = meas.begin(); it != meas.end(); ++it) {
    TrajectoryStateOnSurface updatedTSOS = updator_->update(it->forwardPredictedState(), *it->recHit());
    LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: TSOS for TM " << found << endl;
    if (not updatedTSOS.isValid())
      continue;

    // CHECK if is StereoLayer:
    if (useStereoLayersInTEC_ && (fabs(l2Eta) > 0.8 && fabs(l2Eta) < 1.6)) {
      DetId detid = ((*it).recHit()->hit())->geographicalId();
      if (detid.subdetId() == StripSubdetector::TEC) {
        if (!tTopo->tecIsStereo(detid.rawId()))
          break;  // try another layer
      }
    }

    edm::OwnVector<TrackingRecHit> seedHits;
    seedHits.push_back(*it->recHit()->hit());
    PTrajectoryStateOnDet const& pstate =
        trajectoryStateTransform::persistentState(updatedTSOS, it->recHit()->geographicalId().rawId());
    TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
    LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: number of seedHits: " << seedHits.size() << endl;
    out.push_back(seed);
    found++;
    numSeedsMade++;
    if (found == numOfHitsToTry_)
      break;
  }
  return found;
}
//

double TSGForOI::match_Chi2(const TrajectoryStateOnSurface& tsos1, const TrajectoryStateOnSurface& tsos2) const {
  if (!tsos1.isValid() || !tsos2.isValid())
    return -1.;

  AlgebraicVector5 v(tsos1.localParameters().vector() - tsos2.localParameters().vector());
  AlgebraicSymMatrix55 m(tsos1.localError().matrix() + tsos2.localError().matrix());

  bool ierr = !m.Invert();

  if (ierr) {
    edm::LogInfo("TSGForOI") << "Error inverting covariance matrix";
    return -1;
  }

  double est = ROOT::Math::Similarity(v, m);

  return est;
}

void TSGForOI::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("src", edm::InputTag("hltL2Muons", "UpdatedAtVtx"));
  desc.add<int>("layersToTry", 1);
  desc.add<double>("fixedErrorRescaleFactorForHitless", 2.0);
  desc.add<int>("hitsToTry", 1);
  desc.add<bool>("adjustErrorsDynamicallyForHits", false);
  desc.add<bool>("adjustErrorsDynamicallyForHitless", false);
  desc.add<edm::InputTag>("MeasurementTrackerEvent", edm::InputTag("hltSiStripClusters"));
  desc.add<bool>("UseHitLessSeeds", true);
  desc.add<bool>("UseStereoLayersInTEC", false);
  desc.add<std::string>("estimator", "hltESPChi2MeasurementEstimator100");
  desc.add<double>("maxEtaForTOB", 1.2);
  desc.add<double>("minEtaForTEC", 0.8);
  desc.addUntracked<bool>("debug", true);
  desc.add<double>("fixedErrorRescaleFactorForHits", 2.0);
  desc.add<unsigned int>("maxSeeds", 1);
  desc.add<double>("pT1", 13.0);
  desc.add<double>("pT2", 30.0);
  desc.add<double>("pT3", 70.0);
  desc.add<double>("eta1", 1.0);
  desc.add<double>("eta2", 1.4);
  desc.add<double>("SF1", 3.0);
  desc.add<double>("SF2", 4.0);
  desc.add<double>("SF3", 5.0);
  desc.add<double>("SF4", 7.0);
  desc.add<double>("SF5", 10.0);
  desc.add<double>("tsosDiff", 0.03);
  desc.add<std::string>("propagatorName", "PropagatorWithMaterial");
  descriptions.add("TSGForOI", desc);
}

DEFINE_FWK_MODULE(TSGForOI);