Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:36 

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackAssociator
0004 // Class:      DetIdAssociator
0005 //
0006 /*
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Dmytro Kovalskyi
0015 //         Created:  Fri Apr 21 10:59:41 PDT 2006
0016 //
0017 //
0018 
0019 #include "TrackingTools/TrackAssociator/interface/DetIdAssociator.h"
0020 #include "DetIdInfo.h"
0021 #include "FWCore/Utilities/interface/isFinite.h"
0022 #include <map>
0023 
0024 DetIdAssociator::DetIdAssociator(const int nPhi, const int nEta, const double etaBinSize)
0025     : nPhi_(nPhi),
0026       nEta_(nEta),
0027       lookupMap_(nPhi_ * nEta_, std::pair<unsigned int, unsigned int>(0, 0)),
0028       theMapIsValid_(false),
0029       etaBinSize_(etaBinSize) {
0030   if (nEta_ <= 0 || nPhi_ <= 0)
0031     throw cms::Exception("FatalError") << "incorrect look-up map size. Cannot initialize such a map.";
0032   maxEta_ = etaBinSize_ * nEta_ / 2;
0033   minTheta_ = 2 * atan(exp(-maxEta_));
0034 }
0035 
0036 std::set<DetId> DetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& direction, const int iN) const {
0037   unsigned int n = 0;
0038   if (iN > 0)
0039     n = iN;
0040   return getDetIdsCloseToAPoint(direction, n, n, n, n);
0041 }
0042 
0043 std::set<DetId> DetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& direction,
0044                                                         const unsigned int iNEtaPlus,
0045                                                         const unsigned int iNEtaMinus,
0046                                                         const unsigned int iNPhiPlus,
0047                                                         const unsigned int iNPhiMinus) const {
0048   std::set<DetId> set;
0049   check_setup();
0050   if (!theMapIsValid_)
0051     throw cms::Exception("FatalError") << "map is not valid.";
0052   LogTrace("TrackAssociator") << "(iNEtaPlus, iNEtaMinus, iNPhiPlus, iNPhiMinus): " << iNEtaPlus << ", " << iNEtaMinus
0053                               << ", " << iNPhiPlus << ", " << iNPhiMinus;
0054   LogTrace("TrackAssociator") << "point (eta,phi): " << direction.eta() << "," << direction.phi();
0055   int ieta = iEta(direction);
0056   int iphi = iPhi(direction);
0057   LogTrace("TrackAssociator") << "(ieta,iphi): " << ieta << "," << iphi << "\n";
0058   if (ieta >= 0 && ieta < nEta_ && iphi >= 0 && iphi < nPhi_) {
0059     fillSet(set, ieta, iphi);
0060     // dumpMapContent(ieta,iphi);
0061     // check if any neighbor bin is requested
0062     if (iNEtaPlus + iNEtaMinus + iNPhiPlus + iNPhiMinus > 0) {
0063       LogTrace("TrackAssociator") << "Add neighbors (ieta,iphi): " << ieta << "," << iphi;
0064       // eta
0065       int maxIEta = ieta + iNEtaPlus;
0066       int minIEta = ieta - iNEtaMinus;
0067       if (maxIEta >= nEta_)
0068         maxIEta = nEta_ - 1;
0069       if (minIEta < 0)
0070         minIEta = 0;
0071       // phi
0072       int maxIPhi = iphi + iNPhiPlus;
0073       int minIPhi = iphi - iNPhiMinus;
0074       if (maxIPhi - minIPhi >= nPhi_) {  // all elements in phi
0075         minIPhi = 0;
0076         maxIPhi = nPhi_ - 1;
0077       }
0078       if (minIPhi < 0) {
0079         minIPhi += nPhi_;
0080         maxIPhi += nPhi_;
0081       }
0082       LogTrace("TrackAssociator") << "\tieta (min,max): " << minIEta << "," << maxIEta;
0083       LogTrace("TrackAssociator") << "\tiphi (min,max): " << minIPhi << "," << maxIPhi << "\n";
0084       // dumpMapContent(minIEta,maxIEta,minIPhi,maxIPhi);
0085       for (int i = minIEta; i <= maxIEta; i++)
0086         for (int j = minIPhi; j <= maxIPhi; j++) {
0087           if (i == ieta && j == iphi)
0088             continue;  // already in the set
0089           fillSet(set, i, j % nPhi_);
0090         }
0091     }
0092   }
0093   return set;
0094 }
0095 
0096 std::set<DetId> DetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& point, const double d) const {
0097   return getDetIdsCloseToAPoint(point, d, d, d, d);
0098 }
0099 
0100 std::set<DetId> DetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& point,
0101                                                         const double dThetaPlus,
0102                                                         const double dThetaMinus,
0103                                                         const double dPhiPlus,
0104                                                         const double dPhiMinus) const {
0105   LogTrace("TrackAssociator") << "(dThetaPlus,dThetaMinus,dPhiPlus,dPhiMinus): " << dThetaPlus << ", " << dThetaMinus
0106                               << ", " << dPhiPlus << ", " << dPhiMinus;
0107   unsigned int n = 0;
0108   if (dThetaPlus < 0 || dThetaMinus < 0 || dPhiPlus < 0 || dPhiMinus < 0)
0109     return getDetIdsCloseToAPoint(point, n, n, n, n);
0110   // check that region of interest overlaps with the look-up map
0111   double maxTheta = point.theta() + dThetaPlus;
0112   if (maxTheta > M_PI - minTheta_)
0113     maxTheta = M_PI - minTheta_;
0114   double minTheta = point.theta() - dThetaMinus;
0115   if (minTheta < minTheta_)
0116     minTheta = minTheta_;
0117   if (maxTheta < minTheta_ || minTheta > M_PI - minTheta_)
0118     return std::set<DetId>();
0119 
0120   // take into account non-linear dependence of eta from
0121   // theta in regions with large |eta|
0122   double minEta = -log(tan(maxTheta / 2));
0123   double maxEta = -log(tan(minTheta / 2));
0124   unsigned int iNEtaPlus = abs(int((maxEta - point.eta()) / etaBinSize_));
0125   unsigned int iNEtaMinus = abs(int((point.eta() - minEta) / etaBinSize_));
0126   unsigned int iNPhiPlus = abs(int(dPhiPlus / (2 * M_PI) * nPhi_));
0127   unsigned int iNPhiMinus = abs(int(dPhiMinus / (2 * M_PI) * nPhi_));
0128   // add one more bin in each direction to guaranty that we don't miss anything
0129   return getDetIdsCloseToAPoint(point, iNEtaPlus + 1, iNEtaMinus + 1, iNPhiPlus + 1, iNPhiMinus + 1);
0130 }
0131 
0132 int DetIdAssociator::iEta(const GlobalPoint& point) const { return int(point.eta() / etaBinSize_ + nEta_ / 2); }
0133 
0134 int DetIdAssociator::iPhi(const GlobalPoint& point) const {
0135   return int((double(point.phi()) + M_PI) / (2 * M_PI) * nPhi_);
0136 }
0137 
0138 void DetIdAssociator::buildMap() {
0139   // the map is built in two steps to create a clean memory footprint
0140   // 0) determine how many elements each bin has
0141   // 1) fill the container map
0142   check_setup();
0143   LogTrace("TrackAssociator") << "building map for " << name() << "\n";
0144   // clear the map
0145   if (nEta_ <= 0 || nPhi_ <= 0)
0146     throw cms::Exception("FatalError") << "incorrect look-up map size. Cannot build such a map.";
0147   std::vector<GlobalPoint> pointBuffer;
0148   std::vector<DetId> detIdBuffer;
0149   unsigned int numberOfSubDetectors = getNumberOfSubdetectors();
0150   unsigned totalNumberOfElementsInTheContainer(0);
0151   for (unsigned int step = 0; step < 2; ++step) {
0152     unsigned int numberOfDetIdsOutsideEtaRange = 0;
0153     unsigned int numberOfDetIdsActive = 0;
0154     LogTrace("TrackAssociator") << "Step: " << step;
0155     for (unsigned int subDetectorIndex = 0; subDetectorIndex < numberOfSubDetectors; ++subDetectorIndex) {
0156       getValidDetIds(subDetectorIndex, detIdBuffer);
0157       // This std::move prevents modification
0158       std::vector<DetId> const& validIds = std::move(detIdBuffer);
0159       LogTrace("TrackAssociator") << "Number of valid DetIds for subdetector: " << subDetectorIndex << " is "
0160                                   << validIds.size();
0161       for (std::vector<DetId>::const_iterator id_itr = validIds.begin(); id_itr != validIds.end(); id_itr++) {
0162         std::pair<const_iterator, const_iterator> points = getDetIdPoints(*id_itr, pointBuffer);
0163         LogTrace("TrackAssociatorVerbose") << "Found " << points.second - points.first
0164                                            << " global points to describe geometry of DetId: " << id_itr->rawId();
0165         int etaMax(-1);
0166         int etaMin(-1);
0167         int phiMax(-1);
0168         int phiMin(-1);
0169         // this is a bit overkill, but it should be 100% proof (when debugged :)
0170         for (std::vector<GlobalPoint>::const_iterator iter = points.first; iter != points.second; iter++) {
0171           LogTrace("TrackAssociatorVerbose")
0172               << "\tpoint (rho,phi,z): " << iter->perp() << ", " << iter->phi() << ", " << iter->z();
0173           // FIX ME: this should be a fatal error
0174           if (edm::isNotFinite(iter->mag()) || iter->mag() > 1e5) {  //Detector parts cannot be 1 km away or be NaN
0175             edm::LogWarning("TrackAssociator")
0176                 << "Critical error! Bad detector unit geometry:\n\tDetId:" << id_itr->rawId()
0177                 << "\t mag(): " << iter->mag() << "\n"
0178                 << DetIdInfo::info(*id_itr, nullptr) << "\nSkipped the element";
0179             continue;
0180           }
0181           volume_.addActivePoint(*iter);
0182           int ieta = iEta(*iter);
0183           int iphi = iPhi(*iter);
0184           if (ieta < 0 || ieta >= nEta_) {
0185             LogTrace("TrackAssociator") << "Out of range: DetId:" << id_itr->rawId() << "\t (ieta,iphi): " << ieta
0186                                         << "," << iphi << "\n"
0187                                         << "Point: " << *iter << "\t(eta,phi): " << (*iter).eta() << ","
0188                                         << (*iter).phi() << "\n center: " << getPosition(*id_itr);
0189             continue;
0190           }
0191           if (phiMin < 0) {
0192             // first element
0193             etaMin = ieta;
0194             etaMax = ieta;
0195             phiMin = iphi;
0196             phiMax = iphi;
0197           } else {
0198             // check for discontinuity in phi
0199             int deltaMin = abs(phiMin - iphi);
0200             int deltaMax = abs(phiMax - iphi);
0201             // assume that no single detector element has more than 3.1416 coverage in phi
0202             if (deltaMin > nPhi_ / 2 && phiMin < nPhi_ / 2)
0203               phiMin += nPhi_;
0204             if (deltaMax > nPhi_ / 2) {
0205               if (phiMax < nPhi_ / 2)
0206                 phiMax += nPhi_;
0207               else
0208                 iphi += nPhi_;
0209             }
0210             assert(iphi >= 0);
0211             if (etaMin > ieta)
0212               etaMin = ieta;
0213             if (etaMax < ieta)
0214               etaMax = ieta;
0215             if (phiMin > iphi)
0216               phiMin = iphi;
0217             if (phiMax < iphi)
0218               phiMax = iphi;
0219           }
0220         }
0221         if (etaMax < 0 || phiMax < 0 || etaMin >= nEta_ || phiMin >= nPhi_) {
0222           LogTrace("TrackAssociatorVerbose")
0223               << "Out of range or no geometry: DetId:" << id_itr->rawId() << "\n\teta (min,max): " << etaMin << ","
0224               << etaMax << "\n\tphi (min,max): " << phiMin << "," << phiMax << "\nTower id: " << id_itr->rawId()
0225               << "\n";
0226           numberOfDetIdsOutsideEtaRange++;
0227           continue;
0228         }
0229         numberOfDetIdsActive++;
0230 
0231         LogTrace("TrackAssociatorVerbose") << "DetId (ieta_min,ieta_max,iphi_min,iphi_max): " << id_itr->rawId() << ", "
0232                                            << etaMin << ", " << etaMax << ", " << phiMin << ", " << phiMax;
0233         for (int ieta = etaMin; ieta <= etaMax; ieta++)
0234           for (int iphi = phiMin; iphi <= phiMax; iphi++)
0235             if (step == 0) {
0236               lookupMap_.at(index(ieta, iphi % nPhi_)).second++;
0237               totalNumberOfElementsInTheContainer++;
0238             } else {
0239               container_.at(lookupMap_.at(index(ieta, iphi % nPhi_)).first) = *id_itr;
0240               lookupMap_.at(index(ieta, iphi % nPhi_)).first--;
0241               totalNumberOfElementsInTheContainer--;
0242             }
0243       }
0244     }
0245     LogTrace("TrackAssociator") << "Number of elements outside the allowed range ( |eta|>" << nEta_ / 2 * etaBinSize_
0246                                 << "): " << numberOfDetIdsOutsideEtaRange << "\n";
0247     LogTrace("TrackAssociator") << "Number of active DetId's mapped: " << numberOfDetIdsActive << "\n";
0248     if (step == 0) {
0249       // allocate
0250       container_.resize(totalNumberOfElementsInTheContainer);
0251       // fill the range index in the lookup map to point to the last element in the range
0252       unsigned int index(0);
0253       for (std::vector<std::pair<unsigned int, unsigned int> >::iterator bin = lookupMap_.begin();
0254            bin != lookupMap_.end();
0255            ++bin) {
0256         if (bin->second == 0)
0257           continue;
0258         index += bin->second;
0259         bin->first = index - 1;
0260       }
0261     }
0262   }
0263   if (totalNumberOfElementsInTheContainer != 0)
0264     throw cms::Exception("FatalError")
0265         << "Look-up map filled incorrectly. Structural problem. Get in touch with the developer.";
0266   volume_.determinInnerDimensions();
0267   edm::LogVerbatim("TrackAssociator") << "Fiducial volume for " << name()
0268                                       << " (minR, maxR, minZ, maxZ): " << volume_.minR() << ", " << volume_.maxR()
0269                                       << ", " << volume_.minZ() << ", " << volume_.maxZ();
0270   theMapIsValid_ = true;
0271 }
0272 
0273 std::set<DetId> DetIdAssociator::getDetIdsInACone(const std::set<DetId>& inset,
0274                                                   const std::vector<GlobalPoint>& trajectory,
0275                                                   const double dR) const {
0276   if (selectAllInACone(dR))
0277     return inset;
0278   check_setup();
0279   std::set<DetId> outset;
0280   for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
0281     for (std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end();
0282          point_iter++)
0283       if (nearElement(*point_iter, *id_iter, dR)) {
0284         outset.insert(*id_iter);
0285         break;
0286       }
0287   return outset;
0288 }
0289 
0290 std::vector<DetId> DetIdAssociator::getCrossedDetIds(const std::set<DetId>& inset,
0291                                                      const std::vector<GlobalPoint>& trajectory) const {
0292   check_setup();
0293   std::vector<DetId> output;
0294   std::set<DetId> ids(inset);
0295   for (unsigned int i = 0; i + 1 < trajectory.size(); ++i) {
0296     std::set<DetId>::const_iterator id_iter = ids.begin();
0297     while (id_iter != ids.end()) {
0298       if (crossedElement(trajectory[i], trajectory[i + 1], *id_iter)) {
0299         output.push_back(*id_iter);
0300         ids.erase(id_iter++);
0301       } else
0302         id_iter++;
0303     }
0304   }
0305   return output;
0306 }
0307 
0308 std::vector<DetId> DetIdAssociator::getCrossedDetIds(const std::set<DetId>& inset,
0309                                                      const std::vector<SteppingHelixStateInfo>& trajectory,
0310                                                      const double tolerance) const {
0311   check_setup();
0312   std::vector<DetId> output;
0313   std::set<DetId> ids(inset);
0314   for (unsigned int i = 0; i + 1 < trajectory.size(); ++i) {
0315     std::set<DetId>::const_iterator id_iter = ids.begin();
0316     while (id_iter != ids.end()) {
0317       if (crossedElement(trajectory[i].position(), trajectory[i + 1].position(), *id_iter, tolerance, &trajectory[i])) {
0318         output.push_back(*id_iter);
0319         ids.erase(id_iter++);
0320       } else
0321         id_iter++;
0322     }
0323   }
0324   return output;
0325 }
0326 
0327 void DetIdAssociator::dumpMapContent(int ieta, int iphi) const {
0328   if (!(ieta >= 0 && ieta < nEta_ && iphi >= 0)) {
0329     edm::LogWarning("TrackAssociator") << "ieta or iphi is out of range. Skipped.";
0330     return;
0331   }
0332 
0333   std::vector<GlobalPoint> pointBuffer;
0334   std::set<DetId> set;
0335   fillSet(set, ieta, iphi % nPhi_);
0336   LogTrace("TrackAssociator") << "Map content for cell (ieta,iphi): " << ieta << ", " << iphi % nPhi_;
0337   for (std::set<DetId>::const_iterator itr = set.begin(); itr != set.end(); itr++) {
0338     LogTrace("TrackAssociator") << "\tDetId " << itr->rawId() << ", geometry (x,y,z,rho,eta,phi):";
0339     std::pair<const_iterator, const_iterator> points = getDetIdPoints(*itr, pointBuffer);
0340     for (std::vector<GlobalPoint>::const_iterator point = points.first; point != points.second; point++)
0341       LogTrace("TrackAssociator") << "\t\t" << point->x() << ", " << point->y() << ", " << point->z() << ", "
0342                                   << point->perp() << ", " << point->eta() << ", " << point->phi();
0343   }
0344 }
0345 
0346 void DetIdAssociator::dumpMapContent(int ieta_min, int ieta_max, int iphi_min, int iphi_max) const {
0347   for (int i = ieta_min; i <= ieta_max; i++)
0348     for (int j = iphi_min; j <= iphi_max; j++)
0349       dumpMapContent(i, j);
0350 }
0351 
0352 const FiducialVolume& DetIdAssociator::volume() const {
0353   if (!theMapIsValid_)
0354     throw cms::Exception("FatalError") << "map is not valid.";
0355   return volume_;
0356 }
0357 
0358 std::set<DetId> DetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& direction, const MapRange& mapRange) const {
0359   return getDetIdsCloseToAPoint(
0360       direction, mapRange.dThetaPlus, mapRange.dThetaMinus, mapRange.dPhiPlus, mapRange.dPhiMinus);
0361 }
0362 
0363 bool DetIdAssociator::nearElement(const GlobalPoint& point, const DetId& id, const double distance) const {
0364   GlobalPoint center = getPosition(id);
0365   double deltaPhi(fabs(point.phi() - center.phi()));
0366   if (deltaPhi > M_PI)
0367     deltaPhi = fabs(deltaPhi - M_PI * 2.);
0368   return (point.eta() - center.eta()) * (point.eta() - center.eta()) + deltaPhi * deltaPhi < distance * distance;
0369 }
0370 
0371 void DetIdAssociator::check_setup() const {
0372   if (nEta_ == 0)
0373     throw cms::Exception("FatalError") << "Number of eta bins is not set.\n";
0374   if (nPhi_ == 0)
0375     throw cms::Exception("FatalError") << "Number of phi bins is not set.\n";
0376   // if (ivProp_==0) throw cms::Exception("FatalError") << "Track propagator is not defined\n";
0377   if (etaBinSize_ == 0)
0378     throw cms::Exception("FatalError") << "Eta bin size is not set.\n";
0379 }
0380 
0381 void DetIdAssociator::fillSet(std::set<DetId>& set, unsigned int iEta, unsigned int iPhi) const {
0382   unsigned int i = index(iEta, iPhi);
0383   unsigned int i0 = lookupMap_.at(i).first;
0384   unsigned int size = lookupMap_.at(i).second;
0385   for (i = i0; i < i0 + size; ++i)
0386     set.insert(container_.at(i));
0387 }
0388 
0389 #include "FWCore/PluginManager/interface/ModuleDef.h"
0390 #include "FWCore/Framework/interface/MakerMacros.h"
0391 
0392 #include "FWCore/Utilities/interface/typelookup.h"
0393 TYPELOOKUP_DATA_REG(DetIdAssociator);
0394 
0395 //  LocalWords:  Fiducial
[ 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!