DTSegment/src/DTMeantimerPatternReco.cc

0001 /**
0002  *
0003  * \author Stefano Lacaprara - INFN Legnaro <stefano.lacaprara@pd.infn.it>
0004  * \author Riccardo Bellan - INFN TO <riccardo.bellan@cern.ch>
0005  * \author Piotr Traczyk - SINS Warsaw <ptraczyk@fuw.edu.pl>
0006  */
0007
0008 /* This Class Header */
0009 #include "RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco.h"
0010
0011 /* Collaborating Class Header */
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/ConsumesCollector.h"
0015 #include "DataFormats/Common/interface/OwnVector.h"
0016 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0017 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0018 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0019 #include "Geometry/DTGeometry/interface/DTLayer.h"
0020 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0021 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0022 #include "RecoLocalMuon/DTSegment/src/DTSegmentUpdator.h"
0023 #include "RecoLocalMuon/DTSegment/src/DTSegmentCleaner.h"
0024 #include "RecoLocalMuon/DTSegment/src/DTHitPairForFit.h"
0025 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h"
0026 #include "RecoLocalMuon/DTSegment/src/DTLinearFit.h"
0027
0028 /* C++ Headers */
0029 #include <iterator>
0030 using namespace std;
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032
0033 /* ====================================================================== */
0034
0035 typedef std::vector<std::shared_ptr<DTHitPairForFit>> hitCont;
0036 typedef hitCont::const_iterator hitIter;
0037
0038 /// Constructor
0039 DTMeantimerPatternReco::DTMeantimerPatternReco(const edm::ParameterSet& pset, edm::ConsumesCollector cc)
0040     : DTRecSegment2DBaseAlgo(pset),
0041       theFitter(new DTLinearFit()),
0042       theAlgoName("DTMeantimerPatternReco"),
0043       theDTGeometryToken(cc.esConsumes()) {
0044   theMaxAllowedHits = pset.getParameter<unsigned int>("MaxAllowedHits");  // 100
0045   theAlphaMaxTheta = pset.getParameter<double>("AlphaMaxTheta");          // 0.1 ;
0046   theAlphaMaxPhi = pset.getParameter<double>("AlphaMaxPhi");              // 1.0 ;
0047   theMaxChi2 = pset.getParameter<double>("MaxChi2");                      // 8.0 ;
0048   debug = pset.getUntrackedParameter<bool>("debug");
0049   theUpdator = new DTSegmentUpdator(pset, cc);
0050   theCleaner = new DTSegmentCleaner(pset);
0051 }
0052
0053 /// Destructor
0054 DTMeantimerPatternReco::~DTMeantimerPatternReco() {
0055   delete theUpdator;
0056   delete theCleaner;
0057   delete theFitter;
0058 }
0059
0060 /* Operations */
0061 edm::OwnVector<DTSLRecSegment2D> DTMeantimerPatternReco::reconstruct(const DTSuperLayer* sl,
0062                                                                      const std::vector<DTRecHit1DPair>& pairs) {
0063   edm::OwnVector<DTSLRecSegment2D> result;
0064   std::vector<std::shared_ptr<DTHitPairForFit>> hitsForFit = initHits(sl, pairs);
0065
0066   vector<DTSegmentCand*> candidates = buildSegments(sl, hitsForFit);
0067
0068   vector<DTSegmentCand*>::const_iterator cand = candidates.begin();
0069   while (cand < candidates.end()) {
0070     DTSLRecSegment2D* segment = (**cand);
0071     theUpdator->update(segment, true);
0072
0073     if (debug)
0074       cout << "Reconstructed 2D segments " << *segment << endl;
0075     result.push_back(segment);
0076
0077     delete *(cand++);  // delete the candidate!
0078   }
0079
0080   return result;
0081 }
0082
0083 void DTMeantimerPatternReco::setES(const edm::EventSetup& setup) {
0084   // Get the DT Geometry
0085   theDTGeometry = setup.getHandle(theDTGeometryToken);
0086   theUpdator->setES(setup);
0087 }
0088
0089 vector<std::shared_ptr<DTHitPairForFit>> DTMeantimerPatternReco::initHits(const DTSuperLayer* sl,
0090                                                                           const std::vector<DTRecHit1DPair>& hits) {
0091   hitCont result;
0092   for (vector<DTRecHit1DPair>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0093     result.push_back(std::make_shared<DTHitPairForFit>(*hit, *sl, theDTGeometry));
0094   }
0095   return result;
0096 }
0097
0098 vector<DTSegmentCand*> DTMeantimerPatternReco::buildSegments(
0099     const DTSuperLayer* sl, const std::vector<std::shared_ptr<DTHitPairForFit>>& hits) {
0100   vector<DTSegmentCand*> result;
0101   DTEnums::DTCellSide codes[2] = {DTEnums::Left, DTEnums::Right};
0102
0103   if (debug) {
0104     cout << "buildSegments: " << sl->id() << " nHits " << hits.size() << endl;
0105     for (hitIter hit = hits.begin(); hit != hits.end(); ++hit)
0106       cout << **hit << " wire: " << (*hit)->id() << " DigiTime: " << (*hit)->digiTime() << endl;
0107   }
0108
0109   if (hits.size() > theMaxAllowedHits) {
0110     if (debug) {
0111       cout << "Warning: this SuperLayer " << sl->id() << " has too many hits : " << hits.size() << " max allowed is "
0112            << theMaxAllowedHits << endl;
0113       cout << "Skipping segment reconstruction... " << endl;
0114     }
0115     return result;
0116   }
0117
0118   GlobalPoint IP;
0119   float DAlphaMax;
0120   if ((sl->id()).superlayer() == 2)  // Theta SL
0121     DAlphaMax = theAlphaMaxTheta;
0122   else  // Phi SL
0123     DAlphaMax = theAlphaMaxPhi;
0124
0125   // get two hits in different layers and see if there are other hits
0126   //  compatible with them
0127   for (hitCont::const_iterator firstHit = hits.begin(); firstHit != hits.end(); ++firstHit) {
0128     for (hitCont::const_reverse_iterator lastHit = hits.rbegin(); (*lastHit) != (*firstHit); ++lastHit) {
0129       // a geometrical sensibility cut for the two hits
0130       if (!geometryFilter((*firstHit)->id(), (*lastHit)->id()))
0131         continue;
0132
0133       // create a set of hits for the fit (only the hits between the two selected ones)
0134       hitCont hitsForFit;
0135       for (hitCont::const_iterator tmpHit = firstHit + 1; (*tmpHit) != (*lastHit); tmpHit++)
0136         if ((geometryFilter((*tmpHit)->id(), (*lastHit)->id())) && (geometryFilter((*tmpHit)->id(), (*firstHit)->id())))
0137           hitsForFit.push_back(*tmpHit);
0138
0139       for (int firstLR = 0; firstLR < 2; ++firstLR) {
0140         for (int lastLR = 0; lastLR < 2; ++lastLR) {
0141           // TODO move the global transformation in the DTHitPairForFit class
0142           // when it will be moved I will able to remove the sl from the input parameter
0143           GlobalPoint gposFirst = sl->toGlobal((*firstHit)->localPosition(codes[firstLR]));
0144           GlobalPoint gposLast = sl->toGlobal((*lastHit)->localPosition(codes[lastLR]));
0145           GlobalVector gvec = gposLast - gposFirst;
0146           GlobalVector gvecIP = gposLast - IP;
0147
0148           // difference in angle measured
0149           float DAlpha = fabs(gvec.theta() - gvecIP.theta());
0150           if (DAlpha > DAlphaMax)
0151             continue;
0152
0153           //      if(debug) {
0154           //              cout << "Selected hit pair:" << endl;
0155           //              cout << "  First " << *(*firstHit) << " Layer Id: " << (*firstHit)->id().layerId() << " Side: " << firstLR << " DigiTime: " << (*firstHit)->digiTime() << endl;
0156           //              cout << "  Last "  << *(*lastHit)  << " Layer Id: " << (*lastHit)->id().layerId()  << " Side: " << lastLR << " DigiTime: " << (*lastHit)->digiTime() <<  endl;
0157           //          }
0158
0159           DTSegmentCand::AssPointCont pointSet;
0160           auto segCand = std::make_unique<DTSegmentCand>(pointSet, sl);
0161           segCand->add(*firstHit, codes[firstLR]);
0162           segCand->add(*lastHit, codes[lastLR]);
0163
0164           // run hit adding/segment building
0165           maxfound = 3;
0166           addHits(segCand.get(), hitsForFit, result);
0167         }
0168       }
0169     }
0170   }
0171
0172   // now I have a couple of segment hypotheses, should check for ghosts
0173   if (debug) {
0174     cout << "Result (before cleaning): " << result.size() << endl;
0175     for (vector<DTSegmentCand*>::const_iterator seg = result.begin(); seg != result.end(); ++seg)
0176       cout << *(*seg) << endl;
0177   }
0178
0179   result = theCleaner->clean(result);
0180
0181   if (debug) {
0182     cout << "Result (after cleaning): " << result.size() << endl;
0183     for (vector<DTSegmentCand*>::const_iterator seg = result.begin(); seg != result.end(); ++seg)
0184       cout << *(*seg) << endl;
0185   }
0186
0187   return result;
0188 }
0189
0190 void DTMeantimerPatternReco::addHits(DTSegmentCand* segCand,
0191                                      const vector<std::shared_ptr<DTHitPairForFit>>& hits,
0192                                      vector<DTSegmentCand*>& result) {
0193   double chi2l, chi2r, t0l, t0r;
0194   bool foundSomething = false;
0195
0196   if (debug)
0197     cout << " DTMeantimerPatternReco::addHit " << endl
0198          << "   Picked " << segCand->nHits() << " hits, " << hits.size() << " left." << endl;
0199
0200   if (segCand->nHits() + hits.size() < maxfound)
0201     return;
0202
0203   // loop over the remaining hits
0204   for (hitCont::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0205     //    if (debug) {
0206     //      cout << "     Trying B: " << **hit<< " wire: " << (*hit)->id() << endl;
0207     //      printPattern(assHits,*hit);
0208     //    }
0209
0210     DTSegmentCand::AssPoint lhit(*hit, DTEnums::Left);
0211     DTSegmentCand::AssPoint rhit(*hit, DTEnums::Right);
0212
0213     segCand->add(lhit);
0214     bool left_ok = (fitWithT0(segCand, false) ? true : false);
0215     chi2l = segCand->chi2();
0216     t0l = segCand->t0();
0217     segCand->removeHit(lhit);
0218
0219     segCand->add(rhit);
0220     bool right_ok = (fitWithT0(segCand, false) ? true : false);
0221     chi2r = segCand->chi2();
0222     t0r = segCand->t0();
0223     segCand->removeHit(rhit);
0224
0225     if (debug) {
0226       int nHits = segCand->nHits() + 1;
0227       cout << "    Left:  t0= " << t0l << "  chi2/nHits= " << chi2l << "/" << nHits << "  ok: " << left_ok << endl;
0228       cout << "   Right:  t0= " << t0r << "  chi2/nHits= " << chi2r << "/" << nHits << "  ok: " << right_ok << endl;
0229     }
0230
0231     if (!left_ok && !right_ok)
0232       continue;
0233
0234     foundSomething = true;
0235
0236     // prepare the hit set for the next search, start from the other side
0237     hitCont hitsForFit;
0238     for (hitCont::const_iterator tmpHit = hit + 1; tmpHit != hits.end(); tmpHit++)
0239       if (geometryFilter((*tmpHit)->id(), (*hit)->id()))
0240         hitsForFit.push_back(*tmpHit);
0241
0242     reverse(hitsForFit.begin(), hitsForFit.end());
0243
0244     // choose only one - left or right
0245     if (segCand->nHits() > 3 && left_ok && right_ok) {
0246       if (chi2l < chi2r - 0.1)
0247         right_ok = false;
0248       else if (chi2r < chi2l - 0.1)
0249         left_ok = false;
0250     }
0251
0252     if (left_ok) {
0253       segCand->add(lhit);
0254       addHits(segCand, hitsForFit, result);
0255       segCand->removeHit(lhit);
0256     }
0257
0258     if (right_ok) {
0259       segCand->add(rhit);
0260       addHits(segCand, hitsForFit, result);
0261       segCand->removeHit(rhit);
0262     }
0263   }
0264
0265   if (foundSomething)
0266     return;
0267   // if we didn't find any new hits compatible with the current candidate, we proceed to check and store the candidate
0268
0269   // If we already have a segment with more hits from this hit pair - don't save this one.
0270   if (segCand->nHits() < maxfound)
0271     return;
0272
0273   // Check if semgent Ok, calculate chi2
0274   bool seg_ok = (fitWithT0(segCand, debug) ? true : false);
0275   if (!seg_ok)
0276     return;
0277
0278   if (!segCand->good()) {
0279     //    if (debug) cout << "   Segment not good() - skipping" << endl;
0280     return;
0281   }
0282
0283   if (segCand->nHits() > maxfound)
0284     maxfound = segCand->nHits();
0285   if (debug)
0286     cout << endl << "   Seg t0= " << segCand->t0() << endl << *segCand << endl;
0287
0288   if (checkDoubleCandidates(result, segCand)) {
0289     result.push_back(new DTSegmentCand(*segCand));
0290     if (debug)
0291       cout << "   Result is now " << result.size() << endl;
0292   } else {
0293     if (debug)
0294       cout << "   Exists - skipping" << endl;
0295   }
0296 }
0297
0298 bool DTMeantimerPatternReco::geometryFilter(const DTWireId first, const DTWireId second) const {
0299   //  return true;
0300
0301   const int layerLowerCut[4] = {0, -1, -2, -2};
0302   const int layerUpperCut[4] = {0, 2, 2, 3};
0303   //  const int layerLowerCut[4]={0,-2,-4,-5};
0304   //  const int layerUpperCut[4]={0, 3, 4, 6};
0305
0306   // deal only with hits that are in the same SL
0307   if (first.layerId().superlayerId().superLayer() != second.layerId().superlayerId().superLayer())
0308     return true;
0309
0310   int deltaLayer = abs(first.layerId().layer() - second.layerId().layer());
0311
0312   // drop hits in the same layer
0313   if (!deltaLayer)
0314     return false;
0315
0316   // protection against unexpected layer numbering
0317   if (deltaLayer > 3) {
0318     cout << "*** WARNING! DT Layer numbers differ by more than 3! for hits: " << endl;
0319     cout << "             " << first << endl;
0320     cout << "             " << second << endl;
0321     return false;
0322   }
0323
0324   // accept only hits in cells "not too far away"
0325   int deltaWire = first.wire() - second.wire();
0326   if (second.layerId().layer() % 2 == 0)
0327     deltaWire = -deltaWire;  // yet another trick to get it right...
0328   if ((deltaWire <= layerLowerCut[deltaLayer]) || (deltaWire >= layerUpperCut[deltaLayer]))
0329     return false;
0330
0331   return true;
0332 }
0333
0334 DTSegmentCand* DTMeantimerPatternReco::fitWithT0(DTSegmentCand* seg, const bool fitdebug) {
0335   // perform the 3 parameter fit on the segment candidate
0336   theUpdator->fit(seg, true, fitdebug);
0337   double chi2 = seg->chi2();
0338
0339   // Sanity check - drop segment candidates with a failed 3-par fit.
0340   // (this includes segments with hits after the calculated t0 correction ending up
0341   // beyond the chamber walls or on the other side of the wire)
0342   if (chi2 == -1.)
0343     return nullptr;
0344
0345   // at this point we keep all 3-hit segments that passed the above check
0346   if (seg->nHits() == 3)
0347     return seg;
0348
0349   // for segments with no t0 information we impose a looser chi2 cut
0350   if (seg->t0() == 0) {
0351     if (chi2 < 100.)
0352       return seg;
0353     else
0354       return nullptr;
0355   }
0356
0357   // cut on chi2/ndof of the segment candidate
0358   if ((chi2 / (seg->nHits() - 3) < theMaxChi2))
0359     return seg;
0360   else
0361     return nullptr;
0362 }
0363
0364 bool DTMeantimerPatternReco::checkDoubleCandidates(vector<DTSegmentCand*>& cands, DTSegmentCand* seg) {
0365   for (vector<DTSegmentCand*>::iterator cand = cands.begin(); cand != cands.end(); ++cand) {
0366     if (*(*cand) == *seg)
0367       return false;
0368     if (((*cand)->nHits() >= seg->nHits()) && ((*cand)->chi2ndof() < seg->chi2ndof()))
0369       if ((*cand)->nSharedHitPairs(*seg) > int(seg->nHits() - 2))
0370         return false;
0371   }
0372   return true;
0373 }
0374
0375 void DTMeantimerPatternReco::printPattern(vector<DTSegmentCand::AssPoint>& assHits, const DTHitPairForFit* hit) {
0376   char mark[26] = {". . . . . . . . . . . . "};
0377   int wire[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
0378
0379   for (vector<DTSegmentCand::AssPoint>::const_iterator assHit = assHits.begin(); assHit != assHits.end(); ++assHit) {
0380     int lay =
0381         (((*assHit).first)->id().superlayerId().superLayer() - 1) * 4 + ((*assHit).first)->id().layerId().layer() - 1;
0382     wire[lay] = ((*assHit).first)->id().wire();
0383     if ((*assHit).second == DTEnums::Left)
0384       mark[lay * 2] = 'L';
0385     else
0386       mark[lay * 2] = 'R';
0387   }
0388
0389   int lay = ((*hit).id().superlayerId().superLayer() - 1) * 4 + (*hit).id().layerId().layer() - 1;
0390   wire[lay] = (*hit).id().wire();
0391   mark[lay * 2] = '*';
0392
0393   cout << "   " << mark << endl << "  ";
0394   for (int i = 0; i < 12; i++)
0395     if (wire[i])
0396       cout << setw(2) << wire[i];
0397     else
0398       cout << "  ";
0399   cout << endl;
0400 }