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   static void fillPSetDescription(edm::ParameterSetDescription& desc) {
0054     desc.add<int>("impactParameterType", 0)->setComment("0 = 3D, 1 = 2D");
0055     desc.add<double>("minimumProbability", 0.005);
0056     desc.add<double>("deltaR", -1.0);
0057     desc.add<int>("trackIpSign", 1)->setComment("0 = use both, 1 = positive only, -1 = negative only");
0058     desc.add<unsigned int>("numberOfBTracks", 4);
0059     desc.add<double>("maximumDecayLength", 5.0);
0060     desc.add<double>("maximumDistanceToJetAxis", 0.07);
0061     desc.add<std::string>("trackQualityClass", "any");
0062     desc.add<bool>("useVariableJTA", false);
0063     desc.add<double>("a_dR", -0.001053);
0064     desc.add<double>("b_dR", 0.6263);
0065     desc.add<double>("a_pT", 0.005263);
0066     desc.add<double>("b_pT", 0.3684);
0067     desc.add<double>("min_pT", 120);
0068     desc.add<double>("max_pT", 500);
0069     desc.add<double>("min_pT_dRcut", 0.5);
0070     desc.add<double>("max_pT_dRcut", 0.1);
0071     desc.add<double>("max_pT_trackPTcut", 3);
0072   }
0073
0074   float discriminator(const TagInfoHelper& ti) const override {
0075     const TagInfo& tkip = ti.get<TagInfo>();
0076     const Container& tracks(tkip.selectedTracks());
0077     const std::vector<float>& allProbabilities((tkip.probabilities(m_ipType)));
0078     const std::vector<reco::btag::TrackIPData>& impactParameters((tkip.impactParameterData()));
0079
0080     if (tkip.primaryVertex().isNull())
0081       return 0;
0082
0083     GlobalPoint pv(tkip.primaryVertex()->position().x(),
0084                    tkip.primaryVertex()->position().y(),
0085                    tkip.primaryVertex()->position().z());
0086
0087     std::vector<float> probabilities;
0088     std::vector<float> probabilitiesB;
0089     int i = 0;
0090     for (std::vector<float>::const_iterator it = allProbabilities.begin(); it != allProbabilities.end(); ++it, i++) {
0091       if (fabs(impactParameters[i].distanceToJetAxis.value()) < m_cutMaxDistToAxis &&  // distance to JetAxis
0092           (impactParameters[i].closestToJetAxis - pv).mag() < m_cutMaxDecayLen &&      // max decay len
0093           (m_useAllQualities == true ||
0094            reco::btag::toTrack(tracks[i])->quality(m_trackQuality))  // use selected track qualities
0095       ) {
0096         // Use only positive(or negative) tracks for B
0097         float p = fabs(*it);
0098         if (useVariableJTA_) {
0099           if (tkip.variableJTA(varJTApars)[i]) {
0100             if (m_trackSign > 0 || *it < 0)
0101               probabilities.push_back(p);  //Use all tracks for positive tagger and only negative for negative tagger
0102             if (m_trackSign > 0 && *it >= 0) {
0103               probabilitiesB.push_back(*it);
0104             }  //Use only positive tracks for positive tagger
0105             if (m_trackSign < 0 && *it <= 0) {
0106               probabilitiesB.push_back(-*it);
0107             }  //Use only negative tracks for negative tagger
0108           }
0109         } else if (m_deltaR < 0 ||
0110                    ROOT::Math::VectorUtil::DeltaR((*tkip.jet()).p4().Vect(), (*tracks[i]).momentum()) < m_deltaR) {
0111           //if(m_trackSign>0 || *it >0 ) probabilities.push_back(p); //Use all tracks for positive tagger and only negative for negative tagger
0112           if (m_trackSign > 0 || *it < 0)
0113             probabilities.push_back(p);  //Use all tracks for positive tagger and only negative for negative tagger
0114
0115           if (m_trackSign > 0 && *it >= 0) {
0116             probabilitiesB.push_back(*it);
0117           }  //Use only positive tracks for positive tagger
0118           if (m_trackSign < 0 && *it <= 0) {
0119             probabilitiesB.push_back(-*it);
0120           }  //Use only negative tracks for negative tagger
0121         }
0122       }
0123     }
0124
0125     float all = jetProbability(probabilities);
0126     std::sort(probabilitiesB.begin(), probabilitiesB.end());
0127     if (probabilitiesB.size() > m_nbTracks)
0128       probabilitiesB.resize(m_nbTracks);
0129     float b = jetProbability(probabilitiesB);
0130
0131     return -log(b) / 4 - log(all) / 4;  ///log(all);
0132   }
0133
0134   double jetProbability(const std::vector<float>& v) const {
0135     int ngoodtracks = v.size();
0136     double SumJet = 0.;
0137
0138     for (std::vector<float>::const_iterator q = v.begin(); q != v.end(); q++) {
0139       SumJet += (*q > m_minTrackProb) ? log(*q) : log(m_minTrackProb);
0140     }
0141
0142     double ProbJet;
0143     double Loginvlog = 0;
0144
0145     if (SumJet < 0.) {
0146       if (ngoodtracks >= 2) {
0147         Loginvlog = log(-SumJet);
0148       }
0149       double Prob = 1.;
0150       double lfact = 1.;
0151       for (int l = 1; l != ngoodtracks; l++) {
0152         lfact *= l;
0153         Prob += exp(l * Loginvlog - log(1. * lfact));
0154       }
0155       double LogProb = log(Prob);
0156       ProbJet = std::min(exp(std::max(LogProb + SumJet, -30.)), 1.);
0157     } else {
0158       ProbJet = 1.;
0159     }
0160     if (ProbJet > 1)
0161       std::cout << "ProbJet too high: " << ProbJet << std::endl;
0162
0163     //double LogProbJet=-log(ProbJet);
0164     //  //return 1.-ProbJet;
0165     return ProbJet;
0166   }
0167
0168 private:
0169   bool useVariableJTA_;
0170   reco::btag::variableJTAParameters varJTApars;
0171   double m_minTrackProb;
0172   int m_ipType;
0173   double m_deltaR;
0174   int m_trackSign;
0175   unsigned int m_nbTracks;
0176   double m_cutMaxDecayLen;
0177   double m_cutMaxDistToAxis;
0178   reco::TrackBase::TrackQuality m_trackQuality;
0179   bool m_useAllQualities;
0180 };
0181
0182 #endif  // ImpactParameter_TemplatedJetBProbabilityComputer_h