TrackerSeedGenerator/src/SeedMvaEstimator.cc

0001 #include "RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimator.h"
0002
0003 #include "DataFormats/Math/interface/deltaPhi.h"
0004 #include "DataFormats/Math/interface/deltaR.h"
0005
0006 #include "FWCore/ParameterSet/interface/FileInPath.h"
0007 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0008
0009 SeedMvaEstimator::SeedMvaEstimator(const edm::FileInPath& weightsfile,
0010                                    const std::vector<double>& scale_mean,
0011                                    const std::vector<double>& scale_std,
0012                                    const bool isFromL1,
0013                                    const int minL1Qual)
0014     : scale_mean_(scale_mean), scale_std_(scale_std), isFromL1_(isFromL1), minL1Qual_(minL1Qual) {
0015   gbrForest_ = createGBRForest(weightsfile);
0016 }
0017
0018 SeedMvaEstimator::~SeedMvaEstimator() {}
0019
0020 namespace {
0021   enum inputIndexes {
0022     kTsosErr0,       // 0
0023     kTsosErr2,       // 1
0024     kTsosErr5,       // 2
0025     kTsosDxdz,       // 3
0026     kTsosDydz,       // 4
0027     kTsosQbp,        // 5
0028     kDRdRL1SeedP,    // 6
0029     kDPhidRL1SeedP,  // 7
0030     kLastL1,         // 8
0031
0032     kDRdRL2SeedP = 8,  // 8
0033     kDPhidRL2SeedP,    // 9
0034     kLastL2,           // 10
0035   };
0036 }  // namespace
0037
0038 void SeedMvaEstimator::getL1MuonVariables(const GlobalVector& global_p,
0039                                           const l1t::MuonBxCollection& l1Muons,
0040                                           float& dR2dRL1SeedP,
0041                                           float& dPhidRL1SeedP) const {
0042   for (int ibx = l1Muons.getFirstBX(); ibx <= l1Muons.getLastBX(); ++ibx) {
0043     if (ibx != 0)
0044       continue;  // -- only take when ibx == 0 -- //
0045
0046     for (auto it = l1Muons.begin(ibx); it != l1Muons.end(ibx); it++) {
0047       if (it->hwQual() < minL1Qual_)
0048         continue;
0049
0050       float dR2tmp = reco::deltaR2(it->etaAtVtx(), it->phiAtVtx(), global_p.eta(), global_p.phi());
0051       if (dR2tmp < dR2dRL1SeedP) {
0052         dR2dRL1SeedP = dR2tmp;
0053         dPhidRL1SeedP = reco::deltaPhi(it->phiAtVtx(), global_p.phi());
0054       }
0055     }
0056   }
0057 }
0058
0059 void SeedMvaEstimator::getL2MuonVariables(const GlobalVector& global_p,
0060                                           const reco::RecoChargedCandidateCollection& l2Muons,
0061                                           float& dR2dRL2SeedP,
0062                                           float& dPhidRL2SeedP) const {
0063   for (auto it = l2Muons.begin(); it != l2Muons.end(); it++) {
0064     float dR2tmp = reco::deltaR2(*it, global_p);
0065     if (dR2tmp < dR2dRL2SeedP) {
0066       dR2dRL2SeedP = dR2tmp;
0067       dPhidRL2SeedP = reco::deltaPhi(it->phi(), global_p.phi());
0068     }
0069   }
0070 }
0071
0072 double SeedMvaEstimator::computeMva(const TrajectorySeed& seed,
0073                                     const GlobalVector& global_p,
0074                                     const l1t::MuonBxCollection& l1Muons,
0075                                     const reco::RecoChargedCandidateCollection& l2Muons) const {
0076   static constexpr float initDRdPhi(99999.);
0077   auto kLast = isFromL1_ ? kLastL1 : kLastL2;
0078   float var[kLast];
0079
0080   var[kTsosErr0] = seed.startingState().error(0);
0081   var[kTsosErr2] = seed.startingState().error(2);
0082   var[kTsosErr5] = seed.startingState().error(5);
0083   var[kTsosDxdz] = seed.startingState().parameters().dxdz();
0084   var[kTsosDydz] = seed.startingState().parameters().dydz();
0085   var[kTsosQbp] = seed.startingState().parameters().qbp();
0086
0087   float dR2dRL1SeedP = initDRdPhi;
0088   float dPhidRL1SeedP = initDRdPhi;
0089   getL1MuonVariables(global_p, l1Muons, dR2dRL1SeedP, dPhidRL1SeedP);
0090
0091   var[kDRdRL1SeedP] = std::sqrt(dR2dRL1SeedP);
0092   var[kDPhidRL1SeedP] = dPhidRL1SeedP;
0093
0094   if (!isFromL1_) {
0095     float dR2dRL2SeedP = initDRdPhi;
0096     float dPhidRL2SeedP = initDRdPhi;
0097     getL2MuonVariables(global_p, l2Muons, dR2dRL2SeedP, dPhidRL2SeedP);
0098
0099     var[kDRdRL2SeedP] = std::sqrt(dR2dRL2SeedP);
0100     var[kDPhidRL2SeedP] = dPhidRL2SeedP;
0101   }
0102
0103   for (int iv = 0; iv < kLast; ++iv) {
0104     var[iv] = (var[iv] - scale_mean_.at(iv)) / scale_std_.at(iv);
0105   }
0106
0107   return gbrForest_->GetResponse(var);
0108 }