DTHitPairForFit.cc

Go to the documentation of this file.

/* This Class Header */
#include "RecoLocalMuon/DTSegment/src/DTHitPairForFit.h"

/* Collaborating Class Header */
#include "FWCore/Utilities/interface/Exception.h"

/* C++ Headers */
#include <iostream>
using namespace std;

/* ====================================================================== */

DTHitPairForFit::DTHitPairForFit(const DTRecHit1DPair& pair,
                                 const DTSuperLayer& sl,
                                 const edm::ESHandle<DTGeometry>& dtGeom) {
  theWireId = pair.wireId();
  theDigiTime = pair.digiTime();

  const DTLayer* layer = dtGeom->layer(theWireId.layerId());

  // transform the Local position in Layer-rf in a SL local position
  theLeftPos = sl.toLocal(layer->toGlobal(pair.componentRecHit(DTEnums::Left)->localPosition()));
  theRightPos = sl.toLocal(layer->toGlobal(pair.componentRecHit(DTEnums::Right)->localPosition()));

  // TODO how do I transform an error from local to global?
  theError = pair.componentRecHit(DTEnums::Left)->localPositionError();
  // theError =
  //   layer->surface().toLocal(sl.surface().toGlobal(pair.componentRecHit(DTEnums::Left)->localPositionError()));
}

DTHitPairForFit::~DTHitPairForFit() {}

/* Operations */
LocalPoint DTHitPairForFit::localPosition(DTEnums::DTCellSide s) const {
  if (s == DTEnums::Left)
    return theLeftPos;
  else if (s == DTEnums::Right)
    return theRightPos;
  else {
    throw cms::Exception("DTHitPairForFit") << " localPosition called with undef LR code" << endl;
    return LocalPoint();
  }
}

pair<bool, bool> DTHitPairForFit::isCompatible(const LocalPoint& posIni, const LocalVector& dirIni) const {
  pair<bool, bool> ret;
  LocalPoint segPosAtZLeft = posIni + dirIni * (theLeftPos.z() - posIni.z()) / dirIni.z();
  LocalPoint segPosAtZRight = posIni + dirIni * (theRightPos.z() - posIni.z()) / dirIni.z();
  float dxLeft = fabs(theLeftPos.x() - segPosAtZLeft.x());
  float dxRight = fabs(theRightPos.x() - segPosAtZRight.x());
  float exx = sqrt(theError.xx());
  // if both below 3 sigma, return both
  // if both at 10 sigma or above, return none
  // if one is below N sigma and one above, for 10>=N>=3, match only that one, otherwise none
  if (std::max(dxLeft, dxRight) < 3 * exx) {
    ret = make_pair(true, true);
  } else if (std::min(dxLeft, dxRight) >= 10 * exx) {
    ret = make_pair(false, false);
  } else {
    float sigmasL = floorf(dxLeft / exx), sigmasR = floorf(dxRight / exx);
    ret.first = (sigmasL < sigmasR);
    ret.second = (sigmasR < sigmasL);
  }
  return ret;
}

bool DTHitPairForFit::operator<(const DTHitPairForFit& hit) const {
  //SL if same layer use x() for strict ordering
  if (id() == hit.id())
    return (theLeftPos.x() < hit.localPosition(DTEnums::Left).x());
  return (id() < hit.id());
}

bool DTHitPairForFit::operator==(const DTHitPairForFit& hit) const {
  return (id() == hit.id() && fabs(digiTime() - hit.digiTime()) < 0.1);
}

ostream& operator<<(ostream& out, const DTHitPairForFit& hit) {
  out << hit.leftPos() << " " << hit.rightPos();
  return out;
}