AnalyticalTrackSelector.cc

Go to the documentation of this file.

#include <utility>
#include <vector>
#include <memory>
#include <algorithm>
#include <map>
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"

#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
#include "CommonTools/Utils/interface/StringCutObjectSelector.h"

#include "MultiTrackSelector.h"

using namespace reco;

class dso_hidden AnalyticalTrackSelector final : public MultiTrackSelector {
private:
public:
  explicit AnalyticalTrackSelector(const edm::ParameterSet& cfg);
  ~AnalyticalTrackSelector() override;

private:
  typedef math::XYZPoint Point;
  void run(edm::Event& evt, const edm::EventSetup& es) const override;

  bool copyExtras_;
  bool copyTrajectories_;
  double minEta_;
  double maxEta_;

  edm::EDGetTokenT<std::vector<Trajectory>> srcTraj_;
  edm::EDGetTokenT<TrajTrackAssociationCollection> srcTass_;
};

#include "DataFormats/Common/interface/ValueMap.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <Math/DistFunc.h>
#include "TMath.h"

AnalyticalTrackSelector::AnalyticalTrackSelector(const edm::ParameterSet& cfg) : MultiTrackSelector() {
  //Spoof the pset for each track selector!
  //Size is always 1!!!
  qualityToSet_.reserve(1);
  vtxNumber_.reserve(1);
  vertexCut_.reserve(1);
  res_par_.reserve(1);
  chi2n_par_.reserve(1);
  chi2n_no1Dmod_par_.reserve(1);
  d0_par1_.reserve(1);
  dz_par1_.reserve(1);
  d0_par2_.reserve(1);
  dz_par2_.reserve(1);
  applyAdaptedPVCuts_.reserve(1);
  max_d0_.reserve(1);
  max_z0_.reserve(1);
  nSigmaZ_.reserve(1);
  min_layers_.reserve(1);
  min_3Dlayers_.reserve(1);
  max_lostLayers_.reserve(1);
  min_hits_bypass_.reserve(1);
  applyAbsCutsIfNoPV_.reserve(1);
  max_d0NoPV_.reserve(1);
  max_z0NoPV_.reserve(1);
  preFilter_.reserve(1);
  max_relpterr_.reserve(1);
  min_nhits_.reserve(1);
  max_minMissHitOutOrIn_.reserve(1);
  max_lostHitFraction_.reserve(1);
  min_eta_.reserve(1);
  max_eta_.reserve(1);
  forest_.reserve(1);
  mvaType_.reserve(1);
  useMVA_.reserve(1);

  produces<edm::ValueMap<float>>("MVAVals");
  //foward compatibility
  produces<MVACollection>("MVAValues");
  useAnyMVA_ = false;
  forest_[0] = nullptr;
  if (cfg.exists("useAnyMVA"))
    useAnyMVA_ = cfg.getParameter<bool>("useAnyMVA");

  src_ = consumes<reco::TrackCollection>(cfg.getParameter<edm::InputTag>("src"));
  hSrc_ = consumes<TrackingRecHitCollection>(cfg.getParameter<edm::InputTag>("src"));
  beamspot_ = consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("beamspot"));
  useVertices_ = cfg.getParameter<bool>("useVertices");
  useVtxError_ = cfg.getParameter<bool>("useVtxError");
  if (useVertices_)
    vertices_ = consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("vertices"));
  copyExtras_ = cfg.getUntrackedParameter<bool>("copyExtras", false);
  copyTrajectories_ = cfg.getUntrackedParameter<bool>("copyTrajectories", false);
  if (copyTrajectories_) {
    srcTraj_ = consumes<std::vector<Trajectory>>(cfg.getParameter<edm::InputTag>("src"));
    srcTass_ = consumes<TrajTrackAssociationCollection>(cfg.getParameter<edm::InputTag>("src"));
  }

  qualityToSet_.push_back(TrackBase::undefQuality);
  // parameters for vertex selection
  vtxNumber_.push_back(useVertices_ ? cfg.getParameter<int32_t>("vtxNumber") : 0);
  vertexCut_.push_back(useVertices_ ? cfg.getParameter<std::string>("vertexCut") : "");
  //  parameters for adapted optimal cuts on chi2 and primary vertex compatibility
  res_par_.push_back(cfg.getParameter<std::vector<double>>("res_par"));
  chi2n_par_.push_back(cfg.getParameter<double>("chi2n_par"));
  chi2n_no1Dmod_par_.push_back(cfg.getParameter<double>("chi2n_no1Dmod_par"));
  d0_par1_.push_back(cfg.getParameter<std::vector<double>>("d0_par1"));
  dz_par1_.push_back(cfg.getParameter<std::vector<double>>("dz_par1"));
  d0_par2_.push_back(cfg.getParameter<std::vector<double>>("d0_par2"));
  dz_par2_.push_back(cfg.getParameter<std::vector<double>>("dz_par2"));

  // Boolean indicating if adapted primary vertex compatibility cuts are to be applied.
  applyAdaptedPVCuts_.push_back(cfg.getParameter<bool>("applyAdaptedPVCuts"));
  // Impact parameter absolute cuts.
  max_d0_.push_back(cfg.getParameter<double>("max_d0"));
  max_z0_.push_back(cfg.getParameter<double>("max_z0"));
  nSigmaZ_.push_back(cfg.getParameter<double>("nSigmaZ"));
  // Cuts on numbers of layers with hits/3D hits/lost hits.
  min_layers_.push_back(cfg.getParameter<uint32_t>("minNumberLayers"));
  min_3Dlayers_.push_back(cfg.getParameter<uint32_t>("minNumber3DLayers"));
  max_lostLayers_.push_back(cfg.getParameter<uint32_t>("maxNumberLostLayers"));
  min_hits_bypass_.push_back(cfg.getParameter<uint32_t>("minHitsToBypassChecks"));
  max_relpterr_.push_back(cfg.getParameter<double>("max_relpterr"));
  min_nhits_.push_back(cfg.getParameter<uint32_t>("min_nhits"));
  max_minMissHitOutOrIn_.push_back(
      cfg.existsAs<int32_t>("max_minMissHitOutOrIn") ? cfg.getParameter<int32_t>("max_minMissHitOutOrIn") : 99);
  max_lostHitFraction_.push_back(
      cfg.existsAs<double>("max_lostHitFraction") ? cfg.getParameter<double>("max_lostHitFraction") : 1.0);
  min_eta_.push_back(cfg.getParameter<double>("min_eta"));
  max_eta_.push_back(cfg.getParameter<double>("max_eta"));

  // Flag to apply absolute cuts if no PV passes the selection
  applyAbsCutsIfNoPV_.push_back(cfg.getParameter<bool>("applyAbsCutsIfNoPV"));
  keepAllTracks_.push_back(cfg.exists("keepAllTracks") ? cfg.getParameter<bool>("keepAllTracks") : false);

  setQualityBit_.push_back(false);
  std::string qualityStr = cfg.getParameter<std::string>("qualityBit");

  if (d0_par1_[0].size() != 2 || dz_par1_[0].size() != 2 || d0_par2_[0].size() != 2 || dz_par2_[0].size() != 2) {
    edm::LogError("MisConfiguration") << "vector of size less then 2";
    throw;
  }

  if (cfg.exists("qualityBit")) {
    std::string qualityStr = cfg.getParameter<std::string>("qualityBit");
    if (!qualityStr.empty()) {
      setQualityBit_[0] = true;
      qualityToSet_[0] = TrackBase::qualityByName(cfg.getParameter<std::string>("qualityBit"));
    }
  }

  if (keepAllTracks_[0] && !setQualityBit_[0])
    throw cms::Exception("Configuration")
        << "If you set 'keepAllTracks' to true, you must specify which qualityBit to set.\n";
  if (setQualityBit_[0] && (qualityToSet_[0] == TrackBase::undefQuality))
    throw cms::Exception("Configuration")
        << "You can't set the quality bit " << cfg.getParameter<std::string>("qualityBit")
        << " as it is 'undefQuality' or unknown.\n";
  if (applyAbsCutsIfNoPV_[0]) {
    max_d0NoPV_.push_back(cfg.getParameter<double>("max_d0NoPV"));
    max_z0NoPV_.push_back(cfg.getParameter<double>("max_z0NoPV"));
  } else {  //dummy values
    max_d0NoPV_.push_back(0.);
    max_z0NoPV_.push_back(0.);
  }

  if (useAnyMVA_) {
    bool thisMVA = false;
    if (cfg.exists("useMVA"))
      thisMVA = cfg.getParameter<bool>("useMVA");
    useMVA_.push_back(thisMVA);
    if (thisMVA) {
      double minVal = -1;
      if (cfg.exists("minMVA"))
        minVal = cfg.getParameter<double>("minMVA");
      min_MVA_.push_back(minVal);
      mvaType_.push_back(cfg.exists("mvaType") ? cfg.getParameter<std::string>("mvaType") : "Detached");
      forestLabel_.push_back(cfg.exists("GBRForestLabel") ? cfg.getParameter<std::string>("GBRForestLabel")
                                                          : "MVASelectorIter0");
      useMVAonly_.push_back(cfg.exists("useMVAonly") ? cfg.getParameter<bool>("useMVAonly") : false);
    } else {
      min_MVA_.push_back(-9999.0);
      useMVAonly_.push_back(false);
      mvaType_.push_back("Detached");
      forestLabel_.push_back("MVASelectorIter0");
    }
  } else {
    useMVA_.push_back(false);
    useMVAonly_.push_back(false);
    min_MVA_.push_back(-9999.0);
    mvaType_.push_back("Detached");
    forestLabel_.push_back("MVASelectorIter0");
  }

  std::string alias(cfg.getParameter<std::string>("@module_label"));
  if (copyExtras_) {
    produces<reco::TrackExtraCollection>().setBranchAlias(alias + "TrackExtras");
    produces<TrackingRecHitCollection>().setBranchAlias(alias + "RecHits");
  }
  if (copyTrajectories_) {
    produces<std::vector<Trajectory>>().setBranchAlias(alias + "Trajectories");
    produces<TrajTrackAssociationCollection>().setBranchAlias(alias + "TrajectoryTrackAssociations");
  }
  // TrackCollection refers to TrackingRechit and TrackExtra
  // collections, need to declare its production after them to work
  // around a rare race condition in framework scheduling
  produces<reco::TrackCollection>().setBranchAlias(alias + "Tracks");
}

