Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-07 04:37:39

0001 /**
0002  * \file CSCSegAlgoTC.cc
0003  *
0004  * Last update: 13.02.2015
0005  *
0006  */
0007 
0008 #include "CSCSegAlgoTC.h"
0009 #include "CSCSegFit.h"
0010 
0011 #include "DataFormats/CSCRecHit/interface/CSCSegment.h"
0012 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0013 
0014 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0015 
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 
0019 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0020 
0021 #include <algorithm>
0022 #include <cmath>
0023 #include <iostream>
0024 #include <string>
0025 
0026 CSCSegAlgoTC::CSCSegAlgoTC(const edm::ParameterSet& ps)
0027     : CSCSegmentAlgorithm(ps), sfit_(nullptr), myName("CSCSegAlgoTC") {
0028   debugInfo = ps.getUntrackedParameter<bool>("verboseInfo");
0029 
0030   dRPhiMax = ps.getParameter<double>("dRPhiMax");
0031   dPhiMax = ps.getParameter<double>("dPhiMax");
0032   dRPhiFineMax = ps.getParameter<double>("dRPhiFineMax");
0033   dPhiFineMax = ps.getParameter<double>("dPhiFineMax");
0034   chi2Max = ps.getParameter<double>("chi2Max");
0035   chi2ndfProbMin = ps.getParameter<double>("chi2ndfProbMin");
0036   minLayersApart = ps.getParameter<int>("minLayersApart");
0037   SegmentSorting = ps.getParameter<int>("SegmentSorting");
0038 
0039   LogDebug("CSC") << myName << " has algorithm cuts set to: \n"
0040                   << "--------------------------------------------------------------------\n"
0041                   << "dRPhiMax     = " << dRPhiMax << '\n'
0042                   << "dPhiMax      = " << dPhiMax << '\n'
0043                   << "dRPhiFineMax = " << dRPhiFineMax << '\n'
0044                   << "dPhiFineMax  = " << dPhiFineMax << '\n'
0045                   << "chi2Max      = " << chi2Max << '\n'
0046                   << "chi2ndfProbMin = " << chi2ndfProbMin << '\n'
0047                   << "minLayersApart = " << minLayersApart << '\n'
0048                   << "SegmentSorting = " << SegmentSorting << std::endl;
0049 }
0050 
0051 std::vector<CSCSegment> CSCSegAlgoTC::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
0052   theChamber = aChamber;
0053 
0054   return buildSegments(rechits);
0055 }
0056 
0057 std::vector<CSCSegment> CSCSegAlgoTC::buildSegments(const ChamberHitContainer& urechits) {
0058   //  ChamberHitContainer rechits = urechits;
0059   ChamberHitContainer rechits(urechits);
0060 
0061   //  edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoTC::buildSegments] start building segments in " << theChamber->id();
0062   //  edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoTC::buildSegments] size of rechit container = " << rechits.size();
0063 
0064   if (rechits.size() < 2) {
0065     LogDebug("CSC") << myName << ": " << rechits.size() << "     hit(s) in chamber is not enough to build a segment.\n";
0066     return std::vector<CSCSegment>();
0067   }
0068 
0069   LayerIndex layerIndex(rechits.size());
0070 
0071   for (size_t i = 0; i < rechits.size(); ++i) {
0072     short ilay = rechits[i]->cscDetId().layer();
0073     layerIndex[i] = ilay;
0074     //    edm::LogVerbatim("CSCSegment") << "layerIndex[" << i << "] should   be " << rechits[i]->cscDetId().layer();
0075     //    edm::LogVerbatim("CSCSegment") << "layerIndex[" << i << "] actually is " << layerIndex[i];
0076   }
0077 
0078   double z1 = theChamber->layer(1)->position().z();
0079   double z6 = theChamber->layer(6)->position().z();
0080 
0081   if (z1 > 0.) {
0082     if (z1 > z6) {
0083       reverse(layerIndex.begin(), layerIndex.end());
0084       reverse(rechits.begin(), rechits.end());
0085     }
0086   } else if (z1 < 0.) {
0087     if (z1 < z6) {
0088       reverse(layerIndex.begin(), layerIndex.end());
0089       reverse(rechits.begin(), rechits.end());
0090     }
0091   }
0092 
0093   // Dump rechits after sorting?
0094   //  if (debugInfo) dumpHits(rechits);
0095 
0096   //  if (rechits.size() < 2) {
0097   //    LogDebug("CSC") << myName << ": " << rechits.size() <<
0098   //      "  hit(s) in chamber is not enough to build a segment.\n";
0099   //    return std::vector<CSCSegment>();
0100   //  }
0101 
0102   // We have at least 2 hits. We intend to try all possible pairs of hits to start
0103   // segment building. 'All possible' means each hit lies on different layers in the chamber.
0104   // For now we don't care whether a hit has already been allocated to another segment;
0105   // we'll sort that out after building all possible segments.
0106 
0107   // Choose first hit (as close to IP as possible) h1 and
0108   // second hit (as far from IP as possible) h2
0109   // To do this we iterate over hits in the chamber by layer - pick two layers.
0110   // @@ Require the two layers are at least 3 layers apart. May need tuning?
0111   // Then we iterate over hits within each of these layers and pick h1 and h2 from these.
0112   // If they are 'close enough' we build an empty segment.
0113   // Then try adding hits to this segment.
0114 
0115   // Define buffer for segments we build (at the end we'll sort them somehow, and remove
0116   // those which share hits with better-quality segments.
0117 
0118   std::vector<CSCSegment> segments;
0119 
0120   sfit_ = nullptr;
0121 
0122   ChamberHitContainerCIt ib = rechits.begin();
0123   ChamberHitContainerCIt ie = rechits.end();
0124 
0125   for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
0126     int layer1 = layerIndex[i1 - ib];
0127 
0128     const CSCRecHit2D* h1 = *i1;
0129 
0130     for (ChamberHitContainerCIt i2 = ie - 1; i2 != i1; --i2) {
0131       int layer2 = layerIndex[i2 - ib];
0132 
0133       if (abs(layer2 - layer1) < minLayersApart)
0134         break;
0135 
0136       const CSCRecHit2D* h2 = *i2;
0137 
0138       if (areHitsCloseInLocalX(h1, h2) && areHitsCloseInGlobalPhi(h1, h2)) {
0139         proto_segment.clear();
0140 
0141         const CSCLayer* l1 = theChamber->layer(h1->cscDetId().layer());
0142         GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0143         const CSCLayer* l2 = theChamber->layer(h2->cscDetId().layer());
0144         GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0145         LogDebug("CSCSegment") << "start new segment from hits "
0146                                << "h1: " << gp1 << " - h2: " << gp2 << "\n";
0147         //  edm::LogVerbatim("CSCSegment") << "start new segment from hits " << "h1: " << gp1 << " - h2: " << gp2;
0148         //  edm::LogVerbatim("CSCSegment") << "on layers " << layer1 << " and " << layer2;
0149 
0150         if (!addHit(h1, layer1)) {
0151           LogDebug("CSCSegment") << "  failed to add hit h1\n";
0152           continue;
0153         }
0154 
0155         if (!addHit(h2, layer2)) {
0156           LogDebug("CSCSegment") << "  failed to add hit h2\n";
0157           continue;
0158         }
0159 
0160         if (sfit_)
0161           tryAddingHitsToSegment(rechits, i1, i2);  // can only add hits if there's a segment
0162 
0163         // if a segment has been found push back it into the segment collection
0164         if (proto_segment.empty()) {
0165           LogDebug("CSCSegment") << "No segment found.\n";
0166           //      edm::LogVerbatim("CSCSegment") << "No segment found.";
0167         } else {
0168           //@@ THIS IS GOING TO BE TRICKY - CREATING MANY FITS ON HEAP
0169           //@@ BUT MEMBER VARIABLE JUST POINTS TO MOST RECENT ONE
0170           candidates.push_back(sfit_);  // store the current fit
0171           LogDebug("CSCSegment") << "Found a segment.\n";
0172           //      edm::LogVerbatim("CSCSegment") << "Found a segment.";
0173         }
0174       }
0175     }
0176   }
0177 
0178   //  edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoTC::buildSegments] no. of candidates before pruning = " << candidates.size();
0179 
0180   // We've built all possible segments. Now pick the best, non-overlapping subset.
0181   pruneTheSegments(rechits);
0182 
0183   // Create CSCSegments for the surviving candidates
0184   for (unsigned int i = 0; i < candidates.size(); ++i) {
0185     CSCSegFit* sfit = candidates[i];
0186     //    edm::LogVerbatim("CSCSegment") << "candidate fit " << i+1 << " of " << candidates.size() << " is at " << sfit;
0187     //    if ( !sfit ) {
0188     //      edm::LogVerbatim("CSCSegment") << "stored a null pointer for element " << i+1 << " of " << candidates.size();
0189     //      continue;
0190     //    }
0191     CSCSegment temp(sfit->hits(), sfit->intercept(), sfit->localdir(), sfit->covarianceMatrix(), sfit->chi2());
0192     delete sfit;
0193     segments.push_back(temp);
0194     if (debugInfo)
0195       dumpSegment(temp);
0196   }
0197 
0198   // reset member variables
0199   candidates.clear();
0200   sfit_ = nullptr;
0201 
0202   //  edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoTC::buildSegments] no. of segments returned = " << segments.size();
0203 
0204   return segments;
0205 }
0206 
0207 void CSCSegAlgoTC::tryAddingHitsToSegment(const ChamberHitContainer& rechits,
0208                                           const ChamberHitContainerCIt i1,
0209                                           const ChamberHitContainerCIt i2) {
0210   // Iterate over the layers with hits in the chamber
0211   // Skip the layers containing the segment endpoints
0212   // Test each hit on the other layers to see if it is near the segment
0213   // If it is, see whether there is already a hit on the segment from the same layer
0214   //    - if so, and there are more than 2 hits on the segment, copy the segment,
0215   //      replace the old hit with the new hit. If the new segment chi2 is better
0216   //      then replace the original segment with the new one (by swap)
0217   //    - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory
0218   //      then replace the original segment with the new one (by swap)
0219 
0220   ChamberHitContainerCIt ib = rechits.begin();
0221   ChamberHitContainerCIt ie = rechits.end();
0222 
0223   for (ChamberHitContainerCIt i = ib; i != ie; ++i) {
0224     if (i == i1 || i == i2)
0225       continue;
0226 
0227     int layer = (*i)->cscDetId().layer();
0228     const CSCRecHit2D* h = *i;
0229 
0230     if (isHitNearSegment(h)) {
0231       const CSCLayer* l1 = theChamber->layer(h->cscDetId().layer());
0232       GlobalPoint gp1 = l1->toGlobal(h->localPosition());
0233       LogDebug("CSC") << "    hit at global " << gp1 << " is near segment\n.";
0234 
0235       if (hasHitOnLayer(layer)) {
0236         if (proto_segment.size() <= 2) {
0237           LogDebug("CSC") << "    " << proto_segment.size() << " hits on segment...skip hit on same layer.\n";
0238           continue;
0239         }
0240 
0241         compareProtoSegment(h, layer);
0242       } else
0243         increaseProtoSegment(h, layer);
0244     }  // h & seg close
0245   }  // i
0246 }
0247 
0248 bool CSCSegAlgoTC::addHit(const CSCRecHit2D* aHit, int layer) {
0249   // Return true if hit was added successfully
0250   // (and then parameters are updated).
0251   // Return false if there is already a hit on the same layer, or insert failed.
0252 
0253   bool ok = true;
0254   ChamberHitContainer::const_iterator it;
0255 
0256   for (it = proto_segment.begin(); it != proto_segment.end(); it++)
0257     if (((*it)->cscDetId().layer() == layer) && (aHit != *it))
0258       return false;
0259 
0260   if (ok) {
0261     proto_segment.push_back(aHit);
0262     updateParameters();
0263   }
0264   return ok;
0265 }
0266 
0267 bool CSCSegAlgoTC::replaceHit(const CSCRecHit2D* h, int layer) {
0268   // replace a hit from a layer
0269   ChamberHitContainer::iterator it;
0270   for (it = proto_segment.begin(); it != proto_segment.end();) {
0271     if ((*it)->cscDetId().layer() == layer)
0272       it = proto_segment.erase(it);
0273     else
0274       ++it;
0275   }
0276 
0277   return addHit(h, layer);
0278 }
0279 
0280 void CSCSegAlgoTC::updateParameters(void) {
0281   //@@ DO NOT DELETE EXISTING FIT SINCE WE SAVE IT!!
0282   //  delete sfit_;
0283   sfit_ = new CSCSegFit(theChamber, proto_segment);
0284   sfit_->fit();
0285 }
0286 
0287 float CSCSegAlgoTC::phiAtZ(float z) const {
0288   // Returns a phi in [ 0, 2*pi )
0289   const CSCLayer* l1 = theChamber->layer(1);
0290   GlobalPoint gp = l1->toGlobal(sfit_->intercept());
0291   GlobalVector gv = l1->toGlobal(sfit_->localdir());
0292 
0293   float x = gp.x() + (gv.x() / gv.z()) * (z - gp.z());
0294   float y = gp.y() + (gv.y() / gv.z()) * (z - gp.z());
0295   float phi = atan2(y, x);
0296   if (phi < 0.f)
0297     phi += 2. * M_PI;
0298 
0299   return phi;
0300 }
0301 
0302 void CSCSegAlgoTC::compareProtoSegment(const CSCRecHit2D* h, int layer) {
0303   // Save copy of current fit
0304   CSCSegFit* oldfit = new CSCSegFit(*sfit_);
0305 
0306   bool ok = replaceHit(h, layer);  // possible new fit
0307 
0308   if (ok) {
0309     LogDebug("CSC") << "    hit in same layer as a hit on segment; try replacing old one..."
0310                     << " chi2 new: " << sfit_->chi2() << " old: " << oldfit->chi2() << "\n";
0311   }
0312 
0313   if (ok && (sfit_->chi2() < oldfit->chi2())) {
0314     LogDebug("CSC") << "    segment with replaced hit is better.\n";
0315     delete oldfit;  // new fit is better
0316   } else {
0317     delete sfit_;    // new fit is worse
0318     sfit_ = oldfit;  // restore original fit
0319   }
0320 }
0321 
0322 void CSCSegAlgoTC::increaseProtoSegment(const CSCRecHit2D* h, int layer) {
0323   // save copy of input fit
0324   CSCSegFit* oldfit = new CSCSegFit(*sfit_);
0325 
0326   bool ok = addHit(h, layer);  // possible new fit
0327 
0328   if (ok) {
0329     LogDebug("CSC") << "    hit in new layer: added to segment, new chi2: " << sfit_->chi2() << "\n";
0330   }
0331 
0332   //  int ndf = 2*proto_segment.size() - 4;
0333 
0334   if (ok && ((sfit_->chi2() <= 0) || (sfit_->chi2() / sfit_->ndof() < chi2Max))) {
0335     LogDebug("CSC") << "    segment with added hit is good.\n";
0336     delete oldfit;  // new fit is better
0337   } else {
0338     delete sfit_;    // new fit is worse
0339     sfit_ = oldfit;  // restore original fit
0340   }
0341 }
0342 
0343 bool CSCSegAlgoTC::areHitsCloseInLocalX(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const {
0344   float deltaX = (h1->localPosition() - h2->localPosition()).x();
0345   LogDebug("CSC") << "    Hits at local x= " << h1->localPosition().x() << ", " << h2->localPosition().x()
0346                   << " have separation= " << deltaX;
0347   return (fabs(deltaX) < (dRPhiMax)) ? true : false;  // +v
0348 }
0349 
0350 bool CSCSegAlgoTC::areHitsCloseInGlobalPhi(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const {
0351   const CSCLayer* l1 = theChamber->layer(h1->cscDetId().layer());
0352   GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0353   const CSCLayer* l2 = theChamber->layer(h2->cscDetId().layer());
0354   GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0355 
0356   float h1p = gp1.phi();
0357   float h2p = gp2.phi();
0358   float dphi12 = h1p - h2p;
0359 
0360   // Into range [-pi, pi) (phi() returns values in this range)
0361   if (dphi12 < -M_PI)
0362     dphi12 += 2. * M_PI;
0363   if (dphi12 > M_PI)
0364     dphi12 -= 2. * M_PI;
0365   LogDebug("CSC") << "    Hits at global phi= " << h1p << ", " << h2p << " have separation= " << dphi12;
0366   return (fabs(dphi12) < dPhiMax) ? true : false;  // +v
0367 }
0368 
0369 bool CSCSegAlgoTC::isHitNearSegment(const CSCRecHit2D* h) const {
0370   // Is hit near segment?
0371   // Requires deltaphi and rxy*deltaphi within ranges specified
0372   // in orcarc, or by default, where rxy=sqrt(x**2+y**2) of hit itself.
0373   // Note that to make intuitive cuts on delta(phi) one must work in
0374   // phi range (-pi, +pi] not [0, 2pi
0375 
0376   const CSCLayer* l1 = theChamber->layer(h->cscDetId().layer());
0377   GlobalPoint hp = l1->toGlobal(h->localPosition());
0378 
0379   float hphi = hp.phi();  // in (-pi, +pi]
0380   if (hphi < 0.)
0381     hphi += 2. * M_PI;          // into range [0, 2pi)
0382   float sphi = phiAtZ(hp.z());  // in [0, 2*pi)
0383   float phidif = sphi - hphi;
0384   if (phidif < 0.)
0385     phidif += 2. * M_PI;  // into range [0, 2pi)
0386   if (phidif > M_PI)
0387     phidif -= 2. * M_PI;  // into range (-pi, pi]
0388 
0389   float dRPhi = fabs(phidif) * hp.perp();
0390   LogDebug("CSC") << "    is hit at phi_h= " << hphi << " near segment phi_seg= " << sphi << "? is " << dRPhi << "<"
0391                   << dRPhiFineMax << " ? "
0392                   << " and is |" << phidif << "|<" << dPhiFineMax << " ?";
0393 
0394   return ((dRPhi < dRPhiFineMax) && (fabs(phidif) < dPhiFineMax)) ? true : false;  // +v
0395 }
0396 
0397 bool CSCSegAlgoTC::hasHitOnLayer(int layer) const {
0398   // Is there is already a hit on this layer?
0399   ChamberHitContainerCIt it;
0400 
0401   for (it = proto_segment.begin(); it != proto_segment.end(); it++)
0402     if ((*it)->cscDetId().layer() == layer)
0403       return true;
0404 
0405   return false;
0406 }
0407 
0408 void CSCSegAlgoTC::dumpHits(const ChamberHitContainer& rechits) const {
0409   // Dump positions of RecHit's in each CSCChamber
0410   ChamberHitContainerCIt it;
0411 
0412   for (it = rechits.begin(); it != rechits.end(); it++) {
0413     const CSCLayer* l1 = theChamber->layer((*it)->cscDetId().layer());
0414     GlobalPoint gp1 = l1->toGlobal((*it)->localPosition());
0415 
0416     LogDebug("CSC") << "Global pos.: " << gp1 << ", phi: " << gp1.phi()
0417                     << ". Local position: " << (*it)->localPosition() << ", phi: " << (*it)->localPosition().phi()
0418                     << ". Layer: " << (*it)->cscDetId().layer() << "\n";
0419   }
0420 }
0421 
0422 bool CSCSegAlgoTC::isSegmentGood(std::vector<CSCSegFit*>::iterator seg,
0423                                  const ChamberHitContainer& rechitsInChamber,
0424                                  BoolContainer& used) const {
0425   // Apply any selection cuts to segment
0426 
0427   // 1) Require a minimum no. of hits
0428   //   (@@ THESE VALUES SHOULD BECOME PARAMETERS?)
0429 
0430   // 2) Ensure no hits on segment are already assigned to another segment
0431   //    (typically of higher quality)
0432 
0433   size_t iadd = (rechitsInChamber.size() > 20) ? 1 : 0;
0434 
0435   size_t nhits = (*seg)->nhits();
0436 
0437   if (nhits < 3 + iadd)
0438     return false;
0439 
0440   // Additional part of alternative segment selection scheme: reject
0441   // segments with a chi2 probability of less than chi2ndfProbMin. Relies on list
0442   // being sorted with "SegmentSorting == 2", that is first by nrechits and then
0443   // by chi2prob in subgroups of same no. of rechits.
0444 
0445   if (SegmentSorting == 2) {
0446     double chi2t = (*seg)->chi2();
0447     double ndoft = 2 * nhits - 4;
0448     if (chi2t > 0 && ndoft > 0) {
0449       if (ChiSquaredProbability(chi2t, ndoft) < chi2ndfProbMin) {
0450         return false;
0451       }
0452     } else {
0453       return false;
0454     }
0455   }
0456 
0457   ChamberHitContainer hits_ = (*seg)->hits();
0458 
0459   for (size_t ish = 0; ish < nhits; ++ish) {
0460     ChamberHitContainerCIt ib = rechitsInChamber.begin();
0461 
0462     for (ChamberHitContainerCIt ic = ib; ic != rechitsInChamber.end(); ++ic) {
0463       if ((hits_[ish] == (*ic)) && used[ic - ib])
0464         return false;
0465     }
0466   }
0467 
0468   return true;
0469 }
0470 
0471 void CSCSegAlgoTC::flagHitsAsUsed(std::vector<CSCSegFit*>::iterator seg,
0472                                   const ChamberHitContainer& rechitsInChamber,
0473                                   BoolContainer& used) const {
0474   // Flag hits on segment as used
0475 
0476   ChamberHitContainerCIt ib = rechitsInChamber.begin();
0477   ChamberHitContainer hits = (*seg)->hits();
0478 
0479   for (size_t ish = 0; ish < hits.size(); ++ish) {
0480     for (ChamberHitContainerCIt iu = ib; iu != rechitsInChamber.end(); ++iu)
0481       if (hits[ish] == (*iu))
0482         used[iu - ib] = true;
0483   }
0484 }
0485 
0486 void CSCSegAlgoTC::pruneTheSegments(const ChamberHitContainer& rechitsInChamber) {
0487   // Sort the segment store according to segment 'quality' (chi2/#hits ?) and
0488   // remove any segments which contain hits assigned to higher-quality segments.
0489 
0490   if (candidates.empty())
0491     return;
0492 
0493   // Initialize flags that a given hit has been allocated to a segment
0494   BoolContainer used(rechitsInChamber.size(), false);
0495 
0496   // Sort by chi2/#hits
0497   segmentSort();
0498 
0499   // Select best quality segments, requiring hits are assigned to just one segment
0500   // Want to erase the bad segments, so the iterator must be incremented
0501   // inside the loop, and only when the erase is not called
0502 
0503   std::vector<CSCSegFit*>::iterator is;
0504 
0505   for (is = candidates.begin(); is != candidates.end();) {
0506     bool goodSegment = isSegmentGood(is, rechitsInChamber, used);
0507 
0508     if (goodSegment) {
0509       LogDebug("CSC") << "Accepting segment: ";
0510       flagHitsAsUsed(is, rechitsInChamber, used);
0511       ++is;
0512     } else {
0513       LogDebug("CSC") << "Rejecting segment: ";
0514       delete *is;                 // delete the CSCSegFit*
0515       is = candidates.erase(is);  // erase the element in container
0516     }
0517   }
0518 }
0519 
0520 void CSCSegAlgoTC::segmentSort() {
0521   // The segment collection is sorted according to e.g. chi2/#hits
0522 
0523   for (size_t i = 0; i < candidates.size() - 1; ++i) {
0524     for (size_t j = i + 1; j < candidates.size(); ++j) {
0525       size_t ni = (candidates[i]->hits()).size();
0526       size_t nj = (candidates[j]->hits()).size();
0527 
0528       // Sort criterion: first sort by no. of rechits, then in groups of rechits by chi2prob
0529       if (SegmentSorting == 2) {
0530         if (nj > ni) {  // sort by no. of rechits
0531           CSCSegFit* temp = candidates[j];
0532           candidates[j] = candidates[i];
0533           candidates[i] = temp;
0534         }
0535         // sort by chi2 probability in subgroups with equal nr of rechits
0536         //  if(chi2s[i] != 0. && 2*n2-4 > 0 ) {
0537         if (candidates[i]->chi2() > 0 && candidates[i]->ndof() > 0) {
0538           if (nj == ni && (ChiSquaredProbability(candidates[i]->chi2(), (double)(candidates[i]->ndof())) <
0539                            ChiSquaredProbability(candidates[j]->chi2(), (double)(candidates[j]->ndof())))) {
0540             CSCSegFit* temp = candidates[j];
0541             candidates[j] = candidates[i];
0542             candidates[i] = temp;
0543           }
0544         }
0545       } else if (SegmentSorting == 1) {
0546         if ((candidates[i]->chi2() / ni) > (candidates[j]->chi2() / nj)) {
0547           CSCSegFit* temp = candidates[j];
0548           candidates[j] = candidates[i];
0549           candidates[i] = temp;
0550         }
0551       } else {
0552         LogDebug("CSC") << "No valid segment sorting specified. Algorithm misconfigured! \n";
0553       }
0554     }
0555   }
0556 }
0557 
0558 void CSCSegAlgoTC::dumpSegment(const CSCSegment& seg) const {
0559   edm::LogVerbatim("CSCSegment") << "CSCSegment in " << theChamber->id() << "\nlocal position = " << seg.localPosition()
0560                                  << "\nerror = " << seg.localPositionError()
0561                                  << "\nlocal direction = " << seg.localDirection()
0562                                  << "\nerror =" << seg.localDirectionError() << "\ncovariance matrix"
0563                                  << seg.parametersError() << "chi2/ndf = " << seg.chi2() << "/"
0564                                  << seg.degreesOfFreedom() << "\n#rechits = " << seg.specificRecHits().size()
0565                                  << "\ntime = " << seg.time();
0566 }