DTSegmentUpdator Class Reference

#include <DTSegmentUpdator.h>

Public Member Functions
void	calculateT0corr (DTRecSegment2D *seg) const
void	calculateT0corr (DTRecSegment4D *seg) const
	DTSegmentUpdator (const edm::ParameterSet &config, edm::ConsumesCollector)
	Constructor. More...
bool	fit (DTSegmentCand *seg, bool allow3par, const bool fitdebug) const
void	fit (DTRecSegment2D *seg, bool allow3par, bool block3par) const
void	fit (DTRecSegment4D *seg, bool allow3par) const
void	setES (const edm::EventSetup &setup)
	set the setup More...
void	update (DTRecSegment4D *seg, const bool calcT0, bool allow3par) const
	recompute hits position and refit the segment4D More...
void	update (DTRecSegment2D *seg, bool allow3par) const
	recompute hits position and refit the segment2D More...
	~DTSegmentUpdator ()
	Destructor. More...

Private Member Functions
void	fit (const std::vector< float > &x, const std::vector< float > &y, const std::vector< int > &lfit, const std::vector< double > &dist, const std::vector< float > &sigy, LocalPoint &pos, LocalVector &dir, float &cminf, float &vminf, AlgebraicSymMatrix &covMat, double &chi2, const bool allow3par=false, const bool block3par=false) const
	interface to LinearFit More...
void	rejectBadHits (DTChamberRecSegment2D *) const
void	updateHits (DTRecSegment2D *seg, GlobalPoint &gpos, GlobalVector &gdir, const int step=2) const

Private Attributes
bool	debug
double	intime_cut
bool	perform_delta_rejecting
double	T0_hit_resolution
std::unique_ptr< DTRecHitBaseAlgo >	theAlgo
std::unique_ptr< DTLinearFit >	theFitter
edm::ESHandle< DTGeometry >	theGeom
const edm::ESGetToken < DTGeometry, MuonGeometryRecord >	theGeomToken
bool	vdrift_4parfit

Detailed Description

Perform linear fit and hits update for DT segments. Update a segment by improving the hits thanks to the refined knowledge of impact angle and position (also along the wire) and perform linear fit on improved hits.

Author

Stefano Lacaprara - INFN Legnaro stefa.nosp@m.no.l.nosp@m.acapr.nosp@m.ara@.nosp@m.pd.in.nosp@m.fn.i.nosp@m.t

Riccardo Bellan - INFN TO ricca.nosp@m.rdo..nosp@m.bella.nosp@m.n@ce.nosp@m.rn.ch

M.Meneghelli - INFN BO marco.nosp@m..men.nosp@m.eghel.nosp@m.li@c.nosp@m.ern.c.nosp@m.h

Definition at line 46 of file DTSegmentUpdator.h.

Constructor & Destructor Documentation

DTSegmentUpdator::DTSegmentUpdator	(	const edm::ParameterSet &	config,
		edm::ConsumesCollector	cc
	)

Constructor.

Definition at line 47 of file DTSegmentUpdator.cc.

References gather_cfg::cout, debug, edm::ConsumesCollector::esConsumes(), edm::ParameterSet::exists(), get, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), intime_cut, perform_delta_rejecting, T0_hit_resolution, theAlgo, theGeomToken, vdrift_4parfit, and ~DTSegmentUpdator().

    : theFitter{std::make_unique<DTLinearFit>()},
      theAlgo{DTRecHitAlgoFactory::get()->create(
          config.getParameter<string>("recAlgo"), config.getParameter<ParameterSet>("recAlgoConfig"), cc)},
      theGeomToken(cc.esConsumes()),
      vdrift_4parfit(config.getParameter<bool>("performT0_vdriftSegCorrection")),
      T0_hit_resolution(config.getParameter<double>("hit_afterT0_resolution")),
      perform_delta_rejecting(config.getParameter<bool>("perform_delta_rejecting")),
      debug(config.getUntrackedParameter<bool>("debug", false)) {
  intime_cut = 20.;
  if (config.exists("intime_cut"))
    intime_cut = config.getParameter<double>("intime_cut");

  if (debug)
    cout << "[DTSegmentUpdator] Constructor called" << endl;
}

DTSegmentUpdator::~DTSegmentUpdator ( )

default

Destructor.

Referenced by DTSegmentUpdator().

Member Function Documentation

void DTSegmentUpdator::calculateT0corr

(

DTRecSegment2D *

seg

)

const

Definition at line 619 of file DTSegmentUpdator.cc.

References a, b, gather_cfg::cout, debug, HLT_FULL_cff::distance, TrackingRecHit::geographicalId(), phase1PixelTopology::layer, DTEnums::Left, DTRecSegment2D::setT0(), DTRecSegment2D::setVdrift(), DTRecSegment2D::specificRecHits(), DTLayer::specificTopology(), theFitter, theGeom, vdrift_4parfit, DTTopology::wirePosition(), x, PV3DBase< T, PVType, FrameType >::x(), hit::x, y, and PV3DBase< T, PVType, FrameType >::z().

