---

CSCSegAlgoHitPruning Class Reference

#include <RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.h>

List of all members.

Public Types
typedef std::vector< const CSCRecHit2D * >	ChamberHitContainer
Public Member Functions
	CSCSegAlgoHitPruning (const edm::ParameterSet &ps)
	constructor
std::vector< CSCSegment >	pruneBadHits (const CSCChamber *aChamber, std::vector< CSCSegment > segments)
	clusterize
	~CSCSegAlgoHitPruning ()
	destructor
Private Member Functions
AlgebraicSymMatrix	calculateError (void) const
HepMatrix	derivativeMatrix (void) const
void	fillChiSquared (void)
void	fillLocalDirection (void)
void	fitSlopes (void)
void	flipErrors (AlgebraicSymMatrix &) const
AlgebraicSymMatrix	weightMatrix (void) const
Private Attributes
bool	BrutePruning
double	protoChi2
LocalVector	protoDirection
LocalPoint	protoIntercept
ChamberHitContainer	protoSegment
float	protoSlope_u
float	protoSlope_v
const CSCChamber *	theChamber

---

Detailed Description

Definition at line 23 of file CSCSegAlgoHitPruning.h.

---

Member Typedef Documentation

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

Definition at line 27 of file CSCSegAlgoHitPruning.h.

---

Constructor & Destructor Documentation

CSCSegAlgoHitPruning::CSCSegAlgoHitPruning

(

const edm::ParameterSet &

ps

)

[explicit]

constructor

Definition at line 32 of file CSCSegAlgoHitPruning.cc.

References BrutePruning, and edm::ParameterSet::getUntrackedParameter().

                                                                    {
  BrutePruning           = ps.getUntrackedParameter<bool>("BrutePruning");

}

CSCSegAlgoHitPruning::~CSCSegAlgoHitPruning

(

)

destructor

Definition at line 41 of file CSCSegAlgoHitPruning.cc.

                                           {

}

---

Member Function Documentation

AlgebraicSymMatrix CSCSegAlgoHitPruning::calculateError

(

void

)

const [private]

Definition at line 421 of file CSCSegAlgoHitPruning.cc.

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

Referenced by pruneBadHits().

                                                              {
  
  AlgebraicSymMatrix weights = weightMatrix();
  AlgebraicMatrix A = derivativeMatrix();
  
  // (AT W A)^-1
  // from https://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;
}

HepMatrix CSCSegAlgoHitPruning::derivativeMatrix

(

void

)

const [private]

Definition at line 393 of file CSCSegAlgoHitPruning.cc.

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

Referenced by calculateError().

                                                       {
  
  ChamberHitContainer::const_iterator it;
  int nhits = protoSegment.size();
  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 CSCSegAlgoHitPruning::fillChiSquared

(

void

)

[private]

Definition at line 298 of file CSCSegAlgoHitPruning.cc.

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

Referenced by pruneBadHits().

                                          {
  
  double chsq = 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;
    
    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";
//       std::cout << "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;
}

void CSCSegAlgoHitPruning::fillLocalDirection

(

void

)

[private]

Definition at line 342 of file CSCSegAlgoHitPruning.cc.

References protoDirection, protoIntercept, protoSlope_u, protoSlope_v, funct::sqrt(), theChamber, GeomDet::toGlobal(), and z.

Referenced by pruneBadHits().

                                              {
  // Always enforce direction of segment to point from IP outwards
  // (Incorrect for particles not coming from IP, of course.)
  
  double dxdz = protoSlope_u;
  double dydz = protoSlope_v;
  double dz   = 1./sqrt(1. + dxdz*dxdz + dydz*dydz);
  double dx   = dz*dxdz;
  double dy   = dz*dydz;
  LocalVector localDir(dx,dy,dz);
  
  // localDir may need sign flip to ensure it points outward from IP
  // ptc: Examine its direction and origin in global z: to point outward
  // the localDir should always have same sign as global z...
  
  double globalZpos    = ( theChamber->toGlobal( protoIntercept ) ).z();
  double globalZdir    = ( theChamber->toGlobal( localDir ) ).z();
  double directionSign = globalZpos * globalZdir;
  protoDirection       = (directionSign * localDir).unit();
}

void CSCSegAlgoHitPruning::fitSlopes

(

void

)

[private]

Definition at line 232 of file CSCSegAlgoHitPruning.cc.

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

Referenced by pruneBadHits().

                                     {
  HepMatrix M(4,4,0);
  HepVector B(4,0);
  ChamberHitContainer::const_iterator ih = protoSegment.begin();
  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); 
    // ptc: Local position of hit w.r.t. chamber
    double u = lp.x();
    double v = lp.y();
    double z = lp.z();
    // ptc: Covariance matrix of local errors 
    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";      
//       std::cout<< "CSCSegment::fitSlopes: failed to invert covariance matrix=\n" << IC << "\n"<<std::endl;
    }
    
    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;
  }
  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);
}

void CSCSegAlgoHitPruning::flipErrors

(

AlgebraicSymMatrix &

a

)

const [private]

Definition at line 437 of file CSCSegAlgoHitPruning.cc.

