SeedingTracksConverter.cc

Go to the documentation of this file.

#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"

#include "DataFormats/BTauReco/interface/SeedingTrackFeatures.h"
#include "DataFormats/BTauReco/interface/TrackPairFeatures.h"

#include "RecoBTag/FeatureTools/interface/TrackPairInfoBuilder.h"
#include "RecoBTag/FeatureTools/interface/SeedingTrackInfoBuilder.h"

#include "DataFormats/PatCandidates/interface/Jet.h"

#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/IPTools/interface/IPTools.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "RecoVertex/VertexTools/interface/VertexDistance3D.h"

#include "RecoBTag/TrackProbability/interface/HistogramProbabilityEstimator.h"

#include "RecoBTag/FeatureTools/interface/SeedingTracksConverter.h"

namespace btagbtvdeep {

  void seedingTracksToFeatures(const std::vector<reco::TransientTrack>& selectedTracks,
                               const std::vector<float>& masses,
                               const reco::Jet& jet,
                               const reco::Vertex& pv,
                               HistogramProbabilityEstimator* probabilityEstimator,
                               bool computeProbabilities,
                               std::vector<btagbtvdeep::SeedingTrackFeatures>& seedingT_features_vector)

  {
    GlobalVector jetdirection(jet.px(), jet.py(), jet.pz());
    GlobalPoint pvp(pv.x(), pv.y(), pv.z());

    std::multimap<double, std::pair<btagbtvdeep::SeedingTrackInfoBuilder, std::vector<btagbtvdeep::TrackPairFeatures>>>
        sortedSeedsMap;
    std::multimap<double, btagbtvdeep::TrackPairInfoBuilder> sortedNeighboursMap;

    std::vector<btagbtvdeep::TrackPairFeatures> tp_features_vector;

    sortedSeedsMap.clear();
    seedingT_features_vector.clear();

    std::vector<std::pair<bool, Measurement1D>> absIP3D(selectedTracks.size());
    std::vector<std::pair<bool, Measurement1D>> absIP2D(selectedTracks.size());
    std::vector<bool> absIP3D_filled(selectedTracks.size(), false);
    std::vector<bool> absIP2D_filled(selectedTracks.size(), false);

    unsigned int selTrackCount = 0;

    for (auto const& it : selectedTracks) {
      selTrackCount += 1;
      sortedNeighboursMap.clear();
      tp_features_vector.clear();

      if (reco::deltaR(it.track(), jet) > 0.4)
        continue;

      std::pair<bool, Measurement1D> ip = IPTools::absoluteImpactParameter3D(it, pv);

      absIP3D[selTrackCount - 1] = ip;
      absIP3D_filled[selTrackCount - 1] = true;

      std::pair<double, Measurement1D> jet_dist = IPTools::jetTrackDistance(it, jetdirection, pv);
      TrajectoryStateOnSurface closest =
          IPTools::closestApproachToJet(it.impactPointState(), pv, jetdirection, it.field());
      float length = 999;
      if (closest.isValid())
        length = (closest.globalPosition() - pvp).mag();

      if (!(ip.first && ip.second.value() >= 0.0 && ip.second.significance() >= 1.0 && ip.second.value() <= 9999. &&
            ip.second.significance() <= 9999. && it.track().normalizedChi2() < 5. &&
            std::fabs(it.track().dxy(pv.position())) < 2 && std::fabs(it.track().dz(pv.position())) < 17 &&
            jet_dist.second.value() < 0.07 && length < 5.))
        continue;

      std::pair<bool, Measurement1D> ip2d = IPTools::absoluteTransverseImpactParameter(it, pv);

      absIP2D[selTrackCount - 1] = ip2d;
      absIP2D_filled[selTrackCount - 1] = true;

      btagbtvdeep::SeedingTrackInfoBuilder seedInfo;
      seedInfo.buildSeedingTrackInfo(&(it),
                                     pv,
                                     jet,
                                     masses[selTrackCount - 1],
                                     ip,
                                     ip2d,
                                     jet_dist.second.value(),
                                     length,
                                     probabilityEstimator,
                                     computeProbabilities);

      unsigned int neighbourTrackCount = 0;

      for (auto const& tt : selectedTracks) {
        neighbourTrackCount += 1;

        if (neighbourTrackCount == selTrackCount)
          continue;
        if (std::fabs(pv.z() - tt.track().vz()) > 0.1)
          continue;

        //avoid calling IPs twice
        if (!absIP2D_filled[neighbourTrackCount - 1]) {
          absIP2D[neighbourTrackCount - 1] = IPTools::absoluteTransverseImpactParameter(tt, pv);
          absIP2D_filled[neighbourTrackCount - 1] = true;
        }

        if (!absIP3D_filled[neighbourTrackCount - 1]) {
          absIP3D[neighbourTrackCount - 1] = IPTools::absoluteImpactParameter3D(tt, pv);
          absIP3D_filled[neighbourTrackCount - 1] = true;
        }

        std::pair<bool, Measurement1D> t_ip = absIP3D[neighbourTrackCount - 1];
        std::pair<bool, Measurement1D> t_ip2d = absIP2D[neighbourTrackCount - 1];

        btagbtvdeep::TrackPairInfoBuilder trackPairInfo;
        trackPairInfo.buildTrackPairInfo(&(it), &(tt), pv, masses[neighbourTrackCount - 1], jetdirection, t_ip, t_ip2d);
        sortedNeighboursMap.insert(std::make_pair(trackPairInfo.pca_distance(), trackPairInfo));
      }

      int max_counter = 0;

      for (auto const& im : sortedNeighboursMap) {
        if (max_counter >= 20)
          break;
        btagbtvdeep::TrackPairFeatures tp_features;

        auto const& tp = im.second;

        tp_features.pt = (tp.track_pt() == 0) ? 0 : 1.0 / tp.track_pt();
        tp_features.eta = tp.track_eta();
        tp_features.phi = tp.track_phi();
        tp_features.mass = tp.track_candMass();
        tp_features.dz = logWithOffset(tp.track_dz());
        tp_features.dxy = logWithOffset(tp.track_dxy());
        tp_features.ip3D = log(tp.track_ip3d());
        tp_features.sip3D = log(tp.track_ip3dSig());
        tp_features.ip2D = log(tp.track_ip2d());
        tp_features.sip2D = log(tp.track_ip2dSig());
        tp_features.distPCA = log(tp.pca_distance());
        tp_features.dsigPCA = log(tp.pca_significance());
        tp_features.x_PCAonSeed = tp.pcaSeed_x();
        tp_features.y_PCAonSeed = tp.pcaSeed_y();
        tp_features.z_PCAonSeed = tp.pcaSeed_z();
        tp_features.xerr_PCAonSeed = tp.pcaSeed_xerr();
        tp_features.yerr_PCAonSeed = tp.pcaSeed_yerr();
        tp_features.zerr_PCAonSeed = tp.pcaSeed_zerr();
        tp_features.x_PCAonTrack = tp.pcaTrack_x();
        tp_features.y_PCAonTrack = tp.pcaTrack_y();
        tp_features.z_PCAonTrack = tp.pcaTrack_z();
        tp_features.xerr_PCAonTrack = tp.pcaTrack_xerr();
        tp_features.yerr_PCAonTrack = tp.pcaTrack_yerr();
        tp_features.zerr_PCAonTrack = tp.pcaTrack_zerr();
        tp_features.dotprodTrack = tp.dotprodTrack();
        tp_features.dotprodSeed = tp.dotprodSeed();
        tp_features.dotprodTrackSeed2D = tp.dotprodTrackSeed2D();
        tp_features.dotprodTrackSeed3D = tp.dotprodTrackSeed3D();
        tp_features.dotprodTrackSeedVectors2D = tp.dotprodTrackSeed2DV();
        tp_features.dotprodTrackSeedVectors3D = tp.dotprodTrackSeed3DV();
        tp_features.pvd_PCAonSeed = log(tp.pcaSeed_dist());
        tp_features.pvd_PCAonTrack = log(tp.pcaTrack_dist());
        tp_features.dist_PCAjetAxis = log(tp.pca_jetAxis_dist());
        tp_features.dotprod_PCAjetMomenta = tp.pca_jetAxis_dotprod();
        tp_features.deta_PCAjetDirs = log(tp.pca_jetAxis_dEta());
        tp_features.dphi_PCAjetDirs = tp.pca_jetAxis_dPhi();

        max_counter = max_counter + 1;
        tp_features_vector.push_back(tp_features);
      }

      sortedSeedsMap.insert(std::make_pair(-seedInfo.sip3d_Signed(), std::make_pair(seedInfo, tp_features_vector)));
    }

    int max_counter_seed = 0;

    for (auto const& im : sortedSeedsMap) {
      if (max_counter_seed >= 10)
        break;

      btagbtvdeep::SeedingTrackFeatures seed_features;

      auto const& seed = im.second.first;

      seed_features.nearTracks = im.second.second;
      seed_features.pt = (seed.pt() == 0) ? 0 : 1.0 / seed.pt();
      seed_features.eta = seed.eta();
      seed_features.phi = seed.phi();
      seed_features.mass = seed.mass();
      seed_features.dz = logWithOffset(seed.dz());
      seed_features.dxy = logWithOffset(seed.dxy());
      seed_features.ip3D = log(seed.ip3d());
      seed_features.sip3D = log(seed.sip3d());
      seed_features.ip2D = log(seed.ip2d());
      seed_features.sip2D = log(seed.sip2d());
      seed_features.signedIp3D = logWithOffset(seed.ip3d_Signed());
      seed_features.signedSip3D = logWithOffset(seed.sip3d_Signed());
      seed_features.signedIp2D = logWithOffset(seed.ip2d_Signed());
      seed_features.signedSip2D = logWithOffset(seed.sip2d_Signed());
      seed_features.trackProbability3D = seed.trackProbability3D();
      seed_features.trackProbability2D = seed.trackProbability2D();
      seed_features.chi2reduced = seed.chi2reduced();
      seed_features.nPixelHits = seed.nPixelHits();
      seed_features.nHits = seed.nHits();
      seed_features.jetAxisDistance = log(seed.jetAxisDistance());
      seed_features.jetAxisDlength = log(seed.jetAxisDlength());

      max_counter_seed = max_counter_seed + 1;
      seedingT_features_vector.push_back(seed_features);
    }

    if (sortedSeedsMap.size() < 10) {
      for (unsigned int i = sortedSeedsMap.size(); i < 10; i++) {
        std::vector<btagbtvdeep::TrackPairFeatures> tp_features_zeropad(20);
        btagbtvdeep::SeedingTrackFeatures seed_features_zeropad;
        seed_features_zeropad.nearTracks = tp_features_zeropad;
        seedingT_features_vector.push_back(seed_features_zeropad);
      }
    }
  }

}  // namespace btagbtvdeep