GlobalTrackingTools/plugins/GlobalTrackQualityProducer.cc

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