DTSegment/src/DTCombinatorialPatternReco.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/DTCombinatorialPatternReco.h"
0009
0010 /* Collaborating Class Header */
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "DataFormats/Common/interface/OwnVector.h"
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0017 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0018 #include "Geometry/DTGeometry/interface/DTLayer.h"
0019 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0020 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0021 #include "RecoLocalMuon/DTSegment/src/DTSegmentUpdator.h"
0022 #include "RecoLocalMuon/DTSegment/src/DTSegmentCleaner.h"
0023 #include "RecoLocalMuon/DTSegment/src/DTHitPairForFit.h"
0024 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.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 DTCombinatorialPatternReco::DTCombinatorialPatternReco(const edm::ParameterSet& pset, edm::ConsumesCollector iCC)
0035     : DTRecSegment2DBaseAlgo(pset), theAlgoName("DTCombinatorialPatternReco"), theDTGeometryToken(iCC.esConsumes()) {
0036   theMaxAllowedHits = pset.getParameter<unsigned int>("MaxAllowedHits");  // 100
0037   theAlphaMaxTheta = pset.getParameter<double>("AlphaMaxTheta");          // 0.1 ;
0038   theAlphaMaxPhi = pset.getParameter<double>("AlphaMaxPhi");              // 1.0 ;
0039   debug = pset.getUntrackedParameter<bool>("debug");                      //true;
0040   theUpdator = new DTSegmentUpdator(pset, iCC);
0041   theCleaner = new DTSegmentCleaner(pset);
0042   string theHitAlgoName = pset.getParameter<string>("recAlgo");
0043   usePairs = !(theHitAlgoName == "DTNoDriftAlgo");
0044 }
0045
0046 /// Destructor
0047 DTCombinatorialPatternReco::~DTCombinatorialPatternReco() {
0048   delete theUpdator;
0049   delete theCleaner;
0050 }
0051
0052 /* Operations */
0053 edm::OwnVector<DTSLRecSegment2D> DTCombinatorialPatternReco::reconstruct(const DTSuperLayer* sl,
0054                                                                          const std::vector<DTRecHit1DPair>& pairs) {
0055   edm::OwnVector<DTSLRecSegment2D> result;
0056   vector<std::shared_ptr<DTHitPairForFit>> hitsForFit = initHits(sl, pairs);
0057
0058   vector<DTSegmentCand*> candidates = buildSegments(sl, hitsForFit);
0059
0060   vector<DTSegmentCand*>::const_iterator cand = candidates.begin();
0061   while (cand < candidates.end()) {
0062     DTSLRecSegment2D* segment = (**cand);
0063
0064     theUpdator->update(segment, false);
0065
0066     result.push_back(segment);
0067
0068     if (debug) {
0069       cout << "Reconstructed 2D segments " << result.back() << endl;
0070     }
0071
0072     delete *(cand++);  // delete the candidate!
0073   }
0074
0075   return result;
0076 }
0077
0078 void DTCombinatorialPatternReco::setES(const edm::EventSetup& setup) {
0079   // Get the DT Geometry
0080   theDTGeometry = setup.getHandle(theDTGeometryToken);
0081   theUpdator->setES(setup);
0082 }
0083
0084 vector<std::shared_ptr<DTHitPairForFit>> DTCombinatorialPatternReco::initHits(const DTSuperLayer* sl,
0085                                                                               const std::vector<DTRecHit1DPair>& hits) {
0086   vector<std::shared_ptr<DTHitPairForFit>> result;
0087   for (vector<DTRecHit1DPair>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0088     result.push_back(std::make_shared<DTHitPairForFit>(*hit, *sl, theDTGeometry));
0089   }
0090   return result;
0091 }
0092
0093 vector<DTSegmentCand*> DTCombinatorialPatternReco::buildSegments(
0094     const DTSuperLayer* sl, const std::vector<std::shared_ptr<DTHitPairForFit>>& hits) {
0095   // clear the patterns tried
0096   if (debug) {
0097     cout << "theTriedPattern.size is " << theTriedPattern.size() << "\n";
0098   }
0099   theTriedPattern.clear();
0100
0101   typedef vector<std::shared_ptr<DTHitPairForFit>> hitCont;
0102   typedef hitCont::const_iterator hitIter;
0103   vector<DTSegmentCand*> result;
0104
0105   if (debug) {
0106     cout << "buildSegments: " << sl->id() << " nHits " << hits.size() << endl;
0107     for (vector<std::shared_ptr<DTHitPairForFit>>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit)
0108       cout << **hit << endl;
0109   }
0110
0111   // 10-Mar-2004 SL
0112   // put a protection against heavily populated chambers, for which the segment
0113   // building could lead to infinite memory usage...
0114   if (hits.size() > theMaxAllowedHits) {
0115     if (debug) {
0116       cout << "Warning: this SuperLayer " << sl->id() << " has too many hits : " << hits.size() << " max allowed is "
0117            << theMaxAllowedHits << endl;
0118       cout << "Skipping segment reconstruction... " << endl;
0119     }
0120     return result;
0121   }
0122
0123   /// get two hits in different layers and see if there are other / hits
0124   //  compatible with them
0125   for (hitCont::const_iterator firstHit = hits.begin(); firstHit != hits.end(); ++firstHit) {
0126     for (hitCont::const_reverse_iterator lastHit = hits.rbegin(); (*lastHit) != (*firstHit); ++lastHit) {
0127       //if ( (*lastHit)->id().layerId() == (*firstHit)->id().layerId() ) continue; // hits must be in different layers!
0128       // hits must nor in the same nor in adiacent layers
0129       if (fabs((*lastHit)->id().layerId() - (*firstHit)->id().layerId()) <= 1)
0130         continue;
0131       if (debug) {
0132         cout << "Selected these two hits pair " << endl;
0133         cout << "First " << *(*firstHit) << " Layer Id: " << (*firstHit)->id().layerId() << endl;
0134         cout << "Last " << *(*lastHit) << " Layer Id: " << (*lastHit)->id().layerId() << endl;
0135       }
0136
0137       GlobalPoint IP;
0138       float DAlphaMax;
0139       if ((sl->id()).superlayer() == 2)  // Theta SL
0140         DAlphaMax = theAlphaMaxTheta;
0141       else  // Phi SL
0142         DAlphaMax = theAlphaMaxPhi;
0143
0144       DTEnums::DTCellSide codes[2] = {DTEnums::Right, DTEnums::Left};
0145       for (int firstLR = 0; firstLR < 2; ++firstLR) {
0146         for (int lastLR = 0; lastLR < 2; ++lastLR) {
0147           // TODO move the global transformation in the DTHitPairForFit class
0148           // when it will be moved I will able to remove the sl from the input parameter
0149           GlobalPoint gposFirst = sl->toGlobal((*firstHit)->localPosition(codes[firstLR]));
0150           GlobalPoint gposLast = sl->toGlobal((*lastHit)->localPosition(codes[lastLR]));
0151
0152           GlobalVector gvec = gposLast - gposFirst;
0153           GlobalVector gvecIP = gposLast - IP;
0154
0155           // difference in angle measured
0156           float DAlpha = fabs(gvec.theta() - gvecIP.theta());
0157
0158           // cout << "DAlpha " << DAlpha << endl;
0159           if (DAlpha < DAlphaMax) {
0160             // create a segment hypotesis
0161             // I don't need a true segment, just direction and position
0162             LocalPoint posIni = (*firstHit)->localPosition(codes[firstLR]);
0163             LocalVector dirIni = ((*lastHit)->localPosition(codes[lastLR]) - posIni).unit();
0164
0165             // search for other compatible hits, with or without the L/R solved
0166             vector<DTSegmentCand::AssPoint> assHits = findCompatibleHits(posIni, dirIni, hits);
0167             if (debug)
0168               cout << "compatible hits " << assHits.size() << endl;
0169
0170             // get the best segment with these hits: it's just one!
0171             // (is it correct?)
0172             DTSegmentCand* seg = buildBestSegment(assHits, sl);
0173
0174             if (seg) {
0175               if (debug)
0176                 cout << "segment " << *seg << endl;
0177
0178               // check if the chi2 and #hits are ok
0179               if (!seg->good()) {
0180                 delete seg;
0181               } else {
0182                 // remove duplicated segments (I know, would be better to do it before the
0183                 // fit...)
0184                 if (checkDoubleCandidates(result, seg)) {
0185                   // add to the vector of hypotesis
0186                   result.push_back(seg);
0187                   if (debug)
0188                     cout << "result is now " << result.size() << endl;
0189                 } else {  // delete it!
0190                   delete seg;
0191                   if (debug)
0192                     cout << "already existing" << endl;
0193                 }
0194               }
0195             }
0196           }
0197         }
0198       }
0199     }
0200   }
0201   if (debug) {
0202     for (vector<DTSegmentCand*>::const_iterator seg = result.begin(); seg != result.end(); ++seg)
0203       cout << *(*seg) << endl;
0204   }
0205
0206   // now I have a couple of segment hypotesis, should check for ghost
0207   result = theCleaner->clean(result);
0208   if (debug) {
0209     cout << "result no ghost  " << result.size() << endl;
0210     for (vector<DTSegmentCand*>::const_iterator seg = result.begin(); seg != result.end(); ++seg)
0211       cout << *(*seg) << endl;
0212   }
0213
0214   return result;
0215 }
0216
0217 vector<DTSegmentCand::AssPoint> DTCombinatorialPatternReco::findCompatibleHits(
0218     const LocalPoint& posIni, const LocalVector& dirIni, const vector<std::shared_ptr<DTHitPairForFit>>& hits) {
0219   if (debug)
0220     cout << "Pos: " << posIni << " Dir: " << dirIni << endl;
0221   vector<DTSegmentCand::AssPoint> result;
0222
0223   // counter to early-avoid double counting in hits pattern
0224   TriedPattern tried;
0225   int nCompatibleHits = 0;
0226
0227   typedef vector<std::shared_ptr<DTHitPairForFit>> hitCont;
0228   typedef hitCont::const_iterator hitIter;
0229   for (hitIter hit = hits.begin(); hit != hits.end(); ++hit) {
0230     pair<bool, bool> isCompatible = (*hit)->isCompatible(posIni, dirIni);
0231     if (debug)
0232       cout << "isCompatible " << isCompatible.first << " " << isCompatible.second << endl;
0233
0234     // if only one of the two is compatible, then the LR is assigned,
0235     // otherwise is undefined
0236
0237     DTEnums::DTCellSide lrcode;
0238     if (isCompatible.first && isCompatible.second) {
0239       usePairs ? lrcode = DTEnums::undefLR : lrcode = DTEnums::Left;  // if not usePairs then only use single side
0240       tried.push_back(3);
0241       nCompatibleHits++;
0242     } else if (isCompatible.first) {
0243       lrcode = DTEnums::Left;
0244       tried.push_back(2);
0245       nCompatibleHits++;
0246     } else if (isCompatible.second) {
0247       lrcode = DTEnums::Right;
0248       tried.push_back(1);
0249       nCompatibleHits++;
0250     } else {
0251       tried.push_back(0);
0252       continue;  // neither is compatible
0253     }
0254     result.push_back(DTSegmentCand::AssPoint(*hit, lrcode));
0255   }
0256
0257   // check if too few associated hits or pattern already tried
0258   // TODO: two different if for nCompatibleHits < 3 =>printout and find ->
0259   // printour
0260   if (nCompatibleHits < 3) {
0261     if (debug) {
0262       cout << "Less than 3 hits ";
0263       copy(tried.begin(), tried.end(), ostream_iterator<int>(std::cout));
0264       cout << endl;
0265     }
0266     // No need to insert into the list of patterns, as you will never search for it.
0267     //theTriedPattern.insert(tried);
0268   } else {
0269     // try to insert, return bool if already there
0270     bool isnew = theTriedPattern.insert(tried).second;
0271     if (isnew) {
0272       if (debug) {
0273         cout << "NOT Already tried ";
0274         copy(tried.begin(), tried.end(), ostream_iterator<int>(std::cout));
0275         cout << endl;
0276       }
0277     } else {
0278       if (debug) {
0279         cout << "Already tried ";
0280         copy(tried.begin(), tried.end(), ostream_iterator<int>(std::cout));
0281         cout << endl;
0282       }
0283       // empty the result vector
0284       result.clear();
0285     }
0286   }
0287   return result;
0288 }
0289
0290 DTSegmentCand* DTCombinatorialPatternReco::buildBestSegment(std::vector<DTSegmentCand::AssPoint>& hits,
0291                                                             const DTSuperLayer* sl) {
0292   if (hits.size() < 3) {
0293     //cout << "buildBestSegment: hits " << hits.size()<< endl;
0294     return nullptr;  // a least 3 point
0295   }
0296
0297   // hits with defined LR
0298   vector<DTSegmentCand::AssPoint> points;
0299
0300   // without: I store both L and R, a deque since I need front insertion and
0301   // deletion
0302   deque<std::shared_ptr<DTHitPairForFit>> pointsNoLR;
0303
0304   // first add only the hits with LR assigned
0305   for (vector<DTSegmentCand::AssPoint>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0306     if ((*hit).second != DTEnums::undefLR) {
0307       points.push_back(*hit);
0308     } else {  // then also for the undef'd one
0309       pointsNoLR.push_back((*hit).first);
0310     }
0311   }
0312
0313   if (debug) {
0314     cout << "points " << points.size() << endl;
0315     cout << "pointsNoLR " << pointsNoLR.size() << endl;
0316   }
0317
0318   // build all possible candidates using L/R ambiguity
0319   vector<DTSegmentCand*> candidates;
0320
0321   buildPointsCollection(points, pointsNoLR, candidates, sl);
0322
0323   if (debug)
0324     cout << "candidates " << candidates.size() << endl;
0325
0326   // so now I have build a given number of segments, I should find the best one,
0327   // by #hits and chi2.
0328   vector<DTSegmentCand*>::const_iterator bestCandIter = candidates.end();
0329   double minChi2 = 999999.;
0330   unsigned int maxNumHits = 0;
0331   for (vector<DTSegmentCand*>::const_iterator iter = candidates.begin(); iter != candidates.end(); ++iter) {
0332     if ((*iter)->nHits() == maxNumHits && (*iter)->chi2() < minChi2) {
0333       minChi2 = (*iter)->chi2();
0334       bestCandIter = iter;
0335     } else if ((*iter)->nHits() > maxNumHits) {
0336       maxNumHits = (*iter)->nHits();
0337       minChi2 = (*iter)->chi2();
0338       bestCandIter = iter;
0339     }
0340   }
0341
0342   // delete all candidates but the best one!
0343   for (vector<DTSegmentCand*>::iterator iter = candidates.begin(); iter != candidates.end(); ++iter)
0344     if (iter != bestCandIter)
0345       delete *iter;
0346
0347   // return the best candate if any
0348   if (bestCandIter != candidates.end()) {
0349     return (*bestCandIter);
0350   }
0351   return nullptr;
0352 }
0353
0354 void DTCombinatorialPatternReco::buildPointsCollection(vector<DTSegmentCand::AssPoint>& points,
0355                                                        deque<std::shared_ptr<DTHitPairForFit>>& pointsNoLR,
0356                                                        vector<DTSegmentCand*>& candidates,
0357                                                        const DTSuperLayer* sl) {
0358   if (debug) {
0359     cout << "buildPointsCollection " << endl;
0360     cout << "points: " << points.size() << " NOLR: " << pointsNoLR.size() << endl;
0361   }
0362   if (!pointsNoLR.empty()) {  // still unassociated points!
0363     std::shared_ptr<DTHitPairForFit> unassHit = pointsNoLR.front();
0364     // try with the right
0365     if (debug)
0366       cout << "Right hit" << endl;
0367     points.push_back(DTSegmentCand::AssPoint(unassHit, DTEnums::Right));
0368     pointsNoLR.pop_front();
0369     buildPointsCollection(points, pointsNoLR, candidates, sl);
0370     pointsNoLR.push_front((unassHit));
0371     points.pop_back();
0372
0373     // try with the left
0374     if (debug)
0375       cout << "Left hit" << endl;
0376     points.push_back(DTSegmentCand::AssPoint(unassHit, DTEnums::Left));
0377     pointsNoLR.pop_front();
0378     buildPointsCollection(points, pointsNoLR, candidates, sl);
0379     pointsNoLR.push_front((unassHit));
0380     points.pop_back();
0381   } else {  // all associated
0382
0383     if (debug) {
0384       cout << "The Hits were" << endl;
0385       copy(points.begin(), points.end(), ostream_iterator<DTSegmentCand::AssPoint>(std::cout));
0386       cout << "----" << endl;
0387       cout << "All associated " << endl;
0388     }
0389     DTSegmentCand::AssPointCont pointsSet;
0390
0391     // for (vector<DTSegmentCand::AssPoint>::const_iterator point=points.begin();
0392     //      point!=points.end(); ++point)
0393     pointsSet.insert(points.begin(), points.end());
0394
0395     if (debug) {
0396       cout << "The Hits are" << endl;
0397       copy(pointsSet.begin(), pointsSet.end(), ostream_iterator<DTSegmentCand::AssPoint>(std::cout));
0398       cout << "----" << endl;
0399     }
0400
0401     DTSegmentCand* newCand = new DTSegmentCand(pointsSet, sl);
0402     if (theUpdator->fit(newCand, false, false))
0403       candidates.push_back(newCand);
0404     else
0405       delete newCand;  // bad seg, too few hits
0406   }
0407 }
0408
0409 bool DTCombinatorialPatternReco::checkDoubleCandidates(vector<DTSegmentCand*>& cands, DTSegmentCand* seg) {
0410   for (vector<DTSegmentCand*>::iterator cand = cands.begin(); cand != cands.end(); ++cand)
0411     if (*(*cand) == *seg)
0412       return false;
0413   return true;
0414 }