CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_1/src/RecoMuon/GlobalTrackingTools/plugins/GlobalTrackQualityProducer.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    GlobalTrackingTools
00004 // Class:      GlobalTrackQualityProducer
00005 //
00006 //
00007 // Original Author:  Adam Everett
00008 // $Id: GlobalTrackQualityProducer.cc,v 1.8 2012/12/06 14:45:38 eulisse Exp $
00009 //
00010 //
00011 
00012 // system include files
00013 #include <memory>
00014 
00015 // user include files
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   // service parameters
00033   edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
00034   theService = new MuonServiceProxy(serviceParameters);     
00035   
00036   // TrackRefitter parameters
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   // Take the GLB muon container(s)
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   // reserve some space
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; //normalized inside to /sum(muHits.dimension)
00096       relative_tracker_chi2 = chi.first; //normalized inside to /sum(tkHits.dimension)
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     // Fill the STA-TK match information
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; // a better solution would be to use float above .. m/be not
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   for(int i = 0; i < valuesTkRelChi2.size(); i++) {
00153     LogTrace(theCategory)<<"value " << valuesTkRelChi2[i] ;
00154   }
00155   */
00156 
00157   // create and fill value maps
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   // put value map into event
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     //not used    double estimate = 0.0;
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     //if ( !ok ) continue;
00198     
00199     const TrajectoryStateOnSurface& tsos = (*m).predictedState();
00200 
00201 
00202     if ( tsos.isValid() && rhit->isValid() && rhit->hit()->isValid()
00203          && !edm::isNotFinite(rhit->localPositionError().xx()) //this is paranoia induced by reported case
00204          && !edm::isNotFinite(rhit->localPositionError().xy()) //it's better to track down the origin of bad numbers
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       //LogDebug(theCategory)<<"hitPos " << hitPos;
00221       double error = hitErr.phierr(hitPos);  // error squared
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     //LogTrace(theCategory) << "estimate " << estimate << " TM.est " << m->estimate();
00264     //UNUSED:    double tkDiff = 0.0;
00265     //UNUSED:    double staDiff = 0.0;
00266     if ( hit->isValid() &&  (hit->geographicalId().det()) == DetId::Tracker ) {
00267       tkChi2 += estimate;
00268       //UNUSED:      tkDiff = estimate - m->estimate();
00269       tkNdof += hit->dimension();
00270     }
00271     if ( hit->isValid() &&  (hit->geographicalId().det()) == DetId::Muon ) {
00272       muChi2 += estimate;
00273       //UNUSED      staDiff = estimate - m->estimate();
00274       muNdof += hit->dimension();
00275     }
00276   }
00277   
00278   if (tkNdof < 6 ) tkChi2 = tkChi2; // or should I set it to  a large number ?
00279   else tkChi2 /= (tkNdof-5.);
00280 
00281   if (muNdof < 6 ) muChi2 = muChi2; // or should I set it to  a large number ?
00282   else muChi2 /= (muNdof-5.);
00283 
00284   return std::pair<double,double>(tkChi2,muChi2);
00285        
00286 }
00287 
00288 //
00289 // calculate the tail probability (-ln(P)) of a fit
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 //#include "FWCore/Framework/interface/MakerMacros.h"
00303 //DEFINE_FWK_MODULE(GlobalTrackQualityProducer);