CSCSegAlgoDF.cc

Go to the documentation of this file.

#include "CSCSegAlgoDF.h"
#include "CSCSegFit.h"

#include "Geometry/CSCGeometry/interface/CSCLayer.h"

#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"

#include "CSCSegAlgoPreClustering.h"
#include "CSCSegAlgoShowering.h"

#include <algorithm>
#include <cmath>
#include <iostream>
#include <string>

/* Constructor
 *
 */
CSCSegAlgoDF::CSCSegAlgoDF(const edm::ParameterSet& ps)
    : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoDF"), sfit_(nullptr) {
  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);
}

/* Destructor
 *
 */
CSCSegAlgoDF::~CSCSegAlgoDF() {
  delete preCluster_;
  delete showering_;
}

/* run
 *
 */
std::vector<CSCSegment> CSCSegAlgoDF::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
  // 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;
}

/* 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 loop.
 *
 * 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<CSCSegment> CSCSegAlgoDF::buildSegments(const ChamberHitContainer& _rechits) {
  ChamberHitContainer rechits = _rechits;
  // Clear buffer for segment vector
  std::vector<CSCSegment> segmentInChamber;
  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);

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

  //  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();
  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());
    }
  }

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

    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();

      // We want hit wrt chamber (and local z will be != 0)
      float dz = gp2.z() - gp1.z();
      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(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);

      // 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 (sfit_->chi2() > chi2Max)
        segok = false;

      if (segok) {
        // 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_ = nullptr;  // 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
        flagHitsAsUsed(rechits);
      }
    }
  }
  // Return #7
  return segmentInChamber;
}

/* Method tryAddingHitsToSegment
 *
 * Look at left over hits and try to add them to proto segment by looking how far they
 * are from the segment in terms of the hit error matrix (so how many sigmas away).
 *
 */
void CSCSegAlgoDF::tryAddingHitsToSegment(const ChamberHitContainer& rechits,
                                          const ChamberHitContainerCIt i1,
                                          const ChamberHitContainerCIt i2,
                                          const LayerIndex& layerIndex) {
  /* 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
 */

  //  std::cout << "[CSCSegAlgoDF::tryAddingHitsToSegment] entering"
  //        << " with rechits.size() = " << rechits.size() << std::endl;

  ChamberHitContainerCIt ib = rechits.begin();
  ChamberHitContainerCIt ie = rechits.end();
  closeHits.clear();

  //  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);
        } else {
          closeHits.push_back(h);
        }
      }

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

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

/* isHitNearSegment
 *
 * Compare rechit with expected position from protosegment
 */
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());
  float Hphi = Hgp.phi();
  if (Hphi < 0.)
    Hphi += 2. * M_PI;
  LocalPoint Hlp = theChamber->toLocal(Hgp);
  float z = Hlp.z();

  float LocalX = sfit_->xfit(z);
  float LocalY = sfit_->yfit(z);

  LocalPoint Slp(LocalX, LocalY, z);
  GlobalPoint Sgp = theChamber->toGlobal(Slp);
  float Sphi = Sgp.phi();
  if (Sphi < 0.)
    Sphi += 2. * M_PI;
  float R = sqrt(Sgp.x() * Sgp.x() + Sgp.y() * Sgp.y());

  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;

  float RdeltaPhi = R * deltaPhi;

  if (RdeltaPhi < dRPhiFineMax && deltaPhi < dPhiFineMax)
    return true;

  return false;
}

/* Method addHit
 *
 * Test if can add hit to proto segment. If so, try to add it.
 *
 */
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.

  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)
    if (aHit == (*it))
      return false;

  protoSegment.push_back(aHit);

  return true;
}

/* Method updateParameters
 *      
 * Perform a simple Least Square Fit on proto segment to determine slope and intercept
 *
 */
void CSCSegAlgoDF::updateParameters() {
  // 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;

  return false;
}

/* Method compareProtoSegment
 *      
 * For hit coming from the same layer of an existing hit within the proto segment
 * test if achieve better chi^2 by using this hit than the other
 *
 */
void CSCSegAlgoDF::compareProtoSegment(const CSCRecHit2D* h, int layer) {
  //  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;
  for (it = protoSegment.begin(); it != protoSegment.end();) {
    if ((*it)->cscDetId().layer() == layer) {
      it = protoSegment.erase(it);
    } else {
      ++it;
    }
  }

  //  std::cout << "[CSCSegAlgoDF::compareProtoSegment] about to add hit on layer " << layer
  //        << " with protoSegment.size() = " << protoSegment.size() << std::endl;

  bool ok = addHit(h, layer);

  CSCSegFit* newfit = nullptr;
  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;
  }
}

/* Method flagHitsAsUsed
 *
 * Flag hits which have entered segment building so we don't reuse them.
 * Also flag does which were very close to segment to reduce combinatorics
 */
void CSCSegAlgoDF::flagHitsAsUsed(const ChamberHitContainer& rechitsInChamber) {
  // 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.empty())
    return;
  for (hi = closeHits.begin(); hi != closeHits.end(); ++hi) {
    for (iu = ib; iu != rechitsInChamber.end(); ++iu) {
      if (*hi == *iu)
        usedHits[iu - ib] = true;
    }
  }
}

// Try to clean up segments by quickly looking at residuals
void CSCSegAlgoDF::pruneFromResidual(void) {
  // 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);

    float residual = sfit_->Rdev(lp.x(), lp.y(), lp.z());

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

  // Compute segment parameters
  updateParameters();
}

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