Introduce DeepDoubleX V2

DataFormats/BTauReco/interface/ChargedCandidateFeatures.h

@@ -6,6 +6,9 @@ namespace btagbtvdeep {
   class ChargedCandidateFeatures {
   public:
     float ptrel;
+    float ptrel_noclip;
+    float erel;
+    float etarel;
     float puppiw;
     float vtx_ass;
 
@@ -22,6 +25,14 @@ namespace btagbtvdeep {
 
     float btagPf_trackJetDistVal;
 
+    float drsubjet1;
+    float drsubjet2;
+    float dxy;
+    float dxysig;
+    float dz;
+    float dzsig;
+    float deltaR;
+
     float chi2;
     float quality;
 

DataFormats/BTauReco/interface/DeepDoubleXFeatures.h

@@ -7,6 +7,7 @@
 #include "DataFormats/BTauReco/interface/SecondaryVertexFeatures.h"
 #include "DataFormats/BTauReco/interface/BoostedDoubleSVTagInfoFeatures.h"
 #include "DataFormats/BTauReco/interface/ChargedCandidateFeatures.h"
+#include "DataFormats/BTauReco/interface/NeutralCandidateFeatures.h"
 
 namespace btagbtvdeep {
 
@@ -19,10 +20,12 @@ namespace btagbtvdeep {
     JetFeatures jet_features;
     BoostedDoubleSVTagInfoFeatures tag_info_features;
 
-    std::vector<SecondaryVertexFeatures> sv_features;
-
     std::vector<ChargedCandidateFeatures> c_pf_features;
 
+    std::vector<NeutralCandidateFeatures> n_pf_features;
+
+    std::vector<SecondaryVertexFeatures> sv_features;
+
     std::size_t npv;  // used by deep flavour
 
   private:

DataFormats/BTauReco/interface/NeutralCandidateFeatures.h

@@ -6,9 +6,15 @@ namespace btagbtvdeep {
   class NeutralCandidateFeatures {
   public:
     float ptrel;
+    float ptrel_noclip;
+    float erel;
+
+    float drsubjet1;
+    float drsubjet2;
 
     float puppiw;
     float deltaR;
+    float deltaR_noclip;
     float isGamma;
 
     float hadFrac;

DataFormats/BTauReco/interface/SecondaryVertexFeatures.h

@@ -6,6 +6,7 @@ namespace btagbtvdeep {
   class SecondaryVertexFeatures {
   public:
     float pt;
+    float ptrel;
     float mass;
 
     float deltaR;

DataFormats/BTauReco/src/classes_def.xml

@@ -409,16 +409,19 @@
     <version ClassVersion="3" checksum="386672900"/>
     <version ClassVersion="4" checksum="113477759"/>
   </class>
-  <class name="btagbtvdeep::SecondaryVertexFeatures" ClassVersion="3">
+  <class name="btagbtvdeep::SecondaryVertexFeatures" ClassVersion="4">
+   <version ClassVersion="4" checksum="3987593676"/>
     <version ClassVersion="3" checksum="3610780843"/>
   </class>
   <class name="btagbtvdeep::ShallowTagInfoFeatures" ClassVersion="3">
     <version ClassVersion="3" checksum="273372007"/>
   </class>
-  <class name="btagbtvdeep::NeutralCandidateFeatures" ClassVersion="3">
+  <class name="btagbtvdeep::NeutralCandidateFeatures" ClassVersion="4">
+   <version ClassVersion="4" checksum="3657612864"/>
     <version ClassVersion="3" checksum="938442516"/>
   </class>
-  <class name="btagbtvdeep::ChargedCandidateFeatures" ClassVersion="4">
+  <class name="btagbtvdeep::ChargedCandidateFeatures" ClassVersion="5">
+   <version ClassVersion="5" checksum="3385708923"/>
     <version ClassVersion="4" checksum="2173116472"/>
     <version ClassVersion="3" checksum="1302021141"/>
   </class>
@@ -432,7 +435,8 @@
     <version ClassVersion="3" checksum="2487956635"/>
     <version ClassVersion="4" checksum="3827695028"/>
   </class>
-    <class name="btagbtvdeep::DeepDoubleXFeatures" ClassVersion="3">
+    <class name="btagbtvdeep::DeepDoubleXFeatures" ClassVersion="4">
+     <version ClassVersion="4" checksum="3057931654"/>
     <version ClassVersion="3" checksum="1764960743"/>
   </class>
   <class name="btagbtvdeep::BoostedDoubleSVTagInfoFeatures" ClassVersion="3">

PhysicsTools/PatAlgos/python/recoLayer0/bTagging_cff.py

@@ -206,6 +206,12 @@
   , 'pfMassIndependentDeepDoubleCvLJetTags:probHcc'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
   , 'pfMassIndependentDeepDoubleCvBJetTags:probHbb'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
   , 'pfMassIndependentDeepDoubleCvBJetTags:probHcc'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
+  , 'pfMassIndependentDeepDoubleBvLV2JetTags:probQCD'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
+  , 'pfMassIndependentDeepDoubleBvLV2JetTags:probHbb'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
+  , 'pfMassIndependentDeepDoubleCvLV2JetTags:probQCD'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
+  , 'pfMassIndependentDeepDoubleCvLV2JetTags:probHcc'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
+  , 'pfMassIndependentDeepDoubleCvBV2JetTags:probHbb'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
+  , 'pfMassIndependentDeepDoubleCvBV2JetTags:probHcc'                     : [["pfDeepDoubleXTagInfos"], ['pfBoostedDoubleSVAK8TagInfos', "pfImpactParameterAK8TagInfos", 'pfInclusiveSecondaryVertexFinderAK8TagInfos']]
 }
 
 # meta-taggers are simple arithmetic on top of other taggers, they are stored here

RecoBTag/FeatureTools/interface/ChargedCandidateConverter.h

@@ -30,7 +30,12 @@ namespace btagbtvdeep {
       trackSip3dVal = -trackSip3dVal;
     }
 
+    c_pf_features.deltaR = reco::deltaR(*c_pf, jet);
     c_pf_features.ptrel = catch_infs_and_bound(c_pf->pt() / jet.pt(), 0, -1, 0, -1);
+    c_pf_features.ptrel_noclip = c_pf->pt() / jet.pt();
+    c_pf_features.erel = c_pf->energy() / jet.energy();
+    const float etasign = jet.eta() > 0 ? 1 : -1;
+    c_pf_features.etarel = etasign * (c_pf->eta() - jet.eta());
 
     c_pf_features.btagPf_trackEtaRel = catch_infs_and_bound(track_info.getTrackEtaRel(), 0, -5, 15);
     c_pf_features.btagPf_trackPtRel = catch_infs_and_bound(track_info.getTrackPtRel(), 0, -1, 4);
@@ -45,6 +50,26 @@ namespace btagbtvdeep {
     c_pf_features.btagPf_trackJetDistVal = catch_infs_and_bound(track_info.getTrackJetDistVal(), 0, -20, 1);
 
     c_pf_features.drminsv = catch_infs_and_bound(drminpfcandsv, 0, -1. * jetR, 0, -1. * jetR);
+
+    // subjet related
+    const auto* patJet = dynamic_cast<const pat::Jet*>(&jet);
+    if (!patJet) {
+      throw edm::Exception(edm::errors::InvalidReference) << "Input is not a pat::Jet.";
+    }
+
+    if (patJet->nSubjetCollections() > 0) {
+      auto subjets = patJet->subjets();
+      // sort by pt
+      std::sort(subjets.begin(), subjets.end(), [](const edm::Ptr<pat::Jet>& p1, const edm::Ptr<pat::Jet>& p2) {
+        return p1->pt() > p2->pt();
+      });
+      c_pf_features.drsubjet1 = !subjets.empty() ? reco::deltaR(*c_pf, *subjets.at(0)) : -1;
+      c_pf_features.drsubjet2 = subjets.size() > 1 ? reco::deltaR(*c_pf, *subjets.at(1)) : -1;
+    } else {
+      // AK4 jets don't have subjets
+      c_pf_features.drsubjet1 = -1;
+      c_pf_features.drsubjet2 = -1;
+    }
   }
 
   void packedCandidateToFeatures(const pat::PackedCandidate* c_pf,

RecoBTag/FeatureTools/interface/NeutralCandidateConverter.h

@@ -30,8 +30,31 @@ namespace btagbtvdeep {
                                         const float& drminpfcandsv,
                                         const float& jetR,
                                         NeutralCandidateFeatures& n_pf_features) {
+    const auto* patJet = dynamic_cast<const pat::Jet*>(&jet);
+
+    if (!patJet) {
+      throw edm::Exception(edm::errors::InvalidReference) << "Input is not a pat::Jet.";
+    }
+    // Do Subjets
+    if (patJet->nSubjetCollections() > 0) {
+      auto subjets = patJet->subjets();
+      // sort by pt
+      std::sort(subjets.begin(), subjets.end(), [](const edm::Ptr<pat::Jet>& p1, const edm::Ptr<pat::Jet>& p2) {
+        return p1->pt() > p2->pt();
+      });
+      n_pf_features.drsubjet1 = !subjets.empty() ? reco::deltaR(*n_pf, *subjets.at(0)) : -1;
+      n_pf_features.drsubjet2 = subjets.size() > 1 ? reco::deltaR(*n_pf, *subjets.at(1)) : -1;
+    } else {
+      n_pf_features.drsubjet1 = -1;
+      n_pf_features.drsubjet2 = -1;
+    }
+
+    // Jet relative vars
     n_pf_features.ptrel = catch_infs_and_bound(n_pf->pt() / jet.pt(), 0, -1, 0, -1);
+    n_pf_features.ptrel_noclip = n_pf->pt() / jet.pt();
     n_pf_features.deltaR = catch_infs_and_bound(reco::deltaR(*n_pf, jet), 0, -0.6, 0, -0.6);
+    n_pf_features.deltaR_noclip = reco::deltaR(*n_pf, jet);
+    n_pf_features.erel = n_pf->energy() / jet.energy();
     n_pf_features.isGamma = 0;
     if (std::abs(n_pf->pdgId()) == 22)
       n_pf_features.isGamma = 1;

RecoBTag/FeatureTools/plugins/DeepDoubleXTagInfoProducer.cc

@@ -22,6 +22,7 @@
 #include "DataFormats/BTauReco/interface/DeepDoubleXTagInfo.h"
 
 #include "RecoBTag/FeatureTools/interface/BoostedDoubleSVTagInfoConverter.h"
+#include "RecoBTag/FeatureTools/interface/NeutralCandidateConverter.h"
 #include "RecoBTag/FeatureTools/interface/ChargedCandidateConverter.h"
 #include "RecoBTag/FeatureTools/interface/JetConverter.h"
 #include "RecoBTag/FeatureTools/interface/SecondaryVertexConverter.h"
@@ -119,6 +120,8 @@ void DeepDoubleXTagInfoProducer::produce(edm::Event& iEvent, const edm::EventSet
 
     // reco jet reference (use as much as possible)
     const auto& jet = jets->at(jet_n);
+    const auto* pf_jet = dynamic_cast<const reco::PFJet*>(&jet);
+    const auto* pat_jet = dynamic_cast<const pat::Jet*>(&jet);
 
     edm::RefToBase<reco::Jet> jet_ref(jets, jet_n);
     if (jet.pt() > min_jet_pt_) {
@@ -175,7 +178,8 @@ void DeepDoubleXTagInfoProducer::produce(edm::Event& iEvent, const edm::EventSet
       math::XYZVector jet_dir = jet.momentum().Unit();
       GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
 
-      std::vector<btagbtvdeep::SortingClass<size_t>> c_sorted;
+      std::vector<btagbtvdeep::SortingClass<size_t>> c_sorted, n_sorted;
+      std::vector<int> n_indexes;
 
       // to cache the TrackInfo
       std::map<unsigned int, btagbtvdeep::TrackInfoBuilder> trackinfos;
@@ -205,9 +209,9 @@ void DeepDoubleXTagInfoProducer::produce(edm::Event& iEvent, const edm::EventSet
         }
       }
 
+      std::sort(daughters.begin(), daughters.end(), [](const auto& a, const auto& b) { return a->pt() > b->pt(); });
       for (unsigned int i = 0; i < daughters.size(); i++) {
         auto const* cand = daughters.at(i);
-
         if (cand) {
           // candidates under 950MeV (configurable) are not considered
           // might change if we use also white-listing
@@ -220,48 +224,61 @@ void DeepDoubleXTagInfoProducer::produce(edm::Event& iEvent, const edm::EventSet
                                   trackinfo.getTrackSip2dSig(),
                                   -btagbtvdeep::mindrsvpfcand(svs_unsorted, cand, jet_radius_),
                                   cand->pt() / jet.pt());
+          } else {
+            n_sorted.emplace_back(
+                i, -1, -btagbtvdeep::mindrsvpfcand(svs_unsorted, cand, jet_radius_), cand->pt() / jet.pt());
+            n_indexes.push_back(i);
           }
         }
       }
 
-      // sort collections (open the black-box if you please)
+      // sort collections in added order of priority
       std::sort(c_sorted.begin(), c_sorted.end(), btagbtvdeep::SortingClass<std::size_t>::compareByABCInv);
+      std::sort(n_sorted.begin(), n_sorted.end(), btagbtvdeep::SortingClass<std::size_t>::compareByABCInv);
 
-      std::vector<size_t> c_sortedindices;
+      std::vector<size_t> c_sortedindices, n_sortedindices;
 
       // this puts 0 everywhere and the right position in ind
       c_sortedindices = btagbtvdeep::invertSortingVector(c_sorted);
+      n_sortedindices = btagbtvdeep::invertSortingVector(n_sorted);
 
       // set right size to vectors
       features.c_pf_features.clear();
       features.c_pf_features.resize(c_sorted.size());
+      features.n_pf_features.clear();
+      features.n_pf_features.resize(n_sorted.size());
 
       for (unsigned int i = 0; i < daughters.size(); i++) {
         auto const* cand = daughters.at(i);
         if (cand) {
           // candidates under 950MeV are not considered
           // might change if we use also white-listing
-          if (cand->pt() < 0.95)
+          if (cand->pt() < min_candidate_pt_)
             continue;
+          auto packed_cand = dynamic_cast<const pat::PackedCandidate*>(cand);
+          auto reco_cand = dynamic_cast<const reco::PFCandidate*>(cand);
+
+          // need some edm::Ptr or edm::Ref if reco candidates
+          reco::PFCandidatePtr reco_ptr;
+          if (pf_jet) {
+            reco_ptr = pf_jet->getPFConstituent(i);
+          } else if (pat_jet && reco_cand) {
+            reco_ptr = pat_jet->getPFConstituent(i);
+          }
 
-          // get PUPPI weight from value map
           float puppiw = 1.0;  // fallback value
 
           float drminpfcandsv = btagbtvdeep::mindrsvpfcand(svs_unsorted, cand, jet_radius_);
-
           if (cand->charge() != 0) {
             // is charged candidate
             auto entry = c_sortedindices.at(i);
             // get cached track info
             auto& trackinfo = trackinfos.at(i);
             // get_ref to vector element
             auto& c_pf_features = features.c_pf_features.at(entry);
-            // fill feature structure
-            auto packed_cand = dynamic_cast<const pat::PackedCandidate*>(cand);
-            auto reco_cand = dynamic_cast<const reco::PFCandidate*>(cand);
             if (packed_cand) {
               btagbtvdeep::packedCandidateToFeatures(
-                  packed_cand, jet, trackinfo, drminpfcandsv, static_cast<float>(jet_radius_), c_pf_features);
+                  packed_cand, *pat_jet, trackinfo, drminpfcandsv, static_cast<float>(jet_radius_), c_pf_features);
             } else if (reco_cand) {
               // get vertex association quality
               int pv_ass_quality = 0;  // fallback value
@@ -288,6 +305,19 @@ void DeepDoubleXTagInfoProducer::produce(edm::Event& iEvent, const edm::EventSet
                                                    pv,
                                                    c_pf_features);
             }
+          } else {
+            // is neutral candidate
+            auto entry = n_sortedindices.at(i);
+            // get_ref to vector element
+            auto& n_pf_features = features.n_pf_features.at(entry);
+            // // fill feature structure
+            if (packed_cand) {
+              btagbtvdeep::packedCandidateToFeatures(
+                  packed_cand, *pat_jet, drminpfcandsv, static_cast<float>(jet_radius_), n_pf_features);
+            } else if (reco_cand) {
+              btagbtvdeep::recoCandidateToFeatures(
+                  reco_cand, jet, drminpfcandsv, static_cast<float>(jet_radius_), puppiw, n_pf_features);
+            }
           }
         }
       }

RecoBTag/FeatureTools/src/ChargedCandidateConverter.cc

@@ -27,6 +27,11 @@ namespace btagbtvdeep {
       c_pf_features.chi2 = catch_infs_and_bound(-1, 300, -1, 300);
       c_pf_features.quality = (1 << reco::TrackBase::loose);
     }
+
+    c_pf_features.dxy = catch_infs(c_pf->dxy());
+    c_pf_features.dz = catch_infs(c_pf->dz());
+    c_pf_features.dxysig = c_pf->bestTrack() ? catch_infs(c_pf->dxy() / c_pf->dxyError()) : 0;
+    c_pf_features.dzsig = c_pf->bestTrack() ? catch_infs(c_pf->dz() / c_pf->dzError()) : 0;
   }
 
   void recoCandidateToFeatures(const reco::PFCandidate* c_pf,
@@ -47,6 +52,16 @@ namespace btagbtvdeep {
     const auto& pseudo_track = (c_pf->bestTrack()) ? *c_pf->bestTrack() : reco::Track();
     c_pf_features.chi2 = catch_infs_and_bound(std::floor(pseudo_track.normalizedChi2()), 300, -1, 300);
     c_pf_features.quality = quality_from_pfcand(*c_pf);
+
+    // To be implemented if FatJet tag becomes RECO compatible
+    // const auto *trk =
+    // float dz =
+    // float dxy =
+
+    // c_pf_features.dxy =
+    // c_pf_features.dz =
+    // c_pf_features.dxysig =
+    // c_pf_features.dzsig =
   }
 
 }  // namespace btagbtvdeep

RecoBTag/FeatureTools/src/SecondaryVertexConverter.cc

@@ -16,6 +16,7 @@ namespace btagbtvdeep {
                     const bool flip) {
     math::XYZVector jet_dir = jet.momentum().Unit();
     sv_features.pt = sv.pt();
+    sv_features.ptrel = sv.pt() / jet.pt();
     sv_features.deltaR = catch_infs_and_bound(std::fabs(reco::deltaR(sv, jet_dir)) - 0.5, 0, -2, 0);
     sv_features.mass = sv.mass();
     sv_features.ntracks = sv.numberOfDaughters();