CSCComparatorDigiFitter.cc

Go to the documentation of this file.

#include "L1Trigger/CSCTriggerPrimitives/interface/CSCComparatorDigiFitter.h"
#include "L1Trigger/CSCCommonTrigger/interface/CSCConstants.h"
#include "DataFormats/Math/interface/normalizedPhi.h"

namespace {
  // CSC LCT patterns
  // the number quotes the distance to the center

  // pid=0: no pattern found
  const std::vector<std::vector<int> > pat0delta(CSCConstants::NUM_LAYERS);

  // pid=1: layer-OR trigger
  const std::vector<std::vector<int> > pat1delta{{-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5},
                                                 {-2, -1, 0, 1, 2},
                                                 {0},
                                                 {-2, -1, 0, 1, 2},
                                                 {-4, -3, -2, -1, 0, 1, 2, 3, 4},
                                                 {-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5}};

  // pid=2: right-bending (large)
  const std::vector<std::vector<int> > pat2delta{{3, 4, 5}, {1, 2}, {0}, {-2, -1, 0}, {-4, -3, -2}, {-5, -4, -3}};

  // pid=3: left-bending (large)
  const std::vector<std::vector<int> > pat3delta{{-5, -4, -3}, {-2, -1}, {0}, {0, 1, 2}, {2, 3, 4}, {3, 4, 5}};

  // pid=4: right-bending (medium)
  const std::vector<std::vector<int> > pat4delta{{2, 3, 4}, {1, 2}, {0}, {-2, -1}, {-4, -3, -2}, {-4, -3, -2}

  };

  // pid=5: left-bending (medium)
  const std::vector<std::vector<int> > pat5delta{{-4, -3, -2}, {-2, -1}, {0}, {1, 2}, {2, 3, 4}, {2, 3, 4}

  };

  // pid=6: right-bending (medium)
  const std::vector<std::vector<int> > pat6delta{{1, 2, 3}, {0, 1}, {0}, {-1, 0}, {-2, -1}, {-3, -2, -1}};

  // pid=7: left-bending (medium)
  const std::vector<std::vector<int> > pat7delta{{-3, -2, -1}, {-1, 0}, {0}, {0, 1}, {1, 2}, {1, 2, 3}};

  // pid=8: right-bending (small)
  const std::vector<std::vector<int> > pat8delta{{0, 1, 2}, {0, 1}, {0}, {-1, 0}, {-2, -1, 0}, {-2, -1, 0}};

  // pid=9: left-bending (small)
  const std::vector<std::vector<int> > pat9delta{{-2, -1, 0}, {-1, 0}, {0}, {0, 1}, {0, 1, 2}, {0, 1, 2}};

  // pid=A: straight-through
  const std::vector<std::vector<int> > patAdelta{{-1, 0, 1}, {0}, {0}, {0}, {-1, 0, 1}, {-1, 0, 1}};

  const std::vector<std::vector<std::vector<int> > > patIndexToPatternDelta{pat0delta,
                                                                            pat1delta,
                                                                            pat2delta,
                                                                            pat3delta,
                                                                            pat4delta,
                                                                            pat5delta,
                                                                            pat6delta,
                                                                            pat7delta,
                                                                            pat8delta,
                                                                            pat9delta,
                                                                            patAdelta};
}  // namespace

