EcalTestBeam/test/testFromCMStoTB.cpp

0001 #include <CLHEP/Random/RandFlat.h>
0002 #include <CLHEP/Units/SystemOfUnits.h>
0003 #include "DataFormats/Math/interface/Point3D.h"
0004 #include "Math/GenVector/Rotation3D.h"
0005 #include "Math/GenVector/RotationX.h"
0006 #include "Math/GenVector/RotationY.h"
0007 #include "Math/GenVector/RotationZ.h"
0008
0009 #include <cmath>
0010 #include <vector>
0011
0012 using namespace std;
0013 using CLHEP::deg;
0014
0015 void computeRotation(double &myTheta, double &myPhi, ROOT::Math::Rotation3D &CMStoTB, ROOT::Math::Rotation3D &TBtoCMS) {
0016   // rotation matrix to move from the CMS reference frame to the test beam one
0017
0018   ROOT::Math::Rotation3D *fromCMStoTB = new ROOT::Math::Rotation3D();
0019
0020   // rotation matrix to move from the test beam reference frame to the CMS one
0021
0022   ROOT::Math::Rotation3D *fromTBtoCMS = new ROOT::Math::Rotation3D();
0023
0024   double angle1 = 90. * deg - myPhi;
0025   ROOT::Math::RotationZ *r1 = new ROOT::Math::RotationZ(angle1);
0026   double angle2 = myTheta;
0027   ROOT::Math::RotationX *r2 = new ROOT::Math::RotationX(angle2);
0028   double angle3 = 90. * deg;
0029   ROOT::Math::RotationZ *r3 = new ROOT::Math::RotationZ(angle3);
0030   (*fromCMStoTB) *= (*r3);
0031   (*fromCMStoTB) *= (*r2);
0032   (*fromCMStoTB) *= (*r1);
0033
0034   cout << "Rotation matrix from CMS to test beam frame = " << (*fromCMStoTB) << "built from: \n"
0035        << " the rotation of " << angle1 << " around Z " << (*r1) << "\n"
0036        << " the rotation of " << angle2 << " around X " << (*r2) << "\n"
0037        << " the rotation of " << angle3 << " around Z " << (*r3) << std::endl;
0038
0039   double angle11 = -90. * deg + myPhi;
0040   ROOT::Math::RotationZ *r11 = new ROOT::Math::RotationZ(angle11);
0041   double angle12 = -myTheta;
0042   ROOT::Math::RotationX *r12 = new ROOT::Math::RotationX(angle12);
0043   double angle13 = -90. * deg;
0044   ROOT::Math::RotationZ *r13 = new ROOT::Math::RotationZ(angle13);
0045
0046   (*fromTBtoCMS) *= (*r11);
0047   (*fromTBtoCMS) *= (*r12);
0048   (*fromTBtoCMS) *= (*r13);
0049
0050   cout << "Rotation matrix from test beam to CMS frame = " << (*fromTBtoCMS) << "built from: \n"
0051        << " the rotation of " << angle13 << " around Z " << (*r13) << "\n"
0052        << " the rotation of " << angle12 << " around X " << (*r12) << "\n"
0053        << " the rotation of " << angle11 << " around Z " << (*r11) << std::endl;
0054
0055   ROOT::Math::Rotation3D test = (*fromCMStoTB) * (*fromTBtoCMS);
0056
0057   cout << "Product of the two rotations: " << test << endl;
0058
0059   CMStoTB = (*fromCMStoTB);
0060   TBtoCMS = (*fromTBtoCMS);
0061
0062   delete fromCMStoTB;
0063   delete fromTBtoCMS;
0064   delete r1;
0065   delete r2;
0066   delete r3;
0067 }
0068
0069 void checkTotalRotation(double &myTheta, double &myPhi) {
0070   // rotation matrix to move from the CMS reference frame to the test beam one
0071
0072   double xx = -cos(myTheta) * cos(myPhi);
0073   double xy = -cos(myTheta) * sin(myPhi);
0074   double xz = sin(myTheta);
0075
0076   double yx = sin(myPhi);
0077   double yy = -cos(myPhi);
0078   double yz = 0.;
0079
0080   double zx = sin(myTheta) * cos(myPhi);
0081   double zy = sin(myTheta) * sin(myPhi);
0082   double zz = cos(myTheta);
0083
0084   ROOT::Math::Rotation3D *fromCMStoTB = new ROOT::Math::Rotation3D(xx, xy, xz, yx, yy, yz, zx, zy, zz);
0085
0086   cout << "Total rotation matrix from CMS to test beam frame = " << (*fromCMStoTB) << endl;
0087
0088   // rotation matrix to move from the test beam reference frame to the CMS one
0089
0090   xx = -cos(myTheta) * cos(myPhi);
0091   xy = sin(myPhi);
0092   xz = sin(myTheta) * cos(myPhi);
0093
0094   yx = -cos(myTheta) * sin(myPhi);
0095   yy = -cos(myPhi);
0096   yz = sin(myTheta) * sin(myPhi);
0097
0098   zx = sin(myTheta);
0099   zy = 0.;
0100   zz = cos(myTheta);
0101
0102   ROOT::Math::Rotation3D *fromTBtoCMS = new ROOT::Math::Rotation3D(xx, xy, xz, yx, yy, yz, zx, zy, zz);
0103
0104   cout << "Total rotation matrix from test beam to CMS frame = " << (*fromTBtoCMS) << endl;
0105
0106   ROOT::Math::Rotation3D test = (*fromCMStoTB) * (*fromTBtoCMS);
0107
0108   cout << "Product of the two rotations: " << test << endl;
0109 }
0110
0111 int main() {
0112   // (eta,phi) for a crystal
0113
0114   double myMod = 1.;
0115   double myEta = 0.971226;
0116   double myTheta = 2.0 * atan(exp(-myEta));
0117   double myPhi = 0.115052;
0118
0119   cout << "\n===========================================\n" << endl;
0120   cout << "Input theta = " << myTheta << " phi = " << myPhi << endl;
0121   cout << "\n===========================================\n" << endl;
0122
0123   ROOT::Math::Rotation3D *CMStoTB = new ROOT::Math::Rotation3D();
0124   ROOT::Math::Rotation3D *TBtoCMS = new ROOT::Math::Rotation3D();
0125
0126   checkTotalRotation(myTheta, myPhi);
0127
0128   computeRotation(myTheta, myPhi, (*CMStoTB), (*TBtoCMS));
0129   cout << "\n===========================================\n" << endl;
0130
0131   double xx = 0.;
0132   double yy = 0.;
0133   double zz = 0.;
0134   math::XYZPoint test(xx, yy, zz);
0135   double newTheta = myTheta;
0136   double newPhi = myPhi;
0137
0138   for (int ieta = -1; ieta <= 1; ++ieta) {
0139     for (int iphi = -1; iphi <= 1; ++iphi) {
0140       newTheta = myTheta + (double)ieta;
0141       newPhi = myPhi + (double)iphi;
0142
0143       xx = myMod * sin(newTheta) * cos(newPhi);
0144       yy = myMod * sin(newTheta) * sin(newPhi);
0145       zz = myMod * cos(newTheta);
0146       test.SetCoordinates(xx, yy, zz);
0147
0148       cout << "\n ieta = " << ieta << " iphi = " << iphi << endl;
0149       cout << "\n From CMS to TB \n" << endl;
0150       cout << "Input vector  = " << test << " corresponding to theta = " << newTheta << " phi = " << newPhi << endl;
0151
0152       math::XYZPoint testrot = (*CMStoTB) * test;
0153
0154       cout << "Output vector = " << testrot << " corresponding to theta = " << testrot.theta()
0155            << " phi = " << testrot.phi() << endl;
0156
0157       cout << "\n From TB to CMS \n" << endl;
0158
0159       math::XYZPoint thistest = (*TBtoCMS) * testrot;
0160
0161       cout << "Output vector = " << thistest << " corresponding to theta = " << thistest.theta()
0162            << " phi = " << thistest.phi() << endl;
0163       cout << "\n===========================================\n" << endl;
0164     }
0165   }
0166
0167   cout << "\n===========================================\n" << endl;
0168
0169   for (int ix = -1; ix <= 1; ++ix) {
0170     for (int iy = -1; iy <= 1; ++iy) {
0171       xx = (double)ix * 0.01;
0172       yy = (double)iy * 0.01;
0173       zz = 1.;
0174       test.SetCoordinates(xx, yy, zz);
0175
0176       cout << "\n ix = " << ix << " iy = " << iy << endl;
0177       cout << "\n From TB to CMS \n" << endl;
0178       cout << "Input vector  = " << test << endl;
0179
0180       math::XYZPoint testrot = (*TBtoCMS) * test;
0181
0182       cout << "Output vector = " << testrot << " corresponding to theta = " << testrot.theta()
0183            << " phi = " << testrot.phi() << endl;
0184
0185       cout << "\n From CMS to TB \n" << endl;
0186
0187       math::XYZPoint thistest = (*CMStoTB) * testrot;
0188
0189       cout << "Output vector = " << thistest << endl;
0190       cout << "\n===========================================\n" << endl;
0191     }
0192   }
0193
0194   return 0;
0195 }