BTauReco/interface/IPTagInfo.h

0001 #ifndef DataFormats_BTauReco_IpTagInfo_h
0002 #define DataFormats_BTauReco_IpTagInfo_h
0003
0004 #include "DataFormats/Common/interface/CMS_CLASS_VERSION.h"
0005 #include "DataFormats/BTauReco/interface/RefMacros.h"
0006 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0007 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0008 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0009 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0010 #include "DataFormats/BTauReco/interface/JTATagInfo.h"
0011 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0012 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0013 #include <cmath>
0014 #include <map>
0015 #include <Math/VectorUtil.h>
0016 #include "DataFormats/Math/interface/Vector3D.h"
0017 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019 #include "DataFormats/VertexReco/interface/Vertex.h"
0020
0021 namespace reco {
0022   namespace btag {
0023
0024     inline const reco::Track* toTrack(const reco::TrackBaseRef& t) { return &(*t); }
0025     inline const reco::Track* toTrack(const reco::TrackRef& t) { return &(*t); }
0026     inline const reco::Track* toTrack(const reco::CandidatePtr& c) { return (*c).bestTrack(); }
0027
0028     struct TrackIPData {
0029       GlobalPoint closestToJetAxis;
0030       GlobalPoint closestToGhostTrack;
0031       Measurement1D ip2d;
0032       Measurement1D ip3d;
0033       Measurement1D distanceToJetAxis;
0034       Measurement1D distanceToGhostTrack;
0035       float ghostTrackWeight;
0036     };
0037     struct variableJTAParameters {
0038       double a_dR, b_dR, a_pT, b_pT;
0039       double min_pT, max_pT;
0040       double min_pT_dRcut, max_pT_dRcut;
0041       double max_pT_trackPTcut;
0042     };
0043     enum SortCriteria { IP3DSig = 0, Prob3D, IP2DSig, Prob2D, IP3DValue, IP2DValue };
0044
0045   }  // namespace btag
0046
0047   template <class Container, class Base>
0048   class IPTagInfo : public Base {
0049   public:
0050     typedef Container input_container;
0051     typedef Base base_class;
0052
0053     IPTagInfo(const std::vector<btag::TrackIPData>& ipData,
0054               const std::vector<float>& prob2d,
0055               const std::vector<float>& prob3d,
0056               const Container& selected,
0057               const Base& base,
0058               const edm::Ref<VertexCollection>& pv,
0059               const GlobalVector& axis,
0060               const TrackRef& ghostTrack)
0061         : Base(base),
0062           m_data(ipData),
0063           m_prob2d(prob2d),
0064           m_prob3d(prob3d),
0065           m_selected(selected),
0066           m_pv(pv),
0067           m_axis(axis),
0068           m_ghostTrack(ghostTrack) {}
0069
0070     IPTagInfo() {}
0071
0072     ~IPTagInfo() override {}
0073
0074     /// clone
0075     IPTagInfo* clone(void) const override { return new IPTagInfo(*this); }
0076
0077     /**
0078    Check if probability information is globally available
0079    impact parameters in the collection
0080
0081    Even if true for some tracks it is possible that a -1 probability is returned
0082    if some problem occured
0083   */
0084
0085     virtual bool hasProbabilities() const { return m_data.size() == m_prob3d.size(); }
0086
0087     /**
0088    Vectors of TrackIPData orderd as the selectedTracks()
0089    */
0090     const std::vector<btag::TrackIPData>& impactParameterData() const { return m_data; }
0091
0092     /**
0093    Return the vector of tracks for which the IP information is available
0094    Quality cuts are applied to reject fake tracks
0095   */
0096     const Container& selected() const { return m_selected; }
0097
0098     //legacy name for compatibility
0099     const Container& selectedTracks() const { return m_selected; }
0100
0101     const std::vector<float>& probabilities(int ip) const { return (ip == 0) ? m_prob3d : m_prob2d; }
0102
0103     /**
0104    Return the list of track index sorted by mode
0105    A cut can is specified to select only tracks with
0106    IP value or significance > cut
0107    or
0108    probability < cut
0109    (according to the specified mode)
0110   */
0111     std::vector<size_t> sortedIndexesWithCut(float cut, btag::SortCriteria mode = reco::btag::IP3DSig) const;
0112
0113     /**
0114    variable jet-to track association:
0115    returns vector of bool, indicating for each track whether it passed
0116    the variable JTA.
0117   */
0118     std::vector<bool> variableJTA(const btag::variableJTAParameters& params) const;
0119     static bool passVariableJTA(const btag::variableJTAParameters& params,
0120                                 double jetpt,
0121                                 double trackpt,
0122                                 double jettrackdr);
0123
0124     /**
0125    Return the list of track index sorted by mode
0126   */
0127     std::vector<size_t> sortedIndexes(btag::SortCriteria mode = reco::btag::IP3DSig) const;
0128     Container sorted(const std::vector<size_t>& indexes) const;
0129     Container sortedTracks(const std::vector<size_t>& indexes) const { return sorted(indexes); }
0130
0131     TaggingVariableList taggingVariables(void) const override;
0132
0133     const edm::Ref<VertexCollection>& primaryVertex() const { return m_pv; }
0134
0135     const GlobalVector& axis() const { return m_axis; }
0136     const TrackRef& ghostTrack() const { return m_ghostTrack; }
0137
0138     const Track* selectedTrack(size_t i) const { return reco::btag::toTrack(m_selected[i]); }
0139
0140     // Used by ROOT storage
0141     CMS_CLASS_VERSION(11)
0142
0143   private:
0144     std::vector<btag::TrackIPData> m_data;
0145     std::vector<float> m_prob2d;
0146     std::vector<float> m_prob3d;
0147     Container m_selected;
0148     edm::Ref<VertexCollection> m_pv;
0149     GlobalVector m_axis;
0150     TrackRef m_ghostTrack;
0151   };
0152
0153   //Explicit templates:
0154   //template <> const Track *IPTagInfo<TrackRefVector,JTATagInfo>::selectedTrack(size_t i) const {return &(*m_selected[i]);}
0155   //template <> const Track *IPTagInfo<std::vector<CandidatePtr>,BaseTagInfo>::selectedTrack(size_t i) const {return (*m_selected[i]).bestTrack();}
0156
0157   template <class Container, class Base>
0158   TaggingVariableList IPTagInfo<Container, Base>::taggingVariables(void) const {
0159     TaggingVariableList vars;
0160
0161     math::XYZVector jetDir = Base::jet()->momentum().Unit();
0162     bool havePv = primaryVertex().isNonnull();
0163     GlobalPoint pv;
0164     if (havePv)
0165       pv = GlobalPoint(primaryVertex()->x(), primaryVertex()->y(), primaryVertex()->z());
0166
0167     std::vector<size_t> indexes = sortedIndexes();  // use default criterium
0168     for (std::vector<size_t>::const_iterator it = indexes.begin(); it != indexes.end(); ++it) {
0169       using namespace ROOT::Math;
0170       const Track* track = selectedTrack(*it);
0171       const btag::TrackIPData* data = &m_data[*it];
0172       const math::XYZVector& trackMom = track->momentum();
0173       double trackMag = std::sqrt(trackMom.Mag2());
0174
0175       vars.insert(btau::trackMomentum, trackMag, true);
0176       vars.insert(btau::trackEta, trackMom.Eta(), true);
0177       vars.insert(btau::trackEtaRel, reco::btau::etaRel(jetDir, trackMom), true);
0178       vars.insert(btau::trackPtRel, VectorUtil::Perp(trackMom, jetDir), true);
0179       vars.insert(btau::trackPPar, jetDir.Dot(trackMom), true);
0180       vars.insert(btau::trackDeltaR, VectorUtil::DeltaR(trackMom, jetDir), true);
0181       vars.insert(btau::trackPtRatio, VectorUtil::Perp(trackMom, jetDir) / trackMag, true);
0182       vars.insert(btau::trackPParRatio, jetDir.Dot(trackMom) / trackMag, true);
0183       vars.insert(btau::trackSip3dVal, data->ip3d.value(), true);
0184       vars.insert(btau::trackSip3dSig, data->ip3d.significance(), true);
0185       vars.insert(btau::trackSip2dVal, data->ip2d.value(), true);
0186       vars.insert(btau::trackSip2dSig, data->ip2d.significance(), true);
0187       vars.insert(btau::trackDecayLenVal, havePv ? (data->closestToJetAxis - pv).mag() : -1.0, true);
0188       vars.insert(btau::trackJetDistVal, data->distanceToJetAxis.value(), true);
0189       vars.insert(btau::trackJetDistSig, data->distanceToJetAxis.significance(), true);
0190       vars.insert(btau::trackGhostTrackDistVal, data->distanceToGhostTrack.value(), true);
0191       vars.insert(btau::trackGhostTrackDistSig, data->distanceToGhostTrack.significance(), true);
0192       vars.insert(btau::trackGhostTrackWeight, data->ghostTrackWeight, true);
0193       vars.insert(btau::trackChi2, track->normalizedChi2(), true);
0194       vars.insert(btau::trackNTotalHits, track->hitPattern().numberOfValidHits(), true);
0195       vars.insert(btau::trackNPixelHits, track->hitPattern().numberOfValidPixelHits(), true);
0196     }
0197     vars.finalize();
0198     return vars;
0199   }
0200
0201   template <class Container, class Base>
0202   Container IPTagInfo<Container, Base>::sorted(const std::vector<size_t>& indexes) const {
0203     Container tr;
0204     for (size_t i = 0; i < indexes.size(); i++)
0205       tr.push_back(m_selected[indexes[i]]);
0206     return tr;
0207   }
0208
0209   template <class Container, class Base>
0210   std::vector<bool> IPTagInfo<Container, Base>::variableJTA(const btag::variableJTAParameters& params) const {
0211     std::vector<bool> result;
0212
0213     //Jet parameters
0214     double jetpT = Base::jet()->pt();
0215     math::XYZVector jetDir = Base::jet()->momentum().Unit();
0216
0217     for (size_t i = 0; i < m_selected.size(); i++) {
0218       //Track parameters
0219       const Track* track = selectedTrack(i);
0220       double trackpT = track->pt();
0221       const math::XYZVector& trackMom = track->momentum();
0222
0223       // do the math in passVariableJTA
0224       result.push_back(passVariableJTA(params, jetpT, trackpT, ROOT::Math::VectorUtil::DeltaR(trackMom, jetDir)));
0225     }
0226
0227     return result;
0228   }
0229
0230   template <class Container, class Base>
0231   std::vector<size_t> IPTagInfo<Container, Base>::sortedIndexes(btag::SortCriteria mode) const {
0232     using namespace reco::btag;
0233     float cut = -1e99;
0234     if ((mode == Prob3D || mode == Prob2D))
0235       cut = 1e99;
0236     return sortedIndexesWithCut(cut, mode);
0237   }
0238
0239   template <class Container, class Base>
0240   std::vector<size_t> IPTagInfo<Container, Base>::sortedIndexesWithCut(float cut, btag::SortCriteria mode) const {
0241     std::multimap<float, size_t> sortedIdx;
0242     size_t nSelectedTracks = m_selected.size();
0243     std::vector<size_t> result;
0244     using namespace reco::btag;
0245
0246     //check if probabilities are available
0247     if ((mode == Prob3D || mode == Prob2D) && !hasProbabilities()) {
0248       return result;
0249     }
0250
0251     for (size_t i = 0; i < nSelectedTracks; i++) {
0252       float sortingKey;
0253       switch (mode) {
0254         case IP3DSig:
0255           sortingKey = m_data[i].ip3d.significance();
0256           break;
0257         case IP2DSig:
0258           sortingKey = m_data[i].ip2d.significance();
0259           break;
0260         case IP3DValue:
0261           sortingKey = m_data[i].ip3d.value();
0262           break;
0263         case IP2DValue:
0264           sortingKey = m_data[i].ip2d.value();
0265           break;
0266         case Prob3D:
0267           sortingKey = m_prob3d[i];
0268           break;
0269         case Prob2D:
0270           sortingKey = m_prob2d[i];
0271           break;
0272
0273         default:
0274           sortingKey = i;
0275       }
0276       sortedIdx.insert(std::pair<float, size_t>(sortingKey, i));
0277     }
0278
0279     //Descending:
0280     if (mode == IP3DSig || mode == IP2DSig || mode == IP3DValue || mode == IP2DValue) {
0281       for (std::multimap<float, size_t>::reverse_iterator it = sortedIdx.rbegin(); it != sortedIdx.rend(); it++)
0282         if (it->first >= cut)
0283           result.push_back(it->second);
0284     } else
0285     //Ascending:
0286     {
0287       for (std::multimap<float, size_t>::iterator it = sortedIdx.begin(); it != sortedIdx.end(); it++)
0288         if (it->first <= cut)
0289           result.push_back(it->second);
0290     }
0291     return result;
0292   }
0293
0294   template <class Container, class Base>
0295   bool IPTagInfo<Container, Base>::passVariableJTA(const btag::variableJTAParameters& params,
0296                                                    double jetpT,
0297                                                    double trackpT,
0298                                                    double jettrackdR) {
0299     bool pass = false;
0300
0301     // intermediate pt range (between min_pT and max_pT), apply variable JTA !
0302     if (jetpT > params.min_pT && jetpT < params.max_pT) {
0303       double deltaRfunction_highpt = -jetpT * params.a_dR + params.b_dR;
0304       double ptfunction_highpt = jetpT * params.a_pT + params.b_pT;
0305
0306       if (jettrackdR < deltaRfunction_highpt && trackpT > ptfunction_highpt)
0307         pass = true;
0308
0309       //  cout << "IPTagInfo: passVariableJTA: dR and TrackpT " << jettrackdR << " " << trackpT << endl;
0310
0311       //high pt range, apply fixed default cuts
0312     } else if (jetpT > params.max_pT) {
0313       if (jettrackdR < params.max_pT_dRcut && trackpT > params.max_pT_trackPTcut)
0314         pass = true;
0315
0316       // low pt range, apply fixed default cuts
0317     } else {
0318       if (jettrackdR < params.min_pT_dRcut)
0319         pass = true;
0320     }
0321
0322     return pass;
0323   }
0324 }  // namespace reco
0325 #endif