CombinedSVComputer.cc

Go to the documentation of this file.

#include "RecoBTag/SecondaryVertex/interface/CombinedSVComputer.h"

using namespace reco;

inline edm::ParameterSet CombinedSVComputer::dropDeltaR(const edm::ParameterSet &pset) const {
  edm::ParameterSet psetCopy(pset);
  psetCopy.addParameter<double>("jetDeltaRMax", 99999.0);
  return psetCopy;
}

CombinedSVComputer::CombinedSVComputer(const edm::ParameterSet &params)
    : trackFlip(params.getParameter<bool>("trackFlip")),
      vertexFlip(params.getParameter<bool>("vertexFlip")),
      charmCut(params.getParameter<double>("charmCut")),
      sortCriterium(TrackSorting::getCriterium(params.getParameter<std::string>("trackSort"))),
      trackSelector(params.getParameter<edm::ParameterSet>("trackSelection")),
      trackNoDeltaRSelector(dropDeltaR(params.getParameter<edm::ParameterSet>("trackSelection"))),
      trackPseudoSelector(params.getParameter<edm::ParameterSet>("trackPseudoSelection")),
      pseudoMultiplicityMin(params.getParameter<unsigned int>("pseudoMultiplicityMin")),
      trackMultiplicityMin(params.getParameter<unsigned int>("trackMultiplicityMin")),
      minTrackWeight(params.getParameter<double>("minimumTrackWeight")),
      useTrackWeights(params.getParameter<bool>("useTrackWeights")),
      vertexMassCorrection(params.getParameter<bool>("correctVertexMass")),
      pseudoVertexV0Filter(params.getParameter<edm::ParameterSet>("pseudoVertexV0Filter")),
      trackPairV0Filter(params.getParameter<edm::ParameterSet>("trackPairV0Filter")) {}

inline double CombinedSVComputer::flipValue(double value, bool vertex) const {
  return (vertex ? vertexFlip : trackFlip) ? -value : value;
}

inline CombinedSVComputer::IterationRange CombinedSVComputer::flipIterate(int size, bool vertex) const {
  IterationRange range;
  if (vertex ? vertexFlip : trackFlip) {
    range.begin = size - 1;
    range.end = -1;
    range.increment = -1;
  } else {
    range.begin = 0;
    range.end = size;
    range.increment = +1;
  }

  return range;
}

const btag::TrackIPData &CombinedSVComputer::threshTrack(const CandIPTagInfo &trackIPTagInfo,
                                                         const btag::SortCriteria sort,
                                                         const reco::Jet &jet,
                                                         const GlobalPoint &pv) const {
  const CandIPTagInfo::input_container &tracks = trackIPTagInfo.selectedTracks();
  const std::vector<btag::TrackIPData> &ipData = trackIPTagInfo.impactParameterData();
  std::vector<std::size_t> indices = trackIPTagInfo.sortedIndexes(sort);

  IterationRange range = flipIterate(indices.size(), false);
  TrackKinematics kin;
  range_for(i, range) {
    std::size_t idx = indices[i];
    const btag::TrackIPData &data = ipData[idx];
    const CandidatePtr &track = tracks[idx];

    if (!trackNoDeltaRSelector(track, data, jet, pv))
      continue;

    kin.add(track);
    if (kin.vectorSum().M() > charmCut)
      return data;
  }
  if (trackFlip) {
    static const btag::TrackIPData dummy = {GlobalPoint(),
                                            GlobalPoint(),
                                            Measurement1D(1.0, 1.0),
                                            Measurement1D(1.0, 1.0),
                                            Measurement1D(1.0, 1.0),
                                            Measurement1D(1.0, 1.0),
                                            0.};
    return dummy;
  } else {
    static const btag::TrackIPData dummy = {GlobalPoint(),
                                            GlobalPoint(),
                                            Measurement1D(-1.0, 1.0),
                                            Measurement1D(-1.0, 1.0),
                                            Measurement1D(-1.0, 1.0),
                                            Measurement1D(-1.0, 1.0),
                                            0.};
    return dummy;
  }
}

const btag::TrackIPData &CombinedSVComputer::threshTrack(const TrackIPTagInfo &trackIPTagInfo,
                                                         const btag::SortCriteria sort,
                                                         const reco::Jet &jet,
                                                         const GlobalPoint &pv) const {
  const edm::RefVector<TrackCollection> &tracks = trackIPTagInfo.selectedTracks();
  const std::vector<btag::TrackIPData> &ipData = trackIPTagInfo.impactParameterData();
  std::vector<std::size_t> indices = trackIPTagInfo.sortedIndexes(sort);

  IterationRange range = flipIterate(indices.size(), false);
  TrackKinematics kin;
  range_for(i, range) {
    std::size_t idx = indices[i];
    const btag::TrackIPData &data = ipData[idx];
    const Track &track = *tracks[idx];

    if (!trackNoDeltaRSelector(track, data, jet, pv))
      continue;

    kin.add(track);
    if (kin.vectorSum().M() > charmCut)
      return data;
  }

  if (trackFlip) {
    static const btag::TrackIPData dummy = {GlobalPoint(),
                                            GlobalPoint(),
                                            Measurement1D(1.0, 1.0),
                                            Measurement1D(1.0, 1.0),
                                            Measurement1D(1.0, 1.0),
                                            Measurement1D(1.0, 1.0),
                                            0.};
    return dummy;
  } else {
    static const btag::TrackIPData dummy = {GlobalPoint(),
                                            GlobalPoint(),
                                            Measurement1D(-1.0, 1.0),
                                            Measurement1D(-1.0, 1.0),
                                            Measurement1D(-1.0, 1.0),
                                            Measurement1D(-1.0, 1.0),
                                            0.};
    return dummy;
  }
}

