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) const
void	fit (DTRecSegment2D *seg) const
void	fit (DTRecSegment4D *seg) const
void	setES (const edm::EventSetup &setup)
	set the setup More...
void	update (DTRecSegment4D *seg, const bool calcT0=false) const
	recompute hits position and refit the segment4D More...
void	update (DTRecSegment2D *seg) 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< float > &sigy, LocalPoint &pos, LocalVector &dir, AlgebraicSymMatrix &covMat, double &chi2) const
	interface to LinearFit More...
void	Fit4Var (const std::vector< float > &xfit, const std::vector< float > &yfit, const std::vector< int > &lfit, const std::vector< double > &tfit, const int nptfit, float &cminf, double &vminf, double &chi2fit) const
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 549 of file DTSegmentUpdator.cc.

References Fit4Var(), TrackingRecHit::geographicalId(), DTEnums::Left, DTRecSegment2D::setT0(), DTRecSegment2D::setVdrift(), DTRecSegment2D::specificRecHits(), DTLayer::specificTopology(), theGeom, 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 

  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.;
  double vminf   = 0.;

  if ( nptfit > 2 ) {
    //NB chi2fit is normalized
    Fit4Var(x,y,lc,d_drift,nptfit,cminf,vminf,chi2fit);

    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 544 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

)

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

References dir, DTSegmentCand::good(), DTSegmentCand::hits(), DTSegmentCand::setChi2(), DTSegmentCand::setCovMatrix(), DTSegmentCand::setDirection(), DTSegmentCand::setPosition(), mathSSE::sqrt(), x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::y, and PV3DBase< T, PVType, FrameType >::z().

Referenced by DTMeantimerPatternReco::addHits(), DTCombinatorialPatternReco::buildPointsCollection(), DTCombinatorialExtendedPatternReco::buildPointsCollection(), DTCombinatorialExtendedPatternReco::extendCandidates(), fit(), DTRefitAndCombineReco4D::refitSuperSegments(), and update().

                                                   {
  if (!seg->good()) return false;

  vector<float> x;
  vector<float> y;
  vector<float> sigy;

  DTSegmentCand::AssPointCont hits=seg->hits();
  for (DTSegmentCand::AssPointCont::const_iterator iter=hits.begin();
       iter!=hits.end(); ++iter) {
    LocalPoint pos = (*iter).first->localPosition((*iter).second);
    x.push_back(pos.z()); 
    y.push_back(pos.x());
    sigy.push_back(sqrt((*iter).first->localPositionError().xx()));
  }

  LocalPoint pos;
  LocalVector dir;
  AlgebraicSymMatrix covMat(2);

  double chi2 = 0.;
  fit(x,y,sigy,pos,dir,covMat,chi2);

  seg->setPosition(pos);
  seg->setDirection(dir);

  //cout << "pos " << segPosition<< endl;
  //cout << "dir " << segDirection<< endl;

  seg->setCovMatrix(covMat);
  // cout << "Mat " << covMat << endl;

  seg->setChi2(chi2);
  return true;
}

void DTSegmentUpdator::fit

(

DTRecSegment2D *

seg

)

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

References dir, fit(), TrackingRecHit::geographicalId(), DTRecSegment2D::setChi2(), DTRecSegment2D::setCovMatrix(), DTRecSegment2D::setDirection(), DTRecSegment2D::setPosition(), DTRecSegment2D::specificRecHits(), mathSSE::sqrt(), theGeom, ErrorFrameTransformer::transform(), x, PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), detailsBasic3DVector::y, and PV3DBase< T, PVType, FrameType >::z().

                                                    {
  // WARNING: since this method is called both with a 2D and a 2DPhi as argument
  // seg->geographicalId() can be a superLayerId or a chamberId 

  vector<float> x;
  vector<float> y;
  vector<float> sigy;

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

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

  fit(x,y,sigy,pos,dir,covMat,chi2);

  seg->setPosition(pos);
  seg->setDirection(dir);

  //cout << "pos " << segPosition << endl;
  //cout << "dir " << segDirection << endl;

  seg->setCovMatrix(covMat);
  // cout << "Mat " << mat << endl;

  seg->setChi2(chi2);
}

void DTSegmentUpdator::fit

(

DTRecSegment4D *

seg

)

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