Referenced by calculateT0corr(), DTMeantimerPatternReco4D::reconstruct(), DTCombinatorialPatternReco4D::reconstruct(), and update().

                                                                {
  // WARNING: since this method is called both with a 2D and a 2DPhi as argument
  // seg->geographicalId() can be a superLayerId or a chamberId
  if (debug)
    cout << "[DTSegmentUpdator] CalculateT0corr DTRecSegment4D" << endl;

  vector<double> d_drift;
  vector<float> x;
  vector<float> y;
  vector<int> lc;

  vector<DTRecHit1D> hits = seg->specificRecHits();

  DTWireId wireId;
  int nptfit = 0;

  for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
    // I have to get the hits position (the hit is in the layer rf) in SL frame...
    GlobalPoint glbPos = (theGeom->layer(hit->wireId().layerId()))->toGlobal(hit->localPosition());
    LocalPoint pos = (theGeom->idToDet(seg->geographicalId()))->toLocal(glbPos);

    const DTLayer* layer = theGeom->layer(hit->wireId().layerId());
    float xwire = layer->specificTopology().wirePosition(hit->wireId().wire());
    float distance = fabs(hit->localPosition().x() - xwire);

    int ilc = (hit->lrSide() == DTEnums::Left) ? 1 : -1;

    nptfit++;
    x.push_back(pos.z());
    y.push_back(pos.x());
    lc.push_back(ilc);
    d_drift.push_back(distance);

    // cout << " d_drift "<<distance  <<" npt= " <<npt<<endl;
  }

  double chi2fit = 0.;
  float cminf = 0.;
  float vminf = 0.;
  float a, b;

  if (nptfit > 2) {
    //NB chi2fit is normalized
    theFitter->fit4Var(x, y, lc, d_drift, nptfit, a, b, cminf, vminf, chi2fit, vdrift_4parfit, debug);

    double t0cor = -999.;
    if (cminf > -998.)
      t0cor = -cminf / 0.00543;  // in ns

    //cout << "In calculateT0corr: t0 = " << t0cor << endl;
    //cout << "In calculateT0corr: vminf = " << vminf << endl;
    //cout << "In calculateT0corr: cminf = " << cminf << endl;
    //cout << "In calculateT0corr: chi2 = " << chi2fit << endl;

    seg->setT0(t0cor);      // time  and
    seg->setVdrift(vminf);  //  vdrift correction are recorded in the segment
  }
}

void DTSegmentUpdator::calculateT0corr

(

DTRecSegment4D *

seg

)

const

Definition at line 612 of file DTSegmentUpdator.cc.

References calculateT0corr(), DTRecSegment4D::hasPhi(), DTRecSegment4D::hasZed(), DTRecSegment4D::phiSegment(), and DTRecSegment4D::zSegment().

                                                                {
  if (seg->hasPhi())
    calculateT0corr(seg->phiSegment());
  if (seg->hasZed())
    calculateT0corr(seg->zSegment());
}

bool DTSegmentUpdator::fit	(	DTSegmentCand *	seg,
		bool	allow3par,
		const bool	fitdebug
	)		const

do the linear fit on the hits of the segment candidate and update it. Returns false if the segment candidate is not good()

Definition at line 217 of file DTSegmentUpdator.cc.

References HLT_FULL_cff::chi2, gather_cfg::cout, debug, DeadROC_duringRun::dir, HLT_FULL_cff::distance, DTSegmentCand::hits(), mps_fire::i, DTEnums::Left, DTEnums::Right, DTSegmentCand::setChi2(), DTSegmentCand::setCovMatrix(), DTSegmentCand::setDirection(), DTSegmentCand::setPosition(), DTSegmentCand::sett0(), mathSSE::sqrt(), x, PV3DBase< T, PVType, FrameType >::x(), y, and PV3DBase< T, PVType, FrameType >::z().

