File indexing completed on 2023-03-17 11:17:00
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
0029
0030 std::vector<std::string> variables({"tr_ch_inc",
0031 "sv_ch",
0032 "mu_sip2d",
0033 "mu_delR",
0034 "mu_ptrel",
0035 "mu_pfrac",
0036 "mu_char",
0037 "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
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
0066 float value = -10.;
0067
0068
0069 if (n_ip_info + n_sv_info + n_sm_info + n_se_info > 0) {
0070
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
0080 float mu_delR = 0;
0081 float mu_ptrel = 0;
0082 float mu_pfrac = 0;
0083 int mu_char = 0;
0084
0085
0086
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
0094 float tr_ch_num = 0;
0095 float tr_ch_den = 0;
0096
0097
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
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
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
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
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
0178 if (n_sm_info > 0) {
0179
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
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
0199 if (n_se_info > 0) {
0200
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
0212
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
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
0233
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
0241 value = mvaID->evaluate(inputs);
0242 }
0243
0244
0245 return value;
0246 }
0247
0248
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 }