Skip to content

Commit

Permalink
code-format
Browse files Browse the repository at this point in the history
  • Loading branch information
@drkovalskyi
drkovalskyi committed Dec 21, 2023
1 parent f03ffd2 commit 14e6996
Show file tree
Hide file tree
Showing 6 changed files with 82 additions and 95 deletions.
2 changes: 1 addition & 1 deletion DataFormats/PatCandidates/interface/Muon.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -416,7 +416,7 @@ namespace pat {
/// Muon MVA
float mvaIDValue_;
float softMvaValue_;
float softMvaRun3Value_=0;
float softMvaRun3Value_ = 0;

/// Inverse beta
float inverseBeta_;
Expand Down
2 changes: 1 addition & 1 deletion PhysicsTools/PatAlgos/interface/SoftMuonMvaRun3Estimator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,5 +9,5 @@ namespace pat {
class Muon;

float computeSoftMvaRun3(XGBooster& booster, const Muon& muon);
}
} // namespace pat
#endif
16 changes: 8 additions & 8 deletions PhysicsTools/PatAlgos/interface/XGBooster.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,27 +8,27 @@
#include <xgboost/c_api.h>

namespace pat {
class XGBooster {
class XGBooster {
public:
XGBooster(std::string model_file);
XGBooster(std::string model_file, std::string model_features);

/// Features need to be entered in the order they are used
/// in the model
void addFeature(std::string name);

/// Reset feature values
void reset();

void set(std::string name, float value);

float predict();

private:
std::vector<float> features_;
std::map<std::string,unsigned int> feature_name_to_index_;
BoosterHandle booster_;
std::map<std::string, unsigned int> feature_name_to_index_;
BoosterHandle booster_;
};
}
} // namespace pat