Referenced by DTCombinatorialExtendedPatternReco::buildPointsCollection(), DTCombinatorialPatternReco::buildPointsCollection(), DTCombinatorialExtendedPatternReco::extendCandidates(), fit(), DTMeantimerPatternReco::fitWithT0(), trackingPlots.Iteration::modules(), DTRefitAndCombineReco4D::refitSuperSegments(), and update().

                                                                                        {
  //  if (debug && fitdebug) cout << "[DTSegmentUpdator] Fit DTRecSegment2D" << endl;
  //  if (!seg->good()) return false;

  //  DTSuperLayerId DTid = (DTSuperLayerId)seg->superLayer()->id();
  //  if (DTid.superlayer()==2)
  //    allow3par = 0;

  vector<float> x;
  vector<float> y;
  vector<float> sigy;
  vector<int> lfit;
  vector<double> dist;
  int i = 0;

  x.reserve(8);
  y.reserve(8);
  sigy.reserve(8);
  lfit.reserve(8);
  dist.reserve(8);

  for (DTSegmentCand::AssPointCont::const_iterator iter = seg->hits().begin(); iter != seg->hits().end(); ++iter) {
    LocalPoint pos = (*iter).first->localPosition((*iter).second);
    float xwire =
        (((*iter).first)->localPosition(DTEnums::Left).x() + ((*iter).first)->localPosition(DTEnums::Right).x()) / 2.;
    float distance = pos.x() - xwire;

    if ((*iter).second == DTEnums::Left)
      lfit.push_back(1);
    else
      lfit.push_back(-1);

    dist.push_back(distance);
    sigy.push_back(sqrt((*iter).first->localPositionError().xx()));
    x.push_back(pos.z());
    y.push_back(pos.x());
    i++;
  }

  LocalPoint pos;
  LocalVector dir;
  AlgebraicSymMatrix covMat(2);
  float cminf = 0.;
  float vminf = 0.;
  double chi2 = 0.;
  double t0_corr = 0.;

  fit(x, y, lfit, dist, sigy, pos, dir, cminf, vminf, covMat, chi2, allow3par);
  if (cminf != 0)
    t0_corr = -cminf / 0.00543;  // convert drift distance to time

  if (debug && fitdebug)
    cout << "  DTcand chi2: " << chi2 << "/" << x.size() << "   t0: " << t0_corr << endl;

  seg->setPosition(pos);
  seg->setDirection(dir);
  seg->sett0(t0_corr);
  seg->setCovMatrix(covMat);
  seg->setChi2(chi2);

  // cout << "pos " << pos << endl;
  // cout << "dir " << dir << endl;
  // cout << "Mat " << covMat << endl;

  return true;
}

void DTSegmentUpdator::fit	(	DTRecSegment2D *	seg,
		bool	allow3par,
		bool	block3par
	)		const

ditto for true segment: since the fit is applied on a true segment, by definition the segment is "good", while it's not the case for just candidates

Definition at line 284 of file DTSegmentUpdator.cc.

References HLT_FULL_cff::chi2, gather_cfg::cout, debug, DeadROC_duringRun::dir, HLT_FULL_cff::distance, fit(), TrackingRecHit::geographicalId(), phase1PixelTopology::layer, DTEnums::Left, DTRecSegment2D::setChi2(), DTRecSegment2D::setCovMatrix(), DTRecSegment2D::setDirection(), DTRecSegment2D::setPosition(), DTRecSegment2D::setT0(), DTRecSegment2D::specificRecHits(), DTLayer::specificTopology(), mathSSE::sqrt(), theGeom, ErrorFrameTransformer::transform(), DTTopology::wirePosition(), x, PV3DBase< T, PVType, FrameType >::x(), hit::x, LocalError::xx(), y, and PV3DBase< T, PVType, FrameType >::z().

Referenced by trackingPlots.Iteration::modules().

                                                                                    {
  if (debug)
    cout << "[DTSegmentUpdator] Fit DTRecSegment2D - 3par: " << allow3par << endl;

  vector<float> x;
  vector<float> y;
  vector<float> sigy;
  vector<int> lfit;
  vector<double> dist;
  x.reserve(8);
  y.reserve(8);
  sigy.reserve(8);
  lfit.reserve(8);
  dist.reserve(8);

  //  DTSuperLayerId DTid = (DTSuperLayerId)seg->geographicalId();
  //  if (DTid.superlayer()==2)
  //    allow3par = 0;

  vector<DTRecHit1D> hits = seg->specificRecHits();
  for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
    // I have to get the hits position (the hit is in the layer rf) in SL frame...
    GlobalPoint glbPos = (theGeom->layer(hit->wireId().layerId()))->toGlobal(hit->localPosition());
    LocalPoint pos = (theGeom->idToDet(seg->geographicalId()))->toLocal(glbPos);
    x.push_back(pos.z());
    y.push_back(pos.x());

    const DTLayer* layer = theGeom->layer(hit->wireId().layerId());
    float xwire = layer->specificTopology().wirePosition(hit->wireId().wire());
    float distance = fabs(hit->localPosition().x() - xwire);
    dist.push_back(distance);

    int ilc = (hit->lrSide() == DTEnums::Left) ? 1 : -1;
    lfit.push_back(ilc);

    // Get local error in SL frame
    //RB: is it right in this way?
    ErrorFrameTransformer tran;
    GlobalError glbErr =
        tran.transform(hit->localPositionError(), (theGeom->layer(hit->wireId().layerId()))->surface());
    LocalError slErr = tran.transform(glbErr, (theGeom->idToDet(seg->geographicalId()))->surface());
    sigy.push_back(sqrt(slErr.xx()));
  }

  LocalPoint pos;
  LocalVector dir;
  AlgebraicSymMatrix covMat(2);
  double chi2 = 0.;
  float cminf = 0.;
  float vminf = 0.;
  double t0_corr = 0.;

  fit(x, y, lfit, dist, sigy, pos, dir, cminf, vminf, covMat, chi2, allow3par, block3par);
  if (cminf != 0)
    t0_corr = -cminf / 0.00543;  // convert drift distance to time

  if (debug)
    cout << "   DTSeg2d chi2: " << chi2 << endl;
  if (debug)
    cout << "   DTSeg2d Fit t0: " << t0_corr << endl;
  // cout << "pos " << segPosition << endl;
  // cout << "dir " << segDirection << endl;
  // cout << "Mat " << mat << endl;

  seg->setPosition(pos);
  seg->setDirection(dir);
  seg->setCovMatrix(covMat);
  seg->setChi2(chi2);
  seg->setT0(t0_corr);
}

