CrystalPad.cc

Go to the documentation of this file.

#include "FastSimulation/CaloGeometryTools/interface/CrystalPad.h"

#include <iostream>

std::vector<CLHEP::Hep2Vector> CrystalPad::aVector(4);

CrystalPad::CrystalPad(const CrystalPad& right) {
  corners_ = right.corners_;
  dir_ = right.dir_;
  number_ = right.number_;
  survivalProbability_ = right.survivalProbability_;
  center_ = right.center_;
  epsilon_ = right.epsilon_;
  dummy_ = right.dummy_;
}

CrystalPad& CrystalPad::operator=(const CrystalPad& right) {
  if (this != &right) {  // don't copy into yourself
    corners_ = right.corners_;
    dir_ = right.dir_;
    number_ = right.number_;
    survivalProbability_ = right.survivalProbability_;
    center_ = right.center_;
    epsilon_ = right.epsilon_;
    dummy_ = right.dummy_;
  }
  return *this;
}

CrystalPad::CrystalPad(unsigned number, const std::vector<CLHEP::Hep2Vector>& corners)
    : corners_(corners), dir_(aVector), number_(number), survivalProbability_(1.), center_(0., 0.), epsilon_(0.001) {
  //  std::cout << " Hello " << std::endl;
  if (corners.size() != 4) {
    std::cout << " Try to construct a quadrilateral with " << corners.size() << " points ! " << std::endl;
    dummy_ = true;
  } else {
    dummy_ = false;
    // Set explicity the z to 0 !
    for (unsigned ic = 0; ic < 4; ++ic) {
      dir_[ic] = (corners[(ic + 1) % 4] - corners[ic]).unit();
      center_ += corners_[ic];
    }
    center_ *= 0.25;
  }
  //  std::cout << " End of 1 constructor " << std::endl;
  //  std::cout << " Ncorners " << corners_.size() << std::endl;
  //  std::cout << " Ndirs " << dir_.size() << std::endl;
}

CrystalPad::CrystalPad(unsigned number,
                       int onEcal,
                       const std::vector<XYZPoint>& corners,
                       const XYZPoint& origin,
                       const XYZVector& vec1,
                       const XYZVector& vec2)
    : corners_(aVector), dir_(aVector), number_(number), survivalProbability_(1.), center_(0., 0.), epsilon_(0.001) {
  //  std::cout << " We are in the 2nd constructor " << std::endl;
  if (corners.size() != 4) {
    std::cout << " Try to construct a quadrilateral with " << corners.size() << " points ! " << std::endl;
    dummy_ = true;
  } else {
    dummy_ = false;
    double sign = (onEcal == 1) ? -1. : 1.;

    // the good one in the central
    trans_ = Transform3D((Point)origin,
                         (Point)(origin + vec1),
                         (Point)(origin + vec2),
                         Point(0., 0., 0.),
                         Point(0., 0., sign),
                         Point(0., 1., 0.));
    trans_.GetDecomposition(rotation_, translation_);
    //      std::cout << " Constructor 2; input corners "  << std::endl;
    for (unsigned ic = 0; ic < 4; ++ic) {
      //      std::cout << corners[ic]<< " " ;
      XYZPoint corner = rotation_(corners[ic]) + translation_;
      //      std::cout << corner << std::endl ;
      corners_[ic] = CLHEP::Hep2Vector(corner.X(), corner.Y());
      center_ += corners_[ic];
    }
    for (unsigned ic = 0; ic < 4; ++ic) {
      dir_[ic] = (corners_[(ic + 1) % 4] - corners_[ic]).unit();
    }
    center_ *= 0.25;
  }
  //  std::cout << " End of 2 constructor " << std::endl;
  //  std::cout << " Corners(constructor) " ;
  //  std::cout << corners_[0] << std::endl;
  //  std::cout << corners_[1] << std::endl;
  //  std::cout << corners_[2] << std::endl;
  //  std::cout << corners_[3] << std::endl;
}
CrystalPad::CrystalPad(
    unsigned number, const std::vector<XYZPoint>& corners, const Transform3D& trans, double scaf, bool bothdirections)
    : corners_(aVector),
      dir_(aVector),
      number_(number),
      survivalProbability_(1.),
      center_(0., 0.),
      epsilon_(0.001),
      yscalefactor_(scaf) {
  //  std::cout << " We are in the 2nd constructor " << std::endl;
  if (corners.size() != 4) {
    std::cout << " Try to construct a quadrilateral with " << corners.size() << " points ! " << std::endl;
    dummy_ = true;
  } else {
    dummy_ = false;

    // the good one in the central
    trans_ = trans;
    //      std::cout << " Constructor 2; input corners "  << std::endl;
    trans_.GetDecomposition(rotation_, translation_);
    for (unsigned ic = 0; ic < 4; ++ic) {
      XYZPoint corner = rotation_(corners[ic]) + translation_;
      //      std::cout << corner << std::endl ;
      double xscalefactor = (bothdirections) ? yscalefactor_ : 1.;
      corners_[ic] = CLHEP::Hep2Vector(corner.X() * xscalefactor, corner.Y() * yscalefactor_);
      center_ += corners_[ic];
    }
    for (unsigned ic = 0; ic < 4; ++ic) {
      dir_[ic] = (corners_[(ic + 1) % 4] - corners_[ic]).unit();
    }
    center_ *= 0.25;
  }
}

