1 //--------------------------------------------------------------------------
4 // This software is part of the EvtGen package developed jointly
5 // for the BaBar and CLEO collaborations. If you use all or part
6 // of it, please give an appropriate acknowledgement.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
13 // Description: routines to calculate decay angles.
15 // Modification history:
17 // DJL/RYD September 25, 1996 Module created
19 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenBase/EvtKine.hh"
24 #include "EvtGenBase/EvtConst.hh"
25 #include "EvtGenBase/EvtVector4R.hh"
26 #include "EvtGenBase/EvtVector4C.hh"
27 #include "EvtGenBase/EvtTensor4C.hh"
28 #include "EvtGenBase/EvtdFunction.hh"
29 #include "EvtGenBase/EvtReport.hh"
33 double EvtDecayAngle(const EvtVector4R& p,const EvtVector4R& q,
34 const EvtVector4R& d) {
43 double cost=(pd*mq2-pq*qd)/sqrt((pq*pq-mq2*mp2)*(qd*qd-mq2*md2));
49 double EvtDecayAngleChi(const EvtVector4R& p4_p,const EvtVector4R& p4_d1,
50 const EvtVector4R& p4_d2,const EvtVector4R& p4_h1,
51 const EvtVector4R& p4_h2 ) {
53 EvtVector4R p4_d1p,p4_h1p,p4_h2p,p4_d2p;
56 // boost all vectors parent restframe
58 p4_d1p=boostTo(p4_d1,p4_p);
59 p4_d2p=boostTo(p4_d2,p4_p);
60 p4_h1p=boostTo(p4_h1,p4_p);
61 p4_h2p=boostTo(p4_h2,p4_p);
64 EvtVector4R d1_perp,d1_prime,h1_perp;
69 d1_perp=p4_d1p-(D.dot(p4_d1p)/D.dot(D))*D;
70 h1_perp=p4_h1p-(D.dot(p4_h1p)/D.dot(D))*D;
72 // orthogonal to both D and d1_perp
74 d1_prime=D.cross(d1_perp);
76 d1_perp= d1_perp/d1_perp.d3mag();
77 d1_prime= d1_prime/d1_prime.d3mag();
81 x=d1_perp.dot(h1_perp);
82 y=d1_prime.dot(h1_perp);
84 double chi=atan2(y,x);
86 if (chi<0.0) chi+=EvtConst::twoPi;
94 double EvtDecayPlaneNormalAngle(const EvtVector4R& p,const EvtVector4R& q,
95 const EvtVector4R& d1,const EvtVector4R& d2){
97 EvtVector4C lc=dual(directProd(d1,d2)).cont2(q);
99 EvtVector4R l(real(lc.get(0)),real(lc.get(1)),
100 real(lc.get(2)),real(lc.get(3)));
104 return q.mass()*(p*l)/sqrt(-(pq*pq-p.mass2()*q.mass2())*l.mass2());
110 // Calculate phi using the given 4 vectors (all in the same frame)
111 double EvtDecayAnglePhi( const EvtVector4R& z, const EvtVector4R& p, const
112 EvtVector4R& q, const EvtVector4R& d )
114 double eq = (p * q) / p.mass();
115 double ed = (p * d) / p.mass();
116 double mq = q.mass();
117 double q2 = p.mag2r3(q);
118 double qd = p.dotr3(q,d);
119 double zq = p.dotr3(z,q);
120 double zd = p.dotr3(z,d);
121 double alpha = (eq - mq)/(q2 * mq) * qd - ed/mq;
123 double y = p.scalartripler3(z,q,d) + alpha * p.scalartripler3(z,q,q);
124 double x = (zq * (qd + alpha * q2) - q2 * (zd + alpha * zq)) / sqrt(q2);
126 double phi = atan2(y,x);
128 return phi<0 ? (phi+EvtConst::twoPi) : phi;
131 EvtComplex wignerD( int j, int m1, int m2, double phi,
132 double theta, double gamma )
135 EvtComplex gp(0.0, -phi*m1);
136 EvtComplex gm(0.0, -gamma*m2);
138 return exp( gp ) * EvtdFunction::d(j, m1, m2, theta) * exp( gm );