void CSCComparatorDigiFitter::matchingComparatorDigisLCT(const CSCDetId& ch_id,
                                                         const CSCCorrelatedLCTDigi& stub,
                                                         const CSCComparatorDigiCollection& hCSCComparators) {
  // fetch the CSC comparator digis in this chamber
  for (int iLayer = 1; iLayer <= CSCConstants::NUM_LAYERS; ++iLayer) {
    const CSCDetId layerId(ch_id.endcap(), ch_id.station(), ch_id.ring(), ch_id.chamber(), iLayer);

    // get the digis per layer
    const auto& compRange = hCSCComparators.get(layerId);
    CSCComparatorDigiContainer compDigis;

    for (auto compDigiItr = compRange.first; compDigiItr != compRange.second; compDigiItr++) {
      const auto& compDigi = *compDigiItr;

      //if (stub.getTimeBin() < 4 or stub.getTimeBin() > 8) continue;
      const int stubHalfStrip(compDigi.getHalfStrip());

      // these comparator digis never fit the pattern anyway!
      if (std::abs(stubHalfStrip - stub.getStrip()) > 5)
        continue;

      // check if this comparator digi fits the pattern
      if (comparatorInLCTPattern(stub.getStrip(), stub.getPattern(), iLayer, stubHalfStrip)) {
        compDigis.push_back(compDigi);
      }
    }
    compDigisIds_.emplace_back(layerId, compDigis);
  }
}

void CSCComparatorDigiFitter::getComparatorDigiCoordinates(const CSCDetId& ch_id, const CSCCorrelatedLCTDigi& stub) {
  const auto& cscChamber = cscGeometry_->chamber(ch_id);

  // get the z and phi positions of the comparator digis
  float radius_ = 0.0;

  // loop on all matching digis
  for (const auto& p : compDigisIds_) {
    const auto& detId = p.first;

    float phi_tmp = 0.0;
    float radius_tmp = 0.0;
    float z_tmp = 0.0;

    // ignore layers with no digis
    if (p.second.empty())
      continue;

    // loop on all matching digis in this layer
    for (const auto& hit : p.second) {
      const float fractional_strip = hit.getFractionalStrip();
      const auto& layer_geo = cscChamber->layer(detId.layer())->geometry();
      const float wire = layer_geo->middleWireOfGroup(stub.getKeyWG() + 1);

      // get the phi of each comparator digi
      const LocalPoint& csc_intersect = layer_geo->intersectionOfStripAndWire(fractional_strip, wire);
      const GlobalPoint& csc_gp = cscGeometry_->idToDet(detId)->surface().toGlobal(csc_intersect);
      const float gpphi = csc_gp.phi();

      // normalize phi values according to first one
      if (!phis_.empty() and gpphi > 0 and phis_[0] < 0 and (gpphi - phis_[0]) > M_PI)
        phi_tmp += (gpphi - 2 * M_PI);
      else if (!phis_.empty() and gpphi < 0 and phis_[0] > 0 and (gpphi - phis_[0]) < -M_PI)
        phi_tmp += (gpphi + 2 * M_PI);
      else
        phi_tmp += (csc_gp.phi());

      z_tmp = csc_gp.z();
      radius_tmp += csc_gp.perp();
    }

    //in case there are more than one comparator digis in one layer
    radius_tmp = radius_tmp / (p.second).size();
    radius_ += radius_tmp;

    zs_.push_back(z_tmp);
    ezs_.push_back(0);

    phi_tmp = phi_tmp / (p.second).size();
    phis_.push_back(phi_tmp);
    // assume that the charge is flat distributed across the half-strip
    // this is only approximately valid, but good enough for now
    ephis_.push_back(cscHalfStripWidth(detId) / sqrt(12));
  }
}

