MuonSeedGenerator/src/RPCSeedOverlapper.cc

0001 /**

0002  *  See header file for a description of this class.

0003  *

0004  */
0005
0006 #include "RecoMuon/MuonSeedGenerator/src/RPCSeedOverlapper.h"
0007 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0008 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0009
0010 using namespace std;
0011 using namespace edm;
0012
0013 RPCSeedOverlapper::RPCSeedOverlapper() : rpcGeometry{nullptr} {
0014   isConfigured = false;
0015   isIOset = false;
0016 }
0017
0018 RPCSeedOverlapper::~RPCSeedOverlapper() {}
0019
0020 void RPCSeedOverlapper::configure(const edm::ParameterSet &iConfig) {
0021   isCheckgoodOverlap = iConfig.getParameter<bool>("isCheckgoodOverlap");
0022   isCheckcandidateOverlap = iConfig.getParameter<bool>("isCheckcandidateOverlap");
0023   ShareRecHitsNumberThreshold = iConfig.getParameter<unsigned int>("ShareRecHitsNumberThreshold");
0024   isConfigured = true;
0025 }
0026
0027 void RPCSeedOverlapper::setIO(std::vector<weightedTrajectorySeed> *goodweightedRef,
0028                               std::vector<weightedTrajectorySeed> *candidateweightedRef) {
0029   goodweightedSeedsRef = goodweightedRef;
0030   candidateweightedSeedsRef = candidateweightedRef;
0031   isIOset = true;
0032 }
0033
0034 void RPCSeedOverlapper::unsetIO() { isIOset = false; }
0035
0036 void RPCSeedOverlapper::setGeometry(const RPCGeometry &iGeom) { rpcGeometry = &iGeom; }
0037
0038 void RPCSeedOverlapper::run() {
0039   if (isConfigured == false || isIOset == false || rpcGeometry == nullptr) {
0040     cout << "Configuration or IO is not set yet" << endl;
0041     return;
0042   }
0043   if (isCheckgoodOverlap == true)
0044     CheckOverlap(*rpcGeometry, goodweightedSeedsRef);
0045   if (isCheckcandidateOverlap == true)
0046     CheckOverlap(*rpcGeometry, candidateweightedSeedsRef);
0047 }
0048
0049 void RPCSeedOverlapper::CheckOverlap(const RPCGeometry &rpcGeometry,
0050                                      std::vector<weightedTrajectorySeed> *weightedSeedsRef) {
0051   std::vector<weightedTrajectorySeed> sortweightedSeeds;
0052   std::vector<weightedTrajectorySeed> tempweightedSeeds;
0053   std::vector<TrackingRecHit const *> tempRecHits;
0054
0055   while (!weightedSeedsRef->empty()) {
0056     cout << "Finding the weighted seeds group from " << weightedSeedsRef->size() << " seeds which share some recHits"
0057          << endl;
0058     // Take 1st seed in SeedsRef as referrence and find a collection which always share some recHits with some other

0059     tempRecHits.clear();
0060     tempweightedSeeds.clear();
0061     int N = 0;
0062     for (vector<weightedTrajectorySeed>::iterator itweightedseed = weightedSeedsRef->begin();
0063          itweightedseed != weightedSeedsRef->end();
0064          N++) {
0065       auto const &recHitsRange = itweightedseed->first.recHits();
0066       if (N == 0) {
0067         cout << "Always take the 1st weighted seed to be the referrence." << endl;
0068         for (auto const &hit : recHitsRange) {
0069           cout << "Put its recHits to tempRecHits" << endl;
0070           tempRecHits.push_back(&hit);
0071         }
0072         cout << "Put it to tempweightedSeeds" << endl;
0073         tempweightedSeeds.push_back(*itweightedseed);
0074         cout << "Then erase from weightedSeedsRef->" << endl;
0075         itweightedseed = weightedSeedsRef->erase(itweightedseed);
0076       } else {
0077         cout << "Come to other weighted seed for checking " << itweightedseed->first.nHits() << " recHits from "
0078              << tempRecHits.size() << " temp recHits" << endl;
0079         unsigned int ShareRecHitsNumber = 0;
0080         for (auto const &hit : recHitsRange) {
0081           if (isShareHit(tempRecHits, hit, rpcGeometry))
0082             ShareRecHitsNumber++;
0083         }
0084         if (ShareRecHitsNumber >= ShareRecHitsNumberThreshold) {
0085           cout << "This seed is found to belong to current share group" << endl;
0086           for (auto const &hit : recHitsRange) {
0087             if (!isShareHit(tempRecHits, hit, rpcGeometry)) {
0088               cout << "Put its extra recHits to tempRecHits" << endl;
0089               tempRecHits.push_back(&hit);
0090             }
0091           }
0092           cout << "Put it to tempSeeds" << endl;
0093           tempweightedSeeds.push_back(*itweightedseed);
0094           cout << "Then erase from SeedsRef" << endl;
0095           itweightedseed = weightedSeedsRef->erase(itweightedseed);
0096         } else
0097           itweightedseed++;
0098       }
0099     }
0100     // Find the best weighted seed and kick out those share recHits with it

0101     // The best weighted seed save in sortweightedSeeds, those don't share recHits with it will be push back to weightedSeedsRef for next while loop

0102     weightedTrajectorySeed bestweightedSeed;
0103     vector<weightedTrajectorySeed>::iterator bestweightediter;
0104     // Find the min Spt wrt Pt as the best Seed

0105     double Quality = 1000000;
0106     unsigned NumberofHits = 0;
0107     cout << "Find " << tempweightedSeeds.size() << " seeds into one trajectory group" << endl;
0108     for (vector<weightedTrajectorySeed>::iterator itweightedseed = tempweightedSeeds.begin();
0109          itweightedseed != tempweightedSeeds.end();
0110          itweightedseed++) {
0111       unsigned int nHits = itweightedseed->first.nHits();
0112       //std::vector<float> seed_error = itweightedseed->first.startingState().errorMatrix();

0113       //double Spt = seed_error[1];

0114       double weightedQuality = itweightedseed->second;
0115       cout << "Find a weighted seed with quality " << weightedQuality << endl;
0116       if ((NumberofHits < nHits) || (NumberofHits == nHits && weightedQuality < Quality)) {
0117         NumberofHits = nHits;
0118         Quality = weightedQuality;
0119         bestweightedSeed = *itweightedseed;
0120         bestweightediter = itweightedseed;
0121       }
0122     }
0123     cout << "Best good temp seed's quality is " << Quality << endl;
0124     sortweightedSeeds.push_back(bestweightedSeed);
0125     tempweightedSeeds.erase(bestweightediter);
0126     tempRecHits.clear();
0127
0128     for (auto const &hit : bestweightedSeed.first.recHits()) {
0129       tempRecHits.push_back(&hit);
0130     }
0131
0132     for (vector<weightedTrajectorySeed>::iterator itweightedseed = tempweightedSeeds.begin();
0133          itweightedseed != tempweightedSeeds.end();) {
0134       cout << "Checking the temp weighted seed's " << itweightedseed->first.nHits() << " hits to " << tempRecHits.size()
0135            << " temp recHits" << endl;
0136       bool isShare = false;
0137       for (auto const &hit : itweightedseed->first.recHits()) {
0138         if (isShareHit(tempRecHits, hit, rpcGeometry))
0139           isShare = true;
0140       }
0141
0142       if (isShare == true) {
0143         cout << "Find one temp seed share some recHits with best weighted seed" << endl;
0144         itweightedseed = tempweightedSeeds.erase(itweightedseed);
0145       } else {
0146         cout << "This seed has no relation with best weighted seed" << endl;
0147         weightedSeedsRef->push_back(*itweightedseed);
0148         itweightedseed = tempweightedSeeds.erase(itweightedseed);
0149       }
0150     }
0151   }
0152   // At the end exchange SeedsRef with sortSeeds

0153   weightedSeedsRef->clear();
0154   *weightedSeedsRef = sortweightedSeeds;
0155 }
0156
0157 bool RPCSeedOverlapper::isShareHit(const std::vector<TrackingRecHit const *> &recHits,
0158                                    const TrackingRecHit &hit,
0159                                    const RPCGeometry &rpcGeometry) {
0160   bool istheSame = false;
0161   unsigned int n = 1;
0162   cout << "Checking from " << recHits.size() << " temp recHits" << endl;
0163
0164   LocalPoint lpos1 = hit.localPosition();
0165   DetId RPCId1 = hit.geographicalId();
0166   const GeomDetUnit *rpcroll1 = rpcGeometry.idToDetUnit(RPCId1);
0167   GlobalPoint gpos1 = rpcroll1->toGlobal(lpos1);
0168   cout << "The hit's position: " << gpos1.x() << ", " << gpos1.y() << ", " << gpos1.z() << endl;
0169   for (auto const &recHit : recHits) {
0170     cout << "Checking the " << (n++) << " th recHit from tempRecHits" << endl;
0171     LocalPoint lpos2 = recHit->localPosition();
0172     DetId RPCId2 = recHit->geographicalId();
0173     const GeomDetUnit *rpcroll2 = rpcGeometry.idToDetUnit(RPCId2);
0174     GlobalPoint gpos2 = rpcroll2->toGlobal(lpos2);
0175     cout << "The temp hit's position: " << gpos2.x() << ", " << gpos2.y() << ", " << gpos2.z() << endl;
0176
0177     if ((gpos1.x() == gpos2.x()) && (gpos1.y() == gpos2.y()) && (gpos1.z() == gpos2.z())) {
0178       cout << "This hit is found to be the same" << endl;
0179       istheSame = true;
0180     }
0181   }
0182   return istheSame;
0183 }