AnalyticalTrackSelector::~AnalyticalTrackSelector() {}

void AnalyticalTrackSelector::run(edm::Event& evt, const edm::EventSetup& es) const {
  // storage....
  std::unique_ptr<reco::TrackCollection> selTracks_;
  std::unique_ptr<reco::TrackExtraCollection> selTrackExtras_;
  std::unique_ptr<TrackingRecHitCollection> selHits_;
  std::unique_ptr<std::vector<Trajectory>> selTrajs_;
  std::unique_ptr<std::vector<const Trajectory*>> selTrajPtrs_;
  std::unique_ptr<TrajTrackAssociationCollection> selTTAss_;
  reco::TrackRefProd rTracks_;
  reco::TrackExtraRefProd rTrackExtras_;
  TrackingRecHitRefProd rHits_;
  edm::RefProd<std::vector<Trajectory>> rTrajectories_;
  std::vector<reco::TrackRef> trackRefs_;

  using namespace std;
  using namespace edm;
  using namespace reco;

  Handle<TrackCollection> hSrcTrack;
  Handle<vector<Trajectory>> hTraj;
  Handle<vector<Trajectory>> hTrajP;
  Handle<TrajTrackAssociationCollection> hTTAss;

  // looking for the beam spot
  edm::Handle<reco::BeamSpot> hBsp;
  evt.getByToken(beamspot_, hBsp);
  reco::BeamSpot vertexBeamSpot;
  vertexBeamSpot = *hBsp;

  // Select good primary vertices for use in subsequent track selection
  const reco::VertexCollection dummyVtx;
  const reco::VertexCollection* vtxPtr = &dummyVtx;
  std::vector<Point> points;
  std::vector<float> vterr, vzerr;
  if (useVertices_) {
    vtxPtr = &evt.get(vertices_);
    selectVertices(0, *vtxPtr, points, vterr, vzerr);
    // Debug
    LogDebug("SelectVertex") << points.size() << " good pixel vertices";
  }

  // Get tracks
  evt.getByToken(src_, hSrcTrack);
  // get hits in track..
  Handle<TrackingRecHitCollection> hSrcHits;
  evt.getByToken(hSrc_, hSrcHits);
  const TrackingRecHitCollection& srcHits(*hSrcHits);

  selTracks_ = std::make_unique<TrackCollection>();
  rTracks_ = evt.getRefBeforePut<TrackCollection>();
  if (copyExtras_) {
    selTrackExtras_ = std::make_unique<TrackExtraCollection>();
    selHits_ = std::make_unique<TrackingRecHitCollection>();
    rHits_ = evt.getRefBeforePut<TrackingRecHitCollection>();
    rTrackExtras_ = evt.getRefBeforePut<TrackExtraCollection>();
  }

  if (copyTrajectories_)
    trackRefs_.resize(hSrcTrack->size());

  std::vector<float> mvaVals_(hSrcTrack->size(), -99.f);
  processMVA(evt, es, vertexBeamSpot, *vtxPtr, 0, mvaVals_, true);

  // Loop over tracks
  size_t current = 0;
  for (TrackCollection::const_iterator it = hSrcTrack->begin(), ed = hSrcTrack->end(); it != ed; ++it, ++current) {
    const Track& trk = *it;
    // Check if this track passes cuts

    LogTrace("TrackSelection") << "ready to check track with pt=" << trk.pt();

    float mvaVal = 0;
    if (useAnyMVA_)
      mvaVal = mvaVals_[current];
    bool ok = select(0, vertexBeamSpot, srcHits, trk, points, vterr, vzerr, mvaVal);
    if (!ok) {
      LogTrace("TrackSelection") << "track with pt=" << trk.pt() << " NOT selected";

      if (copyTrajectories_)
        trackRefs_[current] = reco::TrackRef();
      if (!keepAllTracks_[0])
        continue;
    }
    LogTrace("TrackSelection") << "track with pt=" << trk.pt() << " selected";
    selTracks_->push_back(Track(trk));  // clone and store
    if (ok && setQualityBit_[0]) {
      selTracks_->back().setQuality(qualityToSet_[0]);
      if (qualityToSet_[0] == TrackBase::tight) {
        selTracks_->back().setQuality(TrackBase::loose);
      } else if (qualityToSet_[0] == TrackBase::highPurity) {
        selTracks_->back().setQuality(TrackBase::loose);
        selTracks_->back().setQuality(TrackBase::tight);
      }
      if (!points.empty()) {
        if (qualityToSet_[0] == TrackBase::loose) {
          selTracks_->back().setQuality(TrackBase::looseSetWithPV);
        } else if (qualityToSet_[0] == TrackBase::highPurity) {
          selTracks_->back().setQuality(TrackBase::looseSetWithPV);
          selTracks_->back().setQuality(TrackBase::highPuritySetWithPV);
        }
      }
    }
    if (copyExtras_) {
      // TrackExtras
      selTrackExtras_->push_back(TrackExtra(trk.outerPosition(),
                                            trk.outerMomentum(),
                                            trk.outerOk(),
                                            trk.innerPosition(),
                                            trk.innerMomentum(),
                                            trk.innerOk(),
                                            trk.outerStateCovariance(),
                                            trk.outerDetId(),
                                            trk.innerStateCovariance(),
                                            trk.innerDetId(),
                                            trk.seedDirection(),
                                            trk.seedRef()));
      selTracks_->back().setExtra(TrackExtraRef(rTrackExtras_, selTrackExtras_->size() - 1));
      TrackExtra& tx = selTrackExtras_->back();
      tx.setResiduals(trk.residuals());
      // TrackingRecHits
      auto const firstHitIndex = selHits_->size();
      for (trackingRecHit_iterator hit = trk.recHitsBegin(); hit != trk.recHitsEnd(); ++hit) {
        selHits_->push_back((*hit)->clone());
      }
      tx.setHits(rHits_, firstHitIndex, selHits_->size() - firstHitIndex);
      tx.setTrajParams(trk.extra()->trajParams(), trk.extra()->chi2sX5());
    }
    if (copyTrajectories_) {
      trackRefs_[current] = TrackRef(rTracks_, selTracks_->size() - 1);
    }
  }
  if (copyTrajectories_) {
    Handle<vector<Trajectory>> hTraj;
    Handle<TrajTrackAssociationCollection> hTTAss;
    evt.getByToken(srcTass_, hTTAss);
    evt.getByToken(srcTraj_, hTraj);
    selTrajs_ = std::make_unique<std::vector<Trajectory>>();
    rTrajectories_ = evt.getRefBeforePut<vector<Trajectory>>();
    selTTAss_ = std::make_unique<TrajTrackAssociationCollection>(rTrajectories_, rTracks_);
    for (size_t i = 0, n = hTraj->size(); i < n; ++i) {
      Ref<vector<Trajectory>> trajRef(hTraj, i);
      TrajTrackAssociationCollection::const_iterator match = hTTAss->find(trajRef);
      if (match != hTTAss->end()) {
        const Ref<TrackCollection>& trkRef = match->val;
        short oldKey = static_cast<short>(trkRef.key());
        if (trackRefs_[oldKey].isNonnull()) {
          selTrajs_->push_back(Trajectory(*trajRef));
          selTTAss_->insert(Ref<vector<Trajectory>>(rTrajectories_, selTrajs_->size() - 1), trackRefs_[oldKey]);
        }
      }
    }
  }

  evt.put(std::move(selTracks_));
  if (copyExtras_) {
    evt.put(std::move(selTrackExtras_));
    evt.put(std::move(selHits_));
  }
  if (copyTrajectories_) {
    evt.put(std::move(selTrajs_));
    evt.put(std::move(selTTAss_));
  }
}

#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(AnalyticalTrackSelector);