CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
BladeShapeBuilderFromDet.cc
Go to the documentation of this file.
2 
6 
7 #include <iomanip>
8 
9 using namespace std;
10 
11 namespace {
12 
13  pair<DiskSectorBounds*, GlobalVector> computeBounds(const vector<const GeomDet*>& dets, const Plane& plane) {
14  Surface::PositionType tmpPos = dets.front()->surface().position();
15 
16  float rmin(plane.toLocal(tmpPos).perp());
17  float rmax(plane.toLocal(tmpPos).perp());
18  float zmin(plane.toLocal(tmpPos).z());
19  float zmax(plane.toLocal(tmpPos).z());
20  float phimin(plane.toLocal(tmpPos).phi());
21  float phimax(plane.toLocal(tmpPos).phi());
22 
23  for (vector<const GeomDet*>::const_iterator it = dets.begin(); it != dets.end(); it++) {
24  vector<GlobalPoint> corners = BoundingBox().corners((*it)->specificSurface());
25 
26  for (vector<GlobalPoint>::const_iterator i = corners.begin(); i != corners.end(); i++) {
27  float r = plane.toLocal(*i).perp();
28  float z = plane.toLocal(*i).z();
29  float phi = plane.toLocal(*i).phi();
30  rmin = min(rmin, r);
31  rmax = max(rmax, r);
32  zmin = min(zmin, z);
33  zmax = max(zmax, z);
34  if (Geom::phiLess(phi, phimin))
35  phimin = phi;
36  if (Geom::phiLess(phimax, phi))
37  phimax = phi;
38  }
39  // in addition to the corners we have to check the middle of the
40  // det +/- length/2, since the min (max) radius for typical fw
41  // dets is reached there
42 
43  float rdet = (*it)->position().perp();
44  float height = (*it)->surface().bounds().width();
45  rmin = min(rmin, rdet - height / 2.F);
46  rmax = max(rmax, rdet + height / 2.F);
47  }
48 
50  edm::LogError("TkDetLayers") << " BladeShapeBuilderFromDet : "
51  << "Something went wrong with Phi Sorting !";
52  float zPos = (zmax + zmin) / 2.;
53  float phiWin = phimax - phimin;
54  float phiPos = (phimax + phimin) / 2.;
55  float rmed = (rmin + rmax) / 2.;
56  if (phiWin < 0.) {
57  if ((phimin < Geom::pi() / 2.) || (phimax > -Geom::pi() / 2.)) {
58  edm::LogError("TkDetLayers") << " something strange going on, please check " << phimin << " " << phimax << " "
59  << phiWin;
60  }
61  //edm::LogInfo(TkDetLayers) << " Wedge at pi: phi " << phimin << " " << phimax << " " << phiWin
62  // << " " << 2.*Geom::pi()+phiWin << " " ;
63  phiWin += 2. * Geom::pi();
64  phiPos += Geom::pi();
65  }
66 
67  LocalVector localPos(rmed * cos(phiPos), rmed * sin(phiPos), zPos);
68 
69  LogDebug("TkDetLayers") << "localPos in computeBounds: " << localPos << "\n"
70  << "rmin: " << rmin << "\n"
71  << "rmax: " << rmax << "\n"
72  << "zmin: " << zmin << "\n"
73  << "zmax: " << zmax << "\n"
74  << "phiWin: " << phiWin;
75 
76  return make_pair(new DiskSectorBounds(rmin, rmax, zmin, zmax, phiWin), plane.toGlobal(localPos));
77  }
78 
79  Surface::RotationType computeRotation(const vector<const GeomDet*>& dets, const Surface::PositionType& meanPos) {
80  const Plane& plane = dets.front()->surface();
81 
82  GlobalVector xAxis;
83  GlobalVector yAxis;
85 
86  GlobalVector planeXAxis = plane.toGlobal(LocalVector(1, 0, 0));
87  const GlobalPoint& planePosition = plane.position();
88 
89  if (planePosition.x() * planeXAxis.x() + planePosition.y() * planeXAxis.y() > 0.) {
90  yAxis = planeXAxis;
91  } else {
92  yAxis = -planeXAxis;
93  }
94 
95  GlobalVector planeZAxis = plane.toGlobal(LocalVector(0, 0, 1));
96  if (planeZAxis.z() * planePosition.z() > 0.) {
97  zAxis = planeZAxis;
98  } else {
99  zAxis = -planeZAxis;
100  }
101 
102  xAxis = yAxis.cross(zAxis);
103 
104  return Surface::RotationType(xAxis, yAxis);
105  }
106 
107 } // namespace
108 
109 BoundDiskSector* BladeShapeBuilderFromDet::build(const vector<const GeomDet*>& dets) {
110  // find mean position
112  Vector posSum(0, 0, 0);
113  for (vector<const GeomDet*>::const_iterator i = dets.begin(); i != dets.end(); i++) {
114  posSum += (**i).surface().position().basicVector();
115  }
116  Surface::PositionType meanPos(0., 0., posSum.z() / float(dets.size()));
117 
118  // temporary plane - for the computation of bounds
119  Surface::RotationType rotation = computeRotation(dets, meanPos);
120  Plane tmpPlane(meanPos, rotation);
121 
122  auto bo = computeBounds(dets, tmpPlane);
123  GlobalPoint pos = meanPos + bo.second;
124  //edm::LogInfo(TkDetLayers) << "global pos in operator: " << pos ;
125  return new BoundDiskSector(pos, rotation, bo.first);
126 }
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:79
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:110
Local3DVector LocalVector
Definition: LocalVector.h:12
T perp() const
Definition: PV3DBase.h:69
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
T y() const
Definition: PV3DBase.h:60
Log< level::Error, false > LogError
Definition: Plane.h:16
float float float z
LocalPoint toLocal(const GlobalPoint &gp) const
T z() const
Definition: PV3DBase.h:61
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
T min(T a, T b)
Definition: MathUtil.h:58
bool phiLess(float phi1, float phi2)
Definition: VectorUtil.h:18
static BoundDiskSector * build(const std::vector< const GeomDet * > &dets) __attribute__((cold))
static std::vector< GlobalPoint > corners(const Plane &)
Definition: BoundingBox.cc:20
TkRotation< float > RotationType
tuple zAxis
Definition: MetAnalyzer.py:57
constexpr double pi()
Definition: Pi.h:31
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:163
T x() const
Definition: PV3DBase.h:59
const PositionType & position() const
#define LogDebug(id)