CMSSW Reference Manual

MuonTimingExtractor::MuonTimingExtractor ( const edm::ParameterSet & iConfig )

Constructor.

Definition at line 76 of file MuonTimingExtractor.cc.

References edm::ParameterSet::getParameter(), MuonSegmentMatcher_cff::MuonSegmentMatcher, MuonServiceProxy_cff::MuonServiceProxy, theMatcher, and theService.

00077   :
00078   DTSegmentTags_(iConfig.getUntrackedParameter<edm::InputTag>("DTsegments")),
00079   theHitsMin(iConfig.getParameter<int>("HitsMin")),
00080   thePruneCut(iConfig.getParameter<double>("PruneCut")),
00081   useSegmentT0(iConfig.getParameter<bool>("UseSegmentT0")),
00082   debug(iConfig.getParameter<bool>("debug"))
00083 {
00084   // service parameters
00085   ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
00086   // the services
00087   theService = new MuonServiceProxy(serviceParameters);
00088   
00089   ParameterSet matchParameters = iConfig.getParameter<edm::ParameterSet>("MatchParameters");
00090 
00091   theMatcher = new MuonSegmentMatcher(matchParameters, theService);
00092 
00093 }

MuonTimingExtractor::~MuonTimingExtractor ( )

Destructor.

Definition at line 96 of file MuonTimingExtractor.cc.

References theService.

00097 {
00098  
00099   if (theService) delete theService;
00100 
00101 }

reco::MuonTime MuonTimingExtractor::fillTiming	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup,
		reco::TrackRef	muonTrack
	)

Definition at line 110 of file MuonTimingExtractor.cc.

References a, b, GenMuonPlsPt100GeV_cfg::cout, debug, diff, dist(), MuonTimingExtractor::TimeMeasurement::distIP, MuonTimingExtractor::TimeMeasurement::driftCell, lat::endl(), fitT0(), reco::MuonTime::freeInverseBeta, reco::MuonTime::freeInverseBetaErr, TrackingRecHit::geographicalId(), edm::EventSetup::get(), i, GlobalTrackingGeometry::idToDet(), reco::MuonTime::inverseBeta, reco::MuonTime::inverseBetaErr, MuonTimingExtractor::TimeMeasurement::isLeft, MuonTimingExtractor::TimeMeasurement::isPhi, DTEnums::Left, muonGeometry::mag(), MuonSegmentMatcher::matchDT(), reco::MuonTime::nStations, phi, MuonTimingExtractor::TimeMeasurement::posInLayer, GeomDet::position(), edm::ESHandle< T >::product(), Propagator::propagateWithPath(), MuonServiceProxy::propagator(), range, rawFit(), seg, DTRecSegment2D::specificRecHits(), funct::sqrt(), DTChamberId::station(), MuonTimingExtractor::TimeMeasurement::station, GeomDet::surface(), DTRecSegment2D::t0(), theHitsMin, theMatcher, thePruneCut, theService, MuonServiceProxy::theTrackingGeometry, reco::MuonTime::timeAtIpInOut, reco::MuonTime::timeAtIpInOutErr, reco::MuonTime::timeAtIpOutIn, reco::MuonTime::timeAtIpOutInErr, MuonTimingExtractor::TimeMeasurement::timeCorr, GeomDet::toGlobal(), GeomDet::toLocal(), useSegmentT0, PV3DBase< T, PVType, FrameType >::x(), x, y, and z.

Referenced by MuonIdProducer::produce().