References a.

Referenced by pruneBadHits().

                                                                   { 
    
  // 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 
    
} 

std::vector< CSCSegment > CSCSegAlgoHitPruning::pruneBadHits	(	const CSCChamber *	aChamber,
		std::vector< CSCSegment >	segments
	)

clusterize

Definition at line 49 of file CSCSegAlgoHitPruning.cc.

References begin, BrutePruning, calculateError(), CSCSegment::chi2(), ChiSquaredProbability(), e, fillChiSquared(), fillLocalDirection(), fitSlopes(), flipErrors(), it, CSCChamber::layer(), m, CSCSegment::nRecHits(), protoChi2, protoDirection, protoIntercept, protoSegment, radius(), pyDBSRunClass::temp, theChamber, GeomDet::toGlobal(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                                                                   {

  theChamber = aChamber;

  std::vector<CSCSegment>          segments_temp;
  std::vector<ChamberHitContainer> rechits_clusters; 
  
  const float chi2ndfProbMin = 1.0e-4;
  bool use_brute_force = BrutePruning;

  int hit_nr = 0;
  int hit_nr_worst = -1;
  int hit_nr_2ndworst = -1;
  
  for (std::vector<CSCSegment>::iterator it=segments.begin(); it != segments.end(); it++) {
    
    if ( !use_brute_force ) {// find worst hit
      
      float chisq    = (*it).chi2();
      int nhits      = (*it).nRecHits();
      LocalPoint localPos = (*it).localPosition();
      LocalVector segDir = (*it).localDirection();
      const CSCChamber* cscchamber = theChamber;
      float globZ       ;
          
      GlobalPoint globalPosition = cscchamber->toGlobal(localPos);
      globZ = globalPosition.z();
      
      
      if ( ChiSquaredProbability((double)chisq,(double)(2*nhits-4)) < chi2ndfProbMin  ) {

        // find (rough) "residuals" (NOT excluding the hit from the fit - speed!) of hits on segment
        std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
        std::vector<CSCRecHit2D>::const_iterator iRH_worst;
        float xdist_local       = -99999.;

        float xdist_local_worst_sig = -99999.;
        float xdist_local_2ndworst_sig = -99999.;
        float xdist_local_sig       = -99999.;

        hit_nr = 0;
        hit_nr_worst = -1;
        hit_nr_2ndworst = -1;

        for ( std::vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++ ) {
          //mark "worst" hit:
          
          float z_at_target ;
          float radius      ;
          float loc_x_at_target ;
          float loc_y_at_target ;
          float loc_z_at_target ;

          z_at_target  = 0.;
          loc_x_at_target  = 0.;
          loc_y_at_target  = 0.;
          loc_z_at_target  = 0.;
          radius       = 0.;
          
          // set the z target in CMS global coordinates:
          const CSCLayer* csclayerRH = theChamber->layer((*iRH).cscDetId().layer());
          LocalPoint localPositionRH = (*iRH).localPosition();
          GlobalPoint globalPositionRH = csclayerRH->toGlobal(localPositionRH); 
          
          LocalError rerrlocal = (*iRH).localPositionError();  
          float xxerr = rerrlocal.xx();
          
          float target_z     = globalPositionRH.z();  // target z position in cm (z pos of the hit)
          
          loc_x_at_target = localPos.x() + (segDir.x()*( target_z - globZ ));
          loc_y_at_target = localPos.y() + (segDir.y()*( target_z - globZ ));
          loc_z_at_target = target_z;

          // have to transform the segments coordinates back to the local frame... how?!!!!!!!!!!!!
          
          xdist_local  = fabs(localPositionRH.x() - loc_x_at_target);
          xdist_local_sig  = fabs((localPositionRH.x() -loc_x_at_target)/(xxerr));
          
          if( xdist_local_sig > xdist_local_worst_sig ) {
            xdist_local_2ndworst_sig = xdist_local_worst_sig;
            xdist_local_worst_sig    = xdist_local_sig;
            iRH_worst            = iRH;
            hit_nr_2ndworst = hit_nr_worst;
            hit_nr_worst = hit_nr;
          }
          else if(xdist_local_sig > xdist_local_2ndworst_sig) {
            xdist_local_2ndworst_sig = xdist_local_sig;
            hit_nr_2ndworst = hit_nr;
          }
          ++hit_nr;
        }

        // reset worst hit number if certain criteria apply.
        // Criteria: 2nd worst hit must be at least a factor of
        // 1.5 better than the worst in terms of sigma:
        if ( xdist_local_worst_sig / xdist_local_2ndworst_sig < 1.5 ) {
          hit_nr_worst    = -1;
          hit_nr_2ndworst = -1;
        }
      }
    }

    // if worst hit was found, refit without worst hit and select if considerably better than original fit.
    // Can also use brute force: refit all n-1 hit segments and choose one over the n hit if considerably "better"
   
      std::vector< CSCRecHit2D > buffer;
      std::vector< std::vector< CSCRecHit2D > > reduced_segments;
      std::vector< CSCRecHit2D > theseRecHits = (*it).specificRecHits();
      float best_red_seg_prob = 0.0;
      // usefor chi2 1 diff   float best_red_seg_prob = 99999.;
      buffer.clear();
      if( ChiSquaredProbability((double)(*it).chi2(),(double)((2*(*it).nRecHits())-4)) < chi2ndfProbMin  ) {
        
        buffer = theseRecHits;

        // Dirty switch: here one can select to refit all possible subsets or just the one without the 
        // tagged worst hit:
        if( use_brute_force ) { // Brute force method: loop over all possible segments:
          for(uint bi = 0; bi < buffer.size(); bi++) {
            reduced_segments.push_back(buffer);
            reduced_segments[bi].erase(reduced_segments[bi].begin()+(bi),reduced_segments[bi].begin()+(bi+1));
          }
        }
        else { // More elegant but still biased: erase only worst hit
          // Comment: There is not a very strong correlation of the worst hit with the one that one should remove... 
          if( hit_nr_worst >= 0 && hit_nr_worst <= int(buffer.size())  ) {
            // fill segment in buffer, delete worst hit
            buffer.erase(buffer.begin()+(hit_nr_worst),buffer.begin()+(hit_nr_worst+1));
            reduced_segments.push_back(buffer);
          }
          else {
            // only fill segment in array, do not delete anything
            reduced_segments.push_back(buffer);
          }
        }
      }
      
      // Loop over the subsegments and fit (only one segment if "use_brute_force" is false):
      for (uint iSegment=0; iSegment<reduced_segments.size(); iSegment++ ) {
        // loop over hits on given segment and push pointers to hits into protosegment
        protoSegment.clear();
        for (uint m = 0; m<reduced_segments[iSegment].size(); ++m ) {
          protoSegment.push_back(&reduced_segments[iSegment][m]);
        }
        fitSlopes(); 
        fillChiSquared();
        fillLocalDirection();
        // calculate 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);

        // replace n hit segment with n-1 hit segment, if segment probability is 1e3 better:
        if( ( ChiSquaredProbability((double)(*it).chi2(),(double)((2*(*it).nRecHits())-4)) 
              < 
              (1.e-3)*(ChiSquaredProbability((double)temp.chi2(),(double)(2*temp.nRecHits()-4))) )
            && 
            ( (ChiSquaredProbability((double)temp.chi2(),(double)(2*temp.nRecHits()-4))) 
              > best_red_seg_prob 
              )
            &&
            ( (ChiSquaredProbability((double)temp.chi2(),(double)(2*temp.nRecHits()-4))) > 1e-10 )
            ) {
          best_red_seg_prob = ChiSquaredProbability((double)temp.chi2(),(double)(2*temp.nRecHits()-4));
          // exchange current n hit segment (*it) with better n-1 hit segment:
          (*it) = temp;
        }
      }
  }
  
  return segments;
  
}

AlgebraicSymMatrix CSCSegAlgoHitPruning::weightMatrix

(

void

)

const [private]

Definition at line 367 of file CSCSegAlgoHitPruning.cc.

References it, CSCRecHit2D::localPositionError(), matrix, protoSegment, row, 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;
}

