Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:19:47 

0001 #include "L1Trigger/DTTriggerPhase2/interface/HoughGrouping.h"
0002 
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 
0007 #include "DataFormats/Math/interface/CMSUnits.h"
0008 
0009 #include <cmath>
0010 #include <memory>
0011 
0012 #include "TMath.h"
0013 
0014 using namespace std;
0015 using namespace edm;
0016 using namespace cmsdt;
0017 using namespace cms_units::operators;
0018 
0019 namespace {
0020   struct {
0021     bool operator()(ProtoCand a, ProtoCand b) const {
0022       unsigned short int sumhqa = 0;
0023       unsigned short int sumhqb = 0;
0024       unsigned short int sumlqa = 0;
0025       unsigned short int sumlqb = 0;
0026       double sumdista = 0;
0027       double sumdistb = 0;
0028 
0029       for (unsigned short int lay = 0; lay < 8; lay++) {
0030         sumhqa += (unsigned short int)a.isThereHitInLayer_[lay];
0031         sumhqb += (unsigned short int)b.isThereHitInLayer_[lay];
0032         sumlqa += (unsigned short int)a.isThereNeighBourHitInLayer_[lay];
0033         sumlqb += (unsigned short int)b.isThereNeighBourHitInLayer_[lay];
0034         sumdista += a.xDistToPattern_[lay];
0035         sumdistb += b.xDistToPattern_[lay];
0036       }
0037 
0038       if (a.nLayersWithHits_ != b.nLayersWithHits_)
0039         return (a.nLayersWithHits_ > b.nLayersWithHits_);  // number of layers with hits
0040       else if (sumhqa != sumhqb)
0041         return (sumhqa > sumhqb);  // number of hq hits
0042       else if (sumlqa != sumlqb)
0043         return (sumlqa > sumlqb);  // number of lq hits
0044       else if (a.nHitsDiff_ != b.nHitsDiff_)
0045         return (a.nHitsDiff_ < b.nHitsDiff_);  // abs. diff. between SL1 & SL3 hits
0046       else
0047         return (sumdista < sumdistb);  // abs. dist. to digis
0048     }
0049   } const HoughOrdering;
0050 }  // namespace
0051 // ============================================================================
0052 // Constructors and destructor
0053 // ============================================================================
0054 HoughGrouping::HoughGrouping(const ParameterSet& pset, edm::ConsumesCollector& iC)
0055     : MotherGrouping(pset, iC), debug_(pset.getUntrackedParameter<bool>("debug")) {
0056   // Obtention of parameters
0057   if (debug_)
0058     LogDebug("HoughGrouping") << "HoughGrouping: constructor";
0059 
0060   // HOUGH TRANSFORM CONFIGURATION
0061   angletan_ = pset.getParameter<double>("angletan");
0062   anglebinwidth_ = pset.getParameter<double>("anglebinwidth");
0063   posbinwidth_ = pset.getParameter<double>("posbinwidth");
0064 
0065   // MAXIMA SEARCH CONFIGURATION
0066   maxdeltaAngDeg_ = pset.getParameter<double>("maxdeltaAngDeg");
0067   maxdeltaPos_ = pset.getParameter<double>("maxdeltaPos");
0068   upperNumber_ = (unsigned short int)pset.getParameter<int>("UpperNumber");
0069   lowerNumber_ = (unsigned short int)pset.getParameter<int>("LowerNumber");
0070 
0071   // HITS ASSOCIATION CONFIGURATION
0072   maxDistanceToWire_ = pset.getParameter<double>("MaxDistanceToWire");
0073 
0074   // CANDIDATE QUALITY CONFIGURATION
0075   minNLayerHits_ = (unsigned short int)pset.getParameter<int>("minNLayerHits");
0076   minSingleSLHitsMax_ = (unsigned short int)pset.getParameter<int>("minSingleSLHitsMax");
0077   minSingleSLHitsMin_ = (unsigned short int)pset.getParameter<int>("minSingleSLHitsMin");
0078   allowUncorrelatedPatterns_ = pset.getParameter<bool>("allowUncorrelatedPatterns");
0079   minUncorrelatedHits_ = (unsigned short int)pset.getParameter<int>("minUncorrelatedHits");
0080 
0081   dtGeomH = iC.esConsumes<DTGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0082 }
0083 
0084 HoughGrouping::~HoughGrouping() {
0085   if (debug_)
0086     LogDebug("HoughGrouping") << "HoughGrouping: destructor" << endl;
0087 }
0088 
0089 // ============================================================================
0090 // Main methods (initialise, run, finish)
0091 // ============================================================================
0092 void HoughGrouping::initialise(const edm::EventSetup& iEventSetup) {
0093   if (debug_)
0094     LogDebug("HoughGrouping") << "initialise";
0095 
0096   resetAttributes();
0097 
0098   maxrads_ = 0.5 * M_PI - atan(angletan_);
0099   minangle_ = (double)convertDegToRad(anglebinwidth_);
0100   halfanglebins_ = round(maxrads_ / minangle_ + 1);
0101   anglebins_ = (unsigned short int)2 * halfanglebins_;
0102   oneanglebin_ = maxrads_ / halfanglebins_;
0103 
0104   maxdeltaAng_ = maxdeltaAngDeg_ * 2 * M_PI / 360;
0105 
0106   // Initialisation of anglemap. Posmap depends on the size of the chamber.
0107   double phi = 0;
0108   anglemap_ = {};
0109   for (unsigned short int ab = 0; ab < halfanglebins_; ab++) {
0110     anglemap_[ab] = phi;
0111     phi += oneanglebin_;
0112   }
0113 
0114   phi = (M_PI - maxrads_);
0115   for (unsigned short int ab = halfanglebins_; ab < anglebins_; ab++) {
0116     anglemap_[ab] = phi;
0117     phi += oneanglebin_;
0118   }
0119 
0120   if (debug_) {
0121     LogDebug("HoughGrouping")
0122         << "\nHoughGrouping::ResetAttributes - Information from the initialisation of HoughGrouping:";
0123     LogDebug("HoughGrouping") << "ResetAttributes - maxrads: " << maxrads_;
0124     LogDebug("HoughGrouping") << "ResetAttributes - anglebinwidth: " << anglebinwidth_;
0125     LogDebug("HoughGrouping") << "ResetAttributes - minangle: " << minangle_;
0126     LogDebug("HoughGrouping") << "ResetAttributes - halfanglebins: " << halfanglebins_;
0127     LogDebug("HoughGrouping") << "ResetAttributes - anglebins: " << anglebins_;
0128     LogDebug("HoughGrouping") << "ResetAttributes - oneanglebin: " << oneanglebin_;
0129     LogDebug("HoughGrouping") << "ResetAttributes - posbinwidth: " << posbinwidth_;
0130   }
0131 
0132   const MuonGeometryRecord& geom = iEventSetup.get<MuonGeometryRecord>();
0133   dtGeo_ = &geom.get(dtGeomH);
0134 }
0135 
0136 void HoughGrouping::run(edm::Event& iEvent,
0137                         const edm::EventSetup& iEventSetup,
0138                         const DTDigiCollection& digis,
0139                         MuonPathPtrs& outMpath) {
0140   if (debug_)
0141     LogDebug("HoughGrouping") << "\nHoughGrouping::run";
0142 
0143   resetAttributes();
0144 
0145   if (debug_)
0146     LogDebug("HoughGrouping") << "run - Beginning digis' loop...";
0147   LocalPoint wirePosInLay, wirePosInChamber;
0148   GlobalPoint wirePosGlob;
0149   for (DTDigiCollection::DigiRangeIterator dtLayerIdIt = digis.begin(); dtLayerIdIt != digis.end(); dtLayerIdIt++) {
0150     const DTLayer* lay = dtGeo_->layer((*dtLayerIdIt).first);
0151     for (DTDigiCollection::const_iterator digiIt = ((*dtLayerIdIt).second).first;
0152          digiIt != ((*dtLayerIdIt).second).second;
0153          digiIt++) {
0154       if (debug_) {
0155         LogDebug("HoughGrouping") << "\nHoughGrouping::run - Digi number " << idigi_;
0156         LogDebug("HoughGrouping") << "run - Wheel: " << (*dtLayerIdIt).first.wheel();
0157         LogDebug("HoughGrouping") << "run - Chamber: " << (*dtLayerIdIt).first.station();
0158         LogDebug("HoughGrouping") << "run - Sector: " << (*dtLayerIdIt).first.sector();
0159         LogDebug("HoughGrouping") << "run - Superlayer: " << (*dtLayerIdIt).first.superLayer();
0160         LogDebug("HoughGrouping") << "run - Layer: " << (*dtLayerIdIt).first.layer();
0161         LogDebug("HoughGrouping") << "run - Wire: " << (*digiIt).wire();
0162         LogDebug("HoughGrouping") << "run - First wire: " << lay->specificTopology().firstChannel();
0163         LogDebug("HoughGrouping") << "run - Last wire: " << lay->specificTopology().lastChannel();
0164         LogDebug("HoughGrouping") << "run - First wire x: "
0165                                   << lay->specificTopology().wirePosition(lay->specificTopology().firstChannel());
0166         LogDebug("HoughGrouping") << "run - Last wire x: "
0167                                   << lay->specificTopology().wirePosition(lay->specificTopology().lastChannel());
0168         LogDebug("HoughGrouping") << "run - Cell width: " << lay->specificTopology().cellWidth();
0169         LogDebug("HoughGrouping") << "run - Cell height: " << lay->specificTopology().cellHeight();
0170       }
0171       if ((*dtLayerIdIt).first.superLayer() == 2)
0172         continue;
0173 
0174       wirePosInLay = LocalPoint(lay->specificTopology().wirePosition((*digiIt).wire()), 0, 0);
0175       wirePosGlob = lay->toGlobal(wirePosInLay);
0176       wirePosInChamber = lay->chamber()->toLocal(wirePosGlob);
0177 
0178       if ((*dtLayerIdIt).first.superLayer() == 3) {
0179         digimap_[(*dtLayerIdIt).first.layer() + 3][(*digiIt).wire()] = DTPrimitive();
0180         digimap_[(*dtLayerIdIt).first.layer() + 3][(*digiIt).wire()].setTDCTimeStamp((*digiIt).time());
0181         digimap_[(*dtLayerIdIt).first.layer() + 3][(*digiIt).wire()].setChannelId((*digiIt).wire());
0182         digimap_[(*dtLayerIdIt).first.layer() + 3][(*digiIt).wire()].setLayerId((*dtLayerIdIt).first.layer());
0183         digimap_[(*dtLayerIdIt).first.layer() + 3][(*digiIt).wire()].setSuperLayerId((*dtLayerIdIt).first.superLayer());
0184         digimap_[(*dtLayerIdIt).first.layer() + 3][(*digiIt).wire()].setCameraId((*dtLayerIdIt).first.rawId());
0185       } else {
0186         digimap_[(*dtLayerIdIt).first.layer() - 1][(*digiIt).wire()] = DTPrimitive();
0187         digimap_[(*dtLayerIdIt).first.layer() - 1][(*digiIt).wire()].setTDCTimeStamp((*digiIt).time());
0188         digimap_[(*dtLayerIdIt).first.layer() - 1][(*digiIt).wire()].setChannelId((*digiIt).wire());
0189         digimap_[(*dtLayerIdIt).first.layer() - 1][(*digiIt).wire()].setLayerId((*dtLayerIdIt).first.layer());
0190         digimap_[(*dtLayerIdIt).first.layer() - 1][(*digiIt).wire()].setSuperLayerId((*dtLayerIdIt).first.superLayer());
0191         digimap_[(*dtLayerIdIt).first.layer() - 1][(*digiIt).wire()].setCameraId((*dtLayerIdIt).first.rawId());
0192       }
0193 
0194       // Obtaining geometrical info of the chosen chamber
0195       if (xlowlim_ == 0 && xhighlim_ == 0 && zlowlim_ == 0 && zhighlim_ == 0) {
0196         thewheel_ = (*dtLayerIdIt).first.wheel();
0197         thestation_ = (*dtLayerIdIt).first.station();
0198         thesector_ = (*dtLayerIdIt).first.sector();
0199         obtainGeometricalBorders(lay);
0200       }
0201 
0202       if (debug_) {
0203         LogDebug("HoughGrouping") << "run - X position of the cell (chamber frame of reference): "
0204                                   << wirePosInChamber.x();
0205         LogDebug("HoughGrouping") << "run - Y position of the cell (chamber frame of reference): "
0206                                   << wirePosInChamber.y();
0207         LogDebug("HoughGrouping") << "run - Z position of the cell (chamber frame of reference): "
0208                                   << wirePosInChamber.z();
0209       }
0210 
0211       hitvec_.push_back({wirePosInChamber.x() - 1.05, wirePosInChamber.z()});
0212       hitvec_.push_back({wirePosInChamber.x() + 1.05, wirePosInChamber.z()});
0213       nhits_ += 2;
0214 
0215       idigi_++;
0216     }
0217   }
0218 
0219   if (debug_) {
0220     LogDebug("HoughGrouping") << "\nHoughGrouping::run - nhits: " << nhits_;
0221     LogDebug("HoughGrouping") << "run - idigi: " << idigi_;
0222   }
0223 
0224   if (hitvec_.empty()) {
0225     if (debug_)
0226       LogDebug("HoughGrouping") << "run - No digis present in this chamber.";
0227     return;
0228   }
0229 
0230   // Perform the Hough transform of the inputs.
0231   doHoughTransform();
0232 
0233   // Obtain the maxima
0234   maxima_ = getMaximaVector();
0235   resetPosElementsOfLinespace();
0236 
0237   if (maxima_.empty()) {
0238     if (debug_)
0239       LogDebug("HoughGrouping") << "run - No good maxima found in this event.";
0240     return;
0241   }
0242 
0243   DTChamberId TheChambId(thewheel_, thestation_, thesector_);
0244   const DTChamber* TheChamb = dtGeo_->chamber(TheChambId);
0245   std::vector<ProtoCand> cands;
0246 
0247   for (unsigned short int ican = 0; ican < maxima_.size(); ican++) {
0248     if (debug_)
0249       LogDebug("HoughGrouping") << "\nHoughGrouping::run - Candidate number: " << ican;
0250     cands.push_back(associateHits(TheChamb, maxima_.at(ican).first, maxima_.at(ican).second));
0251   }
0252 
0253   // Now we filter them:
0254   orderAndFilter(cands, outMpath);
0255   if (debug_)
0256     LogDebug("HoughGrouping") << "run - now we have our muonpaths! It has " << outMpath.size() << " elements";
0257 
0258   cands.clear();
0259   return;
0260 }
0261 
0262 void HoughGrouping::finish() {
0263   if (debug_)
0264     LogDebug("HoughGrouping") << "finish";
0265   return;
0266 }
0267 
0268 // ============================================================================
0269 // Other methods
0270 // ============================================================================
0271 void HoughGrouping::resetAttributes() {
0272   if (debug_)
0273     LogDebug("HoughGrouping") << "ResetAttributes";
0274   // std::vector's:
0275   maxima_.clear();
0276   hitvec_.clear();
0277 
0278   // Integer-type variables:
0279   spacebins_ = 0;
0280   idigi_ = 0;
0281   nhits_ = 0;
0282   xlowlim_ = 0;
0283   xhighlim_ = 0;
0284   zlowlim_ = 0;
0285   zhighlim_ = 0;
0286   thestation_ = 0;
0287   thesector_ = 0;
0288   thewheel_ = 0;
0289 
0290   // Arrays:
0291   // NOTE: linespace array is treated and reset separately
0292 
0293   // Maps (dictionaries):
0294   posmap_.clear();
0295   for (unsigned short int abslay = 0; abslay < 8; abslay++)
0296     digimap_[abslay].clear();
0297 }
0298 
0299 void HoughGrouping::resetPosElementsOfLinespace() {
0300   if (debug_)
0301     LogDebug("HoughGrouping") << "ResetPosElementsOfLinespace";
0302   for (unsigned short int ab = 0; ab < anglebins_; ab++) {
0303     linespace_[ab].clear();
0304   }
0305   linespace_.clear();
0306 }
0307 
0308 void HoughGrouping::obtainGeometricalBorders(const DTLayer* lay) {
0309   if (debug_)
0310     LogDebug("HoughGrouping") << "ObtainGeometricalBorders";
0311   LocalPoint FirstWireLocal(lay->chamber()->superLayer(1)->layer(1)->specificTopology().wirePosition(
0312                                 lay->chamber()->superLayer(1)->layer(1)->specificTopology().firstChannel()),
0313                             0,
0314                             0);  // TAKING INFO FROM L1 OF SL1 OF THE CHOSEN CHAMBER
0315   GlobalPoint FirstWireGlobal = lay->chamber()->superLayer(1)->layer(1)->toGlobal(FirstWireLocal);
0316   LocalPoint FirstWireLocalCh = lay->chamber()->toLocal(FirstWireGlobal);
0317 
0318   LocalPoint LastWireLocal(lay->chamber()->superLayer(1)->layer(1)->specificTopology().wirePosition(
0319                                lay->chamber()->superLayer(1)->layer(1)->specificTopology().lastChannel()),
0320                            0,
0321                            0);
0322   GlobalPoint LastWireGlobal = lay->chamber()->superLayer(1)->layer(1)->toGlobal(LastWireLocal);
0323   LocalPoint LastWireLocalCh = lay->chamber()->toLocal(LastWireGlobal);
0324 
0325   //   unsigned short int upsl = thestation == 4 ? 2 : 3;
0326   unsigned short int upsl = thestation_ == 4 ? 3 : 3;
0327   if (debug_)
0328     LogDebug("HoughGrouping") << "ObtainGeometricalBorders - uppersuperlayer: " << upsl;
0329 
0330   LocalPoint FirstWireLocalUp(lay->chamber()->superLayer(upsl)->layer(4)->specificTopology().wirePosition(
0331                                   lay->chamber()->superLayer(upsl)->layer(4)->specificTopology().firstChannel()),
0332                               0,
0333                               0);  // TAKING INFO FROM L1 OF SL1 OF THE CHOSEN CHAMBER
0334   GlobalPoint FirstWireGlobalUp = lay->chamber()->superLayer(upsl)->layer(4)->toGlobal(FirstWireLocalUp);
0335   LocalPoint FirstWireLocalChUp = lay->chamber()->toLocal(FirstWireGlobalUp);
0336 
0337   xlowlim_ = FirstWireLocalCh.x() - lay->chamber()->superLayer(1)->layer(1)->specificTopology().cellWidth() / 2;
0338   xhighlim_ = LastWireLocalCh.x() + lay->chamber()->superLayer(1)->layer(1)->specificTopology().cellWidth() / 2;
0339   zlowlim_ = FirstWireLocalChUp.z() - lay->chamber()->superLayer(upsl)->layer(4)->specificTopology().cellHeight() / 2;
0340   zhighlim_ = LastWireLocalCh.z() + lay->chamber()->superLayer(1)->layer(1)->specificTopology().cellHeight() / 2;
0341 
0342   spacebins_ = round(std::abs(xhighlim_ - xlowlim_) / posbinwidth_);
0343 }
0344 
0345 void HoughGrouping::doHoughTransform() {
0346   if (debug_)
0347     LogDebug("HoughGrouping") << "DoHoughTransform";
0348   // First we want to obtain the number of bins in angle that we want. To do so, we will consider at first a maximum angle of
0349   // (in rad.) pi/2 - arctan(0.3) (i.e. ~73º) and a resolution (width of bin angle) of 2º.
0350 
0351   if (debug_) {
0352     LogDebug("HoughGrouping") << "DoHoughTransform - maxrads: " << maxrads_;
0353     LogDebug("HoughGrouping") << "DoHoughTransform - minangle: " << minangle_;
0354     LogDebug("HoughGrouping") << "DoHoughTransform - halfanglebins: " << halfanglebins_;
0355     LogDebug("HoughGrouping") << "DoHoughTransform - anglebins: " << anglebins_;
0356     LogDebug("HoughGrouping") << "DoHoughTransform - oneanglebin: " << oneanglebin_;
0357     LogDebug("HoughGrouping") << "DoHoughTransform - spacebins: " << spacebins_;
0358   }
0359 
0360   double rho = 0, phi = 0, sbx = 0;
0361   double lowinitsb = xlowlim_ + posbinwidth_ / 2;
0362 
0363   // Initialisation
0364   linespace_.resize(anglebins_, std::vector<unsigned short int>(spacebins_));
0365   for (unsigned short int ab = 0; ab < anglebins_; ab++) {
0366     sbx = lowinitsb;
0367     for (unsigned short int sb = 0; sb < spacebins_; sb++) {
0368       posmap_[sb] = sbx;
0369       linespace_[ab][sb] = 0;
0370       sbx += posbinwidth_;
0371     }
0372   }
0373 
0374   // Filling of the double array and actually doing the transform
0375   for (unsigned short int i = 0; i < hitvec_.size(); i++) {
0376     for (unsigned short int ab = 0; ab < anglebins_; ab++) {
0377       phi = anglemap_[ab];
0378       rho = hitvec_.at(i).first * cos(phi) + hitvec_.at(i).second * sin(phi);
0379       sbx = lowinitsb - posbinwidth_ / 2;
0380       for (unsigned short int sb = 0; sb < spacebins_; sb++) {
0381         if (rho < sbx) {
0382           linespace_[ab][sb]++;
0383           break;
0384         }
0385         sbx += posbinwidth_;
0386       }
0387     }
0388   }
0389 }
0390 
0391 PointsInPlane HoughGrouping::getMaximaVector() {
0392   if (debug_)
0393     LogDebug("HoughGrouping") << "getMaximaVector";
0394   PointTuples tmpvec;
0395 
0396   bool flagsearched = false;
0397   unsigned short int numbertolookat = upperNumber_;
0398 
0399   while (!flagsearched) {
0400     for (unsigned short int ab = 0; ab < anglebins_; ab++) {
0401       for (unsigned short int sb = 0; sb < spacebins_; sb++) {
0402         if (linespace_[ab][sb] >= numbertolookat)
0403           tmpvec.push_back({anglemap_[ab], posmap_[sb], linespace_[ab][sb]});
0404       }
0405     }
0406     if (((numbertolookat - 1) < lowerNumber_) || (!tmpvec.empty()))
0407       flagsearched = true;
0408     else
0409       numbertolookat--;
0410   }
0411 
0412   if (tmpvec.empty()) {
0413     if (debug_)
0414       LogDebug("HoughGrouping") << "GetMaximaVector - No maxima could be found";
0415     PointsInPlane finalvec;
0416     return finalvec;
0417   } else {
0418     PointsInPlane finalvec = findTheMaxima(tmpvec);
0419 
0420     // And now obtain the values of m and n of the lines.
0421     for (unsigned short int i = 0; i < finalvec.size(); i++)
0422       finalvec.at(i) = transformPair(finalvec.at(i));
0423     return finalvec;
0424   }
0425 }
0426 
0427 PointInPlane HoughGrouping::transformPair(const PointInPlane& inputpair) {
0428   if (debug_)
0429     LogDebug("HoughGrouping") << "transformPair";
0430   // input: (ang, pos); output: (m, n)
0431   if (inputpair.first == 0)
0432     return {1000, -1000 * inputpair.second};
0433   else
0434     return {-1. / tan(inputpair.first), inputpair.second / sin(inputpair.first)};
0435 }
0436 
0437 PointsInPlane HoughGrouping::findTheMaxima(PointTuples& inputvec) {
0438   if (debug_)
0439     LogDebug("HoughGrouping") << "findTheMaxima";
0440   bool fullyreduced = false;
0441   unsigned short int ind = 0;
0442 
0443   std::vector<unsigned short int> chosenvec;
0444   PointsInPlane resultvec;
0445   PointInPlane finalpair = {};
0446 
0447   if (debug_)
0448     LogDebug("HoughGrouping") << "findTheMaxima - prewhile";
0449   while (!fullyreduced) {
0450     if (debug_) {
0451       LogDebug("HoughGrouping") << "\nHoughGrouping::findTheMaxima - New iteration";
0452       LogDebug("HoughGrouping") << "findTheMaxima - inputvec size: " << inputvec.size();
0453       LogDebug("HoughGrouping") << "findTheMaxima - ind: " << ind;
0454       LogDebug("HoughGrouping") << "findTheMaxima - maximum deltaang: " << maxdeltaAng_
0455                                 << " and maximum deltapos: " << maxdeltaPos_;
0456     }
0457     chosenvec.clear();
0458     //calculate distances and check out the ones that are near
0459     if (debug_)
0460       LogDebug("HoughGrouping") << "findTheMaxima - Ours have " << get<2>(inputvec.at(ind))
0461                                 << " entries, ang.: " << get<0>(inputvec.at(ind))
0462                                 << " and pos.: " << get<1>(inputvec.at(ind));
0463 
0464     for (unsigned short int j = ind + 1; j < inputvec.size(); j++) {
0465       if (getTwoDelta(inputvec.at(ind), inputvec.at(j)).first <= maxdeltaAng_ &&
0466           getTwoDelta(inputvec.at(ind), inputvec.at(j)).second <= maxdeltaPos_) {
0467         chosenvec.push_back(j);
0468         if (debug_)
0469           LogDebug("HoughGrouping") << "findTheMaxima -     - Adding num.  " << j
0470                                     << " with deltaang: " << getTwoDelta(inputvec.at(ind), inputvec.at(j)).first
0471                                     << ", and deltapos: " << getTwoDelta(inputvec.at(ind), inputvec.at(j)).second
0472                                     << " and with " << get<2>(inputvec.at(j))
0473                                     << " entries, ang.: " << get<0>(inputvec.at(j))
0474                                     << " and pos.: " << get<1>(inputvec.at(j));
0475       } else if (debug_)
0476         LogDebug("HoughGrouping") << "findTheMaxima -     - Ignoring num. " << j
0477                                   << " with deltaang: " << getTwoDelta(inputvec.at(ind), inputvec.at(j)).first
0478                                   << ", and deltapos: " << getTwoDelta(inputvec.at(ind), inputvec.at(j)).second
0479                                   << " and with " << get<2>(inputvec.at(j)) << " entries.";
0480     }
0481 
0482     if (debug_)
0483       LogDebug("HoughGrouping") << "findTheMaxima - chosenvecsize: " << chosenvec.size();
0484 
0485     if (chosenvec.empty()) {
0486       if (ind + 1 >= (unsigned short int)inputvec.size())
0487         fullyreduced = true;
0488       if ((get<0>(inputvec.at(ind)) <= maxrads_) || (get<0>(inputvec.at(ind)) >= M_PI - maxrads_))
0489         resultvec.push_back({get<0>(inputvec.at(ind)), get<1>(inputvec.at(ind))});
0490       else if (debug_)
0491         LogDebug("HoughGrouping") << "findTheMaxima -     - Candidate dropped due to an excess in angle";
0492       ind++;
0493       continue;
0494     }
0495 
0496     // Now average them
0497     finalpair = getAveragePoint(inputvec, ind, chosenvec);
0498 
0499     // Erase the ones you used
0500     inputvec.erase(inputvec.begin() + ind);
0501     for (short int j = chosenvec.size() - 1; j > -1; j--) {
0502       if (debug_)
0503         LogDebug("HoughGrouping") << "findTheMaxima - erasing index: " << chosenvec.at(j) - 1;
0504       inputvec.erase(inputvec.begin() + chosenvec.at(j) - 1);
0505     }
0506 
0507     if (debug_)
0508       LogDebug("HoughGrouping") << "findTheMaxima - inputvec size: " << inputvec.size();
0509 
0510     // And add the one you calculated:
0511     if ((finalpair.first <= maxrads_) || (finalpair.first >= M_PI - maxrads_))
0512       resultvec.push_back(finalpair);
0513     else if (debug_)
0514       LogDebug("HoughGrouping") << "findTheMaxima -     - Candidate dropped due to an excess in angle";
0515 
0516     if (ind + 1 >= (unsigned short int)inputvec.size())
0517       fullyreduced = true;
0518     if (debug_)
0519       LogDebug("HoughGrouping") << "findTheMaxima - iteration ends";
0520     ind++;
0521   }
0522   if (debug_)
0523     LogDebug("HoughGrouping") << "findTheMaxima - postwhile";
0524   return resultvec;
0525 }
0526 
0527 PointInPlane HoughGrouping::getTwoDelta(const PointTuple& pair1, const PointTuple& pair2) {
0528   if (debug_)
0529     LogDebug("HoughGrouping") << "getTwoDelta";
0530   return {abs(get<0>(pair1) - get<0>(pair2)), abs(get<1>(pair1) - get<1>(pair2))};
0531 }
0532 
0533 PointInPlane HoughGrouping::getAveragePoint(const PointTuples& inputvec,
0534                                             unsigned short int firstindex,
0535                                             const std::vector<unsigned short int>& indexlist) {
0536   if (debug_)
0537     LogDebug("HoughGrouping") << "getAveragePoint";
0538   std::vector<double> xs;
0539   std::vector<double> ys;
0540   std::vector<double> ws;
0541   xs.push_back(get<0>(inputvec.at(firstindex)));
0542   ys.push_back(get<1>(inputvec.at(firstindex)));
0543   ws.push_back(exp(get<2>(inputvec.at(firstindex))));
0544   for (unsigned short int i = 0; i < indexlist.size(); i++) {
0545     xs.push_back(get<0>(inputvec.at(indexlist.at(i))));
0546     ys.push_back(get<1>(inputvec.at(indexlist.at(i))));
0547     ws.push_back(exp(get<2>(inputvec.at(indexlist.at(i)))));
0548   }
0549   return {TMath::Mean(xs.begin(), xs.end(), ws.begin()), TMath::Mean(ys.begin(), ys.end(), ws.begin())};
0550 }
0551 
0552 ProtoCand HoughGrouping::associateHits(const DTChamber* thechamb, double m, double n) {
0553   if (debug_)
0554     LogDebug("HoughGrouping") << "associateHits";
0555   LocalPoint tmpLocal, AWireLocal, AWireLocalCh, tmpLocalCh, thepoint;
0556   GlobalPoint tmpGlobal, AWireGlobal;
0557   double tmpx = 0;
0558   unsigned short int tmpwire = 0;
0559   unsigned short int abslay = 0;
0560   LATERAL_CASES lat = NONE;
0561   bool isleft = false;
0562   bool isright = false;
0563 
0564   ProtoCand returnPC;
0565   for (auto l = 0; l < NUM_LAYERS_2SL; l++) {
0566     returnPC.isThereHitInLayer_.push_back(false);
0567     returnPC.isThereNeighBourHitInLayer_.push_back(false);
0568     returnPC.xDistToPattern_.push_back(0);
0569     returnPC.dtHits_.push_back(DTPrimitive());
0570   }
0571 
0572   if (debug_)
0573     LogDebug("HoughGrouping") << "associateHits - Beginning SL loop";
0574   for (unsigned short int sl = 1; sl < 3 + 1; sl++) {
0575     if (sl == 2)
0576       continue;
0577     if (debug_)
0578       LogDebug("HoughGrouping") << "associateHits - SL: " << sl;
0579 
0580     for (unsigned short int l = 1; l < 4 + 1; l++) {
0581       if (debug_)
0582         LogDebug("HoughGrouping") << "associateHits - L: " << l;
0583       isleft = false;
0584       isright = false;
0585       lat = NONE;
0586       if (sl == 1)
0587         abslay = l - 1;
0588       else
0589         abslay = l + 3;
0590       AWireLocal = LocalPoint(thechamb->superLayer(sl)->layer(l)->specificTopology().wirePosition(
0591                                   thechamb->superLayer(sl)->layer(l)->specificTopology().firstChannel()),
0592                               0,
0593                               0);
0594       AWireGlobal = thechamb->superLayer(sl)->layer(l)->toGlobal(AWireLocal);
0595       AWireLocalCh = thechamb->toLocal(AWireGlobal);
0596       tmpx = (AWireLocalCh.z() - n) / m;
0597 
0598       if ((tmpx <= xlowlim_) || (tmpx >= xhighlim_)) {
0599         returnPC.dtHits_[abslay] = DTPrimitive();  // empty primitive
0600         continue;
0601       }
0602 
0603       thepoint = LocalPoint(tmpx, 0, AWireLocalCh.z());
0604       tmpwire = thechamb->superLayer(sl)->layer(l)->specificTopology().channel(thepoint);
0605       if (debug_)
0606         LogDebug("HoughGrouping") << "associateHits - Wire number: " << tmpwire;
0607       if (debug_)
0608         LogDebug("HoughGrouping") << "associateHits - First channel in layer: "
0609                                   << thechamb->superLayer(sl)->layer(l)->specificTopology().firstChannel();
0610       if ((digimap_[abslay]).count(tmpwire)) {
0611         // OK, we have a digi, let's choose the laterality, if we can:
0612         tmpLocal = LocalPoint(thechamb->superLayer(sl)->layer(l)->specificTopology().wirePosition(tmpwire), 0, 0);
0613         tmpGlobal = thechamb->superLayer(sl)->layer(l)->toGlobal(tmpLocal);
0614         tmpLocalCh = thechamb->toLocal(tmpGlobal);
0615 
0616         if (abs(tmpLocalCh.x() - thepoint.x()) >= maxDistanceToWire_) {
0617           // The distance where lateralities are not put is 0.03 cm, which is a conservative threshold for the resolution of the cells.
0618           if ((tmpLocalCh.x() - thepoint.x()) > 0)
0619             lat = LEFT;
0620           else
0621             lat = RIGHT;
0622         }
0623 
0624         // Filling info
0625         returnPC.nLayersWithHits_++;
0626         returnPC.isThereHitInLayer_[abslay] = true;
0627         returnPC.isThereNeighBourHitInLayer_[abslay] = true;
0628         if (lat == LEFT)
0629           returnPC.xDistToPattern_[abslay] = abs(tmpx - (tmpLocalCh.x() - 1.05));
0630         else if (lat == RIGHT)
0631           returnPC.xDistToPattern_[abslay] = abs(tmpx - (tmpLocalCh.x() + 1.05));
0632         else
0633           returnPC.xDistToPattern_[abslay] = abs(tmpx - tmpLocalCh.x());
0634         returnPC.dtHits_[abslay] = DTPrimitive(digimap_[abslay][tmpwire]);
0635         returnPC.dtHits_[abslay].setLaterality(lat);
0636       } else {
0637         if (debug_)
0638           LogDebug("HoughGrouping") << "associateHits - No hit in the crossing cell";
0639         if ((digimap_[abslay]).count(tmpwire - 1))
0640           isleft = true;
0641         if ((digimap_[abslay]).count(tmpwire + 1))
0642           isright = true;
0643         if (debug_)
0644           LogDebug("HoughGrouping") << "associateHits - There is in the left: " << (int)isleft;
0645         if (debug_)
0646           LogDebug("HoughGrouping") << "associateHits - There is in the right: " << (int)isright;
0647 
0648         if ((isleft) && (!isright)) {
0649           tmpLocal = LocalPoint(thechamb->superLayer(sl)->layer(l)->specificTopology().wirePosition(tmpwire - 1), 0, 0);
0650           tmpGlobal = thechamb->superLayer(sl)->layer(l)->toGlobal(tmpLocal);
0651           tmpLocalCh = thechamb->toLocal(tmpGlobal);
0652 
0653           // Filling info
0654           returnPC.nLayersWithHits_++;
0655           returnPC.isThereNeighBourHitInLayer_[abslay] = true;
0656           returnPC.xDistToPattern_[abslay] = abs(tmpx - (tmpLocalCh.x() + 1.05));
0657           returnPC.dtHits_[abslay] = DTPrimitive(digimap_[abslay][tmpwire - 1]);
0658           returnPC.dtHits_[abslay].setLaterality(RIGHT);
0659 
0660         } else if ((!isleft) && (isright)) {
0661           tmpLocal = LocalPoint(thechamb->superLayer(sl)->layer(l)->specificTopology().wirePosition(tmpwire + 1), 0, 0);
0662           tmpGlobal = thechamb->superLayer(sl)->layer(l)->toGlobal(tmpLocal);
0663           tmpLocalCh = thechamb->toLocal(tmpGlobal);
0664 
0665           // Filling info
0666           returnPC.nLayersWithHits_++;
0667           returnPC.isThereNeighBourHitInLayer_[abslay] = true;
0668           returnPC.xDistToPattern_[abslay] = abs(tmpx - (tmpLocalCh.x() - 1.05));
0669           returnPC.dtHits_[abslay] = DTPrimitive(digimap_[abslay][tmpwire + 1]);
0670           returnPC.dtHits_[abslay].setLaterality(LEFT);
0671         } else if ((isleft) && (isright)) {
0672           LocalPoint tmpLocal_l =
0673               LocalPoint(thechamb->superLayer(sl)->layer(l)->specificTopology().wirePosition(tmpwire - 1), 0, 0);
0674           GlobalPoint tmpGlobal_l = thechamb->superLayer(sl)->layer(l)->toGlobal(tmpLocal_l);
0675           LocalPoint tmpLocalCh_l = thechamb->toLocal(tmpGlobal_l);
0676 
0677           LocalPoint tmpLocal_r =
0678               LocalPoint(thechamb->superLayer(sl)->layer(l)->specificTopology().wirePosition(tmpwire + 1), 0, 0);
0679           GlobalPoint tmpGlobal_r = thechamb->superLayer(sl)->layer(l)->toGlobal(tmpLocal_r);
0680           LocalPoint tmpLocalCh_r = thechamb->toLocal(tmpGlobal_r);
0681 
0682           // Filling info
0683           returnPC.nLayersWithHits_++;
0684           returnPC.isThereNeighBourHitInLayer_[abslay] = true;
0685 
0686           bool isDistRight = std::abs(thepoint.x() - tmpLocalCh_l.x()) < std::abs(thepoint.x() - tmpLocalCh_r.x());
0687           if (isDistRight) {
0688             returnPC.xDistToPattern_[abslay] = std::abs(tmpx - (tmpLocalCh.x() + 1.05));
0689             returnPC.dtHits_[abslay] = DTPrimitive(digimap_[abslay][tmpwire - 1]);
0690             returnPC.dtHits_[abslay].setLaterality(RIGHT);
0691           } else {
0692             returnPC.xDistToPattern_[abslay] = std::abs(tmpx - (tmpLocalCh.x() - 1.05));
0693             returnPC.dtHits_[abslay] = DTPrimitive(digimap_[abslay][tmpwire + 1]);
0694             returnPC.dtHits_[abslay].setLaterality(LEFT);
0695           }
0696         } else {                                     // case where there are no digis
0697           returnPC.dtHits_[abslay] = DTPrimitive();  // empty primitive
0698         }
0699       }
0700     }
0701   }
0702 
0703   setDifferenceBetweenSL(returnPC);
0704   if (debug_) {
0705     LogDebug("HoughGrouping") << "associateHits - Finishing with the candidate. We have found the following of it:";
0706     LogDebug("HoughGrouping") << "associateHits - # of layers with hits: " << returnPC.nLayersWithHits_;
0707     for (unsigned short int lay = 0; lay < 8; lay++) {
0708       LogDebug("HoughGrouping") << "associateHits - For absolute layer: " << lay;
0709       LogDebug("HoughGrouping") << "associateHits - # of HQ hits: " << returnPC.isThereHitInLayer_[lay];
0710       LogDebug("HoughGrouping") << "associateHits - # of LQ hits: " << returnPC.isThereNeighBourHitInLayer_[lay];
0711     }
0712     LogDebug("HoughGrouping") << "associateHits - Abs. diff. between SL1 and SL3 hits: " << returnPC.nHitsDiff_;
0713     for (unsigned short int lay = 0; lay < 8; lay++) {
0714       LogDebug("HoughGrouping") << "associateHits - For absolute layer: " << lay;
0715       LogDebug("HoughGrouping") << "associateHits - Abs. distance to digi: " << returnPC.xDistToPattern_[lay];
0716     }
0717   }
0718   return returnPC;
0719 }
0720 
0721 void HoughGrouping::setDifferenceBetweenSL(ProtoCand& tupl) {
0722   if (debug_)
0723     LogDebug("HoughGrouping") << "setDifferenceBetweenSL";
0724   short int absres = 0;
0725   for (unsigned short int lay = 0; lay < 8; lay++) {
0726     if (tupl.dtHits_[lay].channelId() > 0) {
0727       if (lay <= 3)
0728         absres++;
0729       else
0730         absres--;
0731     }
0732   }
0733 
0734   if (absres >= 0)
0735     tupl.nHitsDiff_ = absres;
0736   else
0737     tupl.nHitsDiff_ = (unsigned short int)(-absres);
0738 }
0739 
0740 void HoughGrouping::orderAndFilter(std::vector<ProtoCand>& invector, MuonPathPtrs& outMuonPath) {
0741   if (debug_)
0742     LogDebug("HoughGrouping") << "orderAndFilter";
0743   // 0: # of layers with hits.
0744   // 1: # of hits of high quality (the expected line crosses the cell).
0745   // 2: # of hits of low quality (the expected line is in a neighbouring cell).
0746   // 3: absolute diff. between the number of hits in SL1 and SL3.
0747   // 4: absolute distance to all hits of the segment.
0748   // 5: DTPrimitive of the candidate.
0749 
0750   std::vector<unsigned short int> elstoremove;
0751   // ordering:
0752   if (debug_)
0753     LogDebug("HoughGrouping") << "orderAndFilter - First ordering";
0754   std::sort(invector.begin(), invector.end(), HoughOrdering);
0755 
0756   // Now filtering:
0757   unsigned short int ind = 0;
0758   bool filtered = false;
0759   if (debug_)
0760     LogDebug("HoughGrouping") << "orderAndFilter - Entering while";
0761   while (!filtered) {
0762     if (debug_)
0763       LogDebug("HoughGrouping") << "\nHoughGrouping::orderAndFilter - New iteration with ind: " << ind;
0764     elstoremove.clear();
0765     for (unsigned short int i = ind + 1; i < invector.size(); i++) {
0766       if (debug_)
0767         LogDebug("HoughGrouping") << "orderAndFilter - Checking index: " << i;
0768       for (unsigned short int lay = 0; lay < 8; lay++) {
0769         if (debug_)
0770           LogDebug("HoughGrouping") << "orderAndFilter - Checking layer number: " << lay;
0771         if (invector.at(i).dtHits_[lay].channelId() == invector.at(ind).dtHits_[lay].channelId() &&
0772             invector.at(ind).dtHits_[lay].channelId() != -1) {
0773           invector.at(i).nLayersWithHits_--;
0774           invector.at(i).isThereHitInLayer_[lay] = false;
0775           invector.at(i).isThereNeighBourHitInLayer_[lay] = false;
0776           setDifferenceBetweenSL(invector.at(i));
0777           // We check that if its a different laterality, the best candidate of the two of them changes its laterality to not-known (that is, both).
0778           if (invector.at(i).dtHits_[lay].laterality() != invector.at(ind).dtHits_[lay].laterality())
0779             invector.at(ind).dtHits_[lay].setLaterality(NONE);
0780 
0781           invector.at(i).dtHits_[lay] = DTPrimitive();
0782         }
0783       }
0784       if (debug_)
0785         LogDebug("HoughGrouping")
0786             << "orderAndFilter - Finished checking all the layers, now seeing if we should remove the "
0787                "candidate";
0788 
0789       if (!areThereEnoughHits(invector.at(i))) {
0790         if (debug_)
0791           LogDebug("HoughGrouping") << "orderAndFilter - This candidate shall be removed!";
0792         elstoremove.push_back((unsigned short int)i);
0793       }
0794     }
0795 
0796     if (debug_)
0797       LogDebug("HoughGrouping") << "orderAndFilter - We are gonna erase " << elstoremove.size() << " elements";
0798 
0799     for (short int el = (elstoremove.size() - 1); el > -1; el--) {
0800       invector.erase(invector.begin() + elstoremove.at(el));
0801     }
0802 
0803     if (ind + 1 == (unsigned short int)invector.size())
0804       filtered = true;
0805     else
0806       std::sort(invector.begin() + ind + 1, invector.end(), HoughOrdering);
0807     ind++;
0808   }
0809 
0810   // Ultimate filter: if the remaining do not fill the requirements (configurable through pset arguments), they are removed also.
0811   for (short int el = (invector.size() - 1); el > -1; el--) {
0812     if (!areThereEnoughHits(invector.at(el))) {
0813       invector.erase(invector.begin() + el);
0814     }
0815   }
0816 
0817   if (invector.empty()) {
0818     if (debug_)
0819       LogDebug("HoughGrouping") << "OrderAndFilter - We do not have candidates with the minimum hits required.";
0820     return;
0821   } else if (debug_)
0822     LogDebug("HoughGrouping") << "OrderAndFilter - At the end, we have only " << invector.size() << " good paths!";
0823 
0824   // Packing dt primitives
0825   for (unsigned short int i = 0; i < invector.size(); i++) {
0826     DTPrimitivePtrs ptrPrimitive;
0827     unsigned short int tmplowfill = 0;
0828     unsigned short int tmpupfill = 0;
0829     for (unsigned short int lay = 0; lay < 8; lay++) {
0830       auto dtAux = std::make_shared<DTPrimitive>(invector.at(i).dtHits_[lay]);
0831       ptrPrimitive.push_back(std::move(dtAux));
0832       if (debug_) {
0833         LogDebug("HoughGrouping") << "\nHoughGrouping::OrderAndFilter - cameraid: " << ptrPrimitive[lay]->cameraId();
0834         LogDebug("HoughGrouping") << "OrderAndFilter - channelid (GOOD): " << ptrPrimitive[lay]->channelId();
0835         LogDebug("HoughGrouping") << "OrderAndFilter - channelid (AM):   " << ptrPrimitive[lay]->channelId() - 1;
0836       }
0837       // Fixing channel ID to AM conventions...
0838       if (ptrPrimitive[lay]->channelId() != -1)
0839         ptrPrimitive[lay]->setChannelId(ptrPrimitive[lay]->channelId() - 1);
0840 
0841       if (ptrPrimitive[lay]->cameraId() > 0) {
0842         if (lay < 4)
0843           tmplowfill++;
0844         else
0845           tmpupfill++;
0846       }
0847     }
0848 
0849     auto ptrMuonPath = std::make_shared<MuonPath>(ptrPrimitive, tmplowfill, tmpupfill);
0850     outMuonPath.push_back(ptrMuonPath);
0851     if (debug_) {
0852       for (unsigned short int lay = 0; lay < 8; lay++) {
0853         LogDebug("HoughGrouping") << "OrderAndFilter - Final cameraID: "
0854                                   << outMuonPath.back()->primitive(lay)->cameraId();
0855         LogDebug("HoughGrouping") << "OrderAndFilter - Final channelID: "
0856                                   << outMuonPath.back()->primitive(lay)->channelId();
0857         LogDebug("HoughGrouping") << "OrderAndFilter - Final time: "
0858                                   << outMuonPath.back()->primitive(lay)->tdcTimeStamp();
0859       }
0860     }
0861   }
0862   return;
0863 }
0864 
0865 bool HoughGrouping::areThereEnoughHits(const ProtoCand& tupl) {
0866   if (debug_)
0867     LogDebug("HoughGrouping") << "areThereEnoughHits";
0868   unsigned short int numhitssl1 = 0;
0869   unsigned short int numhitssl3 = 0;
0870   for (unsigned short int lay = 0; lay < 8; lay++) {
0871     if ((tupl.dtHits_[lay].channelId() > 0) && (lay < 4))
0872       numhitssl1++;
0873     else if (tupl.dtHits_[lay].channelId() > 0)
0874       numhitssl3++;
0875   }
0876 
0877   if (debug_)
0878     LogDebug("HoughGrouping") << "areThereEnoughHits - Hits in SL1: " << numhitssl1;
0879   if (debug_)
0880     LogDebug("HoughGrouping") << "areThereEnoughHits - Hits in SL3: " << numhitssl3;
0881 
0882   if ((numhitssl1 != 0) && (numhitssl3 != 0)) {  // Correlated candidates
0883     if ((numhitssl1 + numhitssl3) >= minNLayerHits_) {
0884       if (numhitssl1 > numhitssl3) {
0885         return ((numhitssl1 >= minSingleSLHitsMax_) && (numhitssl3 >= minSingleSLHitsMin_));
0886       } else if (numhitssl3 > numhitssl1) {
0887         return ((numhitssl3 >= minSingleSLHitsMax_) && (numhitssl1 >= minSingleSLHitsMin_));
0888       } else
0889         return true;
0890     }
0891   } else if (allowUncorrelatedPatterns_) {  // Uncorrelated candidates
0892     return ((numhitssl1 + numhitssl3) >= minNLayerHits_);
0893   }
0894   return false;
0895 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!