References DTSuperLayerId::chamberId(), DTRecSegment4D::chamberId(), funct::cos(), dir, fit(), DTRecSegment4D::hasPhi(), DTRecSegment4D::hasZed(), DTRecSegment2D::localDirection(), DTRecSegment2D::localPosition(), DTRecSegment2D::parametersError(), DTRecSegment4D::phiSegment(), DTRecSegment4D::setCovMatrix(), DTRecSegment4D::setCovMatrixForZed(), DTRecSegment4D::setDirection(), DTRecSegment4D::setPosition(), DTChamber::superLayer(), DTSLRecSegment2D::superLayerId(), 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().

                                                     {
  // after the update must refit the segments
  if(seg->hasPhi()) fit(seg->phiSegment());
  if(seg->hasZed()) fit(seg->zSegment());

  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< float > &  sigy, LocalPoint &  pos, LocalVector &  dir, AlgebraicSymMatrix &  covMat, double &  chi2  ) const

private

interface to LinearFit

Definition at line 278 of file DTSegmentUpdator.cc.

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

                                                {

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

  // do the fit
  theFitter->fit(x,y,x.size(),sigy,slope,intercept,covss,covii,covsi);
  // 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)

  /* Calculate chi2. */
  chi2 = 0.;
  for(unsigned int i=0; i<x.size() ; ++i) {
    double resid= y[i] - (intercept + slope*x[i]);
    chi2 += (resid/sigy[i])*(resid/sigy[i]);
  }
}

void DTSegmentUpdator::Fit4Var ( const std::vector< float > &  xfit, const std::vector< float > &  yfit, const std::vector< int > &  lfit, const std::vector< double > &  tfit, const int  nptfit, float &  cminf, double &  vminf, double &  chi2fit  ) const

private

Definition at line 607 of file DTSegmentUpdator.cc.

References a, b, alignmentValidation::c1, gather_cfg::cout, debug, delta, j, and vdrift_4parfit.

