Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:24:40

0001 /**
0002  * \file CSCSegAlgoDF.cc
0003  *
0004  *  \author Dominique Fortin - UCR
0005  *
0006  * Last update: 17.02.2015 - Tim Cox - use CSCSegFit for least-squares fit
0007  * @@ Seems to find very few segments so many candidates must be being rejected. 
0008  *
0009  */
0010 
0011 #include "CSCSegAlgoDF.h"
0012 #include "CSCSegFit.h"
0013 
0014 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0015 
0016 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0017 
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 
0021 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0022 
0023 #include "CSCSegAlgoPreClustering.h"
0024 #include "CSCSegAlgoShowering.h"
0025 
0026 #include <algorithm>
0027 #include <cmath>
0028 #include <iostream>
0029 #include <string>
0030 
0031 /* Constructor
0032  *
0033  */
0034 CSCSegAlgoDF::CSCSegAlgoDF(const edm::ParameterSet& ps)
0035     : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoDF"), sfit_(nullptr) {
0036   debug = ps.getUntrackedParameter<bool>("CSCSegmentDebug");
0037   minLayersApart = ps.getParameter<int>("minLayersApart");
0038   minHitsPerSegment = ps.getParameter<int>("minHitsPerSegment");
0039   dRPhiFineMax = ps.getParameter<double>("dRPhiFineMax");
0040   dPhiFineMax = ps.getParameter<double>("dPhiFineMax");
0041   tanThetaMax = ps.getParameter<double>("tanThetaMax");
0042   tanPhiMax = ps.getParameter<double>("tanPhiMax");
0043   chi2Max = ps.getParameter<double>("chi2Max");
0044   preClustering = ps.getUntrackedParameter<bool>("preClustering");
0045   minHitsForPreClustering = ps.getParameter<int>("minHitsForPreClustering");
0046   nHitsPerClusterIsShower = ps.getParameter<int>("nHitsPerClusterIsShower");
0047   Pruning = ps.getUntrackedParameter<bool>("Pruning");
0048   maxRatioResidual = ps.getParameter<double>("maxRatioResidualPrune");
0049 
0050   preCluster_ = new CSCSegAlgoPreClustering(ps);
0051   showering_ = new CSCSegAlgoShowering(ps);
0052 }
0053 
0054 /* Destructor
0055  *
0056  */
0057 CSCSegAlgoDF::~CSCSegAlgoDF() {
0058   delete preCluster_;
0059   delete showering_;
0060 }
0061 
0062 /* run
0063  *
0064  */
0065 std::vector<CSCSegment> CSCSegAlgoDF::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
0066   // Store chamber info in temp memory
0067   theChamber = aChamber;
0068 
0069   int nHits = rechits.size();
0070 
0071   // Segments prior to pruning
0072   std::vector<CSCSegment> segments_temp;
0073 
0074   if (preClustering && nHits > minHitsForPreClustering) {
0075     // This is where the segment origin is in the chamber on avg.
0076     std::vector<CSCSegment> testSegments;
0077     std::vector<ChamberHitContainer> clusteredHits = preCluster_->clusterHits(theChamber, rechits);
0078     // loop over the found clusters:
0079     for (std::vector<ChamberHitContainer>::iterator subrechits = clusteredHits.begin();
0080          subrechits != clusteredHits.end();
0081          ++subrechits) {
0082       // build the subset of segments:
0083       std::vector<CSCSegment> segs = buildSegments((*subrechits));
0084       // add the found subset of segments to the collection of all segments in this chamber:
0085       segments_temp.insert(segments_temp.end(), segs.begin(), segs.end());
0086     }
0087   } else {
0088     std::vector<CSCSegment> segs = buildSegments(rechits);
0089     // add the found subset of segments to the collection of all segments in this chamber:
0090     segments_temp.insert(segments_temp.end(), segs.begin(), segs.end());
0091   }
0092 
0093   return segments_temp;
0094 }
0095 
0096 /* This builds segments by first creating proto-segments from at least 3 hits.
0097  * We intend to try all possible pairs of hits to start segment building. 'All possible' 
0098  * means each hit lies on different layers in the chamber.  Once a hit has been assigned 
0099  * to a segment, we don't consider it again, THAT IS, FOR THE FIRST PASS ONLY !
0100  * In fact, this is one of the possible flaw with the SK algorithms as it sometimes manages
0101  * to build segments with the wrong starting points.  In the DF algorithm, the endpoints
0102  * are tested as the best starting points in a 2nd loop.
0103  *
0104  * Also, a  maximum of   5   segments is allowed in the chamber (and then it just gives up.)
0105  *
0106  * @@ There are 7 return's in this function!
0107  *
0108  */
0109 std::vector<CSCSegment> CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _rechits) {
0110   ChamberHitContainer rechits = _rechits;
0111   // Clear buffer for segment vector
0112   std::vector<CSCSegment> segmentInChamber;
0113   segmentInChamber.clear();
0114 
0115   unsigned nHitInChamber = rechits.size();
0116 
0117   //  std::cout << "[CSCSegAlgoDF::buildSegments] address of chamber = " << theChamber << std::endl;
0118   //  std::cout << "[CSCSegAlgoDF::buildSegments] starting in " << theChamber->id()
0119   //         << " with " << nHitInChamber << " rechits" << std::endl;
0120 
0121   // Return #1 - OK, there aren't enough rechits to build a segment
0122   if (nHitInChamber < 3)
0123     return segmentInChamber;
0124 
0125   LayerIndex layerIndex(nHitInChamber);
0126 
0127   size_t nLayers = 0;
0128   size_t old_layer = 0;
0129   for (size_t i = 0; i < nHitInChamber; ++i) {
0130     size_t this_layer = rechits[i]->cscDetId().layer();
0131     //    std::cout << "[CSCSegAlgoDF::buildSegments] this_layer = " << this_layer << std::endl;
0132     layerIndex[i] = this_layer;
0133     //    std::cout << "[CSCSegAlgoDF::buildSegments] layerIndex[" << i << "] = " << layerIndex[i] << std::endl;
0134     if (this_layer != old_layer) {
0135       old_layer = this_layer;
0136       ++nLayers;
0137     }
0138   }
0139 
0140   //  std::cout << "[CSCSegAlgoDF::buildSegments] layers with rechits = " << nLayers << std::endl;
0141 
0142   // Return #2 - OK, there aren't enough hit layers to build a segment
0143   if (nLayers < 3)
0144     return segmentInChamber;
0145 
0146   double z1 = theChamber->layer(1)->position().z();
0147   double z6 = theChamber->layer(6)->position().z();
0148 
0149   if (z1 > 0.) {
0150     if (z1 > z6) {
0151       reverse(layerIndex.begin(), layerIndex.end());
0152       reverse(rechits.begin(), rechits.end());
0153     }
0154   } else if (z1 < 0.) {
0155     if (z1 < z6) {
0156       reverse(layerIndex.begin(), layerIndex.end());
0157       reverse(rechits.begin(), rechits.end());
0158     }
0159   }
0160 
0161   //  std::cout << "[CSCSegAlgoDF::buildSegments] rechits have been ordered" << std::endl;
0162 
0163   // Showering muon
0164   if (preClustering && int(nHitInChamber) > nHitsPerClusterIsShower && nLayers > 2) {
0165     // std::cout << "[CSCSegAlgoDF::buildSegments] showering block" << std::endl;
0166 
0167     CSCSegment segShower = showering_->showerSeg(theChamber, rechits);
0168 
0169     // Return #3 - OK, this is now 'effectve' rechits
0170     // Make sure have at least 3 hits...
0171     if (segShower.nRecHits() < 3)
0172       return segmentInChamber;
0173 
0174     segmentInChamber.push_back(segShower);
0175     if (debug)
0176       dumpSegment(segShower);
0177 
0178     // Return #4 - OK, only get one try at building a segment from shower
0179     return segmentInChamber;
0180   }
0181 
0182   // Initialize flags that a given hit has been allocated to a segment
0183   BoolContainer used_ini(rechits.size(), false);
0184   usedHits = used_ini;
0185 
0186   ChamberHitContainerCIt ib = rechits.begin();
0187   ChamberHitContainerCIt ie = rechits.end();
0188 
0189   //  std::cout << "[CSCSegAlgoDF::buildSegments] entering rechit loop" << std::endl;
0190 
0191   // Now Loop over hits within the chamber to find 1st seed for segment building
0192   for (ChamberHitContainerCIt i1 = ib; i1 < ie; ++i1) {
0193     if (usedHits[i1 - ib])
0194       continue;
0195 
0196     const CSCRecHit2D* h1 = *i1;
0197     int layer1 = layerIndex[i1 - ib];
0198     const CSCLayer* l1 = theChamber->layer(layer1);
0199     GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0200     LocalPoint lp1 = theChamber->toLocal(gp1);
0201 
0202     // Loop over hits backward to find 2nd seed for segment building
0203     for (ChamberHitContainerCIt i2 = ie - 1; i2 > ib; --i2) {
0204       if (usedHits[i2 - ib])
0205         continue;  // Hit has been used already
0206 
0207       int layer2 = layerIndex[i2 - ib];
0208       if ((layer2 - layer1) < minLayersApart)
0209         continue;
0210 
0211       const CSCRecHit2D* h2 = *i2;
0212       const CSCLayer* l2 = theChamber->layer(layer2);
0213       GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0214       LocalPoint lp2 = theChamber->toLocal(gp2);
0215 
0216       // Clear proto segment so it can be (re)-filled
0217       protoSegment.clear();
0218 
0219       // We want hit wrt chamber (and local z will be != 0)
0220       float dz = gp2.z() - gp1.z();
0221       float slope_u = (lp2.x() - lp1.x()) / dz;
0222       float slope_v = (lp2.y() - lp1.y()) / dz;
0223 
0224       // Test if entrance angle is roughly pointing towards IP
0225       if (fabs(slope_v) > tanThetaMax)
0226         continue;
0227       if (fabs(slope_u) > tanPhiMax)
0228         continue;
0229 
0230       protoSegment.push_back(h1);
0231       protoSegment.push_back(h2);
0232 
0233       //      std::cout << "[CSCSegAlgoDF::buildSegments] about to fit 2 hits on layers "
0234       //        << layer1 << " and " << layer2 << std::endl;
0235 
0236       // protoSegment has just 2 hits - but need to create a CSCSegFit to hold it in case
0237       // we fail to add any more hits
0238       updateParameters();
0239 
0240       // Try adding hits to proto segment
0241       tryAddingHitsToSegment(rechits, i1, i2, layerIndex);
0242 
0243       // Check no. of hits on segment to see if segment is large enough
0244       bool segok = true;
0245       unsigned iadd = 0;
0246 
0247       if (protoSegment.size() < minHitsPerSegment + iadd)
0248         segok = false;
0249 
0250       if (Pruning && segok)
0251         pruneFromResidual();
0252 
0253       // Check if segment satisfies chi2 requirement
0254       if (sfit_->chi2() > chi2Max)
0255         segok = false;
0256 
0257       if (segok) {
0258         // Create an actual CSCSegment - retrieve all info from the current fit
0259         CSCSegment temp(sfit_->hits(), sfit_->intercept(), sfit_->localdir(), sfit_->covarianceMatrix(), sfit_->chi2());
0260         //  std::cout << "[CSCSegAlgoDF::buildSegments] about to delete sfit= = " << sfit_ << std::endl;
0261         delete sfit_;
0262         sfit_ = nullptr;  // avoid possibility of attempting a second delete later
0263 
0264         segmentInChamber.push_back(temp);
0265         if (debug)
0266           dumpSegment(temp);
0267 
0268         // Return #5 - OK, fewer than 3 rechits not on this segment left in chamber
0269         if (nHitInChamber - protoSegment.size() < 3)
0270           return segmentInChamber;
0271         // Return $6 - already have more than 4 segments in this chamber
0272         if (segmentInChamber.size() > 4)
0273           return segmentInChamber;
0274 
0275         // Flag used hits
0276         flagHitsAsUsed(rechits);
0277       }
0278     }
0279   }
0280   // Return #7
0281   return segmentInChamber;
0282 }
0283 
0284 /* Method tryAddingHitsToSegment
0285  *
0286  * Look at left over hits and try to add them to proto segment by looking how far they
0287  * are from the segment in terms of the hit error matrix (so how many sigmas away).
0288  *
0289  */
0290 void CSCSegAlgoDF::tryAddingHitsToSegment(const ChamberHitContainer& rechits,
0291                                           const ChamberHitContainerCIt i1,
0292                                           const ChamberHitContainerCIt i2,
0293                                           const LayerIndex& layerIndex) {
0294   /* Iterate over the layers with hits in the chamber
0295  * Skip the layers containing the segment endpoints on first pass, but then
0296  * try hits on layer containing the segment starting points on 2nd pass
0297  * if segment has >2 hits.  Once a hit is added to a layer, don't replace it 
0298  * until 2nd iteration
0299  */
0300 
0301   //  std::cout << "[CSCSegAlgoDF::tryAddingHitsToSegment] entering"
0302   //        << " with rechits.size() = " << rechits.size() << std::endl;
0303 
0304   ChamberHitContainerCIt ib = rechits.begin();
0305   ChamberHitContainerCIt ie = rechits.end();
0306   closeHits.clear();
0307 
0308   //  int counter1 = 0;
0309   //  int counter2 = 0;
0310 
0311   for (ChamberHitContainerCIt i = ib; i != ie; ++i) {
0312     //    std::cout << "counter1 = " << ++counter1 << std::endl;
0313     if (i == i1 || i == i2)
0314       continue;
0315     if (usedHits[i - ib])
0316       continue;  // Don't use hits already part of a segment.
0317 
0318     //    std::cout << "counter2 = " << ++counter2 << std::endl;
0319     const CSCRecHit2D* h = *i;
0320     int layer = layerIndex[i - ib];
0321     int layer1 = layerIndex[i1 - ib];
0322     int layer2 = layerIndex[i2 - ib];
0323 
0324     //    std::cout << "[CSCSegAlgoDF::tryAddingHitsToSegment] layer, layer1, layer2 = "
0325     //        << layer << ", " << layer1 << ", " << layer2 << std::endl;
0326 
0327     // Low multiplicity case
0328     // only adds hit to protoSegment if no hit on layer already; otherwise adds to closeHits
0329     if (rechits.size() < 9) {
0330       //      std::cout << "low mult" << std::endl;
0331       if (isHitNearSegment(h)) {
0332         if (!hasHitOnLayer(layer)) {
0333           addHit(h, layer);
0334         } else {
0335           closeHits.push_back(h);
0336         }
0337       }
0338 
0339       // High multiplicity case
0340       // only adds hit to protoSegment if no hit on layer already AND then refits; otherwise adds to closeHits
0341     } else {
0342       //      std::cout << "high mult" << std::endl;
0343       if (isHitNearSegment(h)) {
0344         //  std::cout << "near seg" << std::endl;
0345         if (!hasHitOnLayer(layer)) {
0346           //      std::cout << "no hit on layer" << std::endl;
0347           if (addHit(h, layer)) {
0348             //      std::cout << "update fit" << std::endl;
0349             updateParameters();
0350           }
0351           // Don't change the starting points at this stage !!!
0352         } else {
0353           //      std::cout << "already hit on layer" << std::endl;
0354           closeHits.push_back(h);
0355           if (layer != layer1 && layer != layer2)
0356             compareProtoSegment(h, layer);
0357         }
0358       }
0359     }
0360   }
0361 
0362   if (int(protoSegment.size()) < 3)
0363     return;
0364   //  std::cout << "final fit" << std::endl;
0365   updateParameters();
0366 
0367   // 2nd pass to remove biases
0368   // This time, also consider changing the endpoints
0369   for (ChamberHitContainerCIt i = closeHits.begin(); i != closeHits.end(); ++i) {
0370     //    std::cout << "2nd pass" << std::endl;
0371     const CSCRecHit2D* h = *i;
0372     int layer = (*i)->cscDetId().layer();
0373     compareProtoSegment(h, layer);
0374   }
0375 }
0376 
0377 /* isHitNearSegment
0378  *
0379  * Compare rechit with expected position from protosegment
0380  */
0381 bool CSCSegAlgoDF::isHitNearSegment(const CSCRecHit2D* hit) const {
0382   const CSCLayer* layer = theChamber->layer(hit->cscDetId().layer());
0383 
0384   // hit phi position in global coordinates
0385   GlobalPoint Hgp = layer->toGlobal(hit->localPosition());
0386   float Hphi = Hgp.phi();
0387   if (Hphi < 0.)
0388     Hphi += 2. * M_PI;
0389   LocalPoint Hlp = theChamber->toLocal(Hgp);
0390   float z = Hlp.z();
0391 
0392   float LocalX = sfit_->xfit(z);
0393   float LocalY = sfit_->yfit(z);
0394 
0395   LocalPoint Slp(LocalX, LocalY, z);
0396   GlobalPoint Sgp = theChamber->toGlobal(Slp);
0397   float Sphi = Sgp.phi();
0398   if (Sphi < 0.)
0399     Sphi += 2. * M_PI;
0400   float R = sqrt(Sgp.x() * Sgp.x() + Sgp.y() * Sgp.y());
0401 
0402   float deltaPhi = Sphi - Hphi;
0403   if (deltaPhi > 2. * M_PI)
0404     deltaPhi -= 2. * M_PI;
0405   if (deltaPhi < -2. * M_PI)
0406     deltaPhi += 2. * M_PI;
0407   if (deltaPhi < 0.)
0408     deltaPhi = -deltaPhi;
0409 
0410   float RdeltaPhi = R * deltaPhi;
0411 
0412   if (RdeltaPhi < dRPhiFineMax && deltaPhi < dPhiFineMax)
0413     return true;
0414 
0415   return false;
0416 }
0417 
0418 /* Method addHit
0419  *
0420  * Test if can add hit to proto segment. If so, try to add it.
0421  *
0422  */
0423 bool CSCSegAlgoDF::addHit(const CSCRecHit2D* aHit, int layer) {
0424   //  std::cout << "[CSCSegAlgoDF::addHit] on layer " << layer << " to protoSegment.size() = "
0425   //        << protoSegment.size() << std::endl;
0426 
0427   // Return true if hit was added successfully and then parameters are updated.
0428   // Return false if there is already a hit on the same layer, or insert failed.
0429 
0430   if (protoSegment.size() > 5)
0431     return false;  //@@ can only have 6 hits at most
0432 
0433   // Test that we are not trying to add the same hit again
0434   for (ChamberHitContainer::const_iterator it = protoSegment.begin(); it != protoSegment.end(); ++it)
0435     if (aHit == (*it))
0436       return false;
0437 
0438   protoSegment.push_back(aHit);
0439 
0440   return true;
0441 }
0442 
0443 /* Method updateParameters
0444  *      
0445  * Perform a simple Least Square Fit on proto segment to determine slope and intercept
0446  *
0447  */
0448 void CSCSegAlgoDF::updateParameters() {
0449   // Delete existing CSCSegFit, create a new one and make the fit
0450   // Uses internal variables - theChamber & protoSegment
0451 
0452   //  std::cout << "[CSCSegAlgoDF::updateParameters] about to delete sfit_ = " << sfit_ << std::endl;
0453   delete sfit_;
0454   //  std::cout << "[CSCSegAlgoDF::updateParameters] protoSegment.size() = "
0455   //                                 << protoSegment.size() << std::endl;
0456   //  std::cout  << "[CSCSegAlgoDF::updateParameters] theChamber = " << theChamber << std::endl;
0457   sfit_ = new CSCSegFit(theChamber, protoSegment);
0458   //  std::cout << "[CSCSegAlgoDF::updateParameters] new sfit_ = " << sfit_ << std::endl;
0459   sfit_->fit();
0460 }
0461 
0462 /* hasHitOnLayer
0463  *
0464  * Just make sure hit to be added to layer comes from different layer than those in proto segment   
0465  */
0466 bool CSCSegAlgoDF::hasHitOnLayer(int layer) const {
0467   //  std::cout << "[CSCSegAlgoDF::hasHitOnLayer] on layer " << layer << std::endl;
0468 
0469   // Is there already a hit on this layer?
0470   for (ChamberHitContainerCIt it = protoSegment.begin(); it != protoSegment.end(); it++)
0471     if ((*it)->cscDetId().layer() == layer)
0472       return true;
0473 
0474   return false;
0475 }
0476 
0477 /* Method compareProtoSegment
0478  *      
0479  * For hit coming from the same layer of an existing hit within the proto segment
0480  * test if achieve better chi^2 by using this hit than the other
0481  *
0482  */
0483 void CSCSegAlgoDF::compareProtoSegment(const CSCRecHit2D* h, int layer) {
0484   //  std::cout << "[CSCSegAlgoDF::compareProtoSegment] for hit on layer " << layer
0485   //        << " with protoSegment.size() = " << protoSegment.size() << std::endl;
0486 
0487   // Try adding the hit to existing segment, and remove old one existing in same layer
0488   ChamberHitContainer::iterator it;
0489   for (it = protoSegment.begin(); it != protoSegment.end();) {
0490     if ((*it)->cscDetId().layer() == layer) {
0491       it = protoSegment.erase(it);
0492     } else {
0493       ++it;
0494     }
0495   }
0496 
0497   //  std::cout << "[CSCSegAlgoDF::compareProtoSegment] about to add hit on layer " << layer
0498   //        << " with protoSegment.size() = " << protoSegment.size() << std::endl;
0499 
0500   bool ok = addHit(h, layer);
0501 
0502   CSCSegFit* newfit = nullptr;
0503   if (ok) {
0504     newfit = new CSCSegFit(theChamber, protoSegment);
0505     //    std::cout << "[CSCSegAlgoDF::compareProtoSegment] newfit = " << newfit << std::endl;
0506     newfit->fit();
0507   }
0508   if (!ok || (newfit->chi2() > sfit_->chi2())) {
0509     //    std::cout << "[CSCSegAlgoDF::compareProtoSegment] about to delete newfit = " << newfit << std::endl;
0510     delete newfit;  // failed to add a hit or new fit is worse
0511   } else {
0512     //    std::cout << "[CSCSegAlgoDF::compareProtoSegment] about to delete sfit_ = " << sfit_ << std::endl;
0513     delete sfit_;  // new fit is better
0514     sfit_ = newfit;
0515     //    std::cout << "[CSCSegAlgoDF::compareProtoSegment] reset sfit_ = " << sfit_ << std::endl;
0516   }
0517 }
0518 
0519 /* Method flagHitsAsUsed
0520  *
0521  * Flag hits which have entered segment building so we don't reuse them.
0522  * Also flag does which were very close to segment to reduce combinatorics
0523  */
0524 void CSCSegAlgoDF::flagHitsAsUsed(const ChamberHitContainer& rechitsInChamber) {
0525   // Flag hits on segment as used
0526   ChamberHitContainerCIt ib = rechitsInChamber.begin();
0527   ChamberHitContainerCIt hi, iu;
0528 
0529   for (hi = protoSegment.begin(); hi != protoSegment.end(); ++hi) {
0530     for (iu = ib; iu != rechitsInChamber.end(); ++iu) {
0531       if (*hi == *iu)
0532         usedHits[iu - ib] = true;
0533     }
0534   }
0535   // Don't reject hits marked as "nearby" for now.
0536   // So this is bypassed at all times for now !!!
0537   // Perhaps add back to speed up algorithm some more
0538   if (!closeHits.empty())
0539     return;
0540   for (hi = closeHits.begin(); hi != closeHits.end(); ++hi) {
0541     for (iu = ib; iu != rechitsInChamber.end(); ++iu) {
0542       if (*hi == *iu)
0543         usedHits[iu - ib] = true;
0544     }
0545   }
0546 }
0547 
0548 // Try to clean up segments by quickly looking at residuals
0549 void CSCSegAlgoDF::pruneFromResidual(void) {
0550   // Only prune if have at least 5 hits
0551   if (protoSegment.size() < 5)
0552     return;
0553 
0554   // Now Study residuals
0555 
0556   float maxResidual = 0.;
0557   float sumResidual = 0.;
0558   int nHits = 0;
0559   int badIndex = -1;
0560   int j = 0;
0561 
0562   ChamberHitContainer::const_iterator ih;
0563 
0564   for (ih = protoSegment.begin(); ih != protoSegment.end(); ++ih) {
0565     const CSCRecHit2D& hit = (**ih);
0566     const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer());
0567     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0568     LocalPoint lp = theChamber->toLocal(gp);
0569 
0570     float residual = sfit_->Rdev(lp.x(), lp.y(), lp.z());
0571 
0572     sumResidual += residual;
0573     nHits++;
0574     if (residual > maxResidual) {
0575       maxResidual = residual;
0576       badIndex = j;
0577     }
0578     j++;
0579   }
0580 
0581   float corrAvgResidual = (sumResidual - maxResidual) / (nHits - 1);
0582 
0583   // Keep all hits
0584   if (maxResidual / corrAvgResidual < maxRatioResidual)
0585     return;
0586 
0587   // Drop worse hit and recompute segment properties + fill
0588 
0589   ChamberHitContainer newProtoSegment;
0590 
0591   j = 0;
0592   for (ih = protoSegment.begin(); ih != protoSegment.end(); ++ih) {
0593     if (j != badIndex)
0594       newProtoSegment.push_back(*ih);
0595     j++;
0596   }
0597 
0598   protoSegment.clear();
0599 
0600   for (ih = newProtoSegment.begin(); ih != newProtoSegment.end(); ++ih) {
0601     protoSegment.push_back(*ih);
0602   }
0603 
0604   // Compute segment parameters
0605   updateParameters();
0606 }
0607 
0608 void CSCSegAlgoDF::dumpSegment(const CSCSegment& seg) const {
0609   edm::LogVerbatim("CSCSegment") << "CSCSegment in " << theChamber->id() << "\nlocal position = " << seg.localPosition()
0610                                  << "\nerror = " << seg.localPositionError()
0611                                  << "\nlocal direction = " << seg.localDirection()
0612                                  << "\nerror =" << seg.localDirectionError() << "\ncovariance matrix"
0613                                  << seg.parametersError() << "chi2/ndf = " << seg.chi2() << "/"
0614                                  << seg.degreesOfFreedom() << "\n#rechits = " << seg.specificRecHits().size()
0615                                  << "\ntime = " << seg.time();
0616 }