00001 //-------------------------------------------------------------------------- 00002 // 00003 // Environment: 00004 // This software is part of models developed at BES collaboration 00005 // based on the EvtGen framework. If you use all or part 00006 // of it, please give an appropriate acknowledgement. 00007 // 00008 // Copyright Information: See EvtGen/BesCopyright 00009 // Copyright (A) 2006 Ping Rong-Gang @IHEP 00010 // 00011 // Module: EvtDIY.cc 00012 // 00013 // Description: Class to calculate the Euler angles to rotate a system 00014 // 00015 // Modification history: 00016 // 00017 // Ping R.-G. December, 2007 Module created 00018 // 00019 //------------------------------------------------------------------------ 00020 // 00021 #include "EvtGenBase/EvtPatches.hh" 00022 #include <iostream> 00023 #include <math.h> 00024 #include <fstream> 00025 #include <stdio.h> 00026 #include <stdlib.h> 00027 #include "EvtGenBase/EvtEulerAngles.hh" 00028 #include "EvtGenBase/EvtReport.hh" 00029 using namespace std; //::endl; 00030 00031 EvtEulerAngles::~EvtEulerAngles(){} 00032 EvtEulerAngles::EvtEulerAngles(const EvtVector3R & Yaxis, const EvtVector3R & Zaxis){ 00033 _Yaxis=Yaxis; 00034 _Zaxis=Zaxis; 00035 EulerAngles(); 00036 } 00037 00038 EvtEulerAngles::EvtEulerAngles( const EvtVector4R & Pyaxis, const EvtVector4R & Pzaxis){ 00039 for (int i=1;i<4;i++){ 00040 _Yaxis.set(i-1,Pyaxis.get(i)); 00041 _Zaxis.set(i-1,Pzaxis.get(i)); 00042 } 00043 EulerAngles(); 00044 } 00045 00046 EvtEulerAngles::EvtEulerAngles(){ 00047 } 00048 00049 double EvtEulerAngles::getAlpha(){ 00050 return _alpha; 00051 } 00052 00053 double EvtEulerAngles::getBeta(){ 00054 return _beta; 00055 } 00056 00057 double EvtEulerAngles::getGamma(){ 00058 return _gamma; 00059 } 00060 00061 void EvtEulerAngles::EulerAngles(){ 00062 // to calculate Euler angles with y-convention, i.e. R=R(Z, alpha).R(Y,beta).R(Z,gamma) 00063 double pi=3.1415926; 00064 _ry=_Yaxis.d3mag(); 00065 _rz=_Zaxis.d3mag(); 00066 00067 if(_ry==0 ||_rz==0) { 00068 report(ERROR,"") << "Euler angle calculation specified by zero modules of the axis!"<<endl; 00069 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; 00070 ::abort(); 00071 } 00072 double tolerance=1e-10; 00073 bool Y1is0=fabs(_Yaxis.get(0))<tolerance; 00074 bool Y2is0=fabs(_Yaxis.get(1))<tolerance; 00075 bool Y3is0=fabs(_Yaxis.get(2))<tolerance; 00076 bool Z1is0=fabs(_Zaxis.get(0))<tolerance; 00077 bool Z2is0=fabs(_Zaxis.get(1))<tolerance; 00078 bool Z3is0=fabs(_Zaxis.get(2))<tolerance; 00079 00080 if(Y1is0 && Y3is0 && Z1is0 && Z2is0 ){ 00081 _alpha=0; _beta=0; _gamma=0; 00082 return; 00083 } 00084 00085 if( Z1is0 && Z2is0 && !Y2is0){ 00086 _alpha=0; _beta=0; 00087 _gamma=acos(_Yaxis.get(0)/_ry); 00088 if(_Yaxis.get(1)<0) _gamma=2*pi - _gamma; 00089 return; 00090 } 00091 00092 // For general case to calculate Euler angles 00093 // to calculate beta 00094 00095 if(Z1is0 && Z2is0) { 00096 _beta=0; 00097 } else{ _beta =acos(_Zaxis.get(2)/_rz);} 00098 00099 //to calculate alpha 00100 00101 if(_beta==0){ 00102 _alpha=0.0; 00103 }else { 00104 double cosalpha=_Zaxis.get(0)/_rz/sin(_beta); 00105 if(fabs(cosalpha)>1.0) cosalpha=cosalpha/fabs(cosalpha); 00106 _alpha=acos(cosalpha); 00107 if(_Zaxis.get(1)<0.0) _alpha=2*pi - _alpha; 00108 } 00109 00110 //to calculate gamma, alpha=0 and beta=0 case has been calculated, so only alpha !=0 and beta !=0 case left 00111 00112 double singamma=_Yaxis.get(2)/_ry/sin(_beta); 00113 double cosgamma=(-_Yaxis.get(0)/_ry-cos(_alpha)*cos(_beta)*singamma)/sin(_alpha); 00114 if(fabs(singamma)>1.0) singamma=singamma/fabs(singamma); 00115 _gamma=asin(singamma); 00116 if(singamma>0 && cosgamma<0 ) _gamma=pi - _gamma; // _gamma>0 00117 if(singamma<0 && cosgamma<0 ) _gamma=pi - _gamma; //_gamma<0 00118 if(singamma<0 && cosgamma>0 ) _gamma=2*pi + _gamma; //_gamma<0 00119 00120 00121 }