#include <CSCSegAlgoDF.h>

Inheritance diagram for CSCSegAlgoDF:

Public Types
typedef std::deque< bool >	BoolContainer

typedef std::vector< const CSCRecHit2D * >	ChamberHitContainer

typedef std::vector< const CSCRecHit2D * > ::const_iterator	ChamberHitContainerCIt

typedef std::vector< int >	LayerIndex
	Typedefs. More...

Public Member Functions
std::vector< CSCSegment >	buildSegments (const ChamberHitContainer &rechits)

	CSCSegAlgoDF (const edm::ParameterSet &ps)
	Constructor. More...

std::vector< CSCSegment >	run (const CSCChamber *aChamber, const ChamberHitContainer &rechits)

virtual	~CSCSegAlgoDF ()
	Destructor. More...

Public Member Functions inherited from CSCSegmentAlgorithm
	CSCSegmentAlgorithm (const edm::ParameterSet &)
	Constructor. More...

virtual std::vector< CSCSegment >	run (const CSCChamber chamber, const std::vector< const CSCRecHit2D > &rechits)=0

virtual	~CSCSegmentAlgorithm ()
	Destructor. More...

Private Member Functions
bool	addHit (const CSCRecHit2D *hit, int layer)

AlgebraicSymMatrix	calculateError (void) const

void	compareProtoSegment (const CSCRecHit2D *h, int layer)

CLHEP::HepMatrix	derivativeMatrix (void) const

void	flagHitsAsUsed (const ChamberHitContainer &rechitsInChamber)

void	flipErrors (AlgebraicSymMatrix &) const

bool	hasHitOnLayer (int layer) const

bool	isHitNearSegment (const CSCRecHit2D *h) const

void	orderSecondSeed (GlobalPoint gp1, const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2, const ChamberHitContainer &rechits, const LayerIndex &layerIndex)

void	pruneFromResidual ()

void	tryAddingHitsToSegment (const ChamberHitContainer &rechitsInChamber, const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2, const LayerIndex &layerIndex)
	Utility functions. More...

void	updateParameters (void)

AlgebraicSymMatrix	weightMatrix (void) const

Private Attributes
float	chi2Max

ChamberHitContainer	closeHits

bool	debug

double	dPhiFineMax

double	dRPhiFineMax

float	maxRatioResidual

int	minHitsForPreClustering

int	minHitsPerSegment

int	minLayersApart

int	muonsPerChamberMax

const std::string	myName

int	nHitsPerClusterIsShower

float	nSigmaFromSegment

CSCSegAlgoPreClustering *	preCluster_

bool	preClustering

double	protoChi2

LocalVector	protoDirection

LocalPoint	protoIntercept

ChamberHitContainer	protoSegment

float	protoSlope_u

float	protoSlope_v

bool	Pruning

ChamberHitContainer	secondSeedHits

CSCSegAlgoShowering *	showering_

float	tanPhiMax

float	tanThetaMax

bool	testSeg

const CSCChamber *	theChamber

BoolContainer	usedHits

Detailed Description

This is a modified version of the SK algorithm for building endcap muon track segments out of the rechit's in a CSCChamber.

A CSCSegment is a RecSegment4D, and is built from CSCRecHit2D objects, each of which is a RecHit2DLocalPos.

This builds segments by first creating proto-segments from at least 3 hits. We intend to try all possible pairs of hits to start segment building. 'All possible' means each hit lies on different layers in the chamber. Once a hit has been assigned to a segment, we don't consider it again, THAT IS, FOR THE FIRST PASS ONLY ! In fact, this is one of the possible flaw with the SK algorithms as it sometimes manages to build segments with the wrong starting points. In the DF algorithm, the endpoints are tested as the best starting points in a 2nd and 3rd loop.

Another difference with the from the SK algorithm is that rechits can be added to proto segments if they fall within n sigmas of the projected track within a given layer. Hence, a cylinder isn't used as in the SK algorimthm, which allows for pseudo 2D hits built from wire or strip only hits to be used in segment reconstruction.

Also, only a certain muonsPerChamberMax maximum number of segments can be produced in the chamber.

Alternative algorithms can be used for the segment building by writing classes like this, and then selecting which one is actually used via the CSCSegmentBuilder.

Author: Dominique Fortin - UCR

Definition at line 47 of file CSCSegAlgoDF.h.

Member Typedef Documentation

typedef std::deque<bool> CSCSegAlgoDF::BoolContainer

Definition at line 57 of file CSCSegAlgoDF.h.