void DTSegmentUpdator::fit	(	DTRecSegment4D *	seg,
		bool	allow3par
	)		const

ditto for true segment 4D, the fit is done on either projection and then the 4D direction and position is built. Since the fit is applied on a true segment, by definition the segment is "good", while it's not the case for just candidates

Definition at line 116 of file DTSegmentUpdator.cc.

References DTSuperLayerId::chamberId(), DTRecSegment4D::chamberId(), DTRecSegment2D::chi2(), funct::cos(), gather_cfg::cout, debug, DeadROC_duringRun::dir, fit(), DTRecSegment4D::hasPhi(), DTRecSegment4D::hasZed(), intime_cut, DTRecSegment2D::localDirection(), DTRecSegment2D::localPosition(), DTRecSegment2D::parametersError(), DTRecSegment4D::phiSegment(), DTRecSegment4D::setCovMatrix(), DTRecSegment4D::setCovMatrixForZed(), DTRecSegment4D::setDirection(), DTRecSegment4D::setPosition(), DTSLRecSegment2D::superLayerId(), DTRecSegment2D::t0(), theGeom, PV3DBase< T, PVType, FrameType >::theta(), GeomDet::toGlobal(), GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::x(), y, PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), and DTRecSegment4D::zSegment().

Referenced by trackingPlots.Iteration::modules().

                                                                    {
  if (debug)
    cout << "[DTSegmentUpdator] Fit DTRecSegment4D:" << endl;
  // after the update must refit the segments

  if (debug) {
    if (seg->hasPhi())
      cout << "    4D Segment contains a Phi segment. t0= " << seg->phiSegment()->t0()
           << "  chi2= " << seg->phiSegment()->chi2() << endl;
    if (seg->hasZed())
      cout << "    4D Segment contains a Zed segment. t0= " << seg->zSegment()->t0()
           << "  chi2= " << seg->zSegment()->chi2() << endl;
  }

  // If both phi and zed projections are present and the phi segment is in time (segment t0<intime_cut) the 3-par fit is blocked and
  // segments are fit with the 2-par fit. Setting intime_cut to -1 results in the 3-par fit being used always.
  if (seg->hasPhi()) {
    if (seg->hasZed()) {
      if (fabs(seg->phiSegment()->t0()) < intime_cut) {
        fit(seg->phiSegment(), allow3par, true);
        fit(seg->zSegment(), allow3par, true);
      } else {
        fit(seg->phiSegment(), allow3par, false);
        fit(seg->zSegment(), allow3par, false);
      }
    } else
      fit(seg->phiSegment(), allow3par, false);
  } else
    fit(seg->zSegment(), allow3par, false);

  const DTChamber* theChamber = theGeom->chamber(seg->chamberId());

  if (seg->hasPhi() && seg->hasZed()) {
    DTChamberRecSegment2D* segPhi = seg->phiSegment();
    DTSLRecSegment2D* segZed = seg->zSegment();

    // NB Phi seg is already in chamber ref
    LocalPoint posPhiInCh = segPhi->localPosition();
    LocalVector dirPhiInCh = segPhi->localDirection();

    // Zed seg is in SL one
    const DTSuperLayer* zSL = theChamber->superLayer(segZed->superLayerId());
    LocalPoint zPos(segZed->localPosition().x(), (zSL->toLocal(theChamber->toGlobal(segPhi->localPosition()))).y(), 0.);

    LocalPoint posZInCh = theChamber->toLocal(zSL->toGlobal(zPos));

    LocalVector dirZInCh = theChamber->toLocal(zSL->toGlobal(segZed->localDirection()));

    LocalPoint posZAt0 = posZInCh + dirZInCh * (-posZInCh.z()) / cos(dirZInCh.theta());

    // given the actual definition of chamber refFrame, (with z poiniting to IP),
    // the zed component of direction is negative.
    LocalVector dir = LocalVector(dirPhiInCh.x() / fabs(dirPhiInCh.z()), dirZInCh.y() / fabs(dirZInCh.z()), -1.);

    seg->setPosition(LocalPoint(posPhiInCh.x(), posZAt0.y(), 0.));
    seg->setDirection(dir.unit());

    AlgebraicSymMatrix mat(4);

    // set cov matrix
    mat[0][0] = segPhi->parametersError()[0][0];  //sigma (dx/dz)
    mat[0][2] = segPhi->parametersError()[0][1];  //cov(dx/dz,x)
    mat[2][2] = segPhi->parametersError()[1][1];  //sigma (x)

    seg->setCovMatrix(mat);
    seg->setCovMatrixForZed(posZInCh);
  } else if (seg->hasPhi()) {
    DTChamberRecSegment2D* segPhi = seg->phiSegment();

    seg->setPosition(segPhi->localPosition());
    seg->setDirection(segPhi->localDirection());

    AlgebraicSymMatrix mat(4);
    // set cov matrix
    mat[0][0] = segPhi->parametersError()[0][0];  //sigma (dx/dz)
    mat[0][2] = segPhi->parametersError()[0][1];  //cov(dx/dz,x)
    mat[2][2] = segPhi->parametersError()[1][1];  //sigma (x)

    seg->setCovMatrix(mat);
  } else if (seg->hasZed()) {
    DTSLRecSegment2D* segZed = seg->zSegment();

    // Zed seg is in SL one
    GlobalPoint glbPosZ = (theGeom->superLayer(segZed->superLayerId()))->toGlobal(segZed->localPosition());
    LocalPoint posZInCh = (theGeom->chamber(segZed->superLayerId().chamberId()))->toLocal(glbPosZ);

    GlobalVector glbDirZ = (theGeom->superLayer(segZed->superLayerId()))->toGlobal(segZed->localDirection());
    LocalVector dirZInCh = (theGeom->chamber(segZed->superLayerId().chamberId()))->toLocal(glbDirZ);

    LocalPoint posZAt0 = posZInCh + dirZInCh * (-posZInCh.z()) / cos(dirZInCh.theta());

    seg->setPosition(posZAt0);
    seg->setDirection(dirZInCh);

    AlgebraicSymMatrix mat(4);
    // set cov matrix
    seg->setCovMatrix(mat);
    seg->setCovMatrixForZed(posZInCh);
  }
}

