GeomBuilder/src/BaseVolumeHandle.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author N. Amapane - INFN Torino (original developer)
0005  */
0006
0007 #include "BaseVolumeHandle.h"
0008
0009 #include "DataFormats/GeometrySurface/interface/Plane.h"
0010 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0011 #include "DataFormats/GeometrySurface/interface/Cone.h"
0012 #include "DataFormats/GeometryVector/interface/CoordinateSets.h"
0013
0014 #include <cassert>
0015 #include <string>
0016 #include <iterator>
0017 #include <iomanip>
0018 #include <iostream>
0019
0020 using namespace SurfaceOrientation;
0021 using namespace std;
0022 using namespace magneticfield;
0023
0024 BaseVolumeHandle::BaseVolumeHandle(bool expand2Pi, bool debugVal)
0025     : magVolume(nullptr),
0026       masterSector(1),
0027       theRN(0.),
0028       theRMin(0.),
0029       theRMax(0.),
0030       refPlane(nullptr),
0031       expand(expand2Pi),
0032       isIronFlag(false),
0033       debug(debugVal) {}
0034
0035 BaseVolumeHandle::~BaseVolumeHandle() {
0036   if (refPlane != nullptr) {
0037     delete refPlane;
0038     refPlane = nullptr;
0039   }
0040 }
0041
0042 const Surface::GlobalPoint &BaseVolumeHandle::center() const { return center_; }
0043
0044 void BaseVolumeHandle::buildPhiZSurf(double startPhi, double deltaPhi, double zhalf, double rCentr) {
0045   // This is 100% equal for cons and tubs!!!
0046
0047   GlobalVector planeXAxis = refPlane->toGlobal(LocalVector(1, 0, 0));
0048   GlobalVector planeYAxis = refPlane->toGlobal(LocalVector(0, 1, 0));
0049   GlobalVector planeZAxis = refPlane->toGlobal(LocalVector(0, 0, 1));
0050
0051   // Local Y axis of the faces at +-phi.
0052   GlobalVector y_phiplus = refPlane->toGlobal(LocalVector(cos(startPhi + deltaPhi), sin(startPhi + deltaPhi), 0.));
0053   GlobalVector y_phiminus = refPlane->toGlobal(LocalVector(cos(startPhi), sin(startPhi), 0.));
0054
0055   Surface::RotationType rot_Z(planeXAxis, planeYAxis);
0056   Surface::RotationType rot_phiplus(planeZAxis, y_phiplus);
0057   Surface::RotationType rot_phiminus(planeZAxis, y_phiminus);
0058
0059   GlobalPoint pos_zplus(center_.x(), center_.y(), center_.z() + zhalf);
0060   GlobalPoint pos_zminus(center_.x(), center_.y(), center_.z() - zhalf);
0061   // BEWARE: in this case, the origin for phiplus,phiminus surfaces is
0062   // at radius R and not on a plane passing by center_ orthogonal to the radius.
0063   GlobalPoint pos_phiplus(
0064       refPlane->toGlobal(LocalPoint(rCentr * cos(startPhi + deltaPhi), rCentr * sin(startPhi + deltaPhi), 0.)));
0065   GlobalPoint pos_phiminus(refPlane->toGlobal(LocalPoint(rCentr * cos(startPhi), rCentr * sin(startPhi), 0.)));
0066   surfaces[zplus] = new Plane(pos_zplus, rot_Z);
0067   surfaces[zminus] = new Plane(pos_zminus, rot_Z);
0068   surfaces[phiplus] = new Plane(pos_phiplus, rot_phiplus);
0069   surfaces[phiminus] = new Plane(pos_phiminus, rot_phiminus);
0070
0071   if (debug) {
0072     cout << "Actual Center at: " << center_ << " R " << center_.perp() << " phi " << center_.phi() << endl;
0073     cout << "RN            " << theRN << endl;
0074
0075     cout << "pos_zplus    " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl
0076          << "pos_zminus   " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl
0077          << "pos_phiplus  " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl
0078          << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl;
0079
0080     cout << "y_phiplus  " << y_phiplus << endl;
0081     cout << "y_phiminus " << y_phiminus << endl;
0082
0083     cout << "rot_Z    " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl
0084          << "rot_phi+ " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi "
0085          << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl
0086          << "rot_phi- " << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi "
0087          << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl;
0088   }
0089
0090   //   // Check ordering.
0091   if (debug) {
0092     if (pos_zplus.z() < pos_zminus.z()) {
0093       cout << "*** WARNING: pos_zplus < pos_zminus " << endl;
0094     }
0095     if (Geom::Phi<float>(pos_phiplus.phi() - pos_phiminus.phi()) < 0.) {
0096       cout << "*** WARNING: pos_phiplus < pos_phiminus " << endl;
0097     }
0098   }
0099 }
0100
0101 bool BaseVolumeHandle::sameSurface(const Surface &s1, Sides which_side, float tolerance) {
0102   //Check for null comparison
0103   if (&s1 == (surfaces[which_side]).get()) {
0104     if (debug)
0105       cout << "      sameSurface: OK (same ptr)" << endl;
0106     return true;
0107   }
0108
0109   const float maxtilt = 0.999;
0110
0111   const Surface &s2 = *(surfaces[which_side]);
0112   // Try with a plane.
0113   const Plane *p1 = dynamic_cast<const Plane *>(&s1);
0114   if (p1 != nullptr) {
0115     const Plane *p2 = dynamic_cast<const Plane *>(&s2);
0116     if (p2 == nullptr) {
0117       if (debug)
0118         cout << "      sameSurface: different types" << endl;
0119       return false;
0120     }
0121
0122     if ((fabs(p1->normalVector().dot(p2->normalVector())) > maxtilt) &&
0123         (fabs((p1->toLocal(p2->position())).z()) < tolerance)) {
0124       if (debug)
0125         cout << "      sameSurface: OK " << fabs(p1->normalVector().dot(p2->normalVector())) << " "
0126              << fabs((p1->toLocal(p2->position())).z()) << endl;
0127       return true;
0128     } else {
0129       if (debug)
0130         cout << "      sameSurface: not the same: " << p1->normalVector() << p1->position() << endl
0131              << "                                 " << p2->normalVector() << p2->position() << endl
0132              << fabs(p1->normalVector().dot(p2->normalVector())) << " " << (p1->toLocal(p2->position())).z() << endl;
0133       return false;
0134     }
0135   }
0136
0137   // Try with a cylinder.
0138   const Cylinder *cy1 = dynamic_cast<const Cylinder *>(&s1);
0139   if (cy1 != nullptr) {
0140     const Cylinder *cy2 = dynamic_cast<const Cylinder *>(&s2);
0141     if (cy2 == nullptr) {
0142       if (debug)
0143         cout << "      sameSurface: different types" << endl;
0144       return false;
0145     }
0146     // Assume axis is the same!
0147     if (fabs(cy1->radius() - cy2->radius()) < tolerance) {
0148       return true;
0149     } else {
0150       return false;
0151     }
0152   }
0153
0154   // Try with a cone.
0155   const Cone *co1 = dynamic_cast<const Cone *>(&s1);
0156   if (co1 != nullptr) {
0157     const Cone *co2 = dynamic_cast<const Cone *>(&s2);
0158     if (co2 == nullptr) {
0159       if (debug)
0160         cout << "      sameSurface: different types" << endl;
0161       return false;
0162     }
0163     // FIXME
0164     if (fabs(co1->openingAngle() - co2->openingAngle()) < maxtilt &&
0165         (co1->vertex() - co2->vertex()).mag() < tolerance) {
0166       return true;
0167     } else {
0168       return false;
0169     }
0170   }
0171
0172   if (debug)
0173     cout << "      sameSurface: unknown surfaces..." << endl;
0174   return false;
0175 }
0176
0177 bool BaseVolumeHandle::setSurface(const Surface &s1, Sides which_side) {
0178   //Check for null assignment
0179   if (&s1 == (surfaces[which_side]).get()) {
0180     isAssigned[which_side] = true;
0181     return true;
0182   }
0183
0184   if (!sameSurface(s1, which_side)) {
0185     cout << "***ERROR: setSurface: trying to assign a surface that does not match destination surface. Skipping."
0186          << endl;
0187     const Surface &s2 = *(surfaces[which_side]);
0188     //FIXME: Just planes for the time being!!!
0189     const Plane *p1 = dynamic_cast<const Plane *>(&s1);
0190     const Plane *p2 = dynamic_cast<const Plane *>(&s2);
0191     if (p1 != nullptr && p2 != nullptr)
0192       cout << p1->normalVector() << p1->position() << endl << p2->normalVector() << p2->position() << endl;
0193     return false;
0194   }
0195
0196   if (isAssigned[which_side]) {
0197     if (&s1 != (surfaces[which_side]).get()) {
0198       cout << "*** WARNING BaseVolumeHandle::setSurface: trying to reassign a surface to a different surface instance"
0199            << endl;
0200       return false;
0201     }
0202   } else {
0203     surfaces[which_side] = &s1;
0204     isAssigned[which_side] = true;
0205     if (debug)
0206       cout << "     Volume " << name << " # " << copyno << " Assigned: " << (int)which_side << endl;
0207     return true;
0208   }
0209
0210   return false;  // let the compiler be happy
0211 }
0212
0213 const Surface &BaseVolumeHandle::surface(Sides which_side) const { return *(surfaces[which_side]); }
0214
0215 const Surface &BaseVolumeHandle::surface(int which_side) const {
0216   assert(which_side >= 0 && which_side < 6);
0217   return *(surfaces[which_side]);
0218 }