typedef std::vector<const CSCRecHit2D*> CSCSegAlgoDF::ChamberHitContainer

Definition at line 55 of file CSCSegAlgoDF.h.

typedef std::vector<const CSCRecHit2D*>::const_iterator CSCSegAlgoDF::ChamberHitContainerCIt

Definition at line 56 of file CSCSegAlgoDF.h.

typedef std::vector<int> CSCSegAlgoDF::LayerIndex

Typedefs.

Definition at line 54 of file CSCSegAlgoDF.h.

Constructor & Destructor Documentation

CSCSegAlgoDF::CSCSegAlgoDF ( const edm::ParameterSet & ps )

explicit

Constructor.

Definition at line 32 of file CSCSegAlgoDF.cc.

References chi2Max, debug, dPhiFineMax, dRPhiFineMax, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), maxRatioResidual, minHitsForPreClustering, minHitsPerSegment, minLayersApart, nHitsPerClusterIsShower, preCluster_, preClustering, Pruning, showering_, tanPhiMax, and tanThetaMax.

                                                     : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoDF") {
     
   debug                  = ps.getUntrackedParameter<bool>("CSCSegmentDebug");
   minLayersApart         = ps.getParameter<int>("minLayersApart");
   minHitsPerSegment      = ps.getParameter<int>("minHitsPerSegment");
   dRPhiFineMax           = ps.getParameter<double>("dRPhiFineMax");
   dPhiFineMax            = ps.getParameter<double>("dPhiFineMax");
   tanThetaMax            = ps.getParameter<double>("tanThetaMax");
   tanPhiMax              = ps.getParameter<double>("tanPhiMax");    
   chi2Max                = ps.getParameter<double>("chi2Max");  
   preClustering          = ps.getUntrackedParameter<bool>("preClustering");
   minHitsForPreClustering= ps.getParameter<int>("minHitsForPreClustering");
   nHitsPerClusterIsShower= ps.getParameter<int>("nHitsPerClusterIsShower");
   Pruning                = ps.getUntrackedParameter<bool>("Pruning");
   maxRatioResidual       = ps.getParameter<double>("maxRatioResidualPrune");
 
   preCluster_            = new CSCSegAlgoPreClustering( ps );
   showering_             = new CSCSegAlgoShowering( ps );
 }

CSCSegAlgoDF::~CSCSegAlgoDF ( )

virtual

Destructor.

Definition at line 56 of file CSCSegAlgoDF.cc.

References preCluster_, and showering_.

                             {
   delete preCluster_;
   delete showering_;
 }

Member Function Documentation

bool CSCSegAlgoDF::addHit	(	const CSCRecHit2D *	hit,
		int	layer
	)

private

Definition at line 373 of file CSCSegAlgoDF.cc.

References convertSQLiteXML::ok, and protoSegment.

Referenced by compareProtoSegment(), and tryAddingHitsToSegment().

                                                             {
   
   // Return true if hit was added successfully and then parameters are updated.
   // Return false if there is already a hit on the same layer, or insert failed.
   
   bool ok = true;
   
   // Test that we are not trying to add the same hit again
   for ( ChamberHitContainer::const_iterator it = protoSegment.begin(); it != protoSegment.end(); it++ ) 
     if ( aHit == (*it)  ) return false;
   
   protoSegment.push_back(aHit);
 
   return ok;
 }    

std::vector< CSCSegment > CSCSegAlgoDF::buildSegments ( const ChamberHitContainer & rechits )

Build track segments in this chamber (this is where the actual segment-building algorithm hides.)

Definition at line 106 of file CSCSegAlgoDF.cc.

