Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-09-16 03:24:10 

0001 /** \file
0002  *
0003  * \author Stefano Lacaprara - INFN Legnaro <stefano.lacaprara@pd.infn.it>
0004  * \author Riccardo Bellan - INFN TO <riccardo.bellan@cern.ch>
0005  * \       A.Meneguzzo - Padova University  <anna.meneguzzo@pd.infn.it>
0006  * \       M.Pelliccioni - INFN TO <pellicci@cern.ch>
0007  * \       M.Meneghelli - INFN BO <marco.meneghelli@cern.ch>
0008  */
0009 
0010 /* This Class Header */
0011 #include "RecoLocalMuon/DTSegment/src/DTSegmentUpdator.h"
0012 
0013 /* Collaborating Class Header */
0014 
0015 //mene
0016 #include "DataFormats/DTRecHit/interface/DTChamberRecSegment2D.h"
0017 
0018 #include "DataFormats/DTRecHit/interface/DTRecSegment2D.h"
0019 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0020 #include "DataFormats/DTRecHit/interface/DTChamberRecSegment2D.h"
0021 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
0022 #include "DataFormats/DTRecHit/interface/DTRecHit1D.h"
0023 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
0024 
0025 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h"
0026 #include "RecoLocalMuon/DTRecHit/interface/DTRecHitAlgoFactory.h"
0027 #include "RecoLocalMuon/DTSegment/src/DTLinearFit.h"
0028 
0029 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0030 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0031 #include "Geometry/DTGeometry/interface/DTLayer.h"
0032 
0033 #include "FWCore/Utilities/interface/Exception.h"
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #include "FWCore/Framework/interface/EventSetup.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 
0038 /* C++ Headers */
0039 #include <string>
0040 
0041 using namespace std;
0042 using namespace edm;
0043 
0044 /* ====================================================================== */
0045 
0046 /// Constructor
0047 DTSegmentUpdator::DTSegmentUpdator(const ParameterSet& config, edm::ConsumesCollector cc)
0048     : theFitter{std::make_unique<DTLinearFit>()},
0049       theAlgo{DTRecHitAlgoFactory::get()->create(
0050           config.getParameter<string>("recAlgo"), config.getParameter<ParameterSet>("recAlgoConfig"), cc)},
0051       theGeomToken(cc.esConsumes()),
0052       vdrift_4parfit(config.getParameter<bool>("performT0_vdriftSegCorrection")),
0053       T0_hit_resolution(config.getParameter<double>("hit_afterT0_resolution")),
0054       perform_delta_rejecting(config.getParameter<bool>("perform_delta_rejecting")),
0055       debug(config.getUntrackedParameter<bool>("debug", false)) {
0056   intime_cut = 20.;
0057   if (config.exists("intime_cut"))
0058     intime_cut = config.getParameter<double>("intime_cut");
0059 
0060   if (debug)
0061     cout << "[DTSegmentUpdator] Constructor called" << endl;
0062 }
0063 
0064 /// Destructor
0065 DTSegmentUpdator::~DTSegmentUpdator() = default;
0066 
0067 /* Operations */
0068 
0069 void DTSegmentUpdator::setES(const EventSetup& setup) {
0070   theGeom = setup.getHandle(theGeomToken);
0071   theAlgo->setES(setup);
0072 }
0073 
0074 void DTSegmentUpdator::update(DTRecSegment4D* seg, const bool calcT0, bool allow3par) const {
0075   if (debug)
0076     cout << "[DTSegmentUpdator] Starting to update the DTRecSegment4D" << endl;
0077 
0078   const bool hasPhi = seg->hasPhi();
0079   const bool hasZed = seg->hasZed();
0080 
0081   //reject the bad hits (due to delta rays)
0082   if (perform_delta_rejecting && hasPhi)
0083     rejectBadHits(seg->phiSegment());
0084 
0085   int step = (hasPhi && hasZed) ? 3 : 2;
0086   if (calcT0)
0087     step = 4;
0088 
0089   if (debug)
0090     cout << "Step of update is " << step << endl;
0091 
0092   GlobalPoint pos = theGeom->idToDet(seg->geographicalId())->toGlobal(seg->localPosition());
0093   GlobalVector dir = theGeom->idToDet(seg->geographicalId())->toGlobal(seg->localDirection());
0094 
0095   if (calcT0)
0096     calculateT0corr(seg);
0097 
0098   if (hasPhi)
0099     updateHits(seg->phiSegment(), pos, dir, step);
0100   if (hasZed)
0101     updateHits(seg->zSegment(), pos, dir, step);
0102 
0103   fit(seg, allow3par);
0104 }
0105 
0106 void DTSegmentUpdator::update(DTRecSegment2D* seg, bool allow3par) const {
0107   if (debug)
0108     cout << "[DTSegmentUpdator] Starting to update the DTRecSegment2D" << endl;
0109   GlobalPoint pos = (theGeom->idToDet(seg->geographicalId()))->toGlobal(seg->localPosition());
0110   GlobalVector dir = (theGeom->idToDet(seg->geographicalId()))->toGlobal(seg->localDirection());
0111 
0112   updateHits(seg, pos, dir);
0113   fit(seg, allow3par, false);
0114 }
0115 
0116 void DTSegmentUpdator::fit(DTRecSegment4D* seg, bool allow3par) const {
0117   if (debug)
0118     cout << "[DTSegmentUpdator] Fit DTRecSegment4D:" << endl;
0119   // after the update must refit the segments
0120 
0121   if (debug) {
0122     if (seg->hasPhi())
0123       cout << "    4D Segment contains a Phi segment. t0= " << seg->phiSegment()->t0()
0124            << "  chi2= " << seg->phiSegment()->chi2() << endl;
0125     if (seg->hasZed())
0126       cout << "    4D Segment contains a Zed segment. t0= " << seg->zSegment()->t0()
0127            << "  chi2= " << seg->zSegment()->chi2() << endl;
0128   }
0129 
0130   // 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
0131   // segments are fit with the 2-par fit. Setting intime_cut to -1 results in the 3-par fit being used always.
0132   if (seg->hasPhi()) {
0133     if (seg->hasZed()) {
0134       if (fabs(seg->phiSegment()->t0()) < intime_cut) {
0135         fit(seg->phiSegment(), allow3par, true);
0136         fit(seg->zSegment(), allow3par, true);
0137       } else {
0138         fit(seg->phiSegment(), allow3par, false);
0139         fit(seg->zSegment(), allow3par, false);
0140       }
0141     } else
0142       fit(seg->phiSegment(), allow3par, false);
0143   } else
0144     fit(seg->zSegment(), allow3par, false);
0145 
0146   const DTChamber* theChamber = theGeom->chamber(seg->chamberId());
0147 
0148   if (seg->hasPhi() && seg->hasZed()) {
0149     DTChamberRecSegment2D* segPhi = seg->phiSegment();
0150     DTSLRecSegment2D* segZed = seg->zSegment();
0151 
0152     // NB Phi seg is already in chamber ref
0153     LocalPoint posPhiInCh = segPhi->localPosition();
0154     LocalVector dirPhiInCh = segPhi->localDirection();
0155 
0156     // Zed seg is in SL one
0157     const DTSuperLayer* zSL = theChamber->superLayer(segZed->superLayerId());
0158     LocalPoint zPos(segZed->localPosition().x(), (zSL->toLocal(theChamber->toGlobal(segPhi->localPosition()))).y(), 0.);
0159 
0160     LocalPoint posZInCh = theChamber->toLocal(zSL->toGlobal(zPos));
0161 
0162     LocalVector dirZInCh = theChamber->toLocal(zSL->toGlobal(segZed->localDirection()));
0163 
0164     LocalPoint posZAt0 = posZInCh + dirZInCh * (-posZInCh.z()) / cos(dirZInCh.theta());
0165 
0166     // given the actual definition of chamber refFrame, (with z poiniting to IP),
0167     // the zed component of direction is negative.
0168     LocalVector dir = LocalVector(dirPhiInCh.x() / fabs(dirPhiInCh.z()), dirZInCh.y() / fabs(dirZInCh.z()), -1.);
0169 
0170     seg->setPosition(LocalPoint(posPhiInCh.x(), posZAt0.y(), 0.));
0171     seg->setDirection(dir.unit());
0172 
0173     AlgebraicSymMatrix mat(4);
0174 
0175     // set cov matrix
0176     mat[0][0] = segPhi->parametersError()[0][0];  //sigma (dx/dz)
0177     mat[0][2] = segPhi->parametersError()[0][1];  //cov(dx/dz,x)
0178     mat[2][2] = segPhi->parametersError()[1][1];  //sigma (x)
0179 
0180     seg->setCovMatrix(mat);
0181     seg->setCovMatrixForZed(posZInCh);
0182   } else if (seg->hasPhi()) {
0183     DTChamberRecSegment2D* segPhi = seg->phiSegment();
0184 
0185     seg->setPosition(segPhi->localPosition());
0186     seg->setDirection(segPhi->localDirection());
0187 
0188     AlgebraicSymMatrix mat(4);
0189     // set cov matrix
0190     mat[0][0] = segPhi->parametersError()[0][0];  //sigma (dx/dz)
0191     mat[0][2] = segPhi->parametersError()[0][1];  //cov(dx/dz,x)
0192     mat[2][2] = segPhi->parametersError()[1][1];  //sigma (x)
0193 
0194     seg->setCovMatrix(mat);
0195   } else if (seg->hasZed()) {
0196     DTSLRecSegment2D* segZed = seg->zSegment();
0197 
0198     // Zed seg is in SL one
0199     GlobalPoint glbPosZ = (theGeom->superLayer(segZed->superLayerId()))->toGlobal(segZed->localPosition());
0200     LocalPoint posZInCh = (theGeom->chamber(segZed->superLayerId().chamberId()))->toLocal(glbPosZ);
0201 
0202     GlobalVector glbDirZ = (theGeom->superLayer(segZed->superLayerId()))->toGlobal(segZed->localDirection());
0203     LocalVector dirZInCh = (theGeom->chamber(segZed->superLayerId().chamberId()))->toLocal(glbDirZ);
0204 
0205     LocalPoint posZAt0 = posZInCh + dirZInCh * (-posZInCh.z()) / cos(dirZInCh.theta());
0206 
0207     seg->setPosition(posZAt0);
0208     seg->setDirection(dirZInCh);
0209 
0210     AlgebraicSymMatrix mat(4);
0211     // set cov matrix
0212     seg->setCovMatrix(mat);
0213     seg->setCovMatrixForZed(posZInCh);
0214   }
0215 }
0216 
0217 bool DTSegmentUpdator::fit(DTSegmentCand* seg, bool allow3par, const bool fitdebug) const {
0218   //  if (debug && fitdebug) cout << "[DTSegmentUpdator] Fit DTRecSegment2D" << endl;
0219   //  if (!seg->good()) return false;
0220 
0221   //  DTSuperLayerId DTid = (DTSuperLayerId)seg->superLayer()->id();
0222   //  if (DTid.superlayer()==2)
0223   //    allow3par = 0;
0224 
0225   vector<float> x;
0226   vector<float> y;
0227   vector<float> sigy;
0228   vector<int> lfit;
0229   vector<double> dist;
0230   int i = 0;
0231 
0232   x.reserve(8);
0233   y.reserve(8);
0234   sigy.reserve(8);
0235   lfit.reserve(8);
0236   dist.reserve(8);
0237 
0238   for (DTSegmentCand::AssPointCont::const_iterator iter = seg->hits().begin(); iter != seg->hits().end(); ++iter) {
0239     LocalPoint pos = (*iter).first->localPosition((*iter).second);
0240     float xwire =
0241         (((*iter).first)->localPosition(DTEnums::Left).x() + ((*iter).first)->localPosition(DTEnums::Right).x()) / 2.;
0242     float distance = pos.x() - xwire;
0243 
0244     if ((*iter).second == DTEnums::Left)
0245       lfit.push_back(1);
0246     else
0247       lfit.push_back(-1);
0248 
0249     dist.push_back(distance);
0250     sigy.push_back(sqrt((*iter).first->localPositionError().xx()));
0251     x.push_back(pos.z());
0252     y.push_back(pos.x());
0253     i++;
0254   }
0255 
0256   LocalPoint pos;
0257   LocalVector dir;
0258   AlgebraicSymMatrix covMat(2);
0259   float cminf = 0.;
0260   float vminf = 0.;
0261   double chi2 = 0.;
0262   double t0_corr = 0.;
0263 
0264   fit(x, y, lfit, dist, sigy, pos, dir, cminf, vminf, covMat, chi2, allow3par);
0265   if (cminf != 0)
0266     t0_corr = -cminf / 0.00543;  // convert drift distance to time
0267 
0268   if (debug && fitdebug)
0269     cout << "  DTcand chi2: " << chi2 << "/" << x.size() << "   t0: " << t0_corr << endl;
0270 
0271   seg->setPosition(pos);
0272   seg->setDirection(dir);
0273   seg->sett0(t0_corr);
0274   seg->setCovMatrix(covMat);
0275   seg->setChi2(chi2);
0276 
0277   // cout << "pos " << pos << endl;
0278   // cout << "dir " << dir << endl;
0279   // cout << "Mat " << covMat << endl;
0280 
0281   return true;
0282 }
0283 
0284 void DTSegmentUpdator::fit(DTRecSegment2D* seg, bool allow3par, bool block3par) const {
0285   if (debug)
0286     cout << "[DTSegmentUpdator] Fit DTRecSegment2D - 3par: " << allow3par << endl;
0287 
0288   vector<float> x;
0289   vector<float> y;
0290   vector<float> sigy;
0291   vector<int> lfit;
0292   vector<double> dist;
0293   x.reserve(8);
0294   y.reserve(8);
0295   sigy.reserve(8);
0296   lfit.reserve(8);
0297   dist.reserve(8);
0298 
0299   //  DTSuperLayerId DTid = (DTSuperLayerId)seg->geographicalId();
0300   //  if (DTid.superlayer()==2)
0301   //    allow3par = 0;
0302 
0303   vector<DTRecHit1D> hits = seg->specificRecHits();
0304   for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0305     // I have to get the hits position (the hit is in the layer rf) in SL frame...
0306     GlobalPoint glbPos = (theGeom->layer(hit->wireId().layerId()))->toGlobal(hit->localPosition());
0307     LocalPoint pos = (theGeom->idToDet(seg->geographicalId()))->toLocal(glbPos);
0308     x.push_back(pos.z());
0309     y.push_back(pos.x());
0310 
0311     const DTLayer* layer = theGeom->layer(hit->wireId().layerId());
0312     float xwire = layer->specificTopology().wirePosition(hit->wireId().wire());
0313     float distance = fabs(hit->localPosition().x() - xwire);
0314     dist.push_back(distance);
0315 
0316     int ilc = (hit->lrSide() == DTEnums::Left) ? 1 : -1;
0317     lfit.push_back(ilc);
0318 
0319     // Get local error in SL frame
0320     //RB: is it right in this way?
0321     ErrorFrameTransformer tran;
0322     GlobalError glbErr =
0323         tran.transform(hit->localPositionError(), (theGeom->layer(hit->wireId().layerId()))->surface());
0324     LocalError slErr = tran.transform(glbErr, (theGeom->idToDet(seg->geographicalId()))->surface());
0325     sigy.push_back(sqrt(slErr.xx()));
0326   }
0327 
0328   LocalPoint pos;
0329   LocalVector dir;
0330   AlgebraicSymMatrix covMat(2);
0331   double chi2 = 0.;
0332   float cminf = 0.;
0333   float vminf = 0.;
0334   double t0_corr = 0.;
0335 
0336   fit(x, y, lfit, dist, sigy, pos, dir, cminf, vminf, covMat, chi2, allow3par, block3par);
0337   if (cminf != 0)
0338     t0_corr = -cminf / 0.00543;  // convert drift distance to time
0339 
0340   if (debug)
0341     cout << "   DTSeg2d chi2: " << chi2 << endl;
0342   if (debug)
0343     cout << "   DTSeg2d Fit t0: " << t0_corr << endl;
0344   // cout << "pos " << segPosition << endl;
0345   // cout << "dir " << segDirection << endl;
0346   // cout << "Mat " << mat << endl;
0347 
0348   seg->setPosition(pos);
0349   seg->setDirection(dir);
0350   seg->setCovMatrix(covMat);
0351   seg->setChi2(chi2);
0352   seg->setT0(t0_corr);
0353 }
0354 
0355 void DTSegmentUpdator::fit(const vector<float>& x,
0356                            const vector<float>& y,
0357                            const vector<int>& lfit,
0358                            const vector<double>& dist,
0359                            const vector<float>& sigy,
0360                            LocalPoint& pos,
0361                            LocalVector& dir,
0362                            float& cminf,
0363                            float& vminf,
0364                            AlgebraicSymMatrix& covMatrix,
0365                            double& chi2,
0366                            const bool allow3par,
0367                            const bool block3par) const {
0368   float slope = 0.;
0369   float intercept = 0.;
0370   float covss = 0.;
0371   float covii = 0.;
0372   float covsi = 0.;
0373 
0374   cminf = 0;
0375   vminf = 0;
0376 
0377   int leftHits = 0, rightHits = 0;
0378   for (unsigned int i = 0; i < lfit.size(); i++)
0379     if (lfit[i] == 1)
0380       leftHits++;
0381     else
0382       rightHits++;
0383 
0384   theFitter->fit(x, y, x.size(), sigy, slope, intercept, chi2, covss, covii, covsi);
0385 
0386   // If we have at least one left and one right hit we can try the 3 parameter fit (if it is switched on)
0387   // FIXME: currently the covariance matrix from the 2-par fit is kept
0388   if (leftHits && rightHits && (leftHits + rightHits > 3) && allow3par) {
0389     theFitter->fitNpar(3, x, y, lfit, dist, sigy, slope, intercept, cminf, vminf, chi2, debug);
0390     double t0_corr = -cminf / 0.00543;
0391     if (fabs(t0_corr) < intime_cut && block3par) {
0392       theFitter->fit(x, y, x.size(), sigy, slope, intercept, chi2, covss, covii, covsi);
0393       cminf = 0;
0394     }
0395   }
0396 
0397   // cout << "slope " << slope << endl;
0398   // cout << "intercept " << intercept << endl;
0399 
0400   // intercept is the x() in chamber frame when the segment cross the chamber
0401   // plane (at z()=0), the y() is not measured, so let's put the center of the
0402   // chamber.
0403   pos = LocalPoint(intercept, 0., 0.);
0404 
0405   //  slope is dx()/dz(), while dy()/dz() is by definition 0, finally I want the
0406   //  segment to point outward, so opposite to local z
0407   dir = LocalVector(-slope, 0., -1.).unit();
0408 
0409   covMatrix = AlgebraicSymMatrix(2);
0410   covMatrix[0][0] = covss;  // this is var(dy/dz)
0411   covMatrix[1][1] = covii;  // this is var(y)
0412   covMatrix[1][0] = covsi;  // this is cov(dy/dz,y)
0413 }
0414 
0415 // The GlobalPoint and the GlobalVector can be either the glb position and the direction
0416 // of the 2D-segment itself or the glb position and direction of the 4D segment
0417 void DTSegmentUpdator::updateHits(DTRecSegment2D* seg, GlobalPoint& gpos, GlobalVector& gdir, const int step) const {
0418   // it is not necessary to have DTRecHit1D* to modify the obj in the container
0419   // but I have to be carefully, since I cannot make a copy before the iteration!
0420 
0421   vector<DTRecHit1D> toBeUpdatedRecHits = seg->specificRecHits();
0422   vector<DTRecHit1D> updatedRecHits;
0423 
0424   for (vector<DTRecHit1D>::iterator hit = toBeUpdatedRecHits.begin(); hit != toBeUpdatedRecHits.end(); ++hit) {
0425     const DTLayer* layer = theGeom->layer(hit->wireId().layerId());
0426 
0427     LocalPoint segPos = layer->toLocal(gpos);
0428     LocalVector segDir = layer->toLocal(gdir);
0429 
0430     // define impact angle needed by the step 2
0431     const float angle = atan(segDir.x() / -segDir.z());
0432 
0433     // define the local position (extr.) of the segment. Needed by the third step
0434     LocalPoint segPosAtLayer = segPos + segDir * (-segPos.z()) / cos(segDir.theta());
0435 
0436     DTRecHit1D newHit1D = (*hit);
0437     bool ok = true;
0438 
0439     if (step == 2) {
0440       ok = theAlgo->compute(layer, *hit, angle, newHit1D);
0441 
0442     } else if (step == 3) {
0443       LocalPoint hitPos(hit->localPosition().x(), +segPosAtLayer.y(), 0.);
0444       GlobalPoint glbpos = theGeom->layer(hit->wireId().layerId())->toGlobal(hitPos);
0445       newHit1D.setPosition(hitPos);
0446       ok = theAlgo->compute(layer, *hit, angle, glbpos, newHit1D);
0447 
0448     } else if (step == 4) {
0449       //const double vminf = seg->vDrift();   //  vdrift correction are recorded in the segment
0450       double vminf = 0.;
0451       if (vdrift_4parfit)
0452         vminf = seg->vDrift();  // use vdrift recorded in the segment only if vdrift_4parfit=True
0453 
0454       double cminf = 0.;
0455       if (seg->ist0Valid())
0456         cminf = -seg->t0() * 0.00543;
0457 
0458       //cout << "In updateHits: t0 = " << seg->t0() << endl;
0459       //cout << "In updateHits: vminf = " << vminf << endl;
0460       //cout << "In updateHits: cminf = " << cminf << endl;
0461 
0462       const float xwire = layer->specificTopology().wirePosition(hit->wireId().wire());
0463       const float distance = fabs(hit->localPosition().x() - xwire);
0464       const int ilc = (hit->lrSide() == DTEnums::Left) ? 1 : -1;
0465       const double dy_corr = (vminf * ilc * distance - cminf * ilc);
0466 
0467       LocalPoint point(hit->localPosition().x() + dy_corr, +segPosAtLayer.y(), 0.);
0468 
0469       //double final_hit_resol = T0_hit_resolution;
0470       //if(newHit1D.wireId().layerId().superlayerId().superLayer() != 2) final_hit_resol = final_hit_resol * 0.8;
0471       //LocalError error(final_hit_resol * final_hit_resol,0.,0.);
0472       LocalError error(T0_hit_resolution * T0_hit_resolution, 0., 0.);
0473       newHit1D.setPositionAndError(point, error);
0474 
0475       //FIXME: check that the hit is still inside the cell
0476       ok = true;
0477 
0478     } else
0479       throw cms::Exception("DTSegmentUpdator") << " updateHits called with wrong step " << endl;
0480 
0481     if (ok)
0482       updatedRecHits.push_back(newHit1D);
0483     else {
0484       LogError("DTSegmentUpdator") << "DTSegmentUpdator::updateHits failed update" << endl;
0485       throw cms::Exception("DTSegmentUpdator") << "updateHits failed update" << endl;
0486     }
0487   }
0488   seg->update(updatedRecHits);
0489 }
0490 
0491 void DTSegmentUpdator::rejectBadHits(DTChamberRecSegment2D* phiSeg) const {
0492   vector<float> x;
0493   vector<float> y;
0494 
0495   if (debug)
0496     cout << " Inside the segment updator, now loop on hits:   ( x == z_loc , y == x_loc) " << endl;
0497 
0498   vector<DTRecHit1D> hits = phiSeg->specificRecHits();
0499   const size_t N = hits.size();
0500   if (N < 3)
0501     return;
0502 
0503   for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0504     // I have to get the hits position (the hit is in the layer rf) in SL frame...
0505     GlobalPoint glbPos = (theGeom->layer(hit->wireId().layerId()))->toGlobal(hit->localPosition());
0506     LocalPoint pos = (theGeom->idToDet(phiSeg->geographicalId()))->toLocal(glbPos);
0507 
0508     x.push_back(pos.z());
0509     y.push_back(pos.x());
0510   }
0511 
0512   if (debug) {
0513     cout << " end of segment! " << endl;
0514     cout << " size = Number of Hits: " << x.size() << "  " << y.size() << endl;
0515   }
0516 
0517   // Perform the 2 par fit:
0518   float par[2] = {0., 0.};  // q , m
0519 
0520   //variables to perform the fit:
0521   float Sx = 0.;
0522   float Sy = 0.;
0523   float Sx2 = 0.;
0524   float Sy2 = 0.;
0525   float Sxy = 0.;
0526 
0527   for (size_t i = 0; i < N; ++i) {
0528     Sx += x.at(i);
0529     Sy += y.at(i);
0530     Sx2 += x.at(i) * x.at(i);
0531     Sy2 += y.at(i) * y.at(i);
0532     Sxy += x.at(i) * y.at(i);
0533   }
0534 
0535   const float delta = N * Sx2 - Sx * Sx;
0536   par[0] = (Sx2 * Sy - Sx * Sxy) / delta;
0537   par[1] = (N * Sxy - Sx * Sy) / delta;
0538 
0539   if (debug)
0540     cout << "fit 2 parameters done ----> par0: " << par[0] << "  par1: " << par[1] << endl;
0541 
0542   // Calc residuals:
0543   float residuals[N];
0544   float mean_residual = 0.;  //mean of the absolute values of residuals
0545   for (size_t i = 0; i < N; ++i) {
0546     residuals[i] = y.at(i) - par[1] * x.at(i) - par[0];
0547     mean_residual += std::abs(residuals[i]);
0548     if (debug) {
0549       cout << " i: " << i << " y_i " << y.at(i) << " x_i " << x.at(i) << " res_i " << residuals[i];
0550       if (i == N - 1)
0551         cout << endl;
0552     }
0553   }
0554 
0555   if (debug)
0556     cout << " Residuals computed! " << endl;
0557 
0558   mean_residual = mean_residual / (N - 2);
0559   if (debug)
0560     cout << " mean_residual: " << mean_residual << endl;
0561 
0562   int i = 0;
0563 
0564   // Perform bad hit rejecting -- update hits
0565   vector<DTRecHit1D> updatedRecHits;
0566   for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0567     float normResidual = mean_residual > 0 ? std::abs(residuals[i]) / mean_residual : 0;
0568     ++i;
0569     if (normResidual < 1.5) {
0570       DTRecHit1D newHit1D = (*hit);
0571       updatedRecHits.push_back(newHit1D);
0572       if (debug)
0573         cout << " accepted " << i << "th hit"
0574              << "  Irej: " << normResidual << endl;
0575     } else {
0576       if (debug)
0577         cout << " rejected " << i << "th hit"
0578              << "  Irej: " << normResidual << endl;
0579       continue;
0580     }
0581   }
0582 
0583   phiSeg->update(updatedRecHits);
0584 
0585   //final check!
0586   if (debug) {
0587     vector<float> x_upd;
0588     vector<float> y_upd;
0589 
0590     cout << " Check the update action: " << endl;
0591 
0592     vector<DTRecHit1D> hits_upd = phiSeg->specificRecHits();
0593     for (vector<DTRecHit1D>::const_iterator hit = hits_upd.begin(); hit != hits_upd.end(); ++hit) {
0594       // I have to get the hits position (the hit is in the layer rf) in SL frame...
0595       GlobalPoint glbPos = (theGeom->layer(hit->wireId().layerId()))->toGlobal(hit->localPosition());
0596       LocalPoint pos = (theGeom->idToDet(phiSeg->geographicalId()))->toLocal(glbPos);
0597 
0598       x_upd.push_back(pos.z());
0599       y_upd.push_back(pos.x());
0600 
0601       cout << " x_upd: " << pos.z() << "  y_upd: " << pos.x() << endl;
0602     }
0603 
0604     cout << " end of segment! " << endl;
0605     cout << " size = Number of Hits: " << x_upd.size() << "  " << y_upd.size() << endl;
0606 
0607   }  // end debug
0608 
0609   return;
0610 }  //end DTSegmentUpdator::rejectBadHits
0611 
0612 void DTSegmentUpdator::calculateT0corr(DTRecSegment4D* seg) const {
0613   if (seg->hasPhi())
0614     calculateT0corr(seg->phiSegment());
0615   if (seg->hasZed())
0616     calculateT0corr(seg->zSegment());
0617 }
0618 
0619 void DTSegmentUpdator::calculateT0corr(DTRecSegment2D* seg) const {
0620   // WARNING: since this method is called both with a 2D and a 2DPhi as argument
0621   // seg->geographicalId() can be a superLayerId or a chamberId
0622   if (debug)
0623     cout << "[DTSegmentUpdator] CalculateT0corr DTRecSegment4D" << endl;
0624 
0625   vector<double> d_drift;
0626   vector<float> x;
0627   vector<float> y;
0628   vector<int> lc;
0629 
0630   vector<DTRecHit1D> hits = seg->specificRecHits();
0631 
0632   DTWireId wireId;
0633   int nptfit = 0;
0634 
0635   for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0636     // I have to get the hits position (the hit is in the layer rf) in SL frame...
0637     GlobalPoint glbPos = (theGeom->layer(hit->wireId().layerId()))->toGlobal(hit->localPosition());
0638     LocalPoint pos = (theGeom->idToDet(seg->geographicalId()))->toLocal(glbPos);
0639 
0640     const DTLayer* layer = theGeom->layer(hit->wireId().layerId());
0641     float xwire = layer->specificTopology().wirePosition(hit->wireId().wire());
0642     float distance = fabs(hit->localPosition().x() - xwire);
0643 
0644     int ilc = (hit->lrSide() == DTEnums::Left) ? 1 : -1;
0645 
0646     nptfit++;
0647     x.push_back(pos.z());
0648     y.push_back(pos.x());
0649     lc.push_back(ilc);
0650     d_drift.push_back(distance);
0651 
0652     // cout << " d_drift "<<distance  <<" npt= " <<npt<<endl;
0653   }
0654 
0655   double chi2fit = 0.;
0656   float cminf = 0.;
0657   float vminf = 0.;
0658   float a, b;
0659 
0660   if (nptfit > 2) {
0661     //NB chi2fit is normalized
0662     theFitter->fit4Var(x, y, lc, d_drift, nptfit, a, b, cminf, vminf, chi2fit, vdrift_4parfit, debug);
0663 
0664     double t0cor = -999.;
0665     if (cminf > -998.)
0666       t0cor = -cminf / 0.00543;  // in ns
0667 
0668     //cout << "In calculateT0corr: t0 = " << t0cor << endl;
0669     //cout << "In calculateT0corr: vminf = " << vminf << endl;
0670     //cout << "In calculateT0corr: cminf = " << cminf << endl;
0671     //cout << "In calculateT0corr: chi2 = " << chi2fit << endl;
0672 
0673     seg->setT0(t0cor);      // time  and
0674     seg->setVdrift(vminf);  //  vdrift correction are recorded in the segment
0675   }
0676 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!