00001 #include "TrackingTools/DetLayers/interface/RodPlaneBuilderFromDet.h"
00002 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
00003 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
00004
00005 #include <algorithm>
00006
00007 using namespace std;
00008
00009
00010 BoundPlane*
00011 RodPlaneBuilderFromDet::operator()( const vector<const Det*>& dets) const
00012 {
00013
00014 typedef Surface::PositionType::BasicVectorType Vector;
00015 Vector posSum(0,0,0);
00016 for (vector<const Det*>::const_iterator i=dets.begin(); i!=dets.end(); i++) {
00017 posSum += (**i).surface().position().basicVector();
00018 }
00019 Surface::PositionType meanPos( posSum/float(dets.size()));
00020
00021
00022 Surface::RotationType rotation = computeRotation( dets, meanPos);
00023 BoundPlane tmpPlane( meanPos, rotation);
00024 pair<RectangularPlaneBounds,GlobalVector> bo =
00025 computeBounds( dets, tmpPlane);
00026
00027
00028
00029
00030
00031
00032
00033 return new BoundPlane( meanPos+bo.second, rotation, bo.first);
00034 }
00035
00036 pair<RectangularPlaneBounds, GlobalVector>
00037 RodPlaneBuilderFromDet::computeBounds( const vector<const Det*>& dets,
00038 const BoundPlane& plane) const
00039 {
00040
00041 vector<GlobalPoint> corners;
00042 for (vector<const Det*>::const_iterator idet=dets.begin();
00043 idet != dets.end(); idet++) {
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 vector<GlobalPoint> dc = BoundingBox().corners((**idet).specificSurface());
00058 corners.insert( corners.end(),dc.begin(), dc.end() );
00059 }
00060
00061 float xmin(0), xmax(0), ymin(0), ymax(0), zmin(0), zmax(0);
00062 for (vector<GlobalPoint>::const_iterator i=corners.begin();
00063 i!=corners.end(); i++) {
00064 LocalPoint p = plane.toLocal(*i);
00065 if (p.x() < xmin) xmin = p.x();
00066 if (p.x() > xmax) xmax = p.x();
00067 if (p.y() < ymin) ymin = p.y();
00068 if (p.y() > ymax) ymax = p.y();
00069 if (p.z() < zmin) zmin = p.z();
00070 if (p.z() > zmax) zmax = p.z();
00071 }
00072
00073 LocalVector localOffset( (xmin+xmax)/2., (ymin+ymax)/2., (zmin+zmax)/2.);
00074 GlobalVector offset( plane.toGlobal(localOffset));
00075
00076 pair<RectangularPlaneBounds, GlobalVector> result(RectangularPlaneBounds((xmax-xmin)/2, (ymax-ymin)/2, (zmax-zmin)/2), offset);
00077
00078 return result;
00079 }
00080
00081 Surface::RotationType
00082 RodPlaneBuilderFromDet::
00083 computeRotation( const vector<const Det*>& dets,
00084 const Surface::PositionType& meanPos) const
00085 {
00086
00087
00088
00089
00090 const BoundPlane& plane =
00091 dynamic_cast<const BoundPlane&>(dets.front()->surface());
00092
00093
00094 GlobalVector xAxis;
00095 GlobalVector yAxis;
00096 GlobalVector planeYAxis = plane.toGlobal( LocalVector( 0, 1, 0));
00097 if (planeYAxis.z() < 0) yAxis = -planeYAxis;
00098 else yAxis = planeYAxis;
00099
00100 GlobalVector planeXAxis = plane.toGlobal( LocalVector( 1, 0, 0));
00101 GlobalVector n = planeXAxis.cross( planeYAxis);
00102
00103 if (n.x() * meanPos.x() + n.y() * meanPos.y() > 0) {
00104 xAxis = planeXAxis;
00105 }
00106 else {
00107 xAxis = -planeXAxis;
00108 }
00109
00110
00111
00112
00113 return Surface::RotationType( xAxis, yAxis);
00114 }