#include <CSCSegAlgoHitPruning.h>

Public Types
typedef std::vector< const CSCRecHit2D * >	ChamberHitContainer

Public Member Functions
	CSCSegAlgoHitPruning (const edm::ParameterSet &ps)
	constructor More...

std::vector< CSCSegment >	pruneBadHits (const CSCChamber *aChamber, const std::vector< CSCSegment > &segments)
	clusterize More...

	~CSCSegAlgoHitPruning ()
	destructor More...

Private Member Functions
AlgebraicSymMatrix	calculateError (void) const

CLHEP::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::getParameter().

                                                                     {
   BrutePruning           = ps.getParameter<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(), query::result, weightMatrix(), and create_public_pileup_plots::weights.

Referenced by pruneBadHits().

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

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

private

Definition at line 393 of file CSCSegAlgoHitPruning.cc.

References CSCChamber::layer(), 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 CSCSegAlgoHitPruning::fillChiSquared ( void )

private

Definition at line 298 of file CSCSegAlgoHitPruning.cc.

References CSCChamber::layer(), LogDebug, protoChi2, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), findQualityFiles::v, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::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;
     
     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";
 //       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, mathSSE::sqrt(), theChamber, GeomDet::toGlobal(), csvLumiCalc::unit, and detailsBasic3DVector::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 CSCChamber::layer(), LogDebug, AlCaHLTBitMon_ParallelJobs::p, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), findQualityFiles::v, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by pruneBadHits().

                                      {
   CLHEP::HepMatrix M(4,4,0);
   CLHEP::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 
     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";      
 //       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;
   }
   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);
 }

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,
		const std::vector< CSCSegment > &	segments
	)

clusterize

Definition at line 49 of file CSCSegAlgoHitPruning.cc.

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

                                                                                                                            {
 
   theChamber = aChamber;
 
   std::vector<CSCSegment>          segments_temp;
   std::vector<ChamberHitContainer> rechits_clusters; 
   std::vector<CSCSegment> segments = _segments;
   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(size_t 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 (size_t iSegment=0; iSegment<reduced_segments.size(); iSegment++ ) {
         // loop over hits on given segment and push pointers to hits into protosegment
         protoSegment.clear();
         for (size_t 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 makeMuonMisalignmentScenario::matrix, and protoSegment.

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;
 }