00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "MdcRecoUtil/DifFourVector.h"
00023 using std::endl;
00024 using std::cout;
00025
00026
00027 DifFourVector::DifFourVector()
00028 :E(0.0),P(0.0,0.0,0.0)
00029 {}
00030
00031 DifFourVector::DifFourVector
00032 (const DifNumber& m,const DifVector& p)
00033 :E(sqrt(m*m+p*p)),P(p)
00034 {}
00035 DifFourVector::DifFourVector
00036 (const double& m,const DifVector& p)
00037 :E(sqrt(m*m+p*p)),P(p)
00038 {}
00039
00040 DifFourVector::DifFourVector(const DifFourVector& v)
00041 :E(v.E),P(v.P)
00042 {}
00043
00044 HepSymMatrix DifFourVector::errorMatrix(const HepSymMatrix& e)const {
00045 HepSymMatrix temp(4);
00046 temp(1,1)=correlation(E,E,e);
00047 temp(1,2)=correlation(E,P.x,e);
00048 temp(1,3)=correlation(E,P.y,e);
00049 temp(1,4)=correlation(E,P.z,e);
00050 temp(2,2)=correlation(P.x,P.x,e);
00051 temp(2,3)=correlation(P.x,P.y,e);
00052 temp(2,4)=correlation(P.x,P.z,e);
00053 temp(3,3)=correlation(P.y,P.y,e);
00054 temp(3,4)=correlation(P.y,P.z,e);
00055 temp(4,4)=correlation(P.z,P.z,e);
00056 return temp;
00057 }
00058
00059 HepMatrix DifFourVector::jacobian()const{
00060 int npar=E.nPar();
00061 HepMatrix temp(4,npar);
00062 for(int i=1; i<=npar; i++){
00063 temp(1,i)=E.derivative(i);
00064 temp(2,i)=P.x.derivative(i);
00065 temp(3,i)=P.y.derivative(i);
00066 temp(4,i)=P.z.derivative(i);
00067 }
00068 return temp;
00069 }
00070
00071 void
00072 DifFourVector::boostTo(const DifFourVector& pTo)
00073 {
00074 const DifVector xHat(1,0,0);
00075 const DifVector yHat(0,1,0);
00076 const DifVector zHat(0,1,0);
00077 DifVector z=pTo.direction();
00078 DifVector y=zHat-z*(zHat*z);
00079 if(y.length()<0.01) y=xHat-z*(xHat*z);
00080 y.normalize();
00081 DifVector x(cross(y,z));
00082
00083 DifNumber px=P*x;
00084 DifNumber py=P*y;
00085 DifNumber pz=P*z;
00086
00087 DifNumber gamma=pTo.E/pTo.mass();
00088 DifNumber beta=pTo.pMag()/pTo.E;
00089
00090 DifNumber pzP=gamma*(pz-beta*E);
00091 DifNumber eP=gamma*(E-beta*pz);
00092
00093 E=eP;
00094 P=px*x+py*y+pzP*z;
00095
00096 return;
00097
00098 }
00099
00100 void
00101 DifFourVector::boostToMe(std::vector<DifFourVector*>& list)const{
00102 const DifVector xHat(1,0,0);
00103 const DifVector yHat(0,1,0);
00104 const DifVector zHat(0,0,1);
00105 DifVector z=P;
00106 z.normalize();
00107 DifVector y(zHat-z*(zHat*z));
00108 if(y.lengthSq()<0.0001) y=xHat-z*(xHat*z);
00109 y.normalize();
00110 DifVector x(cross(y,z));
00111
00112 DifNumber gamma=E/mass();
00113 DifNumber beta=pMag()/E;
00114
00115 for(int i=0;i<(int)list.size();i++) {
00116 DifFourVector& p4i=*list[i];
00117 DifNumber px=p4i.P*x;
00118 DifNumber py=p4i.P*y;
00119 DifNumber pz=p4i.P*z;
00120 DifNumber e=p4i.E;
00121
00122 DifNumber pzP=gamma*(pz-beta*e);
00123 DifNumber eP=gamma*(e-beta*pz);
00124
00125 p4i.E=eP;
00126 p4i.P=px*x+py*y+pzP*z;
00127
00128 }
00129
00130 }
00131
00132
00133 void
00134 DifFourVector::boostFromMe(std::vector<DifFourVector*>& list)const{
00135 const DifVector xHat(1,0,0);
00136 const DifVector yHat(0,1,0);
00137 const DifVector zHat(0,0,1);
00138 DifVector z=P;
00139 z.normalize();
00140 DifVector y(zHat-z*(zHat*z));
00141 if(y.lengthSq()<0.0001) y=xHat-z*(xHat*z);
00142 y.normalize();
00143 DifVector x(cross(y,z));
00144
00145 DifNumber gamma=E/mass();
00146 DifNumber beta=pMag()/E;
00147
00148 for(int i=0;i<(int)list.size();i++) {
00149 DifFourVector& p4i=*list[i];
00150 DifNumber px=p4i.P*x;
00151 DifNumber py=p4i.P*y;
00152 DifNumber pz=p4i.P*z;
00153 DifNumber e=p4i.E;
00154
00155 DifNumber pzP=gamma*(pz+beta*e);
00156 DifNumber eP=gamma*(e+beta*pz);
00157
00158 p4i.E=eP;
00159 p4i.P=px*x+py*y+pzP*z;
00160
00161 }
00162
00163 }
00164
00165 void
00166 DifFourVector::boostFrom(const DifFourVector& pFrom)
00167 {
00168 const DifVector xHat(1,0,0);
00169 const DifVector yHat(0,1,0);
00170 const DifVector zHat(0,1,0);
00171 DifVector z=pFrom.direction();
00172 DifVector y=zHat-z*(zHat*z);
00173 if(y.length()<0.01) y=xHat-z*(xHat*z);
00174 y.normalize();
00175 DifVector x(cross(y,z));
00176
00177 DifNumber px=P*x;
00178 DifNumber py=P*y;
00179 DifNumber pz=P*z;
00180
00181 DifNumber gamma=pFrom.E/pFrom.mass();
00182 DifNumber beta=pFrom.pMag()/pFrom.E;
00183
00184 DifNumber pzP=gamma*(pz+beta*E);
00185 DifNumber eP=gamma*(E+beta*pz);
00186
00187 E=eP;
00188 P=px*x+py*y+pzP*z;
00189 }
00190
00191
00192 void DifFourVector::print()const {
00193
00194
00195 cout << "SKIP of DifFourVector::print()" <<endl;
00196 }
00197