METProducers/src/MuonTCMETValueMapProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    METProducers
0004 // Class:      MuonTCMETValueMapProducer
0005 //
0006 // Original Author:  Frank Golf
0007 //         Created:  Sun Mar 15 11:33:20 CDT 2009
0008 //
0009 //
0010
0011 //____________________________________________________________________________||
0012 #include "RecoMET/METProducers/interface/MuonTCMETValueMapProducer.h"
0013
0014 #include "RecoMET/METAlgorithms/interface/TCMETAlgo.h"
0015
0016 #include "DataFormats/MuonReco/interface/Muon.h"
0017 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019 #include "DataFormats/TrackReco/interface/TrackBase.h"
0020 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0021 #include "DataFormats/Common/interface/ValueMap.h"
0022 #include "DataFormats/MuonReco/interface/MuonMETCorrectionData.h"
0023 #include "DataFormats/GeometrySurface/interface/Plane.h"
0024 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0025 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0026 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0027 #include "DataFormats/Math/interface/Point3D.h"
0028
0029 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0030 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0031
0032 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034
0035 #include "TVector3.h"
0036 #include "TMath.h"
0037
0038 #include <algorithm>
0039
0040 //____________________________________________________________________________||
0041 typedef math::XYZPoint Point;
0042
0043 //____________________________________________________________________________||
0044 namespace cms {
0045
0046   MuonTCMETValueMapProducer::MuonTCMETValueMapProducer(const edm::ParameterSet& iConfig) {
0047     produces<edm::ValueMap<reco::MuonMETCorrectionData>>("muCorrData");
0048
0049     rfType_ = iConfig.getParameter<int>("rf_type");
0050
0051     nLayers_ = iConfig.getParameter<int>("nLayers");
0052     nLayersTight_ = iConfig.getParameter<int>("nLayersTight");
0053     vertexNdof_ = iConfig.getParameter<int>("vertexNdof");
0054     vertexZ_ = iConfig.getParameter<double>("vertexZ");
0055     vertexRho_ = iConfig.getParameter<double>("vertexRho");
0056     vertexMaxDZ_ = iConfig.getParameter<double>("vertexMaxDZ");
0057     maxpt_eta20_ = iConfig.getParameter<double>("maxpt_eta20");
0058     maxpt_eta25_ = iConfig.getParameter<double>("maxpt_eta25");
0059
0060     // get configuration parameters
0061     std::vector<std::string> algos = iConfig.getParameter<std::vector<std::string>>("trackAlgos");
0062     std::transform(algos.begin(), algos.end(), std::back_inserter(trkAlgos_), [](const std::string& a) {
0063       return reco::TrackBase::algoByName(a);
0064     });
0065     maxd0cut_ = iConfig.getParameter<double>("d0_max");
0066     minpt_ = iConfig.getParameter<double>("pt_min");
0067     maxpt_ = iConfig.getParameter<double>("pt_max");
0068     maxeta_ = iConfig.getParameter<double>("eta_max");
0069     maxchi2_ = iConfig.getParameter<double>("chi2_max");
0070     minhits_ = iConfig.getParameter<double>("nhits_min");
0071     maxPtErr_ = iConfig.getParameter<double>("ptErr_max");
0072
0073     trkQuality_ = iConfig.getParameter<std::vector<int>>("track_quality");
0074     algos = iConfig.getParameter<std::vector<std::string>>("track_algos");
0075     std::transform(algos.begin(), algos.end(), std::back_inserter(trkAlgos_), [](const std::string& a) {
0076       return reco::TrackBase::algoByName(a);
0077     });
0078     maxchi2_tight_ = iConfig.getParameter<double>("chi2_max_tight");
0079     minhits_tight_ = iConfig.getParameter<double>("nhits_min_tight");
0080     maxPtErr_tight_ = iConfig.getParameter<double>("ptErr_max_tight");
0081     usePvtxd0_ = iConfig.getParameter<bool>("usePvtxd0");
0082     d0cuta_ = iConfig.getParameter<double>("d0cuta");
0083     d0cutb_ = iConfig.getParameter<double>("d0cutb");
0084
0085     muon_dptrel_ = iConfig.getParameter<double>("muon_dptrel");
0086     muond0_ = iConfig.getParameter<double>("d0_muon");
0087     muonpt_ = iConfig.getParameter<double>("pt_muon");
0088     muoneta_ = iConfig.getParameter<double>("eta_muon");
0089     muonchi2_ = iConfig.getParameter<double>("chi2_muon");
0090     muonhits_ = iConfig.getParameter<double>("nhits_muon");
0091     muonGlobal_ = iConfig.getParameter<bool>("global_muon");
0092     muonTracker_ = iConfig.getParameter<bool>("tracker_muon");
0093     muonDeltaR_ = iConfig.getParameter<double>("deltaR_muon");
0094     useCaloMuons_ = iConfig.getParameter<bool>("useCaloMuons");
0095     muonMinValidStaHits_ = iConfig.getParameter<int>("muonMinValidStaHits");
0096
0097     response_function = nullptr;
0098     tcmetAlgo_ = new TCMETAlgo();
0099
0100     if (rfType_ == 1)
0101       response_function = tcmetAlgo_->getResponseFunction_fit();
0102     else if (rfType_ == 2)
0103       response_function = tcmetAlgo_->getResponseFunction_mode();
0104
0105     muonToken_ = consumes<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("muonInputTag"));
0106     beamSpotToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotInputTag"));
0107     vertexToken_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexInputTag"));
0108     magFieldToken_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0109   }
0110
0111   //____________________________________________________________________________||
0112   MuonTCMETValueMapProducer::~MuonTCMETValueMapProducer() { delete tcmetAlgo_; }
0113
0114   //____________________________________________________________________________||
0115   void MuonTCMETValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0116     iEvent.getByToken(muonToken_, muons_);
0117     iEvent.getByToken(beamSpotToken_, beamSpot_);
0118
0119     hasValidVertex = false;
0120     if (usePvtxd0_) {
0121       iEvent.getByToken(vertexToken_, vertexHandle_);
0122
0123       if (vertexHandle_.isValid()) {
0124         vertices_ = vertexHandle_.product();
0125         hasValidVertex = isValidVertex();
0126       }
0127     }
0128
0129     edm::ESHandle<MagneticField> theMagField = iSetup.getHandle(magFieldToken_);
0130     bField = theMagField.product();
0131
0132     auto vm_muCorrData = std::make_unique<edm::ValueMap<reco::MuonMETCorrectionData>>();
0133
0134     std::vector<reco::MuonMETCorrectionData> v_muCorrData;
0135
0136     unsigned int nMuons = muons_->size();
0137
0138     for (unsigned int iMu = 0; iMu < nMuons; iMu++) {
0139       const reco::Muon* mu = &(*muons_)[iMu];
0140       double deltax = 0.0;
0141       double deltay = 0.0;
0142
0143       reco::MuonMETCorrectionData muMETCorrData(reco::MuonMETCorrectionData::NotUsed, deltax, deltay);
0144
0145       reco::TrackRef mu_track;
0146       if (mu->isGlobalMuon() || mu->isTrackerMuon() || mu->isCaloMuon())
0147         mu_track = mu->innerTrack();
0148       else {
0149         v_muCorrData.push_back(muMETCorrData);
0150         continue;
0151       }
0152
0153       // figure out depositions muons would make if they were treated as pions
0154       if (isGoodTrack(mu)) {
0155         if (mu_track->pt() < minpt_)
0156           muMETCorrData = reco::MuonMETCorrectionData(reco::MuonMETCorrectionData::TreatedAsPion, deltax, deltay);
0157
0158         else {
0159           int bin_index = response_function->FindBin(mu_track->eta(), mu_track->pt());
0160           double response = response_function->GetBinContent(bin_index);
0161
0162           TVector3 outerTrkPosition = propagateTrack(mu);
0163
0164           deltax = response * mu_track->p() * sin(outerTrkPosition.Theta()) * cos(outerTrkPosition.Phi());
0165           deltay = response * mu_track->p() * sin(outerTrkPosition.Theta()) * sin(outerTrkPosition.Phi());
0166
0167           muMETCorrData = reco::MuonMETCorrectionData(reco::MuonMETCorrectionData::TreatedAsPion, deltax, deltay);
0168         }
0169       }
0170
0171       // figure out muon flag
0172       if (isGoodMuon(mu))
0173         v_muCorrData.push_back(
0174             reco::MuonMETCorrectionData(reco::MuonMETCorrectionData::MuonCandidateValuesUsed, deltax, deltay));
0175
0176       else if (useCaloMuons_ && isGoodCaloMuon(mu, iMu))
0177         v_muCorrData.push_back(
0178             reco::MuonMETCorrectionData(reco::MuonMETCorrectionData::MuonCandidateValuesUsed, deltax, deltay));
0179
0180       else
0181         v_muCorrData.push_back(muMETCorrData);
0182     }
0183
0184     edm::ValueMap<reco::MuonMETCorrectionData>::Filler dataFiller(*vm_muCorrData);
0185
0186     dataFiller.insert(muons_, v_muCorrData.begin(), v_muCorrData.end());
0187     dataFiller.fill();
0188
0189     iEvent.put(std::move(vm_muCorrData), "muCorrData");
0190   }
0191
0192   //____________________________________________________________________________||
0193   bool MuonTCMETValueMapProducer::isGoodMuon(const reco::Muon* muon) {
0194     double d0 = -999;
0195     double nhits = 0;
0196     double chi2 = 999;
0197
0198     // get d0 corrected for beam spot
0199     bool haveBeamSpot = true;
0200     if (!beamSpot_.isValid())
0201       haveBeamSpot = false;
0202
0203     if (muonGlobal_ && !muon->isGlobalMuon())
0204       return false;
0205     if (muonTracker_ && !muon->isTrackerMuon())
0206       return false;
0207
0208     const reco::TrackRef siTrack = muon->innerTrack();
0209     const reco::TrackRef globalTrack = muon->globalTrack();
0210
0211     Point bspot = haveBeamSpot ? beamSpot_->position() : Point(0, 0, 0);
0212     if (siTrack.isNonnull())
0213       nhits = siTrack->numberOfValidHits();
0214     if (globalTrack.isNonnull()) {
0215       d0 = -1 * globalTrack->dxy(bspot);
0216       chi2 = globalTrack->normalizedChi2();
0217     }
0218
0219     if (fabs(d0) > muond0_)
0220       return false;
0221     if (muon->pt() < muonpt_)
0222       return false;
0223     if (fabs(muon->eta()) > muoneta_)
0224       return false;
0225     if (nhits < muonhits_)
0226       return false;
0227     if (chi2 > muonchi2_)
0228       return false;
0229     if (globalTrack->hitPattern().numberOfValidMuonHits() < muonMinValidStaHits_)
0230       return false;
0231
0232     //reject muons with tracker dpt/pt > X
0233     if (!siTrack.isNonnull())
0234       return false;
0235     if (siTrack->ptError() / siTrack->pt() > muon_dptrel_)
0236       return false;
0237
0238     else
0239       return true;
0240   }
0241
0242   //____________________________________________________________________________||
0243   bool MuonTCMETValueMapProducer::isGoodCaloMuon(const reco::Muon* muon, const unsigned int index) {
0244     if (muon->pt() < 10)
0245       return false;
0246
0247     if (!isGoodTrack(muon))
0248       return false;
0249
0250     const reco::TrackRef inputSiliconTrack = muon->innerTrack();
0251     if (!inputSiliconTrack.isNonnull())
0252       return false;
0253
0254     //check if it is in the vicinity of a global or tracker muon
0255     unsigned int nMuons = muons_->size();
0256     for (unsigned int iMu = 0; iMu < nMuons; iMu++) {
0257       if (iMu == index)
0258         continue;
0259
0260       const reco::Muon* mu = &(*muons_)[iMu];
0261
0262       const reco::TrackRef testSiliconTrack = mu->innerTrack();
0263       if (!testSiliconTrack.isNonnull())
0264         continue;
0265
0266       double deltaEta = inputSiliconTrack.get()->eta() - testSiliconTrack.get()->eta();
0267       double deltaPhi = acos(cos(inputSiliconTrack.get()->phi() - testSiliconTrack.get()->phi()));
0268       double deltaR = TMath::Sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
0269
0270       if (deltaR < muonDeltaR_)
0271         return false;
0272     }
0273
0274     return true;
0275   }
0276
0277   //____________________________________________________________________________||
0278   bool MuonTCMETValueMapProducer::isGoodTrack(const reco::Muon* muon) {
0279     double d0;
0280
0281     const reco::TrackRef siTrack = muon->innerTrack();
0282     if (!siTrack.isNonnull())
0283       return false;
0284
0285     if (hasValidVertex) {
0286       //get d0 corrected for primary vertex
0287
0288       const Point pvtx = Point(vertices_->begin()->x(), vertices_->begin()->y(), vertices_->begin()->z());
0289
0290       double dz = siTrack->dz(pvtx);
0291
0292       if (fabs(dz) < vertexMaxDZ_) {
0293         //get d0 corrected for pvtx
0294         d0 = -1 * siTrack->dxy(pvtx);
0295
0296       } else {
0297         // get d0 corrected for beam spot
0298         bool haveBeamSpot = true;
0299         if (!beamSpot_.isValid())
0300           haveBeamSpot = false;
0301
0302         Point bspot = haveBeamSpot ? beamSpot_->position() : Point(0, 0, 0);
0303         d0 = -1 * siTrack->dxy(bspot);
0304       }
0305     } else {
0306       // get d0 corrected for beam spot
0307       bool haveBeamSpot = true;
0308       if (!beamSpot_.isValid())
0309         haveBeamSpot = false;
0310
0311       Point bspot = haveBeamSpot ? beamSpot_->position() : Point(0, 0, 0);
0312       d0 = -1 * siTrack->dxy(bspot);
0313     }
0314
0315     if (std::find(trackAlgos_.begin(), trackAlgos_.end(), siTrack->algo()) != trackAlgos_.end()) {
0316       //1st 4 tracking iterations (pT-dependent d0 cut)
0317
0318       float d0cut = sqrt(std::pow(d0cuta_, 2) + std::pow(d0cutb_ / siTrack->pt(), 2));
0319       if (d0cut > maxd0cut_)
0320         d0cut = maxd0cut_;
0321
0322       if (fabs(d0) > d0cut)
0323         return false;
0324       if (nLayers(siTrack) < nLayers_)
0325         return false;
0326     } else {
0327       //last 2 tracking iterations (tighten chi2, nhits, pt error cuts)
0328
0329       if (siTrack->normalizedChi2() > maxchi2_tight_)
0330         return false;
0331       if (siTrack->numberOfValidHits() < minhits_tight_)
0332         return false;
0333       if ((siTrack->ptError() / siTrack->pt()) > maxPtErr_tight_)
0334         return false;
0335       if (nLayers(siTrack) < nLayersTight_)
0336         return false;
0337     }
0338
0339     if (siTrack->numberOfValidHits() < minhits_)
0340       return false;
0341     if (siTrack->normalizedChi2() > maxchi2_)
0342       return false;
0343     if (fabs(siTrack->eta()) > maxeta_)
0344       return false;
0345     if (siTrack->pt() > maxpt_)
0346       return false;
0347     if ((siTrack->ptError() / siTrack->pt()) > maxPtErr_)
0348       return false;
0349     if (fabs(siTrack->eta()) > 2.5 && siTrack->pt() > maxpt_eta25_)
0350       return false;
0351     if (fabs(siTrack->eta()) > 2.0 && siTrack->pt() > maxpt_eta20_)
0352       return false;
0353
0354     int cut = 0;
0355     for (unsigned int i = 0; i < trkQuality_.size(); i++) {
0356       cut |= (1 << trkQuality_.at(i));
0357     }
0358
0359     if (!((siTrack->qualityMask() & cut) == cut))
0360       return false;
0361
0362     bool isGoodAlgo = false;
0363     if (trkAlgos_.empty())
0364       isGoodAlgo = true;
0365     for (unsigned int i = 0; i < trkAlgos_.size(); i++) {
0366       if (siTrack->algo() == trkAlgos_.at(i))
0367         isGoodAlgo = true;
0368     }
0369
0370     if (!isGoodAlgo)
0371       return false;
0372
0373     return true;
0374   }
0375
0376   //____________________________________________________________________________||
0377   TVector3 MuonTCMETValueMapProducer::propagateTrack(const reco::Muon* muon) {
0378     TVector3 outerTrkPosition;
0379
0380     outerTrkPosition.SetPtEtaPhi(999., -10., 2 * TMath::Pi());
0381
0382     const reco::TrackRef track = muon->innerTrack();
0383
0384     if (!track.isNonnull()) {
0385       return outerTrkPosition;
0386     }
0387
0388     GlobalPoint tpVertex(track->vx(), track->vy(), track->vz());
0389     GlobalVector tpMomentum(track.get()->px(), track.get()->py(), track.get()->pz());
0390     int tpCharge(track->charge());
0391
0392     FreeTrajectoryState fts(tpVertex, tpMomentum, tpCharge, bField);
0393
0394     const double zdist = 314.;
0395
0396     const double radius = 130.;
0397
0398     const double corner = 1.479;
0399
0400     Plane::PlanePointer lendcap = Plane::build(Plane::PositionType(0, 0, -zdist), Plane::RotationType());
0401     Plane::PlanePointer rendcap = Plane::build(Plane::PositionType(0, 0, zdist), Plane::RotationType());
0402
0403     Cylinder::CylinderPointer barrel =
0404         Cylinder::build(Cylinder::PositionType(0, 0, 0), Cylinder::RotationType(), radius);
0405
0406     AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
0407
0408     TrajectoryStateOnSurface tsos;
0409
0410     if (track.get()->eta() < -corner) {
0411       tsos = myAP.propagate(fts, *lendcap);
0412     } else if (fabs(track.get()->eta()) < corner) {
0413       tsos = myAP.propagate(fts, *barrel);
0414     } else if (track.get()->eta() > corner) {
0415       tsos = myAP.propagate(fts, *rendcap);
0416     }
0417
0418     if (tsos.isValid())
0419       outerTrkPosition.SetXYZ(tsos.globalPosition().x(), tsos.globalPosition().y(), tsos.globalPosition().z());
0420
0421     else
0422       outerTrkPosition.SetPtEtaPhi(999., -10., 2 * TMath::Pi());
0423
0424     return outerTrkPosition;
0425   }
0426
0427   //____________________________________________________________________________||
0428   int MuonTCMETValueMapProducer::nLayers(const reco::TrackRef track) {
0429     const reco::HitPattern& p = track->hitPattern();
0430     return p.trackerLayersWithMeasurement();
0431   }
0432
0433   //____________________________________________________________________________||
0434   bool MuonTCMETValueMapProducer::isValidVertex() {
0435     if (vertices_->begin()->isFake())
0436       return false;
0437     if (vertices_->begin()->ndof() < vertexNdof_)
0438       return false;
0439     if (fabs(vertices_->begin()->z()) > vertexZ_)
0440       return false;
0441     if (sqrt(std::pow(vertices_->begin()->x(), 2) + std::pow(vertices_->begin()->y(), 2)) > vertexRho_)
0442       return false;
0443
0444     return true;
0445   }
0446
0447 }  // namespace cms
0448
0449 //____________________________________________________________________________||