ImpactParameter/interface/TemplatedTrackCountingComputer.h

0001 #ifndef ImpactParameter_TemplatedTrackCountingComputer_h
0002 #define ImpactParameter_TemplatedTrackCountingComputer_h
0003
0004 #include "DataFormats/TrackReco/interface/Track.h"
0005 #include "DataFormats/VertexReco/interface/Vertex.h"
0006 #include "DataFormats/BTauReco/interface/TrackCountingTagInfo.h"
0007 #include "DataFormats/BTauReco/interface/IPTagInfo.h"
0008 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0009 #include "Math/GenVector/VectorUtil.h"
0010 #include "RecoBTau/JetTagComputer/interface/JetTagComputer.h"
0011
0012 template <class Container, class Base>
0013 class TemplatedTrackCountingComputer : public JetTagComputer {
0014 public:
0015   using Tokens = void;
0016
0017   typedef reco::IPTagInfo<Container, Base> TagInfo;
0018
0019   TemplatedTrackCountingComputer(const edm::ParameterSet& parameters) {
0020     m_minIP = parameters.getParameter<double>("minimumImpactParameter");
0021     m_useSignedIPSig = parameters.getParameter<bool>("useSignedImpactParameterSig");
0022     m_nthTrack = parameters.getParameter<int>("nthTrack");
0023     m_ipType = parameters.getParameter<int>("impactParameterType");
0024     m_deltaR = parameters.getParameter<double>("deltaR");
0025     m_cutMaxDecayLen = parameters.getParameter<double>("maximumDecayLength");          //used
0026     m_cutMaxDistToAxis = parameters.getParameter<double>("maximumDistanceToJetAxis");  //used
0027     //
0028     // access track quality class; "any" takes everything
0029     //
0030     std::string trackQualityType = parameters.getParameter<std::string>("trackQualityClass");  //used
0031     m_trackQuality = reco::TrackBase::qualityByName(trackQualityType);
0032     m_useAllQualities = false;
0033     if (trackQualityType == "any" || trackQualityType == "Any" || trackQualityType == "ANY")
0034       m_useAllQualities = true;
0035
0036     uses("ipTagInfos");
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   }
0051
0052   float discriminator(const TagInfoHelper& ti) const override {
0053     const TagInfo& tkip = ti.get<TagInfo>();
0054     std::multiset<float> significances = orderedSignificances(tkip);
0055     std::multiset<float>::reverse_iterator nth = significances.rbegin();
0056     for (int i = 0; i < m_nthTrack - 1 && nth != significances.rend(); i++)
0057       nth++;
0058     if (nth != significances.rend())
0059       return *nth;
0060     else
0061       return -100.;
0062   }
0063
0064   static void fillPSetDescription(edm::ParameterSetDescription& desc) {
0065     desc.add<double>("minimumImpactParameter", -1.);
0066     desc.add<bool>("useSignedImpactParameterSig", true);
0067     desc.add<int>("nthTrack", -1);
0068     desc.add<int>("impactParameterType", 0)->setComment("0 = 3D, 1 = 2D");
0069     desc.add<double>("deltaR", -1.0)->setComment("maximum deltaR of track to jet. If -ve just use cut from JTA");
0070     desc.add<double>("maximumDecayLength", 999999.0);
0071     desc.add<double>("maximumDistanceToJetAxis", 999999.0);
0072     desc.add<std::string>("trackQualityClass", "any");
0073     desc.add<bool>("useVariableJTA", false);
0074     desc.add<double>("a_dR", -0.001053);
0075     desc.add<double>("b_dR", 0.6263);
0076     desc.add<double>("a_pT", 0.005263);
0077     desc.add<double>("b_pT", 0.3684);
0078     desc.add<double>("min_pT", 120);
0079     desc.add<double>("max_pT", 500);
0080     desc.add<double>("min_pT_dRcut", 0.5);
0081     desc.add<double>("max_pT_dRcut", 0.1);
0082     desc.add<double>("max_pT_trackPTcut", 3.);
0083   }
0084
0085 protected:
0086   std::multiset<float> orderedSignificances(const TagInfo& tkip) const {
0087     const std::vector<reco::btag::TrackIPData>& impactParameters((tkip.impactParameterData()));
0088     const Container& tracks(tkip.selectedTracks());
0089     std::multiset<float> significances;
0090     int i = 0;
0091     if (tkip.primaryVertex().isNull()) {
0092       return std::multiset<float>();
0093     }
0094
0095     GlobalPoint pv(tkip.primaryVertex()->position().x(),
0096                    tkip.primaryVertex()->position().y(),
0097                    tkip.primaryVertex()->position().z());
0098
0099     for (std::vector<reco::btag::TrackIPData>::const_iterator it = impactParameters.begin();
0100          it != impactParameters.end();
0101          ++it, i++) {
0102       if (fabs(impactParameters[i].distanceToJetAxis.value()) < m_cutMaxDistToAxis &&  // distance to JetAxis
0103           (impactParameters[i].closestToJetAxis - pv).mag() < m_cutMaxDecayLen &&      // max decay len
0104           (m_useAllQualities == true ||
0105            reco::btag::toTrack(tracks[i])->quality(m_trackQuality)) &&       // use selected track qualities
0106           (fabs(((m_ipType == 0) ? it->ip3d : it->ip2d).value()) > m_minIP)  // minimum impact parameter
0107       ) {
0108         //calculate the signed or un-signed significance
0109         float signed_sig = ((m_ipType == 0) ? it->ip3d : it->ip2d).significance();
0110         float unsigned_sig = fabs(signed_sig);
0111         float significance = (m_useSignedIPSig) ? signed_sig : unsigned_sig;
0112
0113         if (useVariableJTA_) {
0114           if (tkip.variableJTA(varJTApars)[i])
0115             significances.insert(significance);
0116         } else  // no using variable JTA, use the default method
0117           if (m_deltaR <= 0 ||
0118               ROOT::Math::VectorUtil::DeltaR((*tkip.jet()).p4().Vect(), (*tracks[i]).momentum()) < m_deltaR)
0119             significances.insert(significance);
0120       }
0121     }
0122
0123     return significances;
0124   }
0125
0126   bool useVariableJTA_;
0127   reco::btag::variableJTAParameters varJTApars;
0128
0129   double m_minIP;
0130   bool m_useSignedIPSig;
0131
0132   int m_nthTrack;
0133   int m_ipType;
0134   double m_deltaR;
0135   double m_cutMaxDecayLen;
0136   double m_cutMaxDistToAxis;
0137   reco::TrackBase::TrackQuality m_trackQuality;
0138   bool m_useAllQualities;
0139 };
0140
0141 #endif  // ImpactParameter_TemplatedTrackCountingComputer_h