00111 {
00112 
00113   //using namespace edm;
00114   using reco::TrackCollection;
00115 
00116   if (debug) 
00117     cout << " *** Muon Timimng Extractor ***" << endl;
00118 
00119   theService->update(iSetup);
00120 
00121   const GlobalTrackingGeometry *theTrackingGeometry = &*theService->trackingGeometry();
00122   
00123   edm::ESHandle<Propagator> propagator;
00124   iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
00125   const Propagator *propag = propagator.product();
00126 
00127   double invbeta=0;
00128   double invbetaerr=0;
00129   int totalWeight=0;
00130   int nStations=0;
00131   vector<TimeMeasurement> tms;
00132     
00133   math::XYZPoint  pos=muonTrack->innerPosition();
00134   math::XYZVector mom=muonTrack->innerMomentum();
00135   GlobalPoint  posp(pos.x(), pos.y(), pos.z());
00136   GlobalVector momv(mom.x(), mom.y(), mom.z());
00137   FreeTrajectoryState muonFTS(posp, momv, (TrackCharge)muonTrack->charge(), theService->magneticField().product());
00138 
00139   // get the DT segments that were used to construct the muon
00140   vector<const DTRecSegment4D*> range = theMatcher->matchDT(*muonTrack,iEvent);
00141 
00142   // create a collection on TimeMeasurements for the track        
00143   for (vector<const DTRecSegment4D*>::iterator rechit = range.begin(); rechit!=range.end();++rechit) {
00144 
00145     // Create the ChamberId
00146     DetId id = (*rechit)->geographicalId();
00147     DTChamberId chamberId(id.rawId());
00148     int station = chamberId.station();
00149 
00150     // use only segments with both phi and theta projections present
00151     if ( (!(*rechit)->hasPhi()) || (!(*rechit)->hasZed()) ) continue;
00152 
00153     // loop over (theta, phi) segments
00154     for (int phi=0; phi<2; phi++) {
00155 
00156       const DTRecSegment2D* segm;
00157       if (phi) segm = dynamic_cast<const DTRecSegment2D*>((*rechit)->phiSegment()); 
00158       else segm = dynamic_cast<const DTRecSegment2D*>((*rechit)->zSegment());
00159       if(segm == 0){cout << "skip this" << endl;  continue; }   
00160       if (!segm->specificRecHits().size()) continue;
00161 
00162       const GeomDet* geomDet = theTrackingGeometry->idToDet(segm->geographicalId());
00163       const vector<DTRecHit1D> hits1d = segm->specificRecHits();
00164 
00165       // store all the hits from the segment
00166       for (vector<DTRecHit1D>::const_iterator hiti=hits1d.begin(); hiti!=hits1d.end(); hiti++) {
00167 
00168         const GeomDet* dtcell = theTrackingGeometry->idToDet(hiti->geographicalId());
00169         TimeMeasurement thisHit;
00170 
00171         std::pair< TrajectoryStateOnSurface, double> tsos;
00172         tsos=propag->propagateWithPath(muonFTS,dtcell->surface());
00173 
00174         double dist;            
00175         if (tsos.first.isValid()) dist = tsos.second+posp.mag(); 
00176           else dist = dtcell->toGlobal(hiti->localPosition()).mag();
00177 
00178         thisHit.driftCell = hiti->geographicalId();
00179         if (hiti->lrSide()==DTEnums::Left) thisHit.isLeft=true; else thisHit.isLeft=false;
00180         thisHit.isPhi = phi;
00181         thisHit.posInLayer = geomDet->toLocal(dtcell->toGlobal(hiti->localPosition())).x();
00182         thisHit.distIP = dist;
00183         thisHit.station = station;
00184         if (useSegmentT0) thisHit.timeCorr=segm->t0();
00185           else thisHit.timeCorr=0.;
00186         tms.push_back(thisHit);
00187       }
00188     } // phi = (0,1)            
00189   } // rechit
00190       
00191   bool modified = false;
00192   vector <double> dstnc, dsegm, dtraj, hitWeight, left;
00193     
00194   // Now loop over the measurements, calculate 1/beta and cut away outliers
00195   do {    
00196 
00197     modified = false;
00198     dstnc.clear();
00199     dsegm.clear();
00200     dtraj.clear();
00201     hitWeight.clear();
00202     left.clear();
00203       
00204     vector <int> hit_idx;
00205     totalWeight=0;
00206     nStations=0;
00207       
00208     // Rebuild segments
00209     for (int sta=1;sta<5;sta++)
00210       for (int phi=0;phi<2;phi++) {
00211         vector <TimeMeasurement> seg;
00212         vector <int> seg_idx;
00213         int tmpos=0;
00214         for (vector<TimeMeasurement>::iterator tm=tms.begin(); tm!=tms.end(); ++tm) {
00215           if ((tm->station==sta) && (tm->isPhi==phi)) {
00216             seg.push_back(*tm);
00217             seg_idx.push_back(tmpos);
00218           }
00219           tmpos++;  
00220         }
00221           
00222         unsigned int segsize = seg.size();
00223           
00224         if (segsize<theHitsMin) continue;
00225 
00226         double a=0, b=0;
00227         vector <double> hitxl,hitxr,hityl,hityr;
00228 
00229         for (vector<TimeMeasurement>::iterator tm=seg.begin(); tm!=seg.end(); ++tm) {
00230  
00231           DetId id = tm->driftCell;
00232           const GeomDet* dtcell = theTrackingGeometry->idToDet(id);
00233           DTChamberId chamberId(id.rawId());
00234           const GeomDet* dtcham = theTrackingGeometry->idToDet(chamberId);
00235 
00236           double celly=dtcham->toLocal(dtcell->position()).z();
00237             
00238           if (tm->isLeft) {
00239             hitxl.push_back(celly);
00240             hityl.push_back(tm->posInLayer);
00241           } else {
00242             hitxr.push_back(celly);
00243             hityr.push_back(tm->posInLayer);
00244           }    
00245         }
00246 
00247         if (!fitT0(a,b,hitxl,hityl,hitxr,hityr)) {
00248           if (debug)
00249             cout << "     t0 = zero, Left hits: " << hitxl.size() << " Right hits: " << hitxr.size() << endl;
00250           continue;
00251         }
00252           
00253         // a segment must have at least one left and one right hit
00254         if ((!hitxl.size()) || (!hityl.size())) continue;
00255 
00256         nStations++;
00257 
00258         int segidx=0;
00259         for (vector<TimeMeasurement>::const_iterator tm=seg.begin(); tm!=seg.end(); ++tm) {
00260 
00261           DetId id = tm->driftCell;
00262           const GeomDet* dtcell = theTrackingGeometry->idToDet(id);
00263           DTChamberId chamberId(id.rawId());
00264           const GeomDet* dtcham = theTrackingGeometry->idToDet(chamberId);
00265 
00266           double layerZ  = dtcham->toLocal(dtcell->position()).z();
00267           double segmLocalPos = b+layerZ*a;
00268           double hitLocalPos = tm->posInLayer;
00269           int hitSide = -tm->isLeft*2+1;
00270           double t0_segm = (-(hitSide*segmLocalPos)+(hitSide*hitLocalPos))/0.00543+tm->timeCorr;
00271             
00272           dstnc.push_back(tm->distIP);
00273           dsegm.push_back(t0_segm);
00274           left.push_back(hitSide);
00275           hitWeight.push_back(((double)seg.size()-2.)/(double)seg.size());
00276           hit_idx.push_back(seg_idx.at(segidx));
00277           segidx++;
00278         }
00279           
00280         totalWeight+=seg.size()-2;
00281 
00282       }
00283 
00284     nStations/=2;
00285     invbeta=0;
00286 
00287     // calculate the value and error of 1/beta from the complete set of 1D hits
00288     if (debug)
00289       cout << " Points for global fit: " << dstnc.size() << endl;
00290 
00291     // inverse beta - weighted average of the contributions from individual hits
00292     for (unsigned int i=0;i<dstnc.size();i++) 
00293       invbeta+=(1.+dsegm.at(i)/dstnc.at(i)*30.)*hitWeight.at(i)/totalWeight;
00294 
00295     double chimax=0.;
00296     vector<TimeMeasurement>::iterator tmmax;
00297     
00298     // the dispersion of inverse beta
00299     double diff;
00300     for (unsigned int i=0;i<dstnc.size();i++) {
00301       diff=(1.+dsegm.at(i)/dstnc.at(i)*30.)-invbeta;
00302       diff=diff*diff*hitWeight.at(i);
00303       invbetaerr+=diff;
00304       if (diff/totalWeight/totalWeight>chimax/100.) { 
00305         tmmax=tms.begin()+hit_idx.at(i);
00306         chimax=diff;
00307       }
00308     }
00309     
00310     invbetaerr=sqrt(invbetaerr/totalWeight);
00311  
00312     // cut away the outliers
00313     if (chimax>thePruneCut) {
00314       tms.erase(tmmax);
00315       modified=true;
00316     }    
00317 
00318     if (debug)
00319       cout << " Measured 1/beta: " << invbeta << " +/- " << invbetaerr << endl;
00320 
00321   } while (modified);
00322 
00323   vector <double> x,y;
00324   double freeBeta, freeBetaErr, freeTime, freeTimeErr;
00325   double vertexTime=0, vertexTimeErr=0, vertexTimeR=0, vertexTimeRErr=0;    
00326 
00327   for (unsigned int i=0;i<dstnc.size();i++) {
00328     x.push_back(dstnc.at(i)/30.);
00329     y.push_back(dsegm.at(i)+dstnc.at(i)/30.);
00330     vertexTime+=dsegm.at(i)*hitWeight.at(i)/totalWeight;
00331     vertexTimeR+=(dsegm.at(i)+2*dstnc.at(i)/30.)*hitWeight.at(i)/totalWeight;
00332   }
00333 
00334   double diff;
00335   for (unsigned int i=0;i<dstnc.size();i++) {
00336     diff=dsegm.at(i)-vertexTime;
00337     vertexTimeErr+=diff*diff*hitWeight.at(i);
00338     diff=dsegm.at(i)+2*dstnc.at(i)/30.-vertexTimeR;
00339     vertexTimeRErr+=diff*diff*hitWeight.at(i);
00340   }
00341   vertexTimeErr=sqrt(vertexTimeErr/totalWeight);
00342   vertexTimeRErr=sqrt(vertexTimeRErr/totalWeight);
00343 
00344   reco::MuonTime muonTime;
00345 
00346   muonTime.inverseBeta=invbeta;                        
00347   muonTime.inverseBetaErr=invbetaerr;  
00348   muonTime.timeAtIpInOut = vertexTime;
00349   muonTime.timeAtIpInOutErr = vertexTimeErr;
00350   muonTime.timeAtIpOutIn = vertexTimeR;
00351   muonTime.timeAtIpOutInErr = vertexTimeRErr;
00352       
00353   rawFit(freeBeta, freeBetaErr, freeTime, freeTimeErr, x, y);
00354 
00355   muonTime.freeInverseBeta = freeBeta;    
00356   muonTime.freeInverseBetaErr = freeBetaErr;      
00357     
00358   muonTime.nStations=nStations;
00359 
00360   if (debug) {
00361     cout << " Free 1/beta: " << freeBeta << " +/- " << freeBetaErr << endl;   
00362     cout << "   Free time: " << freeTime << " +/- " << freeTimeErr << endl;   
00363     cout << " Vertex time: " << vertexTime << " +/- " << freeTimeErr << endl;   
00364     cout << " *** Muon Timing Extractor End ***" << endl;
00365   }
00366   
00367   return(muonTime);
00368 }

