|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
#include "Geometry/TrackerGeometryBuilder/interface/PlaneBuilderForGluedDet.h"
#include "DataFormats/GeometrySurface/interface/Surface.h"
#include "DataFormats/GeometrySurface/interface/BoundingBox.h"
#include "DataFormats/GeometrySurface/interface/MediumProperties.h"
#include "DataFormats/GeometrySurface/interface/OpenBounds.h"
#include <algorithm>
// Warning, remember to assign this pointer to a ReferenceCountingPointer!
PlaneBuilderForGluedDet::ResultType PlaneBuilderForGluedDet::plane(const std::vector<const GeomDetUnit*>& dets) const {
// find mean position
typedef Surface::PositionType::BasicVectorType Vector;
Vector posSum(0, 0, 0);
for (auto det : dets) {
posSum += (*det).surface().position().basicVector();
}
Surface::PositionType meanPos(posSum / float(dets.size()));
Surface::RotationType rotation = dets.front()->surface().rotation();
// Surface::RotationType rotation = computeRotation( dets, meanPos);
Plane tmpPlane = Plane(meanPos, rotation);
// Take the medium properties from the first DetUnit
const MediumProperties& mp = dets.front()->surface().mediumProperties();
MediumProperties newmp(mp.radLen() * 2.0, mp.xi() * 2.0);
std::pair<RectangularPlaneBounds*, GlobalVector> bo = computeRectBounds(dets, tmpPlane);
return new Plane(meanPos + bo.second, rotation, newmp, bo.first);
}
std::pair<RectangularPlaneBounds*, GlobalVector> PlaneBuilderForGluedDet::computeRectBounds(
const std::vector<const GeomDetUnit*>& dets, const Plane& plane) const {
// go over all corners and compute maximum deviations from mean pos.
std::vector<GlobalPoint> corners;
for (auto det : dets) {
const Plane& bplane = dynamic_cast<const Plane&>(det->surface());
std::vector<GlobalPoint> dc = BoundingBox().corners(bplane);
corners.insert(corners.end(), dc.begin(), dc.end());
}
float xmin(0), xmax(0), ymin(0), ymax(0), zmin(0), zmax(0);
for (std::vector<GlobalPoint>::const_iterator i = corners.begin(), cend = corners.end(); i != cend; ++i) {
LocalPoint p = plane.toLocal(*i);
if (p.x() < xmin)
xmin = p.x();
if (p.x() > xmax)
xmax = p.x();
if (p.y() < ymin)
ymin = p.y();
if (p.y() > ymax)
ymax = p.y();
if (p.z() < zmin)
zmin = p.z();
if (p.z() > zmax)
zmax = p.z();
}
LocalVector localOffset((xmin + xmax) / 2., (ymin + ymax) / 2., (zmin + zmax) / 2.);
GlobalVector offset(plane.toGlobal(localOffset));
std::pair<RectangularPlaneBounds*, GlobalVector> result(
new RectangularPlaneBounds((xmax - xmin) / 2, (ymax - ymin) / 2, (zmax - zmin) / 2), offset);
return result;
}
Surface::RotationType PlaneBuilderForGluedDet::computeRotation(const std::vector<GeomDetUnit const*>& dets,
const Surface::PositionType& meanPos) const {
// choose first mono out-pointing rotation
// the rotations of GluedDets coincide with the mono part
// Simply take the x,y of the first Det if z points out,
// or -x, y if it doesn't
const BoundPlane& plane = dynamic_cast<const BoundPlane&>(dets.front()->surface());
//GlobalVector n = plane.normalVector();
GlobalVector xAxis;
GlobalVector yAxis;
GlobalVector planeYAxis = plane.toGlobal(LocalVector(0, 1, 0));
if (planeYAxis.z() < 0)
yAxis = -planeYAxis;
else
yAxis = planeYAxis;
GlobalVector planeXAxis = plane.toGlobal(LocalVector(1, 0, 0));
GlobalVector n = planeXAxis.cross(planeYAxis);
if (n.x() * meanPos.x() + n.y() * meanPos.y() > 0) {
xAxis = planeXAxis;
} else {
xAxis = -planeXAxis;
}
return Surface::RotationType(xAxis, yAxis);
}
|