void DTSegmentUpdator::fit ( const std::vector< float > &  x, const std::vector< float > &  y, const std::vector< int > &  lfit, const std::vector< double > &  dist, const std::vector< float > &  sigy, LocalPoint &  pos, LocalVector &  dir, float &  cminf, float &  vminf, AlgebraicSymMatrix &  covMat, double &  chi2, const bool  allow3par = false, const bool  block3par = false  ) const

private

interface to LinearFit

Definition at line 355 of file DTSegmentUpdator.cc.

References HLT_FULL_cff::chi2, debug, mps_fire::i, intime_cut, slope, theFitter, and Vector3DBase< T, FrameTag >::unit().

Referenced by trackingPlots.Iteration::modules().

                                                       {
  float slope = 0.;
  float intercept = 0.;
  float covss = 0.;
  float covii = 0.;
  float covsi = 0.;

  cminf = 0;
  vminf = 0;

  int leftHits = 0, rightHits = 0;
  for (unsigned int i = 0; i < lfit.size(); i++)
    if (lfit[i] == 1)
      leftHits++;
    else
      rightHits++;

  theFitter->fit(x, y, x.size(), sigy, slope, intercept, chi2, covss, covii, covsi);

  // If we have at least one left and one right hit we can try the 3 parameter fit (if it is switched on)
  // FIXME: currently the covariance matrix from the 2-par fit is kept
  if (leftHits && rightHits && (leftHits + rightHits > 3) && allow3par) {
    theFitter->fitNpar(3, x, y, lfit, dist, sigy, slope, intercept, cminf, vminf, chi2, debug);
    double t0_corr = -cminf / 0.00543;
    if (fabs(t0_corr) < intime_cut && block3par) {
      theFitter->fit(x, y, x.size(), sigy, slope, intercept, chi2, covss, covii, covsi);
      cminf = 0;
    }
  }

  // cout << "slope " << slope << endl;
  // cout << "intercept " << intercept << endl;

  // intercept is the x() in chamber frame when the segment cross the chamber
  // plane (at z()=0), the y() is not measured, so let's put the center of the
  // chamber.
  pos = LocalPoint(intercept, 0., 0.);

  //  slope is dx()/dz(), while dy()/dz() is by definition 0, finally I want the
  //  segment to point outward, so opposite to local z
  dir = LocalVector(-slope, 0., -1.).unit();

  covMatrix = AlgebraicSymMatrix(2);
  covMatrix[0][0] = covss;  // this is var(dy/dz)
  covMatrix[1][1] = covii;  // this is var(y)
  covMatrix[1][0] = covsi;  // this is cov(dy/dz,y)
}

void DTSegmentUpdator::rejectBadHits ( DTChamberRecSegment2D *  phiSeg ) const

private

Definition at line 491 of file DTSegmentUpdator.cc.

