Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:29:12

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