DetLayers/src/RodPlaneBuilderFromDet.cc

0001 #include "TrackingTools/DetLayers/interface/RodPlaneBuilderFromDet.h"
0002 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0003 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
0004
0005 #include <algorithm>
0006
0007 using namespace std;
0008
0009 // Warning, remember to assign this pointer to a ReferenceCountingPointer!
0010 Plane* RodPlaneBuilderFromDet::operator()(const vector<const Det*>& dets) const {
0011   // find mean position
0012   typedef Surface::PositionType::BasicVectorType Vector;
0013   Vector posSum(0, 0, 0);
0014   for (vector<const Det*>::const_iterator i = dets.begin(); i != dets.end(); i++) {
0015     posSum += (**i).surface().position().basicVector();
0016   }
0017   Surface::PositionType meanPos(posSum / float(dets.size()));
0018
0019   // temporary plane - for the computation of bounds
0020   Surface::RotationType rotation = computeRotation(dets, meanPos);
0021   Plane tmpPlane(meanPos, rotation);
0022   auto bo = computeBounds(dets, tmpPlane);
0023
0024   //   LogDebug("DetLayers") << "Creating plane at position " << meanPos
0025   //        << " displaced by " << bo.second ;
0026   //   LogDebug("DetLayers") << "Bounds are (wid/len/thick) " << bo.first.width()
0027   //        << " / " <<  bo.first.length()
0028   //        << " / " <<  bo.first.thickness() ;
0029
0030   return new Plane(meanPos + bo.second, rotation, bo.first);
0031 }
0032
0033 pair<RectangularPlaneBounds*, GlobalVector> RodPlaneBuilderFromDet::computeBounds(const vector<const Det*>& dets,
0034                                                                                   const Plane& plane) const {
0035   // go over all corners and compute maximum deviations from mean pos.
0036   vector<GlobalPoint> corners;
0037   for (vector<const Det*>::const_iterator idet = dets.begin(); idet != dets.end(); idet++) {
0038     /* ---- original implementation. Is it obsolete?
0039     vector<const DetUnit*> detUnits = (**idet).basicComponents();
0040     for (vector<const DetUnit*>::const_iterator detu=detUnits.begin();
0041      detu!=detUnits.end(); detu++) {
0042       vector<GlobalPoint> dc =
0043     BoundingBox().corners((**detu).specificSurface());
0044       corners.insert( corners.end(), dc.begin(), dc.end());
0045     }
0046     ---- */
0047
0048     // temporary implementation (May be the final one if the GluedDet surface
0049     // really contains both the mono and the stereo surfaces
0050     vector<GlobalPoint> dc = BoundingBox().corners((**idet).specificSurface());
0051     corners.insert(corners.end(), dc.begin(), dc.end());
0052   }
0053
0054   float xmin(0), xmax(0), ymin(0), ymax(0), zmin(0), zmax(0);
0055   for (vector<GlobalPoint>::const_iterator i = corners.begin(); i != corners.end(); i++) {
0056     LocalPoint p = plane.toLocal(*i);
0057     if (p.x() < xmin)
0058       xmin = p.x();
0059     if (p.x() > xmax)
0060       xmax = p.x();
0061     if (p.y() < ymin)
0062       ymin = p.y();
0063     if (p.y() > ymax)
0064       ymax = p.y();
0065     if (p.z() < zmin)
0066       zmin = p.z();
0067     if (p.z() > zmax)
0068       zmax = p.z();
0069   }
0070
0071   LocalVector localOffset((xmin + xmax) / 2., (ymin + ymax) / 2., (zmin + zmax) / 2.);
0072   GlobalVector offset(plane.toGlobal(localOffset));
0073
0074   pair<RectangularPlaneBounds*, GlobalVector> result(
0075       new RectangularPlaneBounds((xmax - xmin) / 2, (ymax - ymin) / 2, (zmax - zmin) / 2), offset);
0076
0077   return result;
0078 }
0079
0080 Surface::RotationType RodPlaneBuilderFromDet::computeRotation(const vector<const Det*>& dets,
0081                                                               const Surface::PositionType& meanPos) const {
0082   // choose first mono out-pointing rotation
0083   // the rotations of GluedDets coincide with the mono part
0084   // Simply take the x,y of the first Det if z points out,
0085   // or -x, y if it doesn't
0086   const Plane& plane = dynamic_cast<const Plane&>(dets.front()->surface());
0087   //GlobalVector n = plane.normalVector();
0088
0089   GlobalVector xAxis;
0090   GlobalVector yAxis;
0091   GlobalVector planeYAxis = plane.toGlobal(LocalVector(0, 1, 0));
0092   if (planeYAxis.z() < 0)
0093     yAxis = -planeYAxis;
0094   else
0095     yAxis = planeYAxis;
0096
0097   GlobalVector planeXAxis = plane.toGlobal(LocalVector(1, 0, 0));
0098   GlobalVector n = planeXAxis.cross(planeYAxis);
0099
0100   if (n.x() * meanPos.x() + n.y() * meanPos.y() > 0) {
0101     xAxis = planeXAxis;
0102   } else {
0103     xAxis = -planeXAxis;
0104   }
0105
0106   //   LogDebug("DetLayers") << "Creating rotation with x,y axis "
0107   //        << xAxis << ", " << yAxis ;
0108
0109   return Surface::RotationType(xAxis, yAxis);
0110 }