Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 23:30:49 

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/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   // service parameters
0036   edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
0037   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0038 
0039   // TrackRefitter parameters
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   // Take the GLB muon container(s)
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   //Retrieve tracker topology from geometry
0085   const TrackerTopology* tTopo = &iSetup.getData(tTopoToken_);
0086 
0087   // reserve some space
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;    //normalized inside to /sum(muHits.dimension)
0109       relative_tracker_chi2 = chi.first;  //normalized inside to /sum(tkHits.dimension)
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     // Fill the STA-TK match information
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;  // a better solution would be to use float above .. m/be not
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   for(int i = 0; i < valuesTkRelChi2.size(); i++) {
0170     LogTrace(theCategory)<<"value " << valuesTkRelChi2[i] ;
0171   }
0172   */
0173 
0174   // create and fill value maps
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   // put value map into event
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     //not used    double estimate = 0.0;
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     //if ( !ok ) continue;
0213 
0214     const TrajectoryStateOnSurface& tsos = (*m).predictedState();
0215 
0216     if (tsos.isValid() && rhit->isValid() && rhit->hit()->isValid() &&
0217         !edm::isNotFinite(rhit->localPositionError().xx())     //this is paranoia induced by reported case
0218         && !edm::isNotFinite(rhit->localPositionError().xy())  //it's better to track down the origin of bad numbers
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       //LogDebug(theCategory)<<"hitPos " << hitPos;
0236       double error = hitErr.phierr(hitPos);  // error squared
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     // FIXME FIXME CLONE!!!
0271     // TrackingRecHit::RecHitPointer preciseHit = hit->clone(uptsos);
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     //LogTrace(theCategory) << "estimate " << estimate << " TM.est " << m->estimate();
0279     //UNUSED:    double tkDiff = 0.0;
0280     //UNUSED:    double staDiff = 0.0;
0281     if (hit->isValid() && (hit->geographicalId().det()) == DetId::Tracker) {
0282       tkChi2 += estimate;
0283       //UNUSED:      tkDiff = estimate - m->estimate();
0284       tkNdof += hit->dimension();
0285     }
0286     if (hit->isValid() && (hit->geographicalId().det()) == DetId::Muon) {
0287       muChi2 += estimate;
0288       //UNUSED      staDiff = estimate - m->estimate();
0289       muNdof += hit->dimension();
0290     }
0291   }
0292 
0293   //For tkNdof < 6, should a large number or something else
0294   // be used instead of just tkChi2 directly?
0295   if (tkNdof > 5) {
0296     tkChi2 /= (tkNdof - 5.);
0297   }
0298 
0299   //For muNdof < 6, should a large number or something else
0300   // be used instead of just muChi2 directly?
0301   if (muNdof > 5) {
0302     muChi2 /= (muNdof - 5.);
0303   }
0304 
0305   return std::pair<double, double>(tkChi2, muChi2);
0306 }
0307 
0308 //
0309 // calculate the tail probability (-ln(P)) of a fit
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 //#include "FWCore/Framework/interface/MakerMacros.h"
0381 //DEFINE_FWK_MODULE(GlobalTrackQualityProducer);
[ 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!