PhiBorderFinder.cc

Go to the documentation of this file.

#include "TrackingTools/DetLayers/interface/PhiBorderFinder.h"

PhiBorderFinder::PhiBorderFinder(const std::vector<const Det*>& utheDets) 
  : theNbins(utheDets.size()), 
    isPhiPeriodic_(false), 
    isPhiOverlapping_(false) {
  std::vector<const Det*> theDets = utheDets;
  precomputed_value_sort(theDets.begin(), theDets.end(), DetPhi());
  
  const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|PhiBorderFinder";
  
  double step = 2.*Geom::pi()/theNbins;
  
  LogTrace(metname) << "RecoMuonDetLayers::PhiBorderFinder "
            << "step w: " << step << " # of bins: " << theNbins;
  
  std::vector<double> spread(theNbins);
  std::vector<std::pair<double,double> > phiEdge;
  phiEdge.reserve(theNbins);
  thePhiBorders.reserve(theNbins);
  thePhiBins.reserve(theNbins);
  for ( unsigned int i = 0; i < theNbins; i++ ) {
    thePhiBins.push_back(theDets[i]->position().phi());
    spread.push_back(theDets[i]->position().phi()
             - (theDets[0]->position().phi() + i*step));
    
    LogTrace(metname) << "bin: " << i << " phi bin: " << thePhiBins[i]
              << " spread: " <<  spread[i];
    
    
    ConstReferenceCountingPointer<BoundPlane> plane = 
      dynamic_cast<const BoundPlane*>(&theDets[i]->surface());
    if (plane==nullptr) {
      //FIXME
      throw cms::Exception("UnexpectedState") << ("PhiBorderFinder: det surface is not a BoundPlane");
    }
    
    std::vector<GlobalPoint> dc = 
      BoundingBox().corners(*plane);
    
    float phimin(999.), phimax(-999.);
    for (std::vector<GlobalPoint>::const_iterator pt=dc.begin();
     pt!=dc.end(); pt++) {
      float phi = (*pt).phi();
      //    float z = pt->z();
      if (phi < phimin) phimin = phi;
      if (phi > phimax) phimax = phi;
    }
    if (phimin*phimax < 0. &&           //Handle pi border:
    phimax - phimin > Geom::pi()) { //Assume that the Det is on
      //the shortest side 
      std::swap(phimin,phimax);
    }
    phiEdge.push_back(std::pair<double,double>(phimin,phimax));
  }
  
  for (unsigned int i = 0; i < theNbins; i++) {
    
    // Put the two phi values in the [0,2pi] range
    double firstPhi  = positiveRange(phiEdge[i].first);
    double secondPhi = positiveRange(phiEdge[binIndex(i-1)].second);
    
    // Reformat them in the [-pi,pi] range
    Geom::Phi<double> firstEdge(firstPhi);
    Geom::Phi<double> secondEdge(secondPhi);
    
    // Calculate the mean and format the result in the [-pi,pi] range
    // Get their value in order to perform the mean in the correct way
    Geom::Phi<double> mean((firstEdge.value() + secondEdge.value())/2.);
    
    // Special case: at +-pi there is a discontinuity
    if ( phiEdge[i].first * phiEdge[binIndex(i-1)].second < 0 &&
     fabs(firstPhi-secondPhi) < Geom::pi() ) 
      mean = Geom::pi() - mean;
    
    thePhiBorders.push_back(mean);
  }
  
  for (unsigned int i = 0; i < theNbins; i++) {
    if (Geom::Phi<double>(phiEdge[i].first)
    - Geom::Phi<double>(phiEdge[binIndex(i-1)].second) < 0) {
      isPhiOverlapping_ = true;
      break;
    }
  }
  
  double rms = stat_RMS(spread); 
  if ( rms < 0.01*step) { 
    isPhiPeriodic_ = true;
  }
  
  //Check that everything is proper
  if (thePhiBorders.size() != theNbins || thePhiBins.size() != theNbins) 
    //FIXME
    throw cms::Exception("UnexpectedState") << "PhiBorderFinder: consistency error";
}