References calculateError(), chi2Max, flagHitsAsUsed(), flipErrors(), i, cuy::ib, CSCChamber::layer(), CSCRecHit2D::localPosition(), minHitsPerSegment, minLayersApart, nHitsPerClusterIsShower, CSCSegment::nRecHits(), GeomDet::position(), preClustering, protoChi2, protoDirection, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, pruneFromResidual(), Pruning, HI_PhotonSkim_cff::rechits, showering_, CSCSegAlgoShowering::showerSeg(), mathSSE::sqrt(), tanPhiMax, tanThetaMax, groupFilesInBlocks::temp, theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), tryAddingHitsToSegment(), csvLumiCalc::unit, usedHits, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by run().

                                                                                      {
 
   ChamberHitContainer rechits = _rechits;
   // Clear buffer for segment vector
   std::vector<CSCSegment> segmentInChamber;
   segmentInChamber.clear();
 
   unsigned nHitInChamber = rechits.size();
   if ( nHitInChamber < 3 ) return segmentInChamber;
 
   LayerIndex layerIndex( nHitInChamber );
 
   unsigned nLayers = 0;
   int old_layer = -1;   
   for ( unsigned int i = 0; i < nHitInChamber; i++ ) {    
     int this_layer = rechits[i]->cscDetId().layer();
     layerIndex[i] = this_layer;
     if ( this_layer != old_layer ) {
       old_layer = this_layer;
       nLayers++;   
     }
   }
   
   if ( nLayers < 3 ) return segmentInChamber;
 
   double z1 = theChamber->layer(1)->position().z();
   double z6 = theChamber->layer(6)->position().z();
   
   if ( z1 > 0. ) {
     if ( z1 > z6 ) { 
       reverse( layerIndex.begin(), layerIndex.end() );
       reverse( rechits.begin(),    rechits.end() );
     }    
   }
   else if ( z1 < 0. ) {
     if ( z1 < z6 ) {
       reverse( layerIndex.begin(), layerIndex.end() );
       reverse( rechits.begin(),    rechits.end() );
     }    
   }
 
   // Showering muon
   if ( preClustering && int(nHitInChamber) > nHitsPerClusterIsShower && nLayers > 2 ) {
     CSCSegment segShower = showering_->showerSeg(theChamber, rechits);
 
     // Make sure have at least 3 hits...
     if ( segShower.nRecHits() < 3 ) return segmentInChamber;
 
     segmentInChamber.push_back(segShower);
 
     return segmentInChamber;
   }
 
 
   // Initialize flags that a given hit has been allocated to a segment
   BoolContainer used_ini(rechits.size(), false);
   usedHits = used_ini;
   
   ChamberHitContainerCIt ib = rechits.begin();
   ChamberHitContainerCIt ie = rechits.end();
     
   // Now Loop over hits within the chamber to find 1st seed for segment building
   for ( ChamberHitContainerCIt i1 = ib; i1 < ie; ++i1 ) {
     if ( usedHits[i1-ib] ) continue;
 
     const CSCRecHit2D* h1 = *i1;
     int layer1 = layerIndex[i1-ib];
     const CSCLayer* l1 = theChamber->layer(layer1);
     GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
     LocalPoint lp1 = theChamber->toLocal(gp1);  
            
     // Loop over hits backward to find 2nd seed for segment building
     for ( ChamberHitContainerCIt i2 = ie-1; i2 > ib; --i2 ) {   
 
       if ( usedHits[i2-ib] ) continue;   // Hit has been used already
 
       int layer2 = layerIndex[i2-ib];
       if ( (layer2 - layer1) < minLayersApart ) continue;
 
       const CSCRecHit2D* h2 = *i2;
       const CSCLayer* l2 = theChamber->layer(layer2);
       GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
       LocalPoint lp2 = theChamber->toLocal(gp2);  
 
       // Clear proto segment so it can be (re)-filled 
       protoSegment.clear();
 
       // localPosition is position of hit wrt layer (so local z = 0)
       protoIntercept = h1->localPosition();
 
       // We want hit wrt chamber (and local z will be != 0)
       float dz = gp2.z()-gp1.z();
       protoSlope_u = (lp2.x() - lp1.x())/dz ;
       protoSlope_v = (lp2.y() - lp1.y())/dz ;    
 
       // Test if entrance angle is roughly pointing towards IP
       if (fabs(protoSlope_v) > tanThetaMax) continue;
       if (fabs(protoSlope_u) > tanPhiMax ) continue;
      
       protoSegment.push_back(h1);
       protoSegment.push_back(h2);
     
       // Try adding hits to proto segment
       tryAddingHitsToSegment(rechits, i1, i2, layerIndex); 
     
       // Check no. of hits on segment to see if segment is large enough
       bool segok = true;
       unsigned iadd = 0;
 
       if (protoSegment.size() < minHitsPerSegment+iadd) segok = false;
   
       if ( Pruning && segok ) pruneFromResidual();
 
       // Check if segment satisfies chi2 requirement
       if (protoChi2 > chi2Max) segok = false;
 
       if ( segok ) {
 
         // Fill segment properties
        
         // Local direction
         double dz   = 1./sqrt(1. + protoSlope_u*protoSlope_u + protoSlope_v*protoSlope_v);
         double dx   = dz * protoSlope_u;
         double dy   = dz * protoSlope_v;
         LocalVector localDir(dx,dy,dz);
       
         // localDir may need sign flip to ensure it points outward from IP  
         double globalZpos    = ( theChamber->toGlobal( protoIntercept ) ).z();
         double globalZdir    = ( theChamber->toGlobal( localDir ) ).z();
         double directionSign = globalZpos * globalZdir;
         protoDirection       = (directionSign * localDir).unit();
 
         // Error matrix
         AlgebraicSymMatrix protoErrors = calculateError();     
         // but reorder components to match what's required by TrackingRecHit interface
         // i.e. slopes first, then positions
         flipErrors( protoErrors );
 
         CSCSegment temp(protoSegment, protoIntercept, protoDirection, protoErrors, protoChi2); 
 
         segmentInChamber.push_back(temp); 
 
         if (nHitInChamber-protoSegment.size() < 3) return segmentInChamber; 
         if (segmentInChamber.size() > 4) return segmentInChamber;
 
         // Flag used hits
         flagHitsAsUsed(rechits);
       } 
     } 
   }
   return segmentInChamber;
 }

