TrackerGeometryBuilder/src/PlaneBuilderForGluedDet.cc

0001 #include "Geometry/TrackerGeometryBuilder/interface/PlaneBuilderForGluedDet.h"
0002 #include "DataFormats/GeometrySurface/interface/Surface.h"
0003 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
0004 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0005 #include "DataFormats/GeometrySurface/interface/OpenBounds.h"
0006
0007 #include <algorithm>
0008
0009 // Warning, remember to assign this pointer to a ReferenceCountingPointer!
0010 PlaneBuilderForGluedDet::ResultType PlaneBuilderForGluedDet::plane(const std::vector<const GeomDetUnit*>& dets) const {
0011   // find mean position
0012   typedef Surface::PositionType::BasicVectorType Vector;
0013   Vector posSum(0, 0, 0);
0014   for (auto det : dets) {
0015     posSum += (*det).surface().position().basicVector();
0016   }
0017   Surface::PositionType meanPos(posSum / float(dets.size()));
0018
0019   Surface::RotationType rotation = dets.front()->surface().rotation();
0020   //  Surface::RotationType rotation = computeRotation( dets, meanPos);
0021   Plane tmpPlane = Plane(meanPos, rotation);
0022
0023   // Take the medium properties from the first DetUnit
0024   const MediumProperties& mp = dets.front()->surface().mediumProperties();
0025   MediumProperties newmp(mp.radLen() * 2.0, mp.xi() * 2.0);
0026
0027   std::pair<RectangularPlaneBounds*, GlobalVector> bo = computeRectBounds(dets, tmpPlane);
0028   return new Plane(meanPos + bo.second, rotation, newmp, bo.first);
0029 }
0030
0031 std::pair<RectangularPlaneBounds*, GlobalVector> PlaneBuilderForGluedDet::computeRectBounds(
0032     const std::vector<const GeomDetUnit*>& dets, const Plane& plane) const {
0033   // go over all corners and compute maximum deviations from mean pos.
0034   std::vector<GlobalPoint> corners;
0035   for (auto det : dets) {
0036     const Plane& bplane = dynamic_cast<const Plane&>(det->surface());
0037     std::vector<GlobalPoint> dc = BoundingBox().corners(bplane);
0038     corners.insert(corners.end(), dc.begin(), dc.end());
0039   }
0040
0041   float xmin(0), xmax(0), ymin(0), ymax(0), zmin(0), zmax(0);
0042   for (std::vector<GlobalPoint>::const_iterator i = corners.begin(), cend = corners.end(); i != cend; ++i) {
0043     LocalPoint p = plane.toLocal(*i);
0044     if (p.x() < xmin)
0045       xmin = p.x();
0046     if (p.x() > xmax)
0047       xmax = p.x();
0048     if (p.y() < ymin)
0049       ymin = p.y();
0050     if (p.y() > ymax)
0051       ymax = p.y();
0052     if (p.z() < zmin)
0053       zmin = p.z();
0054     if (p.z() > zmax)
0055       zmax = p.z();
0056   }
0057
0058   LocalVector localOffset((xmin + xmax) / 2., (ymin + ymax) / 2., (zmin + zmax) / 2.);
0059   GlobalVector offset(plane.toGlobal(localOffset));
0060
0061   std::pair<RectangularPlaneBounds*, GlobalVector> result(
0062       new RectangularPlaneBounds((xmax - xmin) / 2, (ymax - ymin) / 2, (zmax - zmin) / 2), offset);
0063
0064   return result;
0065 }
0066
0067 Surface::RotationType PlaneBuilderForGluedDet::computeRotation(const std::vector<GeomDetUnit const*>& dets,
0068                                                                const Surface::PositionType& meanPos) const {
0069   // choose first mono out-pointing rotation
0070   // the rotations of GluedDets coincide with the mono part
0071   // Simply take the x,y of the first Det if z points out,
0072   // or -x, y if it doesn't
0073   const BoundPlane& plane = dynamic_cast<const BoundPlane&>(dets.front()->surface());
0074   //GlobalVector n = plane.normalVector();
0075
0076   GlobalVector xAxis;
0077   GlobalVector yAxis;
0078   GlobalVector planeYAxis = plane.toGlobal(LocalVector(0, 1, 0));
0079   if (planeYAxis.z() < 0)
0080     yAxis = -planeYAxis;
0081   else
0082     yAxis = planeYAxis;
0083
0084   GlobalVector planeXAxis = plane.toGlobal(LocalVector(1, 0, 0));
0085   GlobalVector n = planeXAxis.cross(planeYAxis);
0086
0087   if (n.x() * meanPos.x() + n.y() * meanPos.y() > 0) {
0088     xAxis = planeXAxis;
0089   } else {
0090     xAxis = -planeXAxis;
0091   }
0092
0093   return Surface::RotationType(xAxis, yAxis);
0094 }