Combined/src/CandidateChargeBTagComputer.cc

0001 #include <memory>
0002
0003 #include "RecoBTag/Combined/interface/CandidateChargeBTagComputer.h"
0004
0005 CandidateChargeBTagComputer::Tokens::Tokens(const edm::ParameterSet& parameters, edm::ESConsumesCollector&& cc) {
0006   if (parameters.getParameter<bool>("useCondDB")) {
0007     gbrForest_ = cc.consumes(edm::ESInputTag{"", parameters.getParameter<std::string>("gbrForestLabel")});
0008   }
0009 }
0010
0011 CandidateChargeBTagComputer::CandidateChargeBTagComputer(const edm::ParameterSet& parameters, Tokens tokens)
0012     : weightFile_(parameters.getParameter<edm::FileInPath>("weightFile")),
0013       useAdaBoost_(parameters.getParameter<bool>("useAdaBoost")),
0014       jetChargeExp_(parameters.getParameter<double>("jetChargeExp")),
0015       svChargeExp_(parameters.getParameter<double>("svChargeExp")),
0016       tokens_{tokens} {
0017   uses(0, "pfImpactParameterTagInfos");
0018   uses(1, "pfInclusiveSecondaryVertexFinderCvsLTagInfos");
0019   uses(2, "softPFMuonsTagInfos");
0020   uses(3, "softPFElectronsTagInfos");
0021
0022   mvaID = std::make_unique<TMVAEvaluator>();
0023 }
0024
0025 CandidateChargeBTagComputer::~CandidateChargeBTagComputer() {}
0026
0027 void CandidateChargeBTagComputer::initialize(const JetTagComputerRecord& record) {
0028   // Saving MVA variable names;
0029   // names and order need to be the same as in the training
0030   std::vector<std::string> variables({"tr_ch_inc",
0031                                       "sv_ch",
0032                                       "mu_sip2d",
0033                                       /*"mu_sip3d",*/ "mu_delR",
0034                                       "mu_ptrel",
0035                                       "mu_pfrac",
0036                                       "mu_char",
0037                                       /*"el_sip2d","el_sip3d",*/ "el_delR",
0038                                       "el_ptrel",
0039                                       "el_pfrac",
0040                                       "el_mva",
0041                                       "el_char",
0042                                       "pt1_ch/j_pt",
0043                                       "pt2_ch/j_pt",
0044                                       "pt3_ch/j_pt"});
0045   std::vector<std::string> spectators(0);
0046
0047   if (tokens_.gbrForest_.isInitialized()) {
0048     mvaID->initializeGBRForest(&record.get(tokens_.gbrForest_), variables, spectators, useAdaBoost_);
0049   } else
0050     mvaID->initialize("Color:Silent:Error", "BDT", weightFile_.fullPath(), variables, spectators, true, useAdaBoost_);
0051 }
0052
0053 float CandidateChargeBTagComputer::discriminator(const TagInfoHelper& tagInfo) const {
0054   // get TagInfos
0055   const reco::CandIPTagInfo& ip_info = tagInfo.get<reco::CandIPTagInfo>(0);
0056   const reco::CandSecondaryVertexTagInfo& sv_info = tagInfo.get<reco::CandSecondaryVertexTagInfo>(1);
0057   const reco::CandSoftLeptonTagInfo& sm_info = tagInfo.get<reco::CandSoftLeptonTagInfo>(2);
0058   const reco::CandSoftLeptonTagInfo& se_info = tagInfo.get<reco::CandSoftLeptonTagInfo>(3);
0059
0060   size_t n_ip_info = ip_info.jet()->getJetConstituents().size();
0061   size_t n_sv_info = sv_info.nVertices();
0062   size_t n_sm_info = sm_info.leptons();
0063   size_t n_se_info = se_info.leptons();
0064
0065   // default discriminator value
0066   float value = -10.;
0067
0068   // if no tag info is present, MVA not computed and default discriminator value returned
0069   if (n_ip_info + n_sv_info + n_sm_info + n_se_info > 0) {
0070     // default variable values
0071     float tr_ch_inc = 0;
0072     float pt_ratio1_ch = 0;
0073     float pt_ratio2_ch = 0;
0074     float pt_ratio3_ch = 0;
0075
0076     float sv_ch = 0;
0077
0078     float mu_sip2d = 0;
0079     //float mu_sip3d = 0;
0080     float mu_delR = 0;
0081     float mu_ptrel = 0;
0082     float mu_pfrac = 0;
0083     int mu_char = 0;
0084
0085     //float el_sip2d = 0;
0086     //float el_sip3d = 0;
0087     float el_delR = 0;
0088     float el_ptrel = 0;
0089     float el_pfrac = 0;
0090     float el_mva = 0;
0091     int el_char = 0;
0092
0093     // compute jet-charge
0094     float tr_ch_num = 0;
0095     float tr_ch_den = 0;
0096
0097     // loop over tracks associated to the jet
0098     for (size_t i_ip = 0; i_ip < n_ip_info; ++i_ip) {
0099       const reco::Candidate* trackPtr = ip_info.jet()->getJetConstituents().at(i_ip).get();
0100       if (trackPtr->charge() == 0)
0101         continue;
0102
0103       float tr_ch_weight = pow(trackPtr->pt(), jetChargeExp_);
0104       tr_ch_num += tr_ch_weight * trackPtr->charge();
0105       tr_ch_den += tr_ch_weight;
0106
0107       // find the three higher-pt tracks
0108       if (fabs(pt_ratio1_ch) < trackPtr->pt()) {
0109         pt_ratio3_ch = pt_ratio2_ch;
0110         pt_ratio2_ch = pt_ratio1_ch;
0111         pt_ratio1_ch = trackPtr->pt() * trackPtr->charge();
0112       } else if (fabs(pt_ratio2_ch) < trackPtr->pt()) {
0113         pt_ratio3_ch = pt_ratio2_ch;
0114         pt_ratio2_ch = trackPtr->pt() * trackPtr->charge();
0115       } else if (fabs(pt_ratio3_ch) < trackPtr->pt()) {
0116         pt_ratio3_ch = trackPtr->pt() * trackPtr->charge();
0117       }
0118     }
0119     if (n_ip_info > 0) {
0120       float jet_pt = ip_info.jet()->pt();
0121       if (jet_pt > 0) {
0122         pt_ratio1_ch = pt_ratio1_ch / jet_pt;
0123         pt_ratio2_ch = pt_ratio2_ch / jet_pt;
0124         pt_ratio3_ch = pt_ratio3_ch / jet_pt;
0125       }
0126     }
0127
0128     if (tr_ch_den > 0)
0129       tr_ch_inc = tr_ch_num / tr_ch_den;
0130
0131     // compute secondary vertex charge
0132     if (n_sv_info > 0) {
0133       float jet_pt = sv_info.jet()->pt();
0134
0135       float sv_ch_num = 0;
0136       float sv_ch_den = 0;
0137
0138       // find the selected secondary vertex with highest invariant mass
0139       int vtx_idx = -1;
0140       float max_mass = 0;
0141       for (size_t i_vtx = 0; i_vtx < n_sv_info; ++i_vtx) {
0142         float sv_mass = sv_info.secondaryVertex(i_vtx).p4().mass();
0143         float sv_chi2 = sv_info.secondaryVertex(i_vtx).vertexNormalizedChi2();
0144         float sv_pfrac = sv_info.secondaryVertex(i_vtx).pt() / jet_pt;
0145         float sv_L = sv_info.flightDistance(i_vtx).value();
0146         float sv_sL = sv_info.flightDistance(i_vtx).significance();
0147         float delEta = sv_info.secondaryVertex(i_vtx).momentum().eta() - sv_info.flightDirection(i_vtx).eta();
0148         float delPhi = sv_info.secondaryVertex(i_vtx).momentum().phi() - sv_info.flightDirection(i_vtx).phi();
0149         if (fabs(delPhi) > M_PI)
0150           delPhi = 2 * M_PI - fabs(delPhi);
0151         float sv_delR = sqrt(delEta * delEta + delPhi * delPhi);
0152
0153         if (sv_mass > max_mass && sv_mass > 1.4 && sv_chi2 < 3 && sv_chi2 > 0 && sv_pfrac > 0.25 && sv_L < 2.5 &&
0154             sv_sL > 4 && sv_delR < 0.1) {
0155           max_mass = sv_mass;
0156           vtx_idx = i_vtx;
0157         }
0158       }
0159
0160       if (vtx_idx >= 0) {
0161         // loop over tracks associated to the vertex
0162         size_t n_sv_tracks = sv_info.vertexTracks(vtx_idx).size();
0163         for (size_t i_tr = 0; i_tr < n_sv_tracks; ++i_tr) {
0164           const reco::CandidatePtr trackRef = sv_info.vertexTracks(vtx_idx)[i_tr];
0165           const reco::Track* trackPtr = reco::btag::toTrack(trackRef);
0166
0167           float sv_ch_weight = pow(trackPtr->pt(), svChargeExp_);
0168           sv_ch_num += sv_ch_weight * trackPtr->charge();
0169           sv_ch_den += sv_ch_weight;
0170         }
0171
0172         if (sv_ch_den > 0)
0173           sv_ch = sv_ch_num / sv_ch_den;
0174       }
0175     }
0176
0177     // fill soft muon variables
0178     if (n_sm_info > 0) {
0179       // find the muon with higher transverse momentum
0180       int lep_idx = 0;
0181       float max_pt = 0;
0182       for (size_t i_lep = 0; i_lep < n_sm_info; ++i_lep) {
0183         float lep_pt = sm_info.lepton(0)->pt();
0184         if (lep_pt > max_pt) {
0185           max_pt = lep_pt;
0186           lep_idx = i_lep;
0187         }
0188       }
0189
0190       mu_sip2d = sm_info.properties(lep_idx).sip2d;
0191       //mu_sip3d = sm_info.properties(lep_idx).sip3d;
0192       mu_delR = sm_info.properties(lep_idx).deltaR;
0193       mu_ptrel = sm_info.properties(lep_idx).ptRel;
0194       mu_pfrac = sm_info.properties(lep_idx).ratio;
0195       mu_char = sm_info.properties(lep_idx).charge;
0196     }
0197
0198     // fill soft electron variables
0199     if (n_se_info > 0) {
0200       // find the electron with higher transverse momentum
0201       int lep_idx = 0;
0202       float max_pt = 0;
0203       for (size_t i_lep = 0; i_lep < n_se_info; ++i_lep) {
0204         float lep_pt = se_info.lepton(0)->pt();
0205         if (lep_pt > max_pt) {
0206           max_pt = lep_pt;
0207           lep_idx = i_lep;
0208         }
0209       }
0210
0211       //el_sip2d = se_info.properties(lep_idx).sip2d;
0212       //el_sip3d = se_info.properties(lep_idx).sip3d;
0213       el_delR = se_info.properties(lep_idx).deltaR;
0214       el_ptrel = se_info.properties(lep_idx).ptRel;
0215       el_pfrac = se_info.properties(lep_idx).ratio;
0216       el_mva = se_info.properties(lep_idx).elec_mva;
0217       el_char = se_info.properties(lep_idx).charge;
0218     }
0219
0220     std::map<std::string, float> inputs;
0221     inputs["tr_ch_inc"] = tr_ch_inc;
0222     inputs["pt1_ch/j_pt"] = pt_ratio1_ch;
0223     inputs["pt2_ch/j_pt"] = pt_ratio2_ch;
0224     inputs["pt3_ch/j_pt"] = pt_ratio3_ch;
0225     inputs["sv_ch"] = sv_ch;
0226     inputs["mu_sip2d"] = mu_sip2d;
0227     //inputs["mu_sip3d"] = mu_sip3d;
0228     inputs["mu_delR"] = mu_delR;
0229     inputs["mu_ptrel"] = mu_ptrel;
0230     inputs["mu_pfrac"] = mu_pfrac;
0231     inputs["mu_char"] = mu_char;
0232     //inputs["el_sip2d"] = el_sip2d;
0233     //inputs["el_sip3d"] = el_sip3d;
0234     inputs["el_delR"] = el_delR;
0235     inputs["el_ptrel"] = el_ptrel;
0236     inputs["el_pfrac"] = el_pfrac;
0237     inputs["el_mva"] = el_mva;
0238     inputs["el_char"] = el_char;
0239
0240     // evaluate the MVA
0241     value = mvaID->evaluate(inputs);
0242   }
0243
0244   // return the final discriminator value
0245   return value;
0246 }
0247
0248 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
0249 void CandidateChargeBTagComputer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0250   edm::ParameterSetDescription desc;
0251   desc.add<bool>("useCondDB", false);
0252   desc.add<std::string>("gbrForestLabel", "");
0253   desc.add<edm::FileInPath>("weightFile", edm::FileInPath());
0254   desc.add<bool>("useAdaBoost", true);
0255   desc.add<double>("jetChargeExp", 0.8);
0256   desc.add<double>("svChargeExp", 0.5);
0257   descriptions.addDefault(desc);
0258 }