TaggingVariableList CombinedSVComputer::operator()(const TrackIPTagInfo &ipInfo,
                                                   const SecondaryVertexTagInfo &svInfo) const {
  using namespace ROOT::Math;

  edm::RefToBase<Jet> jet = ipInfo.jet();
  math::XYZVector jetDir = jet->momentum().Unit();
  TaggingVariableList vars;

  TrackKinematics vertexKinematics;

  double vtx_track_ptSum = 0.;
  double vtx_track_ESum = 0.;

  // the following is specific depending on the type of vertex
  int vtx = -1;
  unsigned int numberofvertextracks = 0;

  IterationRange range = flipIterate(svInfo.nVertices(), true);
  range_for(i, range) {
    numberofvertextracks = numberofvertextracks + (svInfo.secondaryVertex(i)).nTracks();

    const Vertex &vertex = svInfo.secondaryVertex(i);
    bool hasRefittedTracks = vertex.hasRefittedTracks();
    for (reco::Vertex::trackRef_iterator track = vertex.tracks_begin(); track != vertex.tracks_end(); track++) {
      double w = vertex.trackWeight(*track);
      if (w < minTrackWeight)
        continue;
      if (hasRefittedTracks) {
        const Track actualTrack = vertex.refittedTrack(*track);
        vertexKinematics.add(actualTrack, w);
        vars.insert(btau::trackEtaRel, reco::btau::etaRel(jetDir, actualTrack.momentum()), true);
        if (vtx < 0)  // calculate this only for the first vertex
        {
          vtx_track_ptSum += std::sqrt(actualTrack.momentum().Perp2());
          vtx_track_ESum += std::sqrt(actualTrack.momentum().Mag2() + ROOT::Math::Square(ParticleMasses::piPlus));
        }
      } else {
        vertexKinematics.add(**track, w);
        vars.insert(btau::trackEtaRel, reco::btau::etaRel(jetDir, (*track)->momentum()), true);
        if (vtx < 0)  // calculate this only for the first vertex
        {
          vtx_track_ptSum += std::sqrt((*track)->momentum().Perp2());
          vtx_track_ESum += std::sqrt((*track)->momentum().Mag2() + ROOT::Math::Square(ParticleMasses::piPlus));
        }
      }
    }

    if (vtx < 0)
      vtx = i;
  }
  if (vtx >= 0) {
    vars.insert(btau::vertexNTracks, numberofvertextracks, true);
    vars.insert(btau::vertexFitProb, (svInfo.secondaryVertex(vtx)).normalizedChi2(), true);
  }

  // after we collected vertex information we let the common code complete the job
  fillCommonVariables(vars, vertexKinematics, ipInfo, svInfo, vtx_track_ptSum, vtx_track_ESum);

  vars.finalize();
  return vars;
}

TaggingVariableList CombinedSVComputer::operator()(const CandIPTagInfo &ipInfo,
                                                   const CandSecondaryVertexTagInfo &svInfo) const {
  using namespace ROOT::Math;

  edm::RefToBase<Jet> jet = ipInfo.jet();
  math::XYZVector jetDir = jet->momentum().Unit();
  TaggingVariableList vars;

  TrackKinematics vertexKinematics;

  double vtx_track_ptSum = 0.;
  double vtx_track_ESum = 0.;

  // the following is specific depending on the type of vertex
  int vtx = -1;
  unsigned int numberofvertextracks = 0;

  IterationRange range = flipIterate(svInfo.nVertices(), true);
  range_for(i, range) {
    numberofvertextracks = numberofvertextracks + (svInfo.secondaryVertex(i)).numberOfSourceCandidatePtrs();

    const reco::VertexCompositePtrCandidate &vertex = svInfo.secondaryVertex(i);
    const std::vector<CandidatePtr> &tracks = vertex.daughterPtrVector();
    for (std::vector<CandidatePtr>::const_iterator track = tracks.begin(); track != tracks.end(); ++track) {
      vertexKinematics.add(*track);
      vars.insert(btau::trackEtaRel, reco::btau::etaRel(jetDir, (*track)->momentum()), true);
      if (vtx < 0)  // calculate this only for the first vertex
      {
        vtx_track_ptSum += std::sqrt((*track)->momentum().Perp2());
        vtx_track_ESum += std::sqrt((*track)->momentum().Mag2() + ROOT::Math::Square(ParticleMasses::piPlus));
      }
    }

    if (vtx < 0)
      vtx = i;
  }
  if (vtx >= 0) {
    vars.insert(btau::vertexNTracks, numberofvertextracks, true);
    vars.insert(btau::vertexFitProb, (svInfo.secondaryVertex(vtx)).vertexNormalizedChi2(), true);
  }

  // after we collected vertex information we let the common code complete the job
  fillCommonVariables(vars, vertexKinematics, ipInfo, svInfo, vtx_track_ptSum, vtx_track_ESum);

  vars.finalize();
  return vars;
}