DTTimingExtractor Class Reference

#include <RecoMuon/MuonIdentification/src/DTTimingExtractor.cc>

Classes
class	TimeMeasurement

Public Member Functions
	DTTimingExtractor (const edm::ParameterSet &, MuonSegmentMatcher *segMatcher)
	Constructor. More...
void	fillTiming (TimeMeasurementSequence &tmSequence, reco::TrackRef muonTrack, const edm::Event &iEvent, const edm::EventSetup &iSetup)
	~DTTimingExtractor ()
	Destructor. More...

Private Member Functions
double	fitT0 (double &a, double &b, const std::vector< double > &xl, const std::vector< double > &yl, const std::vector< double > &xr, const std::vector< double > &yr)

Private Attributes
bool	debug
bool	doWireCorr_
bool	dropTheta_
edm::InputTag	DTSegmentTags_
bool	requireBothProjections_
double	theError_
unsigned int	theHitsMin_
MuonSegmentMatcher *	theMatcher
double	thePruneCut_
std::unique_ptr< MuonServiceProxy >	theService
double	theTimeOffset_
bool	useSegmentT0_

Detailed Description

Extracts timing information associated to a muon track

Description: Produce timing information for a muon track using 1D DT hits from segments used to build the track

Definition at line 54 of file DTTimingExtractor.h.

Constructor & Destructor Documentation

DTTimingExtractor::DTTimingExtractor	(	const edm::ParameterSet &	iConfig,
		MuonSegmentMatcher *	segMatcher
	)

Constructor.

Definition at line 75 of file DTTimingExtractor.cc.

References edm::ParameterSet::getParameter(), theMatcher, and theService.

  :
  theHitsMin_(iConfig.getParameter<int>("HitsMin")),
  thePruneCut_(iConfig.getParameter<double>("PruneCut")),
  theTimeOffset_(iConfig.getParameter<double>("DTTimeOffset")),
  theError_(iConfig.getParameter<double>("HitError")),
  useSegmentT0_(iConfig.getParameter<bool>("UseSegmentT0")),
  doWireCorr_(iConfig.getParameter<bool>("DoWireCorr")),
  dropTheta_(iConfig.getParameter<bool>("DropTheta")),
  requireBothProjections_(iConfig.getParameter<bool>("RequireBothProjections")),
  debug(iConfig.getParameter<bool>("debug"))
{
  edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
  theService = std::make_unique<MuonServiceProxy>(serviceParameters);
  theMatcher = segMatcher;
}

DTTimingExtractor::~DTTimingExtractor

(

)

Destructor.

Definition at line 93 of file DTTimingExtractor.cc.

{
}

Member Function Documentation

