DTSegment/src/DTCombinatorialExtendedPatternReco.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/DTCombinatorialExtendedPatternReco.h"
0009
0010 /* Collaborating Class Header */
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "DataFormats/Common/interface/OwnVector.h"
0014 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0015 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0016 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0017 #include "Geometry/DTGeometry/interface/DTLayer.h"
0018 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0019 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0020 #include "RecoLocalMuon/DTSegment/src/DTSegmentUpdator.h"
0021 #include "RecoLocalMuon/DTSegment/src/DTSegmentCleaner.h"
0022 #include "RecoLocalMuon/DTSegment/src/DTHitPairForFit.h"
0023 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h"
0024 #include "RecoLocalMuon/DTSegment/src/DTSegmentExtendedCand.h"
0025
0026 /* C++ Headers */
0027 #include <iterator>
0028 using namespace std;
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030
0031 /* ====================================================================== */
0032
0033 /// Constructor
0034 DTCombinatorialExtendedPatternReco::DTCombinatorialExtendedPatternReco(const edm::ParameterSet& pset,
0035                                                                        edm::ConsumesCollector cc)
0036     : DTRecSegment2DBaseAlgo(pset),
0037       theAlgoName("DTCombinatorialExtendedPatternReco"),
0038       theDTGeometryToken(cc.esConsumes()) {
0039   theMaxAllowedHits = pset.getParameter<unsigned int>("MaxAllowedHits");  // 100
0040   theAlphaMaxTheta = pset.getParameter<double>("AlphaMaxTheta");          // 0.1 ;
0041   theAlphaMaxPhi = pset.getParameter<double>("AlphaMaxPhi");              // 1.0 ;
0042   debug = pset.getUntrackedParameter<bool>("debug");                      //true;
0043   theUpdator = new DTSegmentUpdator(pset, cc);
0044   theCleaner = new DTSegmentCleaner(pset);
0045   string theHitAlgoName = pset.getParameter<string>("recAlgo");
0046   usePairs = !(theHitAlgoName == "DTNoDriftAlgo");
0047 }
0048
0049 /// Destructor
0050 DTCombinatorialExtendedPatternReco::~DTCombinatorialExtendedPatternReco() {}
0051
0052 /* Operations */
0053 edm::OwnVector<DTSLRecSegment2D> DTCombinatorialExtendedPatternReco::reconstruct(
0054     const DTSuperLayer* sl, const std::vector<DTRecHit1DPair>& pairs) {
0055   if (debug)
0056     cout << "DTCombinatorialExtendedPatternReco::reconstruct" << endl;
0057   theTriedPattern.clear();
0058   edm::OwnVector<DTSLRecSegment2D> result;
0059   vector<std::shared_ptr<DTHitPairForFit>> hitsForFit = initHits(sl, pairs);
0060
0061   vector<DTSegmentCand*> candidates = buildSegments(sl, hitsForFit);
0062
0063   vector<DTSegmentCand*>::const_iterator cand = candidates.begin();
0064   while (cand < candidates.end()) {
0065     DTSLRecSegment2D* segment = (**cand);
0066
0067     theUpdator->update(segment, false);
0068
0069     result.push_back(segment);
0070
0071     if (debug) {
0072       cout << "Reconstructed 2D extended segments " << result.back() << endl;
0073     }
0074
0075     delete *(cand++);  // delete the candidate!
0076   }
0077
0078   return result;
0079 }
0080
0081 void DTCombinatorialExtendedPatternReco::setES(const edm::EventSetup& setup) {
0082   // Get the DT Geometry
0083   theDTGeometry = setup.getHandle(theDTGeometryToken);
0084   theUpdator->setES(setup);
0085 }
0086
0087 void DTCombinatorialExtendedPatternReco::setClusters(const vector<DTSLRecCluster>& clusters) { theClusters = clusters; }
0088
0089 vector<std::shared_ptr<DTHitPairForFit>> DTCombinatorialExtendedPatternReco::initHits(
0090     const DTSuperLayer* sl, const std::vector<DTRecHit1DPair>& hits) {
0091   vector<std::shared_ptr<DTHitPairForFit>> 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*> DTCombinatorialExtendedPatternReco::buildSegments(
0099     const DTSuperLayer* sl, const std::vector<std::shared_ptr<DTHitPairForFit>>& hits) {
0100   typedef vector<std::shared_ptr<DTHitPairForFit>> hitCont;
0101   typedef hitCont::const_iterator hitIter;
0102   vector<DTSegmentCand*> result;
0103
0104   if (debug) {
0105     cout << "DTCombinatorialExtendedPatternReco::buildSegments: " << sl->id() << " nHits " << hits.size() << endl;
0106     for (vector<std::shared_ptr<DTHitPairForFit>>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit)
0107       cout << **hit << endl;
0108   }
0109
0110   // 10-Mar-2004 SL
0111   // put a protection against heavily populated chambers, for which the segment
0112   // building could lead to infinite memory usage...
0113   if (hits.size() > theMaxAllowedHits) {
0114     if (debug) {
0115       cout << "Warning: this SuperLayer " << sl->id() << " has too many hits : " << hits.size() << " max allowed is "
0116            << theMaxAllowedHits << endl;
0117       cout << "Skipping segment reconstruction... " << endl;
0118     }
0119     return result;
0120   }
0121
0122   /// get two hits in different layers and see if there are other / hits
0123   //  compatible with them
0124   for (hitCont::const_iterator firstHit = hits.begin(); firstHit != hits.end(); ++firstHit) {
0125     for (hitCont::const_reverse_iterator lastHit = hits.rbegin(); (*lastHit) != (*firstHit); ++lastHit) {
0126       // hits must nor in the same nor in adiacent layers
0127       if (fabs((*lastHit)->id().layerId() - (*firstHit)->id().layerId()) <= 1)
0128         continue;
0129       if (debug) {
0130         cout << "Selected these two hits pair " << endl;
0131         cout << "First " << *(*firstHit) << " Layer Id: " << (*firstHit)->id().layerId() << endl;
0132         cout << "Last " << *(*lastHit) << " Layer Id: " << (*lastHit)->id().layerId() << endl;
0133       }
0134
0135       GlobalPoint IP;
0136       float DAlphaMax;
0137       if ((sl->id()).superlayer() == 2)  // Theta SL
0138         DAlphaMax = theAlphaMaxTheta;
0139       else  // Phi SL
0140         DAlphaMax = theAlphaMaxPhi;
0141
0142       DTEnums::DTCellSide codes[2] = {DTEnums::Right, DTEnums::Left};
0143       for (int firstLR = 0; firstLR < 2; ++firstLR) {
0144         for (int lastLR = 0; lastLR < 2; ++lastLR) {
0145           // TODO move the global transformation in the DTHitPairForFit class
0146           // when it will be moved I will able to remove the sl from the input parameter
0147           GlobalPoint gposFirst = sl->toGlobal((*firstHit)->localPosition(codes[firstLR]));
0148           GlobalPoint gposLast = sl->toGlobal((*lastHit)->localPosition(codes[lastLR]));
0149
0150           GlobalVector gvec = gposLast - gposFirst;
0151           GlobalVector gvecIP = gposLast - IP;
0152
0153           // difference in angle measured
0154           float DAlpha = fabs(gvec.theta() - gvecIP.theta());
0155
0156           // cout << "DAlpha " << DAlpha << endl;
0157           if (DAlpha < DAlphaMax) {
0158             // create a segment hypotesis
0159             // I don't need a true segment, just direction and position
0160             LocalPoint posIni = (*firstHit)->localPosition(codes[firstLR]);
0161             LocalVector dirIni = ((*lastHit)->localPosition(codes[lastLR]) - posIni).unit();
0162
0163             // search for other compatible hits, with or without the L/R solved
0164             vector<DTSegmentCand::AssPoint> assHits = findCompatibleHits(posIni, dirIni, hits);
0165             if (debug)
0166               cout << "compatible hits " << assHits.size() << endl;
0167
0168             // here return an extended candidate (which _has_ the original
0169             // segment)
0170             DTSegmentExtendedCand* seg = buildBestSegment(assHits, sl);
0171
0172             if (seg) {
0173               if (debug)
0174                 cout << "segment " << *seg << endl;
0175
0176               // check if the chi2 and #hits are ok
0177               if (!seg->good()) {  // good is reimplmented in extended segment
0178                 delete seg;
0179               } else {
0180                 // remove duplicated segments
0181                 if (checkDoubleCandidates(result, seg)) {
0182                   // add to the vector of hypotesis
0183                   // still work with extended segments
0184                   result.push_back(seg);
0185                   if (debug)
0186                     cout << "result is now " << result.size() << endl;
0187                 } else {  // delete it!
0188                   delete seg;
0189                   if (debug)
0190                     cout << "already existing" << endl;
0191                 }
0192               }
0193             }
0194           }
0195         }
0196       }
0197     }
0198   }
0199   if (debug) {
0200     for (vector<DTSegmentCand*>::const_iterator seg = result.begin(); seg != result.end(); ++seg)
0201       cout << *(*seg) << endl;
0202   }
0203
0204   // now I have a couple of segment hypotesis, should check for ghost
0205   // still with extended candidates
0206   result = theCleaner->clean(result);
0207   if (debug) {
0208     cout << "result no ghost  " << result.size() << endl;
0209     for (vector<DTSegmentCand*>::const_iterator seg = result.begin(); seg != result.end(); ++seg)
0210       cout << *(*seg) << endl;
0211   }
0212
0213   // here, finally, I have to return the set of _original_ segments, not the
0214   // extended ones.
0215   return result;
0216 }
0217
0218 vector<DTSegmentCand::AssPoint> DTCombinatorialExtendedPatternReco::findCompatibleHits(
0219     const LocalPoint& posIni, const LocalVector& dirIni, const vector<std::shared_ptr<DTHitPairForFit>>& hits) {
0220   if (debug)
0221     cout << "Pos: " << posIni << " Dir: " << dirIni << endl;
0222   vector<DTSegmentCand::AssPoint> result;
0223
0224   // counter to early-avoid double counting in hits pattern
0225   vector<int> tried;
0226   int nCompatibleHits = 0;
0227
0228   typedef vector<std::shared_ptr<DTHitPairForFit>> hitCont;
0229   typedef hitCont::const_iterator hitIter;
0230   for (hitIter hit = hits.begin(); hit != hits.end(); ++hit) {
0231     pair<bool, bool> isCompatible = (*hit)->isCompatible(posIni, dirIni);
0232     if (debug)
0233       cout << "isCompatible " << isCompatible.first << " " << isCompatible.second << endl;
0234
0235     // if only one of the two is compatible, then the LR is assigned,
0236     // otherwise is undefined
0237
0238     DTEnums::DTCellSide lrcode;
0239     if (isCompatible.first && isCompatible.second) {
0240       usePairs ? lrcode = DTEnums::undefLR : lrcode = DTEnums::Left;  // if not usePairs then only use single side
0241       tried.push_back(3);
0242       nCompatibleHits++;
0243     } else if (isCompatible.first) {
0244       lrcode = DTEnums::Left;
0245       tried.push_back(2);
0246       nCompatibleHits++;
0247     } else if (isCompatible.second) {
0248       lrcode = DTEnums::Right;
0249       tried.push_back(1);
0250       nCompatibleHits++;
0251     } else {
0252       tried.push_back(0);
0253       continue;  // neither is compatible
0254     }
0255     result.push_back(DTSegmentCand::AssPoint(*hit, lrcode));
0256   }
0257
0258   // check if too few associated hits or pattern already tried
0259   if (nCompatibleHits < 3 || find(theTriedPattern.begin(), theTriedPattern.end(), tried) == theTriedPattern.end()) {
0260     theTriedPattern.push_back(tried);
0261   } else {
0262     if (debug) {
0263       vector<vector<int>>::const_iterator t = find(theTriedPattern.begin(), theTriedPattern.end(), tried);
0264       cout << "Already tried";
0265       copy((*t).begin(), (*t).end(), ostream_iterator<int>(std::cout));
0266       cout << endl;
0267     }
0268     // empty the result vector
0269     result.clear();
0270   }
0271   return result;
0272 }
0273
0274 DTSegmentExtendedCand* DTCombinatorialExtendedPatternReco::buildBestSegment(std::vector<DTSegmentCand::AssPoint>& hits,
0275                                                                             const DTSuperLayer* sl) {
0276   if (debug)
0277     cout << "DTCombinatorialExtendedPatternReco::buildBestSegment " << hits.size() << endl;
0278   if (hits.size() < 3) {
0279     //cout << "buildBestSegment: hits " << hits.size()<< endl;
0280     return nullptr;  // a least 3 point
0281   }
0282
0283   // hits with defined LR
0284   vector<DTSegmentCand::AssPoint> points;
0285
0286   // without: I store both L and R, a deque since I need front insertion and
0287   // deletion
0288   deque<std::shared_ptr<DTHitPairForFit>> pointsNoLR;
0289
0290   // first add only the hits with LR assigned
0291   for (vector<DTSegmentCand::AssPoint>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0292     if ((*hit).second != DTEnums::undefLR) {
0293       points.push_back(*hit);
0294     } else {  // then also for the undef'd one
0295       pointsNoLR.push_back((*hit).first);
0296     }
0297   }
0298
0299   if (debug) {
0300     cout << "points " << points.size() << endl;
0301     cout << "pointsNoLR " << pointsNoLR.size() << endl;
0302   }
0303
0304   // build all possible candidates using L/R ambiguity
0305   vector<DTSegmentCand*> candidates;
0306
0307   buildPointsCollection(points, pointsNoLR, candidates, sl);
0308
0309   // here I try to add the external clusters and build a set of "extended
0310   // segment candidate
0311   vector<DTSegmentExtendedCand*> extendedCands = extendCandidates(candidates, sl);
0312   if (debug)
0313     cout << "extended candidates " << extendedCands.size() << endl;
0314
0315   // so now I have build a given number of segments, I should find the best one,
0316   // by #hits and chi2.
0317   vector<DTSegmentExtendedCand*>::const_iterator bestCandIter = extendedCands.end();
0318   double minChi2 = 999999.;
0319   unsigned int maxNumHits = 0;
0320   for (vector<DTSegmentExtendedCand*>::const_iterator iter = extendedCands.begin(); iter != extendedCands.end();
0321        ++iter) {
0322     if ((*iter)->nHits() == maxNumHits && (*iter)->chi2() < minChi2) {
0323       minChi2 = (*iter)->chi2();
0324       bestCandIter = iter;
0325     } else if ((*iter)->nHits() > maxNumHits) {
0326       maxNumHits = (*iter)->nHits();
0327       minChi2 = (*iter)->chi2();
0328       bestCandIter = iter;
0329     }
0330   }
0331
0332   // delete all candidates but the best one!
0333   for (vector<DTSegmentExtendedCand*>::iterator iter = extendedCands.begin(); iter != extendedCands.end(); ++iter)
0334     if (iter != bestCandIter)
0335       delete (*iter);
0336
0337   // return the best candate if any
0338   if (bestCandIter != extendedCands.end()) {
0339     return (*bestCandIter);
0340   }
0341   return nullptr;
0342 }
0343
0344 void DTCombinatorialExtendedPatternReco::buildPointsCollection(vector<DTSegmentCand::AssPoint>& points,
0345                                                                deque<std::shared_ptr<DTHitPairForFit>>& pointsNoLR,
0346                                                                vector<DTSegmentCand*>& candidates,
0347                                                                const DTSuperLayer* sl) {
0348   if (debug) {
0349     cout << "DTCombinatorialExtendedPatternReco::buildPointsCollection " << endl;
0350     cout << "points: " << points.size() << " NOLR: " << pointsNoLR.size() << endl;
0351   }
0352   if (!pointsNoLR.empty()) {  // still unassociated points!
0353     std::shared_ptr<DTHitPairForFit> unassHit = pointsNoLR.front();
0354     // try with the right
0355     if (debug)
0356       cout << "Right hit" << endl;
0357     points.push_back(DTSegmentCand::AssPoint(unassHit, DTEnums::Right));
0358     pointsNoLR.pop_front();
0359     buildPointsCollection(points, pointsNoLR, candidates, sl);
0360     pointsNoLR.push_front((unassHit));
0361     points.pop_back();
0362
0363     // try with the left
0364     if (debug)
0365       cout << "Left hit" << endl;
0366     points.push_back(DTSegmentCand::AssPoint(unassHit, DTEnums::Left));
0367     pointsNoLR.pop_front();
0368     buildPointsCollection(points, pointsNoLR, candidates, sl);
0369     pointsNoLR.push_front((unassHit));
0370     points.pop_back();
0371   } else {  // all associated
0372
0373     if (debug) {
0374       cout << "The Hits were" << endl;
0375       copy(points.begin(), points.end(), ostream_iterator<DTSegmentCand::AssPoint>(std::cout));
0376       cout << "----" << endl;
0377       cout << "All associated " << endl;
0378     }
0379     DTSegmentCand::AssPointCont pointsSet;
0380
0381     // for (vector<DTSegmentCand::AssPoint>::const_iterator point=points.begin();
0382     //      point!=points.end(); ++point)
0383     pointsSet.insert(points.begin(), points.end());
0384
0385     if (debug) {
0386       cout << "The Hits are" << endl;
0387       copy(pointsSet.begin(), pointsSet.end(), ostream_iterator<DTSegmentCand::AssPoint>(std::cout));
0388       cout << "----" << endl;
0389     }
0390
0391     DTSegmentCand* newCand = new DTSegmentCand(pointsSet, sl);
0392     if (theUpdator->fit(newCand, false, false))
0393       candidates.push_back(newCand);
0394     else
0395       delete newCand;  // bad seg, too few hits
0396   }
0397 }
0398
0399 bool DTCombinatorialExtendedPatternReco::checkDoubleCandidates(vector<DTSegmentCand*>& cands, DTSegmentCand* seg) {
0400   for (vector<DTSegmentCand*>::iterator cand = cands.begin(); cand != cands.end(); ++cand)
0401     if (*(*cand) == *seg)
0402       return false;
0403   return true;
0404 }
0405
0406 vector<DTSegmentExtendedCand*> DTCombinatorialExtendedPatternReco::extendCandidates(vector<DTSegmentCand*>& candidates,
0407                                                                                     const DTSuperLayer* sl) {
0408   if (debug)
0409     cout << "extendCandidates " << candidates.size() << endl;
0410   vector<DTSegmentExtendedCand*> result;
0411
0412   // in case of phi SL just return
0413   if (sl->id().superLayer() != 2) {
0414     for (vector<DTSegmentCand*>::const_iterator cand = candidates.begin(); cand != candidates.end(); ++cand) {
0415       DTSegmentExtendedCand* extendedCand = new DTSegmentExtendedCand(*cand);
0416       // and delete the original candidate
0417       delete *cand;
0418       result.push_back(extendedCand);
0419     }
0420     return result;
0421   }
0422
0423   // first I have to select the cluster which are compatible with the actual
0424   // candidate, namely +/-1 sector/station/wheel
0425   vector<DTSegmentExtendedCand::DTSLRecClusterForFit> clustersWithPos;
0426   if (debug)
0427     cout << "AllClustersWithPos " << theClusters.size() << endl;
0428   if (debug)
0429     cout << "SL:   " << sl->id() << endl;
0430   for (vector<DTSLRecCluster>::const_iterator clus = theClusters.begin(); clus != theClusters.end(); ++clus) {
0431     if (debug)
0432       cout << "CLUS: " << (*clus).superLayerId() << endl;
0433     if ((*clus).superLayerId().superLayer() == 2 && closeSL(sl->id(), (*clus).superLayerId())) {
0434       // and then get their pos in the actual SL frame
0435       const DTSuperLayer* clusSl = theDTGeometry->superLayer((*clus).superLayerId());
0436       LocalPoint pos = sl->toLocal(clusSl->toGlobal((*clus).localPosition()));
0437       //LocalError err=sl->toLocal(clusSl->toGlobal((*clus).localPositionError()));
0438       LocalError err = (*clus).localPositionError();
0439       clustersWithPos.push_back(DTSegmentExtendedCand::DTSLRecClusterForFit(*clus, pos, err));
0440     }
0441   }
0442   if (debug)
0443     cout << "closeClustersWithPos " << clustersWithPos.size() << endl;
0444
0445   for (vector<DTSegmentCand*>::const_iterator cand = candidates.begin(); cand != candidates.end(); ++cand) {
0446     // create an extended candidate
0447     DTSegmentExtendedCand* extendedCand = new DTSegmentExtendedCand(*cand);
0448     // and delete the original candidate
0449     delete *cand;
0450     // do this only for theta SL
0451     if (extendedCand->superLayer()->id().superLayer() == 2) {
0452       // first check compatibility between cand and clusForFit
0453       for (vector<DTSegmentExtendedCand::DTSLRecClusterForFit>::const_iterator exClus = clustersWithPos.begin();
0454            exClus != clustersWithPos.end();
0455            ++exClus) {
0456         if (extendedCand->isCompatible(*exClus)) {
0457           if (debug)
0458             cout << "is compatible " << endl;
0459           // add compatible cluster
0460           extendedCand->addClus(*exClus);
0461         }
0462       }
0463       // fit the segment
0464       if (debug)
0465         cout << "extended cands nHits: " << extendedCand->nHits() << endl;
0466       if (theUpdator->fit(extendedCand, false, false)) {
0467         // add to result
0468         result.push_back(extendedCand);
0469       } else {
0470         cout << "Bad fit" << endl;
0471         delete extendedCand;
0472       }
0473     } else {  // Phi SuperLayer
0474       result.push_back(extendedCand);
0475     }
0476   }
0477
0478   return result;
0479 }
0480
0481 bool DTCombinatorialExtendedPatternReco::closeSL(const DTSuperLayerId& id1, const DTSuperLayerId& id2) {
0482   if (id1 == id2)
0483     return false;
0484   if (abs(id1.wheel() - id2.wheel()) > 1)
0485     return false;
0486   // take into account also sector 13 and 14
0487   int sec1 = (id1.sector() == 13) ? 4 : id1.sector();
0488   sec1 = (sec1 == 14) ? 10 : sec1;
0489   int sec2 = (id2.sector() == 13) ? 4 : id2.sector();
0490   sec2 = (sec2 == 14) ? 10 : sec2;
0491   // take into account also sector 1/12
0492   if (abs(sec1 - sec2) > 1 && abs(sec1 - sec2) != 11)
0493     return false;
0494   //if (abs(id1.station()-id2.station())>1 ) return false;
0495   return true;
0496 }