CaloGeometryTools/src/BaseCrystal.cc

0001 //FAMOS headers
0002 #include "FastSimulation/CaloGeometryTools/interface/BaseCrystal.h"
0003
0004 // Data Formats
0005 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0006
0007 typedef ROOT::Math::Plane3D::Point Point;
0008
0009 BaseCrystal::BaseCrystal(const DetId &cell) : cellid_(cell) { subdetn_ = cellid_.subdetId(); }
0010
0011 void BaseCrystal::setCorners(const CaloCellGeometry::CornersVec &vec, const GlobalPoint &pos) {
0012   unsigned ncorners = vec.size();
0013   if (ncorners != 8)
0014     return;
0015
0016   // This is really a pity to have to make the conversion GlobalPoint to XYZPoint, but the latter have many
0017   // useful properties (transformations, conversions....) that aren't implemented
0018   // for the GlobalPoints and GlobalVectors.
0019
0020   if (subdetn_ == EcalBarrel) {
0021     if (pos.z() > 0.) {
0022       for (unsigned ic = 0; ic < 8; ++ic) {
0023         corners_[ic] = XYZPoint(vec[ic].x(), vec[ic].y(), vec[ic].z());
0024       }
0025     } else {
0026       corners_[0] = XYZPoint(vec[2].x(), vec[2].y(), vec[2].z());
0027       corners_[1] = XYZPoint(vec[3].x(), vec[3].y(), vec[3].z());
0028       corners_[2] = XYZPoint(vec[0].x(), vec[0].y(), vec[0].z());
0029       corners_[3] = XYZPoint(vec[1].x(), vec[1].y(), vec[1].z());
0030       corners_[4] = XYZPoint(vec[6].x(), vec[6].y(), vec[6].z());
0031       corners_[5] = XYZPoint(vec[7].x(), vec[7].y(), vec[7].z());
0032       corners_[6] = XYZPoint(vec[4].x(), vec[4].y(), vec[4].z());
0033       corners_[7] = XYZPoint(vec[5].x(), vec[5].y(), vec[5].z());
0034     }
0035   } else if (subdetn_ == EcalEndcap) {
0036     double x = pos.x();
0037     double y = pos.y();
0038     double z = pos.z();
0039     unsigned offset = 0;
0040     int zsign = 1;
0041     if (z > 0) {
0042       if (x > 0 && y > 0)
0043         offset = 1;
0044       else if (x < 0 && y > 0)
0045         offset = 2;
0046       else if (x > 0 && y < 0)
0047         offset = 0;
0048       else if (x < 0 && y < 0)
0049         offset = 3;
0050       zsign = 1;
0051     } else {
0052       if (x > 0 && y > 0)
0053         offset = 3;
0054       else if (x < 0 && y > 0)
0055         offset = 2;
0056       else if (x > 0 && y < 0)
0057         offset = 0;
0058       else if (x < 0 && y < 0)
0059         offset = 1;
0060       zsign = -1;
0061     }
0062     for (unsigned ic = 0; ic < 4; ++ic) {
0063       unsigned i1 = (unsigned)((zsign * ic + offset) % 4);
0064       unsigned i2 = i1 + 4;
0065       corners_[ic] = XYZPoint(vec[i1].x(), vec[i1].y(), vec[i1].z());
0066       corners_[4 + ic] = XYZPoint(vec[i2].x(), vec[i2].y(), vec[i2].z());
0067     }
0068   }
0069   computeBasicProperties();
0070 }
0071
0072 void BaseCrystal::computeBasicProperties() {
0073   //if(corners_.size()==0) return;
0074   center_ = XYZPoint(0., 0., 0.);
0075   for (unsigned ic = 0; ic < 8; ++ic) {
0076     center_ += corners_[ic];
0077   }
0078
0079   center_ *= 0.125;
0080
0081   //  std::cout << " Ncorners ? " << corners_.size() << std::endl;
0082   frontcenter_ = 0.25 * (corners_[0] + corners_[1] + corners_[2] + corners_[3]);
0083   backcenter_ = 0.25 * (corners_[4] + corners_[5] + corners_[6] + corners_[7]);
0084   crystalaxis_ = backcenter_ - frontcenter_;
0085   firstedgedirection_ = -(corners_[1] - corners_[0]).Unit();
0086   fifthedgedirection_ = -(corners_[5] - corners_[4]).Unit();
0087   //  std::cout << " Direction laterales " << std::endl;
0088   for (unsigned il = 0; il < 4; ++il) {
0089     lateraldirection_[il] = -(corners_[(il + 1) % 4] - corners_[il]).Unit();
0090   }
0091
0092   Plane3D frontPlane((Point)corners_[0], (Point)corners_[1], (Point)corners_[2]);
0093   Plane3D backPlane((Point)corners_[4], (Point)corners_[5], (Point)corners_[6]);
0094   for (unsigned i = 0; i < 4; ++i) {
0095     lateralPlane_[i] = Plane3D((Point)corners_[i], (Point)corners_[(i + 1) % 4], (Point)corners_[i + 4]);
0096   }
0097   // Front plane i=4 (UP)
0098   lateralPlane_[4] = frontPlane;
0099   // Back plane i =5 (DOWN)
0100   lateralPlane_[5] = backPlane;
0101
0102   for (unsigned i = 0; i < 6; ++i) {
0103     exitingNormal_[i] = (lateralPlane_[i].Distance(Point(center_.X(), center_.Y(), center_.Z())) < 0.)
0104                             ? lateralPlane_[i].Normal().Unit()
0105                             : -lateralPlane_[i].Normal().Unit();
0106   }
0107 }
0108
0109 void BaseCrystal::getLateralEdges(unsigned i, XYZPoint &a, XYZPoint &b) const {
0110   if (i < 4U)  // i >= 0, since i is unsigned
0111   {
0112     a = corners_[i];
0113     b = corners_[i + 4];
0114   }
0115 }
0116
0117 void BaseCrystal::getFrontSide(XYZPoint &a, XYZPoint &b, XYZPoint &c, XYZPoint &d) const {
0118   a = corners_[0];
0119   b = corners_[1];
0120   c = corners_[2];
0121   d = corners_[3];
0122 }
0123
0124 void BaseCrystal::getFrontSide(std::vector<XYZPoint> &corners) const {
0125   if (corners.size() == 4) {
0126     corners[0] = corners_[0];
0127     corners[1] = corners_[1];
0128     corners[2] = corners_[2];
0129     corners[3] = corners_[3];
0130   }
0131 }
0132
0133 void BaseCrystal::getBackSide(XYZPoint &a, XYZPoint &b, XYZPoint &c, XYZPoint &d) const {
0134   a = corners_[4];
0135   b = corners_[5];
0136   c = corners_[6];
0137   d = corners_[7];
0138 }
0139
0140 void BaseCrystal::getBackSide(std::vector<XYZPoint> &corners) const {
0141   if (corners.size() == 4) {
0142     corners[0] = corners_[4];
0143     corners[1] = corners_[5];
0144     corners[2] = corners_[6];
0145     corners[3] = corners_[7];
0146   }
0147 }
0148
0149 void BaseCrystal::getLateralSide(unsigned i, XYZPoint &a, XYZPoint &b, XYZPoint &c, XYZPoint &d) const {
0150   if (i < 4U)  // i >= 0, since i is unsigned
0151   {
0152     getLateralEdges(i, a, b);
0153     getLateralEdges((i + 1) % 4, c, d);
0154   }
0155 }
0156
0157 void BaseCrystal::getLateralSide(unsigned i, std::vector<XYZPoint> &corners) const {
0158   if (corners.size() == 4 && i < 4U)  // i >= 0, since i is unsigned
0159   {
0160     corners[0] = corners_[i];
0161     corners[1] = corners_[i + 4];
0162     corners[2] = corners_[4 + (i + 1) % 4];
0163     corners[3] = corners_[(i + 1) % 4];
0164   }
0165 }
0166
0167 void BaseCrystal::getDrawingCoordinates(std::vector<float> &x, std::vector<float> &y, std::vector<float> &z) const {
0168   x.clear();
0169   y.clear();
0170   z.clear();
0171
0172   x.push_back(corners_[0].X());
0173   x.push_back(corners_[3].X());
0174   x.push_back(corners_[2].X());
0175   x.push_back(corners_[1].X());
0176   x.push_back(corners_[5].X());
0177   x.push_back(corners_[6].X());
0178   x.push_back(corners_[7].X());
0179   x.push_back(corners_[4].X());
0180   x.push_back(corners_[0].X());
0181   x.push_back(corners_[1].X());
0182   x.push_back(corners_[2].X());
0183   x.push_back(corners_[6].X());
0184   x.push_back(corners_[5].X());
0185   x.push_back(corners_[4].X());
0186   x.push_back(corners_[7].X());
0187   x.push_back(corners_[3].X());
0188
0189   y.push_back(corners_[0].Y());
0190   y.push_back(corners_[3].Y());
0191   y.push_back(corners_[2].Y());
0192   y.push_back(corners_[1].Y());
0193   y.push_back(corners_[5].Y());
0194   y.push_back(corners_[6].Y());
0195   y.push_back(corners_[7].Y());
0196   y.push_back(corners_[4].Y());
0197   y.push_back(corners_[0].Y());
0198   y.push_back(corners_[1].Y());
0199   y.push_back(corners_[2].Y());
0200   y.push_back(corners_[6].Y());
0201   y.push_back(corners_[5].Y());
0202   y.push_back(corners_[4].Y());
0203   y.push_back(corners_[7].Y());
0204   y.push_back(corners_[3].Y());
0205
0206   z.push_back(corners_[0].Z());
0207   z.push_back(corners_[3].Z());
0208   z.push_back(corners_[2].Z());
0209   z.push_back(corners_[1].Z());
0210   z.push_back(corners_[5].Z());
0211   z.push_back(corners_[6].Z());
0212   z.push_back(corners_[7].Z());
0213   z.push_back(corners_[4].Z());
0214   z.push_back(corners_[0].Z());
0215   z.push_back(corners_[1].Z());
0216   z.push_back(corners_[2].Z());
0217   z.push_back(corners_[6].Z());
0218   z.push_back(corners_[5].Z());
0219   z.push_back(corners_[4].Z());
0220   z.push_back(corners_[7].Z());
0221   z.push_back(corners_[3].Z());
0222 }
0223
0224 void BaseCrystal::getSide(const CaloDirection &side, XYZPoint &a, XYZPoint &b, XYZPoint &c, XYZPoint &d) const {
0225   switch (side) {
0226     case UP:
0227       getFrontSide(a, b, c, d);
0228       break;
0229     case DOWN:
0230       getBackSide(a, b, c, d);
0231       break;
0232     default:
0233       getLateralSide(CaloDirectionOperations::Side(side), a, b, c, d);
0234   }
0235 }
0236
0237 void BaseCrystal::print() const {
0238   std::cout << "CellID " << cellid_.rawId() << std::endl;
0239   std::cout << " Corners " << std::endl;
0240   for (unsigned ic = 0; ic < 8; ++ic)
0241     std::cout << corners_[ic] << std::endl;
0242   std::cout << " Center " << center_ << std::endl;
0243   std::cout << " Front Center " << frontcenter_ << std::endl;
0244   std::cout << " Back Center " << backcenter_ << std::endl;
0245   std::cout << " Normales sortantes " << std::endl;
0246   for (unsigned id = 0; id < 6; ++id)
0247     std::cout << exitingNormal_[id] << std::endl;
0248 }
0249
0250 void BaseCrystal::getSide(const CaloDirection &side, std::vector<XYZPoint> &corners) const {
0251   switch (side) {
0252     case UP:
0253       getFrontSide(corners);
0254       break;
0255     case DOWN:
0256       getBackSide(corners);
0257       break;
0258     default:
0259       getLateralSide(CaloDirectionOperations::Side(side), corners);
0260   }
0261 }