DTSegment/src/DTSegmentCand.cc

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  */
0006
0007 /* This Class Header */
0008 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h"
0009
0010 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0011 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0012 #include "DataFormats/DTRecHit/interface/DTChamberRecSegment2D.h"
0013
0014 /* Collaborating Class Header */
0015
0016 /* C++ Headers */
0017
0018 /* ====================================================================== */
0019 const double DTSegmentCand::chi2max = 20.;       // to be tuned!!
0020 const unsigned int DTSegmentCand::nHitsMin = 3;  // to be tuned!!
0021
0022 /// Constructor
0023 DTSegmentCand::DTSegmentCand(AssPointCont& hits, const DTSuperLayer* sl) : theSL(sl), theChi2(-1.), theHits(hits) {}
0024
0025 DTSegmentCand::DTSegmentCand(const AssPointCont& hits,
0026                              LocalPoint& position,
0027                              LocalVector& direction,
0028                              double chi2,
0029                              const AlgebraicSymMatrix& covMat,
0030                              const DTSuperLayer* sl)
0031     : theSL(sl), thePosition(position), theDirection(direction), theChi2(chi2), theCovMatrix(covMat), theHits(hits) {}
0032
0033 /// Destructor
0034 DTSegmentCand::~DTSegmentCand() {}
0035
0036 /* Operations */
0037 bool DTSegmentCand::operator==(const DTSegmentCand& seg) const {
0038   static const double epsilon = 0.00001;
0039   if (nHits() != seg.nHits())
0040     return false;
0041   if (fabs(chi2() - seg.chi2()) > epsilon)
0042     return false;
0043   if (fabs(position().x() - seg.position().x()) > epsilon || fabs(position().y() - seg.position().y()) > epsilon ||
0044       fabs(position().z() - seg.position().z()) > epsilon)
0045     return false;
0046   if (fabs(direction().x() - seg.direction().x()) > epsilon || fabs(direction().y() - seg.direction().y()) > epsilon ||
0047       fabs(direction().z() - seg.direction().z()) > epsilon)
0048     return false;
0049   return true;
0050 }
0051
0052 bool DTSegmentCand::operator<(const DTSegmentCand& seg) const {
0053   if (nHits() == seg.nHits())
0054     return (chi2() > seg.chi2());
0055   return (nHits() < seg.nHits());
0056 }
0057
0058 void DTSegmentCand::add(AssPoint newHit) { theHits.insert(newHit); }
0059
0060 void DTSegmentCand::add(std::shared_ptr<DTHitPairForFit> hit, DTEnums::DTCellSide code) {
0061   AssPoint newHit(hit, code);
0062   theHits.insert(newHit);
0063 }
0064
0065 void DTSegmentCand::removeHit(AssPoint badHit) { theHits.erase(badHit); }
0066
0067 int DTSegmentCand::nSharedHitPairs(const DTSegmentCand& seg) const {
0068   int result = 0;
0069
0070   for (AssPointCont::const_iterator hit = theHits.begin(); hit != theHits.end(); ++hit) {
0071     for (AssPointCont::const_iterator hit2 = seg.hits().begin(); hit2 != seg.hits().end(); ++hit2) {
0072       //  if(result) return result ; // TODO, uncomm this line or move it in another func
0073       if ((*(*hit).first) == (*(*hit2).first)) {
0074         ++result;
0075         continue;
0076       }
0077     }
0078   }
0079   return result;
0080 }
0081
0082 DTSegmentCand::AssPointCont DTSegmentCand::conflictingHitPairs(const DTSegmentCand& seg) const {
0083   AssPointCont result;
0084   const AssPointCont& hits2 = seg.theHits;
0085
0086   //  if (nSharedHitPairs(seg)==0) return result;
0087
0088   AssPointCont::const_iterator hitBegin2 = hits2.begin(), hitEnd2 = hits2.end();
0089   for (AssPointCont::const_iterator hit = theHits.begin(), hitEnd = theHits.end(); hit != hitEnd; ++hit) {
0090     for (AssPointCont::const_iterator hit2 = hitBegin2; hit2 != hitEnd2; ++hit2) {
0091       if ((*(*hit).first) == (*(*hit2).first) && (*hit).second != (*hit2).second) {
0092         result.insert(*hit);
0093         continue;
0094       }
0095     }
0096   }
0097   return result;
0098 }
0099
0100 bool DTSegmentCand::good() const {
0101   // std::cout << NDOF() << "  " << chi2()/NDOF() << "   " << nHits() << std::endl;
0102   if (NDOF() == 0)
0103     return false;
0104   if (chi2() / NDOF() > chi2max || nHits() < nHitsMin)
0105     return false;
0106
0107   if (nHits() == nHitsMin && hitsShareLayer())
0108     return false;
0109
0110   return true;
0111 }
0112
0113 bool DTSegmentCand::hitsShareLayer() const {
0114   const unsigned int hitsSize = theHits.size();
0115   // we don't expect so many 1D hits, if such a segment arrives just drop it
0116   if (hitsSize > 20)
0117     return false;
0118
0119   int layerN[hitsSize];
0120   unsigned int i = 0;
0121   for (DTSegmentCand::AssPointCont::iterator assHit = theHits.begin(); assHit != theHits.end(); ++assHit) {
0122     layerN[i] = (*assHit).first->id().layerId().layer() + 10 * (*assHit).first->id().superlayerId().superlayer();
0123     for (unsigned int j = 0; j < i; j++) {
0124       if (layerN[i] == layerN[j])
0125         return true;
0126     }
0127     i++;
0128   }
0129
0130   return false;
0131 }
0132
0133 int DTSegmentCand::nLayers() const {
0134   // TODO
0135   return 0;
0136 }
0137
0138 DTSegmentCand::operator DTSLRecSegment2D*() const {
0139   LocalPoint seg2Dposition = position();
0140   LocalVector seg2DDirection = direction();
0141   double seg2DChi2 = chi2();
0142   AlgebraicSymMatrix seg2DCovMatrix = covMatrix();
0143
0144   std::vector<DTRecHit1D> hits1D;
0145   for (DTSegmentCand::AssPointCont::iterator assHit = theHits.begin(); assHit != theHits.end(); ++assHit) {
0146     GlobalPoint hitGlobalPos = theSL->toGlobal((*assHit).first->localPosition((*assHit).second));
0147
0148     LocalPoint hitPosInLayer = theSL->layer((*assHit).first->id().layerId())->toLocal(hitGlobalPos);
0149
0150     DTRecHit1D hit(((*assHit).first)->id(),
0151                    (*assHit).second,
0152                    ((*assHit).first)->digiTime(),
0153                    hitPosInLayer,
0154                    ((*assHit).first)->localPositionError());
0155     hits1D.push_back(hit);
0156   }
0157
0158   return new DTSLRecSegment2D(theSL->id(), seg2Dposition, seg2DDirection, seg2DCovMatrix, seg2DChi2, hits1D);
0159 }
0160
0161 DTSegmentCand::operator DTChamberRecSegment2D*() const {
0162   // input position and direction are in sl frame, while must be stored in
0163   // chamber one: so I have to extrapolate the position (along the direction) to
0164   // the chamber reference plane.
0165
0166   LocalPoint posInCh = theSL->chamber()->toLocal(theSL->toGlobal(position()));
0167   LocalVector dirInCh = theSL->chamber()->toLocal(theSL->toGlobal(direction()));
0168
0169   LocalPoint pos = posInCh + dirInCh * posInCh.z() / cos(dirInCh.theta());
0170
0171   double seg2DChi2 = chi2();
0172   AlgebraicSymMatrix seg2DCovMatrix = covMatrix();
0173
0174   std::vector<DTRecHit1D> hits1D;
0175   for (DTSegmentCand::AssPointCont::iterator assHit = theHits.begin(); assHit != theHits.end(); ++assHit) {
0176     GlobalPoint hitGlobalPos = theSL->toGlobal((*assHit).first->localPosition((*assHit).second));
0177
0178     LocalPoint hitPosInLayer = theSL->chamber()
0179                                    ->superLayer((*assHit).first->id().superlayerId())
0180                                    ->layer((*assHit).first->id().layerId())
0181                                    ->toLocal(hitGlobalPos);
0182
0183     DTRecHit1D hit(((*assHit).first)->id(),
0184                    (*assHit).second,
0185                    ((*assHit).first)->digiTime(),
0186                    hitPosInLayer,
0187                    ((*assHit).first)->localPositionError());
0188     hits1D.push_back(hit);
0189   }
0190
0191   return new DTChamberRecSegment2D(theSL->chamber()->id(), pos, dirInCh, seg2DCovMatrix, seg2DChi2, hits1D);
0192
0193   // chamber and Phi SLs' frame are oriented in the same way, only a transaltion,
0194   // so the covariance matrix is the same!
0195 }
0196
0197 bool DTSegmentCand::AssPointLessZ::operator()(const AssPoint& pt1, const AssPoint& pt2) const {
0198   return *(pt1.first) < *(pt2.first);
0199 }
0200
0201 std::ostream& operator<<(std::ostream& out, const DTSegmentCand& seg) {
0202   out << " pos: " << seg.position() << " dir: " << seg.direction() << " chi2/nHits: " << seg.chi2() << "/"
0203       << seg.DTSegmentCand::nHits() << " t0: " << seg.t0();
0204   return out;
0205 }
0206
0207 std::ostream& operator<<(std::ostream& out, const DTSegmentCand::AssPoint& hit) {
0208   // out << "Hits " << (hit.first)->localPosition(DTEnums::Left) <<
0209   //     " " << hit.second  << " Lay " << (hit.first)->layerNumber() << endl;
0210   return out;
0211 }