#endif
4 changes: 2 additions & 2 deletions PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -479,8 +479,8 @@ PATMuonProducer::PATMuonProducer(const edm::ParameterSet& iConfig, PATMuonHeavyO
if (computeSoftMuonMVA_) {
std::string softMvaRun3Model = iConfig.getParameter<string>("softMvaRun3Model");
softMuonMvaRun3Booster_ =
std::make_unique<pat::XGBooster>(edm::FileInPath(softMvaRun3Model + ".model").fullPath(),
edm::FileInPath(softMvaRun3Model + ".features").fullPath());
std::make_unique<pat::XGBooster>(edm::FileInPath(softMvaRun3Model + ".model").fullPath(),
edm::FileInPath(softMvaRun3Model + ".features").fullPath());
}

addTriggerMatching_ = iConfig.getParameter<bool>("addTriggerMatching");
Expand Down
98 changes: 48 additions & 50 deletions PhysicsTools/PatAlgos/src/SoftMuonMvaRun3Estimator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,110 +4,105 @@

typedef std::pair<const reco::MuonChamberMatch*, const reco::MuonSegmentMatch*> MatchPair;

const MatchPair&
getBetterMatch(const MatchPair& match1, const MatchPair& match2){

const MatchPair& getBetterMatch(const MatchPair& match1, const MatchPair& match2) {
// Prefer DT over CSC simply because it's closer to IP
// and will have less multiple scattering (at least for
// RB1 vs ME1/3 case). RB1 & ME1/2 overlap is tiny
if (match2.first->detector() == MuonSubdetId::DT and
match1.first->detector() != MuonSubdetId::DT)
if (match2.first->detector() == MuonSubdetId::DT and match1.first->detector() != MuonSubdetId::DT)
return match2;

// For the rest compare local x match. We expect that
// segments belong to the muon, so the difference in
// local x is a reflection on how well we can measure it
if ( abs(match1.first->x - match1.second->x) >
abs(match2.first->x - match2.second->x) )
if (abs(match1.first->x - match1.second->x) > abs(match2.first->x - match2.second->x))
return match2;

return match1;
}

float dX(const MatchPair& match){
float dX(const MatchPair& match) {
if (match.first and match.second->hasPhi())
return (match.first->x - match.second->x);
else
return 9999.;
}

float pullX(const MatchPair& match){
float pullX(const MatchPair& match) {
if (match.first and match.second->hasPhi())
return dX(match) /
sqrt(pow(match.first->xErr, 2) + pow(match.second->xErr, 2));
return dX(match) / sqrt(pow(match.first->xErr, 2) + pow(match.second->xErr, 2));
else
return 9999.;
}

float pullDxDz(const MatchPair& match){
float pullDxDz(const MatchPair& match) {
if (match.first and match.second->hasPhi())
return (match.first->dXdZ - match.second->dXdZ) /
sqrt(pow(match.first->dXdZErr, 2) + pow(match.second->dXdZErr, 2));
else
return 9999.;
}

float dY(const MatchPair& match){
float dY(const MatchPair& match) {
if (match.first and match.second->hasZed())
return (match.first->y - match.second->y);
else
return 9999.;
}

float pullY(const MatchPair& match){
float pullY(const MatchPair& match) {
if (match.first and match.second->hasZed())
return dY(match) /
sqrt(pow(match.first->yErr, 2) + pow(match.second->yErr, 2));
return dY(match) / sqrt(pow(match.first->yErr, 2) + pow(match.second->yErr, 2));
else
return 9999.;
}

float pullDyDz(const MatchPair& match){
float pullDyDz(const MatchPair& match) {
if (match.first and match.second->hasZed())
return (match.first->dYdZ - match.second->dYdZ) /
sqrt(pow(match.first->dYdZErr, 2) + pow(match.second->dYdZErr, 2));
else
return 9999.;
}

void fillMatchInfoForStation(std::string prefix, pat::XGBooster& booster, const MatchPair& match){
booster.set(prefix + "_dX", dX(match));
booster.set(prefix + "_pullX", pullX(match));
void fillMatchInfoForStation(std::string prefix, pat::XGBooster& booster, const MatchPair& match) {
booster.set(prefix + "_dX", dX(match));
booster.set(prefix + "_pullX", pullX(match));
booster.set(prefix + "_pullDxDz", pullDxDz(match));
booster.set(prefix + "_dY", dY(match));
booster.set(prefix + "_pullY", pullY(match));
booster.set(prefix + "_dY", dY(match));
booster.set(prefix + "_pullY", pullY(match));
booster.set(prefix + "_pullDyDz", pullDyDz(match));
}

void fillMatchInfo(pat::XGBooster& booster, const pat::Muon& muon){
void fillMatchInfo(pat::XGBooster& booster, const pat::Muon& muon) {
// Initiate containter for results
const int n_stations = 2;
std::vector<MatchPair> matches;
for (unsigned int i=0; i < n_stations; ++i)
for (unsigned int i = 0; i < n_stations; ++i)
matches.push_back(std::pair(nullptr, nullptr));

// Find best matches
for (auto& chamberMatch : muon.matches()){
for (auto& chamberMatch : muon.matches()) {
unsigned int station = chamberMatch.station() - 1;
if (station >= n_stations) continue;
if (station >= n_stations)
continue;

// Find best segment match.
// We could consider all segments, but we will restrict to segments
// that match to this candidate better than to other muon candidates
for (auto& segmentMatch : chamberMatch.segmentMatches){
if ( not segmentMatch.isMask(reco::MuonSegmentMatch::BestInStationByDR) ||
not segmentMatch.isMask(reco::MuonSegmentMatch::BelongsToTrackByDR) )
continue;
for (auto& segmentMatch : chamberMatch.segmentMatches) {
if (not segmentMatch.isMask(reco::MuonSegmentMatch::BestInStationByDR) ||
not segmentMatch.isMask(reco::MuonSegmentMatch::BelongsToTrackByDR))
continue;

// Multiple segment matches are possible in different
// chambers that are either overlapping or belong to
// different detectors. We need to select one.
auto match_pair = MatchPair(&chamberMatch, &segmentMatch);

if (matches[station].first)
matches[station] = getBetterMatch(matches[station], match_pair);
matches[station] = getBetterMatch(matches[station], match_pair);
else
matches[station] = match_pair;
matches[station] = match_pair;
}
}

Expand All @@ -117,24 +112,27 @@ void fillMatchInfo(pat::XGBooster& booster, const pat::Muon& muon){
}

float pat::computeSoftMvaRun3(pat::XGBooster& booster, const pat::Muon& muon) {
if (!muon.isTrackerMuon() && !muon.isGlobalMuon()) return 0;
if (!muon.isTrackerMuon() && !muon.isGlobalMuon())
return 0;

fillMatchInfo(booster, muon);
booster.set("pt", muon.pt());
booster.set("eta", muon.eta());
booster.set("trkValidFrac", muon.innerTrack()->validFraction());

booster.set("pt", muon.pt());
booster.set("eta", muon.eta());
booster.set("trkValidFrac", muon.innerTrack()->validFraction());
booster.set("glbTrackProbability", muon.combinedQuality().glbTrackProbability);
booster.set("nLostHitsInner", muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS));
booster.set("nLostHitsOuter", muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS));
booster.set("trkKink", muon.combinedQuality().trkKink);
booster.set("chi2LocalPosition", muon.combinedQuality().chi2LocalPosition);
booster.set("nPixels", muon.innerTrack()->hitPattern().numberOfValidPixelHits());
booster.set("nValidHits", muon.innerTrack()->hitPattern().numberOfValidTrackerHits());
booster.set("nLostHitsOn", muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::TRACK_HITS));
booster.set("glbNormChi2", muon.isGlobalMuon() ? muon.globalTrack()->normalizedChi2() : 9999.);
booster.set("trkLayers", muon.innerTrack()->hitPattern().trackerLayersWithMeasurement());
booster.set("highPurity", muon.innerTrack()->quality(reco::Track::highPurity));

booster.set("nLostHitsInner",
muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS));
booster.set("nLostHitsOuter",
muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS));
booster.set("trkKink", muon.combinedQuality().trkKink);
booster.set("chi2LocalPosition", muon.combinedQuality().chi2LocalPosition);
booster.set("nPixels", muon.innerTrack()->hitPattern().numberOfValidPixelHits());
booster.set("nValidHits", muon.innerTrack()->hitPattern().numberOfValidTrackerHits());
booster.set("nLostHitsOn", muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::TRACK_HITS));
booster.set("glbNormChi2", muon.isGlobalMuon() ? muon.globalTrack()->normalizedChi2() : 9999.);
booster.set("trkLayers", muon.innerTrack()->hitPattern().trackerLayersWithMeasurement());
booster.set("highPurity", muon.innerTrack()->quality(reco::Track::highPurity));

return booster.predict();
}
55 changes: 22 additions & 33 deletions PhysicsTools/PatAlgos/src/XGBooster.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ std::vector<std::string> read_features(const std::string& content) {

while (stream) {
stream >> ch;

if (ch == ']') {
break;
} else if (ch == ',') {
Expand All @@ -43,8 +43,7 @@ std::vector<std::string> read_features(const std::string& content) {
return result;
}

XGBooster::XGBooster(std::string model_file)
{
XGBooster::XGBooster(std::string model_file) {
int status = XGBoosterCreate(nullptr, 0, &booster_);
if (status != 0)
throw std::runtime_error("Failed to create XGBooster");
Expand All @@ -54,59 +53,49 @@ XGBooster::XGBooster(std::string model_file)
XGBoosterSetParam(booster_, "nthread", "1");
}

XGBooster::XGBooster(std::string model_file, std::string model_features):XGBooster(model_file)
{
XGBooster::XGBooster(std::string model_file, std::string model_features) : XGBooster(model_file) {
std::ifstream file(model_features);
if (!file.is_open())
throw std::runtime_error("Failed to open file: " + model_features);

std::string content((std::istreambuf_iterator<char>(file)), std::istreambuf_iterator<char>());
file.close();

std::vector<std::string> features = read_features(content);

for (const auto& feature : features) {
addFeature(feature);
}

}


void XGBooster::reset()
{
std::fill(features_.begin(), features_.end(), std::nan(""));
}

void XGBooster::reset() { std::fill(features_.begin(), features_.end(), std::nan("")); }

void XGBooster::addFeature(std::string name){
void XGBooster::addFeature(std::string name) {
features_.push_back(0);
feature_name_to_index_[name] = features_.size()-1;
feature_name_to_index_[name] = features_.size() - 1;
}

void XGBooster::set(std::string name, float value){
features_.at(feature_name_to_index_[name]) = value;
}
void XGBooster::set(std::string name, float value) { features_.at(feature_name_to_index_[name]) = value; }

float XGBooster::predict()
{
float XGBooster::predict() {
float result(-999.);

// check if all feature values are set properly
for (unsigned int i = 0; i < features_.size(); ++i)
if (std::isnan(features_.at(i))) {
std::string feature_name;
for (const auto& pair : feature_name_to_index_) {
if (pair.second == i) {
feature_name = pair.first;
break;
}
if (pair.second == i) {
feature_name = pair.first;
break;
}
}
throw std::runtime_error("Feature is not set: " + feature_name);
}

DMatrixHandle dvalues;
XGDMatrixCreateFromMat(&features_[0], 1, features_.size(), 9e99, &dvalues);

bst_ulong out_len = 0;
const float* score = nullptr;

Expand All @@ -126,12 +115,12 @@ float XGBooster::predict()

XGDMatrixFree(dvalues);

if (ret==0) {
assert(out_len==1 && "Unexpected prediction format");
if (ret == 0) {
assert(out_len == 1 && "Unexpected prediction format");
result = score[0];
}

reset();
return result;

return result;
}

0 comments on commit 14e6996

Please sign in to comment.