bool CrystalPad::inside(const CLHEP::Hep2Vector& ppoint, bool debug) const {
  //  std::cout << "Inside " << ppoint <<std::endl;
  //  std::cout << "Corners " << corners_.size() << std::endl;
  //  std::cout << corners_[0] << std::endl;
  //  std::cout << corners_[1] << std::endl;
  //  std::cout << corners_[2] << std::endl;
  //  std::cout << corners_[3] << std::endl;
  //  std::cout << " Got the 2D point " << std::endl;
  CLHEP::Hep2Vector pv0(ppoint - corners_[0]);
  CLHEP::Hep2Vector pv2(ppoint - corners_[2]);
  CLHEP::Hep2Vector n1(pv0 - (pv0 * dir_[0]) * dir_[0]);
  CLHEP::Hep2Vector n2(pv2 - (pv2 * dir_[2]) * dir_[2]);

  //  double N1(n1.mag());
  //  double N2(n2.mag());
  double r1(n1 * n2);
  bool inside1(r1 <= 0.);

  if (!inside1)
    return false;

  //  if(debug)
  //    {
  //      std::cout << n1 << std::endl;
  //      std::cout << n2 << std::endl;
  //      std::cout << r1 << std::endl;
  //      std::cout << inside1 << std::endl;
  //    }

  //  bool close1=(N1<epsilon_||N2<epsilon_);
  //
  //  if(!close1&&!inside1) return false;
  //  std::cout << " First calculation " << std::endl;
  CLHEP::Hep2Vector pv1(ppoint - corners_[1]);
  CLHEP::Hep2Vector pv3(ppoint - corners_[3]);
  CLHEP::Hep2Vector n3(pv1 - (pv1 * dir_[1]) * dir_[1]);
  CLHEP::Hep2Vector n4(pv3 - (pv3 * dir_[3]) * dir_[3]);
  //  double N3(n3.mag());
  //  double N4(n4.mag());
  //  bool close2=(N3<epsilon_||N4<epsilon_);
  double r2(n3 * n4);
  bool inside2(r2 <= 0.);
  //  //  std::cout << " pv1 & pv3 " << pv1.mag() << " " << pv3.mag() << std::endl;
  //  //  double tmp=(pv1-(pv1*dir_[1])*dir_[1])*(pv3-(pv3*dir_[3])*dir_[3]);
  //  //  std::cout << " Computed tmp " << tmp << std::endl;
  //  if(debug)
  //    {
  //      std::cout << n3 << std::endl;
  //      std::cout << n4 << std::endl;
  //      std::cout << r2 << std::endl;
  //      std::cout << inside2 << std::endl;
  //    }
  //  if(!close2&&!inside2) return false;
  //  std::cout << " Second calculation " << std::endl;
  //  std::cout << " True " << std::endl;
  //    return (!close1&&!close2||(close2&&inside1||close1&&inside2));

  return inside2;
}

