EcalCommonData/src/DDEcalEndcapTrap.cc

0001 #include "Geometry/EcalCommonData/interface/DDEcalEndcapTrap.h"
0002 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0003 #include <CLHEP/Units/SystemOfUnits.h>
0004 #include <CLHEP/Geometry/Point3D.h>
0005 #include <CLHEP/Geometry/Plane3D.h>
0006 #include <CLHEP/Geometry/Vector3D.h>
0007 #include <CLHEP/Geometry/Transform3D.h>
0008 #include <CLHEP/Vector/EulerAngles.h>
0009
0010 //#define EDM_ML_DEBUG
0011
0012 // Implementation of DDEcalEndcapTrap class
0013
0014 DDEcalEndcapTrap::DDEcalEndcapTrap(const int hand, const double front, const double rear, const double length) {
0015   //
0016   //  Initialise corners of supercrystal.
0017
0018   // Start out with bottom surface on (x,z) plane, front face in (x,y) plane.
0019
0020   double xsign;
0021
0022   if (hand == 2) {
0023     xsign = -1.;
0024   } else {
0025     xsign = 1.;
0026   }
0027
0028   m_hand = hand;
0029   m_front = front;
0030   m_rear = rear;
0031   m_length = length;
0032
0033   int icorner;
0034   icorner = 1;
0035   m_corners[3 * icorner - 3] = xsign * front;
0036   m_corners[3 * icorner - 2] = front;
0037   m_corners[3 * icorner - 1] = 0.;
0038   icorner = 2;
0039   m_corners[3 * icorner - 3] = xsign * front;
0040   m_corners[3 * icorner - 2] = 0.;
0041   m_corners[3 * icorner - 1] = 0.;
0042   icorner = 3;
0043   m_corners[3 * icorner - 3] = 0.;
0044   m_corners[3 * icorner - 2] = 0.;
0045   m_corners[3 * icorner - 1] = 0.;
0046   icorner = 4;
0047   m_corners[3 * icorner - 3] = 0.;
0048   m_corners[3 * icorner - 2] = front;
0049   m_corners[3 * icorner - 1] = 0.;
0050
0051   icorner = 5;
0052   m_corners[3 * icorner - 3] = xsign * rear;
0053   m_corners[3 * icorner - 2] = rear;
0054   m_corners[3 * icorner - 1] = length;
0055   icorner = 6;
0056   m_corners[3 * icorner - 3] = xsign * rear;
0057   m_corners[3 * icorner - 2] = 0.;
0058   m_corners[3 * icorner - 1] = length;
0059   icorner = 7;
0060   m_corners[3 * icorner - 3] = 0.;
0061   m_corners[3 * icorner - 2] = 0.;
0062   m_corners[3 * icorner - 1] = length;
0063   icorner = 8;
0064   m_corners[3 * icorner - 3] = 0.;
0065   m_corners[3 * icorner - 2] = rear;
0066   m_corners[3 * icorner - 1] = length;
0067
0068   calculateCentres();
0069
0070   // Move centre of SC to (0,0,0)
0071
0072   translate();
0073
0074   // Rotate into standard position (face centres on z axis)
0075
0076   //  this->rotate();
0077
0078   calculateCentres();
0079 }
0080
0081 //void DDEcalEndcapTrap::rotate() {
0082 //  //
0083 //  //  Rotate supercrystal to standard position
0084 //  //
0085 //  edm::LogaVerbatim("EcalGeom") << "DDEcalEndcapTrap::rotate() - not yet implemented";
0086 //}
0087
0088 void DDEcalEndcapTrap::rotate(const DDTranslation& frontCentre, const DDTranslation& rearCentre) {
0089   //
0090   //  Rotate supercrystal to bring front and rear face centres to specified points
0091   //
0092   edm::LogWarning("EcalGeom") << "DDEcalEndcapTrap::rotate(DDTranslation,DDTranslation) - not yet implemented";
0093 }
0094
0095 void DDEcalEndcapTrap::rotate(const DDRotationMatrix& rot) {
0096   //
0097   //  Rotate supercrystal by specified rotation about (0,0,0)
0098   //
0099
0100   int icorner;
0101   DDTranslation cc;
0102 #ifdef EDM_ML_DEBUG
0103   edm::LogVerbatim("EcalGeom") << "DDEcalEndcapTrap::rotate - rotation " << rot;
0104 #endif
0105   for (icorner = 1; icorner <= 8; icorner++) {
0106     cc = cornerPos(icorner);
0107 #ifdef EDM_ML_DEBUG
0108     edm::LogVerbatim("EcalGeom") << "   Corner (orig) " << icorner << cc;
0109 #endif
0110     cc = rot * cc;
0111 #ifdef EDM_ML_DEBUG
0112     edm::LogVerbatim("EcalGeom") << "   Corner (rot)  " << icorner << cc;
0113 #endif
0114     cornerPos(icorner, cc);
0115   }
0116   m_rotation = rot * m_rotation;
0117   calculateCentres();
0118 }
0119
0120 void DDEcalEndcapTrap::translate() {
0121 #ifdef EDM_ML_DEBUG
0122   edm::LogVerbatim("EcalGeom") << "DDEcalEndcapTrap::translate() not yet implemented";
0123 #endif
0124   translate(-1. * centrePos());
0125 }
0126
0127 void DDEcalEndcapTrap::translate(const DDTranslation& trans) {
0128   //
0129   //  Translate supercrystal by specified amount
0130   //
0131
0132   DDTranslation tcorner;
0133   for (int icorner = 1; icorner <= 8; icorner++) {
0134     tcorner = cornerPos(icorner) + trans;
0135     cornerPos(icorner, tcorner);
0136   }
0137   calculateCentres();
0138   m_translation = trans + m_translation;
0139 }
0140
0141 void DDEcalEndcapTrap::moveto(const DDTranslation& frontCentre, const DDTranslation& rearCentre) {
0142   //
0143   //  Rotate (about X then about Y) and translate supercrystal to bring axis joining front and rear face centres parallel to line connecting specified points
0144   //
0145
0146   //  Get azimuthal and polar angles of current axis and target axis
0147   double currentTheta = elevationAngle();
0148   double currentPhi = polarAngle();
0149   double targetTheta = elevationAngle(frontCentre - rearCentre);
0150   double targetPhi = polarAngle(frontCentre - rearCentre);
0151
0152   //  Rotate to correct angle (X then Y)
0153 #ifdef EDM_ML_DEBUG
0154   edm::LogVerbatim("EcalGeom") << "moveto: frontCentre " << frontCentre << std::endl
0155                                << "moveto: rearCentre  " << rearCentre << std::endl
0156                                << "moveto: X rotation: " << targetTheta << " " << currentTheta << " "
0157                                << targetTheta - currentTheta << std::endl
0158                                << "moveto: Y rotation: " << targetPhi << " " << currentPhi << " "
0159                                << " " << targetPhi - currentPhi;
0160 #endif
0161   rotateX(targetTheta - currentTheta);
0162   rotateY(targetPhi - currentPhi);
0163
0164   //  Translate SC to final position
0165   DDTranslation targetCentre = 0.5 * (frontCentre + rearCentre);
0166 #ifdef EDM_ML_DEBUG
0167   edm::LogVerbatim("EcalGeom") << "moveto: translation " << targetCentre - centrePos();
0168 #endif
0169   translate(targetCentre - centrePos());
0170 }
0171
0172 void DDEcalEndcapTrap::rotateX(const double angle) {
0173   //
0174   //  Rotate SC through given angle about X axis
0175   //
0176
0177   const CLHEP::HepRotation tmp(CLHEP::Hep3Vector(1., 0., 0.), angle);
0178
0179   rotate(DDRotationMatrix(tmp.xx(), tmp.xy(), tmp.xz(), tmp.yx(), tmp.yy(), tmp.yz(), tmp.zx(), tmp.zy(), tmp.zz()));
0180 }
0181
0182 void DDEcalEndcapTrap::rotateY(const double angle) {
0183   //
0184   //  Rotate SC through given angle about Y axis
0185   //
0186   const CLHEP::HepRotation tmp(CLHEP::Hep3Vector(0., 1., 0.), angle);
0187
0188   rotate(DDRotationMatrix(tmp.xx(), tmp.xy(), tmp.xz(), tmp.yx(), tmp.yy(), tmp.yz(), tmp.zx(), tmp.zy(), tmp.zz()));
0189 }
0190
0191 void DDEcalEndcapTrap::calculateCentres() {
0192   //
0193   //  Calculate crystal centre and front & rear face centres
0194   //
0195
0196   int ixyz, icorner;
0197
0198   for (ixyz = 0; ixyz < 3; ixyz++) {
0199     m_centre[ixyz] = 0;
0200     m_fcentre[ixyz] = 0;
0201     m_rcentre[ixyz] = 0;
0202   }
0203
0204   for (icorner = 1; icorner <= 4; icorner++) {
0205     for (ixyz = 0; ixyz < 3; ixyz++) {
0206       m_centre[ixyz] = m_centre[ixyz] + 0.125 * m_corners[3 * icorner - 3 + ixyz];
0207       m_fcentre[ixyz] = m_fcentre[ixyz] + 0.25 * m_corners[3 * icorner - 3 + ixyz];
0208     }
0209   }
0210   for (icorner = 5; icorner <= 8; icorner++) {
0211     for (ixyz = 0; ixyz < 3; ixyz++) {
0212       m_centre[ixyz] = m_centre[ixyz] + 0.125 * m_corners[3 * icorner - 3 + ixyz];
0213       m_rcentre[ixyz] = m_rcentre[ixyz] + 0.25 * m_corners[3 * icorner - 3 + ixyz];
0214     }
0215   }
0216 }
0217
0218 DDTranslation DDEcalEndcapTrap::cornerPos(const int icorner) {
0219   //
0220   //  Return specified corner as a DDTranslation
0221   //
0222   return DDTranslation(m_corners[3 * icorner - 3], m_corners[3 * icorner - 2], m_corners[3 * icorner - 1]);
0223 }
0224
0225 void DDEcalEndcapTrap::cornerPos(const int icorner, const DDTranslation& cornerxyz) {
0226   //
0227   //  Save position of specified corner.
0228   //
0229   for (int ixyz = 0; ixyz < 3; ixyz++) {
0230     m_corners[3 * icorner - 3 + ixyz] = (0 == ixyz ? cornerxyz.x() : (1 == ixyz ? cornerxyz.y() : cornerxyz.z()));
0231     ;
0232   }
0233 }
0234
0235 DDTranslation DDEcalEndcapTrap::centrePos() {
0236   //
0237   //  Return SC centre as a DDTranslation
0238   //
0239   return DDTranslation(m_centre[0], m_centre[1], m_centre[2]);
0240 }
0241
0242 DDTranslation DDEcalEndcapTrap::fcentrePos() {
0243   //
0244   //  Return SC front face centre as a DDTranslation
0245   //
0246   return DDTranslation(m_fcentre[0], m_fcentre[1], m_fcentre[2]);
0247 }
0248
0249 DDTranslation DDEcalEndcapTrap::rcentrePos() {
0250   //
0251   //  Return SC rear face centre as a DDTranslation
0252   //
0253   return DDTranslation(m_rcentre[0], m_rcentre[1], m_rcentre[2]);
0254 }
0255
0256 double DDEcalEndcapTrap::elevationAngle(const DDTranslation& trans) {
0257   //
0258   //  Return elevation angle (out of x-z plane) of a given translation (seen as a vector from the origin).
0259   //
0260   double sintheta = trans.y() / trans.r();
0261   return asin(sintheta);
0262 }
0263
0264 double DDEcalEndcapTrap::elevationAngle() {
0265   //
0266   //  Return elevation angle (out of x-z plane) of SC in current position.
0267   //
0268   DDTranslation current = fcentrePos() - rcentrePos();
0269   return elevationAngle(current);
0270 }
0271
0272 double DDEcalEndcapTrap::polarAngle(const DDTranslation& trans) {
0273   //
0274   //  Return polar angle (from x to z) of a given translation (seen as a vector from the origin).
0275   //
0276   double tanphi = trans.x() / trans.z();
0277   return atan(tanphi);
0278 }
0279
0280 double DDEcalEndcapTrap::polarAngle() {
0281   //
0282   //  Return elevation angle (out of x-z plane) of SC in current position.
0283   //
0284   DDTranslation current = fcentrePos() - rcentrePos();
0285   return polarAngle(current);
0286 }
0287
0288 void DDEcalEndcapTrap::print() {
0289   //
0290   //  Print SC coordinates for debugging
0291   //
0292   edm::LogVerbatim("EcalGeom") << "Endcap supercrystal";
0293   for (int ic = 1; ic <= 8; ic++) {
0294     DDTranslation cc = cornerPos(ic);
0295     edm::LogVerbatim("EcalGeom") << "Corner " << ic << " " << cc;
0296   }
0297   edm::LogVerbatim("EcalGeom") << "    Centre " << centrePos() << std::endl
0298                                << "   fCentre " << fcentrePos() << std::endl
0299                                << "   rCentre " << rcentrePos();
0300 }