double MuonTimingExtractor::fitT0	(	double &	a,
		double &	b,
		vector< double >	xl,
		vector< double >	yl,
		vector< double >	xr,
		vector< double >	yr
	)			`[private]`

Definition at line 371 of file MuonTimingExtractor.cc.

References i, s, and ss.

Referenced by fillTiming().

00371                                                                                                                             {
00372 
00373   double sx=0,sy=0,sxy=0,sxx=0,ssx=0,ssy=0,s=0,ss=0;
00374 
00375   for (unsigned int i=0; i<xl.size(); i++) {
00376     sx+=xl[i];
00377     sy+=yl[i];
00378     sxy+=xl[i]*yl[i];
00379     sxx+=xl[i]*xl[i];
00380     s++;
00381     ssx+=xl[i];
00382     ssy+=yl[i];
00383     ss++;
00384   } 
00385 
00386   for (unsigned int i=0; i<xr.size(); i++) {
00387     sx+=xr[i];
00388     sy+=yr[i];
00389     sxy+=xr[i]*yr[i];
00390     sxx+=xr[i]*xr[i];
00391     s++;
00392     ssx-=xr[i];
00393     ssy-=yr[i];
00394     ss--;
00395   } 
00396 
00397   double delta = ss*ss*sxx+s*sx*sx+s*ssx*ssx-s*s*sxx-2*ss*sx*ssx;
00398   
00399   double t0_corr=0.;
00400 
00401   if (delta) {
00402     a=(ssy*s*ssx+sxy*ss*ss+sy*sx*s-sy*ss*ssx-ssy*sx*ss-sxy*s*s)/delta;
00403     b=(ssx*sy*ssx+sxx*ssy*ss+sx*sxy*s-sxx*sy*s-ssx*sxy*ss-sx*ssy*ssx)/delta;
00404     t0_corr=(ssx*s*sxy+sxx*ss*sy+sx*sx*ssy-sxx*s*ssy-sx*ss*sxy-ssx*sx*sy)/delta;
00405   }
00406 
00407   // convert drift distance to time
00408   t0_corr/=-0.00543;
00409 
00410   return t0_corr;
00411 }

