SimpleConeBounds.h

Go to the documentation of this file.

#ifndef Geom_SimpleConeBounds_H
#define Geom_SimpleConeBounds_H

#include "DataFormats/GeometryVector/interface/LocalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/Pi.h"
#include "DataFormats/GeometrySurface/interface/LocalError.h"
#include "DataFormats/GeometrySurface/interface/Bounds.h"

#include <cmath>
#include <limits>
#include <algorithm>


class SimpleConeBounds final : public Bounds {
public:

  SimpleConeBounds( float zmin, float rmin_zmin, float rmax_zmin,
            float zmax, float rmin_zmax, float rmax_zmax) : 
    theZmin(zmin), theRminZmin(rmin_zmin), theRmaxZmin(rmax_zmin),
    theZmax(zmax), theRminZmax(rmin_zmax), theRmaxZmax(rmax_zmax) {
    if (theZmin > theZmax) {
      std::swap(theRminZmin, theRminZmax);
      std::swap(theRmaxZmin, theRmaxZmax);
    }
    if (theRminZmin > theRmaxZmin) std::swap(theRminZmin, theRmaxZmin);    
    if (theRminZmax > theRmaxZmax) std::swap(theRminZmax, theRmaxZmax);
  }

  float length()    const override { return theZmax - theZmin;}
  float width()     const override { return 2*std::max(theRmaxZmin,theRmaxZmax);}
  float thickness() const override { return ((theRmaxZmin-theRminZmin)+
                        (theRmaxZmax-theRminZmax))/2.;}
  
  using Bounds::inside;

  bool inside( const Local3DPoint& p) const override {
    float lrmin = (p.z()-theZmin)*(theRminZmax-theRminZmin)/(theZmax-theZmin);
    float lrmax = (p.z()-theZmin)*(theRmaxZmax-theRmaxZmin)/(theZmax-theZmin);
    return p.z()    > theZmin && p.z()    < theZmax &&
           p.perp() > lrmin && p.perp() < lrmax;
  }
    
  bool inside( const Local3DPoint& p, const LocalError& err, float scale) const override {
    // std::cout << "WARNING: SimpleConeBounds::inside(const Local3DPoint&, const LocalError not fully implemented" 
    //        << std::endl;     // FIXME! does not check R.
    SimpleConeBounds tmp(theZmin - sqrt(err.yy())*scale,
             theRminZmin, theRmaxZmin,
             theZmax+ sqrt(err.yy())*scale,
             theRminZmax, theRmaxZmax);
    return tmp.inside(p);
  }

  virtual bool inside( const Local2DPoint& p, const LocalError& err) const {
    return Bounds::inside(p,err);
  }

  Bounds* clone() const override { 
    return new SimpleConeBounds(*this);
  }

  // Extension of the Bounds interface
  Geom::Theta<float> openingAngle() const {
    float theta = atan(((theRmaxZmax+theRminZmax)/2.-
            (theRmaxZmin+theRminZmin)/2.)/length());
    return Geom::Theta<float>(theta < 0 ? theta+Geom::pi() : theta);
  }

  GlobalPoint vertex() const{
    float rAtZmax = (theRmaxZmax+theRminZmax)/2.;
    float rAtZmin = (theRmaxZmin+theRminZmin)/2.;
    float dr = (rAtZmax-rAtZmin);

    if (std::abs(dr) < 0.0001) { // handle degenerate case (cone->cylinder)
      return GlobalPoint(0,0,std::numeric_limits<float>::max());
    } else {
      return GlobalPoint(0,0,(theZmin*rAtZmax - theZmax*rAtZmin)/dr);
    }
  }

private:
  float theZmin;
  float theRminZmin;
  float theRmaxZmin;
  float theZmax;
  float theRminZmax;
  float theRmaxZmax;
};

#endif // Geom_SimpleConeBounds_H