File indexing completed on 2021-02-14 23:30:49
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #include <memory>
0013
0014
0015 #include "FWCore/Framework/interface/ConsumesCollector.h"
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/EDProducer.h"
0018 #include "FWCore/Utilities/interface/isFinite.h"
0019
0020 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0021 #include "RecoMuon/GlobalTrackingTools/plugins/GlobalTrackQualityProducer.h"
0022
0023 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0024 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0025 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0026 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0027
0028 GlobalTrackQualityProducer::GlobalTrackQualityProducer(const edm::ParameterSet& iConfig)
0029 : inputCollection_(iConfig.getParameter<edm::InputTag>("InputCollection")),
0030 inputLinksCollection_(iConfig.getParameter<edm::InputTag>("InputLinksCollection")),
0031 tTopoToken_(esConsumes()),
0032 theService(nullptr),
0033 theGlbRefitter(nullptr),
0034 theGlbMatcher(nullptr) {
0035
0036 edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
0037 theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0038
0039
0040 edm::ConsumesCollector iC = consumesCollector();
0041 edm::ParameterSet refitterParameters = iConfig.getParameter<edm::ParameterSet>("RefitterParameters");
0042 theGlbRefitter = new GlobalMuonRefitter(refitterParameters, theService, iC);
0043
0044 edm::ParameterSet trackMatcherPSet = iConfig.getParameter<edm::ParameterSet>("GlobalMuonTrackMatcher");
0045 theGlbMatcher = new GlobalMuonTrackMatcher(trackMatcherPSet, theService);
0046
0047 double maxChi2 = iConfig.getParameter<double>("MaxChi2");
0048 double nSigma = iConfig.getParameter<double>("nSigma");
0049 theEstimator = new Chi2MeasurementEstimator(maxChi2, nSigma);
0050
0051 glbMuonsToken = consumes<reco::TrackCollection>(inputCollection_);
0052 linkCollectionToken = consumes<reco::MuonTrackLinksCollection>(inputLinksCollection_);
0053
0054 produces<edm::ValueMap<reco::MuonQuality>>();
0055 }
0056
0057 GlobalTrackQualityProducer::~GlobalTrackQualityProducer() {
0058 if (theService)
0059 delete theService;
0060 if (theGlbRefitter)
0061 delete theGlbRefitter;
0062 if (theGlbMatcher)
0063 delete theGlbMatcher;
0064 if (theEstimator)
0065 delete theEstimator;
0066 }
0067
0068 void GlobalTrackQualityProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0069 const std::string theCategory = "Muon|RecoMuon|GlobalTrackQualityProducer";
0070
0071 theService->update(iSetup);
0072
0073 theGlbRefitter->setEvent(iEvent);
0074
0075 theGlbRefitter->setServices(theService->eventSetup());
0076
0077
0078 edm::Handle<reco::TrackCollection> glbMuons;
0079 iEvent.getByToken(glbMuonsToken, glbMuons);
0080
0081 edm::Handle<reco::MuonTrackLinksCollection> linkCollectionHandle;
0082 iEvent.getByToken(linkCollectionToken, linkCollectionHandle);
0083
0084
0085 const TrackerTopology* tTopo = &iSetup.getData(tTopoToken_);
0086
0087
0088 std::vector<reco::MuonQuality> valuesQual;
0089 valuesQual.reserve(glbMuons->size());
0090
0091 int trackIndex = 0;
0092 for (reco::TrackCollection::const_iterator track = glbMuons->begin(); track != glbMuons->end();
0093 ++track, ++trackIndex) {
0094 reco::TrackRef glbRef(glbMuons, trackIndex);
0095 reco::TrackRef staTrack = reco::TrackRef();
0096
0097 std::vector<Trajectory> refitted = theGlbRefitter->refit(*track, 1, tTopo);
0098
0099 LogTrace(theCategory) << "GLBQual N refitted " << refitted.size();
0100
0101 std::pair<double, double> thisKink;
0102 double relative_muon_chi2 = 0.0;
0103 double relative_tracker_chi2 = 0.0;
0104 double glbTrackProbability = 0.0;
0105 if (!refitted.empty()) {
0106 thisKink = kink(refitted.front());
0107 std::pair<double, double> chi = newChi2(refitted.front());
0108 relative_muon_chi2 = chi.second;
0109 relative_tracker_chi2 = chi.first;
0110 glbTrackProbability = trackProbability(refitted.front());
0111 }
0112
0113 LogTrace(theCategory) << "GLBQual: Kink " << thisKink.first << " " << thisKink.second;
0114 LogTrace(theCategory) << "GLBQual: Rel Chi2 " << relative_tracker_chi2 << " " << relative_muon_chi2;
0115 LogTrace(theCategory) << "GLBQual: trackProbability " << glbTrackProbability;
0116
0117
0118 float chi2, d, dist, Rpos;
0119 chi2 = d = dist = Rpos = -1.0;
0120 bool passTight = false;
0121 typedef MuonTrajectoryBuilder::TrackCand TrackCand;
0122 if (linkCollectionHandle.isValid()) {
0123 for (reco::MuonTrackLinksCollection::const_iterator links = linkCollectionHandle->begin();
0124 links != linkCollectionHandle->end();
0125 ++links) {
0126 if (links->trackerTrack().isNull() || links->standAloneTrack().isNull() || links->globalTrack().isNull()) {
0127 edm::LogWarning(theCategory) << "Global muon links to constituent tracks are invalid. There should be no "
0128 "such object. Muon is skipped.";
0129 continue;
0130 }
0131 if (links->globalTrack() == glbRef) {
0132 staTrack = !links->standAloneTrack().isNull() ? links->standAloneTrack() : reco::TrackRef();
0133 TrackCand staCand = TrackCand((Trajectory*)nullptr, links->standAloneTrack());
0134 TrackCand tkCand = TrackCand((Trajectory*)nullptr, links->trackerTrack());
0135 chi2 = theGlbMatcher->match(staCand, tkCand, 0, 0);
0136 d = theGlbMatcher->match(staCand, tkCand, 1, 0);
0137 Rpos = theGlbMatcher->match(staCand, tkCand, 2, 0);
0138 dist = theGlbMatcher->match(staCand, tkCand, 3, 0);
0139 passTight = theGlbMatcher->matchTight(staCand, tkCand);
0140 }
0141 }
0142 }
0143
0144 if (!staTrack.isNull())
0145 LogTrace(theCategory) << "GLBQual: Used UpdatedAtVtx : "
0146 << (iEvent.getStableProvenance(staTrack.id()).productInstanceName() ==
0147 std::string("UpdatedAtVtx"));
0148
0149 float maxFloat01 =
0150 std::numeric_limits<float>::max() * 0.1;
0151 reco::MuonQuality muQual;
0152 if (!staTrack.isNull())
0153 muQual.updatedSta =
0154 iEvent.getStableProvenance(staTrack.id()).productInstanceName() == std::string("UpdatedAtVtx");
0155 muQual.trkKink = thisKink.first > maxFloat01 ? maxFloat01 : thisKink.first;
0156 muQual.glbKink = thisKink.second > maxFloat01 ? maxFloat01 : thisKink.second;
0157 muQual.trkRelChi2 = relative_tracker_chi2 > maxFloat01 ? maxFloat01 : relative_tracker_chi2;
0158 muQual.staRelChi2 = relative_muon_chi2 > maxFloat01 ? maxFloat01 : relative_muon_chi2;
0159 muQual.tightMatch = passTight;
0160 muQual.chi2LocalPosition = dist;
0161 muQual.chi2LocalMomentum = chi2;
0162 muQual.localDistance = d;
0163 muQual.globalDeltaEtaPhi = Rpos;
0164 muQual.glbTrackProbability = glbTrackProbability;
0165 valuesQual.push_back(muQual);
0166 }
0167
0168
0169
0170
0171
0172
0173
0174
0175 auto outQual = std::make_unique<edm::ValueMap<reco::MuonQuality>>();
0176 edm::ValueMap<reco::MuonQuality>::Filler fillerQual(*outQual);
0177 fillerQual.insert(glbMuons, valuesQual.begin(), valuesQual.end());
0178 fillerQual.fill();
0179
0180
0181 iEvent.put(std::move(outQual));
0182 }
0183
0184 std::pair<double, double> GlobalTrackQualityProducer::kink(Trajectory& muon) const {
0185 const std::string theCategory = "Muon|RecoMuon|GlobalTrackQualityProducer";
0186
0187 using namespace std;
0188 using namespace edm;
0189 using namespace reco;
0190
0191 double result = 0.0;
0192 double resultGlb = 0.0;
0193
0194 typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
0195 typedef ConstRecHitPointer RecHit;
0196 typedef std::vector<TrajectoryMeasurement>::const_iterator TMI;
0197
0198 vector<TrajectoryMeasurement> meas = muon.measurements();
0199
0200 for (TMI m = meas.begin(); m != meas.end(); m++) {
0201 TransientTrackingRecHit::ConstRecHitPointer hit = m->recHit();
0202
0203
0204
0205 RecHit rhit = (*m).recHit();
0206 bool ok = false;
0207 if (rhit->isValid()) {
0208 if (DetId::Tracker == rhit->geographicalId().det())
0209 ok = true;
0210 }
0211
0212
0213
0214 const TrajectoryStateOnSurface& tsos = (*m).predictedState();
0215
0216 if (tsos.isValid() && rhit->isValid() && rhit->hit()->isValid() &&
0217 !edm::isNotFinite(rhit->localPositionError().xx())
0218 && !edm::isNotFinite(rhit->localPositionError().xy())
0219 && !edm::isNotFinite(rhit->localPositionError().yy())) {
0220 double phi1 = tsos.globalPosition().phi();
0221 if (phi1 < 0)
0222 phi1 = 2 * M_PI + phi1;
0223
0224 double phi2 = rhit->globalPosition().phi();
0225 if (phi2 < 0)
0226 phi2 = 2 * M_PI + phi2;
0227
0228 double diff = fabs(phi1 - phi2);
0229 if (diff > M_PI)
0230 diff = 2 * M_PI - diff;
0231
0232 GlobalPoint hitPos = rhit->globalPosition();
0233
0234 GlobalError hitErr = rhit->globalPositionError();
0235
0236 double error = hitErr.phierr(hitPos);
0237
0238 double s = (error > 0.0) ? (diff * diff) / error : (diff * diff);
0239
0240 if (ok)
0241 result += s;
0242 resultGlb += s;
0243 }
0244 }
0245
0246 return std::pair<double, double>(result, resultGlb);
0247 }
0248
0249 std::pair<double, double> GlobalTrackQualityProducer::newChi2(Trajectory& muon) const {
0250 const std::string theCategory = "Muon|RecoMuon|GlobalTrackQualityProducer";
0251
0252 using namespace std;
0253 using namespace edm;
0254 using namespace reco;
0255
0256 double muChi2 = 0.0;
0257 double tkChi2 = 0.0;
0258 unsigned int muNdof = 0;
0259 unsigned int tkNdof = 0;
0260
0261 typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
0262 typedef ConstRecHitPointer RecHit;
0263 typedef vector<TrajectoryMeasurement>::const_iterator TMI;
0264
0265 vector<TrajectoryMeasurement> meas = muon.measurements();
0266
0267 for (TMI m = meas.begin(); m != meas.end(); m++) {
0268 TransientTrackingRecHit::ConstRecHitPointer hit = m->recHit();
0269 const TrajectoryStateOnSurface& uptsos = (*m).updatedState();
0270
0271
0272 const auto& preciseHit = hit;
0273 double estimate = 0.0;
0274 if (preciseHit->isValid() && uptsos.isValid()) {
0275 estimate = theEstimator->estimate(uptsos, *preciseHit).second;
0276 }
0277
0278
0279
0280
0281 if (hit->isValid() && (hit->geographicalId().det()) == DetId::Tracker) {
0282 tkChi2 += estimate;
0283
0284 tkNdof += hit->dimension();
0285 }
0286 if (hit->isValid() && (hit->geographicalId().det()) == DetId::Muon) {
0287 muChi2 += estimate;
0288
0289 muNdof += hit->dimension();
0290 }
0291 }
0292
0293
0294
0295 if (tkNdof > 5) {
0296 tkChi2 /= (tkNdof - 5.);
0297 }
0298
0299
0300
0301 if (muNdof > 5) {
0302 muChi2 /= (muNdof - 5.);
0303 }
0304
0305 return std::pair<double, double>(tkChi2, muChi2);
0306 }
0307
0308
0309
0310
0311 double GlobalTrackQualityProducer::trackProbability(Trajectory& track) const {
0312 if (track.ndof() > 0 && track.chiSquared() > 0) {
0313 return -LnChiSquaredProbability(track.chiSquared(), track.ndof());
0314 } else {
0315 return 0.0;
0316 }
0317 }
0318
0319 void GlobalTrackQualityProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0320 edm::ParameterSetDescription desc;
0321 {
0322 edm::ParameterSetDescription psd1;
0323 psd1.setAllowAnything();
0324 desc.add<edm::ParameterSetDescription>("ServiceParameters", psd1);
0325 }
0326 {
0327 edm::ParameterSetDescription psd1;
0328 psd1.setAllowAnything();
0329 desc.add<edm::ParameterSetDescription>("GlobalMuonTrackMatcher", psd1);
0330 }
0331 desc.add<edm::InputTag>("InputCollection", edm::InputTag("globalMuons"));
0332 desc.add<edm::InputTag>("InputLinksCollection", edm::InputTag("globalMuons"));
0333 desc.add<std::string>("BaseLabel", "GLB");
0334 {
0335 edm::ParameterSetDescription descGlbMuonRefitter;
0336 descGlbMuonRefitter.setAllowAnything();
0337 descGlbMuonRefitter.add<edm::InputTag>("DTRecSegmentLabel", edm::InputTag("dt1DRecHits"));
0338 descGlbMuonRefitter.add<edm::InputTag>("CSCRecSegmentLabel", edm::InputTag("csc2DRecHits"));
0339 descGlbMuonRefitter.add<edm::InputTag>("GEMRecHitLabel", edm::InputTag("gemRecHits"));
0340 descGlbMuonRefitter.add<edm::InputTag>("ME0RecHitLabel", edm::InputTag("me0Segments"));
0341 descGlbMuonRefitter.add<edm::InputTag>("RPCRecSegmentLabel", edm::InputTag("rpcRecHits"));
0342
0343 descGlbMuonRefitter.add<std::string>("Fitter", "KFFitterForRefitInsideOut");
0344 descGlbMuonRefitter.add<std::string>("Smoother", "KFSmootherForRefitInsideOut");
0345 descGlbMuonRefitter.add<std::string>("Propagator", "SmartPropagatorAnyRK");
0346 descGlbMuonRefitter.add<std::string>("TrackerRecHitBuilder", "WithAngleAndTemplate");
0347 descGlbMuonRefitter.add<std::string>("MuonRecHitBuilder", "MuonRecHitBuilder");
0348 descGlbMuonRefitter.add<bool>("DoPredictionsOnly", false);
0349 descGlbMuonRefitter.add<std::string>("RefitDirection", "insideOut");
0350 descGlbMuonRefitter.add<bool>("PropDirForCosmics", false);
0351 descGlbMuonRefitter.add<bool>("RefitRPCHits", true);
0352
0353 descGlbMuonRefitter.add<std::vector<int>>("DYTthrs", {10, 10});
0354 descGlbMuonRefitter.add<int>("DYTselector", 1);
0355 descGlbMuonRefitter.add<bool>("DYTupdator", false);
0356 descGlbMuonRefitter.add<bool>("DYTuseAPE", false);
0357 descGlbMuonRefitter.add<bool>("DYTuseThrsParametrization", true);
0358 {
0359 edm::ParameterSetDescription descDYTthrs;
0360 descDYTthrs.add<std::vector<double>>("eta0p8", {1, -0.919853, 0.990742});
0361 descDYTthrs.add<std::vector<double>>("eta1p2", {1, -0.897354, 0.987738});
0362 descDYTthrs.add<std::vector<double>>("eta2p0", {4, -0.986855, 0.998516});
0363 descDYTthrs.add<std::vector<double>>("eta2p2", {1, -0.940342, 0.992955});
0364 descDYTthrs.add<std::vector<double>>("eta2p4", {1, -0.947633, 0.993762});
0365 descGlbMuonRefitter.add<edm::ParameterSetDescription>("DYTthrsParameters", descDYTthrs);
0366 }
0367
0368 descGlbMuonRefitter.add<int>("SkipStation", -1);
0369 descGlbMuonRefitter.add<int>("TrackerSkipSystem", -1);
0370 descGlbMuonRefitter.add<int>("TrackerSkipSection", -1);
0371 descGlbMuonRefitter.add<bool>("RefitFlag", true);
0372
0373 desc.add<edm::ParameterSetDescription>("RefitterParameters", descGlbMuonRefitter);
0374 }
0375 desc.add<double>("nSigma", 3.0);
0376 desc.add<double>("MaxChi2", 100000.0);
0377
0378 descriptions.add("globalTrackQualityProducer", desc);
0379 }
0380
0381