void MuonTimingExtractor::rawFit	(	double &	a,
		double &	da,
		double &	b,
		double &	db,
		const vector< double >	hitsx,
		const vector< double >	hitsy
	)			`[private]`

Definition at line 415 of file MuonTimingExtractor.cc.

References d, i, s, funct::sqrt(), x, and y.

Referenced by fillTiming().

00415                                                                                                                                 {
00416 
00417   double s=0,sx=0,sy=0,x,y;
00418   double sxx=0,sxy=0;
00419 
00420   a=b=0;
00421   if (hitsx.size()==0) return;
00422   if (hitsx.size()==1) {
00423     b=hitsy[0];
00424   } else {
00425     for (unsigned int i = 0; i != hitsx.size(); i++) {
00426       x=hitsx[i];
00427       y=hitsy[i];
00428       sy += y;
00429       sxy+= x*y;
00430       s += 1.;
00431       sx += x;
00432       sxx += x*x;
00433     }
00434 
00435     double d = s*sxx - sx*sx;
00436     b = (sxx*sy- sx*sxy)/ d;
00437     a = (s*sxy - sx*sy) / d;
00438     da = sqrt(sxx/d);
00439     db = sqrt(s/d);
00440   }
00441 }

bool MuonTimingExtractor::debug [private]