AlgebraicSymMatrix CSCSegAlgoDF::calculateError ( void ) const

private

Definition at line 633 of file CSCSegAlgoDF.cc.

References funct::A, derivativeMatrix(), query::result, weightMatrix(), and create_public_pileup_plots::weights.

Referenced by buildSegments().

                                                       {
 
   AlgebraicSymMatrix weights = weightMatrix();
   AlgebraicMatrix A = derivativeMatrix();
 
   // (AT W A)^-1
   // from http://www.phys.ufl.edu/~avery/fitting.html, part I
   int ierr;
   AlgebraicSymMatrix result = weights.similarityT(A);
   result.invert(ierr);
 
   // blithely assuming the inverting never fails...
   return result;
 }

void CSCSegAlgoDF::compareProtoSegment	(	const CSCRecHit2D *	h,
		int	layer
	)

private

Definition at line 519 of file CSCSegAlgoDF.cc.

References addHit(), convertSQLiteXML::ok, protoChi2, protoDirection, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, and updateParameters().

Referenced by tryAddingHitsToSegment().

                                                                       {
   
   // Store old segment first
   double old_protoChi2                  = protoChi2;
   LocalPoint old_protoIntercept         = protoIntercept;
   float old_protoSlope_u                = protoSlope_u;
   float old_protoSlope_v                = protoSlope_v;
   LocalVector old_protoDirection        = protoDirection;
   ChamberHitContainer old_protoSegment  = protoSegment;
  
 
   // Try adding the hit to existing segment, and remove old one existing in same layer
   ChamberHitContainer::iterator it;
   for ( it = protoSegment.begin(); it != protoSegment.end(); ) {
     if ( (*it)->cscDetId().layer() == layer ) {
       it = protoSegment.erase(it);
     } else {
       ++it;
     }
   }
   bool ok = addHit(h, layer);
 
   if (ok) updateParameters();
   
   if ( (protoChi2 > old_protoChi2) || ( !ok ) ) {
     protoChi2       = old_protoChi2;
     protoIntercept  = old_protoIntercept;
     protoSlope_u    = old_protoSlope_u;
     protoSlope_v    = old_protoSlope_v;
     protoDirection  = old_protoDirection;
     protoSegment    = old_protoSegment;
   }
 }

CLHEP::HepMatrix CSCSegAlgoDF::derivativeMatrix ( void ) const

private

Definition at line 606 of file CSCSegAlgoDF.cc.

References CSCRecHit2D::cscDetId(), CSCChamber::layer(), CSCDetId::layer(), CSCRecHit2D::localPosition(), makeMuonMisalignmentScenario::matrix, protoSegment, theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by calculateError().

                                                     {
 
   ChamberHitContainer::const_iterator it;
   int nhits = protoSegment.size();
   CLHEP::HepMatrix matrix(2*nhits, 4);
   int row = 0;
 
   for(it = protoSegment.begin(); it != protoSegment.end(); ++it) {
 
     const CSCRecHit2D& hit = (**it);
     const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer());
     GlobalPoint gp = layer->toGlobal(hit.localPosition());
     LocalPoint lp = theChamber->toLocal(gp);
     float z = lp.z();
     ++row;
     matrix(row, 1) = 1.;
     matrix(row, 3) = z;
     ++row;
     matrix(row, 2) = 1.;
     matrix(row, 4) = z;
   }
   return matrix;
 }

void CSCSegAlgoDF::flagHitsAsUsed ( const ChamberHitContainer & rechitsInChamber )

private

Flag hits on segment as used

Definition at line 559 of file CSCSegAlgoDF.cc.

References closeHits, cuy::ib, protoSegment, and usedHits.