/*
bool 
CrystalPad::globalinside(XYZPoint point) const
{
  //  std::cout << " Global inside " << std::endl;
  //  std::cout << point << " " ;
  ROOT::Math::Rotation3D r;
  XYZVector t;
  point = rotation_(point)+translation_;
  //  std::cout << point << std::endl;
  //  print();
  CLHEP::Hep2Vector ppoint(point.X(),point.Y());
  bool result=inside(ppoint);
  //  std::cout << " Result " << result << std::endl;
  return result;
}
*/

void CrystalPad::print() const {
  std::cout << " Corners " << std::endl;
  std::cout << corners_[0] << std::endl;
  std::cout << corners_[1] << std::endl;
  std::cout << corners_[2] << std::endl;
  std::cout << corners_[3] << std::endl;
}

/*
CLHEP::Hep2Vector 
CrystalPad::localPoint(XYZPoint point) const
{
  point = rotation_(point)+translation_;
  return CLHEP::Hep2Vector(point.X(),point.Y());
}
*/

CLHEP::Hep2Vector& CrystalPad::edge(unsigned iside, int n) { return corners_[(iside + n) % 4]; }

CLHEP::Hep2Vector& CrystalPad::edge(CaloDirection dir) {
  switch (dir) {
    case NORTHWEST:
      return corners_[0];
      break;
    case NORTHEAST:
      return corners_[1];
      break;
    case SOUTHEAST:
      return corners_[2];
      break;
    case SOUTHWEST:
      return corners_[3];
      break;
    default: {
      std::cout << " Serious problem in CrystalPad ! " << dir << std::endl;
      return corners_[0];
    }
  }
  return corners_[0];
}

void CrystalPad::extrems(double& xmin, double& xmax, double& ymin, double& ymax) const {
  xmin = ymin = 999;
  xmax = ymax = -999;
  for (unsigned ic = 0; ic < 4; ++ic) {
    if (corners_[ic].x() < xmin)
      xmin = corners_[ic].x();
    if (corners_[ic].x() > xmax)
      xmax = corners_[ic].x();
    if (corners_[ic].y() < ymin)
      ymin = corners_[ic].y();
    if (corners_[ic].y() > ymax)
      ymax = corners_[ic].y();
  }
}

void CrystalPad::resetCorners() {
  // Find the centre-of-gravity of the Quad (after re-organization)
  center_ = CLHEP::Hep2Vector(0., 0.);
  for (unsigned ic = 0; ic < 4; ++ic)
    center_ += corners_[ic];
  center_ *= 0.25;

  // Rescale the corners to allow for some inaccuracies in
  // in the inside test
  for (unsigned ic = 0; ic < 4; ++ic)
    corners_[ic] += 0.001 * (corners_[ic] - center_);
}

std::ostream& operator<<(std::ostream& ost, CrystalPad& quad) {
  ost << " Number " << quad.getNumber() << std::endl;
  ost << NORTHWEST << quad.edge(NORTHWEST) << std::endl;
  ost << NORTHEAST << quad.edge(NORTHEAST) << std::endl;
  ost << SOUTHEAST << quad.edge(SOUTHEAST) << std::endl;
  ost << SOUTHWEST << quad.edge(SOUTHWEST) << std::endl;

  return ost;
}

void CrystalPad::getDrawingCoordinates(std::vector<float>& x, std::vector<float>& y) const {
  x.clear();
  y.clear();
  x.push_back(corners_[0].x());
  x.push_back(corners_[1].x());
  x.push_back(corners_[2].x());
  x.push_back(corners_[3].x());
  x.push_back(corners_[0].x());
  y.push_back(corners_[0].y());
  y.push_back(corners_[1].y());
  y.push_back(corners_[2].y());
  y.push_back(corners_[3].y());
  y.push_back(corners_[0].y());
}