From c764d245cccd9490c58b4bac5b457f6719e2ac43 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Wed, 18 Feb 2015 18:44:45 +0100 Subject: [PATCH 01/17] remove unused and unwanted file --- .../CSCSegment/src/CSCSegAlgoHitPruning.h | 58 ------------------- 1 file changed, 58 deletions(-) delete mode 100644 RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.h diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.h b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.h deleted file mode 100644 index 8bfad5b9c3c22..0000000000000 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.h +++ /dev/null @@ -1,58 +0,0 @@ -#ifndef CSCSegment_CSCSegAlgoHitPruning_h -#define CSCSegAlgoHitPruning_h -/** - * \file CSCSegAlgoHitPruning.h - * - * \authors: S. Stoynev - NU - * I. Bloch - FNAL - * E. James - FNAL - * D. Fortin - UC Riverside - * - * Go through segments and clean up selection of hits according - * to changes in quality of segment (chi2/d.o.f. probability) - */ - -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include -#include - -#include - -class CSCChamber; - -class CSCSegAlgoHitPruning { - - public: - - typedef std::vector ChamberHitContainer; - - /// constructor - explicit CSCSegAlgoHitPruning(const edm::ParameterSet& ps); - - /// destructor - ~CSCSegAlgoHitPruning(); - - /// clusterize - std::vector pruneBadHits(const CSCChamber* aChamber, const std::vector& segments); - - private: - void fitSlopes(void); - void fillChiSquared(void); - void fillLocalDirection(void); - CLHEP::HepMatrix derivativeMatrix(void) const; - AlgebraicSymMatrix weightMatrix(void) const; - AlgebraicSymMatrix calculateError(void) const; - void flipErrors(AlgebraicSymMatrix&) const; - - const CSCChamber* theChamber; - - ChamberHitContainer protoSegment; - float protoSlope_u; - float protoSlope_v; - LocalPoint protoIntercept; - double protoChi2; - LocalVector protoDirection; - - bool BrutePruning; -}; -#endif From 28ccdd24292e78d51b15d5492ce50eeb5d804b77 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Wed, 18 Feb 2015 18:44:52 +0100 Subject: [PATCH 02/17] remove unused and unwanted file --- .../CSCSegment/src/CSCSegAlgoHitPruning.cc | 458 ------------------ 1 file changed, 458 deletions(-) delete mode 100644 RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.cc diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.cc deleted file mode 100644 index bda2ee00042a6..0000000000000 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.cc +++ /dev/null @@ -1,458 +0,0 @@ -/** - * \file CSCSegAlgoHitPruning.cc - * - * \authors: S. Stoynev - NU - * I. Bloch - FNAL - * E. James - FNAL - * D. Fortin - UC Riverside - * - * See header file for description. - */ - -#include "RecoLocalMuon/CSCSegment/src/CSCSegAlgoHitPruning.h" - -#include "Geometry/CSCGeometry/interface/CSCLayer.h" -#include -#include -#include -#include "CommonTools/Statistics/interface/ChiSquaredProbability.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/MessageLogger/interface/MessageLogger.h" - -#include -#include -#include -#include - - -/* Constructor - * - */ -CSCSegAlgoHitPruning::CSCSegAlgoHitPruning(const edm::ParameterSet& ps) { - BrutePruning = ps.getParameter("BrutePruning"); - -} - - -/* Destructor: - * - */ -CSCSegAlgoHitPruning::~CSCSegAlgoHitPruning(){ - -} - - -/* pruneBadHits - * - */ -std::vector CSCSegAlgoHitPruning::pruneBadHits(const CSCChamber* aChamber, const std::vector& _segments) { - - theChamber = aChamber; - - std::vector segments_temp; - std::vector rechits_clusters; - std::vector 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::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 theseRecHits = (*it).specificRecHits(); - std::vector::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::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 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; - -} - - -/* Method fitSlopes - * - * Perform a Least Square Fit on a segment as per SK algo - * - */ -void CSCSegAlgoHitPruning::fitSlopes() { - 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"<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; -} - - -/* fillLocalDirection - * - */ -void CSCSegAlgoHitPruning::fillLocalDirection() { - // 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(); -} - - -/* weightMatrix - * - */ -AlgebraicSymMatrix CSCSegAlgoHitPruning::weightMatrix() const { - - std::vector::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; -} - - -/* derivativeMatrix - * - */ -CLHEP::HepMatrix CSCSegAlgoHitPruning::derivativeMatrix() const { - - 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; -} - - -/* calculateError - * - */ -AlgebraicSymMatrix CSCSegAlgoHitPruning::calculateError() const { - - 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 CSCSegAlgoHitPruning::flipErrors( AlgebraicSymMatrix& a ) const { - - // 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 - -} - From 2cb7ac1424c5f33288d4fca2e2a477a29d1750eb Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:09:46 +0100 Subject: [PATCH 03/17] remove CLHEP --- RecoLocalMuon/CSCSegment/BuildFile.xml | 1 - 1 file changed, 1 deletion(-) diff --git a/RecoLocalMuon/CSCSegment/BuildFile.xml b/RecoLocalMuon/CSCSegment/BuildFile.xml index 1a9e09d49a8f4..6cb3a8b73d79b 100644 --- a/RecoLocalMuon/CSCSegment/BuildFile.xml +++ b/RecoLocalMuon/CSCSegment/BuildFile.xml @@ -1,5 +1,4 @@ - From d06a52a99c222517f3b9e0a22df1d05855fe0a66 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:10:45 +0100 Subject: [PATCH 04/17] include me42 --- .../python/CSCSegmentAlgorithmDF_cfi.py | 39 +++++++------------ .../python/CSCSegmentAlgorithmSK_cfi.py | 22 ++--------- .../python/CSCSegmentAlgorithmTC_cfi.py | 20 ++-------- 3 files changed, 20 insertions(+), 61 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmDF_cfi.py b/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmDF_cfi.py index 3e33019d78b74..994aa1d0db424 100644 --- a/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmDF_cfi.py +++ b/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmDF_cfi.py @@ -1,8 +1,8 @@ import FWCore.ParameterSet.Config as cms -# The following algorithms looks how far a rechit is from the +# The following DF algorithms looks how far a rechit is from the # proto segment in terms of its error ellipse. This is different -# from the other algorithms which use a cylinder around the proto +# from SK & TC algorithms which use a cylinder around the proto # segment and look for rechits within that cylinder DF_ME1234_1 = cms.PSet( preClustering = cms.untracked.bool(False), @@ -19,7 +19,9 @@ nHitsPerClusterIsShower = cms.int32(20), minLayersApart = cms.int32(2), Pruning = cms.untracked.bool(False), - dYclusBoxMax = cms.double(8.0) + dYclusBoxMax = cms.double(8.0), + maxDPhi = cms.double(999.), + maxDTheta = cms.double(999.) ) DF_ME1234_2 = cms.PSet( preClustering = cms.untracked.bool(False), @@ -36,7 +38,9 @@ nHitsPerClusterIsShower = cms.int32(20), minLayersApart = cms.int32(2), Pruning = cms.untracked.bool(False), - dYclusBoxMax = cms.double(12.0) + dYclusBoxMax = cms.double(12.0), + maxDPhi = cms.double(999.), + maxDTheta = cms.double(999.) ) DF_ME1A = cms.PSet( preClustering = cms.untracked.bool(False), @@ -53,29 +57,14 @@ nHitsPerClusterIsShower = cms.int32(20), minLayersApart = cms.int32(2), Pruning = cms.untracked.bool(False), - dYclusBoxMax = cms.double(8.0) + dYclusBoxMax = cms.double(8.0), + maxDPhi = cms.double(999.), + maxDTheta = cms.double(999.) ) CSCSegAlgoDF = cms.PSet( - chamber_types = cms.vstring('ME1/a', - 'ME1/b', - 'ME1/2', - 'ME1/3', - 'ME2/1', - 'ME2/2', - 'ME3/1', - 'ME3/2', - 'ME4/1'), + chamber_types = cms.vstring('ME1/a', 'ME1/b', 'ME1/2', 'ME1/3', 'ME2/1', 'ME2/2', 'ME3/1', 'ME3/2', 'ME4/1','ME4/2'), algo_name = cms.string('CSCSegAlgoDF'), - algo_psets = cms.VPSet(cms.PSet( - DF_ME1234_1 - ), - cms.PSet( - DF_ME1234_2 - ), - cms.PSet( - DF_ME1A - )), - parameters_per_chamber_type = cms.vint32(3, 1, 2, 2, 1, - 2, 1, 2, 1) + algo_psets = cms.VPSet( cms.PSet(DF_ME1234_1), cms.PSet(DF_ME1234_2), cms.PSet(DF_ME1A) ), + parameters_per_chamber_type = cms.vint32( 3, 1, 2, 2, 1, 2, 1, 2, 1, 2 ) ) diff --git a/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmSK_cfi.py b/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmSK_cfi.py index 4597b4552f216..d7d793badbcd0 100644 --- a/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmSK_cfi.py +++ b/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmSK_cfi.py @@ -23,25 +23,9 @@ dRPhiMax = cms.double(8.0) ) CSCSegAlgoSK = cms.PSet( - chamber_types = cms.vstring('ME1/a', - 'ME1/b', - 'ME1/2', - 'ME1/3', - 'ME2/1', - 'ME2/2', - 'ME3/1', - 'ME3/2', - 'ME4/1'), - # vstring chamber_types = {"ME1/a", "ME1/b", "ME1/2", "ME1/3", "ME2/1", "ME2/2", "ME3/1", "ME3/2", "ME4/1", "ME4/2"} - # vint32 parameters_per_chamber_type = {2, 1, 1, 1, 1, 1, 1, 1, 1, 1} + chamber_types = cms.vstring('ME1/a', 'ME1/b', 'ME1/2', 'ME1/3', 'ME2/1', 'ME2/2', 'ME3/1', 'ME3/2', 'ME4/1','ME4/2'), algo_name = cms.string('CSCSegAlgoSK'), - algo_psets = cms.VPSet(cms.PSet( - SK_ME1234 - ), - cms.PSet( - SK_ME1A - )), - parameters_per_chamber_type = cms.vint32(2, 1, 1, 1, 1, - 1, 1, 1, 1) + algo_psets = cms.VPSet( cms.PSet(SK_ME1234), cms.PSet(SK_ME1A) ), + parameters_per_chamber_type = cms.vint32(2, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) ) diff --git a/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmTC_cfi.py b/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmTC_cfi.py index 48fd0efd0a0e6..864c46b860ed0 100644 --- a/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmTC_cfi.py +++ b/RecoLocalMuon/CSCSegment/python/CSCSegmentAlgorithmTC_cfi.py @@ -28,23 +28,9 @@ dRPhiMax = cms.double(0.6) ) CSCSegAlgoTC = cms.PSet( - chamber_types = cms.vstring('ME1/a', - 'ME1/b', - 'ME1/2', - 'ME1/3', - 'ME2/1', - 'ME2/2', - 'ME3/1', - 'ME3/2', - 'ME4/1'), + chamber_types = cms.vstring('ME1/a', 'ME1/b', 'ME1/2', 'ME1/3', 'ME2/1', 'ME2/2', 'ME3/1', 'ME3/2', 'ME4/1','ME4/2'), algo_name = cms.string('CSCSegAlgoTC'), - algo_psets = cms.VPSet(cms.PSet( - TC_ME1234 - ), - cms.PSet( - TC_ME1A - )), - parameters_per_chamber_type = cms.vint32(2, 1, 1, 1, 1, - 1, 1, 1, 1) + algo_psets = cms.VPSet( cms.PSet(TC_ME1234), cms.PSet(TC_ME1A) ), + parameters_per_chamber_type = cms.vint32( 2, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) ) From 8e503e25a074090c641ab09e3a11291d4e10c546 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:12:20 +0100 Subject: [PATCH 05/17] fix broken formal structure --- RecoLocalMuon/CSCSegment/python/cscSegments_cfi.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/python/cscSegments_cfi.py b/RecoLocalMuon/CSCSegment/python/cscSegments_cfi.py index e3220b81a00e6..f2a28df12c0e1 100644 --- a/RecoLocalMuon/CSCSegment/python/cscSegments_cfi.py +++ b/RecoLocalMuon/CSCSegment/python/cscSegments_cfi.py @@ -4,17 +4,17 @@ from RecoLocalMuon.CSCSegment.CSCSegmentAlgorithmTC_cfi import * from RecoLocalMuon.CSCSegment.CSCSegmentAlgorithmDF_cfi import * from RecoLocalMuon.CSCSegment.CSCSegmentAlgorithmST_cfi import * -# module must be an EDProducer or similar -# ======================================= + cscSegments = cms.EDProducer("CSCSegmentProducer", # Define input inputObjects = cms.InputTag("csc2DRecHits"), # Choice of the building algo: 1 SK, 2 TC, 3 DF, 4 ST, ... algo_type = cms.int32(4), # std::vector - algo_psets = cms.VPSet(cms.PSet( - CSCSegAlgoSK - ), + algo_psets = cms.VPSet( + cms.PSet( + CSCSegAlgoSK + ), cms.PSet( CSCSegAlgoTC ), @@ -23,7 +23,8 @@ ), cms.PSet( CSCSegAlgoST - )) + ) + ) ) From d6f5719db6f30fab8a361b989212c87fbcf56c8a Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:14:25 +0100 Subject: [PATCH 06/17] update to run in 74x --- .../CSCSegment/test/run_on_raw_74x.py | 82 ++++++++++++++++ .../CSCSegment/test/run_on_simdigi.py | 2 +- .../CSCSegment/test/run_on_simdigi_74x.py | 95 ++++++++++++++++++ .../CSCSegment/test/run_on_simdigi_df.py | 88 +++++++++++++++++ .../CSCSegment/test/run_on_simdigi_sk.py | 84 ++++++++++++++++ .../CSCSegment/test/run_on_simdigi_st.py | 97 +++++++++++++++++++ .../CSCSegment/test/run_on_simdigi_tc.py | 81 ++++++++++++++++ 7 files changed, 528 insertions(+), 1 deletion(-) create mode 100644 RecoLocalMuon/CSCSegment/test/run_on_raw_74x.py create mode 100644 RecoLocalMuon/CSCSegment/test/run_on_simdigi_74x.py create mode 100644 RecoLocalMuon/CSCSegment/test/run_on_simdigi_df.py create mode 100644 RecoLocalMuon/CSCSegment/test/run_on_simdigi_sk.py create mode 100644 RecoLocalMuon/CSCSegment/test/run_on_simdigi_st.py create mode 100644 RecoLocalMuon/CSCSegment/test/run_on_simdigi_tc.py diff --git a/RecoLocalMuon/CSCSegment/test/run_on_raw_74x.py b/RecoLocalMuon/CSCSegment/test/run_on_raw_74x.py new file mode 100644 index 0000000000000..c0f26ca2e69e8 --- /dev/null +++ b/RecoLocalMuon/CSCSegment/test/run_on_raw_74x.py @@ -0,0 +1,82 @@ +## Process 100 events in CSC segment builder - Tim Cox - 22.01.2015 +## This version runs in 74X on a Real RAW relval sample in 730 +## Now testing on FullSim+PU TTbar sample in 730 + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("TEST") + +process.load("Configuration/StandardSequences/Geometry_cff") +process.load("Configuration/StandardSequences/MagneticField_cff") +process.load("Configuration/StandardSequences/FrontierConditions_GlobalTag_cff") +process.load("Configuration/StandardSequences/RawToDigi_Data_cff") +process.load("Configuration.StandardSequences.Reconstruction_cff") +process.load("Configuration.StandardSequences.EndOfProcess_cff") + +# --- MATCH GT TO RELEASE AND DATA SAMPLE + +# Tag for 730 Real Data relvals +##process.GlobalTag.globaltag = "GR_P_V43D::All" + +# Tag for 730 Sim relvals +process.GlobalTag.globaltag = "MCRUN2_73_V5::All" + + +# --- NUMBER OF EVENTS --- + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(100) ) + +process.options = cms.untracked.PSet( SkipEvent = cms.untracked.vstring("ProductNotFound") ) +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( +## "/store/relval/CMSSW_7_3_0/SingleMu/RAW/GR_H_V43A_RelVal_zMu2012D-v1/00000/1ED9BE30-8481-E411-8AE9-002618943874.root" + "/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/044157C8-A181-E411-AC04-002354EF3BD2.root" + ) +) + +# -- ACCESSING 'DEEP' PARAMETERS OF THE ALGO IS TRICKY +# THE FOLLOWING FOUND BY EXPLORING CONFIG WITH python -i +# '3' is 4th algo CSCSegAlgoST; '0' and '1' are for ST_ME1234 and ST_ME1A configs +process.cscSegments.algo_psets[3].algo_psets[0].CSCDebug = cms.untracked.bool(True) +process.cscSegments.algo_psets[3].algo_psets[1].CSCDebug = cms.untracked.bool(True) + + +# --- Activate LogVerbatim IN CSCSegment +process.MessageLogger.categories.append("CSCSegment") +process.MessageLogger.destinations = cms.untracked.vstring("cout") +process.MessageLogger.cout = cms.untracked.PSet( + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet( limit = cms.untracked.int32(0) ), + FwkReport = cms.untracked.PSet( limit = cms.untracked.int32(-1) ), + CSCSegment = cms.untracked.PSet( limit = cms.untracked.int32(-1) ) +) + +### --- ACTIVATE LogTrace IN VARIOUS MODULES - NEED TO COMPILE *EACH MODULE* WITH +### scram b -j8 USER_CXXFLAGS="-DEDM_ML_DEBUG" +### LogTrace output goes to cout; all other output to "junk.log" + +#process.load("FWCore.MessageLogger.MessageLogger_cfi") +#process.MessageLogger.categories.append("CSCRecHit") +#process.MessageLogger.categories.append("CSCSegment") +#process.MessageLogger.categories.append("CSCSegAlgoST") + +### module label is something like "muonCSCDigis"... +#process.MessageLogger.debugModules = cms.untracked.vstring("*") +#process.MessageLogger.destinations = cms.untracked.vstring("cout","junk") +#process.MessageLogger.cout = cms.untracked.PSet( +# threshold = cms.untracked.string("DEBUG"), +# default = cms.untracked.PSet( limit = cms.untracked.int32(0) ), +# FwkReport = cms.untracked.PSet( limit = cms.untracked.int32(-1) ), +# CSCRecHit = cms.untracked.PSet( limit = cms.untracked.int32(-1) ), +# CSCSegment = cms.untracked.PSet( limit = cms.untracked.int32(-1) ) +# , CSCSegAlgoST = cms.untracked.PSet( limit = cms.untracked.int32(-1) ) +#) + +# Path and EndPath def +process.unpack = cms.Path(process.muonCSCDigis) +process.reco = cms.Path(process.csc2DRecHits * process.cscSegments) +process.endjob = cms.EndPath(process.endOfProcess) + +# Schedule definition +process.schedule = cms.Schedule(process.unpack, process.reco, process.endjob) diff --git a/RecoLocalMuon/CSCSegment/test/run_on_simdigi.py b/RecoLocalMuon/CSCSegment/test/run_on_simdigi.py index 8a8c425b6a0a9..ffb6724a1017d 100644 --- a/RecoLocalMuon/CSCSegment/test/run_on_simdigi.py +++ b/RecoLocalMuon/CSCSegment/test/run_on_simdigi.py @@ -38,7 +38,7 @@ process.csc2DRecHits.readBadChannels = cms.bool(False) process.csc2DRecHits.CSCUseTimingCorrections = cms.bool(False) -process.csc2DRecHits.CSCUseGasGainCorrection = cms.bool(False) +process.csc2DRecHits.CSCUseGasGainCorrections = cms.bool(False) # Switch input for CSCRecHitD to s i m u l a t e d digis diff --git a/RecoLocalMuon/CSCSegment/test/run_on_simdigi_74x.py b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_74x.py new file mode 100644 index 0000000000000..0baaa96418477 --- /dev/null +++ b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_74x.py @@ -0,0 +1,95 @@ +## Process sim digi events with CSC rechit & segment builders - Tim Cox - 11.02.2015 +## This version runs in 7_4_0_preX on a 7_3_0 simulated data DIGI relval sample. +## Run on 1000 events of a 25ns PU TTbar sample + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("TEST") + +# Geometry access changed after Yana hypernews post 10.02.2015 +# had been using... +##process.load("Configuration.StandardSequences.Geometry_cff") +# which just points to... +##process.load("Configuration.Geometry.GeometryIdeal_cff") +# yana wants... +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") + +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +##process.load("Configuration.StandardSequences.RawToDigi_Data_cff") +process.load("Configuration.StandardSequences.Reconstruction_cff") +process.load("Configuration.StandardSequences.EndOfProcess_cff") + +process.load("CalibMuon.CSCCalibration.CSCChannelMapper_cfi") +process.load("CalibMuon.CSCCalibration.CSCIndexer_cfi") +process.CSCIndexerESProducer.AlgoName = cms.string("CSCIndexerPostls1") +process.CSCChannelMapperESProducer.AlgoName = cms.string("CSCChannelMapperPostls1") + +# --- MATCH GT TO RELEASE AND DATA SAMPLE + +process.GlobalTag.globaltag = "MCRUN2_73_V5::All" + +# --- NUMBER OF EVENTS + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000) ) + +process.options = cms.untracked.PSet( SkipEvent = cms.untracked.vstring("ProductNotFound") ) +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) + +## ttbar+pu is 200 events per file so need 5 for 1000 events +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + "/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/044157C8-A181-E411-AC04-002354EF3BD2.root" + , +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/0A963931-A181-E411-B4C5-0026189438DC.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/145BE1DC-A181-E411-816D-0025905A609E.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/14DDEDD0-A181-E411-9476-0026189438F8.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/18F85F36-A181-E411-8AF4-0025905B85D6.root" + ) +) + +# ME1/1A is u n g a n g e d postls1 + +process.CSCGeometryESModule.useGangedStripsInME1a = False +##process.CSCGeometryESModule.debugV = True +##process.idealForDigiCSCGeometry.useGangedStripsInME1a = False + +# Turn off some flags for CSCRecHitD that are turned ON in default config + +process.csc2DRecHits.readBadChannels = cms.bool(False) +process.csc2DRecHits.CSCUseGasGainCorrections = cms.bool(False) +# Already defaults OFF... +## process.csc2DRecHits.CSCUseTimingCorrections = cms.bool(False) + +# Switch input for CSCRecHitD to s i m u l a t e d digis + +process.csc2DRecHits.wireDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCWireDigi") +process.csc2DRecHits.stripDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCStripDigi") + +# -- ACCESSING "DEEP" PARAMETERS OF THE ALGO IS TRICKY +# THE FOLLOWING FOUND BY EXPLORING CONFIG WITH python -i +# "3" is 4th algo CSCSegAlgoST; "0" and "1" are for ST_ME1234 and ST_ME1A configs +# activate dump of segments +process.cscSegments.algo_psets[3].algo_psets[0].CSCDebug = cms.untracked.bool(True) +process.cscSegments.algo_psets[3].algo_psets[1].CSCDebug = cms.untracked.bool(True) + +# activate the special shower code +process.cscSegments.algo_psets[3].algo_psets[0].useShowering = cms.bool(True) +process.cscSegments.algo_psets[3].algo_psets[1].useShowering = cms.bool(True) + +# --- Activate LogVerbatim IN CSCSegment +process.MessageLogger.categories.append("CSCSegment") +process.MessageLogger.destinations = cms.untracked.vstring("cout") +process.MessageLogger.cout = cms.untracked.PSet( + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet( limit = cms.untracked.int32(0) ), + FwkReport = cms.untracked.PSet( limit = cms.untracked.int32(-1) ), + CSCSegment = cms.untracked.PSet( limit = cms.untracked.int32(-1) ) +) + +# Path and EndPath def +process.reco = cms.Path(process.csc2DRecHits * process.cscSegments) +process.endjob = cms.EndPath(process.endOfProcess) + +# Schedule definition +process.schedule = cms.Schedule(process.reco, process.endjob) diff --git a/RecoLocalMuon/CSCSegment/test/run_on_simdigi_df.py b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_df.py new file mode 100644 index 0000000000000..7b798767cac17 --- /dev/null +++ b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_df.py @@ -0,0 +1,88 @@ +## Process sim digi events with CSC rechit & segment builders - Tim Cox - 16.02.2015 +## This version runs in 7_4_0_preX on a 7_3_0 simulated data DIGI relval sample. +## -- USING OLD ALGO 'DF' +## Run on 100 events of a 25ns PU TTbar sample + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("TEST") + +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +##process.load("Configuration.StandardSequences.RawToDigi_Data_cff") +process.load("Configuration.StandardSequences.Reconstruction_cff") +process.load("Configuration.StandardSequences.EndOfProcess_cff") + +process.load("CalibMuon.CSCCalibration.CSCChannelMapper_cfi") +process.load("CalibMuon.CSCCalibration.CSCIndexer_cfi") +process.CSCIndexerESProducer.AlgoName = cms.string("CSCIndexerPostls1") +process.CSCChannelMapperESProducer.AlgoName = cms.string("CSCChannelMapperPostls1") + +# --- MATCH GT TO RELEASE AND DATA SAMPLE + +process.GlobalTag.globaltag = "MCRUN2_73_V5::All" + +# --- NUMBER OF EVENTS + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(100) ) + +process.options = cms.untracked.PSet( SkipEvent = cms.untracked.vstring("ProductNotFound") ) +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) + +## ttbar+pu is 200 events per file so need 5 for 1000 events +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + "/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/044157C8-A181-E411-AC04-002354EF3BD2.root" + , +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/0A963931-A181-E411-B4C5-0026189438DC.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/145BE1DC-A181-E411-816D-0025905A609E.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/14DDEDD0-A181-E411-9476-0026189438F8.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/18F85F36-A181-E411-8AF4-0025905B85D6.root" + ) +) + +# ME1/1A is u n g a n g e d postls1 + +process.CSCGeometryESModule.useGangedStripsInME1a = False +##process.CSCGeometryESModule.debugV = True +##process.idealForDigiCSCGeometry.useGangedStripsInME1a = False + +# Turn off some flags for CSCRecHitD that are turned ON in default config + +process.csc2DRecHits.readBadChannels = cms.bool(False) +process.csc2DRecHits.CSCUseGasGainCorrections = cms.bool(False) +# Already defaults OFF... +## process.csc2DRecHits.CSCUseTimingCorrections = cms.bool(False) + +# Switch input for CSCRecHitD to s i m u l a t e d digis + +process.csc2DRecHits.wireDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCWireDigi") +process.csc2DRecHits.stripDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCStripDigi") + +# switch to CSCSegAlgoDF +process.cscSegments.algo_type = cms.int32( 3 ) + +# Turn on segment dumps +# "2" is 3rd algo CSCSegAlgoDF; "0","1","2" are for DF_ME1234_1, DF_ME1234_2, DF_ME1A configs +process.cscSegments.algo_psets[2].algo_psets[0].CSCDegmentDebug = cms.untracked.bool(True) +process.cscSegments.algo_psets[2].algo_psets[1].CSCSegmentDebug = cms.untracked.bool(True) +process.cscSegments.algo_psets[2].algo_psets[2].CSCDegmentDebug = cms.untracked.bool(True) + +# --- Activate LogVerbatim IN CSCSegment + +process.MessageLogger.categories.append("CSCSegment") +process.MessageLogger.destinations = cms.untracked.vstring("cout") +process.MessageLogger.cout = cms.untracked.PSet( + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet( limit = cms.untracked.int32(0) ), + FwkReport = cms.untracked.PSet( limit = cms.untracked.int32(-1) ), + CSCSegment = cms.untracked.PSet( limit = cms.untracked.int32(-1) ) +) + +# Path and EndPath def +process.reco = cms.Path(process.csc2DRecHits * process.cscSegments) +process.endjob = cms.EndPath(process.endOfProcess) + +# Schedule definition +process.schedule = cms.Schedule(process.reco, process.endjob) diff --git a/RecoLocalMuon/CSCSegment/test/run_on_simdigi_sk.py b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_sk.py new file mode 100644 index 0000000000000..20c8f853ae809 --- /dev/null +++ b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_sk.py @@ -0,0 +1,84 @@ +## Process sim digi events with CSC rechit & segment builders - Tim Cox - 11.02.2015 +## This version runs in 7_4_0_preX on a 7_3_0 simulated data DIGI relval sample. +## -- USING OLD ALGO 'SK' -- +## Run on 100 events of a 25ns PU TTbar sample + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("TEST") + +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +##process.load("Configuration.StandardSequences.RawToDigi_Data_cff") +process.load("Configuration.StandardSequences.Reconstruction_cff") +process.load("Configuration.StandardSequences.EndOfProcess_cff") + +process.load("CalibMuon.CSCCalibration.CSCChannelMapper_cfi") +process.load("CalibMuon.CSCCalibration.CSCIndexer_cfi") +process.CSCIndexerESProducer.AlgoName = cms.string("CSCIndexerPostls1") +process.CSCChannelMapperESProducer.AlgoName = cms.string("CSCChannelMapperPostls1") + +# --- MATCH GT TO RELEASE AND DATA SAMPLE + +process.GlobalTag.globaltag = "MCRUN2_73_V5::All" + +# --- NUMBER OF EVENTS + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(100) ) + +process.options = cms.untracked.PSet( SkipEvent = cms.untracked.vstring("ProductNotFound") ) +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) + +## ttbar+pu is 200 events per file so need 5 for 1000 events +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + "/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/044157C8-A181-E411-AC04-002354EF3BD2.root" + , +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/0A963931-A181-E411-B4C5-0026189438DC.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/145BE1DC-A181-E411-816D-0025905A609E.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/14DDEDD0-A181-E411-9476-0026189438F8.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/18F85F36-A181-E411-8AF4-0025905B85D6.root" + ) +) + +# ME1/1A is u n g a n g e d postls1 + +process.CSCGeometryESModule.useGangedStripsInME1a = False +##process.CSCGeometryESModule.debugV = True +##process.idealForDigiCSCGeometry.useGangedStripsInME1a = False + +# Turn off some flags for CSCRecHitD that are turned ON in default config + +process.csc2DRecHits.readBadChannels = cms.bool(False) +process.csc2DRecHits.CSCUseGasGainCorrections = cms.bool(False) +# Already defaults OFF... +## process.csc2DRecHits.CSCUseTimingCorrections = cms.bool(False) + +# Switch input for CSCRecHitD to s i m u l a t e d digis + +process.csc2DRecHits.wireDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCWireDigi") +process.csc2DRecHits.stripDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCStripDigi") + +# -- ACCESSING "DEEP" PARAMETERS OF THE ALGO IS TRICKY +# THE FOLLOWING FOUND BY EXPLORING CONFIG WITH python -i + +# switch to CSCSegAlgoSK +process.cscSegments.algo_type = cms.int32(1) + +# --- Activate LogVerbatim IN CSCSegment +process.MessageLogger.categories.append("CSCSegment") +process.MessageLogger.destinations = cms.untracked.vstring("cout") +process.MessageLogger.cout = cms.untracked.PSet( + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet( limit = cms.untracked.int32(0) ), + FwkReport = cms.untracked.PSet( limit = cms.untracked.int32(-1) ), + CSCSegment = cms.untracked.PSet( limit = cms.untracked.int32(-1) ) +) + +# Path and EndPath def +process.reco = cms.Path(process.csc2DRecHits * process.cscSegments) +process.endjob = cms.EndPath(process.endOfProcess) + +# Schedule definition +process.schedule = cms.Schedule(process.reco, process.endjob) diff --git a/RecoLocalMuon/CSCSegment/test/run_on_simdigi_st.py b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_st.py new file mode 100644 index 0000000000000..944799a916908 --- /dev/null +++ b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_st.py @@ -0,0 +1,97 @@ +## Process sim digi events with CSC rechit & segment builders - Tim Cox - 11.02.2015 +## This version runs in 7_4_0_preX on a 7_3_0 simulated data DIGI relval sample. +## -- USING DEFAULT ALGO "ST" +## Run on 100 events of a 25ns PU TTbar sample + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("TEST") + +# Geometry access changed after Yana hypernews post 10.02.2015 +# had been using... +##process.load("Configuration.StandardSequences.Geometry_cff") +# which just points to... +##process.load("Configuration.Geometry.GeometryIdeal_cff") +# yana wants... +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") + +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +##process.load("Configuration.StandardSequences.RawToDigi_Data_cff") +process.load("Configuration.StandardSequences.Reconstruction_cff") +process.load("Configuration.StandardSequences.EndOfProcess_cff") + +process.load("CalibMuon.CSCCalibration.CSCChannelMapper_cfi") +process.load("CalibMuon.CSCCalibration.CSCIndexer_cfi") +process.CSCIndexerESProducer.AlgoName = cms.string("CSCIndexerPostls1") +process.CSCChannelMapperESProducer.AlgoName = cms.string("CSCChannelMapperPostls1") + +# --- MATCH GT TO RELEASE AND DATA SAMPLE + +process.GlobalTag.globaltag = "MCRUN2_73_V5::All" + +# --- NUMBER OF EVENTS + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(100) ) + +process.options = cms.untracked.PSet( SkipEvent = cms.untracked.vstring("ProductNotFound") ) +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) + +## ttbar+pu is 200 events per file so need 5 for 1000 events +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + "/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/044157C8-A181-E411-AC04-002354EF3BD2.root" + , +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/0A963931-A181-E411-B4C5-0026189438DC.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/145BE1DC-A181-E411-816D-0025905A609E.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/14DDEDD0-A181-E411-9476-0026189438F8.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/18F85F36-A181-E411-8AF4-0025905B85D6.root" + ) +) + +# ME1/1A is u n g a n g e d postls1 + +process.CSCGeometryESModule.useGangedStripsInME1a = False +##process.CSCGeometryESModule.debugV = True +##process.idealForDigiCSCGeometry.useGangedStripsInME1a = False + +# Turn off some flags for CSCRecHitD that are turned ON in default config + +process.csc2DRecHits.readBadChannels = cms.bool(False) +process.csc2DRecHits.CSCUseGasGainCorrections = cms.bool(False) +# Already defaults OFF... +## process.csc2DRecHits.CSCUseTimingCorrections = cms.bool(False) + +# Switch input for CSCRecHitD to s i m u l a t e d digis + +process.csc2DRecHits.wireDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCWireDigi") +process.csc2DRecHits.stripDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCStripDigi") + +# -- ACCESSING "DEEP" PARAMETERS OF THE ALGO IS TRICKY +# THE FOLLOWING FOUND BY EXPLORING CONFIG WITH python -i + +# activate dump of segments +# "3" is 4th algo CSCSegAlgoST; "0" and "1" are for ST_ME1234 and ST_ME1A configs +process.cscSegments.algo_psets[3].algo_psets[0].CSCDebug = cms.untracked.bool(True) +process.cscSegments.algo_psets[3].algo_psets[1].CSCDebug = cms.untracked.bool(True) + +# activate the special shower code +process.cscSegments.algo_psets[3].algo_psets[0].useShowering = cms.bool(True) +process.cscSegments.algo_psets[3].algo_psets[1].useShowering = cms.bool(True) + +# --- Activate LogVerbatim IN CSCSegment +process.MessageLogger.categories.append("CSCSegment") +process.MessageLogger.destinations = cms.untracked.vstring("cout") +process.MessageLogger.cout = cms.untracked.PSet( + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet( limit = cms.untracked.int32(0) ), + FwkReport = cms.untracked.PSet( limit = cms.untracked.int32(-1) ), + CSCSegment = cms.untracked.PSet( limit = cms.untracked.int32(-1) ) +) + +# Path and EndPath def +process.reco = cms.Path(process.csc2DRecHits * process.cscSegments) +process.endjob = cms.EndPath(process.endOfProcess) + +# Schedule definition +process.schedule = cms.Schedule(process.reco, process.endjob) diff --git a/RecoLocalMuon/CSCSegment/test/run_on_simdigi_tc.py b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_tc.py new file mode 100644 index 0000000000000..99b98adb4d12c --- /dev/null +++ b/RecoLocalMuon/CSCSegment/test/run_on_simdigi_tc.py @@ -0,0 +1,81 @@ +## Process sim digi events with CSC rechit & segment builders - Tim Cox - 12.02.2015 +## This version runs in 7_4_0_preX on a 7_3_0 simulated data DIGI relval sample. +## -- USING OLD ALGO 'TC' +## Run on 100 events of a 25ns PU TTbar sample + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("TEST") + +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +##process.load("Configuration.StandardSequences.RawToDigi_Data_cff") +process.load("Configuration.StandardSequences.Reconstruction_cff") +process.load("Configuration.StandardSequences.EndOfProcess_cff") + +process.load("CalibMuon.CSCCalibration.CSCChannelMapper_cfi") +process.load("CalibMuon.CSCCalibration.CSCIndexer_cfi") +process.CSCIndexerESProducer.AlgoName = cms.string("CSCIndexerPostls1") +process.CSCChannelMapperESProducer.AlgoName = cms.string("CSCChannelMapperPostls1") + +# --- MATCH GT TO RELEASE AND DATA SAMPLE + +process.GlobalTag.globaltag = "MCRUN2_73_V5::All" + +# --- NUMBER OF EVENTS + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(100) ) + +process.options = cms.untracked.PSet( SkipEvent = cms.untracked.vstring("ProductNotFound") ) +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) + +## ttbar+pu is 200 events per file so need 5 for 1000 events +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + "/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/044157C8-A181-E411-AC04-002354EF3BD2.root" + , +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/0A963931-A181-E411-B4C5-0026189438DC.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/145BE1DC-A181-E411-816D-0025905A609E.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/14DDEDD0-A181-E411-9476-0026189438F8.root", +"/store/relval/CMSSW_7_3_0/RelValTTbar_13/GEN-SIM-DIGI-RAW-HLTDEBUG/PU25ns_MCRUN2_73_V7_71XGENSIM-v1/00000/18F85F36-A181-E411-8AF4-0025905B85D6.root" + ) +) + +# ME1/1A is u n g a n g e d postls1 + +process.CSCGeometryESModule.useGangedStripsInME1a = False +##process.CSCGeometryESModule.debugV = True +##process.idealForDigiCSCGeometry.useGangedStripsInME1a = False + +# Turn off some flags for CSCRecHitD that are turned ON in default config + +process.csc2DRecHits.readBadChannels = cms.bool(False) +process.csc2DRecHits.CSCUseGasGainCorrections = cms.bool(False) +# Already defaults OFF... +## process.csc2DRecHits.CSCUseTimingCorrections = cms.bool(False) + +# Switch input for CSCRecHitD to s i m u l a t e d digis + +process.csc2DRecHits.wireDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCWireDigi") +process.csc2DRecHits.stripDigiTag = cms.InputTag("simMuonCSCDigis","MuonCSCStripDigi") + +# switch to CSCSegAlgoTC +process.cscSegments.algo_type = cms.int32( 2 ) + +# --- Activate LogVerbatim IN CSCSegment +process.MessageLogger.categories.append("CSCSegment") +process.MessageLogger.destinations = cms.untracked.vstring("cout") +process.MessageLogger.cout = cms.untracked.PSet( + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet( limit = cms.untracked.int32(0) ), + FwkReport = cms.untracked.PSet( limit = cms.untracked.int32(-1) ), + CSCSegment = cms.untracked.PSet( limit = cms.untracked.int32(-1) ) +) + +# Path and EndPath def +process.reco = cms.Path(process.csc2DRecHits * process.cscSegments) +process.endjob = cms.EndPath(process.endOfProcess) + +# Schedule definition +process.schedule = cms.Schedule(process.reco, process.endjob) From c70b154c52c351b9eb89d78d942b61e72c95d213 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:15:45 +0100 Subject: [PATCH 07/17] segment fit using smatrix factored out of algorithms --- RecoLocalMuon/CSCSegment/src/CSCSegFit.cc | 476 ++++++++++++++++++++++ RecoLocalMuon/CSCSegment/src/CSCSegFit.h | 127 ++++++ 2 files changed, 603 insertions(+) create mode 100644 RecoLocalMuon/CSCSegment/src/CSCSegFit.cc create mode 100644 RecoLocalMuon/CSCSegment/src/CSCSegFit.h diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegFit.cc b/RecoLocalMuon/CSCSegment/src/CSCSegFit.cc new file mode 100644 index 0000000000000..713c485e6a064 --- /dev/null +++ b/RecoLocalMuon/CSCSegment/src/CSCSegFit.cc @@ -0,0 +1,476 @@ +// CSCSegFit.cc +// Last mod: 03.02.2015 + +#include "RecoLocalMuon/CSCSegment/src/CSCSegFit.h" + +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" + +#include "Geometry/CSCGeometry/interface/CSCLayer.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + + +void CSCSegFit::fit(void) { + if ( fitdone() ) return; // don't redo fit unnecessarily + short n = nhits(); + switch ( n ) { + case 1: + edm::LogVerbatim("CSCSegFit") << "[CSCSegFit::fit] - cannot fit just 1 hit!!"; + break; + case 2: + fit2(); + break; + case 3: + case 4: + case 5: + case 6: + fitlsq(); + break; + default: + edm::LogVerbatim("CSCSegFit") << "[CSCSegFit::fit] - cannot fit more than 6 hits!!"; + } +} + +void CSCSegFit::fit2(void) { + + // Just join the two points + // Equation of straight line between (x1, y1) and (x2, y2) in xy-plane is + // y = mx + c + // with m = (y2-y1)/(x2-x1) + // and c = (y1*x2-x2*y1)/(x2-x1) + + CSCSetOfHits::const_iterator ih = hits_.begin(); + int il1 = (*ih)->cscDetId().layer(); + const CSCRecHit2D& h1 = (**ih); + ++ih; + int il2 = (*ih)->cscDetId().layer(); + const CSCRecHit2D& h2 = (**ih); + + // Skip if on same layer, but should be impossible :) + if (il1 == il2) { + edm::LogVerbatim("CSCSegFit") << "[CSCSegFit:fit]2 - 2 hits on same layer!!"; + return; + } + + const CSCLayer* layer1 = chamber()->layer(il1); + const CSCLayer* layer2 = chamber()->layer(il2); + + GlobalPoint h1glopos = layer1->toGlobal(h1.localPosition()); + GlobalPoint h2glopos = layer2->toGlobal(h2.localPosition()); + + // We want hit wrt chamber (and local z will be != 0) + LocalPoint h1pos = chamber()->toLocal(h1glopos); + LocalPoint h2pos = chamber()->toLocal(h2glopos); + + float dz = h2pos.z()-h1pos.z(); + + uslope_ = ( h2pos.x() - h1pos.x() ) / dz ; + vslope_ = ( h2pos.y() - h1pos.y() ) / dz ; + + float uintercept = ( h1pos.x()*h2pos.z() - h2pos.x()*h1pos.z() ) / dz; + float vintercept = ( h1pos.y()*h2pos.z() - h2pos.y()*h1pos.z() ) / dz; + intercept_ = LocalPoint( uintercept, vintercept, 0.); + + setOutFromIP(); + + //@@ NOT SURE WHAT IS SENSIBLE FOR THESE... + chi2_ = 0.; + ndof_ = 0; + + fitdone_ = true; +} + + +void CSCSegFit::fitlsq(void) { + + // Linear least-squares fit to up to 6 CSC rechits, one per layer in a CSC. + // Comments adapted from mine in original CSCSegAlgoSK algorithm. + + // Fit to the local x, y rechit coordinates in z projection + // The strip measurement controls the precision of x + // The wire measurement controls the precision of y. + // Typical precision: u (strip, sigma~200um), v (wire, sigma~1cm) + + // Set up the normal equations for the least-squares fit as a matrix equation + + // We have a vector of measurements m, which is a 2n x 1 dim matrix + // The transpose mT is (u1, v1, u2, v2, ..., un, vn) where + // ui is the strip-associated measurement and + // vi is the wire-associated measurement + // for a given rechit i. + + // The fit is to + // u = u0 + uz * z + // v = v0 + vz * z + // where u0, uz, v0, vz are the parameters to be obtained from the fit. + + // These are contained in a vector p which is a 4x1 dim matrix, and + // its transpose pT is (u0, v0, uz, vz). Note the ordering! + + // The covariance matrix for each pair of measurements is 2 x 2 and + // the inverse of this is the error matrix E. + // The error matrix for the whole set of n measurements is a diagonal + // matrix with diagonal elements the individual 2 x 2 error matrices + // (because the inverse of a diagonal matrix is a diagonal matrix + // with each element the inverse of the original.) + + // In function 'weightMatrix()', the variable 'matrix' is filled with this + // block-diagonal overall covariance matrix. Then 'matrix' is inverted to the + // block-diagonal error matrix, and returned. + + // Define the matrix A as + // 1 0 z1 0 + // 0 1 0 z1 + // 1 0 z2 0 + // 0 1 0 z2 + // .. .. .. .. + // 1 0 zn 0 + // 0 1 0 zn + + // This matrix A is set up and returned by function 'derivativeMatrix()'. + + // Then the normal equations are described by the matrix equation + // + // (AT E A)p = (AT E)m + // + // where AT is the transpose of A. + + // Call the combined matrix on the LHS, M, and that on the RHS, B: + // M p = B + + // We solve this for the parameter vector, p. + // The elements of M and B then involve sums over the hits + + // The covariance matrix of the parameters is obtained by + // (AT E A)^-1 calculated in 'covarianceMatrix()'. + + + // NOTE + // We need local position of a RecHit w.r.t. the CHAMBER + // and the RecHit itself only knows its local position w.r.t. + // the LAYER, so we must explicitly transform global position. + + + SMatrix4 M; // 4x4, init to 0 + SVector4 B; // 4x1, init to 0; + + CSCSetOfHits::const_iterator ih = hits_.begin(); + + for (ih = hits_.begin(); ih != hits_.end(); ++ih) { + const CSCRecHit2D& hit = (**ih); + const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer()); + GlobalPoint gp = layer->toGlobal(hit.localPosition()); + LocalPoint lp = chamber()->toLocal(gp); + + // Local position of hit w.r.t. chamber + double u = lp.x(); + double v = lp.y(); + double z = lp.z(); + + // Covariance matrix of local errors + SMatrixSym2 IC; // 2x2, init to 0 + + IC(0,0) = hit.localPositionError().xx(); + IC(1,1) = hit.localPositionError().yy(); + //@@ NOT SURE WHICH OFF-DIAGONAL ELEMENT MUST BE DEFINED BUT (1,0) WORKS + //@@ (and SMatrix enforces symmetry) + IC(1,0) = hit.localPositionError().xy(); + // IC(0,1) = IC(1,0); + + // Invert covariance matrix (and trap if it fails!) + bool ok = IC.Invert(); + if ( !ok ) { + edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::fit] Failed to invert covariance matrix: \n" << IC; + // return ok; //@@ SHOULD PASS THIS BACK TO CALLER? + } + + M(0,0) += IC(0,0); + M(0,1) += IC(0,1); + M(0,2) += IC(0,0) * z; + M(0,3) += IC(0,1) * z; + B(0) += u * IC(0,0) + v * IC(0,1); + + M(1,0) += IC(1,0); + M(1,1) += IC(1,1); + M(1,2) += IC(1,0) * z; + M(1,3) += IC(1,1) * z; + B(1) += u * IC(1,0) + v * IC(1,1); + + M(2,0) += IC(0,0) * z; + M(2,1) += IC(0,1) * z; + M(2,2) += IC(0,0) * z * z; + M(2,3) += IC(0,1) * z * z; + B(2) += ( u * IC(0,0) + v * IC(0,1) ) * z; + + M(3,0) += IC(1,0) * z; + M(3,1) += IC(1,1) * z; + M(3,2) += IC(1,0) * z * z; + M(3,3) += IC(1,1) * z * z; + B(3) += ( u * IC(1,0) + v * IC(1,1) ) * z; + + } + + SVector4 p; + bool ok = M.Invert(); + if (!ok ){ + edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::fit] Failed to invert matrix: \n" << M; + // return ok; //@@ SHOULD PASS THIS BACK TO CALLER? + } + else { + p = M * B; + } + + // LogTrace("CSCSegFit") << "[CSCSegFit::fit] p = " + // << p(0) << ", " << p(1) << ", " << p(2) << ", " << p(3); + + // fill member variables (note origin has local z = 0) + // intercept_ + intercept_ = LocalPoint(p(0), p(1), 0.); + + // localdir_ - set so segment points outwards from IP + uslope_ = p(2); + vslope_ = p(3); + setOutFromIP(); + + // calculate chi2 of fit + setChi2( ); + + // flag fit has been done + fitdone_ = true; + +} + + + +void CSCSegFit::setChi2(void) { + + double chsq = 0.; + + CSCSetOfHits::const_iterator ih; + for (ih = hits_.begin(); ih != hits_.end(); ++ih) { + + const CSCRecHit2D& hit = (**ih); + const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer()); + GlobalPoint gp = layer->toGlobal(hit.localPosition()); + LocalPoint lp = chamber()->toLocal(gp); + + double u = lp.x(); + double v = lp.y(); + double z = lp.z(); + + double du = intercept_.x() + uslope_ * z - u; + double dv = intercept_.y() + vslope_ * z - v; + + // LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] u, v, z = " << u << ", " << v << ", " << z; + + SMatrixSym2 IC; // 2x2, init to 0 + + IC(0,0) = hit.localPositionError().xx(); + // IC(0,1) = hit.localPositionError().xy(); + IC(1,0) = hit.localPositionError().xy(); + IC(1,1) = hit.localPositionError().yy(); + // IC(1,0) = IC(0,1); + + // LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] IC before = \n" << IC; + + // Invert covariance matrix + bool ok = IC.Invert(); + if (!ok ){ + edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::setChi2] Failed to invert covariance matrix: \n" << IC; + // return ok; + } + // LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] IC after = \n" << IC; + chsq += du*du*IC(0,0) + 2.*du*dv*IC(0,1) + dv*dv*IC(1,1); + } + + // fill member variables + chi2_ = chsq; + ndof_ = 2.*hits_.size() - 4; + + // LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] chi2 = " << chi2_ << "/" << ndof_ << " dof"; + +} + + + + +CSCSegFit::SMatrixSym12 CSCSegFit::weightMatrix() { + + bool ok = true; + + SMatrixSym12 matrix = ROOT::Math::SMatrixIdentity(); // 12x12, init to 1's on diag + + int row = 0; + + for (CSCSetOfHits::const_iterator it = hits_.begin(); it != hits_.end(); ++it) { + + const CSCRecHit2D& hit = (**it); + +// Note scaleXError allows rescaling the x error if necessary + + matrix(row, row) = scaleXError()*hit.localPositionError().xx(); + matrix(row, row+1) = hit.localPositionError().xy(); + ++row; + matrix(row, row-1) = hit.localPositionError().xy(); + matrix(row, row) = hit.localPositionError().yy(); + ++row; + } + + ok = matrix.Invert(); // invert in place + if ( !ok ) { + edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::weightMatrix] Failed to invert matrix: \n" << matrix; + // return ok; //@@ SHOULD PASS THIS BACK TO CALLER? + } + return matrix; +} + + + + +CSCSegFit::SMatrix12by4 CSCSegFit::derivativeMatrix() { + + SMatrix12by4 matrix; // 12x4, init to 0 + int row = 0; + + for( CSCSetOfHits::const_iterator it = hits_.begin(); it != hits_.end(); ++it) { + + const CSCRecHit2D& hit = (**it); + const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer()); + GlobalPoint gp = layer->toGlobal(hit.localPosition()); + LocalPoint lp = chamber()->toLocal(gp); + float z = lp.z(); + + matrix(row, 0) = 1.; + matrix(row, 2) = z; + ++row; + matrix(row, 1) = 1.; + matrix(row, 3) = z; + ++row; + } + return matrix; +} + + +void CSCSegFit::setOutFromIP() { + // Set direction of segment to point from IP outwards + // (Incorrect for particles not coming from IP, of course.) + + double dxdz = uslope_; + double dydz = vslope_; + double dz = 1./sqrt(1. + dxdz*dxdz + dydz*dydz); + double dx = dz*dxdz; + double dy = dz*dydz; + LocalVector localDir(dx,dy,dz); + + // localDir sometimes needs a sign flip + // Examine its direction and origin in global z: to point outward + // the localDir should always have same sign as global z... + + double globalZpos = ( chamber()->toGlobal( intercept_ ) ).z(); + double globalZdir = ( chamber()->toGlobal( localDir ) ).z(); + double directionSign = globalZpos * globalZdir; + localdir_ = (directionSign * localDir ).unit(); +} + + + +AlgebraicSymMatrix CSCSegFit::covarianceMatrix() { + + SMatrixSym12 weights = weightMatrix(); + SMatrix12by4 A = derivativeMatrix(); + // LogTrace("CSCSegFit") << "[CSCSegFit::covarianceMatrix] weights matrix W: \n" << weights; + // LogTrace("CSCSegFit") << "[CSCSegFit::covarianceMatrix] derivatives matrix A: \n" << A; + + // (AT W A)^-1 + // e.g. See http://www.phys.ufl.edu/~avery/fitting.html, part I + + bool ok; + SMatrixSym4 result = ROOT::Math::SimilarityT(A, weights); + // LogTrace("CSCSegFit") << "[CSCSegFit::covarianceMatrix] (AT W A): \n" << result; + ok = result.Invert(); // inverts in place + if ( !ok ) { + edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::calculateError] Failed to invert matrix: \n" << result; + // return ok; //@@ SHOULD PASS THIS BACK TO CALLER? + } + // LogTrace("CSCSegFit") << "[CSCSegFit::covarianceMatrix] (AT W A)^-1: \n" << result; + + // reorder components to match TrackingRecHit interface (CSCSegment isa TrackingRecHit) + // i.e. slopes first, then positions + AlgebraicSymMatrix flipped = flipErrors( result ); + + return flipped; +} + + +AlgebraicSymMatrix CSCSegFit::flipErrors( const SMatrixSym4& a ) { + + // The CSCSegment needs the error matrix re-arranged to match + // parameters in order (uz, vz, u0, v0) + // where uz, vz = slopes, u0, v0 = intercepts + + // LogTrace("CSCSegFit") << "[CSCSegFit::flipErrors] input: \n" << a; + + AlgebraicSymMatrix hold(4, 0. ); + + for ( short j=0; j!=4; ++j) { + for (short i=0; i!=4; ++i) { + hold(i+1,j+1) = a(i,j); // SMatrix counts from 0, AlgebraicMatrix from 1 + } + } + + // LogTrace("CSCSegFit") << "[CSCSegFit::flipErrors] after copy:"; + // LogTrace("CSCSegFit") << "(" << hold(1,1) << " " << hold(1,2) << " " << hold(1,3) << " " << hold(1,4); + // LogTrace("CSCSegFit") << " " << hold(2,1) << " " << hold(2,2) << " " << hold(2,3) << " " << hold(2,4); + // LogTrace("CSCSegFit") << " " << hold(3,1) << " " << hold(3,2) << " " << hold(3,3) << " " << hold(3,4); + // LogTrace("CSCSegFit") << " " << hold(4,1) << " " << hold(4,2) << " " << hold(4,3) << " " << hold(4,4) << ")"; + + // errors on slopes into upper left + hold(1,1) = a(2,2); + hold(1,2) = a(2,3); + hold(2,1) = a(3,2); + hold(2,2) = a(3,3); + + // errors on positions into lower right + hold(3,3) = a(0,0); + hold(3,4) = a(0,1); + hold(4,3) = a(1,0); + hold(4,4) = a(1,1); + + // must also interchange off-diagonal elements of off-diagonal 2x2 submatrices + hold(4,1) = a(1,2); + hold(3,2) = a(0,3); + hold(2,3) = a(3,0); // = a(0,3) + hold(1,4) = a(2,1); // = a(1,2) + + // LogTrace("CSCSegFit") << "[CSCSegFit::flipErrors] after flip:"; + // LogTrace("CSCSegFit") << "(" << hold(1,1) << " " << hold(1,2) << " " << hold(1,3) << " " << hold(1,4); + // LogTrace("CSCSegFit") << " " << hold(2,1) << " " << hold(2,2) << " " << hold(2,3) << " " << hold(2,4); + // LogTrace("CSCSegFit") << " " << hold(3,1) << " " << hold(3,2) << " " << hold(3,3) << " " << hold(3,4); + // LogTrace("CSCSegFit") << " " << hold(4,1) << " " << hold(4,2) << " " << hold(4,3) << " " << hold(4,4) << ")"; + + return hold; +} + +float CSCSegFit::xfit( float z ) const { + //@@ ADD THIS TO EACH ACCESSOR OF FIT RESULTS? + // if ( !fitdone() ) fit(); + return intercept_.x() + uslope_ * z; +} + +float CSCSegFit::yfit( float z ) const { + return intercept_.y() + vslope_ * z; +} + +float CSCSegFit::xdev( float x, float z ) const { + return intercept_.x() + uslope_ * z - x; +} + +float CSCSegFit::ydev( float y, float z ) const { + return intercept_.y() + vslope_ * z - y; +} + +float CSCSegFit::Rdev(float x, float y, float z) const { + return sqrt ( xdev(x,z)*xdev(x,z) + ydev(y,z)*ydev(y,z) ); +} + diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegFit.h b/RecoLocalMuon/CSCSegment/src/CSCSegFit.h new file mode 100644 index 0000000000000..97f11d41a84a0 --- /dev/null +++ b/RecoLocalMuon/CSCSegment/src/CSCSegFit.h @@ -0,0 +1,127 @@ +#ifndef CSCSegment_CSCSegFit_h +#define CSCSegment_CSCSegFit_h + +// CSCSegFit.h - Segment fitting factored out of CSC segment builders - Tim Cox +// Last mod: 03.02.2015 + +/* This as an object which is initialized by a set of rechits (2 to 6) in a + * specific CSC and has the functionality to make a least squares fit to a + * straight line in 2-dim for those rechits. + * The covariance matrix and chi2 of the fit are calculated. + * The original code made use of CLHEP matrices but this version uses + * ROOT SMatrices because they are multithreading compatible. + * Because of this, the no. of rechits that can be handled is limited to + * a maximum of 6, one per layer of a CSC. This means maximum dimensions + * can be specified at compile time and hence satisfies SMatrix constraints. + * For 2 hits of course there is no fit - just draw a straight line between them. + * Details of the algorithm are in the .cc file + * + */ + +#include +#include + +#include +#include +#include + +#include + +class CSCSegFit { + +public: + +// TYPES + + typedef std::vector CSCSetOfHits; + + // 12 x12 Symmetric + typedef ROOT::Math::SMatrix > SMatrixSym12; + + // 12 x 4 + typedef ROOT::Math::SMatrix SMatrix12by4; + + // 4 x 4 General + Symmetric + typedef ROOT::Math::SMatrix SMatrix4; + typedef ROOT::Math::SMatrix > SMatrixSym4; + + // 2 x 2 Symmetric + typedef ROOT::Math::SMatrix > SMatrixSym2; + + // 4-dim vector + typedef ROOT::Math::SVector SVector4; + + + // PUBLIC FUNCTIONS + + //@@ WANT OBJECT TO CACHE THE SET OF HITS SO CANNOT PASS BY REF + CSCSegFit( const CSCChamber* csc, CSCSetOfHits hits) : + chamber_( csc ), hits_( hits ), scaleXError_( 1.0 ), fitdone_( false ) {} + + virtual ~CSCSegFit() {} + + // Least-squares fit + void fit( void ); // fill uslope_, vslope_, intercept_ @@ FKA fitSlopes() + // Calculate covariance matrix of fitted parameters + AlgebraicSymMatrix covarianceMatrix(void); + + // Change scale factor of rechit x error + // - expert use only! + void setScaleXError ( double factor ) { scaleXError_ = factor; } + + // Fit values + float xfit( float z ) const; + float yfit( float z ) const; + + // Deviations from fit for given input (local w.r.t. chamber) + float xdev( float x, float z ) const; + float ydev ( float y, float z ) const; + float Rdev( float x, float y, float z ) const; + + // Other public functions are accessors + CSCSetOfHits hits(void) const { return hits_; } + double scaleXError(void) const { return scaleXError_; } + size_t nhits(void) const { return hits_.size(); } + double chi2(void) const { return chi2_; } + int ndof(void) const { return ndof_; } + LocalPoint intercept() const { return intercept_;} + LocalVector localdir() const { return localdir_;} + const CSCChamber* chamber() const { return chamber_; } + bool fitdone() const { return fitdone_; } + + private: + + // PRIVATE FUNCTIONS + + void fit2(void); // fit for 2 hits + void fitlsq(void); // least-squares fit for 3-6 hits + void setChi2(void); // fill chi2_ & ndof_ @@ FKA fillChiSquared() + + + protected: + + // PROTECTED FUNCTIONS - derived class needs access + + // Set segment direction 'out' from IP + void setOutFromIP(void); // fill localdir_ @@ FKA fillLocalDirection() + + SMatrix12by4 derivativeMatrix(void); + SMatrixSym12 weightMatrix(void); + AlgebraicSymMatrix flipErrors(const SMatrixSym4&); + + // PROTECTED MEMBER VARIABLES - derived class needs access + + const CSCChamber* chamber_; + CSCSetOfHits hits_; //@@ FKA protoSegment + float uslope_; //@@ FKA protoSlope_u + float vslope_; //@@ FKA protoSlope_v + LocalPoint intercept_; //@@ FKA protoIntercept + LocalVector localdir_; //@@ FKA protoDirection + double chi2_; //@@ FKA protoChi2 + int ndof_; //@@ FKA protoNDF, which was double!! + double scaleXError_; + bool fitdone_; +}; + +#endif + From e0e2c0de63fcd0c53aa7ce88914348423cf698df Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:16:35 +0100 Subject: [PATCH 08/17] conditioned version of cscsegfit factored out of cscsegalgost --- RecoLocalMuon/CSCSegment/src/CSCCondSegFit.cc | 295 ++++++++++++++++++ RecoLocalMuon/CSCSegment/src/CSCCondSegFit.h | 72 +++++ 2 files changed, 367 insertions(+) create mode 100644 RecoLocalMuon/CSCSegment/src/CSCCondSegFit.cc create mode 100644 RecoLocalMuon/CSCSegment/src/CSCCondSegFit.h diff --git a/RecoLocalMuon/CSCSegment/src/CSCCondSegFit.cc b/RecoLocalMuon/CSCSegment/src/CSCCondSegFit.cc new file mode 100644 index 0000000000000..709c00a643ad4 --- /dev/null +++ b/RecoLocalMuon/CSCSegment/src/CSCCondSegFit.cc @@ -0,0 +1,295 @@ +// CSCCondSegFit.cc +// Last mod: 23.01.2015 + +#include "RecoLocalMuon/CSCSegment/src/CSCCondSegFit.h" + +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "Geometry/CSCGeometry/interface/CSCLayer.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include + +void CSCCondSegFit::fit( bool condpass1, bool condpass2 ) { + + // See base class CSCSegFit::fit for detailed explanation of algorithm. + // This is exactly the same, except it adds two conditioning passes + // which can improve the robustness of the calculations. + + lex_.clear(); // Vector of the error matrix (only xx) + + if(condpass1 && !condpass2){ + correctTheCovX(); + if(lex_.size()!=hits_.size()){ + LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] lex_.size()!=hits_.size() ALARM! Please inform CSC DPG " << std::endl; + } + } + + SMatrix4 M; // 4x4, init to 0 + SVector4 B; // 4x1, init to 0; + + CSCSetOfHits::const_iterator ih = hits_.begin(); + + for (ih = hits_.begin(); ih != hits_.end(); ++ih) { + const CSCRecHit2D& hit = (**ih); + const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer()); + GlobalPoint gp = layer->toGlobal(hit.localPosition()); + LocalPoint lp = chamber()->toLocal(gp); + + // Local position of hit w.r.t. chamber + double u = lp.x(); + double v = lp.y(); + double z = lp.z(); + + // Covariance matrix of local errors + SMatrixSym2 IC; // 2x2, init to 0 + + // Adjust covariance matrix elements to improve numerical stability? + if(condpass1 && !condpass2){ + IC(0,0) = lex_.at(ih-hits_.begin()); + } + else { + IC(0,0) = hit.localPositionError().xx(); + } + IC(1,1) = hit.localPositionError().yy(); + //@@ NOT SURE WHICH OFF-DIAGONAL ELEMENT MUST BE DEFINED BUT (1,0) WORKS + //@@ (and SMatrix enforces symmetry) + IC(1,0) = hit.localPositionError().xy(); + // IC(0,1) = IC(1,0); + + // Correct the covariance matrix to improve numerical stability? + if(condpass2){ + correctTheCovMatrix(IC); + } + + // Invert covariance matrix (and trap if it fails!) + bool ok = IC.Invert(); + if ( !ok ) { + LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] Failed to invert covariance matrix: \n" << IC; + // return ok; //@@ SHOULD PASS THIS BACK TO CALLER? + } + + M(0,0) += IC(0,0); + M(0,1) += IC(0,1); + M(0,2) += IC(0,0) * z; + M(0,3) += IC(0,1) * z; + B(0) += u * IC(0,0) + v * IC(0,1); + + M(1,0) += IC(1,0); + M(1,1) += IC(1,1); + M(1,2) += IC(1,0) * z; + M(1,3) += IC(1,1) * z; + B(1) += u * IC(1,0) + v * IC(1,1); + + M(2,0) += IC(0,0) * z; + M(2,1) += IC(0,1) * z; + M(2,2) += IC(0,0) * z * z; + M(2,3) += IC(0,1) * z * z; + B(2) += ( u * IC(0,0) + v * IC(0,1) ) * z; + + M(3,0) += IC(1,0) * z; + M(3,1) += IC(1,1) * z; + M(3,2) += IC(1,0) * z * z; + M(3,3) += IC(1,1) * z * z; + B(3) += ( u * IC(1,0) + v * IC(1,1) ) * z; + + } + + SVector4 p; + bool ok = M.Invert(); + if (!ok ){ + LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] Failed to invert matrix: \n" << M; + // return ok; //@@ SHOULD PASS THIS BACK TO CALLER? + } + else { + p = M * B; + } + + // LogTrace("CSCSegFit") << "[CSCSegFit::fit] p = " + // << p(0) << ", " << p(1) << ", " << p(2) << ", " << p(3); + + // fill member variables (note origin has local z = 0) + intercept_ = LocalPoint(p(0), p(1), 0.); + + // localdir_ + uslope_ = p(2); + vslope_ = p(3); + setOutFromIP(); + + // calculate chi2 of fit + setChi2( condpass1, condpass2 ); + +} + + + +void CSCCondSegFit::setChi2( bool condpass1, bool condpass2 ) { + + double chsq = 0.; + + CSCSetOfHits::const_iterator ih; + for (ih = hits_.begin(); ih != hits_.end(); ++ih) { + + const CSCRecHit2D& hit = (**ih); + const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer()); + GlobalPoint gp = layer->toGlobal(hit.localPosition()); + LocalPoint lp = chamber()->toLocal(gp); + + double u = lp.x(); + double v = lp.y(); + double z = lp.z(); + + double du = intercept_.x() + uslope_ * z - u; + double dv = intercept_.y() + vslope_ * z - v; + + // LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] u, v, z = " << u << ", " << v << ", " << z; + + SMatrixSym2 IC; // 2x2, init to 0 + + if(condpass1 && !condpass2){ + IC(0,0) = lex_.at(ih-hits_.begin()); + } + else{ + IC(0,0) = hit.localPositionError().xx(); + } + // IC(0,1) = hit.localPositionError().xy(); + IC(1,0) = hit.localPositionError().xy(); + IC(1,1) = hit.localPositionError().yy(); + // IC(1,0) = IC(0,1); + + // LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] IC before = \n" << IC; + + /// Correct the cov matrix + if(condpass2){ + correctTheCovMatrix(IC); + } + + // Invert covariance matrix + bool ok = IC.Invert(); + if (!ok ){ + LogTrace("CSCSegment|CSC") << "[CSCSegFit::setChi2] Failed to invert covariance matrix: \n" << IC; + // return ok; + } + // LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] IC after = \n" << IC; + chsq += du*du*IC(0,0) + 2.*du*dv*IC(0,1) + dv*dv*IC(1,1); + } + + // fill member variables + chi2_ = chsq; + ndof_ = 2.*hits_.size() - 4; + + // LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] chi2 = " << chi2_ << "/" << ndof_ << " dof"; + +} + + + +void CSCCondSegFit::correctTheCovX(void){ + std::vector uu, vv, zz; // Vectors of coordinates + + lex_.clear(); // x component of local error for each hit + + double sum_U_err=0.0; + double sum_Z_U_err=0.0; + double sum_Z2_U_err=0.0; + double sum_U_U_err=0.0; + double sum_UZ_U_err=0.0; + std::vector chiUZind; + std::vector::iterator chiContribution; + double chiUZ=0.0; + CSCSetOfHits::const_iterator ih = hits_.begin(); + for (ih = hits_.begin(); ih != hits_.end(); ++ih) { + const CSCRecHit2D& hit = (**ih); + lex_.push_back(hit.localPositionError().xx()); + + const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer()); + GlobalPoint gp = layer->toGlobal(hit.localPosition()); + LocalPoint lp = chamber()->toLocal(gp); + // Local position of hit w.r.t. chamber + double u = lp.x(); + double v = lp.y(); + double z = lp.z(); + uu.push_back(u); + vv.push_back(v); + zz.push_back(z); + // Prepare the sums for the standard linear fit + sum_U_err += 1./lex_.back(); + sum_Z_U_err += z/lex_.back(); + sum_Z2_U_err += (z*z)/lex_.back(); + sum_U_U_err += u/lex_.back(); + sum_UZ_U_err += (u*z)/lex_.back(); + } + + /// Make a one dimensional fit in U-Z plane (i.e. chamber local x-z) + /// U=U0+UZ*Z fit parameters + + double denom=sum_U_err*sum_Z2_U_err-pow(sum_Z_U_err,2); + double U0=(sum_Z2_U_err*sum_U_U_err-sum_Z_U_err*sum_UZ_U_err)/denom; + double UZ=(sum_U_err*sum_UZ_U_err-sum_Z_U_err*sum_U_U_err)/denom; + + /// Calculate the fit line trace + /// Calculate one dimensional chi^2 and normalize the errors if needed + + for(unsigned i=0; i=chi2Norm_){ + double chi2uCorrection = chiUZ/chi2Norm_; + for(unsigned i=0; icondNumberCorr2))|| + ((diag2>condNumberCorr2)&&(detCov + +#include + +class CSCCondSegFit : public CSCSegFit { + +public: + + CSCCondSegFit( const edm::ParameterSet& ps, const CSCChamber* csc, const CSCSetOfHits& hits) : + CSCSegFit( csc, hits ), + worstHit_( 0 ), + chi2Norm_ ( ps.getParameter("NormChi2Cut2D") ), + condSeed1_ ( ps.getParameter("SeedSmall") ), + condSeed2_ ( ps.getParameter("SeedBig") ), + covToAnyNumber_ ( ps.getParameter("ForceCovariance") ), + covToAnyNumberAll_ ( ps.getParameter("ForceCovarianceAll") ), + covAnyNumber_ ( ps.getParameter("Covariance") ) {} + + ~CSCCondSegFit() {} + + // The fit - override base class version with this version + // which passes in bool flags for up to two extra conditioning passes + void fit( bool condpass1 = false, bool condpass2 = false ); // fill uslope_, vslope_, intercept_ + + int worstHit( void ) { return worstHit_; } + + private: + + // Rest can all be private since we don't plan on more derived classes + + // PRIVATE MEMBER FUNCTIONS + + // Calculate chi2 - override base class version with this version + // which passes in bool flags for up to two extra conditioning passes + void setChi2( bool condpass1, bool condpass2 ); // fill chi2_ & ndof_ + void correctTheCovMatrix(CSCSegFit::SMatrixSym2& IC); + void correctTheCovX(void); + + + // EXTRA MEMBER VARIABLES + + int worstHit_; //@@ FKA maxContrIndex + + // Parameters related to adjustment for numerical robustness + std::vector lex_; //@@ FKA e_Cxx; LOCAL ERROR x COMPONENT FOR EACH HIT + double chi2Norm_; //@@ FKA chi2Norm_2D_ + + // PSet values that might reasonably be accessed directly + // since used ONLY in correctTheCovMatrix: + //@@ the comments on following parameters don't help me understand them + double condSeed1_, condSeed2_; /// The correction parameters + bool covToAnyNumber_; /// Allow to use any number for covariance (by hand) + bool covToAnyNumberAll_; /// Allow to use any number for covariance for all RecHits + double covAnyNumber_; /// The number to force the Covariance + +}; + +#endif + From 0adc7501c2cdb95fb88c9c27481ed5ca460adaa5 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:18:28 +0100 Subject: [PATCH 09/17] logdebug to logverbatim for convenience when dumping --- RecoLocalMuon/CSCSegment/src/CSCSegmentBuilder.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegmentBuilder.cc b/RecoLocalMuon/CSCSegment/src/CSCSegmentBuilder.cc index 44f5b416eebdb..a1df7d49a2b97 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegmentBuilder.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegmentBuilder.cc @@ -53,7 +53,7 @@ CSCSegmentBuilder::CSCSegmentBuilder(const edm::ParameterSet& ps) : geom_(0) { for (size_t j=0; j create(algoName, segAlgoPSet[algoToType[j]-1]); - LogDebug("CSCSegment|CSC")<< "using algorithm #" << algoToType[j] << " for chamber type " << chType[j]; + edm::LogVerbatim("CSCSegment|CSC")<< "using algorithm #" << algoToType[j] << " for chamber type " << chType[j]; } } From 2f701997404a6701ef415c19ba2b00002982c3e0 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:19:31 +0100 Subject: [PATCH 10/17] back to stream again --- RecoLocalMuon/CSCSegment/src/CSCSegmentProducer.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegmentProducer.h b/RecoLocalMuon/CSCSegment/src/CSCSegmentProducer.h index 1c82fe8bd221f..2db9a09c82c65 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegmentProducer.h +++ b/RecoLocalMuon/CSCSegment/src/CSCSegmentProducer.h @@ -4,12 +4,11 @@ /** \class CSCSegmentProducer * Produces a collection of CSCSegment's in endcap muon CSCs. * - * \author M. Sani */ -#include +#include "FWCore/Framework/interface/ConsumesCollector.h" #include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -17,7 +16,7 @@ class CSCSegmentBuilder; -class CSCSegmentProducer : public edm::EDProducer { +class CSCSegmentProducer : public edm::stream::EDProducer<> { public: /// Constructor explicit CSCSegmentProducer(const edm::ParameterSet&); From 09df29d80b8f06e6d0940016fbfa8ab37fa6b9e0 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Thu, 19 Feb 2015 12:20:54 +0100 Subject: [PATCH 11/17] factored out segment fit which now uses smatrix not clhep --- RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.cc | 464 ++++------- RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h | 25 +- RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc | 495 ++---------- RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.h | 26 +- RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc | 702 ++++------------- RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.h | 87 +-- .../CSCSegment/src/CSCSegAlgoShowering.cc | 380 ++------- .../CSCSegment/src/CSCSegAlgoShowering.h | 19 +- RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc | 723 ++++-------------- RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.h | 53 +- 10 files changed, 671 insertions(+), 2303 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.cc index ff365045a2bba..02babf7bdeff9 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.cc @@ -2,13 +2,16 @@ * \file CSCSegAlgoDF.cc * * \author Dominique Fortin - UCR + * + * Last update: 17.02.2015 - Tim Cox - use CSCSegFit for least-squares fit + * @@ Seems to find very few segments so many candidates must be being rejected. + * */ #include "CSCSegAlgoDF.h" +#include "CSCSegFit.h" #include "Geometry/CSCGeometry/interface/CSCLayer.h" -// For clhep Matrix::solve -#include "DataFormats/CLHEP/interface/AlgebraicObjects.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" @@ -29,7 +32,8 @@ /* Constructor * */ -CSCSegAlgoDF::CSCSegAlgoDF(const edm::ParameterSet& ps) : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoDF") { +CSCSegAlgoDF::CSCSegAlgoDF(const edm::ParameterSet& ps) + : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoDF"), sfit_(0) { debug = ps.getUntrackedParameter("CSCSegmentDebug"); minLayersApart = ps.getParameter("minLayersApart"); @@ -101,7 +105,10 @@ std::vector CSCSegAlgoDF::run(const CSCChamber* aChamber, const Cham * to build segments with the wrong starting points. In the DF algorithm, the endpoints * are tested as the best starting points in a 2nd loop. * - * Also, only a certain muonsPerChamberMax maximum number of segments can be produced in the chamber + * Also, a maximum of 5 segments is allowed in the chamber (and then it just gives up.) + * + * @@ There are 7 return's in this function! + * */ std::vector CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _rechits) { @@ -111,21 +118,32 @@ std::vector CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _ segmentInChamber.clear(); unsigned nHitInChamber = rechits.size(); + + // std::cout << "[CSCSegAlgoDF::buildSegments] address of chamber = " << theChamber << std::endl; + // std::cout << "[CSCSegAlgoDF::buildSegments] starting in " << theChamber->id() + // << " with " << nHitInChamber << " rechits" << std::endl; + + // Return #1 - OK, there aren't enough rechits to build a segment 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(); + size_t nLayers = 0; + size_t old_layer = 0; + for ( size_t i = 0; i < nHitInChamber; ++i ) { + size_t this_layer = rechits[i]->cscDetId().layer(); + // std::cout << "[CSCSegAlgoDF::buildSegments] this_layer = " << this_layer << std::endl; layerIndex[i] = this_layer; + // std::cout << "[CSCSegAlgoDF::buildSegments] layerIndex[" << i << "] = " << layerIndex[i] << std::endl; if ( this_layer != old_layer ) { old_layer = this_layer; - nLayers++; + ++nLayers; } } + + // std::cout << "[CSCSegAlgoDF::buildSegments] layers with rechits = " << nLayers << std::endl; + // Return #2 - OK, there aren't enough hit layers to build a segment if ( nLayers < 3 ) return segmentInChamber; double z1 = theChamber->layer(1)->position().z(); @@ -144,26 +162,35 @@ std::vector CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _ } } + // std::cout << "[CSCSegAlgoDF::buildSegments] rechits have been ordered" << std::endl; + // Showering muon if ( preClustering && int(nHitInChamber) > nHitsPerClusterIsShower && nLayers > 2 ) { + + // std::cout << "[CSCSegAlgoDF::buildSegments] showering block" << std::endl; + CSCSegment segShower = showering_->showerSeg(theChamber, rechits); + // Return #3 - OK, this is now 'effectve' rechits // Make sure have at least 3 hits... if ( segShower.nRecHits() < 3 ) return segmentInChamber; segmentInChamber.push_back(segShower); + if (debug) dumpSegment( segShower ); + // Return #4 - OK, only get one try at building a segment from shower 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(); - + + // std::cout << "[CSCSegAlgoDF::buildSegments] entering rechit loop" << std::endl; + // 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; @@ -190,21 +217,25 @@ std::vector CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _ // 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 ; + float slope_u = (lp2.x() - lp1.x())/dz ; + float slope_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; + if (fabs(slope_v) > tanThetaMax) continue; + if (fabs(slope_u) > tanPhiMax ) continue; protoSegment.push_back(h1); protoSegment.push_back(h2); - + + // std::cout << "[CSCSegAlgoDF::buildSegments] about to fit 2 hits on layers " + // << layer1 << " and " << layer2 << std::endl; + + // protoSegment has just 2 hits - but need to create a CSCSegFit to hold it in case + // we fail to add any more hits + updateParameters(); + // Try adding hits to proto segment tryAddingHitsToSegment(rechits, i1, i2, layerIndex); @@ -217,35 +248,23 @@ std::vector CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _ if ( Pruning && segok ) pruneFromResidual(); // Check if segment satisfies chi2 requirement - if (protoChi2 > chi2Max) segok = false; + if ( sfit_->chi2() > 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); + // Create an actual CSCSegment - retrieve all info from the current fit + CSCSegment temp(sfit_->hits(), sfit_->intercept(), sfit_->localdir(), + sfit_->covarianceMatrix(), sfit_->chi2()); + // std::cout << "[CSCSegAlgoDF::buildSegments] about to delete sfit= = " << sfit_ << std::endl; + delete sfit_; + sfit_ = 0; // avoid possibility of attempting a second delete later segmentInChamber.push_back(temp); + if (debug) dumpSegment( temp ); + // Return #5 - OK, fewer than 3 rechits not on this segment left in chamber if (nHitInChamber-protoSegment.size() < 3) return segmentInChamber; + // Return $6 - already have more than 4 segments in this chamber if (segmentInChamber.size() > 4) return segmentInChamber; // Flag used hits @@ -253,6 +272,7 @@ std::vector CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _ } } } + // Return #7 return segmentInChamber; } @@ -266,7 +286,7 @@ std::vector CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _ void CSCSegAlgoDF::tryAddingHitsToSegment( const ChamberHitContainer& rechits, const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2, - const LayerIndex& layerIndex ) { + const LayerIndex& layerIndex) { /* Iterate over the layers with hits in the chamber * Skip the layers containing the segment endpoints on first pass, but then @@ -274,23 +294,37 @@ void CSCSegAlgoDF::tryAddingHitsToSegment( const ChamberHitContainer& rechits, * if segment has >2 hits. Once a hit is added to a layer, don't replace it * until 2nd iteration */ + + +// std::cout << "[CSCSegAlgoDF::tryAddingHitsToSegment] entering" +// << " with rechits.size() = " << rechits.size() << std::endl; ChamberHitContainerCIt ib = rechits.begin(); ChamberHitContainerCIt ie = rechits.end(); closeHits.clear(); - for ( ChamberHitContainerCIt i = ib; i != ie; ++i ) { + // int counter1 = 0; + // int counter2 = 0; + + for ( ChamberHitContainerCIt i = ib; i != ie; ++i ) { + // std::cout << "counter1 = " << ++counter1 << std::endl; if (i == i1 || i == i2 ) continue; if ( usedHits[i-ib] ) continue; // Don't use hits already part of a segment. + // std::cout << "counter2 = " << ++counter2 << std::endl; const CSCRecHit2D* h = *i; int layer = layerIndex[i-ib]; int layer1 = layerIndex[i1-ib]; int layer2 = layerIndex[i2-ib]; + // std::cout << "[CSCSegAlgoDF::tryAddingHitsToSegment] layer, layer1, layer2 = " + // << layer << ", " << layer1 << ", " << layer2 << std::endl; + // Low multiplicity case + // only adds hit to protoSegment if no hit on layer already; otherwise adds to closeHits if (rechits.size() < 9) { + // std::cout << "low mult" << std::endl; if ( isHitNearSegment( h ) ) { if ( !hasHitOnLayer(layer) ) { addHit(h, layer); @@ -300,13 +334,20 @@ void CSCSegAlgoDF::tryAddingHitsToSegment( const ChamberHitContainer& rechits, } // High multiplicity case + // only adds hit to protoSegment if no hit on layer already AND then refits; otherwise adds to closeHits } else { + // std::cout << "high mult" << std::endl; if ( isHitNearSegment( h ) ) { + // std::cout << "near seg" << std::endl; if ( !hasHitOnLayer(layer) ) { - addHit(h, layer); - updateParameters(); + // std::cout << "no hit on layer" << std::endl; + if ( addHit(h, layer) ) { + // std::cout << "update fit" << std::endl; + updateParameters(); + } // Don't change the starting points at this stage !!! } else { + // std::cout << "already hit on layer" << std::endl; closeHits.push_back(h); if (layer != layer1 && layer != layer2 ) compareProtoSegment(h, layer); } @@ -315,12 +356,13 @@ void CSCSegAlgoDF::tryAddingHitsToSegment( const ChamberHitContainer& rechits, } if ( int(protoSegment.size()) < 3) return; - + // std::cout << "final fit" << std::endl; updateParameters(); // 2nd pass to remove biases // This time, also consider changing the endpoints for ( ChamberHitContainerCIt i = closeHits.begin() ; i != closeHits.end(); ++i ) { + // std::cout << "2nd pass" << std::endl; const CSCRecHit2D* h = *i; int layer = (*i)->cscDetId().layer(); compareProtoSegment(h, layer); @@ -331,33 +373,34 @@ void CSCSegAlgoDF::tryAddingHitsToSegment( const ChamberHitContainer& rechits, /* isHitNearSegment * - * Compare rechit with expected position from proto_segment + * Compare rechit with expected position from protosegment */ -bool CSCSegAlgoDF::isHitNearSegment( const CSCRecHit2D* hit) const { +bool CSCSegAlgoDF::isHitNearSegment( const CSCRecHit2D* hit ) const { const CSCLayer* layer = theChamber->layer(hit->cscDetId().layer()); // hit phi position in global coordinates GlobalPoint Hgp = layer->toGlobal(hit->localPosition()); - double Hphi = Hgp.phi(); + float Hphi = Hgp.phi(); if (Hphi < 0.) Hphi += 2.*M_PI; LocalPoint Hlp = theChamber->toLocal(Hgp); - double z = Hlp.z(); + float z = Hlp.z(); + + float LocalX = sfit_->xfit(z); + float LocalY = sfit_->yfit(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(); + float Sphi = Sgp.phi(); if (Sphi < 0.) Sphi += 2.*M_PI; - double R = sqrt(Sgp.x()*Sgp.x() + Sgp.y()*Sgp.y()); + float R = sqrt(Sgp.x()*Sgp.x() + Sgp.y()*Sgp.y()); - double deltaPhi = Sphi - Hphi; + float 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; + float RdeltaPhi = R * deltaPhi; if (RdeltaPhi < dRPhiFineMax && deltaPhi < dPhiFineMax ) return true; @@ -371,19 +414,23 @@ bool CSCSegAlgoDF::isHitNearSegment( const CSCRecHit2D* hit) const { * */ bool CSCSegAlgoDF::addHit(const CSCRecHit2D* aHit, int layer) { + + + // std::cout << "[CSCSegAlgoDF::addHit] on layer " << layer << " to protoSegment.size() = " + // << protoSegment.size() << std::endl; // 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; + if ( protoSegment.size() > 5 ) return false; //@@ can only have 6 hits at most // Test that we are not trying to add the same hit again - for ( ChamberHitContainer::const_iterator it = protoSegment.begin(); it != protoSegment.end(); it++ ) + for ( ChamberHitContainer::const_iterator it = protoSegment.begin(); it != protoSegment.end(); ++it ) if ( aHit == (*it) ) return false; protoSegment.push_back(aHit); - return ok; + return true; } @@ -394,114 +441,29 @@ bool CSCSegAlgoDF::addHit(const CSCRecHit2D* aHit, int layer) { */ void CSCSegAlgoDF::updateParameters() { - // 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; + // Delete existing CSCSegFit, create a new one and make the fit + // Uses internal variables - theChamber & protoSegment + + // std::cout << "[CSCSegAlgoDF::updateParameters] about to delete sfit_ = " << sfit_ << std::endl; + delete sfit_; + // std::cout << "[CSCSegAlgoDF::updateParameters] protoSegment.size() = " + // << protoSegment.size() << std::endl; + // std::cout << "[CSCSegAlgoDF::updateParameters] theChamber = " << theChamber << std::endl; + sfit_ = new CSCSegFit( theChamber, protoSegment ); + // std::cout << "[CSCSegAlgoDF::updateParameters] new sfit_ = " << sfit_ << std::endl; + sfit_->fit(); } - /* hasHitOnLayer * * Just make sure hit to be added to layer comes from different layer than those in proto segment */ bool CSCSegAlgoDF::hasHitOnLayer(int layer) const { + + // std::cout << "[CSCSegAlgoDF::hasHitOnLayer] on layer " << layer << std::endl; + + // Is there already a hit on this layer? for ( ChamberHitContainerCIt it = protoSegment.begin(); it != protoSegment.end(); it++ ) if ( (*it)->cscDetId().layer() == layer ) return true; @@ -517,15 +479,10 @@ bool CSCSegAlgoDF::hasHitOnLayer(int layer) const { * */ void CSCSegAlgoDF::compareProtoSegment(const CSCRecHit2D* h, int layer) { - - // 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; - + + + // std::cout << "[CSCSegAlgoDF::compareProtoSegment] for hit on layer " << layer + // << " with protoSegment.size() = " << protoSegment.size() << std::endl; // Try adding the hit to existing segment, and remove old one existing in same layer ChamberHitContainer::iterator it; @@ -536,17 +493,27 @@ void CSCSegAlgoDF::compareProtoSegment(const CSCRecHit2D* h, int layer) { ++it; } } + + // std::cout << "[CSCSegAlgoDF::compareProtoSegment] about to add hit on layer " << layer + // << " with protoSegment.size() = " << protoSegment.size() << std::endl; + 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; + CSCSegFit* newfit = 0; + if ( ok ) { + newfit = new CSCSegFit( theChamber, protoSegment ); + // std::cout << "[CSCSegAlgoDF::compareProtoSegment] newfit = " << newfit << std::endl; + newfit->fit(); + } + if ( !ok || (newfit->chi2() > sfit_->chi2()) ) { + // std::cout << "[CSCSegAlgoDF::compareProtoSegment] about to delete newfit = " << newfit << std::endl; + delete newfit; // failed to add a hit or new fit is worse + } + else { + // std::cout << "[CSCSegAlgoDF::compareProtoSegment] about to delete sfit_ = " << sfit_ << std::endl; + delete sfit_; // new fit is better + sfit_ = newfit; + // std::cout << "[CSCSegAlgoDF::compareProtoSegment] reset sfit_ = " << sfit_ << std::endl; } } @@ -580,101 +547,12 @@ void CSCSegAlgoDF::flagHitsAsUsed(const ChamberHitContainer& rechitsInChamber) { } -AlgebraicSymMatrix CSCSegAlgoDF::weightMatrix() const { - - std::vector::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; -} - - -CLHEP::HepMatrix CSCSegAlgoDF::derivativeMatrix() const { - - 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; -} - -/* calculateError - * - */ -AlgebraicSymMatrix CSCSegAlgoDF::calculateError() const { - - 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::flipErrors( AlgebraicSymMatrix& a ) const { - - // 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 - -} - - // Try to clean up segments by quickly looking at residuals -void CSCSegAlgoDF::pruneFromResidual(){ +void CSCSegAlgoDF::pruneFromResidual(void){ // Only prune if have at least 5 hits if ( protoSegment.size() < 5 ) return ; - // Now Study residuals float maxResidual = 0.; @@ -692,14 +570,7 @@ void CSCSegAlgoDF::pruneFromResidual(){ 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); + float residual = sfit_->Rdev(lp.x(), lp.y(), lp.z()); sumResidual += residual; nHits++; @@ -732,38 +603,21 @@ void CSCSegAlgoDF::pruneFromResidual(){ protoSegment.push_back(*ih); } - // Update segment parameters + // Compute segment parameters updateParameters(); } - -/* - * Order the hits such that 2nd one is closest in x,y to first seed hit in global coordinates - */ -void CSCSegAlgoDF::orderSecondSeed( GlobalPoint gp1, - const ChamberHitContainerCIt i1, - const ChamberHitContainerCIt i2, - const ChamberHitContainer& rechits, - const LayerIndex& layerIndex ) { - - 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::dumpSegment( const CSCSegment& seg ) const { + + edm::LogVerbatim("CSCSegment") << "CSCSegment in " << theChamber->id() + << "\nlocal position = " << seg.localPosition() + << "\nerror = " << seg.localPositionError() + << "\nlocal direction = " << seg.localDirection() + << "\nerror =" << seg.localDirectionError() + << "\ncovariance matrix" + << seg.parametersError() + << "chi2/ndf = " << seg.chi2() << "/" << seg.degreesOfFreedom() + << "\n#rechits = " << seg.specificRecHits().size() + << "\ntime = " << seg.time(); } diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h index 23a827d0fa567..6a3be86c4070c 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h @@ -43,6 +43,7 @@ class CSCSegAlgoPreClustering; class CSCSegAlgoShowering; +class CSCSegFit; class CSCSegAlgoDF : public CSCSegmentAlgorithm { @@ -102,28 +103,14 @@ class CSCSegAlgoDF : public CSCSegmentAlgorithm { /** * Prune bad segment from the worse hit based on residuals */ - void pruneFromResidual(); + void pruneFromResidual(void); - - /** - * Order the hits on the 2nd layer for seed building - */ - void orderSecondSeed( GlobalPoint gp1, - const ChamberHitContainerCIt i1, - const ChamberHitContainerCIt i2, - const ChamberHitContainer& rechits, - const LayerIndex& layerIndex ); - - bool isHitNearSegment(const CSCRecHit2D* h) const; bool addHit(const CSCRecHit2D* hit, int layer); void updateParameters(void); bool hasHitOnLayer(int layer) const; void compareProtoSegment(const CSCRecHit2D* h, int layer); - CLHEP::HepMatrix derivativeMatrix(void) const; - AlgebraicSymMatrix weightMatrix(void) const; - AlgebraicSymMatrix calculateError(void) const; - void flipErrors(AlgebraicSymMatrix&) const; + void dumpSegment( const CSCSegment& seg ) const; // Member variables const std::string myName; @@ -134,11 +121,6 @@ class CSCSegAlgoDF : public CSCSegmentAlgorithm { ChamberHitContainer protoSegment; ChamberHitContainer secondSeedHits; - float protoSlope_u; - float protoSlope_v; - LocalPoint protoIntercept; - double protoChi2; - LocalVector protoDirection; // input from .cfi file bool debug; @@ -160,6 +142,7 @@ class CSCSegAlgoDF : public CSCSegmentAlgorithm { CSCSegAlgoPreClustering* preCluster_; CSCSegAlgoShowering* showering_; + CSCSegFit* sfit_; }; diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc index d531f1062a6be..06dab0a4bfdaf 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc @@ -1,13 +1,11 @@ /** * \file CSCSegAlgoSK.cc * - * \author M. Sani */ #include "CSCSegAlgoSK.h" +#include "CSCSegFit.h" #include "Geometry/CSCGeometry/interface/CSCLayer.h" -// For clhep Matrix::solve -#include "DataFormats/CLHEP/interface/AlgebraicObjects.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -18,8 +16,8 @@ #include #include -CSCSegAlgoSK::CSCSegAlgoSK(const edm::ParameterSet& ps) : CSCSegmentAlgorithm(ps), - myName("CSCSegAlgoSK") { +CSCSegAlgoSK::CSCSegAlgoSK(const edm::ParameterSet& ps) + : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoSK"), sfit_(0) { debugInfo = ps.getUntrackedParameter("verboseInfo"); @@ -77,11 +75,8 @@ std::vector CSCSegAlgoSK::buildSegments(const ChamberHitContainer& u } } - if (debugInfo) { - // dump after sorting - dumpHits(rechits); - } - + // if (debugInfo) dumpHits(rechits); + if (rechits.size() < 2) { LogDebug("CSC") << myName << ": " << rechits.size() << " hit(s) in chamber is not enough to build a segment.\n"; @@ -106,6 +101,9 @@ std::vector CSCSegAlgoSK::buildSegments(const ChamberHitContainer& u // Define buffer for segments we build std::vector segments; + + // This is going to point to fits to hits, and its content will be used to create a CSCSegment + sfit_ = 0; ChamberHitContainerCIt ib = rechits.begin(); ChamberHitContainerCIt ie = rechits.end(); @@ -144,7 +142,8 @@ std::vector CSCSegAlgoSK::buildSegments(const ChamberHitContainer& u GlobalPoint gp2 = l2->toGlobal(h2->localPosition()); LogDebug("CSC") << "start new segment from hits " << "h1: " << gp1 << " - h2: " << gp2 << "\n"; - + + //@@ TRY ADDING A HIT - AND FIT if (!addHit(h1, layer1)) { LogDebug("CSC") << " failed to add hit h1\n"; continue; @@ -155,7 +154,8 @@ std::vector CSCSegAlgoSK::buildSegments(const ChamberHitContainer& u continue; } - tryAddingHitsToSegment(rechits, used, layerIndex, i1, i2); + // Can only add hits if already have a segment + if ( sfit_ ) tryAddingHitsToSegment(rechits, used, layerIndex, i1, i2); // Check no. of hits on segment, and if enough flag them as used // and store the segment @@ -169,18 +169,13 @@ std::vector CSCSegAlgoSK::buildSegments(const ChamberHitContainer& u LogDebug("CSC") << "No segment has been found !!!\n"; } else { - - // calculate error matrix... - AlgebraicSymMatrix errors = calculateError(); - - // but reorder components to match what's required by TrackingRecHit interface - // i.e. slopes first, then positions - - flipErrors( errors ); - - CSCSegment temp(proto_segment, theOrigin, theDirection, errors, theChi2); - + // Create an actual CSCSegment - retrieve all info from the fit + CSCSegment temp(sfit_->hits(), sfit_->intercept(), + sfit_->localdir(), sfit_->covarianceMatrix(), sfit_->chi2() ); + delete sfit_; + sfit_ = 0; LogDebug("CSC") << "Found a segment !!!\n"; + if ( debugInfo ) dumpSegment( temp ); segments.push_back(temp); } } @@ -308,11 +303,18 @@ bool CSCSegAlgoSK::isHitNearSegment(const CSCRecHit2D* h) const { float CSCSegAlgoSK::phiAtZ(float z) const { + if ( !sfit_ ) { + edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoSK::phiAtZ] Segment fit undefined"; + return 0.; + } + // Returns a phi in [ 0, 2*pi ) const CSCLayer* l1 = theChamber->layer((*(proto_segment.begin()))->cscDetId().layer()); - GlobalPoint gp = l1->toGlobal(theOrigin); - GlobalVector gv = l1->toGlobal(theDirection); - + GlobalPoint gp = l1->toGlobal(sfit_->intercept()); + GlobalVector gv = l1->toGlobal(sfit_->localdir()); + + LogTrace("CSCSegment") << "[CSCSegAlgoSK::phiAtZ] Global intercept = " << gp << ", direction = " << gv; + float x = gp.x() + (gv.x()/gv.z())*(z - gp.z()); float y = gp.y() + (gv.y()/gv.z())*(z - gp.z()); float phi = atan2(y, x); @@ -326,13 +328,13 @@ void CSCSegAlgoSK::dumpHits(const ChamberHitContainer& rechits) const { // Dump positions of RecHit's in each CSCChamber ChamberHitContainerCIt it; - edm::LogInfo("CSC") << "CSCChamber rechit dump.\n"; + edm::LogInfo("CSCSegment") << "CSCChamber rechit dump.\n"; for(it=rechits.begin(); it!=rechits.end(); it++) { const CSCLayer* l1 = theChamber->layer((*it)->cscDetId().layer()); GlobalPoint gp1 = l1->toGlobal((*it)->localPosition()); - edm::LogInfo("CSC") << "Global pos.: " << gp1 << ", phi: " << gp1.phi() << ". Local position: " + edm::LogInfo("CSCSegment") << "Global pos.: " << gp1 << ", phi: " << gp1.phi() << ". Local position: " << (*it)->localPosition() << ", phi: " << (*it)->localPosition().phi() << ". Layer: " << (*it)->cscDetId().layer() << "\n"; @@ -343,7 +345,7 @@ bool CSCSegAlgoSK::isSegmentGood(const ChamberHitContainer& rechitsInChamber) co // If the chamber has 20 hits or fewer, require at least 3 hits on segment // If the chamber has >20 hits require at least 4 hits - //@@ THESE VALUES SHOULD BECOME PARAMETERS? + bool ok = false; unsigned int iadd = ( rechitsInChamber.size() > 20 )? 1 : 0; @@ -379,295 +381,24 @@ bool CSCSegAlgoSK::addHit(const CSCRecHit2D* aHit, int layer) { // Return false if there is already a hit on the same layer, or insert failed. ChamberHitContainer::const_iterator it; - + for(it = proto_segment.begin(); it != proto_segment.end(); it++) if (((*it)->cscDetId().layer() == layer) && (aHit != (*it))) - return false; + return false; proto_segment.push_back(aHit); + + // make a fit updateParameters(); - return true; } -void CSCSegAlgoSK::updateParameters() { - - // Note that we need local position of a RecHit w.r.t. the CHAMBER - // and the RecHit itself only knows its local position w.r.t. - // the LAYER, so need to explicitly transform to global position. - - // no. of hits in the RecHitsOnSegment - // By construction this is the no. of layers with hits - // since we allow just one hit per layer in a segment. - - int nh = proto_segment.size(); - - // First hit added to a segment must always fail here - if (nh < 2) - return; - - if (nh == 2) { - - // Once we have two hits we can calculate a straight line - // (or rather, a straight line for each projection in xz and yz.) - ChamberHitContainer::const_iterator ih = proto_segment.begin(); - int il1 = (*ih)->cscDetId().layer(); - const CSCRecHit2D& h1 = (**ih); - ++ih; - int il2 = (*ih)->cscDetId().layer(); - const CSCRecHit2D& h2 = (**ih); - - //@@ Skip if on same layer, but should be impossible - if (il1 == il2) - return; - - const CSCLayer* layer1 = theChamber->layer(il1); - const CSCLayer* layer2 = theChamber->layer(il2); - - GlobalPoint h1glopos = layer1->toGlobal(h1.localPosition()); - GlobalPoint h2glopos = layer2->toGlobal(h2.localPosition()); - - // localPosition is position of hit wrt layer (so local z = 0) - theOrigin = h1.localPosition(); - - // We want hit wrt chamber (and local z will be != 0) - LocalPoint h1pos = theChamber->toLocal(h1glopos); - LocalPoint h2pos = theChamber->toLocal(h2glopos); - - float dz = h2pos.z()-h1pos.z(); - uz = (h2pos.x()-h1pos.x())/dz ; - vz = (h2pos.y()-h1pos.y())/dz ; - - theChi2 = 0.; - } - else if (nh > 2) { - - // When we have more than two hits then we can fit projections to straight lines - fitSlopes(); - fillChiSquared(); - } // end of 'if' testing no. of hits - - fillLocalDirection(); -} - -void CSCSegAlgoSK::fitSlopes() { - - // Update parameters of fit - // ptc 13-Aug-02: This does a linear least-squares fit - // to the hits associated with the segment, in the z projection. - - // In principle perhaps one could fit the strip and wire - // measurements (u, v respectively), to - // u = u0 + uz * z - // v = v0 + vz * z - // where u0, uz, v0, vz are the parameters resulting from the fit. - // But what is actually done is fit to the local x, y coordinates - // of the RecHits. However the strip measurement controls the precision - // of x, and the wire measurement controls that of y. - // Precision in local coordinate: - // u (strip, sigma~200um), v (wire, sigma~1cm) - - // I have verified that this code agrees with the formulation given - // on p246-247 of 'Data analysis techniques for high-energy physics - // experiments' by Bock, Grote, Notz & Regler, and that on p111-113 - // of 'Statistics' by Barlow. - - // Formulate the matrix equation representing the least-squares fit - // We have a vector of measurements m, which is a 2n x 1 dim matrix - // The transpose mT is (u1, v1, u2, v2, ..., un, vn) - // where ui is the strip-associated measurement and vi is the - // wire-associated measurement for a given RecHit i. - // The fit is to - // u = u0 + uz * z - // v = v0 + vz * z - // where u0, uz, v0, vz are the parameters to be obtained from the fit. - // These are contained in a vector p which is a 4x1 dim matrix, and - // its transpose pT is (u0, v0, uz, vz). Note the ordering! - // The covariance matrix for each pair of measurements is 2 x 2 and - // the inverse of this is the error matrix E. - // The error matrix for the whole set of n measurements is a diagonal - // matrix with diagonal elements the individual 2 x 2 error matrices - // (because the inverse of a diagonal matrix is a diagonal matrix - // with each element the inverse of the original.) - - // The matrix 'matrix' in method 'CSCSegment::weightMatrix()' is this - // block-diagonal overall covariance matrix. It is inverted to the - // block-diagonal error matrix right before it is returned. - - // Use the matrix A defined by - // 1 0 z1 0 - // 0 1 0 z1 - // 1 0 z2 0 - // 0 1 0 z2 - // .. .. .. .. - // 1 0 zn 0 - // 0 1 0 zn - - // The matrix A is returned by 'CSCSegment::derivativeMatrix()'. - - // Then the normal equations are encapsulated in the matrix equation - // - // (AT E A)p = (AT E)m - // - // where AT is the transpose of A. - // We'll call the combined matrix on the LHS, M, and that on the RHS, B: - // M p = B - - // We solve this for the parameter vector, p. - // The elements of M and B then involve sums over the hits - - // The error matrix of the parameters is obtained by - // (AT E A)^-1 calculated in 'calculateError()'. - - // The 4 values in p can be accessed from 'CSCSegment::parameters()' - // in the order uz, vz, u0, v0. - - // NOTE 1 - // Do the #hits = 2 case separately. - // (I hope they're not on the same layer! They should not be, by construction.) - - // NOTE 2 - // We need local position of a RecHit w.r.t. the CHAMBER - // and the RecHit itself only knows its local position w.r.t. - // the LAYER, so we must explicitly transform global position. - - CLHEP::HepMatrix M(4,4,0); - CLHEP::HepVector B(4,0); - - ChamberHitContainer::const_iterator ih = proto_segment.begin(); - - for (ih = proto_segment.begin(); ih != proto_segment.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,1) = IC(1,2); // since Cov is symmetric - IC(2,2) = hit.localPositionError().yy(); - - // 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"; - - //@@ NOW WHAT TO DO? Exception? Return? Ignore? - //@@ Implement throw - } - - // ptc: Note that IC is no longer Cov but Cov^-1 - 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; - } - - // Solve the matrix equation using CLHEP's 'solve' - //@@ ptc: CAN solve FAIL?? UNCLEAR FROM (LACK OF) CLHEP DOC - CLHEP::HepVector p = solve(M, B); - - // Update member variables uz, vz, theOrigin - // Note it has local z = 0 - theOrigin = LocalPoint(p(1), p(2), 0.); - uz = p(3); - vz = p(4); -} - -void CSCSegAlgoSK::fillChiSquared() { - - // The chi-squared is (m-Ap)T E (m-Ap) - // where T denotes transpose. - // This collapses to a simple sum over contributions from each - // pair of measurements. - float u0 = theOrigin.x(); - float v0 = theOrigin.y(); - double chsq = 0.; - - ChamberHitContainer::const_iterator ih; - for (ih = proto_segment.begin(); ih != proto_segment.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); // FIX !! - - double hu = lp.x(); - double hv = lp.y(); - double hz = lp.z(); - - double du = u0 + uz * hz - hu; - double dv = v0 + vz * hz - hv; - - CLHEP::HepMatrix IC(2,2); - IC(1,1) = hit.localPositionError().xx(); - IC(1,2) = hit.localPositionError().xy(); - IC(2,1) = IC(1,2); - IC(2,2) = hit.localPositionError().yy(); - - // Invert covariance matrix - int ierr = 0; - IC.invert(ierr); - if (ierr != 0) { - LogDebug("CSC") << "CSCSegment::fillChiSquared: failed to invert covariance matrix=\n" << IC << "\n"; - - // @@ NOW WHAT TO DO? Exception? Return? Ignore? - } - - chsq += du*du*IC(1,1) + 2.*du*dv*IC(1,2) + dv*dv*IC(2,2); - } - theChi2 = chsq; -} +void CSCSegAlgoSK::updateParameters(void ) { -void CSCSegAlgoSK::fillLocalDirection() { - // Always enforce direction of segment to point from IP outwards - // (Incorrect for particles not coming from IP, of course.) - - double dxdz = uz; - double dydz = vz; - 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( theOrigin ) ).z(); - double globalZdir = ( theChamber->toGlobal( localDir ) ).z(); - double directionSign = globalZpos * globalZdir; - - theDirection = (directionSign * localDir).unit(); - + // Delete input CSCSegFit, create a new one and make the fit + delete sfit_; + sfit_ = new CSCSegFit( theChamber, proto_segment ); + sfit_->fit(); } bool CSCSegAlgoSK::hasHitOnLayer(int layer) const { @@ -693,142 +424,68 @@ bool CSCSegAlgoSK::replaceHit(const CSCRecHit2D* h, int layer) { ++it; } - return addHit(h, layer); + return addHit(h, layer); } void CSCSegAlgoSK::compareProtoSegment(const CSCRecHit2D* h, int layer) { - // compare the chi2 of two segments - double oldChi2 = theChi2; - LocalPoint oldOrigin = theOrigin; - LocalVector oldDirection = theDirection; - ChamberHitContainer oldSegment = proto_segment; + // Copy the input CSCSegFit + CSCSegFit* oldfit = new CSCSegFit( *sfit_ ); + // May create a new fit bool ok = replaceHit(h, layer); if (ok) { - LogDebug("CSC") << " hit in same layer as a hit on segment; try replacing old one..." - << " chi2 new: " << theChi2 << " old: " << oldChi2 << "\n"; + LogDebug("CSCSegment") << " hit in same layer as a hit on segment; try replacing old one..." + << " chi2 new: " << sfit_->chi2() << " old: " << oldfit->chi2() << "\n"; } - if ((theChi2 < oldChi2) && (ok)) { + if ( ( sfit_->chi2() < oldfit->chi2() ) && ok ) { LogDebug("CSC") << " segment with replaced hit is better.\n"; } else { - proto_segment = oldSegment; - theChi2 = oldChi2; - theOrigin = oldOrigin; - theDirection = oldDirection; + // keep original fit + delete sfit_; // now the new fit + sfit_ = oldfit; // reset to the original input fit } } void CSCSegAlgoSK::increaseProtoSegment(const CSCRecHit2D* h, int layer) { - double oldChi2 = theChi2; - LocalPoint oldOrigin = theOrigin; - LocalVector oldDirection = theDirection; - ChamberHitContainer oldSegment = proto_segment; - + // Copy input fit + CSCSegFit* oldfit = new CSCSegFit( *sfit_ ); + + // Creates a new fit bool ok = addHit(h, layer); if (ok) { - LogDebug("CSC") << " hit in new layer: added to segment, new chi2: " - << theChi2 << "\n"; + LogDebug("CSCSegment") << " hit in new layer: added to segment, new chi2: " + << sfit_->chi2() << "\n"; } - int ndf = 2*proto_segment.size() - 4; - - if (ok && ((ndf <= 0) || (theChi2/ndf < chi2Max))) { - LogDebug("CSC") << " segment with added hit is good.\n" ; // FIX + // int ndf = 2*proto_segment.size() - 4; + + //@@ TEST ON ndof<=0 IS JUST TO ACCEPT nhits=2 CASE?? + if ( ok && ( (sfit_->ndof() <= 0) || (sfit_->chi2()/sfit_->ndof() < chi2Max)) ) { + LogDebug("CSCSegment") << " segment with added hit is good.\n" ; } else { - proto_segment = oldSegment; - theChi2 = oldChi2; - theOrigin = oldOrigin; - theDirection = oldDirection; + // reset to original fit + delete sfit_; + sfit_ = oldfit; } } -CLHEP::HepMatrix CSCSegAlgoSK::derivativeMatrix() const { - - ChamberHitContainer::const_iterator it; - int nhits = proto_segment.size(); - CLHEP::HepMatrix matrix(2*nhits, 4); - int row = 0; - - for(it = proto_segment.begin(); it != proto_segment.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; -} - - -AlgebraicSymMatrix CSCSegAlgoSK::weightMatrix() const { - - std::vector::const_iterator it; - int nhits = proto_segment.size(); - AlgebraicSymMatrix matrix(2*nhits, 0); - int row = 0; - - for (it = proto_segment.begin(); it != proto_segment.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; -} - -AlgebraicSymMatrix CSCSegAlgoSK::calculateError() const { - - 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 CSCSegAlgoSK::flipErrors( AlgebraicSymMatrix& a ) const { - - // 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 - +void CSCSegAlgoSK::dumpSegment( const CSCSegment& seg ) const { + + edm::LogVerbatim("CSCSegment") << "CSCSegment in " << theChamber->id() + << "\nlocal position = " << seg.localPosition() + << "\nerror = " << seg.localPositionError() + << "\nlocal direction = " << seg.localDirection() + << "\nerror =" << seg.localDirectionError() + << "\ncovariance matrix" + << seg.parametersError() + << "chi2/ndf = " << seg.chi2() << "/" << seg.degreesOfFreedom() + << "\n#rechits = " << seg.specificRecHits().size() + << "\ntime = " << seg.time(); } diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.h b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.h index 650b8b9f87657..b66d07168389a 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.h +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.h @@ -20,18 +20,18 @@ * Ported to C++ and improved: Rick.Wilkinson@cern.ch
* Reimplemented in terms of layer index, and bug fix: Tim.Cox@cern.ch
* Ported to CMSSW 2006-04-03: Matteo.Sani@cern.ch
+ * Factored out segment fitter Tim.Cox@cern.ch Feb-2015
* - * \author M. Sani */ #include #include -//#include - #include #include +class CSCSegFit; + class CSCSegAlgoSK : public CSCSegmentAlgorithm { public: @@ -92,7 +92,7 @@ class CSCSegAlgoSK : public CSCSegmentAlgorithm { */ void tryAddingHitsToSegment(const ChamberHitContainer& rechitsInChamber, const BoolContainer& used, const LayerIndex& layerIndex, - const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2); + const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2); /** * Return true if segment is 'good'. @@ -108,22 +108,12 @@ class CSCSegAlgoSK : public CSCSegmentAlgorithm { /// Utility functions bool addHit(const CSCRecHit2D* hit, int layer); void updateParameters(void); - void fitSlopes(void); - void fillChiSquared(void); - /** - * Always enforce direction of segment to point from IP outwards - * (Incorrect for particles not coming from IP, of course.) - */ - void fillLocalDirection(void); float phiAtZ(float z) const; bool hasHitOnLayer(int layer) const; bool replaceHit(const CSCRecHit2D* h, int layer); void compareProtoSegment(const CSCRecHit2D* h, int layer); void increaseProtoSegment(const CSCRecHit2D* h, int layer); - CLHEP::HepMatrix derivativeMatrix(void) const; - AlgebraicSymMatrix weightMatrix(void) const; - AlgebraicSymMatrix calculateError(void) const; - void flipErrors(AlgebraicSymMatrix&) const; + void dumpSegment( const CSCSegment& seg ) const; // Member variables // ================ @@ -132,10 +122,6 @@ class CSCSegAlgoSK : public CSCSegmentAlgorithm { ChamberHitContainer proto_segment; const std::string myName; - double theChi2; - LocalPoint theOrigin; - LocalVector theDirection; - float uz, vz; float windowScale; float dRPhiMax ; float dPhiMax; @@ -145,6 +131,8 @@ class CSCSegAlgoSK : public CSCSegmentAlgorithm { float wideSeg; int minLayersApart; bool debugInfo; + + CSCSegFit* sfit_; }; #endif diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc index 1b8315cf28417..63630f10f60f3 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc @@ -1,3 +1,4 @@ + /** * \file CSCSegAlgoST.cc * @@ -5,16 +6,17 @@ * I. Bloch - FNAL * E. James - FNAL * A. Sakharov - WSU + * T. Cox - UC Davis - segment fit factored out of entangled code - Jan 2015 */ - + #include "CSCSegAlgoST.h" +#include "CSCCondSegFit.h" #include "CSCSegAlgoShowering.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include "Geometry/CSCGeometry/interface/CSCLayer.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "CommonTools/Statistics/interface/ChiSquaredProbability.h" @@ -29,7 +31,7 @@ * */ CSCSegAlgoST::CSCSegAlgoST(const edm::ParameterSet& ps) : - CSCSegmentAlgorithm(ps), myName("CSCSegAlgoST"), showering_(0) { + CSCSegmentAlgorithm(ps), myName_("CSCSegAlgoST"), ps_(ps), showering_(0) { debug = ps.getUntrackedParameter("CSCDebug"); // minLayersApart = ps.getParameter("minLayersApart"); @@ -64,24 +66,14 @@ CSCSegAlgoST::CSCSegAlgoST(const edm::ParameterSet& ps) : useShowering = ps.getParameter("useShowering"); if (useShowering) showering_ = new CSCSegAlgoShowering( ps ); - /// Correct the Error Matrix - correctCov_ = ps.getParameter("CorrectTheErrors"); - chi2Norm_2D_ = ps.getParameter("NormChi2Cut2D"); - chi2Norm_3D_ = ps.getParameter("NormChi2Cut3D"); - prePrun_ = ps.getParameter("prePrun"); - prePrunLimit_ = ps.getParameter("prePrunLimit"); - - condSeed1_ = ps.getParameter("SeedSmall"); - condSeed2_ = ps.getParameter("SeedBig"); - covToAnyNumber_ = ps.getParameter("ForceCovariance"); - covToAnyNumberAll_ = ps.getParameter("ForceCovarianceAll"); - covAnyNumber_ = ps.getParameter("Covariance"); - passCondNumber=false; - passCondNumber_2=false; - protoChiUCorrection=1.0; - maxContrIndex=0; - protoNDF = 1.; + /// Improve the covariance matrix? + adjustCovariance_ = ps.getParameter("CorrectTheErrors"); + chi2Norm_3D_ = ps.getParameter("NormChi2Cut3D"); + prePrun_ = ps.getParameter("prePrun"); + prePrunLimit_ = ps.getParameter("prePrunLimit"); + + if (debug) edm::LogVerbatim("CSCSegment") << "CSCSegAlgoST: with factored conditioned segment fit"; } /* Destructor @@ -122,7 +114,7 @@ std::vector CSCSegAlgoST::run(const CSCChamber* aChamber, const Cham a_yweightPenaltyThreshold[3][1] = yweightPenaltyThreshold * 7.60; a_yweightPenaltyThreshold[3][2] = yweightPenaltyThreshold * 20.40; a_yweightPenaltyThreshold[4][1] = yweightPenaltyThreshold * 6.75; - a_yweightPenaltyThreshold[4][2] = yweightPenaltyThreshold * 20.40; //@@ ADD ME4/2 WITH RING 2 VALUE + a_yweightPenaltyThreshold[4][2] = yweightPenaltyThreshold * 20.40; if(preClustering) { // run a pre-clusterer on the given rechits to split clearly-separated segment seeds: @@ -155,7 +147,7 @@ std::vector CSCSegAlgoST::run(const CSCChamber* aChamber, const Cham segments.swap(segments_temp); // segments_temp needed?!?! } - //@@ Ganged strips in ME1/1A? + // Ganged strips in ME1/1A? if ( ("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips() ){ // if ( aChamber->specs()->gangedStrips() ){ findDuplicates(segments); @@ -171,7 +163,7 @@ std::vector CSCSegAlgoST::run(const CSCChamber* aChamber, const Cham segments.swap(segments_temp); // segments_temp needed?!?! } - //@@ Ganged strips in ME1/1A? + // Ganged strips in ME1/1A? if ( ("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips() ){ // if ( aChamber->specs()->gangedStrips() ){ findDuplicates(segments); @@ -342,55 +334,61 @@ std::vector CSCSegAlgoST::prune_bad_hits(const CSCChamber* aChamber, for(size_t m = 0; msetScaleXError( 1.0 ); + segfit->fit(condpass1, condpass2); + // Attempt to handle numerical instability of the fit; // The same as in the build method; // Preprune is not applied; - if(correctCov_){ - if(protoChi2/protoNDF>chi2Norm_3D_){ - passCondNumber = true; - doSlopesAndChi2(); + if( adjustCovariance() ){ + if(segfit->chi2()/segfit->ndof()>chi2Norm_3D_){ + condpass1 = true; + segfit->fit(condpass1, condpass2); } - if((protoChiUCorrection<1.00005)&&(protoChi2/protoNDF>chi2Norm_3D_)){ - passCondNumber_2=true; - doSlopesAndChi2(); + if( (segfit->scaleXError() < 1.00005)&&(segfit->chi2()/segfit->ndof()>chi2Norm_3D_) ){ + condpass2 = true; + segfit->fit(condpass1, condpass2); } } - fillLocalDirection(); - // calculate error matrix - // AlgebraicSymMatrix protoErrors = calculateError(); - SMatrixSym4 protoErrors = calculateError(); - // but reorder components to match what's required by TrackingRecHit interface - // i.e. slopes first, then positions - AlgebraicSymMatrix temp2 = flipErrors( protoErrors ); //@@ TO SATISFY CSCSegment INTERFACE - - CSCSegment temp(protoSegment, protoIntercept, protoDirection, temp2, protoChi2); - - // replace n hit segment with n-1 hit segment, if segment probability is BPMinImprovement better: - if( ( ChiSquaredProbability((double)(*it).chi2(),(double)((2*(*it).nRecHits())-4)) - < - (1./BPMinImprovement)*(ChiSquaredProbability((double)temp.chi2(),(double)(2*temp.nRecHits()-4))) ) // was (1.e-3) 081202 - - && - ( (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)); - // The alternative n-1 segment is much cleaner. If this segment - // has >= minHitsPerSegment hits exchange current n hit segment (*it) - // with better n-1 hit segment: + + // calculate error matrix + // AlgebraicSymMatrix temp2 = segfit->covarianceMatrix(); + + // build an actual segment + CSCSegment temp(protoSegment, segfit->intercept(), segfit->localdir(), + segfit->covarianceMatrix(), segfit->chi2() ); + + // and finished with this fit + delete segfit; + + // n-hit segment is (*it) + // (n-1)-hit segment is temp + // replace n-hit segment with (n-1)-hit segment if segment probability is BPMinImprovement better + double oldchi2 = (*it).chi2(); + double olddof = 2 * (*it).nRecHits() - 4; + double newchi2 = temp.chi2(); + double newdof = 2 * temp.nRecHits() - 4; + if( ( ChiSquaredProbability(oldchi2,olddof) < (1./BPMinImprovement)* + ChiSquaredProbability(newchi2,newdof) ) + && ( ChiSquaredProbability(newchi2,newdof) > best_red_seg_prob ) + && ( ChiSquaredProbability(newchi2,newdof) > 1e-10 ) + ) { + best_red_seg_prob = ChiSquaredProbability(newchi2,newdof); + // The (n-1)- segment is better than the n-hit segment. + // If it has at least minHitsPerSegment hits replace the n-hit segment + // with this better (n-1)-hit segment: if( temp.nRecHits() >= minHitsPerSegment ) { (*it) = temp; } } - } - } + } // end of loop over subsegments (iSegment) + + } // end loop over segments (it) return segments; @@ -535,7 +533,7 @@ std::vector< std::vector > CSCSegAlgoST::chainHits(const CSC seeds.push_back(temp); usedCluster.push_back(false); } - //@@ Only ME1/1A can have ganged strips so no need to test name + // Only ME1/1A can have ganged strips so no need to test name bool gangedME11a = false; if ( ("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips() ){ // if ( aChamber->specs()->gangedStrips() ){ @@ -665,8 +663,8 @@ double CSCSegAlgoST::theWeight(double coordinate_1, double coordinate_2, double } /* - * This algorithm is based on the Minimum Spanning Tree (ST) approach - * for building endcap muon track segments out of the rechit's in a CSCChamber. + * This algorithm uses a Minimum Spanning Tree (ST) approach to build + * endcap muon track segments from the rechits in the 6 layers of a CSC. */ std::vector CSCSegAlgoST::buildSegments(const ChamberHitContainer& rechits) { @@ -814,21 +812,23 @@ std::vector CSCSegAlgoST::buildSegments(const ChamberHitContainer& r //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Showering muon - returns nothing if chi2 == -1 (see comment in SegAlgoShowering) //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - if (useShowering) { - CSCSegment segShower = showering_->showerSeg(theChamber, rechits); - - // Make sure have at least 3 hits... - if ( segShower.nRecHits() < 3 ) return segmentInChamber; - if ( segShower.chi2() == -1 ) return segmentInChamber; - - segmentInChamber.push_back(segShower); - return segmentInChamber; - - } else{ - LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::buildSegments] No. of rechits in the cluster/chamber > " << UpperLimit << - " ... Segment finding in the cluster/chamber canceled!"; - return segmentInChamber; - } + if (useShowering) { + CSCSegment segShower = showering_->showerSeg(theChamber, rechits); + if( debug ) dumpSegment( segShower ); + // Make sure have at least 3 hits... + if ( segShower.nRecHits() < 3 ) return segmentInChamber; + if ( segShower.chi2() == -1 ) return segmentInChamber; + segmentInChamber.push_back(segShower); + return segmentInChamber; + } + else { + // LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::buildSegments] No. of rechits in the cluster/chamber > " + // << UpperLimit << " ... Segment finding in the cluster/chamber canceled!"; + CSCDetId id = rechits[0]->cscDetId(); + edm::LogVerbatim("CSCSegment|CSC") << "[CSCSegAlgoST::buildSegments] No. of rechits in the cluster/chamber > " + << UpperLimit << " ... Segment finding canceled in CSC" << id; + return segmentInChamber; + } } // Find out which station, ring and chamber we are in @@ -1552,35 +1552,38 @@ std::vector CSCSegAlgoST::buildSegments(const ChamberHitContainer& r for(unsigned int iSegment=0; iSegmentchi2Norm_3D_ - // considered as that one potentially suffering from - // numerical instability in fit. - if(correctCov_){ - // Call the fit with prefitting option; - // First fit a straight line to X-Z coordinates - // and calculate chi^2 (chiUZ in correctTheCovX(void)) for X-Z fit; - // Scale up errors in X if chiUZ too big (default 20); - // Refit XY-Z with the scaled up X errors - if(protoChi2/protoNDF>chi2Norm_3D_){ - passCondNumber = true; - doSlopesAndChi2(); + + // Create new fitter object + CSCCondSegFit* segfit = new CSCCondSegFit( pset(), chamber(), protoSegment ); + condpass1 = false; + condpass2 = false; + segfit->setScaleXError( 1.0 ); + segfit->fit(condpass1, condpass2); + + // Attempt to handle numerical instability of the fit. + // (Any segment with chi2/dof > chi2Norm_3D_ is considered + // as potentially suffering from numerical instability in fit.) + if( adjustCovariance() ){ + // Call the fit with prefitting option: + // First fit a straight line to X-Z coordinates and calculate chi2 + // This is done in CSCCondSegFit::correctTheCovX() + // Scale up errors in X if this chi2 is too big (default 'too big' is >20); + // Then refit XY-Z with the scaled-up X errors + if(segfit->chi2()/segfit->ndof()>chi2Norm_3D_){ + condpass1 = true; + segfit->fit(condpass1, condpass2); } - if(protoChiUCorrection<1.00005){ + if(segfit->scaleXError() < 1.00005){ LogTrace("CSCWeirdSegment") << "[CSCSegAlgoST::buildSegments] Segment ErrXX scaled and refit " << std::endl; - if(protoChi2/protoNDF>chi2Norm_3D_){ + if(segfit->chi2()/segfit->ndof()>chi2Norm_3D_){ // Call the fit with direct adjustment of condition number; // If the attempt to improve fit by scaling up X error fails // call the procedure to make the condition number of M compatible with // the precision of X and Y measurements; // Achieved by decreasing abs value of the Covariance LogTrace("CSCWeirdSegment") << "[CSCSegAlgoST::buildSegments] Segment ErrXY changed to match cond. number and refit " << std::endl; - passCondNumber_2=true; - doSlopesAndChi2(); + condpass2 = true; + segfit->fit(condpass1, condpass2); } } // Call the pre-pruning procedure; @@ -1588,27 +1591,33 @@ std::vector CSCSegAlgoST::buildSegments(const ChamberHitContainer& r // while scale factor for X errors is too big. // Prune the recHit inducing the biggest contribution into X-Z chi^2 // and refit; - if(prePrun_ && (sqrt(protoChiUCorrection)>prePrunLimit_) && - (protoSegment.size()>3)){ - LogTrace("CSCWeirdSegment") << "[CSCSegAlgoST::buildSegments] Scale factor protoChiUCorrection too big, pre-Prune and refit " << std::endl; - protoSegment.erase(protoSegment.begin()+(maxContrIndex), - protoSegment.begin()+(maxContrIndex+1)); - doSlopesAndChi2(); + if(prePrun_ && (sqrt(segfit->scaleXError())>prePrunLimit_) && + (segfit->nhits()>3)){ + LogTrace("CSCWeirdSegment") << "[CSCSegAlgoST::buildSegments] Scale factor chi2uCorrection too big, pre-Prune and refit " << std::endl; + protoSegment.erase(protoSegment.begin() + segfit->worstHit(), + protoSegment.begin() + segfit->worstHit()+1 ); + + // Need to create new fitter object to repeat fit with fewer hits + // Original code maintained current values of condpass1, condpass2, scaleXError - calc in CorrectTheCovX() + //@@ DO THE SAME THING HERE, BUT IS THAT CORRECT?! It does make a difference. + double tempcorr = segfit->scaleXError(); // save current value + delete segfit; + segfit = new CSCCondSegFit( pset(), chamber(), protoSegment ); + segfit->setScaleXError( tempcorr ); // reset to previous value (rather than init to 1) + segfit->fit(condpass1, condpass2); } } + + // calculate covariance matrix + // AlgebraicSymMatrix temp2 = segfit->covarianceMatrix(); - fillLocalDirection(); - // calculate error matrix - // AlgebraicSymMatrix protoErrors = calculateError(); - SMatrixSym4 protoErrors = calculateError(); - // but reorder components to match what's required by TrackingRecHit interface - // i.e. slopes first, then positions - AlgebraicSymMatrix temp2 = flipErrors( protoErrors ); //@@ TO SATISFY CSCSegment INTERFACE - - CSCSegment temp(protoSegment, protoIntercept, protoDirection, temp2, protoChi2); - - dumpSegment( temp ); + // build an actual CSC segment + CSCSegment temp(protoSegment, segfit->intercept(), segfit->localdir(), + segfit->covarianceMatrix(), segfit->chi2() ); + delete segfit; + if( debug ) dumpSegment( temp ); + segmentInChamber.push_back(temp); } return segmentInChamber; @@ -1721,456 +1730,8 @@ void CSCSegAlgoST::ChooseSegments2(int best_seg) { } } } -//Method doSlopesAndChi2 -// fitSlopes() and fillChiSquared() are always called one after the other -// In fact the code is duplicated in the two functions (as we need 2 loops) - -// it is much better to fix that at some point -void CSCSegAlgoST::doSlopesAndChi2(){ - fitSlopes(); - fillChiSquared(); -} -/* Method fitSlopes - * - * Perform a Least Square Fit on a segment as per SK algo - * - */ -void CSCSegAlgoST::fitSlopes() { - e_Cxx.clear(); /// Vector of the error matrix (only xx) - if(passCondNumber && !passCondNumber_2){ - correctTheCovX(); - if(e_Cxx.size()!=protoSegment.size()){ - LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::fitSlopes] e_Cxx.size()!=protoSegment.size() ALARM! Please inform CSC DPG " << std::endl; - } - } - - SMatrix4 M; // 4x4, init to 0 - SVector4 B; // 4x1, init to 0; //@@ 4x1 OR 1x4 ?? - - 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); - // Local position of hit w.r.t. chamber - double u = lp.x(); - double v = lp.y(); - double z = lp.z(); - // Covariance matrix of local errors - SMatrixSym2 IC; // 2x2, init to 0 - if(passCondNumber&& !passCondNumber_2){ - IC(0,0) = e_Cxx.at(ih-protoSegment.begin()); - } - else{ - IC(0,0) = hit.localPositionError().xx(); - } - // IC(0,1) = hit.localPositionError().xy(); - IC(1,0) = hit.localPositionError().xy(); - IC(1,1) = hit.localPositionError().yy(); - // IC(1,0) = IC(0,1); // since Cov is symmetric - /// Correct the cov matrix - if(passCondNumber_2){ - correctTheCovMatrix(IC); - } - // Invert covariance matrix (and trap if it fails!) - bool ok = IC.Invert(); - if ( !ok ) { - LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::fitSlopes] Failed to invert covariance matrix: \n" << IC; - // return ok; - } - - M(0,0) += IC(0,0); - M(0,1) += IC(0,1); - M(0,2) += IC(0,0) * z; - M(0,3) += IC(0,1) * z; - B(0) += u * IC(0,0) + v * IC(0,1); - - M(1,0) += IC(1,0); - M(1,1) += IC(1,1); - M(1,2) += IC(1,0) * z; - M(1,3) += IC(1,1) * z; - B(1) += u * IC(1,0) + v * IC(1,1); - - M(2,0) += IC(0,0) * z; - M(2,1) += IC(0,1) * z; - M(2,2) += IC(0,0) * z * z; - M(2,3) += IC(0,1) * z * z; - B(2) += ( u * IC(0,0) + v * IC(0,1) ) * z; - - M(3,0) += IC(1,0) * z; - M(3,1) += IC(1,1) * z; - M(3,2) += IC(1,0) * z * z; - M(3,3) += IC(1,1) * z * z; - B(3) += ( u * IC(1,0) + v * IC(1,1) ) * z; - - } - - SVector4 p; - bool ok = M.Invert(); - if (!ok ){ - LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::fitSlopes[ Failed to invert matrix: \n" << M; - // return ok; - } - else { - p = M * B; - } - - // Update member variables - // Note that origin has local z = 0 - - protoIntercept = LocalPoint(p(0), p(1), 0.); - protoSlope_u = p(2); - protoSlope_v = p(3); - - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::fillSlopes] p = " - // << p(0) << ", " << p(1) << ", " << p(2) << ", " << p(3); - -} -/* Method fillChiSquared - * - * Determine Chi^2 for the proto wire segment - * - */ -void CSCSegAlgoST::fillChiSquared() { - - 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; - - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::fillChiSquared] u, v, z = " << u << ", " << v << ", " << z; - - SMatrixSym2 IC; // 2x2, init to 0 - - if(passCondNumber&& !passCondNumber_2){ - IC(0,0) = e_Cxx.at(ih-protoSegment.begin()); - } - else{ - IC(0,0) = hit.localPositionError().xx(); - } - // IC(0,1) = hit.localPositionError().xy(); - IC(1,0) = hit.localPositionError().xy(); - IC(1,1) = hit.localPositionError().yy(); - // IC(1,0) = IC(0,1); - - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::fillChiSquared] IC before = \n" << IC; - - /// Correct the cov matrix - if(passCondNumber_2){ - correctTheCovMatrix(IC); - } - - // Invert covariance matrix - bool ok = IC.Invert(); - if (!ok ){ - LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::fillChiSquared] Failed to invert covariance matrix: \n" << IC; - // return ok; - } - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::fillChiSquared] IC after = \n" << IC; - chsq += du*du*IC(0,0) + 2.*du*dv*IC(0,1) + dv*dv*IC(1,1); - } - protoChi2 = chsq; - protoNDF = 2.*protoSegment.size() - 4; - - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::fillChiSquared] chi2 = " << protoChi2 << "/" << protoNDF << " dof"; - -} - -/* fillLocalDirection - * - */ -void CSCSegAlgoST::fillLocalDirection() { - // 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 - // 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(); -} - -/* weightMatrix - * - */ -// AlgebraicSymMatrix CSCSegAlgoST::weightMatrix() const { -CSCSegAlgoST::SMatrixSym12 CSCSegAlgoST::weightMatrix() const { - - std::vector::const_iterator it; - - bool ok = true; - - SMatrixSym12 matrix = ROOT::Math::SMatrixIdentity(); // 12x12, init to 1's on diag - - int row = 0; - - for (it = protoSegment.begin(); it != protoSegment.end(); ++it) { - - const CSCRecHit2D& hit = (**it); - - matrix(row, row) = protoChiUCorrection*hit.localPositionError().xx(); - matrix(row, row+1) = hit.localPositionError().xy(); - ++row; - matrix(row, row-1) = hit.localPositionError().xy(); - matrix(row, row) = hit.localPositionError().yy(); - ++row; - } - - ok = matrix.Invert(); // invert in place - if ( !ok ) { - LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::weightMatrix] Failed to invert matrix: \n" << matrix; - // return ok; - } - return matrix; -} - - -/* derivativeMatrix - * - */ -// CLHEP::HepMatrix CSCSegAlgoST::derivativeMatrix() const { -CSCSegAlgoST::SMatrix12by4 CSCSegAlgoST::derivativeMatrix() const { - - ChamberHitContainer::const_iterator it; - - SMatrix12by4 matrix; // 12x4, init to 0 - 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(); - - matrix(row, 0) = 1.; - matrix(row, 2) = z; - ++row; - matrix(row, 1) = 1.; - matrix(row, 3) = z; - ++row; - } - return matrix; -} - - -/* calculateError - * - */ - -CSCSegAlgoST::SMatrixSym4 CSCSegAlgoST::calculateError() const { - - SMatrixSym12 weights = weightMatrix(); - SMatrix12by4 A = derivativeMatrix(); - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::calculateError] weights matrix W: \n" << weights; - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::calculateError] derivatives matrix A: \n" << A; - - // (AT W A)^-1 - // from http://www.phys.ufl.edu/~avery/fitting.html, part I - - bool ok; - SMatrixSym4 result = ROOT::Math::SimilarityT(A, weights); - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::calculateError] (AT W A): \n" << result; - ok = result.Invert(); // inverts in place - if ( !ok ) { - LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::calculateError] Failed to invert matrix: \n" << result; - // return ok; - } - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::calculateError] (AT W A)^-1: \n" << result; - // blithely assuming the inverting never fails... - return result; -} -AlgebraicSymMatrix CSCSegAlgoST::flipErrors( const SMatrixSym4& a ) const { - - // The CSCSegment needs the error matrix re-arranged to match - // parameters in order (uz, vz, u0, v0) where uz, vz = slopes, u0, v0 = intercepts - - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::flipErrors] input: \n" << a; - - AlgebraicSymMatrix hold(4, 0. ); - - for ( short j=0; j!=4; ++j) { - for (short i=0; i!=4; ++i) { - hold(i+1,j+1) = a(i,j); // SMatrix counts from 0, AlgebraicMatrix from 1 - } - } - - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::flipErrors] after copy:"; - // LogTrace("CSCSegAlgoST") << "(" << hold(1,1) << " " << hold(1,2) << " " << hold(1,3) << " " << hold(1,4); - // LogTrace("CSCSegAlgoST") << " " << hold(2,1) << " " << hold(2,2) << " " << hold(2,3) << " " << hold(2,4); - // LogTrace("CSCSegAlgoST") << " " << hold(3,1) << " " << hold(3,2) << " " << hold(3,3) << " " << hold(3,4); - // LogTrace("CSCSegAlgoST") << " " << hold(4,1) << " " << hold(4,2) << " " << hold(4,3) << " " << hold(4,4) << ")"; - - // errors on slopes into upper left - hold(1,1) = a(2,2); - hold(1,2) = a(2,3); - hold(2,1) = a(3,2); - hold(2,2) = a(3,3); - - // errors on positions into lower right - hold(3,3) = a(0,0); - hold(3,4) = a(0,1); - hold(4,3) = a(1,0); - hold(4,4) = a(1,1); - - // must also interchange off-diagonal elements of off-diagonal 2x2 submatrices - hold(4,1) = a(1,2); - hold(3,2) = a(0,3); - hold(2,3) = a(3,0); // = a(0,3) - hold(1,4) = a(2,1); // = a(1,2) - - // LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::flipErrors] after flip:"; - // LogTrace("CSCSegAlgoST") << "(" << hold(1,1) << " " << hold(1,2) << " " << hold(1,3) << " " << hold(1,4); - // LogTrace("CSCSegAlgoST") << " " << hold(2,1) << " " << hold(2,2) << " " << hold(2,3) << " " << hold(2,4); - // LogTrace("CSCSegAlgoST") << " " << hold(3,1) << " " << hold(3,2) << " " << hold(3,3) << " " << hold(3,4); - // LogTrace("CSCSegAlgoST") << " " << hold(4,1) << " " << hold(4,2) << " " << hold(4,3) << " " << hold(4,4) << ")"; - - return hold; -} - -void CSCSegAlgoST::correctTheCovX(void){ - std::vector uu, vv, zz; /// Vectors of coordinates - //std::vector e_Cxx; - e_Cxx.clear(); - double sum_U_err=0.0; - double sum_Z_U_err=0.0; - double sum_Z2_U_err=0.0; - double sum_U_U_err=0.0; - double sum_UZ_U_err=0.0; - std::vector chiUZind; - std::vector::iterator chiContribution; - double chiUZ=0.0; - ChamberHitContainer::const_iterator ih = protoSegment.begin(); - for (ih = protoSegment.begin(); ih != protoSegment.end(); ++ih) { - const CSCRecHit2D& hit = (**ih); - e_Cxx.push_back(hit.localPositionError().xx()); - - const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer()); - GlobalPoint gp = layer->toGlobal(hit.localPosition()); - LocalPoint lp = theChamber->toLocal(gp); - // Local position of hit w.r.t. chamber - double u = lp.x(); - double v = lp.y(); - double z = lp.z(); - uu.push_back(u); - vv.push_back(v); - zz.push_back(z); - // Prepare the sums for the standard linear fit - sum_U_err += 1./e_Cxx.back(); - sum_Z_U_err += z/e_Cxx.back(); - sum_Z2_U_err += (z*z)/e_Cxx.back(); - sum_U_U_err += u/e_Cxx.back(); - sum_UZ_U_err += (u*z)/e_Cxx.back(); - } - - /// Make a one dimensional fit in U-Z plane (i.e. chamber local x-z) - /// U=U0+UZ*Z fit parameters - - double denom=sum_U_err*sum_Z2_U_err-pow(sum_Z_U_err,2); - double U0=(sum_Z2_U_err*sum_U_U_err-sum_Z_U_err*sum_UZ_U_err)/denom; - double UZ=(sum_U_err*sum_UZ_U_err-sum_Z_U_err*sum_U_U_err)/denom; - - /// Calculate the fit line trace - /// Calculate one dimentional chi^2 and normilize the errors if needed - - for(unsigned i=0; i=chi2Norm_2D_){ - protoChiUCorrection = chiUZ/chi2Norm_2D_; - for(unsigned i=0; iprePrunLimit_){ - chiContribution=max_element(chiUZind.begin(),chiUZind.end()); - maxContrIndex = chiContribution - chiUZind.begin(); - /* - for(unsigned i=0; icondNumberCorr2))|| - ((diag2>condNumberCorr2)&&(detCov & segments ){ // this is intended for ME1/1a only - we have ghost segments because of the strips ganging // this function finds them (first the rechits by sharesInput() ) @@ -2199,7 +1760,10 @@ void CSCSegAlgoST::findDuplicates(std::vector & segments ){ } void CSCSegAlgoST::dumpSegment( const CSCSegment& seg ) const { - LogTrace("CSCSegment") << "CSCSegment in " << chamber()->id() + + // Only called if pset value 'CSCDebug' is set in config + + edm::LogVerbatim("CSCSegment") << "CSCSegment in " << chamber()->id() << "\nlocal position = " << seg.localPosition() << "\nerror = " << seg.localPositionError() << "\nlocal direction = " << seg.localDirection() @@ -2207,7 +1771,7 @@ void CSCSegAlgoST::dumpSegment( const CSCSegment& seg ) const { << "\ncovariance matrix" << seg.parametersError() << "chi2/ndf = " << seg.chi2() << "/" << seg.degreesOfFreedom() - << "\n#rechits = " << seg.specificRecHits().size(); + << "\n#rechits = " << seg.specificRecHits().size() + << "\ntime = " << seg.time(); } - diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.h b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.h index 166340936c775..a4b6adba2b02b 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.h +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.h @@ -10,27 +10,21 @@ * A CSCSegment is a RecSegment4D, and is built from * CSCRecHit2D objects, each of which is a RecHit2DLocalPos.
* - * This builds segments consisting of at least 3 hits. It is allowed for segments to have - * a common (only one) rechit. + * This builds segments consisting of at least 3 hits. + * Segments can share a common rechit, but only one. * - * The program is under construction/testing, as it has been for many years?! - * As of Aug-2014 replace CLHEP matrices by ROOT::SMatrix in a minimal way. - * * \authors S. Stoynev - NWU * I. Bloch - FNAL * E. James - FNAL * A. Sakharov - WSU (extensive revision to handle weird segments) + * ... ... ... * T. Cox - UC Davis (struggling to handle this monster) * */ #include #include - -#include -#include -#include - +#include #include #include @@ -46,25 +40,9 @@ class CSCSegAlgoST : public CSCSegmentAlgorithm { typedef std::vector < std::vector > Segments; typedef std::deque BoolContainer; - // 12 x12 Symmetric - typedef ROOT::Math::SMatrix > SMatrixSym12; - - // 12 x 4 - typedef ROOT::Math::SMatrix SMatrix12by4; - - // 4 x 4 General + Symmetric - typedef ROOT::Math::SMatrix SMatrix4; - typedef ROOT::Math::SMatrix > SMatrixSym4; - - // 2 x 2 Symmetric - typedef ROOT::Math::SMatrix > SMatrixSym2; - - // 4-dim vector - typedef ROOT::Math::SVector SVector4; - - /// Constructor explicit CSCSegAlgoST(const edm::ParameterSet& ps); + /// Destructor virtual ~CSCSegAlgoST(); @@ -104,6 +82,12 @@ class CSCSegAlgoST : public CSCSegmentAlgorithm { private: + // Retrieve pset + const edm::ParameterSet& pset(void) const { return ps_;} + + // Adjust covariance matrix? + bool adjustCovariance(void) { return adjustCovariance_;} + /// Utility functions double theWeight(double coordinate_1, double coordinate_2, double coordinate_3, float layer_1, float layer_2, float layer_3); @@ -118,31 +102,20 @@ class CSCSegAlgoST : public CSCSegmentAlgorithm { void ChooseSegments3(int best_seg); void ChooseSegments3(std::vector< ChamberHitContainer > & best_segments, std::vector< float > & best_weight, int best_seg); // - void fitSlopes(void); - void fillChiSquared(void); - void fillLocalDirection(void); - void doSlopesAndChi2(void); // Find duplicates in ME1/1a, if it has ganged strips (i.e. pre-LS1) void findDuplicates(std::vector & segments ); bool isGoodToMerge(bool isME11a, ChamberHitContainer & newChain, ChamberHitContainer & oldChain); - // THE FOLLOWING BLITHELY PASS MATRICES AROUND LIKE CONFETTI - // (And that must be fixed!) - - SMatrix12by4 derivativeMatrix(void) const; - SMatrixSym12 weightMatrix(void) const; - SMatrixSym4 calculateError(void) const; - AlgebraicSymMatrix flipErrors(const SMatrixSym4&) const; - void correctTheCovMatrix(SMatrixSym2& IC); - void correctTheCovX(void); - void dumpSegment( const CSCSegment& seg ) const; const CSCChamber* chamber() const {return theChamber;} // Member variables - const std::string myName; + const std::string myName_; + const edm::ParameterSet ps_; + CSCSegAlgoShowering* showering_; + const CSCChamber* theChamber; Segments GoodSegments; @@ -182,15 +155,7 @@ class CSCSegAlgoST : public CSCSegmentAlgorithm { std::vector< float > weight_noL5_B; std::vector< float > weight_noL6_B; - //ibl - ChamberHitContainer protoSegment; - float protoSlope_u; - float protoSlope_v; - LocalPoint protoIntercept; - double protoChi2; - double protoNDF; - LocalVector protoDirection; // input from .cfi file bool debug; @@ -221,26 +186,18 @@ class CSCSegAlgoST : public CSCSegmentAlgorithm { double curvePenaltyThreshold; double curvePenalty; - CSCSegAlgoShowering* showering_; - /// Correct the Error Matrix - bool correctCov_; /// Allow to correct the error matrix - double protoChiUCorrection; - std::vector e_Cxx; - double chi2Norm_2D_; /// Chi^2 normalization for the corrected fit + + bool adjustCovariance_; /// Flag whether to 'improve' covariance matrix + + bool condpass1, condpass2; + double chi2Norm_3D_; /// Chi^2 normalization for the initial fit - unsigned maxContrIndex; /// The index of the worst x RecHit in Chi^2-X method + bool prePrun_; /// Allow to prune a (rechit in a) segment in segment buld method - /// once it passed through Chi^2-X and protoChiUCorrection is big. + /// once it passed through Chi^2-X and chi2uCorrection is big. double prePrunLimit_; /// The upper limit of protoChiUCorrection to apply prePrun - /// Correct the error matrix for the condition number - double condSeed1_, condSeed2_; /// The correction parameters - bool covToAnyNumber_; /// Allow to use any number for covariance (by hand) - bool covToAnyNumberAll_; /// Allow to use any number for covariance for all RecHits - double covAnyNumber_; /// The number to force the Covariance - bool passCondNumber; /// Passage the condition number calculations - bool passCondNumber_2; /// Passage the condition number calculations }; #endif diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.cc index 18ab9e4c1e1c9..3eb1f29748d20 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.cc @@ -1,13 +1,12 @@ /** * \file CSCSegAlgoShowering.cc * - * \author: D. Fortin - UC Riverside - * \modified: J. Babb - UC Riverside - * - * See header file for description. + * Last update: 17.02.2015 + * */ #include "RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.h" +#include "RecoLocalMuon/CSCSegment/src/CSCSegFit.h" #include "Geometry/CSCGeometry/interface/CSCLayer.h" #include @@ -25,7 +24,7 @@ /* Constructor * */ -CSCSegAlgoShowering::CSCSegAlgoShowering(const edm::ParameterSet& ps) { +CSCSegAlgoShowering::CSCSegAlgoShowering(const edm::ParameterSet& ps) : sfit_(0) { // debug = ps.getUntrackedParameter("CSCSegmentDebug"); dRPhiFineMax = ps.getParameter("dRPhiFineMax"); dPhiFineMax = ps.getParameter("dPhiFineMax"); @@ -65,7 +64,7 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha } // Loop over hits to find center-of-mass position in each layer - for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); it++ ) { + for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); ++it ) { const CSCRecHit2D& hit = (**it); const CSCDetId id = hit.cscDetId(); int l_id = id.layer(); @@ -73,7 +72,7 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha GlobalPoint gp = layer->toGlobal(hit.localPosition()); LocalPoint lp = theChamber->toLocal(gp); - n[l_id -1]++; + ++n[l_id -1]; x[l_id -1] += lp.x(); y[l_id -1] += lp.y(); gz[l_id -1] += gp.z(); @@ -115,9 +114,9 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha // by projecting the vector std::vector layerPoints; - for (unsigned i = 1; i < 7; i++) { + for (size_t i = 0; i!=6; ++i) { // Get the layer z coordinates in global frame - const CSCLayer* layer = theChamber->layer(i); + const CSCLayer* layer = theChamber->layer(i+1); LocalPoint temp(0., 0., 0.); GlobalPoint gp = layer->toGlobal(temp); float layer_Z = gp.z(); @@ -138,7 +137,7 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha std::vector r_closest; std::vector id; - for (unsigned i = 0; i < 6; ++i ) { + for (size_t i = 0; i!=6; ++i ) { id.push_back(-1); r_closest.push_back(9999.); } @@ -146,10 +145,11 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha int idx = 0; // Loop over all hits and find hit closest to com for that layer. - for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); it++ ) { + for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); ++it ) { const CSCRecHit2D& hit = (**it); int layId = hit.cscDetId().layer(); - const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer()); + + const CSCLayer* layer = theChamber->layer(layId); GlobalPoint gp = layer->toGlobal(hit.localPosition()); LocalPoint lp = theChamber->toLocal(gp); @@ -157,12 +157,12 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha float d_y = lp.y() - layerPoints[layId-1].y(); LocalPoint diff(d_x, d_y, 0.); - + if ( fabs(diff.mag() ) < r_closest[layId-1] ) { r_closest[layId-1] = fabs(diff.mag()); id[layId-1] = idx; } - idx++; + ++idx; } // Now fill vector of rechits closest to center of mass: @@ -170,13 +170,13 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha idx = 0; // Loop over all hits and find hit closest to com for that layer. - for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); it++ ) { + for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); ++it ) { const CSCRecHit2D& hit = (**it); int layId = hit.cscDetId().layer(); if ( idx == id[layId-1] )protoSegment.push_back(*it); - idx++; + ++idx; } // Reorder hits in protosegment @@ -191,42 +191,26 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha } } - - // Compute the segment properties + // Fit the protosegment updateParameters(); - // Clean up protosegment if there is one very bad hit on segment + // If there is one very bad hit on segment, remove it and refit if (protoSegment.size() > 4) pruneFromResidual(); - // Look for better hits near segment + // If any hit on a layer is closer to segment than original, replace it and refit for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); it++ ) { - const CSCRecHit2D* h = *it; int layer = h->cscDetId().layer(); - - if ( isHitNearSegment( h ) ) compareProtoSegment(h, layer); + if ( isHitNearSegment( h ) ) compareProtoSegment( h, layer ); } - // Prune worse hit if necessary - if ( protoSegment.size() > 5 ) pruneFromResidual(); + // Check again for a bad hit, and remove and refit if necessary + if ( sfit_->nhits() > 5 ) pruneFromResidual( ); - // Update the parameters - updateParameters(); + // Does the fitted line point to the IP? + // If it doesn't, the algorithm has probably failed i.e. that's life! - // 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; - LocalVector protoDirection = (directionSign * localDir).unit(); - - // Check to make sure the fitting didn't fuck things up - GlobalVector protoGlobalDir = theChamber->toGlobal( protoDirection ); + GlobalVector protoGlobalDir = theChamber->toGlobal( sfit_->localdir() ); double protoTheta = protoGlobalDir.theta(); double protoPhi = protoGlobalDir.phi(); double simTheta = gvCOM.theta(); @@ -236,22 +220,24 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha float dPhi = fabs(protoPhi - simPhi); // float dR = sqrt(dEta*dEta + dPhi*dPhi); - // Error matrix - AlgebraicSymMatrix protoErrors = calculateError(); - // but reorder components to match what's required by TrackingRecHit interface - // i.e. slopes first, then positions - flipErrors( protoErrors ); - - //Flag the segment with chi12 of -1 of the segment isn't pointing toward origin + // Flag the segment with chi2=-1 of the segment isn't pointing toward origin + // i.e. flag that the algorithm has probably just failed (I presume we expect + // a segment to point to the IP if the algorithm is successful!) + + double theFlag = -1.; if (dTheta > maxDTheta || dPhi > maxDPhi) { - CSCSegment temp(protoSegment, protoIntercept, protoDirection, protoErrors, -1); - return temp; } else { - CSCSegment temp(protoSegment, protoIntercept, protoDirection, protoErrors, protoChi2); - return temp; + theFlag = sfit_->chi2(); // If it points to IP, just pass fit chi2 as usual } + // Create an actual CSCSegment - retrieve all info from the fit + CSCSegment temp(sfit_->hits(), sfit_->intercept(), + sfit_->localdir(), sfit_->covarianceMatrix(), theFlag ); + delete sfit_; + sfit_ = 0; + + return temp; } @@ -261,7 +247,7 @@ CSCSegment CSCSegAlgoShowering::showerSeg( const CSCChamber* aChamber, const Cha * * Compare rechit with expected position from proto_segment */ -bool CSCSegAlgoShowering::isHitNearSegment( const CSCRecHit2D* hit) const { +bool CSCSegAlgoShowering::isHitNearSegment( const CSCRecHit2D* hit ) const { const CSCLayer* layer = theChamber->layer(hit->cscDetId().layer()); @@ -272,8 +258,8 @@ bool CSCSegAlgoShowering::isHitNearSegment( const CSCRecHit2D* hit) const { LocalPoint Hlp = theChamber->toLocal(Hgp); double z = Hlp.z(); - double LocalX = protoIntercept.x() + protoSlope_u * z; - double LocalY = protoIntercept.y() + protoSlope_v * z; + double LocalX = sfit_->xfit(z); + double LocalY = sfit_->yfit(z); LocalPoint Slp(LocalX, LocalY, z); GlobalPoint Sgp = theChamber->toGlobal(Slp); double Sphi = Sgp.phi(); @@ -300,8 +286,8 @@ bool CSCSegAlgoShowering::isHitNearSegment( const CSCRecHit2D* hit) const { */ bool CSCSegAlgoShowering::addHit(const CSCRecHit2D* aHit, int layer) { - // 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. + // Return true if hit was added successfully + // Return false if there is already a hit on the same layer bool ok = true; @@ -315,115 +301,6 @@ bool CSCSegAlgoShowering::addHit(const CSCRecHit2D* aHit, int layer) { } -/* Method updateParameters - * - * Perform a simple Least Square Fit on proto segment to determine slope and intercept - * - */ -void CSCSegAlgoShowering::updateParameters() { - - // 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; -} - - - /* Method compareProtoSegment * * For hit coming from the same layer of an existing hit within the proto segment @@ -432,16 +309,7 @@ void CSCSegAlgoShowering::updateParameters() { */ void CSCSegAlgoShowering::compareProtoSegment(const CSCRecHit2D* h, int layer) { - // 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 + // Try adding the hit to existing segment, and removing one from the same layer ChamberHitContainer::iterator it; for ( it = protoSegment.begin(); it != protoSegment.end(); ) { if ( (*it)->cscDetId().layer() == layer ) { @@ -451,123 +319,39 @@ void CSCSegAlgoShowering::compareProtoSegment(const CSCRecHit2D* h, int layer) { } } 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; - } -} - -AlgebraicSymMatrix CSCSegAlgoShowering::weightMatrix() const { - - std::vector::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; -} - - -CLHEP::HepMatrix CSCSegAlgoShowering::derivativeMatrix() const { - - 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; + if (ok) { + CSCSegFit* newfit = new CSCSegFit(theChamber, protoSegment); + newfit->fit(); + if ( newfit->chi2() > sfit_->chi2() ) { + // new fit is worse: revert to old fit + delete newfit; + } + else { + // new fit is better + delete sfit_; + sfit_ = newfit; + } } - return matrix; } -/* calculateError - * - */ -AlgebraicSymMatrix CSCSegAlgoShowering::calculateError() const { - AlgebraicSymMatrix weights = weightMatrix(); - AlgebraicMatrix A = derivativeMatrix(); +// Look for a hit with a large deviation from fit - // (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); +void CSCSegAlgoShowering::pruneFromResidual(void){ - // blithely assuming the inverting never fails... - return result; -} - -void CSCSegAlgoShowering::flipErrors( AlgebraicSymMatrix& a ) const { - - // 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 - -} - -// Try to clean up segments by quickly looking at residuals -void CSCSegAlgoShowering::pruneFromResidual(){ + //@@ THIS FUNCTION HAS 3 RETURNS PATHS! // Only prune if have at least 5 hits - if ( protoSegment.size() < 5 ) return ; - + 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; + ChamberHitContainer::iterator ih; + ChamberHitContainer::iterator ibad = protoSegment.end(); for ( ih = protoSegment.begin(); ih != protoSegment.end(); ++ih ) { const CSCRecHit2D& hit = (**ih); @@ -579,18 +363,17 @@ void CSCSegAlgoShowering::pruneFromResidual(){ double v = lp.y(); double z = lp.z(); - double du = protoIntercept.x() + protoSlope_u * z - u; - double dv = protoIntercept.y() + protoSlope_v * z - v; + // double du = segfit->xdev( u, z ); + // double dv = segfit->ydev( v, z ); - float residual = sqrt(du*du + dv*dv); + float residual = sfit_->Rdev(u, v, z); // == sqrt(du*du + dv*dv) sumResidual += residual; - nHits++; + ++nHits; if ( residual > maxResidual ) { maxResidual = residual; - badIndex = j; + ibad = ih; } - j++; } float corrAvgResidual = (sumResidual - maxResidual)/(nHits -1); @@ -598,28 +381,23 @@ void CSCSegAlgoShowering::pruneFromResidual(){ // Keep all hits if ( maxResidual/corrAvgResidual < maxRatioResidual ) return; + // Drop worst hit + if( ibad != protoSegment.end() ) protoSegment.erase(ibad); - // 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 + // Make a new fit updateParameters(); + return; } +void CSCSegAlgoShowering::updateParameters(void) { + // Create fit for the hits in the protosegment & run it + delete sfit_; + sfit_ = new CSCSegFit( theChamber, protoSegment ); + sfit_->fit(); +} + + diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.h b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.h index 394ad9c6b4adc..b105b3e00d125 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.h +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.h @@ -6,6 +6,7 @@ * * \author: D. Fortin - UC Riverside * \modified by J. Babb - UC Riverside + * Updated Tim Cox - UC Davis Feb-2015 - factored out segment fit * * Handle case where too many hits are reconstructed in the chamber, even after preclustering * for normal segment reconstruction to properly handle these. @@ -18,9 +19,10 @@ #include #include #include - #include +class CSCSegFit; + class CSCSegAlgoShowering { public: @@ -41,26 +43,16 @@ class CSCSegAlgoShowering { /// Utility functions bool isHitNearSegment(const CSCRecHit2D* h) const; bool addHit(const CSCRecHit2D* hit, int layer); - void updateParameters(void); bool hasHitOnLayer(int layer) const; void compareProtoSegment(const CSCRecHit2D* h, int layer); - CLHEP::HepMatrix derivativeMatrix(void) const; - AlgebraicSymMatrix weightMatrix(void) const; - AlgebraicSymMatrix calculateError(void) const; - void flipErrors(AlgebraicSymMatrix&) const; - - void pruneFromResidual(); + void pruneFromResidual(void); + void updateParameters(void); // Member variables const std::string myName; const CSCChamber* theChamber; ChamberHitContainer protoSegment; - float protoSlope_u; - float protoSlope_v; - LocalPoint protoIntercept; - double protoChi2; - LocalVector protoDirection; // input from .cfi file bool debug; @@ -75,6 +67,7 @@ class CSCSegAlgoShowering { float maxDTheta; float maxDPhi; + CSCSegFit* sfit_; // current fit }; #endif diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc index e5c75daf3d291..da74bd2c2b3b4 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc @@ -1,17 +1,16 @@ /** * \file CSCSegAlgoTC.cc * - * \author M. Sani - * + * Last update: 13.02.2015 + * */ #include "CSCSegAlgoTC.h" +#include "CSCSegFit.h" #include "DataFormats/CSCRecHit/interface/CSCSegment.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" -// For clhep Matrix::solve -#include "DataFormats/CLHEP/interface/AlgebraicObjects.h" #include "Geometry/CSCGeometry/interface/CSCLayer.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -24,8 +23,8 @@ #include #include -CSCSegAlgoTC::CSCSegAlgoTC(const edm::ParameterSet& ps) : CSCSegmentAlgorithm(ps), - myName("CSCSegAlgoTC") { +CSCSegAlgoTC::CSCSegAlgoTC(const edm::ParameterSet& ps) + : CSCSegmentAlgorithm(ps), sfit_(0), myName("CSCSegAlgoTC") { debugInfo = ps.getUntrackedParameter("verboseInfo"); @@ -52,22 +51,31 @@ CSCSegAlgoTC::CSCSegAlgoTC(const edm::ParameterSet& ps) : CSCSegmentAlgorithm(ps std::vector CSCSegAlgoTC::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits) { theChamber = aChamber; + return buildSegments(rechits); } -std::vector CSCSegAlgoTC::buildSegments(const ChamberHitContainer& _rechits) { +std::vector CSCSegAlgoTC::buildSegments(const ChamberHitContainer& urechits) { - // Reimplementation of original algorithm of CSCSegmentizer, Mar-06 + // ChamberHitContainer rechits = urechits; + ChamberHitContainer rechits(urechits); + + // edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoTC::buildSegments] start building segments in " << theChamber->id(); + // edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoTC::buildSegments] size of rechit container = " << rechits.size(); + + if (rechits.size() < 2) { + LogDebug("CSC") << myName << ": " << rechits.size() << + " hit(s) in chamber is not enough to build a segment.\n"; + return std::vector(); + } - LogDebug("CSC") << "*********************************************"; - LogDebug("CSC") << "Start segment building in the new chamber: " << theChamber->specs()->chamberTypeName(); - LogDebug("CSC") << "*********************************************"; - ChamberHitContainer rechits = _rechits; LayerIndex layerIndex(rechits.size()); - - for(unsigned int i = 0; i < rechits.size(); i++) { - - layerIndex[i] = rechits[i]->cscDetId().layer(); + + for ( size_t i = 0; i < rechits.size(); ++i ) { + short ilay = rechits[i]->cscDetId().layer(); + layerIndex[i] = ilay; + // edm::LogVerbatim("CSCSegment") << "layerIndex[" << i << "] should be " << rechits[i]->cscDetId().layer(); + // edm::LogVerbatim("CSCSegment") << "layerIndex[" << i << "] actually is " << layerIndex[i]; } double z1 = theChamber->layer(1)->position().z(); @@ -86,16 +94,14 @@ std::vector CSCSegAlgoTC::buildSegments(const ChamberHitContainer& _ } } - if (debugInfo) { - // dump after sorting - dumpHits(rechits); - } + // Dump rechits after sorting? + // if (debugInfo) dumpHits(rechits); - if (rechits.size() < 2) { - LogDebug("CSC") << myName << ": " << rechits.size() << - " hit(s) in chamber is not enough to build a segment.\n"; - return std::vector(); - } + // if (rechits.size() < 2) { + // LogDebug("CSC") << myName << ": " << rechits.size() << + // " hit(s) in chamber is not enough to build a segment.\n"; + // return std::vector(); + // } // We have at least 2 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. @@ -115,6 +121,8 @@ std::vector CSCSegAlgoTC::buildSegments(const ChamberHitContainer& _ std::vector segments; + + sfit_ = 0; ChamberHitContainerCIt ib = rechits.begin(); ChamberHitContainerCIt ie = rechits.end(); @@ -122,6 +130,7 @@ std::vector CSCSegAlgoTC::buildSegments(const ChamberHitContainer& _ for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) { int layer1 = layerIndex[i1-ib]; + const CSCRecHit2D* h1 = *i1; for (ChamberHitContainerCIt i2 = ie-1; i2 != i1; --i2) { @@ -141,67 +150,69 @@ std::vector CSCSegAlgoTC::buildSegments(const ChamberHitContainer& _ GlobalPoint gp1 = l1->toGlobal(h1->localPosition()); const CSCLayer* l2 = theChamber->layer(h2->cscDetId().layer()); GlobalPoint gp2 = l2->toGlobal(h2->localPosition()); - LogDebug("CSC") << "start new segment from hits " << "h1: " << gp1 << " - h2: " << gp2 << "\n"; + LogDebug("CSCSegment") << "start new segment from hits " << "h1: " << gp1 << " - h2: " << gp2 << "\n"; + // edm::LogVerbatim("CSCSegment") << "start new segment from hits " << "h1: " << gp1 << " - h2: " << gp2; + // edm::LogVerbatim("CSCSegment") << "on layers " << layer1 << " and " << layer2; - if (!addHit(h1, layer1)) { - LogDebug("CSC") << " failed to add hit h1\n"; + if ( !addHit(h1, layer1) ) { + LogDebug("CSCSegment") << " failed to add hit h1\n"; continue; } - if (!addHit(h2, layer2)) { - LogDebug("CSC") << " failed to add hit h2\n"; + if ( !addHit(h2, layer2) ) { + LogDebug("CSCSegment") << " failed to add hit h2\n"; continue; } - tryAddingHitsToSegment(rechits, i1, i2); // changed seg + if ( sfit_ ) tryAddingHitsToSegment(rechits, i1, i2); // can only add hits if there's a segment // if a segment has been found push back it into the segment collection if (proto_segment.empty()) { - - LogDebug("CSC") << "No segment has been found !!!\n"; + LogDebug("CSCSegment") << "No segment found.\n"; + // edm::LogVerbatim("CSCSegment") << "No segment found."; } else { - - // calculate error matrix - AlgebraicSymMatrix error_matrix = calculateError(); - - // but reorder components to match what's required by TrackingRecHit interface - // i.e. slopes first, then positions - flipErrors( error_matrix ); - - candidates.push_back(proto_segment); - origins.push_back(theOrigin); - directions.push_back(theDirection); - errors.push_back(error_matrix); - chi2s.push_back(theChi2); - LogDebug("CSC") << "Found a segment !!!\n"; + //@@ THIS IS GOING TO BE TRICKY - CREATING MANY FITS ON HEAP + //@@ BUT MEMBER VARIABLE JUST POINTS TO MOST RECENT ONE + candidates.push_back( sfit_ ); // store the current fit + LogDebug("CSCSegment") << "Found a segment.\n"; + // edm::LogVerbatim("CSCSegment") << "Found a segment."; } } } } + + // edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoTC::buildSegments] no. of candidates before pruning = " << candidates.size(); // We've built all possible segments. Now pick the best, non-overlapping subset. pruneTheSegments(rechits); - // Copy the selected proto segments into the CSCSegment vector - for(unsigned int i=0; i < candidates.size(); i++) { - - CSCSegment temp(candidates[i], origins[i], directions[i], errors[i], chi2s[i]); + // Create CSCSegments for the surviving candidates + for(unsigned int i = 0; i < candidates.size(); ++i ) { + CSCSegFit*sfit = candidates[i]; + // edm::LogVerbatim("CSCSegment") << "candidate fit " << i+1 << " of " << candidates.size() << " is at " << sfit; + // if ( !sfit ) { + // edm::LogVerbatim("CSCSegment") << "stored a null pointer for element " << i+1 << " of " << candidates.size(); + // continue; + // } + CSCSegment temp(sfit->hits(), sfit->intercept(), sfit->localdir(), + sfit->covarianceMatrix(), sfit->chi2() ); + delete sfit; segments.push_back(temp); + if (debugInfo) dumpSegment( temp ); } - + + // reset member variables candidates.clear(); - origins.clear(); - directions.clear(); - errors.clear(); - chi2s.clear(); - - // Give the segments to the CSCChamber + sfit_ = 0; + + // edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoTC::buildSegments] no. of segments returned = " << segments.size(); + return segments; } void CSCSegAlgoTC::tryAddingHitsToSegment(const ChamberHitContainer& rechits, - const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2) { + const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2) { // Iterate over the layers with hits in the chamber // Skip the layers containing the segment endpoints @@ -212,7 +223,7 @@ void CSCSegAlgoTC::tryAddingHitsToSegment(const ChamberHitContainer& rechits, // then replace the original segment with the new one (by swap) // - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory // then replace the original segment with the new one (by swap) - + ChamberHitContainerCIt ib = rechits.begin(); ChamberHitContainerCIt ie = rechits.end(); @@ -250,6 +261,7 @@ bool CSCSegAlgoTC::addHit(const CSCRecHit2D* aHit, int layer) { // 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; ChamberHitContainer::const_iterator it; @@ -265,7 +277,7 @@ bool CSCSegAlgoTC::addHit(const CSCRecHit2D* aHit, int layer) { } bool CSCSegAlgoTC::replaceHit(const CSCRecHit2D* h, int layer) { - + // replace a hit from a layer ChamberHitContainer::iterator it; for (it = proto_segment.begin(); it != proto_segment.end();) { @@ -275,295 +287,25 @@ bool CSCSegAlgoTC::replaceHit(const CSCRecHit2D* h, int layer) { ++it; } - return addHit(h, layer); -} - -void CSCSegAlgoTC::updateParameters() { - - // Note that we need local position of a RecHit w.r.t. the CHAMBER - // and the RecHit itself only knows its local position w.r.t. - // the LAYER, so need to explicitly transform to global position. - - // no. of hits in the RecHitsOnSegment - // By construction this is the no. of layers with hitsna parte รจ da inserirsi tra le Contrade aperte ad accettare quello che - - - // since we allow just one hit per layer in a segment. - - int nh = proto_segment.size(); - - // First hit added to a segment must always fail here - if (nh < 2) - return; - - if (nh == 2) { - - // Once we have two hits we can calculate a straight line - // (or rather, a straight line for each projection in xz and yz.) - ChamberHitContainer::const_iterator ih = proto_segment.begin(); - int il1 = (*ih)->cscDetId().layer(); - const CSCRecHit2D& h1 = (**ih); - ++ih; - int il2 = (*ih)->cscDetId().layer(); - const CSCRecHit2D& h2 = (**ih); - - //@@ Skip if on same layer, but should be impossible - if (il1 == il2) - return; - - const CSCLayer* layer1 = theChamber->layer(il1); - const CSCLayer* layer2 = theChamber->layer(il2); - - GlobalPoint h1glopos = layer1->toGlobal(h1.localPosition()); - GlobalPoint h2glopos = layer2->toGlobal(h2.localPosition()); - - // localPosition is position of hit wrt layer (so local z = 0) - theOrigin = h1.localPosition(); - - // We want hit wrt chamber (and local z will be != 0) - LocalPoint h1pos = theChamber->toLocal(h1glopos); // FIX !! - LocalPoint h2pos = theChamber->toLocal(h2glopos); // FIX !! - - float dz = h2pos.z()-h1pos.z(); - uz = (h2pos.x()-h1pos.x())/dz ; - vz = (h2pos.y()-h1pos.y())/dz ; - - theChi2 = 0.; - } - else if (nh > 2) { - - // When we have more than two hits then we can fit projections to straight lines - fitSlopes(); - fillChiSquared(); - } // end of 'if' testing no. of hits - - fillLocalDirection(); -} + return addHit(h, layer); -void CSCSegAlgoTC::fitSlopes() { - - // Update parameters of fit - // ptc 13-Aug-02: This does a linear least-squares fit - // to the hits associated with the segment, in the z projection. - - // In principle perhaps one could fit the strip and wire - // measurements (u, v respectively), to - // u = u0 + uz * z - // v = v0 + vz * z - // where u0, uz, v0, vz are the parameters resulting from the fit. - // But what is actually done is fit to the local x, y coordinates - // of the RecHits. However the strip measurement controls the precision - // of x, and the wire measurement controls that of y. - // Precision in local coordinate: - // u (strip, sigma~200um), v (wire, sigma~1cm) - - // I have verified that this code agrees with the formulation given - // on p246-247 of 'Data analysis techniques for high-energy physics - // experiments' by Bock, Grote, Notz & Regler, and that on p111-113 - // of 'Statistics' by Barlow. - - // Formulate the matrix equation representing the least-squares fit - // We have a vector of measurements m, which is a 2n x 1 dim matrix - // The transpose mT is (u1, v1, u2, v2, ..., un, vn) - // where ui is the strip-associated measurement and vi is the - // wire-associated measurement for a given RecHit i. - // The fit is to - // u = u0 + uz * z - // v = v0 + vz * z - // where u0, uz, v0, vz are the parameters to be obtained from the fit. - // These are contained in a vector p which is a 4x1 dim matrix, and - // its transpose pT is (u0, v0, uz, vz). Note the ordering! - // The covariance matrix for each pair of measurements is 2 x 2 and - // the inverse of this is the error matrix E. - // The error matrix for the whole set of n measurements is a diagonal - // matrix with diagonal elements the individual 2 x 2 error matrices - // (because the inverse of a diagonal matrix is a diagonal matrix - // with each element the inverse of the original.) - - // The matrix 'matrix' in method 'CSCSegment::weightMatrix()' is this - // block-diagonal overall covariance matrix. It is inverted to the - // block-diagonal error matrix right before it is returned. - - // Use the matrix A defined by - // 1 0 z1 0 - // 0 1 0 z1 - // 1 0 z2 0 - // 0 1 0 z2 - // .. .. .. .. - // 1 0 zn 0 - // 0 1 0 zn - - // The matrix A is returned by 'CSCSegment::derivativeMatrix()'. - - // Then the normal equations are encapsulated in the matrix equation - // - // (AT E A)p = (AT E)m - // - // where AT is the transpose of A. - // We'll call the combined matrix on the LHS, M, and that on the RHS, B: - // M p = B - - // We solve this for the parameter vector, p. - // The elements of M and B then involve sums over the hits - - // The 4 values in p are returned by 'CSCSegment::parameters()' - // in the order uz, vz, u0, v0. - // The error matrix of the parameters is obtained by - // (AT E A)^-1 - // calculated in 'CSCSegment::parametersErrors()'. - - // NOTE 1 - // It does the #hits = 2 case separately. - // (I hope they're not on the same layer! They should not be, by construction.) - - // NOTE 2 - // We need local position of a RecHit w.r.t. the CHAMBER - // and the RecHit itself only knows its local position w.r.t. - // the LAYER, so we must explicitly transform global position. - - CLHEP::HepMatrix M(4,4,0); - CLHEP::HepVector B(4,0); - - ChamberHitContainer::const_iterator ih = proto_segment.begin(); - - for (ih = proto_segment.begin(); ih != proto_segment.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); // FIX !! - - // 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 MUST BE CHECKED IF COMAPTIBLE - CLHEP::HepMatrix IC(2,2); - IC(1,1) = hit.localPositionError().xx(); - IC(1,2) = hit.localPositionError().xy(); - IC(2,1) = IC(1,2); // since Cov is symmetric - IC(2,2) = hit.localPositionError().yy(); - - // 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"; - - // @@ NOW WHAT TO DO? Exception? Return? Ignore? - } - - // ptc: Note that IC is no longer Cov but Cov^-1 - 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; - } - - // Solve the matrix equation using CLHEP's 'solve' - //@@ ptc: CAN solve FAIL?? UNCLEAR FROM (LACK OF) CLHEP DOC - CLHEP::HepVector p = solve(M, B); - - // Update member variables uz, vz, theOrigin - theOrigin = LocalPoint(p(1), p(2), 0.); - uz = p(3); - vz = p(4); } -void CSCSegAlgoTC::fillChiSquared() { - - // The chi-squared is (m-Ap)T E (m-Ap) - // where T denotes transpose. - // This collapses to a simple sum over contributions from each - // pair of measurements. - float u0 = theOrigin.x(); - float v0 = theOrigin.y(); - double chsq = 0.; +void CSCSegAlgoTC::updateParameters(void) { - ChamberHitContainer::const_iterator ih; - for (ih = proto_segment.begin(); ih != proto_segment.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); // FIX !! - - double hu = lp.x(); - double hv = lp.y(); - double hz = lp.z(); - - double du = u0 + uz * hz - hu; - double dv = v0 + vz * hz - hv; - - CLHEP::HepMatrix IC(2,2); - IC(1,1) = hit.localPositionError().xx(); - IC(1,2) = hit.localPositionError().xy(); - IC(2,1) = IC(1,2); - IC(2,2) = hit.localPositionError().yy(); - - // Invert covariance matrix - int ierr = 0; - IC.invert(ierr); - if (ierr != 0) { - LogDebug("CSC") << "CSCSegment::fillChiSquared: failed to invert covariance matrix=\n" << IC << "\n"; - - // @@ NOW WHAT TO DO? Exception? Return? Ignore? - } - - chsq += du*du*IC(1,1) + 2.*du*dv*IC(1,2) + dv*dv*IC(2,2); - } - theChi2 = chsq; -} + //@@ DO NOT DELETE EXISTING FIT SINCE WE SAVE IT!! + // delete sfit_; + sfit_ = new CSCSegFit( theChamber, proto_segment ); + sfit_->fit(); -void CSCSegAlgoTC::fillLocalDirection() { - - // Always enforce direction of segment to point from IP outwards - // (Incorrect for particles not coming from IP, of course.) - - double dxdz = uz; - double dydz = vz; - 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( theOrigin ) ).z(); - double globalZdir = ( theChamber->toGlobal( localDir ) ).z(); - double directionSign = globalZpos * globalZdir; - - theDirection = (directionSign * localDir).unit(); } float CSCSegAlgoTC::phiAtZ(float z) const { // Returns a phi in [ 0, 2*pi ) const CSCLayer* l1 = theChamber->layer(1); - GlobalPoint gp = l1->toGlobal(theOrigin); - GlobalVector gv = l1->toGlobal(theDirection); + GlobalPoint gp = l1->toGlobal(sfit_->intercept()); + GlobalVector gv = l1->toGlobal(sfit_->localdir()); float x = gp.x() + (gv.x()/gv.z())*(z - gp.z()); float y = gp.y() + (gv.y()/gv.z())*(z - gp.z()); @@ -575,55 +317,48 @@ float CSCSegAlgoTC::phiAtZ(float z) const { } void CSCSegAlgoTC::compareProtoSegment(const CSCRecHit2D* h, int layer) { - - // compare the chi2 of two segments - double oldChi2 = theChi2; - LocalPoint oldOrigin = theOrigin; - LocalVector oldDirection = theDirection; - ChamberHitContainer oldSegment = proto_segment; - - bool ok = replaceHit(h, layer); + + // Save copy of current fit + CSCSegFit* oldfit = new CSCSegFit( *sfit_ ); + + bool ok = replaceHit(h, layer); // possible new fit if (ok) { LogDebug("CSC") << " hit in same layer as a hit on segment; try replacing old one..." - << " chi2 new: " << theChi2 << " old: " << oldChi2 << "\n"; + << " chi2 new: " << sfit_->chi2() << " old: " << oldfit->chi2() << "\n"; } - if ((theChi2 < oldChi2) && (ok)) { + if ( ok && (sfit_->chi2() < oldfit->chi2()) ) { LogDebug("CSC") << " segment with replaced hit is better.\n"; + delete oldfit; // new fit is better } else { - proto_segment = oldSegment; - theChi2 = oldChi2; - theOrigin = oldOrigin; - theDirection = oldDirection; + delete sfit_; // new fit is worse + sfit_ = oldfit; // restore original fit } } void CSCSegAlgoTC::increaseProtoSegment(const CSCRecHit2D* h, int layer) { - double oldChi2 = theChi2; - LocalPoint oldOrigin = theOrigin; - LocalVector oldDirection = theDirection; - ChamberHitContainer oldSegment = proto_segment; - - bool ok = addHit(h, layer); + // save copy of input fit + CSCSegFit* oldfit = new CSCSegFit( *sfit_ ); + + bool ok = addHit(h, layer); // possible new fit if (ok) { LogDebug("CSC") << " hit in new layer: added to segment, new chi2: " - << theChi2 << "\n"; + << sfit_->chi2() << "\n"; } - int ndf = 2*proto_segment.size() - 4; + // int ndf = 2*proto_segment.size() - 4; - if (ok && ((ndf <= 0) || (theChi2/ndf < chi2Max))) { + if (ok && ((sfit_->chi2() <= 0) || (sfit_->chi2()/sfit_->ndof() < chi2Max))) { LogDebug("CSC") << " segment with added hit is good.\n" ; + delete oldfit; // new fit is better } else { - proto_segment = oldSegment; - theChi2 = oldChi2; - theOrigin = oldOrigin; - theDirection = oldDirection; + delete sfit_; // new fit is worse + sfit_ = oldfit; // restore original fit } } @@ -664,19 +399,19 @@ bool CSCSegAlgoTC::isHitNearSegment(const CSCRecHit2D* h) const { // in orcarc, or by default, where rxy=sqrt(x**2+y**2) of hit itself. // Note that to make intuitive cuts on delta(phi) one must work in // phi range (-pi, +pi] not [0, 2pi - + const CSCLayer* l1 = theChamber->layer(h->cscDetId().layer()); GlobalPoint hp = l1->toGlobal(h->localPosition()); - float hphi = hp.phi(); // in (-pi, +pi] + float hphi = hp.phi(); // in (-pi, +pi] if (hphi < 0.) - hphi += 2.*M_PI; // into range [0, 2pi) - float sphi = phiAtZ(hp.z()); // in [0, 2*pi) + hphi += 2.*M_PI; // into range [0, 2pi) + float sphi = phiAtZ(hp.z()); // in [0, 2*pi) float phidif = sphi-hphi; if (phidif < 0.) - phidif += 2.*M_PI; // into range [0, 2pi) + phidif += 2.*M_PI; // into range [0, 2pi) if (phidif > M_PI) - phidif -= 2.*M_PI; // into range (-pi, pi] + phidif -= 2.*M_PI; // into range (-pi, pi] float dRPhi = fabs(phidif)*hp.perp(); LogDebug("CSC") << " is hit at phi_h= " << hphi << " near segment phi_seg= " << sphi @@ -717,7 +452,8 @@ void CSCSegAlgoTC::dumpHits(const ChamberHitContainer& rechits) const { } -bool CSCSegAlgoTC::isSegmentGood(std::vector::iterator seg, std::vector::iterator chi2, + +bool CSCSegAlgoTC::isSegmentGood(std::vector::iterator seg, const ChamberHitContainer& rechitsInChamber, BoolContainer& used) const { // Apply any selection cuts to segment @@ -727,34 +463,41 @@ bool CSCSegAlgoTC::isSegmentGood(std::vector::iterator seg, // 2) Ensure no hits on segment are already assigned to another segment // (typically of higher quality) + + size_t iadd = (rechitsInChamber.size() > 20 )? 1 : 0; - unsigned int iadd = (rechitsInChamber.size() > 20 )? 1 : 0; - - if (seg->size() < 3 + iadd) + size_t nhits = (*seg)->nhits(); + + if (nhits < 3 + iadd) return false; // Additional part of alternative segment selection scheme: reject // segments with a chi2 probability of less than chi2ndfProbMin. Relies on list // being sorted with "SegmentSorting == 2", that is first by nrechits and then - // by chi2prob in subgroups of same nr of rechits. + // by chi2prob in subgroups of same no. of rechits. if( SegmentSorting == 2 ){ - if( (*chi2) != 0 && ((2*seg->size())-4) >0 ) { - if ( ChiSquaredProbability((*chi2),(double)(2*seg->size()-4)) < chi2ndfProbMin ) { + double chi2t = (*seg)->chi2(); + double ndoft = 2*nhits - 4 ; + if( chi2t > 0 && ndoft > 0 ) { + if ( ChiSquaredProbability(chi2t,ndoft) < chi2ndfProbMin ) { return false; } } - if((*chi2) == 0 ) return false; + else { + return false; + } } + ChamberHitContainer hits_ = (*seg)->hits(); - for(unsigned int ish = 0; ish < seg->size(); ++ish) { + for(size_t ish = 0; ish < nhits; ++ish) { ChamberHitContainerCIt ib = rechitsInChamber.begin(); for(ChamberHitContainerCIt ic = ib; ic != rechitsInChamber.end(); ++ic) { - if(((*seg)[ish] == (*ic)) && used[ic-ib]) + if((hits_[ish] == (*ic)) && used[ic-ib]) return false; } } @@ -762,17 +505,17 @@ bool CSCSegAlgoTC::isSegmentGood(std::vector::iterator seg, return true; } -void CSCSegAlgoTC::flagHitsAsUsed(std::vector::iterator seg, const ChamberHitContainer& rechitsInChamber, - BoolContainer& used) const { +void CSCSegAlgoTC::flagHitsAsUsed(std::vector::iterator seg, + const ChamberHitContainer& rechitsInChamber, BoolContainer& used) const { // Flag hits on segment as used ChamberHitContainerCIt ib = rechitsInChamber.begin(); + ChamberHitContainer hits = (*seg)->hits(); - for(unsigned int ish = 0; ish < seg->size(); ish++) { - + for(size_t ish = 0; ish < hits.size(); ++ish) { for(ChamberHitContainerCIt iu = ib; iu != rechitsInChamber.end(); ++iu) - if((*seg)[ish] == (*iu)) + if( hits[ish] == (*iu)) used[iu-ib] = true; } } @@ -781,7 +524,7 @@ void CSCSegAlgoTC::pruneTheSegments(const ChamberHitContainer& rechitsInChamber) // Sort the segment store according to segment 'quality' (chi2/#hits ?) and // remove any segments which contain hits assigned to higher-quality segments. - + if (candidates.empty()) return; @@ -792,214 +535,82 @@ void CSCSegAlgoTC::pruneTheSegments(const ChamberHitContainer& rechitsInChamber) segmentSort(); // Select best quality segments, requiring hits are assigned to just one segment - // Because I want to erase the bad segments, the iterator must be incremented + // Want to erase the bad segments, so the iterator must be incremented // inside the loop, and only when the erase is not called - std::vector::iterator is; - std::vector::iterator ichi = chi2s.begin(); - std::vector::iterator iErr = errors.begin(); - std::vector::iterator iOrig = origins.begin(); - std::vector::iterator iDir = directions.begin(); + std::vector::iterator is; for (is = candidates.begin(); is != candidates.end(); ) { - bool goodSegment = isSegmentGood(is, ichi, rechitsInChamber, used); + bool goodSegment = isSegmentGood(is, rechitsInChamber, used); if (goodSegment) { LogDebug("CSC") << "Accepting segment: "; - flagHitsAsUsed(is, rechitsInChamber, used); ++is; - ++ichi; - ++iErr; - ++iOrig; - ++iDir; } else { LogDebug("CSC") << "Rejecting segment: "; - is = candidates.erase(is); - ichi = chi2s.erase(ichi); - iErr = errors.erase(iErr); - iOrig = origins.erase(iOrig); - iDir = directions.erase(iDir); + delete *is; // delete the CSCSegFit* + is = candidates.erase(is); // erase the element in container } } } void CSCSegAlgoTC::segmentSort() { - // The segment collection is sorted according chi2/#hits + // The segment collection is sorted according to e.g. chi2/#hits - for(unsigned int i=0; ihits()).size(); + size_t nj = (candidates[j]->hits()).size(); - int n1 = candidates[i].size(); - int n2 = candidates[j].size(); - - if( SegmentSorting == 2 ){ // Sort criterion: first sort by Nr of rechits, then in groups of rechits by chi2prob: - if ( n2 > n1 ) { // sort by nr of rechits - ChamberHitContainer temp = candidates[j]; + // Sort criterion: first sort by no. of rechits, then in groups of rechits by chi2prob + if( SegmentSorting == 2 ){ + if ( nj > ni ) { // sort by no. of rechits + CSCSegFit* temp = candidates[j]; candidates[j] = candidates[i]; candidates[i] = temp; - - double temp1 = chi2s[j]; - chi2s[j] = chi2s[i]; - chi2s[i] = temp1; - - AlgebraicSymMatrix temp2 = errors[j]; - errors[j] = errors[i]; - errors[i] = temp2; - - LocalPoint temp3 = origins[j]; - origins[j] = origins[i]; - origins[i] = temp3; - - LocalVector temp4 = directions[j]; - directions[j] = directions[i]; - directions[i] = temp4; } // sort by chi2 probability in subgroups with equal nr of rechits - if(chi2s[i] != 0. && 2*n2-4 > 0 ) { - if( n2 == n1 && (ChiSquaredProbability( chi2s[i],(double)(2*n1-4)) < ChiSquaredProbability(chi2s[j],(double)(2*n2-4))) ){ - ChamberHitContainer temp = candidates[j]; + // if(chi2s[i] != 0. && 2*n2-4 > 0 ) { + if( candidates[i]->chi2() > 0 && candidates[i]->ndof() > 0 ) { + if( nj == ni && + ( ChiSquaredProbability( candidates[i]->chi2(),(double)(candidates[i]->ndof()) ) < + ChiSquaredProbability( candidates[j]->chi2(),(double)(candidates[j]->ndof()) ) ) + ){ + CSCSegFit* temp = candidates[j]; candidates[j] = candidates[i]; candidates[i] = temp; - - double temp1 = chi2s[j]; - chi2s[j] = chi2s[i]; - chi2s[i] = temp1; - - AlgebraicSymMatrix temp2 = errors[j]; - errors[j] = errors[i]; - errors[i] = temp2; - - LocalPoint temp3 = origins[j]; - origins[j] = origins[i]; - origins[i] = temp3; - - LocalVector temp4 = directions[j]; - directions[j] = directions[i]; - directions[i] = temp4; } } } else if( SegmentSorting == 1 ){ - if ((chi2s[i]/n1) > (chi2s[j]/n2)) { - - ChamberHitContainer temp = candidates[j]; + if ((candidates[i]->chi2()/ni) > (candidates[j]->chi2()/nj)) { + CSCSegFit* temp = candidates[j]; candidates[j] = candidates[i]; candidates[i] = temp; - - double temp1 = chi2s[j]; - chi2s[j] = chi2s[i]; - chi2s[i] = temp1; - - AlgebraicSymMatrix temp2 = errors[j]; - errors[j] = errors[i]; - errors[i] = temp2; - - LocalPoint temp3 = origins[j]; - origins[j] = origins[i]; - origins[i] = temp3; - - LocalVector temp4 = directions[j]; - directions[j] = directions[i]; - directions[i] = temp4; } } else{ - LogDebug("CSC") << "No valid segment sorting specified - BAD !!!\n"; + LogDebug("CSC") << "No valid segment sorting specified. Algorithm misconfigured! \n"; } } } } -AlgebraicSymMatrix CSCSegAlgoTC::calculateError() const { - - 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 CSCSegAlgoTC::derivativeMatrix() const { - - ChamberHitContainer::const_iterator it; - int nhits = proto_segment.size(); - CLHEP::HepMatrix matrix(2*nhits, 4); - int row = 0; - - for(it = proto_segment.begin(); it != proto_segment.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); // FIX - float z = lp.z(); - ++row; - matrix(row, 1) = 1.; - matrix(row, 3) = z; - ++row; - matrix(row, 2) = 1.; - matrix(row, 4) = z; - } - return matrix; -} - - -AlgebraicSymMatrix CSCSegAlgoTC::weightMatrix() const { - - std::vector::const_iterator it; - int nhits = proto_segment.size(); - AlgebraicSymMatrix matrix(2*nhits, 0); - int row = 0; - - for (it = proto_segment.begin(); it != proto_segment.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; +void CSCSegAlgoTC::dumpSegment( const CSCSegment& seg ) const { + + edm::LogVerbatim("CSCSegment") << "CSCSegment in " << theChamber->id() + << "\nlocal position = " << seg.localPosition() + << "\nerror = " << seg.localPositionError() + << "\nlocal direction = " << seg.localDirection() + << "\nerror =" << seg.localDirectionError() + << "\ncovariance matrix" + << seg.parametersError() + << "chi2/ndf = " << seg.chi2() << "/" << seg.degreesOfFreedom() + << "\n#rechits = " << seg.specificRecHits().size() + << "\ntime = " << seg.time(); } - -void CSCSegAlgoTC::flipErrors( AlgebraicSymMatrix& a ) const { - - // 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 - -} - diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.h b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.h index dbfd30653028a..7e91d7d7408b2 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.h +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.h @@ -18,8 +18,8 @@ * * Ported to CMSSW 2006-04-03: Matteo.Sani@cern.ch
* - * \author M. Sani - * + * Replaced least-squares fit by CSCSegFit - Feb 2015
+ * */ #include @@ -29,21 +29,11 @@ #include #include +class CSCSegFit; class CSCSegAlgoTC : public CSCSegmentAlgorithm { public: - // Tim tried using map as basic container of all (space-point) RecHit's in a chamber: - // The 'key' is a pseudo-layer number (1-6 but with 1 always closest to IP). - // The 'value' is a vector of the RecHit's on that layer. - // Using the layer number like this removes the need to sort in global z. - // Instead we just have to ensure the layer index is correctly adjusted - // to enforce the requirement that 'layer 1' is closest in the chamber - // to the IP. - // However a map is very painful to use when handling the original 'SK' algorithm, - // particularly when we need to flag a hit as 'used'. Because of this I am now - // favouring a pair of vectors, vector and vector for the layer id. - /// Typedefs typedef std::vector LayerIndex; typedef std::vector ChamberHitContainer; @@ -80,10 +70,6 @@ class CSCSegAlgoTC : public CSCSegmentAlgorithm { bool replaceHit(const CSCRecHit2D* h, int layer); void compareProtoSegment(const CSCRecHit2D* h, int layer); void increaseProtoSegment(const CSCRecHit2D* h, int layer); - AlgebraicSymMatrix calculateError() const; - CLHEP::HepMatrix derivativeMatrix() const; - AlgebraicSymMatrix weightMatrix() const; - void flipErrors(AlgebraicSymMatrix&) const; /** * Return true if the difference in (local) x of two hits is < dRPhiMax @@ -100,7 +86,7 @@ class CSCSegAlgoTC : public CSCSegmentAlgorithm { /** * Return true if hit is near segment. * 'Near' means deltaphi and rxy*deltaphi are within ranges - * specified by orcarc parameters dPhiFineMax and dRPhiFineMax, + * specified by config parameters dPhiFineMax and dRPhiFineMax, * where rxy = sqrt(x**2+y**2) of the hit in global coordinates. */ bool isHitNearSegment(const CSCRecHit2D* h) const; @@ -122,15 +108,14 @@ class CSCSegAlgoTC : public CSCSegmentAlgorithm { * In this algorithm, this means it shares no hits with any other segment. * If "SegmentSort=2" also require a minimal chi2 probability of "chi2ndfProbMin". */ - bool isSegmentGood(std::vector::iterator is, - std::vector::iterator ichi, + bool isSegmentGood(std::vector::iterator is, const ChamberHitContainer& rechitsInChamber, BoolContainer& used) const; /** * Flag hits on segment as used */ - void flagHitsAsUsed(std::vector::iterator is, + void flagHitsAsUsed(std::vector::iterator is, const ChamberHitContainer& rechitsInChamber, BoolContainer& used) const; /** @@ -141,29 +126,26 @@ class CSCSegAlgoTC : public CSCSegmentAlgorithm { /** * Sort criterion for segment quality, for use in pruneTheSegments. */ - void segmentSort(); + void segmentSort(void); float phiAtZ(float z) const; - void fillLocalDirection(); - void fillChiSquared(); - void fitSlopes(); - void updateParameters(); + + void updateParameters(void); + + void dumpSegment( const CSCSegment& seg ) const; /// Member variables // ================ const CSCChamber* theChamber; - std::vector candidates; - std::vector origins; - std::vector directions; - std::vector errors; - std::vector chi2s; ChamberHitContainer proto_segment; - double theChi2; - LocalPoint theOrigin; - LocalVector theDirection; - float uz, vz; + + // Pointer to most recent candidate fit + CSCSegFit* sfit_; + + // Store pointers to set of candidate fits + std::vector candidates; /** max segment chi squared */ @@ -209,6 +191,7 @@ class CSCSegAlgoTC : public CSCSegmentAlgorithm { */ const std::string myName; bool debugInfo; + }; #endif From 88abc6253508d39ca13beb76fc4bc644d3316136 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Fri, 20 Feb 2015 12:02:46 +0100 Subject: [PATCH 12/17] remove unused variables --- RecoLocalMuon/CSCSegment/test/CSCSegmentReader.h | 7 +++---- RecoLocalMuon/CSCSegment/test/CSCSegmentVisualise.h | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/test/CSCSegmentReader.h b/RecoLocalMuon/CSCSegment/test/CSCSegmentReader.h index 2bd0f0cff0b27..50a266e25402f 100644 --- a/RecoLocalMuon/CSCSegment/test/CSCSegmentReader.h +++ b/RecoLocalMuon/CSCSegment/test/CSCSegmentReader.h @@ -4,8 +4,7 @@ /** \class CSCSegmentReader * Basic analyzer class which accesses CSCSegment * and plot efficiency of the builder - * - * \author M. Sani + * [This hasn't been used in many years but I am loath to delete code.] */ #include @@ -86,8 +85,8 @@ class CSCSegmentReader : public edm::EDAnalyzer { int simhit; int maxNhits; int minNhits; - int n6hitSegmentMC[9]; - int n6hitSegmentReco[9]; + // int n6hitSegmentMC[9]; + // int n6hitSegmentReco[9]; int near_segment; }; diff --git a/RecoLocalMuon/CSCSegment/test/CSCSegmentVisualise.h b/RecoLocalMuon/CSCSegment/test/CSCSegmentVisualise.h index 42d5e659ebbbd..71c1f58a3b3c8 100644 --- a/RecoLocalMuon/CSCSegment/test/CSCSegmentVisualise.h +++ b/RecoLocalMuon/CSCSegment/test/CSCSegmentVisualise.h @@ -58,7 +58,7 @@ class CSCSegmentVisualise : public edm::EDAnalyzer { int idxHisto; int minRechitChamber; int maxRechitChamber; - double maxPhi, maxTheta; + // double maxPhi, maxTheta; int ievt; }; From 6bce0fdb5905fcc465b2b83b01f9203a203ce2c3 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Fri, 20 Feb 2015 12:04:43 +0100 Subject: [PATCH 13/17] fix mem leak due to too much clean up of code --- RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc index 06dab0a4bfdaf..64a65818cc713 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoSK.cc @@ -242,11 +242,9 @@ void CSCSegAlgoSK::tryAddingHitsToSegment(const ChamberHitContainer& rechits, } bool CSCSegAlgoSK::areHitsCloseInLocalX(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const { - float h1x = h1->localPosition().x(); - float h2x = h2->localPosition().x(); - float deltaX = (h1->localPosition()-h2->localPosition()).x(); - LogDebug("CSC") << " Hits at local x= " << h1x << ", " - << h2x << " have separation= " << deltaX; + float deltaX = ( h1->localPosition() - h2->localPosition() ).x(); + LogDebug("CSC") << " Hits at local x= " << h1->localPosition().x() << ", " + << h2->localPosition().x() << " have separation= " << deltaX; return (fabs(deltaX) < (dRPhiMax * windowScale))? true:false; // +v } @@ -442,6 +440,7 @@ void CSCSegAlgoSK::compareProtoSegment(const CSCRecHit2D* h, int layer) { if ( ( sfit_->chi2() < oldfit->chi2() ) && ok ) { LogDebug("CSC") << " segment with replaced hit is better.\n"; + delete oldfit; // new fit is better } else { // keep original fit @@ -468,6 +467,7 @@ void CSCSegAlgoSK::increaseProtoSegment(const CSCRecHit2D* h, int layer) { //@@ TEST ON ndof<=0 IS JUST TO ACCEPT nhits=2 CASE?? if ( ok && ( (sfit_->ndof() <= 0) || (sfit_->chi2()/sfit_->ndof() < chi2Max)) ) { LogDebug("CSCSegment") << " segment with added hit is good.\n" ; + delete oldfit; // new fit is better } else { // reset to original fit From f1e4d9d093d1d826589557e1fd65370039d842d8 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Fri, 20 Feb 2015 12:06:09 +0100 Subject: [PATCH 14/17] remove unused variables --- RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc index da74bd2c2b3b4..8dcc0470158d6 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoTC.cc @@ -363,11 +363,9 @@ void CSCSegAlgoTC::increaseProtoSegment(const CSCRecHit2D* h, int layer) { } bool CSCSegAlgoTC::areHitsCloseInLocalX(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const { - float h1x = h1->localPosition().x(); - float h2x = h2->localPosition().x(); float deltaX = (h1->localPosition()-h2->localPosition()).x(); - LogDebug("CSC") << " Hits at local x= " << h1x << ", " - << h2x << " have separation= " << deltaX; + LogDebug("CSC") << " Hits at local x= " << h1->localPosition().x() << ", " + << h2->localPosition().x() << " have separation= " << deltaX; return (fabs(deltaX) < (dRPhiMax))? true:false; // +v } From 1b299aa26073f0eed08d4438e153b5df14e2b26c Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Fri, 20 Feb 2015 12:06:25 +0100 Subject: [PATCH 15/17] remove unused variables --- RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc index 63630f10f60f3..687cbde9675ee 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoST.cc @@ -1,6 +1,6 @@ /** - * \file CSCSegAlgoST.cc + * \File CSCSegAlgoST.cc * * \authors: S. Stoynev - NU * I. Bloch - FNAL @@ -232,14 +232,11 @@ std::vector CSCSegAlgoST::prune_bad_hits(const CSCChamber* aChamber, //float z_at_target ; //float radius ; - float loc_x_at_target ; - //float loc_y_at_target ; - //float loc_z_at_target ; + 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: From 9bdb2517b2af41f2025cad27328d8b397956b273 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Fri, 20 Feb 2015 12:06:43 +0100 Subject: [PATCH 16/17] remove unused variables --- RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h index 6a3be86c4070c..1b2a6eab3837d 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoDF.h @@ -24,7 +24,7 @@ * 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.
+ * chamber. [Seems to be hardwired rather than using this variable?]
* * Alternative algorithms can be used for the segment building * by writing classes like this, and then selecting which one is actually @@ -126,13 +126,13 @@ class CSCSegAlgoDF : public CSCSegmentAlgorithm { bool debug; bool preClustering; int minHitsForPreClustering; - bool testSeg; + // bool testSeg; bool Pruning; int minLayersApart; int nHitsPerClusterIsShower; - float nSigmaFromSegment; + // float nSigmaFromSegment; int minHitsPerSegment; - int muonsPerChamberMax; + // int muonsPerChamberMax; double dRPhiFineMax; double dPhiFineMax; float tanPhiMax; From 5751c57927402b7c2729ee78245fe12b55f05f20 Mon Sep 17 00:00:00 2001 From: Tim Cox Date: Fri, 20 Feb 2015 12:07:12 +0100 Subject: [PATCH 17/17] fix misnamed header guard --- RecoLocalMuon/CSCSegment/src/CSCSegAlgoPreClustering.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoPreClustering.h b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoPreClustering.h index d2a6804b918c2..4edbdf32661f5 100644 --- a/RecoLocalMuon/CSCSegment/src/CSCSegAlgoPreClustering.h +++ b/RecoLocalMuon/CSCSegment/src/CSCSegAlgoPreClustering.h @@ -1,5 +1,5 @@ #ifndef CSCSegment_CSCSegAlgoPreClustering_h -#define CSCSegAlgoPreClustering_h +#define CSCSegment_CSCSegAlgoPreClustering_h /** * \file CSCSegAlgoPreClustering.h *