References funct::abs(), gather_cfg::cout, debug, CommonMethods::delta(), TrackingRecHit::geographicalId(), mps_fire::i, N, DTRecSegment2D::specificRecHits(), theGeom, DTRecSegment2D::update(), x, PV3DBase< T, PVType, FrameType >::x(), y, and PV3DBase< T, PVType, FrameType >::z().

Referenced by update().

                                                                        {
  vector<float> x;
  vector<float> y;

  if (debug)
    cout << " Inside the segment updator, now loop on hits:   ( x == z_loc , y == x_loc) " << endl;

  vector<DTRecHit1D> hits = phiSeg->specificRecHits();
  const size_t N = hits.size();
  if (N < 3)
    return;

  for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
    // I have to get the hits position (the hit is in the layer rf) in SL frame...
    GlobalPoint glbPos = (theGeom->layer(hit->wireId().layerId()))->toGlobal(hit->localPosition());
    LocalPoint pos = (theGeom->idToDet(phiSeg->geographicalId()))->toLocal(glbPos);

    x.push_back(pos.z());
    y.push_back(pos.x());
  }

  if (debug) {
    cout << " end of segment! " << endl;
    cout << " size = Number of Hits: " << x.size() << "  " << y.size() << endl;
  }

  // Perform the 2 par fit:
  float par[2] = {0., 0.};  // q , m

  //variables to perform the fit:
  float Sx = 0.;
  float Sy = 0.;
  float Sx2 = 0.;
  float Sy2 = 0.;
  float Sxy = 0.;

  for (size_t i = 0; i < N; ++i) {
    Sx += x.at(i);
    Sy += y.at(i);
    Sx2 += x.at(i) * x.at(i);
    Sy2 += y.at(i) * y.at(i);
    Sxy += x.at(i) * y.at(i);
  }

  const float delta = N * Sx2 - Sx * Sx;
  par[0] = (Sx2 * Sy - Sx * Sxy) / delta;
  par[1] = (N * Sxy - Sx * Sy) / delta;

  if (debug)
    cout << "fit 2 parameters done ----> par0: " << par[0] << "  par1: " << par[1] << endl;

  // Calc residuals:
  float residuals[N];
  float mean_residual = 0.;  //mean of the absolute values of residuals
  for (size_t i = 0; i < N; ++i) {
    residuals[i] = y.at(i) - par[1] * x.at(i) - par[0];
    mean_residual += std::abs(residuals[i]);
    if (debug) {
      cout << " i: " << i << " y_i " << y.at(i) << " x_i " << x.at(i) << " res_i " << residuals[i];
      if (i == N - 1)
        cout << endl;
    }
  }

  if (debug)
    cout << " Residuals computed! " << endl;

  mean_residual = mean_residual / (N - 2);
  if (debug)
    cout << " mean_residual: " << mean_residual << endl;

  int i = 0;

  // Perform bad hit rejecting -- update hits
  vector<DTRecHit1D> updatedRecHits;
  for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
    float normResidual = mean_residual > 0 ? std::abs(residuals[i]) / mean_residual : 0;
    ++i;
    if (normResidual < 1.5) {
      DTRecHit1D newHit1D = (*hit);
      updatedRecHits.push_back(newHit1D);
      if (debug)
        cout << " accepted " << i << "th hit"
             << "  Irej: " << normResidual << endl;
    } else {
      if (debug)
        cout << " rejected " << i << "th hit"
             << "  Irej: " << normResidual << endl;
      continue;
    }
  }

  phiSeg->update(updatedRecHits);

  //final check!
  if (debug) {
    vector<float> x_upd;
    vector<float> y_upd;

    cout << " Check the update action: " << endl;

    vector<DTRecHit1D> hits_upd = phiSeg->specificRecHits();
    for (vector<DTRecHit1D>::const_iterator hit = hits_upd.begin(); hit != hits_upd.end(); ++hit) {
      // I have to get the hits position (the hit is in the layer rf) in SL frame...
      GlobalPoint glbPos = (theGeom->layer(hit->wireId().layerId()))->toGlobal(hit->localPosition());
      LocalPoint pos = (theGeom->idToDet(phiSeg->geographicalId()))->toLocal(glbPos);

      x_upd.push_back(pos.z());
      y_upd.push_back(pos.x());

      cout << " x_upd: " << pos.z() << "  y_upd: " << pos.x() << endl;
    }

    cout << " end of segment! " << endl;
    cout << " size = Number of Hits: " << x_upd.size() << "  " << y_upd.size() << endl;

  }  // end debug

  return;
}  //end DTSegmentUpdator::rejectBadHits

void DTSegmentUpdator::setES

(

const edm::EventSetup &

setup

)

set the setup

Definition at line 69 of file DTSegmentUpdator.cc.

References edm::EventSetup::getHandle(), theAlgo, theGeom, and theGeomToken.