---

Member Data Documentation

bool CSCSegAlgoHitPruning::BrutePruning [private]

Definition at line 56 of file CSCSegAlgoHitPruning.h.

Referenced by CSCSegAlgoHitPruning(), and pruneBadHits().

double CSCSegAlgoHitPruning::protoChi2 [private]

Definition at line 53 of file CSCSegAlgoHitPruning.h.

Referenced by fillChiSquared(), and pruneBadHits().

LocalVector CSCSegAlgoHitPruning::protoDirection [private]

Definition at line 54 of file CSCSegAlgoHitPruning.h.

Referenced by fillLocalDirection(), and pruneBadHits().

LocalPoint CSCSegAlgoHitPruning::protoIntercept [private]

Definition at line 52 of file CSCSegAlgoHitPruning.h.

Referenced by fillChiSquared(), fillLocalDirection(), fitSlopes(), and pruneBadHits().

ChamberHitContainer CSCSegAlgoHitPruning::protoSegment [private]

Definition at line 49 of file CSCSegAlgoHitPruning.h.

Referenced by derivativeMatrix(), fillChiSquared(), fitSlopes(), pruneBadHits(), and weightMatrix().

float CSCSegAlgoHitPruning::protoSlope_u [private]

Definition at line 50 of file CSCSegAlgoHitPruning.h.

Referenced by fillChiSquared(), fillLocalDirection(), and fitSlopes().

float CSCSegAlgoHitPruning::protoSlope_v [private]

Definition at line 51 of file CSCSegAlgoHitPruning.h.

Referenced by fillChiSquared(), fillLocalDirection(), and fitSlopes().

const CSCChamber* CSCSegAlgoHitPruning::theChamber [private]

Definition at line 47 of file CSCSegAlgoHitPruning.h.

Referenced by derivativeMatrix(), fillChiSquared(), fillLocalDirection(), fitSlopes(), and pruneBadHits().

---

The documentation for this class was generated from the following files:

• RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.h
• RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.cc

---