Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <memory>
00014
00015
00016 #include "FWCore/Framework/interface/Frameworkfwd.h"
00017 #include "FWCore/Framework/interface/EDProducer.h"
00018 #include "FWCore/Utilities/interface/isFinite.h"
00019
00020 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
00021 #include "RecoMuon/GlobalTrackingTools/plugins/GlobalTrackQualityProducer.h"
00022
00023 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00024 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00025 #include "DataFormats/MuonReco/interface/MuonQuality.h"
00026 #include "DataFormats/MuonReco/interface/MuonTrackLinks.h"
00027 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00028
00029 GlobalTrackQualityProducer::GlobalTrackQualityProducer(const edm::ParameterSet& iConfig):
00030 inputCollection_(iConfig.getParameter<edm::InputTag>("InputCollection")),inputLinksCollection_(iConfig.getParameter<edm::InputTag>("InputLinksCollection")),theService(0),theGlbRefitter(0),theGlbMatcher(0)
00031 {
00032
00033 edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
00034 theService = new MuonServiceProxy(serviceParameters);
00035
00036
00037 edm::ParameterSet refitterParameters = iConfig.getParameter<edm::ParameterSet>("RefitterParameters");
00038 theGlbRefitter = new GlobalMuonRefitter(refitterParameters, theService);
00039
00040 edm::ParameterSet trackMatcherPSet = iConfig.getParameter<edm::ParameterSet>("GlobalMuonTrackMatcher");
00041 theGlbMatcher = new GlobalMuonTrackMatcher(trackMatcherPSet,theService);
00042
00043 double maxChi2 = iConfig.getParameter<double>("MaxChi2");
00044 double nSigma = iConfig.getParameter<double>("nSigma");
00045 theEstimator = new Chi2MeasurementEstimator(maxChi2,nSigma);
00046
00047 produces<edm::ValueMap<reco::MuonQuality> >();
00048 }
00049
00050 GlobalTrackQualityProducer::~GlobalTrackQualityProducer() {
00051 if (theService) delete theService;
00052 if (theGlbRefitter) delete theGlbRefitter;
00053 if (theGlbMatcher) delete theGlbMatcher;
00054 if (theEstimator) delete theEstimator;
00055 }
00056
00057 void
00058 GlobalTrackQualityProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00059 {
00060 const std::string theCategory = "Muon|RecoMuon|GlobalTrackQualityProducer";
00061
00062 theService->update(iSetup);
00063
00064 theGlbRefitter->setEvent(iEvent);
00065
00066 theGlbRefitter->setServices(theService->eventSetup());
00067
00068
00069 edm::Handle<reco::TrackCollection> glbMuons;
00070 iEvent.getByLabel(inputCollection_,glbMuons);
00071
00072 edm::Handle<reco::MuonTrackLinksCollection> linkCollectionHandle;
00073 iEvent.getByLabel(inputLinksCollection_, linkCollectionHandle);
00074
00075
00076 std::vector<reco::MuonQuality> valuesQual;
00077 valuesQual.reserve(glbMuons->size());
00078
00079 int trackIndex = 0;
00080 for (reco::TrackCollection::const_iterator track = glbMuons->begin(); track!=glbMuons->end(); ++track , ++trackIndex) {
00081 reco::TrackRef glbRef(glbMuons,trackIndex);
00082 reco::TrackRef staTrack = reco::TrackRef();
00083
00084 std::vector<Trajectory> refitted=theGlbRefitter->refit(*track,1);
00085
00086 LogTrace(theCategory)<<"GLBQual N refitted " << refitted.size();
00087
00088 std::pair<double,double> thisKink;
00089 double relative_muon_chi2 = 0.0;
00090 double relative_tracker_chi2 = 0.0;
00091 double glbTrackProbability = 0.0;
00092 if(refitted.size()>0) {
00093 thisKink = kink(refitted.front()) ;
00094 std::pair<double,double> chi = newChi2(refitted.front());
00095 relative_muon_chi2 = chi.second;
00096 relative_tracker_chi2 = chi.first;
00097 glbTrackProbability = trackProbability(refitted.front());
00098 }
00099
00100 LogTrace(theCategory)<<"GLBQual: Kink " << thisKink.first << " " << thisKink.second;
00101 LogTrace(theCategory)<<"GLBQual: Rel Chi2 " << relative_tracker_chi2 << " " << relative_muon_chi2;
00102 LogTrace(theCategory)<<"GLBQual: trackProbability " << glbTrackProbability;
00103
00104
00105 float chi2, d, dist, Rpos;
00106 chi2 = d = dist = Rpos = -1.0;
00107 bool passTight = false;
00108 typedef MuonTrajectoryBuilder::TrackCand TrackCand;
00109 if ( linkCollectionHandle.isValid() ) {
00110 for ( reco::MuonTrackLinksCollection::const_iterator links = linkCollectionHandle->begin();
00111 links != linkCollectionHandle->end(); ++links )
00112 {
00113 if ( links->trackerTrack().isNull() ||
00114 links->standAloneTrack().isNull() ||
00115 links->globalTrack().isNull() )
00116 {
00117 edm::LogWarning(theCategory) << "Global muon links to constituent tracks are invalid. There should be no such object. Muon is skipped.";
00118 continue;
00119 }
00120 if (links->globalTrack() == glbRef) {
00121 staTrack = !links->standAloneTrack().isNull() ? links->standAloneTrack() : reco::TrackRef();
00122 TrackCand staCand = TrackCand((Trajectory*)(0),links->standAloneTrack());
00123 TrackCand tkCand = TrackCand((Trajectory*)(0),links->trackerTrack());
00124 chi2 = theGlbMatcher->match(staCand,tkCand,0,0);
00125 d = theGlbMatcher->match(staCand,tkCand,1,0);
00126 Rpos = theGlbMatcher->match(staCand,tkCand,2,0);
00127 dist = theGlbMatcher->match(staCand,tkCand,3,0);
00128 passTight = theGlbMatcher->matchTight(staCand,tkCand);
00129 }
00130 }
00131 }
00132
00133 if(!staTrack.isNull()) LogTrace(theCategory)<<"GLBQual: Used UpdatedAtVtx : " << (iEvent.getProvenance(staTrack.id()).productInstanceName() == std::string("UpdatedAtVtx"));
00134
00135 float maxFloat01 = std::numeric_limits<float>::max()*0.1;
00136 reco::MuonQuality muQual;
00137 if(!staTrack.isNull()) muQual.updatedSta = iEvent.getProvenance(staTrack.id()).productInstanceName() == std::string("UpdatedAtVtx");
00138 muQual.trkKink = thisKink.first > maxFloat01 ? maxFloat01 : thisKink.first;
00139 muQual.glbKink = thisKink.second > maxFloat01 ? maxFloat01 : thisKink.second;
00140 muQual.trkRelChi2 = relative_tracker_chi2 > maxFloat01 ? maxFloat01 : relative_tracker_chi2;
00141 muQual.staRelChi2 = relative_muon_chi2 > maxFloat01 ? maxFloat01 : relative_muon_chi2;
00142 muQual.tightMatch = passTight;
00143 muQual.chi2LocalPosition = dist;
00144 muQual.chi2LocalMomentum = chi2;
00145 muQual.localDistance = d;
00146 muQual.globalDeltaEtaPhi = Rpos;
00147 muQual.glbTrackProbability = glbTrackProbability;
00148 valuesQual.push_back(muQual);
00149 }
00150
00151
00152
00153
00154
00155
00156
00157
00158 std::auto_ptr<edm::ValueMap<reco::MuonQuality> > outQual(new edm::ValueMap<reco::MuonQuality>());
00159 edm::ValueMap<reco::MuonQuality>::Filler fillerQual(*outQual);
00160 fillerQual.insert(glbMuons, valuesQual.begin(), valuesQual.end());
00161 fillerQual.fill();
00162
00163
00164 iEvent.put(outQual);
00165 }
00166
00167 std::pair<double,double> GlobalTrackQualityProducer::kink(Trajectory& muon) const {
00168
00169 const std::string theCategory = "Muon|RecoMuon|GlobalTrackQualityProducer";
00170
00171 using namespace std;
00172 using namespace edm;
00173 using namespace reco;
00174
00175 double result = 0.0;
00176 double resultGlb = 0.0;
00177
00178
00179 typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
00180 typedef ConstRecHitPointer RecHit;
00181 typedef std::vector<TrajectoryMeasurement>::const_iterator TMI;
00182
00183
00184 vector<TrajectoryMeasurement> meas = muon.measurements();
00185
00186 for ( TMI m = meas.begin(); m != meas.end(); m++ ) {
00187 TransientTrackingRecHit::ConstRecHitPointer hit = m->recHit();
00188
00189
00190
00191 RecHit rhit = (*m).recHit();
00192 bool ok = false;
00193 if ( rhit->isValid() ) {
00194 if(DetId::Tracker == rhit->geographicalId().det()) ok = true;
00195 }
00196
00197
00198
00199 const TrajectoryStateOnSurface& tsos = (*m).predictedState();
00200
00201
00202 if ( tsos.isValid() && rhit->isValid() && rhit->hit()->isValid()
00203 && !edm::isNotFinite(rhit->localPositionError().xx())
00204 && !edm::isNotFinite(rhit->localPositionError().xy())
00205 && !edm::isNotFinite(rhit->localPositionError().yy())
00206 ) {
00207
00208 double phi1 = tsos.globalPosition().phi();
00209 if ( phi1 < 0 ) phi1 = 2*M_PI + phi1;
00210
00211 double phi2 = rhit->globalPosition().phi();
00212 if ( phi2 < 0 ) phi2 = 2*M_PI + phi2;
00213
00214 double diff = fabs(phi1 - phi2);
00215 if ( diff > M_PI ) diff = 2*M_PI - diff;
00216
00217 GlobalPoint hitPos = rhit->globalPosition();
00218
00219 GlobalError hitErr = rhit->globalPositionError();
00220
00221 double error = hitErr.phierr(hitPos);
00222
00223 double s = ( error > 0.0 ) ? (diff*diff)/error : (diff*diff);
00224
00225 if(ok) result += s;
00226 resultGlb += s;
00227 }
00228
00229 }
00230
00231
00232 return std::pair<double,double>(result,resultGlb);
00233
00234 }
00235
00236 std::pair<double,double> GlobalTrackQualityProducer::newChi2(Trajectory& muon) const {
00237 const std::string theCategory = "Muon|RecoMuon|GlobalTrackQualityProducer";
00238
00239 using namespace std;
00240 using namespace edm;
00241 using namespace reco;
00242
00243 double muChi2 = 0.0;
00244 double tkChi2 = 0.0;
00245 unsigned int muNdof = 0;
00246 unsigned int tkNdof = 0;
00247
00248 typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
00249 typedef ConstRecHitPointer RecHit;
00250 typedef vector<TrajectoryMeasurement>::const_iterator TMI;
00251
00252 vector<TrajectoryMeasurement> meas = muon.measurements();
00253
00254 for ( TMI m = meas.begin(); m != meas.end(); m++ ) {
00255 TransientTrackingRecHit::ConstRecHitPointer hit = m->recHit();
00256 const TrajectoryStateOnSurface& uptsos = (*m).updatedState();
00257 TransientTrackingRecHit::RecHitPointer preciseHit = hit->clone(uptsos);
00258 double estimate = 0.0;
00259 if (preciseHit->isValid() && uptsos.isValid()) {
00260 estimate = theEstimator->estimate(uptsos, *preciseHit ).second;
00261 }
00262
00263
00264
00265
00266 if ( hit->isValid() && (hit->geographicalId().det()) == DetId::Tracker ) {
00267 tkChi2 += estimate;
00268
00269 tkNdof += hit->dimension();
00270 }
00271 if ( hit->isValid() && (hit->geographicalId().det()) == DetId::Muon ) {
00272 muChi2 += estimate;
00273
00274 muNdof += hit->dimension();
00275 }
00276 }
00277
00278 if (tkNdof < 6 ) tkChi2 = tkChi2;
00279 else tkChi2 /= (tkNdof-5.);
00280
00281 if (muNdof < 6 ) muChi2 = muChi2;
00282 else muChi2 /= (muNdof-5.);
00283
00284 return std::pair<double,double>(tkChi2,muChi2);
00285
00286 }
00287
00288
00289
00290
00291 double
00292 GlobalTrackQualityProducer::trackProbability(Trajectory& track) const {
00293
00294 if ( track.ndof() > 0 && track.chiSquared() > 0 ) {
00295 return -LnChiSquaredProbability(track.chiSquared(), track.ndof());
00296 } else {
00297 return 0.0;
00298 }
00299
00300 }
00301
00302
00303