Referenced by calculateT0corr().

                                                      { 

  const double sigma = 0.0295;// errors can be inserted .just load them/that is the usual TB resolution value for DT chambers 
  double aminf = 0.;
  double bminf = 0.;
  int nppar = 0;
  double sx = 0.;
  double  sx2 = 0.;
  double sy = 0.;
  double sxy = 0.;
  double sl = 0.;
  double sl2 = 0.;
  double sly = 0.;
  double slx = 0.;
  double st = 0.;
  double st2 = 0.;
  double slt = 0.;
  double sltx = 0.;
  double slty = 0.;
  double chi2fitN2 = -1. ;
  double chi2fit3 = -1.;
  double chi2fitN3 = -1. ;
  double chi2fitN4 = -1.;
  float bminf3 = bminf;
  float aminf3 = aminf;
  float cminf3 = cminf;
  int nppar2 = 0;
  int nppar3 = 0;
  int nppar4 = 0;

  cminf = -999.;
  vminf = 0.;

  for (int j=0; j<nptfit; j++){
    sx  = sx + xfit[j];       
    sy  = sy + yfit[j];
    sx2 = sx2 + xfit[j]*xfit[j];
    sxy = sxy + xfit[j]*yfit[j];
    sl  = sl + lfit[j];       
    sl2 = sl2 + lfit[j]*lfit[j];
    sly = sly + lfit[j]*yfit[j];
    slx = slx + lfit[j]*xfit[j];
    st = st + tfit[j];
    st2 = st2 + tfit[j] * tfit[j];
    slt = slt + lfit[j] * tfit[j];
    sltx = sltx + lfit[j] * tfit[j]*xfit[j];
    slty = slty + lfit[j] * tfit[j]*yfit[j];

  } //end loop

  const double delta = nptfit*sx2 - sx*sx;

  double a = 0.;
  double b = 0.;           

  if (delta!=0){   //
    a = (sx2*sy - sx*sxy)/delta;
    b = (nptfit*sxy - sx*sy)/delta;

    //  cout << " NPAR=2 : slope = "<<b<< "    intercept = "<<a <<endl;
    for (int j=0; j<nptfit; j++){
      const double ypred = a + b*xfit[j];
      const double dy = (yfit[j] - ypred)/sigma;
      chi2fit = chi2fit + dy*dy;
    } //end loop chi2
  }

  bminf = b;
  aminf = a;

  nppar = 2; 
  nppar2 = nppar; 

  chi2fitN2 = chi2fit/(nptfit-2);

  // cout << "dt0 = 0chi2fit = " << chi2fit << "  slope = "<<b<<endl;

  if (nptfit >= 3) {

    const double d1 = sy;
    const double d2 = sxy;
    const double d3 = sly;
    const double c1 = sl;
    const double c2 = slx;
    const double c3 = sl2;
    const double b1 = sx;
    const double b2 = sx2;
    const double b3 = slx;
    const double a1 = nptfit;
    const double a2 = sx;
    const double a3 = sl;

    //these parameters are not used in the 4-variables fit
    const double b4 = b2*a1-b1*a2;
    const double c4 = c2*a1-c1*a2;
    const double d4 = d2*a1-d1*a2;
    const double b5 = a1*b3-a3*b1;
    const double c5 = a1*c3-a3*c1;
    const double d5 = a1*d3-d1*a3;
    const double a6 = slt;
    const double b6 = sltx;
    const double c6 = st;
    const double v6 = st2;  
    const double d6 = slty;

    if (((c5*b4-c4*b5)*b4*a1)!=0) {
      nppar = 3;
      chi2fit = 0.;
      cminf = (d5*b4-d4*b5)/(c5*b4-c4*b5);
      bminf = d4/b4 -cminf *c4/b4;
      aminf = (d1/a1 -cminf*c1/a1 -bminf*b1/a1);

      for (int j=0; j<nptfit; j++){
        const double ypred = aminf + bminf*xfit[j];
        const double dy = (yfit[j]-cminf*lfit[j] - ypred)/sigma;
        chi2fit = chi2fit + dy*dy;

      } //end loop chi2
      chi2fit3 = chi2fit;
      if (nptfit>3)
        chi2fitN3 = chi2fit /(nptfit-3);

    }
    else {
      cminf = -999.;
      bminf = b;
      aminf = a;
      chi2fit3 = chi2fit;
      chi2fitN3 = chi2fit /(nptfit-2);
    }

    bminf3 = bminf;
    aminf3 = aminf;
    cminf3 = cminf;
    nppar3 = nppar;

    if (debug) {
      cout << "dt0= 0 : slope 2 = " << b << " pos in  = " << a << " chi2fitN2 = " << chi2fitN2
       << " nppar = " << nppar2 << " nptfit = " << nptfit << endl;
      cout << "dt0 = 0 : slope 3 = " << bminf << " pos out = " << aminf << " chi2fitN3 = "
       << chi2fitN3 << " nppar = " << nppar3 << " T0_ev ns = " << cminf/0.00543 << endl;
    } 

    //***********************************
    //     cout << " vdrift_4parfit "<< vdrift_4parfit<<endl;
    if( nptfit>=5) { 
      const double det = (a1*a1*(b2*v6 - b6*b6) - a1*(a2*a2*v6 - 2*a2*a6*b6 + a6*a6*b2 + b2*c6*c6 + b3*(b3*v6 - 2*b6*c6))
              + a2*a2*c6*c6 + 2*a2*(a3*(b3*v6 - b6*c6) - a6*b3*c6) + a3*a3*(b6*b6 - b2*v6)
              + a6*(2*a3*(b2*c6 - b3*b6) + a6*b3*b3)); 

      // the dv/vdrift correction may be computed  under vdrift_4parfit request;
      if (det != 0) { 
        nppar = 4;
        chi2fit = 0.;
        // computation of   a, b, c e v
        aminf = (a1*(a2*(b6*d6 - v6*d2) + a6*(b6*d2 - b2*d6) + d1*(b2*v6 - b6*b6)) - a2*(b3*(c6*d6 - v6*d3)
                 + c6*(b6*d3 - c6*d2)) + a3*(b2*(c6*d6 - v6*d3) + b3*(v6*d2 - b6*d6) + b6*(b6*d3 - c6*d2))
                 + a6*(b2*c6*d3 + b3*(b3*d6 - b6*d3 - c6*d2)) - d1*(b2*c6*c6 + b3*(b3*v6 - 2*b6*c6)))/det;
        bminf = - (a1*a1*(b6*d6 - v6*d2) - a1*(a2*(a6*d6 - v6*d1) - a6*a6*d2 + a6*b6*d1 + b3*(c6*d6 - v6*d3)
                 + c6*(b6*d3 - c6*d2)) + a2*(a3*(c6*d6 - v6*d3) + c6*(a6*d3 - c6*d1)) + a3*a3*(v6*d2 - b6*d6)
                 + a3*(a6*(b3*d6 + b6*d3 - 2*c6*d2) - d1*(b3*v6 - b6*c6)) - a6*b3*(a6*d3 - c6*d1))/det;
        cminf = -(a1*(b2*(c6*d6 - v6*d3) + b3*(v6*d2 - b6*d6) + b6*(b6*d3 - c6*d2)) + a2*a2*(v6*d3 - c6*d6)
                 + a2*(a3*(b6*d6 - v6*d2) + a6*(b3*d6 - 2*b6*d3 + c6*d2) - d1*(b3*v6 - b6*c6))
                 + a3*(d1*(b2*v6 - b6*b6) - a6*(b2*d6 - b6*d2)) + a6*(a6*(b2*d3 - b3*d2) - d1*(b2*c6 - b3*b6)))/det;
        vminf = - (a1*a1*(b2*d6 - b6*d2) - a1*(a2*a2*d6 - a2*(a6*d2 + b6*d1) + a6*b2*d1 + b2*c6*d3
                 + b3*(b3*d6 - b6*d3 - c6*d2)) + a2*a2*c6*d3 + a2*(a3*(2*b3*d6 - b6*d3 - c6*d2) - b3*(a6*d3 + c6*d1))
                 + a3*a3*(b6*d2 - b2*d6) + a3*(a6*(b2*d3 - b3*d2) + d1*(b2*c6 - b3*b6)) + a6*b3*b3*d1)/det;

        //  chi 2
        for (int j=0; j<nptfit; j++) {
          const double ypred = aminf + bminf*xfit[j];
          const double dy = (yfit[j]+vminf*lfit[j]*tfit[j]-cminf*lfit[j] -ypred)/sigma; 
          chi2fit = chi2fit + dy*dy;

        } //end loop chi2
        if (nptfit<=nppar){ 
          chi2fitN4=-1;
          //        cout << "nptfit " << nptfit << " nppar " << nppar << endl;
        }
        else{
          chi2fitN4= chi2fit / (nptfit-nppar); 
        }
      }
      else {
        vminf = 0.;

        if (nptfit <= nppar) chi2fitN4=-1;
        else chi2fitN4  = chi2fit / (nptfit-nppar); 
      }

      if (fabs(vminf) >= 0.29) {
        // for safety and for code construction..dont accept correction on dv/vdrift greater then 0.09
        vminf = 0.;
        cminf = cminf3;
        aminf = aminf3;
        bminf = bminf3;
        nppar = 3;
        chi2fit = chi2fit3;
      }

    }  //end if vdrift

     if(!vdrift_4parfit){         //if not required explicitly leave the t0 and track step as at step 3
                                  // just update vdrift value vmin for storing in the segments for monitoring
       cminf = cminf3;
       aminf = aminf3;
       bminf = bminf3;
       nppar = 3;
       chi2fit = chi2fit3;
     }

    nppar4 = nppar;

  }  //end nptfit >=3

  if (debug) {
    cout << "   dt0= 0 : slope 4  = " << bminf << " pos out = " << aminf <<" chi2fitN4 = " << chi2fitN4
     << "  nppar= " << nppar4 << " T0_ev ns= " << cminf/0.00543 <<" delta v = " << vminf <<endl;
    cout << nptfit << " nptfit " << " end  chi2fit = " << chi2fit/ (nptfit-nppar ) << " T0_ev ns= " << cminf/0.00543 << " delta v = " << vminf <<endl;
  }

  if ( fabs(vminf) >= 0.09 && debug ) {  //checks only vdrift less then 10 % accepted
    cout << "vminf gt 0.09 det=  " << endl;
    cout << "dt0= 0 : slope 4 = "<< bminf << " pos out = " << aminf << " chi2fitN4 = " << chi2fitN4
     << " T0_ev ns = " << cminf/0.00543 << " delta v = "<< vminf << endl;
    cout << "dt0 = 0 : slope 2 = "<< b << " pos in = " << a <<" chi2fitN2 = " << chi2fitN2
     << " nppar = " << nppar-1 << " nptfit = " << nptfit <<endl;
    cout << "dt0 = 0 : slope 3 = " << bminf << " pos out = " << aminf << " chi2fitN3 = "
     << chi2fitN3 << " T0_ev ns = " << cminf/0.00543 << endl;
    cout << nptfit   <<" nptfit "<< "   end  chi2fit = " << chi2fit << "T0_ev ns= " << cminf/0.00543 << "delta v = "<< vminf <<endl;        
  }

  if (nptfit != nppar) chi2fit = chi2fit / (nptfit-nppar);
}

void DTSegmentUpdator::rejectBadHits ( DTChamberRecSegment2D *  phiSeg ) const

private

Definition at line 406 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 = false
	)		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 segment" << 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);
}

void DTSegmentUpdator::update

(

DTRecSegment2D *

seg

)

const

recompute hits position and refit the segment2D

Definition at line 100 of file DTSegmentUpdator.cc.

References 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().

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

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

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

private

Definition at line 321 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 114 of file DTSegmentUpdator.h.

Referenced by DTSegmentUpdator(), Fit4Var(), rejectBadHits(), and update().

bool DTSegmentUpdator::perform_delta_rejecting

private

Definition at line 113 of file DTSegmentUpdator.h.

Referenced by update().

double DTSegmentUpdator::T0_hit_resolution

private

Definition at line 112 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 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 111 of file DTSegmentUpdator.h.

Referenced by Fit4Var(), and updateHits().