Referenced by buildSegments().

                                                                              {
   
   // Flag hits on segment as used
   ChamberHitContainerCIt ib = rechitsInChamber.begin();
   ChamberHitContainerCIt hi, iu;
   
   for ( hi = protoSegment.begin(); hi != protoSegment.end(); ++hi ) {
     for ( iu = ib; iu != rechitsInChamber.end(); ++iu ) {
       if (*hi == *iu) usedHits[iu-ib] = true;
     }
   }
   // Don't reject hits marked as "nearby" for now.
   // So this is bypassed at all times for now !!!
   // Perhaps add back to speed up algorithm some more
   if (closeHits.size() > 0) return;  
   for ( hi = closeHits.begin(); hi != closeHits.end(); ++hi ) {
     for ( iu = ib; iu != rechitsInChamber.end(); ++iu ) {
       if (*hi == *iu) usedHits[iu-ib] = true;
     }
   }
 
 }

void CSCSegAlgoDF::flipErrors ( AlgebraicSymMatrix & a ) const

private

Definition at line 648 of file CSCSegAlgoDF.cc.

References a.

Referenced by buildSegments().

                                                            {
 
   // The CSCSegment needs the error matrix re-arranged
 
   AlgebraicSymMatrix hold( a );
 
   // errors on slopes into upper left
   a(1,1) = hold(3,3);
   a(1,2) = hold(3,4);
   a(2,1) = hold(4,3);
   a(2,2) = hold(4,4);
 
   // errors on positions into lower right
   a(3,3) = hold(1,1);
   a(3,4) = hold(1,2);
   a(4,3) = hold(2,1);
   a(4,4) = hold(2,2);
 
   // off-diagonal elements remain unchanged
 
 }

bool CSCSegAlgoDF::hasHitOnLayer ( int layer ) const

private

Definition at line 503 of file CSCSegAlgoDF.cc.

References protoSegment.

Referenced by tryAddingHitsToSegment().

                                                 {
   
   // Is there already a hit on this layer?
   for ( ChamberHitContainerCIt it = protoSegment.begin(); it != protoSegment.end(); it++ )
     if ( (*it)->cscDetId().layer() == layer ) return true;
   
   return false;
 }

bool CSCSegAlgoDF::isHitNearSegment ( const CSCRecHit2D * h ) const

private

Definition at line 336 of file CSCSegAlgoDF.cc.

References CSCRecHit2D::cscDetId(), SiPixelRawToDigiRegional_cfi::deltaPhi, dPhiFineMax, dRPhiFineMax, CSCChamber::layer(), CSCDetId::layer(), CSCRecHit2D::localPosition(), M_PI, PV3DBase< T, PVType, FrameType >::phi(), protoIntercept, protoSlope_u, protoSlope_v, dttmaxenums::R, mathSSE::sqrt(), theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by tryAddingHitsToSegment().

                                                                  {
 
   const CSCLayer* layer = theChamber->layer(hit->cscDetId().layer());
 
   // hit phi position in global coordinates
   GlobalPoint Hgp = layer->toGlobal(hit->localPosition());
   double Hphi = Hgp.phi();                                
   if (Hphi < 0.) Hphi += 2.*M_PI;
   LocalPoint Hlp = theChamber->toLocal(Hgp);
   double z = Hlp.z();  
 
   double LocalX = protoIntercept.x() + protoSlope_u * z;
   double LocalY = protoIntercept.y() + protoSlope_v * z;
   LocalPoint Slp(LocalX, LocalY, z);
   GlobalPoint Sgp = theChamber->toGlobal(Slp); 
   double Sphi = Sgp.phi();
   if (Sphi < 0.) Sphi += 2.*M_PI;
   double R = sqrt(Sgp.x()*Sgp.x() + Sgp.y()*Sgp.y());
   
   double deltaPhi = Sphi - Hphi;
   if (deltaPhi >  2.*M_PI) deltaPhi -= 2.*M_PI;
   if (deltaPhi < -2.*M_PI) deltaPhi += 2.*M_PI;
   if (deltaPhi < 0.) deltaPhi = -deltaPhi; 
 
   double RdeltaPhi = R * deltaPhi;
 
   if (RdeltaPhi < dRPhiFineMax && deltaPhi < dPhiFineMax ) return true;
 
   return false;
 }

void CSCSegAlgoDF::orderSecondSeed	(	GlobalPoint	gp1,
		const ChamberHitContainerCIt	i1,
		const ChamberHitContainerCIt	i2,
		const ChamberHitContainer &	rechits,
		const LayerIndex &	layerIndex
	)

