Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-08 03:36:26 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!