00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <math.h>
00023 #include "EvtGenBase/EvtKine.hh"
00024 #include "EvtGenBase/EvtConst.hh"
00025 #include "EvtGenBase/EvtVector4R.hh"
00026 #include "EvtGenBase/EvtVector4C.hh"
00027 #include "EvtGenBase/EvtTensor4C.hh"
00028 #include "EvtGenBase/EvtdFunction.hh"
00029 #include "EvtGenBase/EvtReport.hh"
00030
00031
00032
00033 double EvtDecayAngle(const EvtVector4R& p,const EvtVector4R& q,
00034 const EvtVector4R& d) {
00035
00036 double pd=p*d;
00037 double pq=p*q;
00038 double qd=q*d;
00039 double mp2=p.mass2();
00040 double mq2=q.mass2();
00041 double md2=d.mass2();
00042
00043 double cost=(pd*mq2-pq*qd)/sqrt((pq*pq-mq2*mp2)*(qd*qd-mq2*md2));
00044
00045 return cost;
00046
00047 }
00048
00049 double EvtDecayAngleChi(const EvtVector4R& p4_p,const EvtVector4R& p4_d1,
00050 const EvtVector4R& p4_d2,const EvtVector4R& p4_h1,
00051 const EvtVector4R& p4_h2 ) {
00052
00053 EvtVector4R p4_d1p,p4_h1p,p4_h2p,p4_d2p;
00054
00055
00056
00057
00058 p4_d1p=boostTo(p4_d1,p4_p);
00059 p4_d2p=boostTo(p4_d2,p4_p);
00060 p4_h1p=boostTo(p4_h1,p4_p);
00061 p4_h2p=boostTo(p4_h2,p4_p);
00062
00063
00064 EvtVector4R d1_perp,d1_prime,h1_perp;
00065 EvtVector4R D;
00066
00067 D=p4_d1p+p4_d2p;
00068
00069 d1_perp=p4_d1p-(D.dot(p4_d1p)/D.dot(D))*D;
00070 h1_perp=p4_h1p-(D.dot(p4_h1p)/D.dot(D))*D;
00071
00072
00073
00074 d1_prime=D.cross(d1_perp);
00075
00076 d1_perp= d1_perp/d1_perp.d3mag();
00077 d1_prime= d1_prime/d1_prime.d3mag();
00078
00079 double x,y;
00080
00081 x=d1_perp.dot(h1_perp);
00082 y=d1_prime.dot(h1_perp);
00083
00084 double chi=atan2(y,x);
00085
00086 if (chi<0.0) chi+=EvtConst::twoPi;
00087
00088 return chi;
00089
00090 }
00091
00092
00093
00094 double EvtDecayPlaneNormalAngle(const EvtVector4R& p,const EvtVector4R& q,
00095 const EvtVector4R& d1,const EvtVector4R& d2){
00096
00097 EvtVector4C lc=dual(directProd(d1,d2)).cont2(q);
00098
00099 EvtVector4R l(real(lc.get(0)),real(lc.get(1)),
00100 real(lc.get(2)),real(lc.get(3)));
00101
00102 double pq=p*q;
00103
00104 return q.mass()*(p*l)/sqrt(-(pq*pq-p.mass2()*q.mass2())*l.mass2());
00105
00106
00107 }
00108
00109
00110
00111 double EvtDecayAnglePhi( const EvtVector4R& z, const EvtVector4R& p, const
00112 EvtVector4R& q, const EvtVector4R& d )
00113 {
00114 double eq = (p * q) / p.mass();
00115 double ed = (p * d) / p.mass();
00116 double mq = q.mass();
00117 double q2 = p.mag2r3(q);
00118 double qd = p.dotr3(q,d);
00119 double zq = p.dotr3(z,q);
00120 double zd = p.dotr3(z,d);
00121 double alpha = (eq - mq)/(q2 * mq) * qd - ed/mq;
00122
00123 double y = p.scalartripler3(z,q,d) + alpha * p.scalartripler3(z,q,q);
00124 double x = (zq * (qd + alpha * q2) - q2 * (zd + alpha * zq)) / sqrt(q2);
00125
00126 double phi = atan2(y,x);
00127
00128 return phi<0 ? (phi+EvtConst::twoPi) : phi;
00129 }
00130
00131
00132 EvtComplex wignerD( int j, int m1, int m2, double phi, double theta, double
00133 gamma )
00134 {
00135 EvtComplex gp(0.0, -phi*m1);
00136 EvtComplex gm(0.0, -gamma*m2);
00137
00138
00139 return exp( gp ) * EvtdFunction::d(2*j, 2*m1,2*m2, theta) * exp( gm );
00140 }
00141
00142