CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_6/src/RecoMuon/MuonIdentification/src/DTTimingExtractor.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    MuonIdentification
00004 // Class:      DTTimingExtractor
00005 // 
00011 //
00012 // Original Author:  Traczyk Piotr
00013 //         Created:  Thu Oct 11 15:01:28 CEST 2007
00014 // $Id: DTTimingExtractor.cc,v 1.15 2011/02/24 15:41:53 farrell3 Exp $
00015 //
00016 //
00017 
00018 #include "RecoMuon/MuonIdentification/interface/DTTimingExtractor.h"
00019 
00020 
00021 // user include files
00022 #include "FWCore/Framework/interface/Frameworkfwd.h"
00023 #include "FWCore/Framework/interface/EDProducer.h"
00024 
00025 #include "FWCore/Framework/interface/Event.h"
00026 #include "FWCore/Framework/interface/MakerMacros.h"
00027 
00028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00029 
00030 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
00031 
00032 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00033 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00034 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00035 #include "Geometry/DTGeometry/interface/DTLayer.h"
00036 #include "Geometry/DTGeometry/interface/DTChamber.h"
00037 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
00038 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
00039 #include "DataFormats/DTRecHit/interface/DTRecSegment2D.h"
00040 
00041 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00042 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00043 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00044 
00045 #include "DataFormats/MuonReco/interface/Muon.h"
00046 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00047 #include "DataFormats/DTRecHit/interface/DTRecSegment2D.h"
00048 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00049 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
00050 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00051 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
00052 
00053 #include "DataFormats/TrackReco/interface/Track.h"
00054 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00055 #include "DataFormats/TrackReco/interface/TrackExtra.h"
00056 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
00057 
00058 
00059 // system include files
00060 #include <memory>
00061 #include <vector>
00062 #include <string>
00063 #include <iostream>
00064 
00065 namespace edm {
00066   class ParameterSet;
00067   class EventSetup;
00068   class InputTag;
00069 }
00070 
00071 class MuonServiceProxy;
00072 
00073 //
00074 // constructors and destructor
00075 //
00076 DTTimingExtractor::DTTimingExtractor(const edm::ParameterSet& iConfig)
00077   :
00078   DTSegmentTags_(iConfig.getParameter<edm::InputTag>("DTsegments")),
00079   theHitsMin_(iConfig.getParameter<int>("HitsMin")),
00080   thePruneCut_(iConfig.getParameter<double>("PruneCut")),
00081   theTimeOffset_(iConfig.getParameter<double>("DTTimeOffset")),
00082   theError_(iConfig.getParameter<double>("HitError")),
00083   useSegmentT0_(iConfig.getParameter<bool>("UseSegmentT0")),
00084   doWireCorr_(iConfig.getParameter<bool>("DoWireCorr")),
00085   dropTheta_(iConfig.getParameter<bool>("DropTheta")),
00086   requireBothProjections_(iConfig.getParameter<bool>("RequireBothProjections")),
00087   debug(iConfig.getParameter<bool>("debug"))
00088 {
00089   edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
00090   theService = new MuonServiceProxy(serviceParameters);
00091   
00092   edm::ParameterSet matchParameters = iConfig.getParameter<edm::ParameterSet>("MatchParameters");
00093 
00094   theMatcher = new MuonSegmentMatcher(matchParameters, theService);
00095 }
00096 
00097 
00098 DTTimingExtractor::~DTTimingExtractor()
00099 {
00100   if (theService) delete theService;
00101   if (theMatcher) delete theMatcher;
00102 }
00103 
00104 
00105 //
00106 // member functions
00107 //
00108 
00109 // ------------ method called to produce the data  ------------
00110 void
00111 DTTimingExtractor::fillTiming(TimeMeasurementSequence &tmSequence, reco::TrackRef muonTrack, const edm::Event& iEvent, const edm::EventSetup& iSetup)
00112 {
00113 
00114 //  using reco::TrackCollection;
00115 
00116   if (debug) 
00117     std::cout << " *** Muon Timimng Extractor ***" << std::endl;
00118 
00119   theService->update(iSetup);
00120 
00121   const GlobalTrackingGeometry *theTrackingGeometry = &*theService->trackingGeometry();
00122 
00123   // get the DT geometry
00124   edm::ESHandle<DTGeometry> theDTGeom;
00125   iSetup.get<MuonGeometryRecord>().get(theDTGeom);
00126   
00127   edm::ESHandle<Propagator> propagator;
00128   iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
00129   const Propagator *propag = propagator.product();
00130 
00131   double invbeta=0;
00132   double invbetaerr=0;
00133   double totalWeightInvbeta=0;
00134   double totalWeightVertex=0;
00135   std::vector<TimeMeasurement> tms;
00136 
00137   math::XYZPoint  pos=muonTrack->innerPosition();
00138   math::XYZVector mom=muonTrack->innerMomentum();
00139   GlobalPoint  posp(pos.x(), pos.y(), pos.z());
00140   GlobalVector momv(mom.x(), mom.y(), mom.z());
00141   FreeTrajectoryState muonFTS(posp, momv, (TrackCharge)muonTrack->charge(), theService->magneticField().product());
00142 
00143   // get the DT segments that were used to construct the muon
00144   std::vector<const DTRecSegment4D*> range = theMatcher->matchDT(*muonTrack,iEvent);
00145 
00146   // create a collection on TimeMeasurements for the track        
00147   for (std::vector<const DTRecSegment4D*>::iterator rechit = range.begin(); rechit!=range.end();++rechit) {
00148 
00149     // Create the ChamberId
00150     DetId id = (*rechit)->geographicalId();
00151     DTChamberId chamberId(id.rawId());
00152     int station = chamberId.station();
00153 
00154     // use only segments with both phi and theta projections present (optional)
00155     bool bothProjections = ( ((*rechit)->hasPhi()) && ((*rechit)->hasZed()) );
00156     
00157     if (requireBothProjections_ && !bothProjections) continue;
00158 
00159     // loop over (theta, phi) segments
00160     for (int phi=0; phi<2; phi++) {
00161 
00162       if (dropTheta_ && !phi) continue;
00163 
00164       const DTRecSegment2D* segm;
00165       if (phi) segm = dynamic_cast<const DTRecSegment2D*>((*rechit)->phiSegment()); 
00166         else segm = dynamic_cast<const DTRecSegment2D*>((*rechit)->zSegment());
00167 
00168       if(segm == 0) continue;
00169       if (!segm->specificRecHits().size()) continue;
00170 
00171       const GeomDet* geomDet = theTrackingGeometry->idToDet(segm->geographicalId());
00172       const std::vector<DTRecHit1D> hits1d = segm->specificRecHits();
00173 
00174       // store all the hits from the segment
00175       for (std::vector<DTRecHit1D>::const_iterator hiti=hits1d.begin(); hiti!=hits1d.end(); hiti++) {
00176 
00177         const GeomDet* dtcell = theTrackingGeometry->idToDet(hiti->geographicalId());
00178         TimeMeasurement thisHit;
00179 
00180         std::pair< TrajectoryStateOnSurface, double> tsos;
00181         tsos=propag->propagateWithPath(muonFTS,dtcell->surface());
00182 
00183         double dist;            
00184         double dist_straight = dtcell->toGlobal(hiti->localPosition()).mag(); 
00185         if (tsos.first.isValid()) { 
00186           dist = tsos.second+posp.mag(); 
00187 //        std::cout << "Propagate distance: " << dist << " ( innermost: " << posp.mag() << ")" << std::endl; 
00188         } else { 
00189           dist = dist_straight;
00190 //        std::cout << "Geom distance: " << dist << std::endl; 
00191         }
00192 
00193         thisHit.driftCell = hiti->geographicalId();
00194         if (hiti->lrSide()==DTEnums::Left) thisHit.isLeft=true; else thisHit.isLeft=false;
00195         thisHit.isPhi = phi;
00196         thisHit.posInLayer = geomDet->toLocal(dtcell->toGlobal(hiti->localPosition())).x();
00197         thisHit.distIP = dist;
00198         thisHit.station = station;
00199         if (useSegmentT0_ && segm->ist0Valid()) thisHit.timeCorr=segm->t0();
00200         else thisHit.timeCorr=0.;
00201         thisHit.timeCorr += theTimeOffset_;
00202           
00203         // signal propagation along the wire correction for unmached theta or phi segment hits
00204         if (doWireCorr_ && !bothProjections && tsos.first.isValid()) {
00205           const DTLayer* layer = theDTGeom->layer(hiti->wireId());
00206           float propgL = layer->toLocal( tsos.first.globalPosition() ).y();
00207           float wirePropCorr = propgL/24.4*0.00543; // signal propagation speed along the wire
00208           if (thisHit.isLeft) wirePropCorr=-wirePropCorr;
00209           thisHit.posInLayer += wirePropCorr;
00210           const DTSuperLayer *sl = layer->superLayer();
00211           float tofCorr = sl->surface().position().mag()-tsos.first.globalPosition().mag();
00212           tofCorr = (tofCorr/29.979)*0.00543;
00213           if (thisHit.isLeft) tofCorr=-tofCorr;
00214           thisHit.posInLayer += tofCorr;
00215         } else {
00216           // non straight-line path correction
00217           float slCorr = (dist_straight-dist)/29.979*0.00543;
00218           if (thisHit.isLeft) slCorr=-slCorr;
00219           thisHit.posInLayer += slCorr;
00220         }
00221 
00222         tms.push_back(thisHit);
00223       }
00224     } // phi = (0,1)            
00225   } // rechit
00226       
00227   bool modified = false;
00228   std::vector <double> dstnc, dsegm, dtraj, hitWeightVertex, hitWeightInvbeta, left;
00229     
00230   // Now loop over the measurements, calculate 1/beta and cut away outliers
00231   do {    
00232 
00233     modified = false;
00234     dstnc.clear();
00235     dsegm.clear();
00236     dtraj.clear();
00237     hitWeightVertex.clear();
00238     hitWeightInvbeta.clear();
00239     left.clear();
00240       
00241     std::vector <int> hit_idx;
00242     totalWeightInvbeta=0;
00243     totalWeightVertex=0;
00244       
00245     // Rebuild segments
00246     for (int sta=1;sta<5;sta++)
00247       for (int phi=0;phi<2;phi++) {
00248         std::vector <TimeMeasurement> seg;
00249         std::vector <int> seg_idx;
00250         int tmpos=0;
00251         for (std::vector<TimeMeasurement>::iterator tm=tms.begin(); tm!=tms.end(); ++tm) {
00252           if ((tm->station==sta) && (tm->isPhi==phi)) {
00253             seg.push_back(*tm);
00254             seg_idx.push_back(tmpos);
00255           }
00256           tmpos++;  
00257         }
00258 
00259         unsigned int segsize = seg.size();
00260         if (segsize<theHitsMin_) continue;
00261 
00262         double a=0, b=0;
00263         std::vector <double> hitxl,hitxr,hityl,hityr;
00264 
00265         for (std::vector<TimeMeasurement>::iterator tm=seg.begin(); tm!=seg.end(); ++tm) {
00266  
00267           DetId id = tm->driftCell;
00268           const GeomDet* dtcell = theTrackingGeometry->idToDet(id);
00269           DTChamberId chamberId(id.rawId());
00270           const GeomDet* dtcham = theTrackingGeometry->idToDet(chamberId);
00271 
00272           double celly=dtcham->toLocal(dtcell->position()).z();
00273             
00274           if (tm->isLeft) {
00275             hitxl.push_back(celly);
00276             hityl.push_back(tm->posInLayer);
00277           } else {
00278             hitxr.push_back(celly);
00279             hityr.push_back(tm->posInLayer);
00280           }    
00281         }
00282 
00283         if (!fitT0(a,b,hitxl,hityl,hitxr,hityr)) {
00284           if (debug)
00285             std::cout << "     t0 = zero, Left hits: " << hitxl.size() << " Right hits: " << hitxr.size() << std::endl;
00286           continue;
00287         }
00288           
00289         // a segment must have at least one left and one right hit
00290         if ((!hitxl.size()) || (!hityl.size())) continue;
00291 
00292         int segidx=0;
00293         for (std::vector<TimeMeasurement>::const_iterator tm=seg.begin(); tm!=seg.end(); ++tm) {
00294 
00295           DetId id = tm->driftCell;
00296           const GeomDet* dtcell = theTrackingGeometry->idToDet(id);
00297           DTChamberId chamberId(id.rawId());
00298           const GeomDet* dtcham = theTrackingGeometry->idToDet(chamberId);
00299 
00300           double layerZ  = dtcham->toLocal(dtcell->position()).z();
00301           double segmLocalPos = b+layerZ*a;
00302           double hitLocalPos = tm->posInLayer;
00303           int hitSide = -tm->isLeft*2+1;
00304           double t0_segm = (-(hitSide*segmLocalPos)+(hitSide*hitLocalPos))/0.00543+tm->timeCorr;
00305             
00306           dstnc.push_back(tm->distIP);
00307           dsegm.push_back(t0_segm);
00308           left.push_back(hitSide);
00309           hitWeightInvbeta.push_back(((double)seg.size()-2.)*tm->distIP*tm->distIP/((double)seg.size()*30.*30.*theError_*theError_));
00310           hitWeightVertex.push_back(((double)seg.size()-2.)/((double)seg.size()*theError_*theError_));
00311           hit_idx.push_back(seg_idx.at(segidx));
00312           segidx++;
00313           totalWeightInvbeta+=((double)seg.size()-2.)*tm->distIP*tm->distIP/((double)seg.size()*30.*30.*theError_*theError_);
00314           totalWeightVertex+=((double)seg.size()-2.)/((double)seg.size()*theError_*theError_);
00315         }
00316       }
00317 
00318     if (totalWeightInvbeta==0) break;        
00319 
00320     // calculate the value and error of 1/beta from the complete set of 1D hits
00321     if (debug)
00322       std::cout << " Points for global fit: " << dstnc.size() << std::endl;
00323 
00324     // inverse beta - weighted average of the contributions from individual hits
00325     invbeta=0;
00326     for (unsigned int i=0;i<dstnc.size();i++) 
00327       invbeta+=(1.+dsegm.at(i)/dstnc.at(i)*30.)*hitWeightInvbeta.at(i)/totalWeightInvbeta;
00328 
00329     double chimax=0.;
00330     std::vector<TimeMeasurement>::iterator tmmax;
00331     
00332     // the dispersion of inverse beta
00333     double diff;
00334     for (unsigned int i=0;i<dstnc.size();i++) {
00335       diff=(1.+dsegm.at(i)/dstnc.at(i)*30.)-invbeta;
00336       diff=diff*diff*hitWeightInvbeta.at(i);
00337       invbetaerr+=diff;
00338       if (diff>chimax) { 
00339         tmmax=tms.begin()+hit_idx.at(i);
00340         chimax=diff;
00341       }
00342     }
00343     
00344     invbetaerr=sqrt(invbetaerr/totalWeightInvbeta); 
00345  
00346     // cut away the outliers
00347     if (chimax>thePruneCut_) {
00348       tms.erase(tmmax);
00349       modified=true;
00350     }    
00351 
00352     if (debug)
00353       std::cout << " Measured 1/beta: " << invbeta << " +/- " << invbetaerr << std::endl;
00354 
00355   } while (modified);
00356 
00357   for (unsigned int i=0;i<dstnc.size();i++) {
00358     tmSequence.dstnc.push_back(dstnc.at(i));
00359     tmSequence.local_t0.push_back(dsegm.at(i));
00360     tmSequence.weightInvbeta.push_back(hitWeightInvbeta.at(i));
00361     tmSequence.weightVertex.push_back(hitWeightVertex.at(i));
00362   }
00363 
00364   tmSequence.totalWeightInvbeta=totalWeightInvbeta;
00365   tmSequence.totalWeightVertex=totalWeightVertex;
00366 
00367 }
00368 
00369 double
00370 DTTimingExtractor::fitT0(double &a, double &b, std::vector<double> xl, std::vector<double> yl, std::vector<double> xr, std::vector<double> yr ) {
00371 
00372   double sx=0,sy=0,sxy=0,sxx=0,ssx=0,ssy=0,s=0,ss=0;
00373 
00374   for (unsigned int i=0; i<xl.size(); i++) {
00375     sx+=xl[i];
00376     sy+=yl[i];
00377     sxy+=xl[i]*yl[i];
00378     sxx+=xl[i]*xl[i];
00379     s++;
00380     ssx+=xl[i];
00381     ssy+=yl[i];
00382     ss++;
00383   } 
00384 
00385   for (unsigned int i=0; i<xr.size(); i++) {
00386     sx+=xr[i];
00387     sy+=yr[i];
00388     sxy+=xr[i]*yr[i];
00389     sxx+=xr[i]*xr[i];
00390     s++;
00391     ssx-=xr[i];
00392     ssy-=yr[i];
00393     ss--;
00394   } 
00395 
00396   double delta = ss*ss*sxx+s*sx*sx+s*ssx*ssx-s*s*sxx-2*ss*sx*ssx;
00397   
00398   double t0_corr=0.;
00399 
00400   if (delta) {
00401     a=(ssy*s*ssx+sxy*ss*ss+sy*sx*s-sy*ss*ssx-ssy*sx*ss-sxy*s*s)/delta;
00402     b=(ssx*sy*ssx+sxx*ssy*ss+sx*sxy*s-sxx*sy*s-ssx*sxy*ss-sx*ssy*ssx)/delta;
00403     t0_corr=(ssx*s*sxy+sxx*ss*sy+sx*sx*ssy-sxx*s*ssy-sx*ss*sxy-ssx*sx*sy)/delta;
00404   }
00405 
00406   // convert drift distance to time
00407   t0_corr/=-0.00543;
00408 
00409   return t0_corr;
00410 }
00411 
00412 
00413 //define this as a plug-in
00414 //DEFINE_FWK_MODULE(DTTimingExtractor);