EcalCommonData/src/DDEcalEndcapTrapX.cc

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