Referenced by DTRefitAndCombineReco4D::setES(), DTMeantimerPatternReco4D::setES(), DTCombinatorialPatternReco4D::setES(), DTCombinatorialPatternReco::setES(), DTCombinatorialExtendedPatternReco::setES(), and DTMeantimerPatternReco::setES().

                                                    {
  theGeom = setup.getHandle(theGeomToken);
  theAlgo->setES(setup);
}

void DTSegmentUpdator::update	(	DTRecSegment4D *	seg,
		const bool	calcT0,
		bool	allow3par
	)		const

recompute hits position and refit the segment4D

Definition at line 74 of file DTSegmentUpdator.cc.

References calculateT0corr(), gather_cfg::cout, debug, DeadROC_duringRun::dir, fit(), TrackingRecHit::geographicalId(), DTRecSegment4D::hasPhi(), DTRecSegment4D::hasZed(), DTRecSegment4D::localDirection(), DTRecSegment4D::localPosition(), perform_delta_rejecting, DTRecSegment4D::phiSegment(), rejectBadHits(), theGeom, updateHits(), and DTRecSegment4D::zSegment().

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), MatrixUtil.Steps::overwrite(), DTRefitAndCombineReco4D::reconstruct(), DTMeantimerPatternReco4D::reconstruct(), DTCombinatorialPatternReco4D::reconstruct(), DTCombinatorialPatternReco::reconstruct(), DTCombinatorialExtendedPatternReco::reconstruct(), and DTMeantimerPatternReco::reconstruct().

                                                                                          {
  if (debug)
    cout << "[DTSegmentUpdator] Starting to update the DTRecSegment4D" << endl;

  const bool hasPhi = seg->hasPhi();
  const bool hasZed = seg->hasZed();

  //reject the bad hits (due to delta rays)
  if (perform_delta_rejecting && hasPhi)
    rejectBadHits(seg->phiSegment());

  int step = (hasPhi && hasZed) ? 3 : 2;
  if (calcT0)
    step = 4;

  if (debug)
    cout << "Step of update is " << step << endl;

  GlobalPoint pos = theGeom->idToDet(seg->geographicalId())->toGlobal(seg->localPosition());
  GlobalVector dir = theGeom->idToDet(seg->geographicalId())->toGlobal(seg->localDirection());

  if (calcT0)
    calculateT0corr(seg);

  if (hasPhi)
    updateHits(seg->phiSegment(), pos, dir, step);
  if (hasZed)
    updateHits(seg->zSegment(), pos, dir, step);

  fit(seg, allow3par);
}

void DTSegmentUpdator::update	(	DTRecSegment2D *	seg,
		bool	allow3par
	)		const

recompute hits position and refit the segment2D

Definition at line 106 of file DTSegmentUpdator.cc.

References gather_cfg::cout, debug, DeadROC_duringRun::dir, fit(), TrackingRecHit::geographicalId(), DTRecSegment2D::localDirection(), DTRecSegment2D::localPosition(), theGeom, and updateHits().

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                                       {
  if (debug)
    cout << "[DTSegmentUpdator] Starting to update the DTRecSegment2D" << endl;
  GlobalPoint pos = (theGeom->idToDet(seg->geographicalId()))->toGlobal(seg->localPosition());
  GlobalVector dir = (theGeom->idToDet(seg->geographicalId()))->toGlobal(seg->localDirection());

  updateHits(seg, pos, dir);
  fit(seg, allow3par, false);
}

void DTSegmentUpdator::updateHits ( DTRecSegment2D *  seg, GlobalPoint &  gpos, GlobalVector &  gdir, const int  step = 2  ) const

private

Definition at line 417 of file DTSegmentUpdator.cc.

References angle(), funct::cos(), HLT_FULL_cff::distance, relativeConstraints::error, Exception, DTRecSegment2D::ist0Valid(), phase1PixelTopology::layer, DTEnums::Left, convertSQLiteXML::ok, point, DTRecSegment2D::specificRecHits(), DTLayer::specificTopology(), DTRecSegment2D::t0(), T0_hit_resolution, theAlgo, theGeom, PV3DBase< T, PVType, FrameType >::theta(), GeomDet::toLocal(), DTRecSegment2D::update(), DTRecSegment2D::vDrift(), vdrift_4parfit, DTTopology::wirePosition(), PV3DBase< T, PVType, FrameType >::x(), hit::x, and PV3DBase< T, PVType, FrameType >::z().

