BTauReco/interface/TrackProbabilityTagInfo.h

0001 #ifndef BTauReco_BJetTagTrackProbability_h
0002 #define BTauReco_BJetTagTrackProbability_h
0003
0004 #include "DataFormats/BTauReco/interface/RefMacros.h"
0005 #include "DataFormats/TrackReco/interface/Track.h"
0006 #include "DataFormats/BTauReco/interface/JTATagInfo.h"
0007 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0008
0009 namespace reco {
0010
0011   class TrackProbabilityTagInfo : public JTATagInfo {
0012   public:
0013     TrackProbabilityTagInfo(const std::vector<double>& probability2d,
0014                             const std::vector<double>& probability3d,
0015                             const std::vector<int>& trackOrder2d,
0016                             const std::vector<int>& trackOrder3d,
0017                             const JetTracksAssociationRef& jtaRef)
0018         : JTATagInfo(jtaRef),
0019           m_probability2d(probability2d),
0020           m_probability3d(probability3d),
0021           m_trackOrder2d(trackOrder2d),
0022           m_trackOrder3d(trackOrder3d) {}
0023
0024     TrackProbabilityTagInfo() {}
0025
0026     ~TrackProbabilityTagInfo() override {}
0027
0028     int factorial(int n) const {
0029       if (n < 2)
0030         return 1;
0031       else
0032         return n * factorial(n - 1);
0033     }
0034
0035     virtual float probability(size_t n, int ip) const {
0036       if (ip == 0) {
0037         if (n < m_probability3d.size())
0038           return m_probability3d[n];
0039       } else {
0040         if (n < m_probability2d.size())
0041           return m_probability2d[n];
0042       }
0043       return -10.;
0044     }
0045
0046     virtual float jetProbability(int ip, float minTrackProb) const {
0047       const std::vector<double>* vp;
0048       if (ip == 0)
0049         vp = &m_probability3d;
0050       else
0051         vp = &m_probability2d;
0052       const std::vector<double>& v = *vp;
0053
0054       int ngoodtracks = v.size();
0055       double SumJet = 0.;
0056
0057       for (std::vector<double>::const_iterator q = v.begin(); q != v.end(); q++) {
0058         SumJet += (*q > minTrackProb) ? log(*q) : log(minTrackProb);
0059       }
0060
0061       double ProbJet;
0062       double Loginvlog = 0;
0063
0064       if (SumJet < 0.) {
0065         if (ngoodtracks >= 2) {
0066           Loginvlog = log(-SumJet);
0067         }
0068         double Prob = 1.;
0069         for (int l = 1; l != ngoodtracks; l++) {
0070           Prob += exp(l * Loginvlog - log(1. * factorial(l)));
0071         }
0072         double LogProb = log(Prob);
0073         ProbJet = std::min(exp(std::max(LogProb + SumJet, -30.)), 1.);
0074       } else {
0075         ProbJet = 1.;
0076       }
0077       if (ProbJet > 1)
0078         std::cout << "ProbJet too high: " << ProbJet << std::endl;
0079
0080       //double LogProbJet=-log(ProbJet);
0081       //return 1.-ProbJet;
0082       return -log10(ProbJet) / 4.;
0083     }
0084
0085     /**
0086   Recompute discriminator
0087   ipType = 0 means 3d impact parameter
0088   ipType = 1 means transverse impact parameter
0089
0090   minProb is the minimum probability allowed  for a single track. Tracks with lower probability
0091   are considered with a probability = minProb.
0092  */
0093     virtual float discriminator(int ipType, float minProb) const { return jetProbability(ipType, minProb); }
0094
0095     virtual int selectedTracks(int ipType) const {
0096       if (ipType == 0)
0097         return m_probability3d.size();
0098       else
0099         return m_probability2d.size();
0100     }
0101     virtual int trackIndex(size_t n, int ip) const {
0102       if (ip == 0) {
0103         if (n < m_probability3d.size())
0104           return m_trackOrder3d[n];
0105       } else {
0106         if (n < m_probability2d.size())
0107           return m_trackOrder2d[n];
0108       }
0109       return 0;
0110     }
0111
0112     virtual const Track& track(size_t n, int ipType) const { return *tracks()[trackIndex(n, ipType)]; }
0113
0114     TrackProbabilityTagInfo* clone() const override { return new TrackProbabilityTagInfo(*this); }
0115
0116   private:
0117     std::vector<double> m_probability2d;  //
0118     std::vector<double> m_probability3d;  // create a smarter container instead of
0119     std::vector<int> m_trackOrder2d;      // this pair of vectors.
0120     std::vector<int> m_trackOrder3d;      //
0121   };
0122
0123   //typedef edm::ExtCollection< TrackProbabilityTagInfo,JetTagCollection> TrackProbabilityExtCollection;
0124   //typedef edm::OneToOneAssociation<JetTagCollection, TrackProbabilityTagInfo> TrackProbabilityExtCollection;
0125
0126   DECLARE_EDM_REFS(TrackProbabilityTagInfo)
0127
0128 }  // namespace reco
0129 #endif