Definition at line 88 of file MuonTimingExtractor.h.

Referenced by fillTiming().

edm::InputTag MuonTimingExtractor::DTSegmentTags_ [private]

Definition at line 84 of file MuonTimingExtractor.h.

unsigned int MuonTimingExtractor::theHitsMin [private]

Definition at line 85 of file MuonTimingExtractor.h.

Referenced by fillTiming().

MuonSegmentMatcher* MuonTimingExtractor::theMatcher [private]

Definition at line 92 of file MuonTimingExtractor.h.

Referenced by fillTiming(), and MuonTimingExtractor().

double MuonTimingExtractor::thePruneCut [private]

Definition at line 86 of file MuonTimingExtractor.h.

Referenced by fillTiming().

MuonServiceProxy* MuonTimingExtractor::theService [private]

Definition at line 90 of file MuonTimingExtractor.h.

Referenced by fillTiming(), MuonTimingExtractor(), and ~MuonTimingExtractor().

bool MuonTimingExtractor::useSegmentT0 [private]

Definition at line 87 of file MuonTimingExtractor.h.

Referenced by fillTiming().

MuonTimingExtractor Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Classes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
reco::MuonTime	fillTiming (edm::Event &, const edm::EventSetup &, reco::TrackRef muonTrack)
	MuonTimingExtractor (const edm::ParameterSet &)
	Constructor.
	~MuonTimingExtractor ()
	Destructor.
Private Member Functions
double	fitT0 (double &a, double &b, vector< double > xl, vector< double > yl, vector< double > xr, vector< double > yr)
void	rawFit (double &a, double &da, double &b, double &db, const vector< double > hitsx, const vector< double > hitsy)
Private Attributes
bool	debug
edm::InputTag	DTSegmentTags_
unsigned int	theHitsMin
MuonSegmentMatcher *	theMatcher
double	thePruneCut
MuonServiceProxy *	theService
bool	useSegmentT0
Classes
class	TimeMeasurement