TkDetLayers/src/ForwardDiskSectorBuilderFromDet.cc

0001 #include "ForwardDiskSectorBuilderFromDet.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0005 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0006 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0007
0008 using namespace std;
0009
0010 // Warning, remember to assign this pointer to a ReferenceCountingPointer!
0011 BoundDiskSector* ForwardDiskSectorBuilderFromDet::operator()(const vector<const GeomDet*>& dets) const {
0012   // check that the dets are all at about the same radius and z
0013   float rcheck = dets.front()->surface().position().perp();
0014   float zcheck = dets.front()->surface().position().z();
0015   for (vector<const GeomDet*>::const_iterator i = dets.begin(); i != dets.end(); i++) {
0016     float rdiff = (**i).surface().position().perp() - rcheck;
0017     if (std::abs(rdiff) > 1.)
0018       edm::LogError("TkDetLayers") << " ForwardDiskSectorBuilderFromDet: Trying to build Petal Wedge from "
0019                                    << "Dets at different radii !! Delta_r = " << rdiff;
0020     float zdiff = zcheck - (**i).surface().position().z();
0021     if (std::abs(zdiff) > 0.8)
0022       edm::LogError("TkDetLayers") << " ForwardDiskSectorBuilderFromDet: Trying to build Petal Wedge from "
0023                                    << "Dets at different z positions !! Delta_z = " << zdiff;
0024   }
0025
0026   auto bo = computeBounds(dets);
0027
0028   Surface::PositionType pos(bo.second.x(), bo.second.y(), bo.second.z());
0029   Surface::RotationType rot = computeRotation(dets, pos);
0030   return new BoundDiskSector(pos, rot, bo.first);
0031 }
0032
0033 pair<DiskSectorBounds*, GlobalVector> ForwardDiskSectorBuilderFromDet::computeBounds(
0034     const vector<const GeomDet*>& dets) const {
0035   // go over all corners and compute maximum deviations
0036   float rmin = (**(dets.begin())).surface().position().perp();
0037   float rmax(rmin);
0038   float zmin((**(dets.begin())).surface().position().z());
0039   float zmax(zmin);
0040   float phimin((**(dets.begin())).surface().position().phi());
0041   float phimax(phimin);
0042
0043   for (vector<const GeomDet*>::const_iterator idet = dets.begin(); idet != dets.end(); idet++) {
0044     vector<const GeomDet*> detUnits = (**idet).components();
0045     if (!detUnits.empty()) {
0046       for (vector<const GeomDet*>::const_iterator detu = detUnits.begin(); detu != detUnits.end(); detu++) {
0047         // edm::LogInfo(TkDetLayers) << " Builder: Position of detUnit :"<< (**detu).position() ;
0048         vector<GlobalPoint> corners = computeTrapezoidalCorners(*detu);
0049         for (vector<GlobalPoint>::const_iterator i = corners.begin(); i != corners.end(); i++) {
0050           float r = i->perp();
0051           float z = i->z();
0052           float phi = i->phi();
0053           rmin = min(rmin, r);
0054           rmax = max(rmax, r);
0055           zmin = min(zmin, z);
0056           zmax = max(zmax, z);
0057           if (Geom::phiLess(phi, phimin))
0058             phimin = phi;
0059           if (Geom::phiLess(phimax, phi))
0060             phimax = phi;
0061         }
0062         // in addition to the corners we have to check the middle of the
0063         // det +/- length/2, since the min (max) radius for typical fw
0064         // dets is reached there
0065         float rdet = (**detu).position().perp();
0066         float len = (**detu).surface().bounds().length();
0067         float width = (**detu).surface().bounds().width();
0068
0069         GlobalVector xAxis = (**detu).toGlobal(LocalVector(1, 0, 0));
0070         GlobalVector yAxis = (**detu).toGlobal(LocalVector(0, 1, 0));
0071         GlobalVector perpDir = GlobalVector((**detu).position() - GlobalPoint(0, 0, (**detu).position().z()));
0072
0073         double xAxisCos = xAxis.unit().dot(perpDir.unit());
0074         double yAxisCos = yAxis.unit().dot(perpDir.unit());
0075
0076         if (std::abs(xAxisCos) > std::abs(yAxisCos)) {
0077           rmin = min(rmin, rdet - width / 2.F);
0078           rmax = max(rmax, rdet + width / 2.F);
0079         } else {
0080           rmin = min(rmin, rdet - len / 2.F);
0081           rmax = max(rmax, rdet + len / 2.F);
0082         }
0083       }
0084     } else {
0085       vector<GlobalPoint> corners = computeTrapezoidalCorners(*idet);
0086       for (vector<GlobalPoint>::const_iterator i = corners.begin(); i != corners.end(); i++) {
0087         float r = i->perp();
0088         float z = i->z();
0089         float phi = i->phi();
0090         rmin = min(rmin, r);
0091         rmax = max(rmax, r);
0092         zmin = min(zmin, z);
0093         zmax = max(zmax, z);
0094         if (Geom::phiLess(phi, phimin))
0095           phimin = phi;
0096         if (Geom::phiLess(phimax, phi))
0097           phimax = phi;
0098       }
0099       // in addition to the corners we have to check the middle of the
0100       // det +/- length/2, since the min (max) radius for typical fw
0101       // dets is reached there
0102
0103       float rdet = (**idet).position().perp();
0104       float len = (**idet).surface().bounds().length();
0105       float width = (**idet).surface().bounds().width();
0106
0107       GlobalVector xAxis = (**idet).toGlobal(LocalVector(1, 0, 0));
0108       GlobalVector yAxis = (**idet).toGlobal(LocalVector(0, 1, 0));
0109       GlobalVector perpDir = GlobalVector((**idet).position() - GlobalPoint(0, 0, (**idet).position().z()));
0110
0111       double xAxisCos = xAxis.unit().dot(perpDir.unit());
0112       double yAxisCos = yAxis.unit().dot(perpDir.unit());
0113
0114       if (std::abs(xAxisCos) > std::abs(yAxisCos)) {
0115         rmin = min(rmin, rdet - width / 2.F);
0116         rmax = max(rmax, rdet + width / 2.F);
0117       } else {
0118         rmin = min(rmin, rdet - len / 2.F);
0119         rmax = max(rmax, rdet + len / 2.F);
0120       }
0121     }
0122   }
0123
0124   if (!Geom::phiLess(phimin, phimax))
0125     edm::LogError("TkDetLayers") << " ForwardDiskSectorBuilderFromDet : "
0126                                  << "Something went wrong with Phi Sorting !";
0127   float zPos = (zmax + zmin) / 2.;
0128   float phiWin = phimax - phimin;
0129   float phiPos = (phimax + phimin) / 2.;
0130   float rmed = (rmin + rmax) / 2.;
0131   if (phiWin < 0.) {
0132     if ((phimin < Geom::pi() / 2.) || (phimax > -Geom::pi() / 2.)) {
0133       edm::LogError("TkDetLayers") << " Debug: something strange going on, please check ";
0134     }
0135     // edm::LogInfo(TkDetLayers) << " Wedge at pi: phi " << phimin << " " << phimax << " " << phiWin
0136     //   << " " << 2.*Geom::pi()+phiWin << " " ;
0137     phiWin += 2. * Geom::pi();
0138     phiPos += Geom::pi();
0139   }
0140   GlobalVector pos(rmed * cos(phiPos), rmed * sin(phiPos), zPos);
0141   return make_pair(new DiskSectorBounds(rmin, rmax, zmin - zPos, zmax - zPos, phiWin), pos);
0142 }
0143
0144 Surface::RotationType ForwardDiskSectorBuilderFromDet::computeRotation(const vector<const GeomDet*>& dets,
0145                                                                        Surface::PositionType pos) const {
0146   GlobalVector yAxis = (GlobalVector(pos.x(), pos.y(), 0.)).unit();
0147
0148   GlobalVector zAxis(0., 0., 1.);
0149   GlobalVector xAxis = yAxis.cross(zAxis);
0150
0151   return Surface::RotationType(xAxis, yAxis);
0152 }
0153
0154 vector<GlobalPoint> ForwardDiskSectorBuilderFromDet::computeTrapezoidalCorners(const GeomDet* det) const {
0155   const Plane& plane(det->specificSurface());
0156
0157   const TrapezoidalPlaneBounds* myBounds(static_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0158
0159   /*
0160   if (myBounds == 0) {
0161     string errmsg="ForwardDiskSectorBuilderFromDet: problems with dynamic cast to trapezoidal bounds for DetUnits";
0162     throw DetLayerException(errmsg);
0163     edm::LogError("TkDetLayers") << errmsg ;
0164   }
0165   */
0166
0167   // array<const float, 4>
0168   auto const& parameters = (*myBounds).parameters();
0169
0170   if (parameters[0] == 0) {
0171     edm::LogError("TkDetLayers") << "ForwardDiskSectorBuilder: something weird going on !";
0172     edm::LogError("TkDetLayers") << " Trapezoidal parameters of GeomDet (L2/L1/T/H): ";
0173     for (int i = 0; i < 4; i++)
0174       edm::LogError("TkDetLayers") << "  " << 2. * parameters[i] << "\n";
0175   }
0176
0177   float hbotedge = parameters[0];
0178   float htopedge = parameters[1];
0179   float hapothem = parameters[3];
0180   float hthick = parameters[2];
0181
0182   vector<GlobalPoint> corners;
0183
0184   corners.push_back(plane.toGlobal(LocalPoint(-htopedge, hapothem, hthick)));
0185   corners.push_back(plane.toGlobal(LocalPoint(-htopedge, hapothem, -hthick)));
0186   corners.push_back(plane.toGlobal(LocalPoint(htopedge, hapothem, hthick)));
0187   corners.push_back(plane.toGlobal(LocalPoint(htopedge, hapothem, -hthick)));
0188   corners.push_back(plane.toGlobal(LocalPoint(hbotedge, -hapothem, hthick)));
0189   corners.push_back(plane.toGlobal(LocalPoint(hbotedge, -hapothem, -hthick)));
0190   corners.push_back(plane.toGlobal(LocalPoint(-hbotedge, -hapothem, hthick)));
0191   corners.push_back(plane.toGlobal(LocalPoint(-hbotedge, -hapothem, -hthick)));
0192
0193   return corners;
0194 }