void CSCComparatorDigiFitter::fit(const CSCDetId& ch_id,
                                  const CSCCorrelatedLCTDigi& stub,
                                  const CSCComparatorDigiCollection& hCSCComparators,
                                  std::vector<float>& fit_phi_layers,
                                  std::vector<float>& fit_z_layers,
                                  float& keyRadius) {
  // clear fit results
  fit_phi_layers.clear();
  fit_z_layers.clear();
  keyRadius = 0;

  // first, match the comparator digis to the LCT
  matchingComparatorDigisLCT(ch_id, stub, hCSCComparators);

  // second, get the coordinates
  getComparatorDigiCoordinates(ch_id, stub);

  // get radius of the stub from key layer
  const CSCDetId key_id(ch_id.endcap(), ch_id.station(), ch_id.ring(), ch_id.chamber(), CSCConstants::KEY_CLCT_LAYER);
  const float fractional_strip = stub.getFractionalStrip();

  const auto& cscChamber = cscGeometry_->chamber(ch_id);
  const auto& layer_geo = cscChamber->layer(CSCConstants::KEY_CLCT_LAYER)->geometry();

  // LCT::getKeyWG() also starts from 0
  const float wire = layer_geo->middleWireOfGroup(stub.getKeyWG() + 1);
  const LocalPoint& csc_intersect = layer_geo->intersectionOfStripAndWire(fractional_strip, wire);
  const GlobalPoint& csc_gp = cscGeometry_->idToDet(key_id)->surface().toGlobal(csc_intersect);

  // get radius from key layer
  if (useKeyRadius_)
    radius_ = radius_ / phis_.size();
  else
    radius_ = csc_gp.perp();

  float alpha = -99., beta = 0.;
  // do a fit to the comparator digis
  calculateSlopeIntercept(alpha, beta);
  if (phis_.size() <= 2 or std::abs(alpha) >= 99) {
    alpha = csc_gp.phi();
    beta = 0.0;
  }

  // determine the pitch of the chamber
  // option to discretize the pitch
  float stripPhiPitch = layer_geo->stripPhiPitch();
  if (nStripBits_)
    stripPhiPitch = stripPhiPitch / nStripBits_;

  // get the fit results
  // option to discretize the fitted phi value
  keyRadius = radius_;

  for (int i = 0; i < CSCConstants::NUM_LAYERS; i++) {
    const float fit_z = cscChamber->layer(i + 1)->centerOfStrip(20).z();
    const float fit_phi = normalizedPhi(alpha + beta * fit_z);
    fit_z_layers.push_back(fit_z);
    fit_phi_layers.push_back(fit_phi);
    if (nStripBits_)
      fit_phi_layers.push_back((std::floor(fit_phi / stripPhiPitch) + 0.5) * stripPhiPitch);
  }
}

void CSCComparatorDigiFitter::calculateSlopeIntercept(float& alpha, float& beta) {
  // if there are at least 3 hits in the chamber, do a linear fit to the
  // comparator digi positions with the chi2 method
  if (phis_.size() >= 3) {
    float Sxx = 0, Sxy = 0, Sx = 0, Sy = 0, S = 0;
    for (unsigned i = 0; i < phis_.size(); ++i) {
      float sigma2_inv = 1. / ephis_[i] * ephis_[i];
      Sxx += zs_[i] * zs_[i] * sigma2_inv;
      Sxy += zs_[i] * phis_[i] * sigma2_inv;
      Sx += zs_[i] * zs_[i] * sigma2_inv;
      Sy += phis_[i] * sigma2_inv;
      S += sigma2_inv;
    }
    float delta = S * Sxx - Sx * Sx;
    alpha = (Sxx * Sy - Sx * Sxy) / delta;
    beta = (S * Sxy - Sx * Sy) / delta;
  } else {
    alpha = -99;
    beta = 0.0;
  }
}

float CSCComparatorDigiFitter::cscHalfStripWidth(const CSCDetId& id) const {
  // what is the chamber type?
  int index = id.iChamberType() - 1;

  // calculate the half strip width of this chamber
  return degrees_[index] * M_PI / 180. / (2. * strips_[index]);
}

bool CSCComparatorDigiFitter::comparatorInLCTPattern(int keyStrip, int pattern, int layer, int halfStrip) const {
  // get the (sub)pattern
  const std::vector<int>& subpat = patIndexToPatternDelta[pattern].at(layer - 1);

  // due to comparator digi time extension in the CLCT processor we need to
  // search a bigger region around the key HS. +/-1, 0 should be sufficient
  const int halfStripDelta = halfStrip - keyStrip;
  return (std::find(subpat.begin(), subpat.end(), halfStripDelta + 1) != subpat.end() or
          std::find(subpat.begin(), subpat.end(), halfStripDelta) != subpat.end() or
          std::find(subpat.begin(), subpat.end(), halfStripDelta - 1) != subpat.end());
}