DetLayers/src/MTDDiskSectorBuilderFromDet.cc

0001 //#define EDM_ML_DEBUG
0002
0003 #include "MTDDiskSectorBuilderFromDet.h"
0004
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0007 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
0008
0009 #include <iomanip>
0010
0011 using namespace std;
0012
0013 namespace {
0014
0015   pair<DiskSectorBounds*, GlobalVector> computeBounds(const vector<const GeomDet*>& dets) {
0016     // go over all corners and compute maximum deviations
0017     float rmin(dets.front()->surface().position().perp());
0018     float rmax(rmin);
0019     float zmin(dets.front()->surface().position().z());
0020     float zmax(zmin);
0021     float phimin(dets.front()->surface().position().phi());
0022     float phimax(phimin);
0023
0024     for (auto const& idet : dets) {
0025       vector<GlobalPoint> corners = BoundingBox().corners(idet->specificSurface());
0026       for (auto const& i : corners) {
0027         float r = i.perp();
0028         float z = i.z();
0029         float phi = 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     }
0040
0041     if (!Geom::phiLess(phimin, phimax))
0042       edm::LogError("MTDDetLayers") << " MTDDiskSectorBuilderFromDet : "
0043                                     << "Something went wrong with Phi Sorting !";
0044     float zPos = (zmax + zmin) / 2.;
0045     float phiWin = phimax - phimin;
0046     float phiPos = (phimax + phimin) / 2.;
0047     float rmed = (rmin + rmax) / 2.;
0048     if (phiWin < 0.) {
0049       if ((phimin < Geom::pi() / 2.) || (phimax > -Geom::pi() / 2.)) {
0050         edm::LogError("MTDDetLayers") << " something strange going on, please check " << phimin << " " << phimax << " "
0051                                       << phiWin;
0052       }
0053       phiWin += 2. * Geom::pi();
0054       phiPos += Geom::pi();
0055     }
0056
0057     GlobalVector pos(rmed * cos(phiPos), rmed * sin(phiPos), zPos);
0058
0059     LogTrace("MTDDetLayers") << "MTDDiskSectorBuilderFromDet::computeBounds sector at: " << std::fixed << pos << "\n"
0060                              << "zmin    : " << std::setw(14) << zmin << "\n"
0061                              << "zmax    : " << std::setw(14) << zmax << "\n"
0062                              << "rmin    : " << std::setw(14) << rmin << "\n"
0063                              << "rmax    : " << std::setw(14) << rmax << "\n"
0064                              << "phi ref : " << std::setw(14) << phiPos << "\n"
0065                              << "phi win : " << std::setw(14) << phiWin;
0066
0067     return make_pair(new DiskSectorBounds(rmin, rmax, zmin - zPos, zmax - zPos, phiWin), pos);
0068   }
0069
0070   Surface::RotationType computeRotation(const vector<const GeomDet*>& dets, const Surface::PositionType pos) {
0071     GlobalVector yAxis = (GlobalVector(pos.x(), pos.y(), 0.)).unit();
0072
0073     GlobalVector zAxis(0., 0., 1.);
0074     GlobalVector xAxis = yAxis.cross(zAxis);
0075
0076     return Surface::RotationType(xAxis, yAxis);
0077   }
0078
0079 }  // namespace
0080
0081 BoundDiskSector* MTDDiskSectorBuilderFromDet::operator()(const vector<const GeomDet*>& dets) const {
0082   // check that the dets are all at about the same z
0083   float zcheck = dets.front()->surface().position().z();
0084   constexpr double tol(0.5);  // minimal safety check on z position of modules within a sector, width ~ 10 mm
0085   for (auto const& idet : dets) {
0086     float zdiff = zcheck - (*idet).surface().position().z();
0087     if (std::abs(zdiff) > tol) {
0088       edm::LogError("MTDDetLayers")
0089           << " MTDDiskSectorBuilderFromDet: Trying to build sector from Dets at different z positions !! Delta_z = "
0090           << zdiff;
0091     }
0092   }
0093
0094   auto bo = computeBounds(dets);
0095
0096   Surface::PositionType pos(bo.second.x(), bo.second.y(), bo.second.z());
0097   Surface::RotationType rot = computeRotation(dets, pos);
0098   return new BoundDiskSector(pos, rot, bo.first);
0099 }