private

Order the hits on the 2nd layer for seed building

Definition at line 744 of file CSCSegAlgoDF.cc.

References secondSeedHits.

                                                                           {
 
   secondSeedHits.clear();
 
   //ChamberHitContainerCIt ib = rechits.begin();
   ChamberHitContainerCIt ie = rechits.end();
 
   //  int layer1 = layerIndex[i1-ib];
   //  int layer2 = layerIndex[i2-ib];
 
   // Now fill vector of rechits closest to center of mass:
   // secondSeedHitsIdx.clear() = 0;
 
   // Loop over all hits and find hit closest to 1st seed.
   for ( ChamberHitContainerCIt i2 = ie-1; i2 > i1; --i2 ) { 
 
   }
 
 }

void CSCSegAlgoDF::pruneFromResidual ( )

private

Prune bad segment from the worse hit based on residuals

Definition at line 672 of file CSCSegAlgoDF.cc.

References CSCRecHit2D::cscDetId(), j, CSCChamber::layer(), CSCDetId::layer(), CSCRecHit2D::localPosition(), maxRatioResidual, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, mathSSE::sqrt(), theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), updateParameters(), findQualityFiles::v, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by buildSegments().

                                     {
 
   // Only prune if have at least 5 hits 
   if ( protoSegment.size() < 5 ) return ;
 
 
   // Now Study residuals
       
   float maxResidual = 0.;
   float sumResidual = 0.;
   int nHits = 0;
   int badIndex = -1;
   int j = 0;
 
 
   ChamberHitContainer::const_iterator ih;
 
   for ( ih = protoSegment.begin(); ih != protoSegment.end(); ++ih ) {
     const CSCRecHit2D& hit = (**ih);
     const CSCLayer* layer  = theChamber->layer(hit.cscDetId().layer());
     GlobalPoint gp         = layer->toGlobal(hit.localPosition());
     LocalPoint lp          = theChamber->toLocal(gp);
 
     double u = lp.x();
     double v = lp.y();
     double z = lp.z();
 
     double du = protoIntercept.x() + protoSlope_u * z - u;
     double dv = protoIntercept.y() + protoSlope_v * z - v;
 
     float residual = sqrt(du*du + dv*dv);
 
     sumResidual += residual;
     nHits++;
     if ( residual > maxResidual ) {
       maxResidual = residual;
       badIndex = j;
     }
     j++;
   }
 
   float corrAvgResidual = (sumResidual - maxResidual)/(nHits -1);
 
   // Keep all hits 
   if ( maxResidual/corrAvgResidual < maxRatioResidual ) return;
 
 
   // Drop worse hit and recompute segment properties + fill
 
   ChamberHitContainer newProtoSegment;
 
   j = 0;
   for ( ih = protoSegment.begin(); ih != protoSegment.end(); ++ih ) {
     if ( j != badIndex ) newProtoSegment.push_back(*ih);
     j++;
   }
   
   protoSegment.clear();
 
   for ( ih = newProtoSegment.begin(); ih != newProtoSegment.end(); ++ih ) {
     protoSegment.push_back(*ih);
   }
 
   // Update segment parameters
   updateParameters();
 
 }

std::vector< CSCSegment > CSCSegAlgoDF::run	(	const CSCChamber *	aChamber,
		const ChamberHitContainer &	rechits
	)

Here we must implement the algorithm

Definition at line 65 of file CSCSegAlgoDF.cc.

References buildSegments(), CSCSegAlgoPreClustering::clusterHits(), minHitsForPreClustering, preCluster_, preClustering, and theChamber.

                                                                                                       {
 
   // Store chamber info in temp memory
   theChamber = aChamber; 
 
   int nHits = rechits.size();
 
   // Segments prior to pruning
   std::vector<CSCSegment> segments_temp;  
 
   if ( preClustering && nHits > minHitsForPreClustering ) {
     // This is where the segment origin is in the chamber on avg.
     std::vector<CSCSegment> testSegments;
     std::vector<ChamberHitContainer> clusteredHits = preCluster_->clusterHits(theChamber, rechits);
     // loop over the found clusters:
     for (std::vector<ChamberHitContainer>::iterator subrechits = clusteredHits.begin(); subrechits != clusteredHits.end(); ++subrechits ) {
       // build the subset of segments:
       std::vector<CSCSegment> segs = buildSegments( (*subrechits) );
       // add the found subset of segments to the collection of all segments in this chamber:
       segments_temp.insert( segments_temp.end(), segs.begin(), segs.end() );
     }
   } else {
     std::vector<CSCSegment> segs = buildSegments( rechits );
     // add the found subset of segments to the collection of all segments in this chamber:
     segments_temp.insert( segments_temp.end(), segs.begin(), segs.end() ); 
   }
 
   return segments_temp; 
 }