Referenced by update().

                                                                                                                  {
  // it is not necessary to have DTRecHit1D* to modify the obj in the container
  // but I have to be carefully, since I cannot make a copy before the iteration!

  vector<DTRecHit1D> toBeUpdatedRecHits = seg->specificRecHits();
  vector<DTRecHit1D> updatedRecHits;

  for (vector<DTRecHit1D>::iterator hit = toBeUpdatedRecHits.begin(); hit != toBeUpdatedRecHits.end(); ++hit) {
    const DTLayer* layer = theGeom->layer(hit->wireId().layerId());

    LocalPoint segPos = layer->toLocal(gpos);
    LocalVector segDir = layer->toLocal(gdir);

    // define impact angle needed by the step 2
    const float angle = atan(segDir.x() / -segDir.z());

    // define the local position (extr.) of the segment. Needed by the third step
    LocalPoint segPosAtLayer = segPos + segDir * (-segPos.z()) / cos(segDir.theta());

    DTRecHit1D newHit1D = (*hit);
    bool ok = true;

    if (step == 2) {
      ok = theAlgo->compute(layer, *hit, angle, newHit1D);

    } else if (step == 3) {
      LocalPoint hitPos(hit->localPosition().x(), +segPosAtLayer.y(), 0.);
      GlobalPoint glbpos = theGeom->layer(hit->wireId().layerId())->toGlobal(hitPos);
      newHit1D.setPosition(hitPos);
      ok = theAlgo->compute(layer, *hit, angle, glbpos, newHit1D);

    } else if (step == 4) {
      //const double vminf = seg->vDrift();   //  vdrift correction are recorded in the segment
      double vminf = 0.;
      if (vdrift_4parfit)
        vminf = seg->vDrift();  // use vdrift recorded in the segment only if vdrift_4parfit=True

      double cminf = 0.;
      if (seg->ist0Valid())
        cminf = -seg->t0() * 0.00543;

      //cout << "In updateHits: t0 = " << seg->t0() << endl;
      //cout << "In updateHits: vminf = " << vminf << endl;
      //cout << "In updateHits: cminf = " << cminf << endl;

      const float xwire = layer->specificTopology().wirePosition(hit->wireId().wire());
      const float distance = fabs(hit->localPosition().x() - xwire);
      const int ilc = (hit->lrSide() == DTEnums::Left) ? 1 : -1;
      const double dy_corr = (vminf * ilc * distance - cminf * ilc);

      LocalPoint point(hit->localPosition().x() + dy_corr, +segPosAtLayer.y(), 0.);

      //double final_hit_resol = T0_hit_resolution;
      //if(newHit1D.wireId().layerId().superlayerId().superLayer() != 2) final_hit_resol = final_hit_resol * 0.8;
      //LocalError error(final_hit_resol * final_hit_resol,0.,0.);
      LocalError error(T0_hit_resolution * T0_hit_resolution, 0., 0.);
      newHit1D.setPositionAndError(point, error);

      //FIXME: check that the hit is still inside the cell
      ok = true;

    } else
      throw cms::Exception("DTSegmentUpdator") << " updateHits called with wrong step " << endl;

    if (ok)
      updatedRecHits.push_back(newHit1D);
    else {
      LogError("DTSegmentUpdator") << "DTSegmentUpdator::updateHits failed update" << endl;
      throw cms::Exception("DTSegmentUpdator") << "updateHits failed update" << endl;
    }
  }
  seg->update(updatedRecHits);
}

Member Data Documentation

bool DTSegmentUpdator::debug

private

Definition at line 114 of file DTSegmentUpdator.h.

Referenced by calculateT0corr(), DTSegmentUpdator(), fit(), runTauIdMVA.TauIDEmbedder::loadMVA_WPs_run2_2017(), rejectBadHits(), runTauIdMVA.TauIDEmbedder::runTauID(), and update().

double DTSegmentUpdator::intime_cut

private

Definition at line 110 of file DTSegmentUpdator.h.

Referenced by DTSegmentUpdator(), and fit().

bool DTSegmentUpdator::perform_delta_rejecting

private

Definition at line 113 of file DTSegmentUpdator.h.

Referenced by DTSegmentUpdator(), and update().

double DTSegmentUpdator::T0_hit_resolution

private

Definition at line 112 of file DTSegmentUpdator.h.

Referenced by DTSegmentUpdator(), and updateHits().

std::unique_ptr<DTRecHitBaseAlgo> DTSegmentUpdator::theAlgo

private

Definition at line 86 of file DTSegmentUpdator.h.

Referenced by DTSegmentUpdator(), setES(), and updateHits().

std::unique_ptr<DTLinearFit> DTSegmentUpdator::theFitter

private

Definition at line 85 of file DTSegmentUpdator.h.

Referenced by calculateT0corr(), and fit().

edm::ESHandle<DTGeometry> DTSegmentUpdator::theGeom

private

Definition at line 87 of file DTSegmentUpdator.h.

Referenced by calculateT0corr(), fit(), rejectBadHits(), setES(), update(), and updateHits().

const edm::ESGetToken<DTGeometry, MuonGeometryRecord> DTSegmentUpdator::theGeomToken

private

Definition at line 88 of file DTSegmentUpdator.h.

Referenced by DTSegmentUpdator(), and setES().

bool DTSegmentUpdator::vdrift_4parfit

private

Definition at line 111 of file DTSegmentUpdator.h.

Referenced by calculateT0corr(), DTSegmentUpdator(), and updateHits().