Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:17 

0001 #include "RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimatorPhase2.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 #include <cmath>
0010 
0011 using namespace std;
0012 
0013 namespace Phase2 {
0014   enum inputIndexesPhase2 {
0015     kTsosErr0,           // 0
0016     kTsosErr2,           // 1
0017     kTsosErr5,           // 2
0018     kTsosDxdz,           // 3
0019     kTsosDydz,           // 4
0020     kTsosQbp,            // 5
0021     kDRdRL1TkMuSeedP,    // 6
0022     kDRdPhiL1TkMuSeedP,  // 7
0023     kExpd2HitL1Tk1,      // 8
0024     kExpd2HitL1Tk2,      // 9
0025     kExpd2HitL1Tk3,      // 10
0026     kLast                // 11
0027   };
0028 }
0029 
0030 SeedMvaEstimatorPhase2::SeedMvaEstimatorPhase2(const edm::FileInPath& weightsfile,
0031                                                const std::vector<double>& scale_mean,
0032                                                const std::vector<double>& scale_std)
0033     : scale_mean_(scale_mean), scale_std_(scale_std) {
0034   gbrForest_ = createGBRForest(weightsfile);
0035 }
0036 
0037 SeedMvaEstimatorPhase2::~SeedMvaEstimatorPhase2() {}
0038 
0039 void SeedMvaEstimatorPhase2::getL1TTVariables(const TrajectorySeed& seed,
0040                                               const GlobalVector& global_p,
0041                                               const GlobalPoint& global_x,
0042                                               const edm::Handle<l1t::TrackerMuonCollection>& h_L1TkMu,
0043                                               float& DRL1TkMu,
0044                                               float& DPhiL1TkMu) const {
0045   for (auto L1TkMu = h_L1TkMu->begin(); L1TkMu != h_L1TkMu->end(); ++L1TkMu) {
0046     auto TkRef = L1TkMu->trkPtr();
0047     float DRL1TkMu_tmp = reco::deltaR(*TkRef, global_p);
0048     float DPhiL1TkMu_tmp = reco::deltaPhi(TkRef->phi(), global_p.phi());
0049     if (DRL1TkMu_tmp < DRL1TkMu) {
0050       DRL1TkMu = DRL1TkMu_tmp;
0051       DPhiL1TkMu = DPhiL1TkMu_tmp;
0052     }
0053   }
0054 }
0055 
0056 vector<LayerTSOS> SeedMvaEstimatorPhase2::getTsosOnPixels(const TTTrack<Ref_Phase2TrackerDigi_>& l1tk,
0057                                                           const edm::ESHandle<MagneticField>& magfieldH,
0058                                                           const Propagator& propagatorAlong,
0059                                                           const GeometricSearchTracker& geomTracker) const {
0060   vector<LayerTSOS> v_tsos = {};
0061 
0062   std::vector<BarrelDetLayer const*> const& bpix = geomTracker.pixelBarrelLayers();
0063   std::vector<ForwardDetLayer const*> const& nfpix = geomTracker.negPixelForwardLayers();
0064   std::vector<ForwardDetLayer const*> const& pfpix = geomTracker.posPixelForwardLayers();
0065 
0066   int chargeTk = l1tk.rInv() > 0. ? 1 : -1;
0067   GlobalPoint gpos = l1tk.POCA();
0068   GlobalVector gmom = l1tk.momentum();
0069 
0070   FreeTrajectoryState fts = FreeTrajectoryState(gpos, gmom, chargeTk, magfieldH.product());
0071 
0072   for (auto it = bpix.begin(); it != bpix.end(); ++it) {
0073     TrajectoryStateOnSurface tsos = propagatorAlong.propagate(fts, (**it).specificSurface());
0074     if (!tsos.isValid())
0075       continue;
0076 
0077     auto z0 = std::abs(tsos.globalPosition().z() - (**it).specificSurface().position().z());
0078     auto deltaZ = 0.5f * (**it).surface().bounds().length();
0079     deltaZ += 0.5f * (**it).surface().bounds().thickness() * std::abs(tsos.globalDirection().z()) /
0080               tsos.globalDirection().perp();
0081     bool is_compatible = (z0 < deltaZ);
0082 
0083     if (is_compatible) {
0084       v_tsos.push_back(make_pair((const DetLayer*)(*it), tsos));
0085     }
0086   }
0087   for (auto it = nfpix.begin(); it != nfpix.end(); ++it) {
0088     TrajectoryStateOnSurface tsos = propagatorAlong.propagate(fts, (**it).specificSurface());
0089     if (!tsos.isValid())
0090       continue;
0091 
0092     auto r2 = tsos.localPosition().perp2();
0093     float deltaR = 0.5f * (**it).surface().bounds().thickness() * tsos.localDirection().perp() /
0094                    std::abs(tsos.localDirection().z());
0095     auto ri2 = std::max((**it).specificSurface().innerRadius() - deltaR, 0.f);
0096     ri2 *= ri2;
0097     auto ro2 = (**it).specificSurface().outerRadius() + deltaR;
0098     ro2 *= ro2;
0099     bool is_compatible = ((r2 > ri2) && (r2 < ro2));
0100 
0101     if (is_compatible) {
0102       v_tsos.push_back(make_pair((const DetLayer*)(*it), tsos));
0103     }
0104   }
0105   for (auto it = pfpix.begin(); it != pfpix.end(); ++it) {
0106     TrajectoryStateOnSurface tsos = propagatorAlong.propagate(fts, (**it).specificSurface());
0107     if (!tsos.isValid())
0108       continue;
0109 
0110     auto r2 = tsos.localPosition().perp2();
0111     float deltaR = 0.5f * (**it).surface().bounds().thickness() * tsos.localDirection().perp() /
0112                    std::abs(tsos.localDirection().z());
0113     auto ri2 = std::max((**it).specificSurface().innerRadius() - deltaR, 0.f);
0114     ri2 *= ri2;
0115     auto ro2 = (**it).specificSurface().outerRadius() + deltaR;
0116     ro2 *= ro2;
0117     bool is_compatible = ((r2 > ri2) && (r2 < ro2));
0118     if (is_compatible) {
0119       v_tsos.push_back(make_pair((const DetLayer*)(*it), tsos));
0120     }
0121   }
0122   return v_tsos;
0123 }
0124 
0125 // FIX HERE
0126 vector<pair<LayerHit, LayerTSOS> > SeedMvaEstimatorPhase2::getHitTsosPairs(
0127     const TrajectorySeed& seed,
0128     const edm::Handle<l1t::TrackerMuonCollection>& L1TkMuonHandle,
0129     const edm::ESHandle<MagneticField>& magfieldH,
0130     const Propagator& propagatorAlong,
0131     const GeometricSearchTracker& geomTracker) const {
0132   vector<pair<LayerHit, LayerTSOS> > out = {};
0133   float av_dr_min = 999.;
0134 
0135   // -- loop on L1TkMu
0136   for (auto L1TkMu = L1TkMuonHandle->begin(); L1TkMu != L1TkMuonHandle->end(); ++L1TkMu) {
0137     const vector<LayerTSOS> v_tsos = getTsosOnPixels(*L1TkMu->trkPtr(), magfieldH, propagatorAlong, geomTracker);
0138 
0139     // -- loop on recHits
0140     vector<int> v_tsos_skip(v_tsos.size(), 0);
0141     vector<pair<LayerHit, LayerTSOS> > hitTsosPair = {};
0142     int ihit = 0;
0143     for (const auto& hit : seed.recHits()) {
0144       // -- look for closest tsos by absolute distance
0145       int the_tsos = -99999;
0146       float dr_min = 20.;
0147       for (auto i = 0U; i < v_tsos.size(); ++i) {
0148         if (v_tsos_skip.at(i))
0149           continue;
0150         float dr = (v_tsos.at(i).second.globalPosition() - hit.globalPosition()).mag();
0151         if (dr < dr_min) {
0152           dr_min = dr;
0153           the_tsos = i;
0154           v_tsos_skip.at(i) = 1;
0155         }
0156       }
0157       if (the_tsos > -1) {
0158         const DetLayer* thelayer = geomTracker.idToLayer(hit.geographicalId());
0159         hitTsosPair.push_back(make_pair(make_pair(thelayer, &hit), v_tsos.at(the_tsos)));
0160       }
0161       ihit++;
0162     }  // loop on recHits
0163 
0164     // -- find tsos for all recHits?
0165     if ((int)hitTsosPair.size() == ihit) {
0166       float av_dr = 0.;
0167       for (auto it = hitTsosPair.begin(); it != hitTsosPair.end(); ++it) {
0168         auto hit = it->first.second;
0169         auto tsos = it->second.second;
0170         av_dr += (hit->globalPosition() - tsos.globalPosition()).mag();
0171       }
0172       av_dr = av_dr > 0 ? av_dr / (float)hitTsosPair.size() : 1.e6;
0173       if (av_dr < av_dr_min) {
0174         av_dr_min = av_dr;
0175         out = hitTsosPair;
0176       }
0177     }
0178   }  // loop on L1TkMu
0179   return out;
0180 }
0181 
0182 void SeedMvaEstimatorPhase2::getHitL1TkVariables(const TrajectorySeed& seed,
0183                                                  const edm::Handle<l1t::TrackerMuonCollection>& L1TkMuonHandle,
0184                                                  const edm::ESHandle<MagneticField>& magfieldH,
0185                                                  const Propagator& propagatorAlong,
0186                                                  const GeometricSearchTracker& geomTracker,
0187                                                  float& expd2HitL1Tk1,
0188                                                  float& expd2HitL1Tk2,
0189                                                  float& expd2HitL1Tk3) const {
0190   vector<pair<LayerHit, LayerTSOS> > hitTsosPair =
0191       getHitTsosPairs(seed, L1TkMuonHandle, magfieldH, propagatorAlong, geomTracker);
0192 
0193   bool found = (!hitTsosPair.empty());
0194 
0195   if (found) {
0196     float l1x1 = 99999.;
0197     float l1y1 = 99999.;
0198     float l1z1 = 99999.;
0199     float hitx1 = 99999.;
0200     float hity1 = 99999.;
0201     float hitz1 = 99999.;
0202 
0203     float l1x2 = 99999.;
0204     float l1y2 = 99999.;
0205     float l1z2 = 99999.;
0206     float hitx2 = 99999.;
0207     float hity2 = 99999.;
0208     float hitz2 = 99999.;
0209 
0210     float l1x3 = 99999.;
0211     float l1y3 = 99999.;
0212     float l1z3 = 99999.;
0213     float hitx3 = 99999.;
0214     float hity3 = 99999.;
0215     float hitz3 = 99999.;
0216 
0217     if (hitTsosPair.size() > 1) {
0218       auto hit1 = hitTsosPair.at(0).first.second;
0219       auto tsos1 = hitTsosPair.at(0).second.second;
0220 
0221       l1x1 = tsos1.globalPosition().x();
0222       l1y1 = tsos1.globalPosition().y();
0223       l1z1 = tsos1.globalPosition().z();
0224       hitx1 = hit1->globalPosition().x();
0225       hity1 = hit1->globalPosition().y();
0226       hitz1 = hit1->globalPosition().z();
0227 
0228       auto hit2 = hitTsosPair.at(1).first.second;
0229       auto tsos2 = hitTsosPair.at(1).second.second;
0230 
0231       l1x2 = tsos2.globalPosition().x();
0232       l1y2 = tsos2.globalPosition().y();
0233       l1z2 = tsos2.globalPosition().z();
0234       hitx2 = hit2->globalPosition().x();
0235       hity2 = hit2->globalPosition().y();
0236       hitz2 = hit2->globalPosition().z();
0237     }
0238     if (hitTsosPair.size() > 2) {
0239       auto hit3 = hitTsosPair.at(2).first.second;
0240       auto tsos3 = hitTsosPair.at(2).second.second;
0241 
0242       l1x3 = tsos3.globalPosition().x();
0243       l1y3 = tsos3.globalPosition().y();
0244       l1z3 = tsos3.globalPosition().z();
0245       hitx3 = hit3->globalPosition().x();
0246       hity3 = hit3->globalPosition().y();
0247       hitz3 = hit3->globalPosition().z();
0248     }
0249 
0250     // If hit == 99999. >> cannot find hit info >> distance btw hit~tsos is large >> expd2HitL1Tk ~ 0
0251     if (hitx1 != 99999.) {
0252       expd2HitL1Tk1 = exp(-(pow((l1x1 - hitx1), 2) + pow((l1y1 - hity1), 2) + pow((l1z1 - hitz1), 2)));
0253     } else {
0254       expd2HitL1Tk1 = 0.;
0255     }
0256 
0257     if (hitx2 != 99999.) {
0258       expd2HitL1Tk2 = exp(-(pow((l1x2 - hitx2), 2) + pow((l1y2 - hity2), 2) + pow((l1z2 - hitz2), 2)));
0259     } else {
0260       expd2HitL1Tk2 = 0.;
0261     }
0262 
0263     if (hitx3 != 99999.) {
0264       expd2HitL1Tk3 = exp(-(pow((l1x3 - hitx3), 2) + pow((l1y3 - hity3), 2) + pow((l1z3 - hitz3), 2)));
0265     } else {
0266       expd2HitL1Tk3 = 0.;
0267     }
0268   }
0269 }
0270 
0271 double SeedMvaEstimatorPhase2::computeMva(const TrajectorySeed& seed,
0272                                           const GlobalVector& global_p,
0273                                           const GlobalPoint& global_x,
0274                                           const edm::Handle<l1t::TrackerMuonCollection>& h_L1TkMu,
0275                                           const edm::ESHandle<MagneticField>& magfieldH,
0276                                           const Propagator& propagatorAlong,
0277                                           const GeometricSearchTracker& geomTracker) const {
0278   float Phase2var[Phase2::kLast]{};
0279 
0280   Phase2var[Phase2::kTsosErr0] = seed.startingState().error(0);
0281   Phase2var[Phase2::kTsosErr2] = seed.startingState().error(2);
0282   Phase2var[Phase2::kTsosErr5] = seed.startingState().error(5);
0283   Phase2var[Phase2::kTsosDxdz] = seed.startingState().parameters().dxdz();
0284   Phase2var[Phase2::kTsosDydz] = seed.startingState().parameters().dydz();
0285   Phase2var[Phase2::kTsosQbp] = seed.startingState().parameters().qbp();
0286 
0287   float initDRdPhi = 99999.;
0288 
0289   float DRL1TkMuSeedP = initDRdPhi;
0290   float DPhiL1TkMuSeedP = initDRdPhi;
0291   getL1TTVariables(seed, global_p, global_x, h_L1TkMu, DRL1TkMuSeedP, DPhiL1TkMuSeedP);
0292   Phase2var[Phase2::kDRdRL1TkMuSeedP] = DRL1TkMuSeedP;
0293   Phase2var[Phase2::kDRdPhiL1TkMuSeedP] = DPhiL1TkMuSeedP;
0294 
0295   float expd2HitL1Tk1 = initDRdPhi;
0296   float expd2HitL1Tk2 = initDRdPhi;
0297   float expd2HitL1Tk3 = initDRdPhi;
0298   getHitL1TkVariables(
0299       seed, h_L1TkMu, magfieldH, propagatorAlong, geomTracker, expd2HitL1Tk1, expd2HitL1Tk2, expd2HitL1Tk3);
0300   Phase2var[Phase2::kExpd2HitL1Tk1] = expd2HitL1Tk1;
0301   Phase2var[Phase2::kExpd2HitL1Tk2] = expd2HitL1Tk2;
0302   Phase2var[Phase2::kExpd2HitL1Tk3] = expd2HitL1Tk3;
0303 
0304   for (int iv = 0; iv < Phase2::kLast; ++iv) {
0305     Phase2var[iv] = (Phase2var[iv] - scale_mean_.at(iv)) / scale_std_.at(iv);
0306   }
0307 
0308   return gbrForest_->GetResponse(Phase2var);
0309 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!