void DTTimingExtractor::fillTiming	(	TimeMeasurementSequence &	tmSequence,
		reco::TrackRef	muonTrack,
		const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

Definition at line 104 of file DTTimingExtractor.cc.

References a, b, gather_cfg::cout, debug, DTTimingExtractor::TimeMeasurement::distIP, doWireCorr_, DTTimingExtractor::TimeMeasurement::driftCell, dropTheta_, TimeMeasurementSequence::dstnc, fitT0(), TrackingRecHit::geographicalId(), edm::EventSetup::get(), mps_fire::i, GlobalTrackingGeometry::idToDet(), reco::if(), DTTimingExtractor::TimeMeasurement::isLeft, DTTimingExtractor::TimeMeasurement::isPhi, DTRecSegment2D::ist0Valid(), DTGeometry::layer(), DTEnums::Left, TimeMeasurementSequence::local_t0, PV3DBase< T, PVType, FrameType >::mag(), mag(), MuonSegmentMatcher::matchDT(), phi, DTTimingExtractor::TimeMeasurement::posInLayer, GloballyPositioned< T >::position(), GeomDet::position(), edm::ESHandle< T >::product(), Propagator::propagateWithPath(), PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi::propagator, requireBothProjections_, DTRecSegment2D::specificRecHits(), mathSSE::sqrt(), DTChamberId::station(), relativeConstraints::station, DTTimingExtractor::TimeMeasurement::station, DTLayer::superLayer(), GeomDet::surface(), DTRecSegment2D::t0(), theError_, theHitsMin_, theMatcher, thePruneCut_, theService, theTimeOffset_, MuonServiceProxy::theTrackingGeometry, DTTimingExtractor::TimeMeasurement::timeCorr, GeomDet::toGlobal(), GeomDet::toLocal(), TimeMeasurementSequence::totalWeightInvbeta, TimeMeasurementSequence::totalWeightTimeVtx, useSegmentT0_, TimeMeasurementSequence::weightInvbeta, TimeMeasurementSequence::weightTimeVtx, PV3DBase< T, PVType, FrameType >::x(), y, and z.

{
  if (debug) 
    std::cout << " *** DT Timimng Extractor ***" << std::endl;

  theService->update(iSetup);

  const GlobalTrackingGeometry *theTrackingGeometry = &*theService->trackingGeometry();

  // get the DT geometry
  edm::ESHandle<DTGeometry> theDTGeom;
  iSetup.get<MuonGeometryRecord>().get(theDTGeom);

  // get the propagator  
  edm::ESHandle<Propagator> propagator;
  iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
  const Propagator *propag = propagator.product();

  math::XYZPoint  pos=muonTrack->innerPosition();
  math::XYZVector mom=muonTrack->innerMomentum();
  GlobalPoint  posp(pos.x(), pos.y(), pos.z());
  GlobalVector momv(mom.x(), mom.y(), mom.z());
  FreeTrajectoryState muonFTS(posp, momv, (TrackCharge)muonTrack->charge(), theService->magneticField().product());

  // get the DT segments that were used to construct the muon
  std::vector<const DTRecSegment4D*> range = theMatcher->matchDT(*muonTrack,iEvent);
  
  if (debug) 
    std::cout << " The muon track matches " << range.size() << " segments." << std::endl;
  
  // create a collection on TimeMeasurements for the track        
  std::vector<TimeMeasurement> tms;
  for (std::vector<const DTRecSegment4D*>::iterator rechit = range.begin(); rechit!=range.end();++rechit) {

    // Create the ChamberId
    DetId id = (*rechit)->geographicalId();
    DTChamberId chamberId(id.rawId());
    int station = chamberId.station();
    if (debug) std::cout << "Matched DT segment in station " << station << std::endl;

    // use only segments with both phi and theta projections present (optional)
    bool bothProjections = ( ((*rechit)->hasPhi()) && ((*rechit)->hasZed()) );
    if (requireBothProjections_ && !bothProjections) continue;

    // loop over (theta, phi) segments
    for (int phi=0; phi<2; phi++) {

      if (dropTheta_ && !phi) continue;

      const DTRecSegment2D* segm;
      if (phi) { 
        segm = dynamic_cast<const DTRecSegment2D*>((*rechit)->phiSegment()); 
        if (debug) std::cout << " *** Segment t0: " << segm->t0() << std::endl;
      }    
        else segm = dynamic_cast<const DTRecSegment2D*>((*rechit)->zSegment());

      if(segm == 0) continue;
      if (!segm->specificRecHits().size()) continue;

      const GeomDet* geomDet = theTrackingGeometry->idToDet(segm->geographicalId());
      const std::vector<DTRecHit1D> hits1d = segm->specificRecHits();

      // store all the hits from the segment
      for (std::vector<DTRecHit1D>::const_iterator hiti=hits1d.begin(); hiti!=hits1d.end(); hiti++) {

    const GeomDet* dtcell = theTrackingGeometry->idToDet(hiti->geographicalId());
    TimeMeasurement thisHit;

    std::pair< TrajectoryStateOnSurface, double> tsos;
    tsos=propag->propagateWithPath(muonFTS,dtcell->surface());

        double dist;            
        double dist_straight = dtcell->toGlobal(hiti->localPosition()).mag(); 
    if (tsos.first.isValid()) { 
      dist = tsos.second+posp.mag(); 
//    std::cout << "Propagate distance: " << dist << " ( innermost: " << posp.mag() << ")" << std::endl; 
    } else { 
      dist = dist_straight;
//    std::cout << "Geom distance: " << dist << std::endl; 
    }

    thisHit.driftCell = hiti->geographicalId();
    if (hiti->lrSide()==DTEnums::Left) thisHit.isLeft=true; else thisHit.isLeft=false;
    thisHit.isPhi = phi;
    thisHit.posInLayer = geomDet->toLocal(dtcell->toGlobal(hiti->localPosition())).x();
    thisHit.distIP = dist;
    thisHit.station = station;
    if (useSegmentT0_ && segm->ist0Valid()) thisHit.timeCorr=segm->t0();
    else thisHit.timeCorr=0.;
    thisHit.timeCorr += theTimeOffset_;
      
    // signal propagation along the wire correction for unmached theta or phi segment hits
    if (doWireCorr_ && !bothProjections && tsos.first.isValid()) {
      const DTLayer* layer = theDTGeom->layer(hiti->wireId());
      float propgL = layer->toLocal( tsos.first.globalPosition() ).y();
      float wirePropCorr = propgL/24.4*0.00543; // signal propagation speed along the wire
      if (thisHit.isLeft) wirePropCorr=-wirePropCorr;
      thisHit.posInLayer += wirePropCorr;
      const DTSuperLayer *sl = layer->superLayer();
      float tofCorr = sl->surface().position().mag()-tsos.first.globalPosition().mag();
      tofCorr = (tofCorr/29.979)*0.00543;
      if (thisHit.isLeft) tofCorr=-tofCorr;
      thisHit.posInLayer += tofCorr;
    } else {
          // non straight-line path correction
          float slCorr = (dist_straight-dist)/29.979*0.00543;
      if (thisHit.isLeft) slCorr=-slCorr;
      thisHit.posInLayer += slCorr;
    }

        if (debug) std::cout << " dist: " << dist << "  t0: " << thisHit.posInLayer << std::endl;
 
    tms.push_back(thisHit);
      }
    } // phi = (0,1)            
  } // rechit
      
  bool modified = false;
  std::vector <double> dstnc, local_t0, hitWeightTimeVtx, hitWeightInvbeta, left;
  double totalWeightInvbeta=0;
  double totalWeightTimeVtx=0;
    
  // Now loop over the measurements, calculate 1/beta and time at vertex and cut away outliers
  do {    

    modified = false;
    dstnc.clear();
    local_t0.clear();
    hitWeightTimeVtx.clear();
    hitWeightInvbeta.clear();
    left.clear();
      
    std::vector <int> hit_idx;
    totalWeightInvbeta=0;
    totalWeightTimeVtx=0;
      
    // Rebuild segments
    for (int sta=1;sta<5;sta++)
      for (int phi=0;phi<2;phi++) {
        std::vector <TimeMeasurement> seg;
        std::vector <int> seg_idx;
    int tmpos=0;
    for (std::vector<TimeMeasurement>::iterator tm=tms.begin(); tm!=tms.end(); ++tm) {
      if ((tm->station==sta) && (tm->isPhi==phi)) {
        seg.push_back(*tm);
        seg_idx.push_back(tmpos);
      }
      tmpos++;  
    }

    unsigned int segsize = seg.size();
    if (segsize<theHitsMin_) continue;

    double a=0, b=0;
        std::vector <double> hitxl,hitxr,hityl,hityr;

    for (std::vector<TimeMeasurement>::iterator tm=seg.begin(); tm!=seg.end(); ++tm) {
 
      DetId id = tm->driftCell;
      const GeomDet* dtcell = theTrackingGeometry->idToDet(id);
      DTChamberId chamberId(id.rawId());
      const GeomDet* dtcham = theTrackingGeometry->idToDet(chamberId);

      double celly=dtcham->toLocal(dtcell->position()).z();
            
      if (tm->isLeft) {
        hitxl.push_back(celly);
        hityl.push_back(tm->posInLayer);
      } else {
        hitxr.push_back(celly);
        hityr.push_back(tm->posInLayer);
      }    
    }

    if (!fitT0(a,b,hitxl,hityl,hitxr,hityr)) {
      if (debug)
        std::cout << "     t0 = zero, Left hits: " << hitxl.size() << " Right hits: " << hitxr.size() << std::endl;
      continue;
    }

    // a segment must have at least one left and one right hit
    if ((!hitxl.size()) || (!hityl.size())) continue;

    int segidx=0;
    for (std::vector<TimeMeasurement>::const_iterator tm=seg.begin(); tm!=seg.end(); ++tm) {

      DetId id = tm->driftCell;
      const GeomDet* dtcell = theTrackingGeometry->idToDet(id);
      DTChamberId chamberId(id.rawId());
      const GeomDet* dtcham = theTrackingGeometry->idToDet(chamberId);

      double layerZ  = dtcham->toLocal(dtcell->position()).z();
      double segmLocalPos = b+layerZ*a;
      double hitLocalPos = tm->posInLayer;
      int hitSide = -tm->isLeft*2+1;
      double t0_segm = (-(hitSide*segmLocalPos)+(hitSide*hitLocalPos))/0.00543+tm->timeCorr;
      
      if (debug) std::cout << "   Segm hit.  dstnc: " << tm->distIP << "   t0: " << t0_segm << std::endl;
            
      dstnc.push_back(tm->distIP);
      local_t0.push_back(t0_segm);
      left.push_back(hitSide);
      hitWeightInvbeta.push_back(((double)seg.size()-2.)*tm->distIP*tm->distIP/((double)seg.size()*30.*30.*theError_*theError_));
          hitWeightTimeVtx.push_back(((double)seg.size()-2.)/((double)seg.size()*theError_*theError_));
      hit_idx.push_back(seg_idx.at(segidx));
      segidx++;
      totalWeightInvbeta+=((double)seg.size()-2.)*tm->distIP*tm->distIP/((double)seg.size()*30.*30.*theError_*theError_);
      totalWeightTimeVtx+=((double)seg.size()-2.)/((double)seg.size()*theError_*theError_);
    }
      }

    if (totalWeightInvbeta==0) break;        

    // calculate the value and error of 1/beta and timeVtx from the complete set of 1D hits
    if (debug)
      std::cout << " Points for global fit: " << dstnc.size() << std::endl;

    double invbeta=0,invbetaErr=0;
    double timeVtx=0,timeVtxErr=0;

    for (unsigned int i=0;i<dstnc.size();i++) {
      invbeta+=(1.+local_t0.at(i)/dstnc.at(i)*30.)*hitWeightInvbeta.at(i)/totalWeightInvbeta;
      timeVtx+=local_t0.at(i)*hitWeightTimeVtx.at(i)/totalWeightTimeVtx;
    }
    
    double chimax=0.;
    std::vector<TimeMeasurement>::iterator tmmax;
    
    // Calculate the inv beta and time at vertex dispersion
    double diff_ibeta,diff_tvtx;
    for (unsigned int i=0;i<dstnc.size();i++) {
      diff_ibeta=(1.+local_t0.at(i)/dstnc.at(i)*30.)-invbeta;
      diff_ibeta=diff_ibeta*diff_ibeta*hitWeightInvbeta.at(i);
      diff_tvtx=local_t0.at(i)-timeVtx;
      diff_tvtx=diff_tvtx*diff_tvtx*hitWeightTimeVtx.at(i);
      invbetaErr+=diff_ibeta;
      timeVtxErr+=diff_tvtx;
      
      // decide if we cut away time at vertex outliers or inverse beta outliers
      // currently not configurable.
      if (diff_tvtx>chimax) { 
    tmmax=tms.begin()+hit_idx.at(i);
    chimax=diff_tvtx;
      }
    }
    
    // cut away the outliers
    if (chimax>thePruneCut_) {
      tms.erase(tmmax);
      modified=true;
    }    

    if (debug) {
      double cf = 1./(dstnc.size()-1);
      invbetaErr=sqrt(invbetaErr/totalWeightInvbeta*cf); 
      timeVtxErr=sqrt(timeVtxErr/totalWeightTimeVtx*cf); 
      std::cout << " Measured 1/beta: " << invbeta << " +/- " << invbetaErr << std::endl;
      std::cout << " Measured time: " << timeVtx << " +/- " << timeVtxErr << std::endl;
    }  

  } while (modified);

  for (unsigned int i=0;i<dstnc.size();i++) {
    tmSequence.dstnc.push_back(dstnc.at(i));
    tmSequence.local_t0.push_back(local_t0.at(i));
    tmSequence.weightInvbeta.push_back(hitWeightInvbeta.at(i));
    tmSequence.weightTimeVtx.push_back(hitWeightTimeVtx.at(i));
  }

  tmSequence.totalWeightInvbeta=totalWeightInvbeta;
  tmSequence.totalWeightTimeVtx=totalWeightTimeVtx;

}

double DTTimingExtractor::fitT0 ( double &  a, double &  b, const std::vector< double > &  xl, const std::vector< double > &  yl, const std::vector< double > &  xr, const std::vector< double > &  yr  )

private

Definition at line 380 of file DTTimingExtractor.cc.

References delta, mps_fire::i, alignCSCRings::s, fftjetcommon_cfi::sx, and fftjetcommon_cfi::sy.

Referenced by fillTiming().

                                                                                                                                                                  {

  double sx=0,sy=0,sxy=0,sxx=0,ssx=0,ssy=0,s=0,ss=0;

  for (unsigned int i=0; i<xl.size(); i++) {
    sx+=xl[i];
    sy+=yl[i];
    sxy+=xl[i]*yl[i];
    sxx+=xl[i]*xl[i];
    s++;
    ssx+=xl[i];
    ssy+=yl[i];
    ss++;
  } 

  for (unsigned int i=0; i<xr.size(); i++) {
    sx+=xr[i];
    sy+=yr[i];
    sxy+=xr[i]*yr[i];
    sxx+=xr[i]*xr[i];
    s++;
    ssx-=xr[i];
    ssy-=yr[i];
    ss--;
  } 

  double delta = ss*ss*sxx+s*sx*sx+s*ssx*ssx-s*s*sxx-2*ss*sx*ssx;
  
  double t0_corr=0.;

  if (delta) {
    a=(ssy*s*ssx+sxy*ss*ss+sy*sx*s-sy*ss*ssx-ssy*sx*ss-sxy*s*s)/delta;
    b=(ssx*sy*ssx+sxx*ssy*ss+sx*sxy*s-sxx*sy*s-ssx*sxy*ss-sx*ssy*ssx)/delta;
    t0_corr=(ssx*s*sxy+sxx*ss*sy+sx*sx*ssy-sxx*s*ssy-sx*ss*sxy-ssx*sx*sy)/delta;
  }

  // convert drift distance to time
  t0_corr/=-0.00543;

  return t0_corr;
}

Member Data Documentation

bool DTTimingExtractor::debug

private

Definition at line 91 of file DTTimingExtractor.h.

Referenced by fillTiming().

bool DTTimingExtractor::doWireCorr_

private

Definition at line 88 of file DTTimingExtractor.h.

Referenced by fillTiming().

bool DTTimingExtractor::dropTheta_

private

Definition at line 89 of file DTTimingExtractor.h.

Referenced by fillTiming().

edm::InputTag DTTimingExtractor::DTSegmentTags_

private

Definition at line 82 of file DTTimingExtractor.h.

bool DTTimingExtractor::requireBothProjections_

private

Definition at line 90 of file DTTimingExtractor.h.

Referenced by fillTiming().

double DTTimingExtractor::theError_

private

Definition at line 86 of file DTTimingExtractor.h.

Referenced by fillTiming().

unsigned int DTTimingExtractor::theHitsMin_

private

Definition at line 83 of file DTTimingExtractor.h.

Referenced by fillTiming().

MuonSegmentMatcher* DTTimingExtractor::theMatcher

private

Definition at line 94 of file DTTimingExtractor.h.

Referenced by DTTimingExtractor(), and fillTiming().

double DTTimingExtractor::thePruneCut_

private

Definition at line 84 of file DTTimingExtractor.h.

Referenced by fillTiming().

std::unique_ptr<MuonServiceProxy> DTTimingExtractor::theService

private

Definition at line 93 of file DTTimingExtractor.h.

Referenced by DTTimingExtractor(), and fillTiming().

double DTTimingExtractor::theTimeOffset_

private

Definition at line 85 of file DTTimingExtractor.h.

Referenced by fillTiming().

bool DTTimingExtractor::useSegmentT0_

private

Definition at line 87 of file DTTimingExtractor.h.

Referenced by fillTiming().