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

Classes
class	TimeMeasurement
Public Member Functions
	CSCTimingExtractor (const edm::ParameterSet &)
	Constructor.
void	fillTiming (TimeMeasurementSequence &tmSequence, reco::TrackRef muonTrack, const edm::Event &iEvent, const edm::EventSetup &iSetup)
	~CSCTimingExtractor ()
	Destructor.
Private Attributes
edm::InputTag	CSCSegmentTags_
bool	debug
unsigned int	theHitsMin_
MuonSegmentMatcher *	theMatcher
double	thePruneCut_
MuonServiceProxy *	theService
double	theTimeOffset_

Detailed Description

Extracts timing information associated to a muon track

Description: <one line="" class="" summary>="">

Definition at line 55 of file CSCTimingExtractor.h.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 70 of file CSCTimingExtractor.cc.

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

  :
  CSCSegmentTags_(iConfig.getParameter<edm::InputTag>("CSCsegments")),
  thePruneCut_(iConfig.getParameter<double>("PruneCut")),
  theTimeOffset_(iConfig.getParameter<double>("CSCTimeOffset")),
  debug(iConfig.getParameter<bool>("debug"))
{
  edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
  theService = new MuonServiceProxy(serviceParameters);
  
  edm::ParameterSet matchParameters = iConfig.getParameter<edm::ParameterSet>("MatchParameters");

  theMatcher = new MuonSegmentMatcher(matchParameters, theService);
}

CSCTimingExtractor::~CSCTimingExtractor ( )

Destructor.

Definition at line 86 of file CSCTimingExtractor.cc.

References theMatcher, and theService.

{
  if (theService) delete theService;
  if (theMatcher) delete theMatcher;
}

Member Function Documentation

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

Definition at line 99 of file CSCTimingExtractor.cc.

References gather_cfg::cout, debug, diffTreeTool::diff, CSCTimingExtractor::TimeMeasurement::distIP, TimeMeasurementSequence::dstnc, edm::EventSetup::get(), i, GlobalTrackingGeometry::idToDet(), TimeMeasurementSequence::local_t0, mag(), MuonServiceProxy::magneticField(), MuonSegmentMatcher::matchCSC(), pos, edm::ESHandle< T >::product(), Propagator::propagateWithPath(), LargeD0_PixelPairStep_cff::propagator, mathSSE::sqrt(), CSCTimingExtractor::TimeMeasurement::station, relativeConstraints::station, CSCDetId::station(), GeomDet::surface(), theMatcher, thePruneCut_, theService, theTimeOffset_, MuonServiceProxy::theTrackingGeometry, CSCTimingExtractor::TimeMeasurement::timeCorr, GeomDet::toGlobal(), TimeMeasurementSequence::totalWeight, MuonServiceProxy::trackingGeometry(), MuonServiceProxy::update(), and TimeMeasurementSequence::weight.

Referenced by MuonTimingFiller::fillTiming().

{

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

  theService->update(iSetup);

  const GlobalTrackingGeometry *theTrackingGeometry = &*theService->trackingGeometry();
  
  edm::ESHandle<Propagator> propagator;
  iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
  const Propagator *propag = propagator.product();

  double invbeta=0;
  double invbetaerr=0;
  int totalWeight=0;
  std::vector<TimeMeasurement> tms;

  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 CSC segments that were used to construct the muon
  std::vector<const CSCSegment*> range = theMatcher->matchCSC(*muonTrack,iEvent);

  // create a collection on TimeMeasurements for the track        
  for (std::vector<const CSCSegment*>::iterator rechit = range.begin(); rechit!=range.end();++rechit) {

    // Create the ChamberId
    DetId id = (*rechit)->geographicalId();
    CSCDetId chamberId(id.rawId());
    int station = chamberId.station();

    if (!(*rechit)->specificRecHits().size()) continue;

    const std::vector<CSCRecHit2D> hits2d = (*rechit)->specificRecHits();

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

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

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

      double dist;            
      if (tsos.first.isValid()) dist = tsos.second+posp.mag(); 
        else dist = cscDet->toGlobal(hiti->localPosition()).mag();

      thisHit.distIP = dist;
      thisHit.station = station;
      thisHit.timeCorr = hiti->tpeak()+theTimeOffset_;
      tms.push_back(thisHit);
      
//      std::cout << " CSC Hit. Dist= " << dist << "    Time= " << thisHit.timeCorr 
//           << "   invBeta= " << (1.+thisHit.timeCorr/dist*30.) << std::endl;
    }

  } // rechit
      
  bool modified = false;
  std::vector <double> dstnc, dsegm, dtraj, hitWeight;

  // Now loop over the measurements, calculate 1/beta and cut away outliers
  do {    

    modified = false;
    dstnc.clear();
    dsegm.clear();
    dtraj.clear();
    hitWeight.clear();
      
    totalWeight=0;
      
        for (std::vector<TimeMeasurement>::iterator tm=tms.begin(); tm!=tms.end(); ++tm) {
          dstnc.push_back(tm->distIP);
          dsegm.push_back(tm->timeCorr);
          hitWeight.push_back(1.);
          totalWeight++;
        }
          
    if (totalWeight==0) break;        

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

    // inverse beta - weighted average of the contributions from individual hits
    invbeta=0;
    for (unsigned int i=0;i<dstnc.size();i++) 
      invbeta+=(1.+dsegm.at(i)/dstnc.at(i)*30.)*hitWeight.at(i)/totalWeight;

    double chimax=0.;
    std::vector<TimeMeasurement>::iterator tmmax;
    
    // the dispersion of inverse beta
    double diff;
    for (unsigned int i=0;i<dstnc.size();i++) {
      diff=(1.+dsegm.at(i)/dstnc.at(i)*30.)-invbeta;
      diff=diff*diff*hitWeight.at(i);
      invbetaerr+=diff;
      if (diff/totalWeight>chimax) { 
        tmmax=tms.begin()+i;
        chimax=diff;
      }
    }
    
    invbetaerr=sqrt(invbetaerr/totalWeight); 
 
    // cut away the outliers
    if (chimax>thePruneCut_) {
      tms.erase(tmmax);
      modified=true;
    }    

    if (debug)
      std::cout << " Measured 1/beta: " << invbeta << " +/- " << invbetaerr << std::endl;

  } while (modified);

  // std::cout << " *** FINAL Measured 1/beta: " << invbeta << " +/- " << invbetaerr << std::endl;

  for (unsigned int i=0;i<dstnc.size();i++) {
    tmSequence.dstnc.push_back(dstnc.at(i));
    tmSequence.local_t0.push_back(dsegm.at(i));
    tmSequence.weight.push_back(hitWeight.at(i));
  }

  tmSequence.totalWeight=totalWeight;

}