DetLayers/src/PhiBorderFinder.cc

0001 #include "TrackingTools/DetLayers/interface/PhiBorderFinder.h"
0002
0003 PhiBorderFinder::PhiBorderFinder(const std::vector<const Det*>& utheDets)
0004     : theNbins(utheDets.size()), isPhiPeriodic_(false), isPhiOverlapping_(false) {
0005   std::vector<const Det*> theDets = utheDets;
0006   precomputed_value_sort(theDets.begin(), theDets.end(), DetPhi());
0007
0008   const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|PhiBorderFinder";
0009
0010   double step = 2. * Geom::pi() / theNbins;
0011
0012   LogTrace(metname) << "RecoMuonDetLayers::PhiBorderFinder "
0013                     << "step w: " << step << " # of bins: " << theNbins;
0014
0015   std::vector<double> spread(theNbins);
0016   std::vector<std::pair<double, double> > phiEdge;
0017   phiEdge.reserve(theNbins);
0018   thePhiBorders.reserve(theNbins);
0019   thePhiBins.reserve(theNbins);
0020   for (unsigned int i = 0; i < theNbins; i++) {
0021     thePhiBins.push_back(theDets[i]->position().phi());
0022     spread.push_back(theDets[i]->position().phi() - (theDets[0]->position().phi() + i * step));
0023
0024     LogTrace(metname) << "bin: " << i << " phi bin: " << thePhiBins[i] << " spread: " << spread[i];
0025
0026     ConstReferenceCountingPointer<BoundPlane> plane = dynamic_cast<const BoundPlane*>(&theDets[i]->surface());
0027     if (plane == nullptr) {
0028       //FIXME
0029       throw cms::Exception("UnexpectedState") << ("PhiBorderFinder: det surface is not a BoundPlane");
0030     }
0031
0032     std::vector<GlobalPoint> dc = BoundingBox().corners(*plane);
0033
0034     float phimin(999.), phimax(-999.);
0035     for (std::vector<GlobalPoint>::const_iterator pt = dc.begin(); pt != dc.end(); pt++) {
0036       float phi = (*pt).phi();
0037       //    float z = pt->z();
0038       if (phi < phimin)
0039         phimin = phi;
0040       if (phi > phimax)
0041         phimax = phi;
0042     }
0043     if (phimin * phimax < 0. &&          //Handle pi border:
0044         phimax - phimin > Geom::pi()) {  //Assume that the Det is on
0045       //the shortest side
0046       std::swap(phimin, phimax);
0047     }
0048     phiEdge.push_back(std::pair<double, double>(phimin, phimax));
0049   }
0050
0051   for (unsigned int i = 0; i < theNbins; i++) {
0052     // Put the two phi values in the [0,2pi] range
0053     double firstPhi = positiveRange(phiEdge[i].first);
0054     double secondPhi = positiveRange(phiEdge[binIndex(i - 1)].second);
0055
0056     // Reformat them in the [-pi,pi] range
0057     Geom::Phi<double> firstEdge(firstPhi);
0058     Geom::Phi<double> secondEdge(secondPhi);
0059
0060     // Calculate the mean and format the result in the [-pi,pi] range
0061     // Get their value in order to perform the mean in the correct way
0062     Geom::Phi<double> mean((firstEdge.value() + secondEdge.value()) / 2.);
0063
0064     // Special case: at +-pi there is a discontinuity
0065     if (phiEdge[i].first * phiEdge[binIndex(i - 1)].second < 0 && fabs(firstPhi - secondPhi) < Geom::pi())
0066       mean = Geom::pi() - mean;
0067
0068     thePhiBorders.push_back(mean);
0069   }
0070
0071   for (unsigned int i = 0; i < theNbins; i++) {
0072     if (Geom::Phi<double>(phiEdge[i].first) - Geom::Phi<double>(phiEdge[binIndex(i - 1)].second) < 0) {
0073       isPhiOverlapping_ = true;
0074       break;
0075     }
0076   }
0077
0078   double rms = stat_RMS(spread);
0079   if (rms < 0.01 * step) {
0080     isPhiPeriodic_ = true;
0081   }
0082
0083   //Check that everything is proper
0084   if (thePhiBorders.size() != theNbins || thePhiBins.size() != theNbins)
0085     //FIXME
0086     throw cms::Exception("UnexpectedState") << "PhiBorderFinder: consistency error";
0087 }