ImpactParameter/interface/TemplatedJetBProbabilityComputer.h

0001 #ifndef ImpactParameter_TemplatedJetBProbabilityComputer_h
0002 #define ImpactParameter_TemplatedJetBProbabilityComputer_h
0003
0004 #include "DataFormats/TrackReco/interface/Track.h"
0005 #include "DataFormats/BTauReco/interface/TrackProbabilityTagInfo.h"
0006 #include "DataFormats/BTauReco/interface/TrackIPTagInfo.h"
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "Math/GenVector/VectorUtil.h"
0009 #include "RecoBTau/JetTagComputer/interface/JetTagComputer.h"
0010 #include <algorithm>
0011 #include <iostream>
0012
0013 template <class Container, class Base>
0014 class TemplatedJetBProbabilityComputer : public JetTagComputer {
0015 public:
0016   using Tokens = void;
0017
0018   typedef reco::IPTagInfo<Container, Base> TagInfo;
0019
0020   TemplatedJetBProbabilityComputer(const edm::ParameterSet& parameters) {
0021     m_ipType = parameters.getParameter<int>("impactParameterType");
0022     m_minTrackProb = parameters.getParameter<double>("minimumProbability");
0023     m_deltaR = parameters.getParameter<double>("deltaR");
0024     m_trackSign = parameters.getParameter<int>("trackIpSign");
0025     m_nbTracks = parameters.getParameter<unsigned int>("numberOfBTracks");
0026     m_cutMaxDecayLen = parameters.getParameter<double>("maximumDecayLength");
0027     m_cutMaxDistToAxis = parameters.getParameter<double>("maximumDistanceToJetAxis");
0028
0029     //
0030     // access track quality class; "any" takes everything
0031     //
0032     std::string trackQualityType = parameters.getParameter<std::string>("trackQualityClass");  //used
0033     m_trackQuality = reco::TrackBase::qualityByName(trackQualityType);
0034     m_useAllQualities = false;
0035     if (trackQualityType == "any" || trackQualityType == "Any" || trackQualityType == "ANY")
0036       m_useAllQualities = true;
0037
0038     useVariableJTA_ = parameters.getParameter<bool>("useVariableJTA");
0039     if (useVariableJTA_)
0040       varJTApars = {parameters.getParameter<double>("a_dR"),
0041                     parameters.getParameter<double>("b_dR"),
0042                     parameters.getParameter<double>("a_pT"),
0043                     parameters.getParameter<double>("b_pT"),
0044                     parameters.getParameter<double>("min_pT"),
0045                     parameters.getParameter<double>("max_pT"),
0046                     parameters.getParameter<double>("min_pT_dRcut"),
0047                     parameters.getParameter<double>("max_pT_dRcut"),
0048                     parameters.getParameter<double>("max_pT_trackPTcut")};
0049
0050     uses("ipTagInfos");
0051   }
0052
0053   float discriminator(const TagInfoHelper& ti) const override {
0054     const TagInfo& tkip = ti.get<TagInfo>();
0055     const Container& tracks(tkip.selectedTracks());
0056     const std::vector<float>& allProbabilities((tkip.probabilities(m_ipType)));
0057     const std::vector<reco::btag::TrackIPData>& impactParameters((tkip.impactParameterData()));
0058
0059     if (tkip.primaryVertex().isNull())
0060       return 0;
0061
0062     GlobalPoint pv(tkip.primaryVertex()->position().x(),
0063                    tkip.primaryVertex()->position().y(),
0064                    tkip.primaryVertex()->position().z());
0065
0066     std::vector<float> probabilities;
0067     std::vector<float> probabilitiesB;
0068     int i = 0;
0069     for (std::vector<float>::const_iterator it = allProbabilities.begin(); it != allProbabilities.end(); ++it, i++) {
0070       if (fabs(impactParameters[i].distanceToJetAxis.value()) < m_cutMaxDistToAxis &&  // distance to JetAxis
0071           (impactParameters[i].closestToJetAxis - pv).mag() < m_cutMaxDecayLen &&      // max decay len
0072           (m_useAllQualities == true ||
0073            reco::btag::toTrack(tracks[i])->quality(m_trackQuality))  // use selected track qualities
0074       ) {
0075         // Use only positive(or negative) tracks for B
0076         float p = fabs(*it);
0077         if (useVariableJTA_) {
0078           if (tkip.variableJTA(varJTApars)[i]) {
0079             if (m_trackSign > 0 || *it < 0)
0080               probabilities.push_back(p);  //Use all tracks for positive tagger and only negative for negative tagger
0081             if (m_trackSign > 0 && *it >= 0) {
0082               probabilitiesB.push_back(*it);
0083             }  //Use only positive tracks for positive tagger
0084             if (m_trackSign < 0 && *it <= 0) {
0085               probabilitiesB.push_back(-*it);
0086             }  //Use only negative tracks for negative tagger
0087           }
0088         } else if (m_deltaR < 0 ||
0089                    ROOT::Math::VectorUtil::DeltaR((*tkip.jet()).p4().Vect(), (*tracks[i]).momentum()) < m_deltaR) {
0090           //if(m_trackSign>0 || *it >0 ) probabilities.push_back(p); //Use all tracks for positive tagger and only negative for negative tagger
0091           if (m_trackSign > 0 || *it < 0)
0092             probabilities.push_back(p);  //Use all tracks for positive tagger and only negative for negative tagger
0093
0094           if (m_trackSign > 0 && *it >= 0) {
0095             probabilitiesB.push_back(*it);
0096           }  //Use only positive tracks for positive tagger
0097           if (m_trackSign < 0 && *it <= 0) {
0098             probabilitiesB.push_back(-*it);
0099           }  //Use only negative tracks for negative tagger
0100         }
0101       }
0102     }
0103
0104     float all = jetProbability(probabilities);
0105     std::sort(probabilitiesB.begin(), probabilitiesB.end());
0106     if (probabilitiesB.size() > m_nbTracks)
0107       probabilitiesB.resize(m_nbTracks);
0108     float b = jetProbability(probabilitiesB);
0109
0110     return -log(b) / 4 - log(all) / 4;  ///log(all);
0111   }
0112
0113   double jetProbability(const std::vector<float>& v) const {
0114     int ngoodtracks = v.size();
0115     double SumJet = 0.;
0116
0117     for (std::vector<float>::const_iterator q = v.begin(); q != v.end(); q++) {
0118       SumJet += (*q > m_minTrackProb) ? log(*q) : log(m_minTrackProb);
0119     }
0120
0121     double ProbJet;
0122     double Loginvlog = 0;
0123
0124     if (SumJet < 0.) {
0125       if (ngoodtracks >= 2) {
0126         Loginvlog = log(-SumJet);
0127       }
0128       double Prob = 1.;
0129       double lfact = 1.;
0130       for (int l = 1; l != ngoodtracks; l++) {
0131         lfact *= l;
0132         Prob += exp(l * Loginvlog - log(1. * lfact));
0133       }
0134       double LogProb = log(Prob);
0135       ProbJet = std::min(exp(std::max(LogProb + SumJet, -30.)), 1.);
0136     } else {
0137       ProbJet = 1.;
0138     }
0139     if (ProbJet > 1)
0140       std::cout << "ProbJet too high: " << ProbJet << std::endl;
0141
0142     //double LogProbJet=-log(ProbJet);
0143     //  //return 1.-ProbJet;
0144     return ProbJet;
0145   }
0146
0147 private:
0148   bool useVariableJTA_;
0149   reco::btag::variableJTAParameters varJTApars;
0150   double m_minTrackProb;
0151   int m_ipType;
0152   double m_deltaR;
0153   int m_trackSign;
0154   unsigned int m_nbTracks;
0155   double m_cutMaxDecayLen;
0156   double m_cutMaxDistToAxis;
0157   reco::TrackBase::TrackQuality m_trackQuality;
0158   bool m_useAllQualities;
0159 };
0160
0161 #endif  // ImpactParameter_TemplatedJetBProbabilityComputer_h