TkDetLayers/src/BladeShapeBuilderFromDet.cc

0001 #include "BladeShapeBuilderFromDet.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0005 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
0006
0007 #include <iomanip>
0008
0009 using namespace std;
0010
0011 namespace {
0012
0013   pair<DiskSectorBounds*, GlobalVector> computeBounds(const vector<const GeomDet*>& dets, const Plane& plane) {
0014     Surface::PositionType tmpPos = dets.front()->surface().position();
0015
0016     float rmin(plane.toLocal(tmpPos).perp());
0017     float rmax(plane.toLocal(tmpPos).perp());
0018     float zmin(plane.toLocal(tmpPos).z());
0019     float zmax(plane.toLocal(tmpPos).z());
0020     float phimin(plane.toLocal(tmpPos).phi());
0021     float phimax(plane.toLocal(tmpPos).phi());
0022
0023     for (vector<const GeomDet*>::const_iterator it = dets.begin(); it != dets.end(); it++) {
0024       vector<GlobalPoint> corners = BoundingBox().corners((*it)->specificSurface());
0025
0026       for (vector<GlobalPoint>::const_iterator i = corners.begin(); i != corners.end(); i++) {
0027         float r = plane.toLocal(*i).perp();
0028         float z = plane.toLocal(*i).z();
0029         float phi = plane.toLocal(*i).phi();
0030         rmin = min(rmin, r);
0031         rmax = max(rmax, r);
0032         zmin = min(zmin, z);
0033         zmax = max(zmax, z);
0034         if (Geom::phiLess(phi, phimin))
0035           phimin = phi;
0036         if (Geom::phiLess(phimax, phi))
0037           phimax = phi;
0038       }
0039       // in addition to the corners we have to check the middle of the
0040       // det +/- length/2, since the min (max) radius for typical fw
0041       // dets is reached there
0042
0043       float rdet = (*it)->position().perp();
0044       float height = (*it)->surface().bounds().width();
0045       rmin = min(rmin, rdet - height / 2.F);
0046       rmax = max(rmax, rdet + height / 2.F);
0047     }
0048
0049     if (!Geom::phiLess(phimin, phimax))
0050       edm::LogError("TkDetLayers") << " BladeShapeBuilderFromDet : "
0051                                    << "Something went wrong with Phi Sorting !";
0052     float zPos = (zmax + zmin) / 2.;
0053     float phiWin = phimax - phimin;
0054     float phiPos = (phimax + phimin) / 2.;
0055     float rmed = (rmin + rmax) / 2.;
0056     if (phiWin < 0.) {
0057       if ((phimin < Geom::pi() / 2.) || (phimax > -Geom::pi() / 2.)) {
0058         edm::LogError("TkDetLayers") << " something strange going on, please check " << phimin << " " << phimax << " "
0059                                      << phiWin;
0060       }
0061       //edm::LogInfo(TkDetLayers) << " Wedge at pi: phi " << phimin << " " << phimax << " " << phiWin
0062       //     << " " << 2.*Geom::pi()+phiWin << " " ;
0063       phiWin += 2. * Geom::pi();
0064       phiPos += Geom::pi();
0065     }
0066
0067     LocalVector localPos(rmed * cos(phiPos), rmed * sin(phiPos), zPos);
0068
0069     LogDebug("TkDetLayers") << "localPos in computeBounds: " << localPos << "\n"
0070                             << "rmin:   " << rmin << "\n"
0071                             << "rmax:   " << rmax << "\n"
0072                             << "zmin:   " << zmin << "\n"
0073                             << "zmax:   " << zmax << "\n"
0074                             << "phiWin: " << phiWin;
0075
0076     return make_pair(new DiskSectorBounds(rmin, rmax, zmin, zmax, phiWin), plane.toGlobal(localPos));
0077   }
0078
0079   Surface::RotationType computeRotation(const vector<const GeomDet*>& dets, const Surface::PositionType& meanPos) {
0080     const Plane& plane = dets.front()->surface();
0081
0082     GlobalVector xAxis;
0083     GlobalVector yAxis;
0084     GlobalVector zAxis;
0085
0086     GlobalVector planeXAxis = plane.toGlobal(LocalVector(1, 0, 0));
0087     const GlobalPoint& planePosition = plane.position();
0088
0089     if (planePosition.x() * planeXAxis.x() + planePosition.y() * planeXAxis.y() > 0.) {
0090       yAxis = planeXAxis;
0091     } else {
0092       yAxis = -planeXAxis;
0093     }
0094
0095     GlobalVector planeZAxis = plane.toGlobal(LocalVector(0, 0, 1));
0096     if (planeZAxis.z() * planePosition.z() > 0.) {
0097       zAxis = planeZAxis;
0098     } else {
0099       zAxis = -planeZAxis;
0100     }
0101
0102     xAxis = yAxis.cross(zAxis);
0103
0104     return Surface::RotationType(xAxis, yAxis);
0105   }
0106
0107 }  // namespace
0108
0109 BoundDiskSector* BladeShapeBuilderFromDet::build(const vector<const GeomDet*>& dets) {
0110   // find mean position
0111   typedef Surface::PositionType::BasicVectorType Vector;
0112   Vector posSum(0, 0, 0);
0113   for (vector<const GeomDet*>::const_iterator i = dets.begin(); i != dets.end(); i++) {
0114     posSum += (**i).surface().position().basicVector();
0115   }
0116   Surface::PositionType meanPos(0., 0., posSum.z() / float(dets.size()));
0117
0118   // temporary plane - for the computation of bounds
0119   Surface::RotationType rotation = computeRotation(dets, meanPos);
0120   Plane tmpPlane(meanPos, rotation);
0121
0122   auto bo = computeBounds(dets, tmpPlane);
0123   GlobalPoint pos = meanPos + bo.second;
0124   //edm::LogInfo(TkDetLayers) << "global pos in operator: " << pos ;
0125   return new BoundDiskSector(pos, rotation, bo.first);
0126 }