void CSCSegAlgoDF::tryAddingHitsToSegment	(	const ChamberHitContainer &	rechitsInChamber,
		const ChamberHitContainerCIt	i1,
		const ChamberHitContainerCIt	i2,
		const LayerIndex &	layerIndex
	)

private

Utility functions.

Try adding non-used hits to segment
Skip the layers containing the segment endpoints on first 2 passes, but then
try hits on layer containing the segment starting points on 2nd and/or 3rd pass
if segment has >2 hits. Test each hit on the other layers to see if it is near
the segment using rechit error matrix.
If it is, see whether there is already a hit on the segment from the same layer

if so, and there are more than 2 hits on the segment, copy the segment,
replace the old hit with the new hit. If the new segment chi2 is better
then replace the original segment with the new one
if not, copy the segment, add the hit if it's within a certain range.

Definition at line 266 of file CSCSegAlgoDF.cc.

References addHit(), closeHits, compareProtoSegment(), h, hasHitOnLayer(), i, cuy::ib, isHitNearSegment(), protoSegment, updateParameters(), and usedHits.

Referenced by buildSegments().

                                                                           {
   
 /* Iterate over the layers with hits in the chamber
  * Skip the layers containing the segment endpoints on first pass, but then
  * try hits on layer containing the segment starting points on 2nd pass
  * if segment has >2 hits.  Once a hit is added to a layer, don't replace it 
  * until 2nd iteration
  */  
   
   ChamberHitContainerCIt ib = rechits.begin();
   ChamberHitContainerCIt ie = rechits.end();
   closeHits.clear();
 
   for ( ChamberHitContainerCIt i = ib; i != ie; ++i ) {
 
     if (i == i1 || i == i2 ) continue;
     if ( usedHits[i-ib] ) continue;   // Don't use hits already part of a segment.
 
     const CSCRecHit2D* h = *i;
     int layer = layerIndex[i-ib];
     int layer1 = layerIndex[i1-ib];
     int layer2 = layerIndex[i2-ib];
 
     // Low multiplicity case
     if (rechits.size() < 9) {
       if ( isHitNearSegment( h ) ) {
         if ( !hasHitOnLayer(layer) ) {
           addHit(h, layer);
         } else {
           closeHits.push_back(h);
         }
       }
 
     // High multiplicity case
     } else { 
       if ( isHitNearSegment( h ) ) {
         if ( !hasHitOnLayer(layer) ) {
           addHit(h, layer);
           updateParameters();
         // Don't change the starting points at this stage !!!
         } else {
           closeHits.push_back(h);
           if (layer != layer1 && layer != layer2 ) compareProtoSegment(h, layer);
         }
       }
     }
   }
  
   if ( int(protoSegment.size()) < 3) return;
 
   updateParameters();
 
   // 2nd pass to remove biases 
   // This time, also consider changing the endpoints
   for ( ChamberHitContainerCIt i = closeHits.begin() ; i != closeHits.end(); ++i ) {      
     const CSCRecHit2D* h = *i;      
     int layer = (*i)->cscDetId().layer();     
     compareProtoSegment(h, layer); 
   } 
 
 }

void CSCSegAlgoDF::updateParameters ( void )

private

Definition at line 395 of file CSCSegAlgoDF.cc.

