BaseCrystal.cc

Go to the documentation of this file.

//FAMOS headers
#include "FastSimulation/CaloGeometryTools/interface/BaseCrystal.h"

// Data Formats
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"

typedef ROOT::Math::Plane3D::Point Point;

BaseCrystal::BaseCrystal(const DetId &cell) : cellid_(cell) { subdetn_ = cellid_.subdetId(); }

void BaseCrystal::setCorners(const CaloCellGeometry::CornersVec &vec, const GlobalPoint &pos) {
  unsigned ncorners = vec.size();
  if (ncorners != 8)
    return;

  // This is really a pity to have to make the conversion GlobalPoint to XYZPoint, but the latter have many
  // useful properties (transformations, conversions....) that aren't implemented
  // for the GlobalPoints and GlobalVectors.

  if (subdetn_ == EcalBarrel) {
    if (pos.z() > 0.) {
      for (unsigned ic = 0; ic < 8; ++ic) {
        corners_[ic] = XYZPoint(vec[ic].x(), vec[ic].y(), vec[ic].z());
      }
    } else {
      corners_[0] = XYZPoint(vec[2].x(), vec[2].y(), vec[2].z());
      corners_[1] = XYZPoint(vec[3].x(), vec[3].y(), vec[3].z());
      corners_[2] = XYZPoint(vec[0].x(), vec[0].y(), vec[0].z());
      corners_[3] = XYZPoint(vec[1].x(), vec[1].y(), vec[1].z());
      corners_[4] = XYZPoint(vec[6].x(), vec[6].y(), vec[6].z());
      corners_[5] = XYZPoint(vec[7].x(), vec[7].y(), vec[7].z());
      corners_[6] = XYZPoint(vec[4].x(), vec[4].y(), vec[4].z());
      corners_[7] = XYZPoint(vec[5].x(), vec[5].y(), vec[5].z());
    }
  } else if (subdetn_ == EcalEndcap) {
    double x = pos.x();
    double y = pos.y();
    double z = pos.z();
    unsigned offset = 0;
    int zsign = 1;
    if (z > 0) {
      if (x > 0 && y > 0)
        offset = 1;
      else if (x < 0 && y > 0)
        offset = 2;
      else if (x > 0 && y < 0)
        offset = 0;
      else if (x < 0 && y < 0)
        offset = 3;
      zsign = 1;
    } else {
      if (x > 0 && y > 0)
        offset = 3;
      else if (x < 0 && y > 0)
        offset = 2;
      else if (x > 0 && y < 0)
        offset = 0;
      else if (x < 0 && y < 0)
        offset = 1;
      zsign = -1;
    }
    for (unsigned ic = 0; ic < 4; ++ic) {
      unsigned i1 = (unsigned)((zsign * ic + offset) % 4);
      unsigned i2 = i1 + 4;
      corners_[ic] = XYZPoint(vec[i1].x(), vec[i1].y(), vec[i1].z());
      corners_[4 + ic] = XYZPoint(vec[i2].x(), vec[i2].y(), vec[i2].z());
    }
  }
  computeBasicProperties();
}

void BaseCrystal::computeBasicProperties() {
  //if(corners_.size()==0) return;
  center_ = XYZPoint(0., 0., 0.);
  for (unsigned ic = 0; ic < 8; ++ic) {
    center_ += corners_[ic];
  }

  center_ *= 0.125;

  //  std::cout << " Ncorners ? " << corners_.size() << std::endl;
  frontcenter_ = 0.25 * (corners_[0] + corners_[1] + corners_[2] + corners_[3]);
  backcenter_ = 0.25 * (corners_[4] + corners_[5] + corners_[6] + corners_[7]);
  crystalaxis_ = backcenter_ - frontcenter_;
  firstedgedirection_ = -(corners_[1] - corners_[0]).Unit();
  fifthedgedirection_ = -(corners_[5] - corners_[4]).Unit();
  //  std::cout << " Direction laterales " << std::endl;
  for (unsigned il = 0; il < 4; ++il) {
    lateraldirection_[il] = -(corners_[(il + 1) % 4] - corners_[il]).Unit();
  }

  Plane3D frontPlane((Point)corners_[0], (Point)corners_[1], (Point)corners_[2]);
  Plane3D backPlane((Point)corners_[4], (Point)corners_[5], (Point)corners_[6]);
  for (unsigned i = 0; i < 4; ++i) {
    lateralPlane_[i] = Plane3D((Point)corners_[i], (Point)corners_[(i + 1) % 4], (Point)corners_[i + 4]);
  }
  // Front plane i=4 (UP)
  lateralPlane_[4] = frontPlane;
  // Back plane i =5 (DOWN)
  lateralPlane_[5] = backPlane;

  for (unsigned i = 0; i < 6; ++i) {
    exitingNormal_[i] = (lateralPlane_[i].Distance(Point(center_.X(), center_.Y(), center_.Z())) < 0.)
                            ? lateralPlane_[i].Normal().Unit()
                            : -lateralPlane_[i].Normal().Unit();
  }
}

void BaseCrystal::getLateralEdges(unsigned i, XYZPoint &a, XYZPoint &b) const {
  if (i < 4U)  // i >= 0, since i is unsigned
  {
    a = corners_[i];
    b = corners_[i + 4];
  }
}

