File indexing completed on 2023-03-17 11:17:10
0001 #ifndef RecoBTag_SecondaryVertex_CombinedSVComputer_h
0002 #define RecoBTag_SecondaryVertex_CombinedSVComputer_h
0003
0004 #include <iostream>
0005 #include <cstddef>
0006 #include <string>
0007 #include <cmath>
0008 #include <vector>
0009
0010 #include <Math/VectorUtil.h>
0011
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/Utilities/interface/Exception.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 #include "DataFormats/Math/interface/Vector3D.h"
0017 #include "DataFormats/Math/interface/LorentzVector.h"
0018 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0019 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0020 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0021 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0022 #include "DataFormats/TrackReco/interface/Track.h"
0023 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0024 #include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"
0025 #include "DataFormats/BTauReco/interface/CandSecondaryVertexTagInfo.h"
0026 #include "DataFormats/BTauReco/interface/TaggingVariable.h"
0027 #include "DataFormats/BTauReco/interface/VertexTypes.h"
0028 #include "DataFormats/BTauReco/interface/ParticleMasses.h"
0029 #include "DataFormats/JetReco/interface/PFJet.h"
0030 #include "DataFormats/PatCandidates/interface/Jet.h"
0031
0032 #include "RecoBTag/SecondaryVertex/interface/TrackSelector.h"
0033 #include "RecoBTag/SecondaryVertex/interface/V0Filter.h"
0034 #include "RecoBTag/SecondaryVertex/interface/TrackSorting.h"
0035 #include "RecoBTag/SecondaryVertex/interface/TrackKinematics.h"
0036
0037 #define range_for(i, x) for (int i = (x).begin; i != (x).end; i += (x).increment)
0038
0039 class CombinedSVComputer {
0040 public:
0041 explicit CombinedSVComputer(const edm::ParameterSet ¶ms);
0042 virtual ~CombinedSVComputer() = default;
0043 virtual reco::TaggingVariableList operator()(const reco::TrackIPTagInfo &ipInfo,
0044 const reco::SecondaryVertexTagInfo &svInfo) const;
0045 virtual reco::TaggingVariableList operator()(const reco::CandIPTagInfo &ipInfo,
0046 const reco::CandSecondaryVertexTagInfo &svInfo) const;
0047
0048 struct IterationRange {
0049 int begin, end, increment;
0050 };
0051 double flipValue(double value, bool vertex) const;
0052 IterationRange flipIterate(int size, bool vertex) const;
0053
0054 edm::ParameterSet dropDeltaR(const edm::ParameterSet &pset) const;
0055
0056 const reco::btag::TrackIPData &threshTrack(const reco::CandIPTagInfo &trackIPTagInfo,
0057 const reco::btag::SortCriteria sort,
0058 const reco::Jet &jet,
0059 const GlobalPoint &pv) const;
0060 const reco::btag::TrackIPData &threshTrack(const reco::TrackIPTagInfo &trackIPTagInfo,
0061 const reco::btag::SortCriteria sort,
0062 const reco::Jet &jet,
0063 const GlobalPoint &pv) const;
0064 template <class SVTI, class IPTI>
0065 void fillCommonVariables(reco::TaggingVariableList &vars,
0066 reco::TrackKinematics &vertexKinematics,
0067 const IPTI &ipInfo,
0068 const SVTI &svInfo,
0069 double &vtx_track_ptSum,
0070 double &vtx_track_ESum) const;
0071
0072 private:
0073 bool trackFlip;
0074 bool vertexFlip;
0075 double charmCut;
0076 reco::btag::SortCriteria sortCriterium;
0077 reco::TrackSelector trackSelector;
0078 reco::TrackSelector trackNoDeltaRSelector;
0079 reco::TrackSelector trackPseudoSelector;
0080 unsigned int pseudoMultiplicityMin;
0081 unsigned int trackMultiplicityMin;
0082 double minTrackWeight;
0083 bool useTrackWeights;
0084 bool vertexMassCorrection;
0085 reco::V0Filter pseudoVertexV0Filter;
0086 reco::V0Filter trackPairV0Filter;
0087 std::vector<reco::btau::TaggingVariableName> taggingVariables;
0088 };
0089
0090 template <class SVTI, class IPTI>
0091 void CombinedSVComputer::fillCommonVariables(reco::TaggingVariableList &vars,
0092 reco::TrackKinematics &vertexKinematics,
0093 const IPTI &ipInfo,
0094 const SVTI &svInfo,
0095 double &vtx_track_ptSum,
0096 double &vtx_track_ESum) const {
0097 using namespace ROOT::Math;
0098 using namespace reco;
0099
0100 typedef typename IPTI::input_container Container;
0101 typedef typename Container::value_type TrackRef;
0102
0103 edm::RefToBase<Jet> jet = ipInfo.jet();
0104 math::XYZVector jetDir = jet->momentum().Unit();
0105 bool havePv = ipInfo.primaryVertex().isNonnull();
0106 GlobalPoint pv;
0107 if (havePv)
0108 pv = GlobalPoint(ipInfo.primaryVertex()->x(), ipInfo.primaryVertex()->y(), ipInfo.primaryVertex()->z());
0109
0110 btag::Vertices::VertexType vtxType = btag::Vertices::NoVertex;
0111
0112 vars.insert(btau::jetPt, jet->pt(), true);
0113 vars.insert(btau::jetEta, jet->eta(), true);
0114 vars.insert(btau::jetAbsEta, fabs(jet->eta()), true);
0115
0116 if (ipInfo.selectedTracks().size() < trackMultiplicityMin)
0117 return;
0118
0119 vars.insert(btau::jetNTracks, ipInfo.selectedTracks().size(), true);
0120
0121 TrackKinematics allKinematics;
0122 TrackKinematics trackJetKinematics;
0123
0124 double jet_track_ESum = 0.;
0125
0126 int vtx = -1;
0127
0128 IterationRange range = flipIterate(svInfo.nVertices(), true);
0129 range_for(i, range) {
0130 if (vtx < 0)
0131 vtx = i;
0132 }
0133
0134 if (vtx >= 0) {
0135 vtxType = btag::Vertices::RecoVertex;
0136 vars.insert(btau::flightDistance1dVal, flipValue(svInfo.flightDistance(vtx, 1).value(), true), true);
0137 vars.insert(btau::flightDistance1dSig, flipValue(svInfo.flightDistance(vtx, 1).significance(), true), true);
0138 vars.insert(btau::flightDistance2dVal, flipValue(svInfo.flightDistance(vtx, 2).value(), true), true);
0139 vars.insert(btau::flightDistance2dSig, flipValue(svInfo.flightDistance(vtx, 2).significance(), true), true);
0140 vars.insert(btau::flightDistance3dVal, flipValue(svInfo.flightDistance(vtx, 3).value(), true), true);
0141 vars.insert(btau::flightDistance3dSig, flipValue(svInfo.flightDistance(vtx, 3).significance(), true), true);
0142 vars.insert(btau::vertexJetDeltaR, Geom::deltaR(svInfo.flightDirection(vtx), vertexFlip ? -jetDir : jetDir), true);
0143 vars.insert(btau::jetNSecondaryVertices, svInfo.nVertices(), true);
0144 }
0145
0146 std::vector<std::size_t> indices = ipInfo.sortedIndexes(sortCriterium);
0147 const std::vector<reco::btag::TrackIPData> &ipData = ipInfo.impactParameterData();
0148
0149 const Container &tracks = ipInfo.selectedTracks();
0150 std::vector<TrackRef> pseudoVertexTracks;
0151
0152 std::vector<TrackRef> trackPairV0Test(2);
0153 range = flipIterate(indices.size(), false);
0154 range_for(i, range) {
0155 std::size_t idx = indices[i];
0156 const reco::btag::TrackIPData &data = ipData[idx];
0157 const TrackRef &track = tracks[idx];
0158
0159 jet_track_ESum += std::sqrt(track->momentum().Mag2() + ROOT::Math::Square(ParticleMasses::piPlus));
0160
0161
0162
0163
0164
0165 if (!trackSelector(track, data, *jet, pv))
0166 continue;
0167
0168
0169 allKinematics.add(track);
0170
0171
0172 if (vtxType == btag::Vertices::NoVertex && trackPseudoSelector(track, data, *jet, pv)) {
0173 pseudoVertexTracks.push_back(track);
0174 vertexKinematics.add(track);
0175 }
0176
0177
0178 trackPairV0Test[0] = track;
0179 bool ok = true;
0180 range_for(j, range) {
0181 if (i == j)
0182 continue;
0183
0184 std::size_t pairIdx = indices[j];
0185 const reco::btag::TrackIPData &pairTrackData = ipData[pairIdx];
0186 const TrackRef &pairTrack = tracks[pairIdx];
0187
0188 if (!trackSelector(pairTrack, pairTrackData, *jet, pv))
0189 continue;
0190
0191 trackPairV0Test[1] = pairTrack;
0192 if (!trackPairV0Filter(trackPairV0Test)) {
0193 ok = false;
0194 break;
0195 }
0196 }
0197 if (!ok)
0198 continue;
0199
0200 trackJetKinematics.add(track);
0201
0202
0203 math::XYZVector trackMom = track->momentum();
0204 double trackMag = std::sqrt(trackMom.Mag2());
0205
0206 vars.insert(btau::trackSip3dVal, flipValue(data.ip3d.value(), false), true);
0207 vars.insert(btau::trackSip3dSig, flipValue(data.ip3d.significance(), false), true);
0208 vars.insert(btau::trackSip2dVal, flipValue(data.ip2d.value(), false), true);
0209 vars.insert(btau::trackSip2dSig, flipValue(data.ip2d.significance(), false), true);
0210 vars.insert(btau::trackJetDistVal, data.distanceToJetAxis.value(), true);
0211
0212
0213
0214
0215 vars.insert(btau::trackDecayLenVal, havePv ? (data.closestToJetAxis - pv).mag() : -1.0, true);
0216
0217 vars.insert(btau::trackMomentum, trackMag, true);
0218 vars.insert(btau::trackEta, trackMom.Eta(), true);
0219
0220
0221 double perp_trackMom_jetDir = VectorUtil::Perp2(trackMom, jetDir);
0222 perp_trackMom_jetDir = (perp_trackMom_jetDir > 0. ? std::sqrt(perp_trackMom_jetDir) : 0.);
0223
0224 vars.insert(btau::trackPtRel, perp_trackMom_jetDir, true);
0225 vars.insert(btau::trackPPar, jetDir.Dot(trackMom), true);
0226 vars.insert(btau::trackDeltaR, VectorUtil::DeltaR(trackMom, jetDir), true);
0227 vars.insert(btau::trackPtRatio, perp_trackMom_jetDir / trackMag, true);
0228 vars.insert(btau::trackPParRatio, jetDir.Dot(trackMom) / trackMag, true);
0229 }
0230
0231 if (vtxType == btag::Vertices::NoVertex && vertexKinematics.numberOfTracks() >= pseudoMultiplicityMin &&
0232 pseudoVertexV0Filter(pseudoVertexTracks)) {
0233 vtxType = btag::Vertices::PseudoVertex;
0234 for (typename std::vector<TrackRef>::const_iterator trkIt = pseudoVertexTracks.begin();
0235 trkIt != pseudoVertexTracks.end();
0236 ++trkIt) {
0237 vars.insert(btau::trackEtaRel, reco::btau::etaRel(jetDir, (*trkIt)->momentum()), true);
0238 vtx_track_ptSum += std::sqrt((*trkIt)->momentum().Perp2());
0239 vtx_track_ESum += std::sqrt((*trkIt)->momentum().Mag2() + ROOT::Math::Square(ParticleMasses::piPlus));
0240 }
0241 }
0242
0243 vars.insert(btau::vertexCategory, vtxType, true);
0244
0245 vars.insert(btau::trackJetPt, trackJetKinematics.vectorSum().Pt(), true);
0246 vars.insert(btau::trackSumJetDeltaR, VectorUtil::DeltaR(allKinematics.vectorSum(), jetDir), true);
0247 vars.insert(btau::trackSumJetEtRatio, allKinematics.vectorSum().Et() / ipInfo.jet()->et(), true);
0248
0249 vars.insert(btau::trackSip3dSigAboveCharm,
0250 flipValue(threshTrack(ipInfo, reco::btag::IP3DSig, *jet, pv).ip3d.significance(), false),
0251 true);
0252 vars.insert(btau::trackSip3dValAboveCharm,
0253 flipValue(threshTrack(ipInfo, reco::btag::IP3DSig, *jet, pv).ip3d.value(), false),
0254 true);
0255 vars.insert(btau::trackSip2dSigAboveCharm,
0256 flipValue(threshTrack(ipInfo, reco::btag::IP2DSig, *jet, pv).ip2d.significance(), false),
0257 true);
0258 vars.insert(btau::trackSip2dValAboveCharm,
0259 flipValue(threshTrack(ipInfo, reco::btag::IP2DSig, *jet, pv).ip2d.value(), false),
0260 true);
0261
0262 if (vtxType != btag::Vertices::NoVertex) {
0263 math::XYZTLorentzVector allSum = useTrackWeights ? allKinematics.weightedVectorSum() : allKinematics.vectorSum();
0264 math::XYZTLorentzVector vertexSum =
0265 useTrackWeights ? vertexKinematics.weightedVectorSum() : vertexKinematics.vectorSum();
0266
0267 if (vtxType != btag::Vertices::RecoVertex) {
0268 vars.insert(btau::vertexNTracks, vertexKinematics.numberOfTracks(), true);
0269 vars.insert(btau::vertexJetDeltaR, VectorUtil::DeltaR(vertexSum, jetDir), true);
0270 }
0271
0272 double vertexMass = vertexSum.M();
0273 if (vtxType == btag::Vertices::RecoVertex && vertexMassCorrection) {
0274 const GlobalVector &dir = svInfo.flightDirection(vtx);
0275 double vertexPt2 = math::XYZVector(dir.x(), dir.y(), dir.z()).Cross(vertexSum).Mag2() / dir.mag2();
0276 vertexMass = std::sqrt(vertexMass * vertexMass + vertexPt2) + std::sqrt(vertexPt2);
0277 }
0278 vars.insert(btau::vertexMass, vertexMass, true);
0279
0280 double varPi = (vertexMass / 5.2794) *
0281 (vtx_track_ESum / jet_track_ESum);
0282 vars.insert(btau::massVertexEnergyFraction, varPi / (varPi + 0.04), true);
0283 double varB = (std::sqrt(5.2794) * vtx_track_ptSum) / (vertexMass * std::sqrt(jet->pt()));
0284 vars.insert(btau::vertexBoostOverSqrtJetPt, varB * varB / (varB * varB + 10.), true);
0285
0286 if (allKinematics.numberOfTracks()) {
0287 vars.insert(btau::vertexEnergyRatio, vertexSum.E() / allSum.E(), true);
0288 } else {
0289 vars.insert(btau::vertexEnergyRatio, 1, true);
0290 }
0291 }
0292
0293 reco::PFJet const *pfJet = dynamic_cast<reco::PFJet const *>(&*jet);
0294 pat::Jet const *patJet = dynamic_cast<pat::Jet const *>(&*jet);
0295 if (pfJet != nullptr) {
0296 vars.insert(btau::chargedHadronEnergyFraction, pfJet->chargedHadronEnergyFraction(), true);
0297 vars.insert(btau::neutralHadronEnergyFraction, pfJet->neutralHadronEnergyFraction(), true);
0298 vars.insert(btau::photonEnergyFraction, pfJet->photonEnergyFraction(), true);
0299 vars.insert(btau::electronEnergyFraction, pfJet->electronEnergyFraction(), true);
0300 vars.insert(btau::muonEnergyFraction, pfJet->muonEnergyFraction(), true);
0301 vars.insert(btau::chargedHadronMultiplicity, pfJet->chargedHadronMultiplicity(), true);
0302 vars.insert(btau::neutralHadronMultiplicity, pfJet->neutralHadronMultiplicity(), true);
0303 vars.insert(btau::photonMultiplicity, pfJet->photonMultiplicity(), true);
0304 vars.insert(btau::electronMultiplicity, pfJet->electronMultiplicity(), true);
0305 vars.insert(btau::muonMultiplicity, pfJet->muonMultiplicity(), true);
0306 vars.insert(
0307 btau::hadronMultiplicity, pfJet->chargedHadronMultiplicity() + pfJet->neutralHadronMultiplicity(), true);
0308 vars.insert(btau::hadronPhotonMultiplicity,
0309 pfJet->chargedHadronMultiplicity() + pfJet->neutralHadronMultiplicity() + pfJet->photonMultiplicity(),
0310 true);
0311 vars.insert(btau::totalMultiplicity,
0312 pfJet->chargedHadronMultiplicity() + pfJet->neutralHadronMultiplicity() + pfJet->photonMultiplicity() +
0313 pfJet->electronMultiplicity() + pfJet->muonMultiplicity(),
0314 true);
0315
0316 } else if (patJet != nullptr && patJet->isPFJet()) {
0317 vars.insert(btau::chargedHadronEnergyFraction, patJet->chargedHadronEnergyFraction(), true);
0318 vars.insert(btau::neutralHadronEnergyFraction, patJet->neutralHadronEnergyFraction(), true);
0319 vars.insert(btau::photonEnergyFraction, patJet->photonEnergyFraction(), true);
0320 vars.insert(btau::electronEnergyFraction, patJet->electronEnergyFraction(), true);
0321 vars.insert(btau::muonEnergyFraction, patJet->muonEnergyFraction(), true);
0322 vars.insert(btau::chargedHadronMultiplicity, patJet->chargedHadronMultiplicity(), true);
0323 vars.insert(btau::neutralHadronMultiplicity, patJet->neutralHadronMultiplicity(), true);
0324 vars.insert(btau::photonMultiplicity, patJet->photonMultiplicity(), true);
0325 vars.insert(btau::electronMultiplicity, patJet->electronMultiplicity(), true);
0326 vars.insert(btau::muonMultiplicity, patJet->muonMultiplicity(), true);
0327 vars.insert(
0328 btau::hadronMultiplicity, patJet->chargedHadronMultiplicity() + patJet->neutralHadronMultiplicity(), true);
0329 vars.insert(
0330 btau::hadronPhotonMultiplicity,
0331 patJet->chargedHadronMultiplicity() + patJet->neutralHadronMultiplicity() + patJet->photonMultiplicity(),
0332 true);
0333 vars.insert(btau::totalMultiplicity,
0334 patJet->chargedHadronMultiplicity() + patJet->neutralHadronMultiplicity() +
0335 patJet->photonMultiplicity() + patJet->electronMultiplicity() + patJet->muonMultiplicity(),
0336 true);
0337
0338 } else {
0339 vars.insert(btau::chargedHadronEnergyFraction, 0., true);
0340 vars.insert(btau::neutralHadronEnergyFraction, 0., true);
0341 vars.insert(btau::photonEnergyFraction, 0., true);
0342 vars.insert(btau::electronEnergyFraction, 0., true);
0343 vars.insert(btau::muonEnergyFraction, 0., true);
0344 vars.insert(btau::chargedHadronMultiplicity, 0, true);
0345 vars.insert(btau::neutralHadronMultiplicity, 0, true);
0346 vars.insert(btau::photonMultiplicity, 0, true);
0347 vars.insert(btau::electronMultiplicity, 0, true);
0348 vars.insert(btau::muonMultiplicity, 0, true);
0349 vars.insert(btau::hadronMultiplicity, 0, true);
0350 vars.insert(btau::hadronPhotonMultiplicity, 0, true);
0351 vars.insert(btau::totalMultiplicity, 0, true);
0352 }
0353 }
0354
0355 #endif