DTSegment/src/DTSegmentCleaner.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/DTSegmentCleaner.h"
0009
0010 /* Collaborating Class Header */
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 /* C++ Headers */
0015 using namespace std;
0016
0017 /* ====================================================================== */
0018
0019 /* Constructor */
0020 DTSegmentCleaner::DTSegmentCleaner(const edm::ParameterSet& pset) {
0021   nSharedHitsMax = pset.getParameter<int>("nSharedHitsMax");
0022
0023   nUnSharedHitsMin = pset.getParameter<int>("nUnSharedHitsMin");
0024
0025   segmCleanerMode = pset.getParameter<int>("segmCleanerMode");
0026
0027   if ((segmCleanerMode != 1) && (segmCleanerMode != 2) && (segmCleanerMode != 3))
0028     edm::LogError("Muon|RecoLocalMuon|DTSegmentCleaner")
0029         << "Wrong segmCleanerMode! It must be 1,2 or 3. The default is 1";
0030 }
0031
0032 /* Destructor */
0033 DTSegmentCleaner::~DTSegmentCleaner() {}
0034
0035 /* Operations */
0036 vector<DTSegmentCand*> DTSegmentCleaner::clean(const std::vector<DTSegmentCand*>& inputCands) const {
0037   if (inputCands.size() < 2)
0038     return inputCands;
0039   //cout << "[DTSegmentCleaner] # of candidates: " << inputCands.size() << endl;
0040   vector<DTSegmentCand*> result = solveConflict(inputCands);
0041
0042   //cout << "[DTSegmentCleaner] to ghostbuster: " << result.size() << endl;
0043   result = ghostBuster(result);
0044
0045   return result;
0046 }
0047
0048 vector<DTSegmentCand*> DTSegmentCleaner::solveConflict(const std::vector<DTSegmentCand*>& inputCands) const {
0049   vector<DTSegmentCand*> result;
0050
0051   vector<DTSegmentCand*> ghosts;
0052
0053   for (vector<DTSegmentCand*>::const_iterator cand = inputCands.begin(); cand != inputCands.end(); ++cand) {
0054     for (vector<DTSegmentCand*>::const_iterator cand2 = cand + 1; cand2 != inputCands.end(); ++cand2) {
0055       DTSegmentCand::AssPointCont confHits = (*cand)->conflictingHitPairs(*(*cand2));
0056
0057       if (!confHits.empty()) {
0058         ///treatment of LR ambiguity cases: 1 chooses the best chi2
0059         ///                                 2 chooses the smaller angle
0060         ///                                 3 keeps both candidates
0061         if ((confHits.size()) == ((*cand)->nHits()) && (confHits.size()) == ((*cand2)->nHits()) &&
0062             (fabs((*cand)->chi2() - (*cand2)->chi2()) < 0.1)) {  // cannot choose on the basis of # of hits or chi2
0063
0064           if (segmCleanerMode == 2) {  // mode 2: choose on the basis of the angle
0065
0066             DTSegmentCand* badCand = nullptr;
0067             if ((*cand)->superLayer()->id().superlayer() != 2) {  // we are in the phi view
0068
0069               LocalVector dir1 = (*cand)->direction();
0070               LocalVector dir2 = (*cand2)->direction();
0071               float phi1 = (atan((dir1.x()) / (dir1.z())));
0072               float phi2 = (atan((dir2.x()) / (dir2.z())));
0073
0074               badCand = (fabs(phi1) > fabs(phi2)) ? (*cand) : (*cand2);
0075
0076             } else {  // we are in the theta view: choose the most pointing one
0077
0078               GlobalPoint IP;
0079
0080               GlobalVector cand1GlobDir = (*cand)->superLayer()->toGlobal((*cand)->direction());
0081               GlobalPoint cand1GlobPos = (*cand)->superLayer()->toGlobal((*cand)->position());
0082               GlobalVector cand1GlobVecIP = cand1GlobPos - IP;
0083               float DAlpha1 = fabs(cand1GlobDir.theta() - cand1GlobVecIP.theta());
0084
0085               GlobalVector cand2GlobDir = (*cand2)->superLayer()->toGlobal((*cand2)->direction());
0086               GlobalPoint cand2GlobPos = (*cand2)->superLayer()->toGlobal((*cand2)->position());
0087               GlobalVector cand2GlobVecIP = cand2GlobPos - IP;
0088               float DAlpha2 = fabs(cand2GlobDir.theta() - cand2GlobVecIP.theta());
0089
0090               badCand = (DAlpha1 > DAlpha2) ? (*cand) : (*cand2);
0091             }
0092
0093             for (DTSegmentCand::AssPointCont::const_iterator cHit = confHits.begin(); cHit != confHits.end(); ++cHit) {
0094               badCand->removeHit(*cHit);
0095             }
0096
0097           } else {  // mode 3: keep both candidates
0098             continue;
0099           }
0100
0101         } else {  // mode 1: take > # hits or best chi2
0102           //cout << " Choose based on hits/chi2. " << endl << (**cand) << " vs " << (**cand2) << endl;
0103
0104           // if one segment is in-time and the other not then keep both
0105           if (((*cand)->t0() * (*cand2)->t0() != 0) || ((*cand)->t0() == (*cand2)->t0())) {
0106             DTSegmentCand* badCand = (**cand) < (**cand2) ? (*cand) : (*cand2);
0107             //cout << " Remove hits in " << (*badCand) << endl;
0108             for (DTSegmentCand::AssPointCont::const_iterator cHit = confHits.begin(); cHit != confHits.end(); ++cHit)
0109               badCand->removeHit(*cHit);
0110           }
0111         }
0112       }
0113     }
0114   }
0115
0116   vector<DTSegmentCand*>::const_iterator cand = inputCands.begin();
0117   while (cand < inputCands.end()) {
0118     if ((*cand)->good())
0119       result.push_back(*cand);
0120     else {
0121       vector<DTSegmentCand*>::const_iterator badCand = cand;
0122       delete *badCand;
0123     }
0124     ++cand;
0125   }
0126   return result;
0127 }
0128
0129 vector<DTSegmentCand*> DTSegmentCleaner::ghostBuster(const std::vector<DTSegmentCand*>& inputCands) const {
0130   vector<DTSegmentCand*> ghosts;
0131   for (vector<DTSegmentCand*>::const_iterator cand = inputCands.begin(); cand != inputCands.end(); ++cand) {
0132     for (vector<DTSegmentCand*>::const_iterator cand2 = cand + 1; cand2 != inputCands.end(); ++cand2) {
0133       unsigned int nSharedHits = (*cand)->nSharedHitPairs(*(*cand2));
0134       ///cout << "Sharing " << (**cand) << " " << (**cand2) << " " << nSharedHits
0135       //       << " (first or second) " << ((**cand)<(**cand2)) << endl;
0136       if ((nSharedHits == ((*cand)->nHits())) && (nSharedHits == ((*cand2)->nHits())) &&
0137           (fabs((*cand)->chi2() - (*cand2)->chi2()) < 0.1) && (segmCleanerMode == 3)) {
0138         continue;
0139       }
0140
0141       if (((*cand)->nHits() == 3 || (*cand2)->nHits() == 3) && (fabs((*cand)->chi2() - (*cand2)->chi2()) < 0.0001)) {
0142         continue;
0143       }
0144
0145       // if one segment is in-time and the other not then keep both
0146       if (((*cand2)->nHits() == (*cand)->nHits()) && ((*cand)->t0() * (*cand2)->t0() == 0) &&
0147           ((*cand)->t0() != (*cand2)->t0())) {
0148         continue;
0149       }
0150
0151       // remove the worst segment if too many shared hits or too few unshared
0152       if ((int)nSharedHits >= nSharedHitsMax || (int)((*cand)->nHits() - nSharedHits) <= nUnSharedHitsMin ||
0153           (int)((*cand2)->nHits() - nSharedHits) <= nUnSharedHitsMin) {
0154         if ((**cand) < (**cand2)) {
0155           ghosts.push_back(*cand);
0156         } else {
0157           ghosts.push_back(*cand2);
0158         }
0159         continue;
0160       }
0161     }
0162   }
0163
0164   vector<DTSegmentCand*> result;
0165   for (vector<DTSegmentCand*>::const_iterator cand = inputCands.begin(); cand != inputCands.end(); ++cand) {
0166     bool isGhost = false;
0167     for (vector<DTSegmentCand*>::const_iterator ghost = ghosts.begin(); ghost != ghosts.end(); ++ghost) {
0168       if ((*cand) == (*ghost)) {
0169         isGhost = true;
0170         break;
0171       }
0172     }
0173     if (!isGhost)
0174       result.push_back(*cand);
0175     else
0176       delete *cand;
0177   }
0178   // cout << "No Ghosts ------" << endl;
0179   // for (vector<DTSegmentCand*>::iterator cand=result.begin();
0180   //      cand!=result.end(); ++cand) {
0181   //   cout << "cand " << *cand << " nH " <<(*cand)->nHits() << " chi2 " << (*cand)->chi2() << endl;
0182   // }
0183   // cout << "----------------" << endl;
0184
0185   return result;
0186 }