DTSegmentUpdator Class Reference

#include <DTSegmentUpdator.h>

Public Member Functions
void	calculateT0corr (DTRecSegment2D *seg) const
void	calculateT0corr (DTRecSegment4D *seg) const
	DTSegmentUpdator (const edm::ParameterSet &config)
	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
bool	perform_delta_rejecting
double	T0_hit_resolution
DTRecHitBaseAlgo *	theAlgo
DTLinearFit *	theFitter
edm::ESHandle< DTGeometry >	theGeom
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 40 of file DTSegmentUpdator.h.

Constructor & Destructor Documentation

DTSegmentUpdator::DTSegmentUpdator

(

const edm::ParameterSet &

config

)

Constructor.

Definition at line 47 of file DTSegmentUpdator.cc.

References gather_cfg::cout, debug, reco::get(), edm::ParameterSet::getParameter(), and theAlgo.

                                                             :
  theFitter(new DTLinearFit()) ,
  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)) 
{  
  string theAlgoName = config.getParameter<string>("recAlgo");
  theAlgo = DTRecHitAlgoFactory::get()->create(theAlgoName, 
                                               config.getParameter<ParameterSet>("recAlgoConfig"));

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

DTSegmentUpdator::~DTSegmentUpdator

(

)

Destructor.

Definition at line 63 of file DTSegmentUpdator.cc.

References theFitter.

                                    {
  delete theFitter;
}

Member Function Documentation

void DTSegmentUpdator::calculateT0corr

(

DTRecSegment2D *

seg

)

const

Definition at line 626 of file DTSegmentUpdator.cc.

References a, b, gather_cfg::cout, debug, DTLinearFit::fit4Var(), TrackingRecHit::geographicalId(), DTEnums::Left, DTRecSegment2D::setT0(), DTRecSegment2D::setVdrift(), DTRecSegment2D::specificRecHits(), DTLayer::specificTopology(), theFitter, theGeom, vdrift_4parfit, DTTopology::wirePosition(), x, PV3DBase< T, PVType, FrameType >::x(), hit::x, detailsBasic3DVector::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 621 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 215 of file DTSegmentUpdator.cc.

References gather_cfg::cout, debug, dir, DTSegmentCand::good(), DTSegmentCand::hits(), i, getDQMSummary::iter, DTEnums::Left, DTEnums::Right, DTSegmentCand::setChi2(), DTSegmentCand::setCovMatrix(), DTSegmentCand::setDirection(), DTSegmentCand::setPosition(), DTSegmentCand::sett0(), mathSSE::sqrt(), x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::y, and PV3DBase< T, PVType, FrameType >::z().

Referenced by DTCombinatorialPatternReco::buildPointsCollection(), DTCombinatorialExtendedPatternReco::buildPointsCollection(), DTCombinatorialExtendedPatternReco::extendCandidates(), fit(), DTMeantimerPatternReco::fitWithT0(), 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);

  DTSegmentCand::AssPointCont hits=seg->hits();
  for (DTSegmentCand::AssPointCont::const_iterator iter=hits.begin(); iter!=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 << endl;
  if (debug && fitdebug) cout << "  DTcand   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 280 of file DTSegmentUpdator.cc.

References gather_cfg::cout, debug, dir, fit(), TrackingRecHit::geographicalId(), 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(), detailsBasic3DVector::y, and PV3DBase< T, PVType, FrameType >::z().

                                                                                    {

  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 109 of file DTSegmentUpdator.cc.

References DTSuperLayerId::chamberId(), DTRecSegment4D::chamberId(), DTRecSegment2D::chi2(), funct::cos(), gather_cfg::cout, debug, dir, fit(), DTRecSegment4D::hasPhi(), DTRecSegment4D::hasZed(), 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(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), zPos, and DTRecSegment4D::zSegment().

                                                                     {
  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;
  }
  
  
  // use the 2-par fit UNLESS the 4D segment only has the Phi projection or is out of time by more than 20ns
  // - in these cases use the 3-par fit
  if(seg->hasPhi()) {
    if(seg->hasZed()) {

      // fit in-time Phi segments with the 2par fit and out-of-time segments with the 3par fit
      if (fabs(seg->phiSegment()->t0())<40.) {
        fit(seg->phiSegment(),allow3par,1);
        fit(seg->zSegment(),allow3par,1);      
      } else {
        fit(seg->phiSegment(),allow3par,0);
        fit(seg->zSegment(),allow3par,0);      
      }

    } else fit(seg->phiSegment(),allow3par,0);
  } else fit(seg->zSegment(),allow3par,0);
 
  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 349 of file DTSegmentUpdator.cc.

References debug, DTLinearFit::fit(), DTLinearFit::fitNpar(), i, slope, theFitter, and Vector3DBase< T, FrameTag >::unit().

                                                        {

  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)<20. && 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 486 of file DTSegmentUpdator.cc.

References gather_cfg::cout, debug, delta, TrackingRecHit::geographicalId(), i, N, DTRecSegment2D::specificRecHits(), theGeom, DTRecSegment2D::update(), x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::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();
  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.;

  size_t N =  x.size();

  if (N == 0)
    return;
    
  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];
    
  for(size_t i = 0; i < N;++i)
    residuals[i] = 0;
    
  for(size_t i = 0; i < N;++i){
    residuals[i] = y.at(i) - par[1]*x.at(i) - par[0];
    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;
            
  // Perform bad hit rejecting -- update hits
  vector<DTRecHit1D> updatedRecHits;
    
  float mean_residual = 0.; //mean of the absolute values of residuals
    
  for (size_t i = 0; i < N; ++i)
    mean_residual += fabs(residuals[i]);
    
  mean_residual = mean_residual/(N - 2);    
    
  if(debug) cout << " mean_residual: "<< mean_residual << endl;
    
  int i = 0;
    
  for (vector<DTRecHit1D>::const_iterator hit=hits.begin();
       hit!=hits.end(); ++hit) {
        
    DTRecHit1D newHit1D = (*hit);

    float normResidual = mean_residual > 0 ? fabs(residuals[i])/mean_residual : 0;
    if(normResidual < 1.5){
                    
      updatedRecHits.push_back(newHit1D);
      if(debug) cout << " accepted "<< i+1 << "th hit" <<"  Irej: " << normResidual << endl;
      ++i;
    }
    else {
      if(debug) cout << " rejected "<< i+1 << "th hit" <<"  Irej: " << normResidual << endl;
      ++i;
      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::get(), DTRecHitBaseAlgo::setES(), theAlgo, and theGeom.

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

                                                   {
  setup.get<MuonGeometryRecord>().get(theGeom);
  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, dir, fit(), TrackingRecHit::geographicalId(), DTRecSegment4D::hasPhi(), DTRecSegment4D::hasZed(), DTRecSegment4D::localDirection(), DTRecSegment4D::localPosition(), perform_delta_rejecting, DTRecSegment4D::phiSegment(), rejectBadHits(), relval_parameters_module::step, theGeom, updateHits(), and DTRecSegment4D::zSegment().

Referenced by progressbar.ProgressBar::__next__(), relval_steps.Matrix::__setitem__(), relval_steps.Steps::__setitem__(), Vispa.Gui.VispaWidget.VispaWidget::autosize(), Vispa.Views.LineDecayView.LineDecayContainer::createObject(), Vispa.Views.LineDecayView.LineDecayContainer::deselectAllObjects(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::deselectAllWidgets(), Vispa.Gui.VispaWidget.VispaWidget::enableAutosizing(), progressbar.ProgressBar::finish(), Vispa.Gui.MenuWidget.MenuWidget::leaveEvent(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseMoveEvent(), Vispa.Gui.MenuWidget.MenuWidget::mouseMoveEvent(), Vispa.Views.LineDecayView.LineDecayContainer::mouseMoveEvent(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseReleaseEvent(), Vispa.Views.LineDecayView.LineDecayContainer::objectMoved(), relval_steps.Steps::overwrite(), DTMeantimerPatternReco4D::reconstruct(), DTCombinatorialPatternReco4D::reconstruct(), DTRefitAndCombineReco4D::reconstruct(), DTCombinatorialPatternReco::reconstruct(), DTCombinatorialExtendedPatternReco::reconstruct(), DTMeantimerPatternReco::reconstruct(), Vispa.Views.LineDecayView.LineDecayContainer::removeObject(), Vispa.Gui.ConnectableWidget.ConnectableWidget::removePorts(), Vispa.Gui.FindDialog.FindDialog::reset(), Vispa.Gui.PortConnection.PointToPointConnection::select(), Vispa.Gui.VispaWidget.VispaWidget::select(), Vispa.Views.LineDecayView.LineDecayContainer::select(), Vispa.Gui.VispaWidget.VispaWidget::setText(), Vispa.Gui.VispaWidget.VispaWidget::setTitle(), Vispa.Gui.ZoomableWidget.ZoomableWidget::setZoom(), Vispa.Views.LineDecayView.LineDecayContainer::setZoom(), and Vispa.Gui.PortConnection.PointToPointConnection::updateConnection().

                                                                                          {

  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 100 of file DTSegmentUpdator.cc.

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

Referenced by progressbar.ProgressBar::__next__(), relval_steps.Matrix::__setitem__(), relval_steps.Steps::__setitem__(), Vispa.Gui.VispaWidget.VispaWidget::autosize(), Vispa.Views.LineDecayView.LineDecayContainer::createObject(), Vispa.Views.LineDecayView.LineDecayContainer::deselectAllObjects(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::deselectAllWidgets(), Vispa.Gui.VispaWidget.VispaWidget::enableAutosizing(), progressbar.ProgressBar::finish(), Vispa.Gui.MenuWidget.MenuWidget::leaveEvent(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseMoveEvent(), Vispa.Gui.MenuWidget.MenuWidget::mouseMoveEvent(), Vispa.Views.LineDecayView.LineDecayContainer::mouseMoveEvent(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseReleaseEvent(), Vispa.Views.LineDecayView.LineDecayContainer::objectMoved(), relval_steps.Steps::overwrite(), Vispa.Views.LineDecayView.LineDecayContainer::removeObject(), Vispa.Gui.ConnectableWidget.ConnectableWidget::removePorts(), Vispa.Gui.FindDialog.FindDialog::reset(), Vispa.Gui.PortConnection.PointToPointConnection::select(), Vispa.Gui.VispaWidget.VispaWidget::select(), Vispa.Views.LineDecayView.LineDecayContainer::select(), Vispa.Gui.VispaWidget.VispaWidget::setText(), Vispa.Gui.VispaWidget.VispaWidget::setTitle(), Vispa.Gui.ZoomableWidget.ZoomableWidget::setZoom(), Vispa.Views.LineDecayView.LineDecayContainer::setZoom(), and Vispa.Gui.PortConnection.PointToPointConnection::updateConnection().

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

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

private

Definition at line 409 of file DTSegmentUpdator.cc.

References angle(), DTRecHitBaseAlgo::compute(), funct::cos(), relativeConstraints::error, edm::hlt::Exception, DTRecSegment2D::ist0Valid(), 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 112 of file DTSegmentUpdator.h.

Referenced by calculateT0corr(), DTSegmentUpdator(), fit(), rejectBadHits(), and update().

bool DTSegmentUpdator::perform_delta_rejecting

private

Definition at line 111 of file DTSegmentUpdator.h.

Referenced by update().

double DTSegmentUpdator::T0_hit_resolution

private

Definition at line 110 of file DTSegmentUpdator.h.

Referenced by updateHits().

DTRecHitBaseAlgo* DTSegmentUpdator::theAlgo

private

Definition at line 83 of file DTSegmentUpdator.h.

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

DTLinearFit* DTSegmentUpdator::theFitter

private

Definition at line 82 of file DTSegmentUpdator.h.

Referenced by calculateT0corr(), fit(), and ~DTSegmentUpdator().

edm::ESHandle<DTGeometry> DTSegmentUpdator::theGeom

private

Definition at line 84 of file DTSegmentUpdator.h.

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

bool DTSegmentUpdator::vdrift_4parfit

private

Definition at line 109 of file DTSegmentUpdator.h.

Referenced by calculateT0corr(), and updateHits().