ImpactParameter/interface/TemplatedJetProbabilityComputer.h

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