void BaseCrystal::getFrontSide(XYZPoint &a, XYZPoint &b, XYZPoint &c, XYZPoint &d) const {
  a = corners_[0];
  b = corners_[1];
  c = corners_[2];
  d = corners_[3];
}

void BaseCrystal::getFrontSide(std::vector<XYZPoint> &corners) const {
  if (corners.size() == 4) {
    corners[0] = corners_[0];
    corners[1] = corners_[1];
    corners[2] = corners_[2];
    corners[3] = corners_[3];
  }
}

void BaseCrystal::getBackSide(XYZPoint &a, XYZPoint &b, XYZPoint &c, XYZPoint &d) const {
  a = corners_[4];
  b = corners_[5];
  c = corners_[6];
  d = corners_[7];
}

void BaseCrystal::getBackSide(std::vector<XYZPoint> &corners) const {
  if (corners.size() == 4) {
    corners[0] = corners_[4];
    corners[1] = corners_[5];
    corners[2] = corners_[6];
    corners[3] = corners_[7];
  }
}

void BaseCrystal::getLateralSide(unsigned i, XYZPoint &a, XYZPoint &b, XYZPoint &c, XYZPoint &d) const {
  if (i < 4U)  // i >= 0, since i is unsigned
  {
    getLateralEdges(i, a, b);
    getLateralEdges((i + 1) % 4, c, d);
  }
}

void BaseCrystal::getLateralSide(unsigned i, std::vector<XYZPoint> &corners) const {
  if (corners.size() == 4 && i < 4U)  // i >= 0, since i is unsigned
  {
    corners[0] = corners_[i];
    corners[1] = corners_[i + 4];
    corners[2] = corners_[4 + (i + 1) % 4];
    corners[3] = corners_[(i + 1) % 4];
  }
}

void BaseCrystal::getDrawingCoordinates(std::vector<float> &x, std::vector<float> &y, std::vector<float> &z) const {
  x.clear();
  y.clear();
  z.clear();

  x.push_back(corners_[0].X());
  x.push_back(corners_[3].X());
  x.push_back(corners_[2].X());
  x.push_back(corners_[1].X());
  x.push_back(corners_[5].X());
  x.push_back(corners_[6].X());
  x.push_back(corners_[7].X());
  x.push_back(corners_[4].X());
  x.push_back(corners_[0].X());
  x.push_back(corners_[1].X());
  x.push_back(corners_[2].X());
  x.push_back(corners_[6].X());
  x.push_back(corners_[5].X());
  x.push_back(corners_[4].X());
  x.push_back(corners_[7].X());
  x.push_back(corners_[3].X());

  y.push_back(corners_[0].Y());
  y.push_back(corners_[3].Y());
  y.push_back(corners_[2].Y());
  y.push_back(corners_[1].Y());
  y.push_back(corners_[5].Y());
  y.push_back(corners_[6].Y());
  y.push_back(corners_[7].Y());
  y.push_back(corners_[4].Y());
  y.push_back(corners_[0].Y());
  y.push_back(corners_[1].Y());
  y.push_back(corners_[2].Y());
  y.push_back(corners_[6].Y());
  y.push_back(corners_[5].Y());
  y.push_back(corners_[4].Y());
  y.push_back(corners_[7].Y());
  y.push_back(corners_[3].Y());

  z.push_back(corners_[0].Z());
  z.push_back(corners_[3].Z());
  z.push_back(corners_[2].Z());
  z.push_back(corners_[1].Z());
  z.push_back(corners_[5].Z());
  z.push_back(corners_[6].Z());
  z.push_back(corners_[7].Z());
  z.push_back(corners_[4].Z());
  z.push_back(corners_[0].Z());
  z.push_back(corners_[1].Z());
  z.push_back(corners_[2].Z());
  z.push_back(corners_[6].Z());
  z.push_back(corners_[5].Z());
  z.push_back(corners_[4].Z());
  z.push_back(corners_[7].Z());
  z.push_back(corners_[3].Z());
}

void BaseCrystal::getSide(const CaloDirection &side, XYZPoint &a, XYZPoint &b, XYZPoint &c, XYZPoint &d) const {
  switch (side) {
    case UP:
      getFrontSide(a, b, c, d);
      break;
    case DOWN:
      getBackSide(a, b, c, d);
      break;
    default:
      getLateralSide(CaloDirectionOperations::Side(side), a, b, c, d);
  }
}

void BaseCrystal::print() const {
  std::cout << "CellID " << cellid_.rawId() << std::endl;
  std::cout << " Corners " << std::endl;
  for (unsigned ic = 0; ic < 8; ++ic)
    std::cout << corners_[ic] << std::endl;
  std::cout << " Center " << center_ << std::endl;
  std::cout << " Front Center " << frontcenter_ << std::endl;
  std::cout << " Back Center " << backcenter_ << std::endl;
  std::cout << " Normales sortantes " << std::endl;
  for (unsigned id = 0; id < 6; ++id)
    std::cout << exitingNormal_[id] << std::endl;
}

void BaseCrystal::getSide(const CaloDirection &side, std::vector<XYZPoint> &corners) const {
  switch (side) {
    case UP:
      getFrontSide(corners);
      break;
    case DOWN:
      getBackSide(corners);
      break;
    default:
      getLateralSide(CaloDirectionOperations::Side(side), corners);
  }
}