References CSCRecHit2D::cscDetId(), CSCChamber::layer(), CSCDetId::layer(), CSCRecHit2D::localPosition(), CSCRecHit2D::localPositionError(), LogDebug, AlCaHLTBitMon_ParallelJobs::p, protoChi2, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), findQualityFiles::v, PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), LocalError::xy(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by compareProtoSegment(), pruneFromResidual(), and tryAddingHitsToSegment().

                                     {
 
   // Compute slope from Least Square Fit    
   CLHEP::HepMatrix M(4,4,0);
   CLHEP::HepVector B(4,0);
 
   ChamberHitContainer::const_iterator ih;
   
   for (ih = protoSegment.begin(); ih != protoSegment.end(); ++ih) {
     
     const CSCRecHit2D& hit = (**ih);
     const CSCLayer* layer  = theChamber->layer(hit.cscDetId().layer());
     GlobalPoint gp         = layer->toGlobal(hit.localPosition());
     LocalPoint  lp         = theChamber->toLocal(gp); 
     
     double u = lp.x();
     double v = lp.y();
     double z = lp.z();
     
     // ptc: Covariance matrix of local errors 
     CLHEP::HepMatrix IC(2,2);
     IC(1,1) = hit.localPositionError().xx();
     IC(1,2) = hit.localPositionError().xy();
     IC(2,2) = hit.localPositionError().yy();
     IC(2,1) = IC(1,2); // since Cov is symmetric
     
     // ptc: Invert covariance matrix (and trap if it fails!)
     int ierr = 0;
     IC.invert(ierr); // inverts in place
     if (ierr != 0) {
       LogDebug("CSC") << "CSCSegment::fitSlopes: failed to invert covariance matrix=\n" << IC << "\n";      
     }
     
     M(1,1) += IC(1,1);
     M(1,2) += IC(1,2);
     M(1,3) += IC(1,1) * z;
     M(1,4) += IC(1,2) * z;
     B(1)   += u * IC(1,1) + v * IC(1,2);
     
     M(2,1) += IC(2,1);
     M(2,2) += IC(2,2);
     M(2,3) += IC(2,1) * z;
     M(2,4) += IC(2,2) * z;
     B(2)   += u * IC(2,1) + v * IC(2,2);
     
     M(3,1) += IC(1,1) * z;
     M(3,2) += IC(1,2) * z;
     M(3,3) += IC(1,1) * z * z;
     M(3,4) += IC(1,2) * z * z;
     B(3)   += ( u * IC(1,1) + v * IC(1,2) ) * z;
     
     M(4,1) += IC(2,1) * z;
     M(4,2) += IC(2,2) * z;
     M(4,3) += IC(2,1) * z * z;
     M(4,4) += IC(2,2) * z * z;
     B(4)   += ( u * IC(2,1) + v * IC(2,2) ) * z;
   }
   
   CLHEP::HepVector p = solve(M, B);
   
   // Update member variables 
   // Note that origin has local z = 0
 
   protoIntercept = LocalPoint(p(1), p(2), 0.);
   protoSlope_u = p(3);
   protoSlope_v = p(4);
 
 
   // Determine Chi^2 for the proto wire segment
   
   double chsq = 0.;
   
   for (ih = protoSegment.begin(); ih != protoSegment.end(); ++ih) {
     
     const CSCRecHit2D& hit = (**ih);
     const CSCLayer* layer  = theChamber->layer(hit.cscDetId().layer());
     GlobalPoint gp         = layer->toGlobal(hit.localPosition());
     LocalPoint lp          = theChamber->toLocal(gp);
     
     double u = lp.x();
     double v = lp.y();
     double z = lp.z();
     
     double du = protoIntercept.x() + protoSlope_u * z - u;
     double dv = protoIntercept.y() + protoSlope_v * z - v;
     
     CLHEP::HepMatrix IC(2,2);
     IC(1,1) = hit.localPositionError().xx();
     IC(1,2) = hit.localPositionError().xy();
     IC(2,2) = hit.localPositionError().yy();
     IC(2,1) = IC(1,2);
     
     // Invert covariance matrix
     int ierr = 0;
     IC.invert(ierr);
     if (ierr != 0) {
       LogDebug("CSC") << "CSCSegment::fillChiSquared: failed to invert covariance matrix=\n" << IC << "\n";      
     }
     chsq += du*du*IC(1,1) + 2.*du*dv*IC(1,2) + dv*dv*IC(2,2);
   }
   protoChi2 = chsq;
 }

AlgebraicSymMatrix CSCSegAlgoDF::weightMatrix ( void ) const

private

Definition at line 583 of file CSCSegAlgoDF.cc.

References CSCRecHit2D::localPositionError(), makeMuonMisalignmentScenario::matrix, protoSegment, LocalError::xx(), LocalError::xy(), and LocalError::yy().

Referenced by calculateError().

                                                     {
 
   std::vector<const CSCRecHit2D*>::const_iterator it;
   int nhits = protoSegment.size();
   AlgebraicSymMatrix matrix(2*nhits, 0);
   int row = 0;
 
   for (it = protoSegment.begin(); it != protoSegment.end(); ++it) {
 
     const CSCRecHit2D& hit = (**it);
     ++row;
     matrix(row, row)   = hit.localPositionError().xx();
     matrix(row, row+1) = hit.localPositionError().xy();
     ++row;
     matrix(row, row-1) = hit.localPositionError().xy();
     matrix(row, row)   = hit.localPositionError().yy();
   }
   int ierr;
   matrix.invert(ierr);
   return matrix;
 }