|
|
|
#include <CSCSegAlgoDF.h>
This is a modified version of the SK algorithm for building endcap muon track segments out of the rechit's in a CSCChamber.
A CSCSegment is a RecSegment4D, and is built from CSCRecHit2D objects, each of which is a RecHit2DLocalPos.
This builds segments by first creating proto-segments from at least 3 hits. We intend to try all possible pairs of hits to start segment building. 'All possible' means each hit lies on different layers in the chamber. Once a hit has been assigned to a segment, we don't consider it again, THAT IS, FOR THE FIRST PASS ONLY ! In fact, this is one of the possible flaw with the SK algorithms as it sometimes manages to build segments with the wrong starting points. In the DF algorithm, the endpoints are tested as the best starting points in a 2nd and 3rd loop.
Another difference with the from the SK algorithm is that rechits can be added to proto segments if they fall within n sigmas of the projected track within a given layer. Hence, a cylinder isn't used as in the SK algorimthm, which allows for pseudo 2D hits built from wire or strip only hits to be used in segment reconstruction.
Also, only a certain muonsPerChamberMax maximum number of segments can be produced in the chamber.
Alternative algorithms can be used for the segment building by writing classes like this, and then selecting which one is actually used via the CSCSegmentBuilder.
Definition at line 47 of file CSCSegAlgoDF.h.
| typedef std::deque<bool> CSCSegAlgoDF::BoolContainer |
Definition at line 57 of file CSCSegAlgoDF.h.
| typedef std::vector<const CSCRecHit2D*> CSCSegAlgoDF::ChamberHitContainer |
Definition at line 55 of file CSCSegAlgoDF.h.
| typedef std::vector<const CSCRecHit2D*>::const_iterator CSCSegAlgoDF::ChamberHitContainerCIt |
Definition at line 56 of file CSCSegAlgoDF.h.
| typedef std::vector<int> CSCSegAlgoDF::LayerIndex |
Typedefs.
Definition at line 54 of file CSCSegAlgoDF.h.
| CSCSegAlgoDF::CSCSegAlgoDF | ( | const edm::ParameterSet & | ps | ) | [explicit] |
Constructor.
Definition at line 32 of file CSCSegAlgoDF.cc.
References chi2Max, debug, dPhiFineMax, dRPhiFineMax, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), maxRatioResidual, minHitsForPreClustering, minHitsPerSegment, minLayersApart, nHitsPerClusterIsShower, preCluster_, preClustering, Pruning, showering_, tanPhiMax, and tanThetaMax.
: CSCSegmentAlgorithm(ps), myName("CSCSegAlgoDF") { debug = ps.getUntrackedParameter<bool>("CSCSegmentDebug"); minLayersApart = ps.getParameter<int>("minLayersApart"); minHitsPerSegment = ps.getParameter<int>("minHitsPerSegment"); dRPhiFineMax = ps.getParameter<double>("dRPhiFineMax"); dPhiFineMax = ps.getParameter<double>("dPhiFineMax"); tanThetaMax = ps.getParameter<double>("tanThetaMax"); tanPhiMax = ps.getParameter<double>("tanPhiMax"); chi2Max = ps.getParameter<double>("chi2Max"); preClustering = ps.getUntrackedParameter<bool>("preClustering"); minHitsForPreClustering= ps.getParameter<int>("minHitsForPreClustering"); nHitsPerClusterIsShower= ps.getParameter<int>("nHitsPerClusterIsShower"); Pruning = ps.getUntrackedParameter<bool>("Pruning"); maxRatioResidual = ps.getParameter<double>("maxRatioResidualPrune"); preCluster_ = new CSCSegAlgoPreClustering( ps ); showering_ = new CSCSegAlgoShowering( ps ); }
| CSCSegAlgoDF::~CSCSegAlgoDF | ( | ) | [virtual] |
Destructor.
Definition at line 56 of file CSCSegAlgoDF.cc.
References preCluster_, and showering_.
{
delete preCluster_;
delete showering_;
}
| bool CSCSegAlgoDF::addHit | ( | const CSCRecHit2D * | hit, |
| int | layer | ||
| ) | [private] |
Definition at line 372 of file CSCSegAlgoDF.cc.
References convertSQLiteXML::ok, and protoSegment.
Referenced by compareProtoSegment(), and tryAddingHitsToSegment().
{
// Return true if hit was added successfully and then parameters are updated.
// Return false if there is already a hit on the same layer, or insert failed.
bool ok = true;
// Test that we are not trying to add the same hit again
for ( ChamberHitContainer::const_iterator it = protoSegment.begin(); it != protoSegment.end(); it++ )
if ( aHit == (*it) ) return false;
protoSegment.push_back(aHit);
return ok;
}
| std::vector< CSCSegment > CSCSegAlgoDF::buildSegments | ( | ChamberHitContainer | rechits | ) |
Build track segments in this chamber (this is where the actual segment-building algorithm hides.)
Definition at line 106 of file CSCSegAlgoDF.cc.
References calculateError(), chi2Max, flagHitsAsUsed(), flipErrors(), i, CSCChamber::layer(), CSCRecHit2D::localPosition(), minHitsPerSegment, minLayersApart, nHitsPerClusterIsShower, CSCSegment::nRecHits(), GeomDet::position(), preClustering, protoChi2, protoDirection, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, pruneFromResidual(), Pruning, showering_, CSCSegAlgoShowering::showerSeg(), mathSSE::sqrt(), tanPhiMax, tanThetaMax, groupFilesInBlocks::temp, theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), tryAddingHitsToSegment(), csvLumiCalc::unit, usedHits, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), and z.
Referenced by run().
{
// Clear buffer for segment vector
std::vector<CSCSegment> segmentInChamber;
segmentInChamber.clear();
unsigned nHitInChamber = rechits.size();
if ( nHitInChamber < 3 ) return segmentInChamber;
LayerIndex layerIndex( nHitInChamber );
unsigned nLayers = 0;
int old_layer = -1;
for ( unsigned int i = 0; i < nHitInChamber; i++ ) {
int this_layer = rechits[i]->cscDetId().layer();
layerIndex[i] = this_layer;
if ( this_layer != old_layer ) {
old_layer = this_layer;
nLayers++;
}
}
if ( nLayers < 3 ) return segmentInChamber;
double z1 = theChamber->layer(1)->position().z();
double z6 = theChamber->layer(6)->position().z();
if ( z1 > 0. ) {
if ( z1 > z6 ) {
reverse( layerIndex.begin(), layerIndex.end() );
reverse( rechits.begin(), rechits.end() );
}
}
else if ( z1 < 0. ) {
if ( z1 < z6 ) {
reverse( layerIndex.begin(), layerIndex.end() );
reverse( rechits.begin(), rechits.end() );
}
}
// Showering muon
if ( preClustering && int(nHitInChamber) > nHitsPerClusterIsShower && nLayers > 2 ) {
CSCSegment segShower = showering_->showerSeg(theChamber, rechits);
// Make sure have at least 3 hits...
if ( segShower.nRecHits() < 3 ) return segmentInChamber;
segmentInChamber.push_back(segShower);
return segmentInChamber;
}
// Initialize flags that a given hit has been allocated to a segment
BoolContainer used_ini(rechits.size(), false);
usedHits = used_ini;
ChamberHitContainerCIt ib = rechits.begin();
ChamberHitContainerCIt ie = rechits.end();
// Now Loop over hits within the chamber to find 1st seed for segment building
for ( ChamberHitContainerCIt i1 = ib; i1 < ie; ++i1 ) {
if ( usedHits[i1-ib] ) continue;
const CSCRecHit2D* h1 = *i1;
int layer1 = layerIndex[i1-ib];
const CSCLayer* l1 = theChamber->layer(layer1);
GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
LocalPoint lp1 = theChamber->toLocal(gp1);
// Loop over hits backward to find 2nd seed for segment building
for ( ChamberHitContainerCIt i2 = ie-1; i2 > ib; --i2 ) {
if ( usedHits[i2-ib] ) continue; // Hit has been used already
int layer2 = layerIndex[i2-ib];
if ( (layer2 - layer1) < minLayersApart ) continue;
const CSCRecHit2D* h2 = *i2;
const CSCLayer* l2 = theChamber->layer(layer2);
GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
LocalPoint lp2 = theChamber->toLocal(gp2);
// Clear proto segment so it can be (re)-filled
protoSegment.clear();
// localPosition is position of hit wrt layer (so local z = 0)
protoIntercept = h1->localPosition();
// We want hit wrt chamber (and local z will be != 0)
float dz = gp2.z()-gp1.z();
protoSlope_u = (lp2.x() - lp1.x())/dz ;
protoSlope_v = (lp2.y() - lp1.y())/dz ;
// Test if entrance angle is roughly pointing towards IP
if (fabs(protoSlope_v) > tanThetaMax) continue;
if (fabs(protoSlope_u) > tanPhiMax ) continue;
protoSegment.push_back(h1);
protoSegment.push_back(h2);
// Try adding hits to proto segment
tryAddingHitsToSegment(rechits, i1, i2, layerIndex);
// Check no. of hits on segment to see if segment is large enough
bool segok = true;
unsigned iadd = 0;
if (protoSegment.size() < minHitsPerSegment+iadd) segok = false;
if ( Pruning && segok ) pruneFromResidual();
// Check if segment satisfies chi2 requirement
if (protoChi2 > chi2Max) segok = false;
if ( segok ) {
// Fill segment properties
// Local direction
double dz = 1./sqrt(1. + protoSlope_u*protoSlope_u + protoSlope_v*protoSlope_v);
double dx = dz * protoSlope_u;
double dy = dz * protoSlope_v;
LocalVector localDir(dx,dy,dz);
// localDir may need sign flip to ensure it points outward from IP
double globalZpos = ( theChamber->toGlobal( protoIntercept ) ).z();
double globalZdir = ( theChamber->toGlobal( localDir ) ).z();
double directionSign = globalZpos * globalZdir;
protoDirection = (directionSign * localDir).unit();
// Error matrix
AlgebraicSymMatrix protoErrors = calculateError();
// but reorder components to match what's required by TrackingRecHit interface
// i.e. slopes first, then positions
flipErrors( protoErrors );
CSCSegment temp(protoSegment, protoIntercept, protoDirection, protoErrors, protoChi2);
segmentInChamber.push_back(temp);
if (nHitInChamber-protoSegment.size() < 3) return segmentInChamber;
if (segmentInChamber.size() > 4) return segmentInChamber;
// Flag used hits
flagHitsAsUsed(rechits);
}
}
}
return segmentInChamber;
}
| AlgebraicSymMatrix CSCSegAlgoDF::calculateError | ( | void | ) | const [private] |
Definition at line 632 of file CSCSegAlgoDF.cc.
References funct::A, derivativeMatrix(), query::result, weightMatrix(), and create_public_pileup_plots::weights.
Referenced by buildSegments().
{
AlgebraicSymMatrix weights = weightMatrix();
AlgebraicMatrix A = derivativeMatrix();
// (AT W A)^-1
// from https://www.phys.ufl.edu/~avery/fitting.html, part I
int ierr;
AlgebraicSymMatrix result = weights.similarityT(A);
result.invert(ierr);
// blithely assuming the inverting never fails...
return result;
}
| void CSCSegAlgoDF::compareProtoSegment | ( | const CSCRecHit2D * | h, |
| int | layer | ||
| ) | [private] |
Definition at line 518 of file CSCSegAlgoDF.cc.
References addHit(), convertSQLiteXML::ok, protoChi2, protoDirection, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, and updateParameters().
Referenced by tryAddingHitsToSegment().
{
// Store old segment first
double old_protoChi2 = protoChi2;
LocalPoint old_protoIntercept = protoIntercept;
float old_protoSlope_u = protoSlope_u;
float old_protoSlope_v = protoSlope_v;
LocalVector old_protoDirection = protoDirection;
ChamberHitContainer old_protoSegment = protoSegment;
// Try adding the hit to existing segment, and remove old one existing in same layer
ChamberHitContainer::iterator it;
for ( it = protoSegment.begin(); it != protoSegment.end(); ) {
if ( (*it)->cscDetId().layer() == layer ) {
it = protoSegment.erase(it);
} else {
++it;
}
}
bool ok = addHit(h, layer);
if (ok) updateParameters();
if ( (protoChi2 > old_protoChi2) || ( !ok ) ) {
protoChi2 = old_protoChi2;
protoIntercept = old_protoIntercept;
protoSlope_u = old_protoSlope_u;
protoSlope_v = old_protoSlope_v;
protoDirection = old_protoDirection;
protoSegment = old_protoSegment;
}
}
| CLHEP::HepMatrix CSCSegAlgoDF::derivativeMatrix | ( | void | ) | const [private] |
Definition at line 605 of file CSCSegAlgoDF.cc.
References CSCRecHit2D::cscDetId(), CSCDetId::layer(), CSCChamber::layer(), CSCRecHit2D::localPosition(), makeMuonMisalignmentScenario::matrix, protoSegment, theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::z(), and z.
Referenced by calculateError().
{
ChamberHitContainer::const_iterator it;
int nhits = protoSegment.size();
CLHEP::HepMatrix matrix(2*nhits, 4);
int row = 0;
for(it = protoSegment.begin(); it != protoSegment.end(); ++it) {
const CSCRecHit2D& hit = (**it);
const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer());
GlobalPoint gp = layer->toGlobal(hit.localPosition());
LocalPoint lp = theChamber->toLocal(gp);
float z = lp.z();
++row;
matrix(row, 1) = 1.;
matrix(row, 3) = z;
++row;
matrix(row, 2) = 1.;
matrix(row, 4) = z;
}
return matrix;
}
| void CSCSegAlgoDF::flagHitsAsUsed | ( | const ChamberHitContainer & | rechitsInChamber | ) | [private] |
Flag hits on segment as used
Definition at line 558 of file CSCSegAlgoDF.cc.
References closeHits, protoSegment, and usedHits.
Referenced by buildSegments().
{
// Flag hits on segment as used
ChamberHitContainerCIt ib = rechitsInChamber.begin();
ChamberHitContainerCIt hi, iu;
for ( hi = protoSegment.begin(); hi != protoSegment.end(); ++hi ) {
for ( iu = ib; iu != rechitsInChamber.end(); ++iu ) {
if (*hi == *iu) usedHits[iu-ib] = true;
}
}
// Don't reject hits marked as "nearby" for now.
// So this is bypassed at all times for now !!!
// Perhaps add back to speed up algorithm some more
if (closeHits.size() > 0) return;
for ( hi = closeHits.begin(); hi != closeHits.end(); ++hi ) {
for ( iu = ib; iu != rechitsInChamber.end(); ++iu ) {
if (*hi == *iu) usedHits[iu-ib] = true;
}
}
}
| void CSCSegAlgoDF::flipErrors | ( | AlgebraicSymMatrix & | a | ) | const [private] |
Definition at line 647 of file CSCSegAlgoDF.cc.
References a.
Referenced by buildSegments().
{
// The CSCSegment needs the error matrix re-arranged
AlgebraicSymMatrix hold( a );
// errors on slopes into upper left
a(1,1) = hold(3,3);
a(1,2) = hold(3,4);
a(2,1) = hold(4,3);
a(2,2) = hold(4,4);
// errors on positions into lower right
a(3,3) = hold(1,1);
a(3,4) = hold(1,2);
a(4,3) = hold(2,1);
a(4,4) = hold(2,2);
// off-diagonal elements remain unchanged
}
| bool CSCSegAlgoDF::hasHitOnLayer | ( | int | layer | ) | const [private] |
Definition at line 502 of file CSCSegAlgoDF.cc.
References protoSegment.
Referenced by tryAddingHitsToSegment().
{
// Is there already a hit on this layer?
for ( ChamberHitContainerCIt it = protoSegment.begin(); it != protoSegment.end(); it++ )
if ( (*it)->cscDetId().layer() == layer ) return true;
return false;
}
| bool CSCSegAlgoDF::isHitNearSegment | ( | const CSCRecHit2D * | h | ) | const [private] |
Definition at line 335 of file CSCSegAlgoDF.cc.
References CSCRecHit2D::cscDetId(), SiPixelRawToDigiRegional_cfi::deltaPhi, dPhiFineMax, dRPhiFineMax, CSCDetId::layer(), CSCChamber::layer(), CSCRecHit2D::localPosition(), M_PI, PV3DBase< T, PVType, FrameType >::phi(), protoIntercept, protoSlope_u, protoSlope_v, dttmaxenums::R, mathSSE::sqrt(), theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), and z.
Referenced by tryAddingHitsToSegment().
{
const CSCLayer* layer = theChamber->layer(hit->cscDetId().layer());
// hit phi position in global coordinates
GlobalPoint Hgp = layer->toGlobal(hit->localPosition());
double Hphi = Hgp.phi();
if (Hphi < 0.) Hphi += 2.*M_PI;
LocalPoint Hlp = theChamber->toLocal(Hgp);
double z = Hlp.z();
double LocalX = protoIntercept.x() + protoSlope_u * z;
double LocalY = protoIntercept.y() + protoSlope_v * z;
LocalPoint Slp(LocalX, LocalY, z);
GlobalPoint Sgp = theChamber->toGlobal(Slp);
double Sphi = Sgp.phi();
if (Sphi < 0.) Sphi += 2.*M_PI;
double R = sqrt(Sgp.x()*Sgp.x() + Sgp.y()*Sgp.y());
double deltaPhi = Sphi - Hphi;
if (deltaPhi > 2.*M_PI) deltaPhi -= 2.*M_PI;
if (deltaPhi < -2.*M_PI) deltaPhi += 2.*M_PI;
if (deltaPhi < 0.) deltaPhi = -deltaPhi;
double RdeltaPhi = R * deltaPhi;
if (RdeltaPhi < dRPhiFineMax && deltaPhi < dPhiFineMax ) return true;
return false;
}
| void CSCSegAlgoDF::orderSecondSeed | ( | GlobalPoint | gp1, |
| const ChamberHitContainerCIt | i1, | ||
| const ChamberHitContainerCIt | i2, | ||
| const ChamberHitContainer & | rechits, | ||
| LayerIndex | layerIndex | ||
| ) | [private] |
Order the hits on the 2nd layer for seed building
Definition at line 743 of file CSCSegAlgoDF.cc.
References secondSeedHits.
{
secondSeedHits.clear();
//ChamberHitContainerCIt ib = rechits.begin();
ChamberHitContainerCIt ie = rechits.end();
// int layer1 = layerIndex[i1-ib];
// int layer2 = layerIndex[i2-ib];
// Now fill vector of rechits closest to center of mass:
// secondSeedHitsIdx.clear() = 0;
// Loop over all hits and find hit closest to 1st seed.
for ( ChamberHitContainerCIt i2 = ie-1; i2 > i1; --i2 ) {
}
}
| void CSCSegAlgoDF::pruneFromResidual | ( | ) | [private] |
Prune bad segment from the worse hit based on residuals
Definition at line 671 of file CSCSegAlgoDF.cc.
References CSCRecHit2D::cscDetId(), j, CSCDetId::layer(), CSCChamber::layer(), CSCRecHit2D::localPosition(), maxRatioResidual, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, mathSSE::sqrt(), theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), updateParameters(), v, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), and z.
Referenced by buildSegments().
{
// Only prune if have at least 5 hits
if ( protoSegment.size() < 5 ) return ;
// Now Study residuals
float maxResidual = 0.;
float sumResidual = 0.;
int nHits = 0;
int badIndex = -1;
int j = 0;
ChamberHitContainer::const_iterator ih;
for ( ih = protoSegment.begin(); ih != protoSegment.end(); ++ih ) {
const CSCRecHit2D& hit = (**ih);
const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer());
GlobalPoint gp = layer->toGlobal(hit.localPosition());
LocalPoint lp = theChamber->toLocal(gp);
double u = lp.x();
double v = lp.y();
double z = lp.z();
double du = protoIntercept.x() + protoSlope_u * z - u;
double dv = protoIntercept.y() + protoSlope_v * z - v;
float residual = sqrt(du*du + dv*dv);
sumResidual += residual;
nHits++;
if ( residual > maxResidual ) {
maxResidual = residual;
badIndex = j;
}
j++;
}
float corrAvgResidual = (sumResidual - maxResidual)/(nHits -1);
// Keep all hits
if ( maxResidual/corrAvgResidual < maxRatioResidual ) return;
// Drop worse hit and recompute segment properties + fill
ChamberHitContainer newProtoSegment;
j = 0;
for ( ih = protoSegment.begin(); ih != protoSegment.end(); ++ih ) {
if ( j != badIndex ) newProtoSegment.push_back(*ih);
j++;
}
protoSegment.clear();
for ( ih = newProtoSegment.begin(); ih != newProtoSegment.end(); ++ih ) {
protoSegment.push_back(*ih);
}
// Update segment parameters
updateParameters();
}
| std::vector< CSCSegment > CSCSegAlgoDF::run | ( | const CSCChamber * | aChamber, |
| ChamberHitContainer | rechits | ||
| ) |
Here we must implement the algorithm
Definition at line 65 of file CSCSegAlgoDF.cc.
References buildSegments(), CSCSegAlgoPreClustering::clusterHits(), minHitsForPreClustering, preCluster_, preClustering, and theChamber.
{
// Store chamber info in temp memory
theChamber = aChamber;
int nHits = rechits.size();
// Segments prior to pruning
std::vector<CSCSegment> segments_temp;
if ( preClustering && nHits > minHitsForPreClustering ) {
// This is where the segment origin is in the chamber on avg.
std::vector<CSCSegment> testSegments;
std::vector<ChamberHitContainer> clusteredHits = preCluster_->clusterHits(theChamber, rechits);
// loop over the found clusters:
for (std::vector<ChamberHitContainer>::iterator subrechits = clusteredHits.begin(); subrechits != clusteredHits.end(); ++subrechits ) {
// build the subset of segments:
std::vector<CSCSegment> segs = buildSegments( (*subrechits) );
// add the found subset of segments to the collection of all segments in this chamber:
segments_temp.insert( segments_temp.end(), segs.begin(), segs.end() );
}
} else {
std::vector<CSCSegment> segs = buildSegments( rechits );
// add the found subset of segments to the collection of all segments in this chamber:
segments_temp.insert( segments_temp.end(), segs.begin(), segs.end() );
}
return segments_temp;
}
| void CSCSegAlgoDF::tryAddingHitsToSegment | ( | const ChamberHitContainer & | rechitsInChamber, |
| const ChamberHitContainerCIt | i1, | ||
| const ChamberHitContainerCIt | i2, | ||
| LayerIndex | layerIndex | ||
| ) | [private] |
Utility functions.
Try adding non-used hits to segment
Skip the layers containing the segment endpoints on first 2 passes, but then
try hits on layer containing the segment starting points on 2nd and/or 3rd pass
if segment has >2 hits. Test each hit on the other layers to see if it is near
the segment using rechit error matrix.
If it is, see whether there is already a hit on the segment from the same layer
Definition at line 265 of file CSCSegAlgoDF.cc.
References addHit(), closeHits, compareProtoSegment(), h, hasHitOnLayer(), i, isHitNearSegment(), protoSegment, updateParameters(), and usedHits.
Referenced by buildSegments().
{
/* Iterate over the layers with hits in the chamber
* Skip the layers containing the segment endpoints on first pass, but then
* try hits on layer containing the segment starting points on 2nd pass
* if segment has >2 hits. Once a hit is added to a layer, don't replace it
* until 2nd iteration
*/
ChamberHitContainerCIt ib = rechits.begin();
ChamberHitContainerCIt ie = rechits.end();
closeHits.clear();
for ( ChamberHitContainerCIt i = ib; i != ie; ++i ) {
if (i == i1 || i == i2 ) continue;
if ( usedHits[i-ib] ) continue; // Don't use hits already part of a segment.
const CSCRecHit2D* h = *i;
int layer = layerIndex[i-ib];
int layer1 = layerIndex[i1-ib];
int layer2 = layerIndex[i2-ib];
// Low multiplicity case
if (rechits.size() < 9) {
if ( isHitNearSegment( h ) ) {
if ( !hasHitOnLayer(layer) ) {
addHit(h, layer);
} else {
closeHits.push_back(h);
}
}
// High multiplicity case
} else {
if ( isHitNearSegment( h ) ) {
if ( !hasHitOnLayer(layer) ) {
addHit(h, layer);
updateParameters();
// Don't change the starting points at this stage !!!
} else {
closeHits.push_back(h);
if (layer != layer1 && layer != layer2 ) compareProtoSegment(h, layer);
}
}
}
}
if ( int(protoSegment.size()) < 3) return;
updateParameters();
// 2nd pass to remove biases
// This time, also consider changing the endpoints
for ( ChamberHitContainerCIt i = closeHits.begin() ; i != closeHits.end(); ++i ) {
const CSCRecHit2D* h = *i;
int layer = (*i)->cscDetId().layer();
compareProtoSegment(h, layer);
}
}
| void CSCSegAlgoDF::updateParameters | ( | void | ) | [private] |
Definition at line 394 of file CSCSegAlgoDF.cc.
References CSCRecHit2D::cscDetId(), CSCDetId::layer(), CSCChamber::layer(), CSCRecHit2D::localPosition(), CSCRecHit2D::localPositionError(), LogDebug, AlCaHLTBitMon_ParallelJobs::p, protoChi2, protoIntercept, protoSegment, protoSlope_u, protoSlope_v, theChamber, GeomDet::toGlobal(), GeomDet::toLocal(), v, PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), LocalError::xy(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), PV3DBase< T, PVType, FrameType >::z(), and z.
Referenced by compareProtoSegment(), pruneFromResidual(), and tryAddingHitsToSegment().
{
// Compute slope from Least Square Fit
CLHEP::HepMatrix M(4,4,0);
CLHEP::HepVector B(4,0);
ChamberHitContainer::const_iterator ih;
for (ih = protoSegment.begin(); ih != protoSegment.end(); ++ih) {
const CSCRecHit2D& hit = (**ih);
const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer());
GlobalPoint gp = layer->toGlobal(hit.localPosition());
LocalPoint lp = theChamber->toLocal(gp);
double u = lp.x();
double v = lp.y();
double z = lp.z();
// ptc: Covariance matrix of local errors
CLHEP::HepMatrix IC(2,2);
IC(1,1) = hit.localPositionError().xx();
IC(1,2) = hit.localPositionError().xy();
IC(2,2) = hit.localPositionError().yy();
IC(2,1) = IC(1,2); // since Cov is symmetric
// ptc: Invert covariance matrix (and trap if it fails!)
int ierr = 0;
IC.invert(ierr); // inverts in place
if (ierr != 0) {
LogDebug("CSC") << "CSCSegment::fitSlopes: failed to invert covariance matrix=\n" << IC << "\n";
}
M(1,1) += IC(1,1);
M(1,2) += IC(1,2);
M(1,3) += IC(1,1) * z;
M(1,4) += IC(1,2) * z;
B(1) += u * IC(1,1) + v * IC(1,2);
M(2,1) += IC(2,1);
M(2,2) += IC(2,2);
M(2,3) += IC(2,1) * z;
M(2,4) += IC(2,2) * z;
B(2) += u * IC(2,1) + v * IC(2,2);
M(3,1) += IC(1,1) * z;
M(3,2) += IC(1,2) * z;
M(3,3) += IC(1,1) * z * z;
M(3,4) += IC(1,2) * z * z;
B(3) += ( u * IC(1,1) + v * IC(1,2) ) * z;
M(4,1) += IC(2,1) * z;
M(4,2) += IC(2,2) * z;
M(4,3) += IC(2,1) * z * z;
M(4,4) += IC(2,2) * z * z;
B(4) += ( u * IC(2,1) + v * IC(2,2) ) * z;
}
CLHEP::HepVector p = solve(M, B);
// Update member variables
// Note that origin has local z = 0
protoIntercept = LocalPoint(p(1), p(2), 0.);
protoSlope_u = p(3);
protoSlope_v = p(4);
// Determine Chi^2 for the proto wire segment
double chsq = 0.;
for (ih = protoSegment.begin(); ih != protoSegment.end(); ++ih) {
const CSCRecHit2D& hit = (**ih);
const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer());
GlobalPoint gp = layer->toGlobal(hit.localPosition());
LocalPoint lp = theChamber->toLocal(gp);
double u = lp.x();
double v = lp.y();
double z = lp.z();
double du = protoIntercept.x() + protoSlope_u * z - u;
double dv = protoIntercept.y() + protoSlope_v * z - v;
CLHEP::HepMatrix IC(2,2);
IC(1,1) = hit.localPositionError().xx();
IC(1,2) = hit.localPositionError().xy();
IC(2,2) = hit.localPositionError().yy();
IC(2,1) = IC(1,2);
// Invert covariance matrix
int ierr = 0;
IC.invert(ierr);
if (ierr != 0) {
LogDebug("CSC") << "CSCSegment::fillChiSquared: failed to invert covariance matrix=\n" << IC << "\n";
}
chsq += du*du*IC(1,1) + 2.*du*dv*IC(1,2) + dv*dv*IC(2,2);
}
protoChi2 = chsq;
}
| AlgebraicSymMatrix CSCSegAlgoDF::weightMatrix | ( | void | ) | const [private] |
Definition at line 582 of file CSCSegAlgoDF.cc.
References CSCRecHit2D::localPositionError(), makeMuonMisalignmentScenario::matrix, protoSegment, LocalError::xx(), LocalError::xy(), and LocalError::yy().
Referenced by calculateError().
{
std::vector<const CSCRecHit2D*>::const_iterator it;
int nhits = protoSegment.size();
AlgebraicSymMatrix matrix(2*nhits, 0);
int row = 0;
for (it = protoSegment.begin(); it != protoSegment.end(); ++it) {
const CSCRecHit2D& hit = (**it);
++row;
matrix(row, row) = hit.localPositionError().xx();
matrix(row, row+1) = hit.localPositionError().xy();
++row;
matrix(row, row-1) = hit.localPositionError().xy();
matrix(row, row) = hit.localPositionError().yy();
}
int ierr;
matrix.invert(ierr);
return matrix;
}
float CSCSegAlgoDF::chi2Max [private] |
Definition at line 158 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), and CSCSegAlgoDF().
ChamberHitContainer CSCSegAlgoDF::closeHits [private] |
Definition at line 133 of file CSCSegAlgoDF.h.
Referenced by flagHitsAsUsed(), and tryAddingHitsToSegment().
bool CSCSegAlgoDF::debug [private] |
Definition at line 144 of file CSCSegAlgoDF.h.
Referenced by CSCSegAlgoDF().
double CSCSegAlgoDF::dPhiFineMax [private] |
Definition at line 155 of file CSCSegAlgoDF.h.
Referenced by CSCSegAlgoDF(), and isHitNearSegment().
double CSCSegAlgoDF::dRPhiFineMax [private] |
Definition at line 154 of file CSCSegAlgoDF.h.
Referenced by CSCSegAlgoDF(), and isHitNearSegment().
float CSCSegAlgoDF::maxRatioResidual [private] |
Definition at line 159 of file CSCSegAlgoDF.h.
Referenced by CSCSegAlgoDF(), and pruneFromResidual().
int CSCSegAlgoDF::minHitsForPreClustering [private] |
Definition at line 146 of file CSCSegAlgoDF.h.
Referenced by CSCSegAlgoDF(), and run().
int CSCSegAlgoDF::minHitsPerSegment [private] |
Definition at line 152 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), and CSCSegAlgoDF().
int CSCSegAlgoDF::minLayersApart [private] |
Definition at line 149 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), and CSCSegAlgoDF().
int CSCSegAlgoDF::muonsPerChamberMax [private] |
Definition at line 153 of file CSCSegAlgoDF.h.
const std::string CSCSegAlgoDF::myName [private] |
Definition at line 129 of file CSCSegAlgoDF.h.
int CSCSegAlgoDF::nHitsPerClusterIsShower [private] |
Definition at line 150 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), and CSCSegAlgoDF().
float CSCSegAlgoDF::nSigmaFromSegment [private] |
Definition at line 151 of file CSCSegAlgoDF.h.
CSCSegAlgoPreClustering* CSCSegAlgoDF::preCluster_ [private] |
Definition at line 161 of file CSCSegAlgoDF.h.
Referenced by CSCSegAlgoDF(), run(), and ~CSCSegAlgoDF().
bool CSCSegAlgoDF::preClustering [private] |
Definition at line 145 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), CSCSegAlgoDF(), and run().
double CSCSegAlgoDF::protoChi2 [private] |
Definition at line 140 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), compareProtoSegment(), and updateParameters().
LocalVector CSCSegAlgoDF::protoDirection [private] |
Definition at line 141 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), and compareProtoSegment().
LocalPoint CSCSegAlgoDF::protoIntercept [private] |
Definition at line 139 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), compareProtoSegment(), isHitNearSegment(), pruneFromResidual(), and updateParameters().
Definition at line 135 of file CSCSegAlgoDF.h.
Referenced by addHit(), buildSegments(), compareProtoSegment(), derivativeMatrix(), flagHitsAsUsed(), hasHitOnLayer(), pruneFromResidual(), tryAddingHitsToSegment(), updateParameters(), and weightMatrix().
float CSCSegAlgoDF::protoSlope_u [private] |
Definition at line 137 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), compareProtoSegment(), isHitNearSegment(), pruneFromResidual(), and updateParameters().
float CSCSegAlgoDF::protoSlope_v [private] |
Definition at line 138 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), compareProtoSegment(), isHitNearSegment(), pruneFromResidual(), and updateParameters().
bool CSCSegAlgoDF::Pruning [private] |
Definition at line 148 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), and CSCSegAlgoDF().
Definition at line 136 of file CSCSegAlgoDF.h.
Referenced by orderSecondSeed().
CSCSegAlgoShowering* CSCSegAlgoDF::showering_ [private] |
Definition at line 162 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), CSCSegAlgoDF(), and ~CSCSegAlgoDF().
float CSCSegAlgoDF::tanPhiMax [private] |
Definition at line 156 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), and CSCSegAlgoDF().
float CSCSegAlgoDF::tanThetaMax [private] |
Definition at line 157 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), and CSCSegAlgoDF().
bool CSCSegAlgoDF::testSeg [private] |
Definition at line 147 of file CSCSegAlgoDF.h.
const CSCChamber* CSCSegAlgoDF::theChamber [private] |
Definition at line 130 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), derivativeMatrix(), isHitNearSegment(), pruneFromResidual(), run(), and updateParameters().
BoolContainer CSCSegAlgoDF::usedHits [private] |
Definition at line 131 of file CSCSegAlgoDF.h.
Referenced by buildSegments(), flagHitsAsUsed(), and tryAddingHitsToSegment().
1.7.3