]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //-------------------------------------------------------------------------- |
2 | // | |
3 | // Environment: | |
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. | |
7 | // | |
8 | // Copyright Information: See EvtGen/COPYRIGHT | |
9 | // Copyright (C) 1998 Caltech, UCSB | |
10 | // | |
11 | // Module: EvtKine.cc | |
12 | // | |
13 | // Description: routines to calculate decay angles. | |
14 | // | |
15 | // Modification history: | |
16 | // | |
17 | // DJL/RYD September 25, 1996 Module created | |
18 | // | |
19 | //------------------------------------------------------------------------ | |
20 | // | |
21 | #include "EvtGenBase/EvtPatches.hh" | |
22 | #include <math.h> | |
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" | |
30 | ||
31 | ||
32 | ||
33 | double EvtDecayAngle(const EvtVector4R& p,const EvtVector4R& q, | |
34 | const EvtVector4R& d) { | |
35 | ||
36 | double pd=p*d; | |
37 | double pq=p*q; | |
38 | double qd=q*d; | |
39 | double mp2=p.mass2(); | |
40 | double mq2=q.mass2(); | |
41 | double md2=d.mass2(); | |
42 | ||
43 | double cost=(pd*mq2-pq*qd)/sqrt((pq*pq-mq2*mp2)*(qd*qd-mq2*md2)); | |
44 | ||
45 | return cost; | |
46 | ||
47 | } | |
48 | ||
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 ) { | |
52 | ||
53 | EvtVector4R p4_d1p,p4_h1p,p4_h2p,p4_d2p; | |
54 | ||
55 | ||
56 | // boost all vectors parent restframe | |
57 | ||
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); | |
62 | ||
63 | ||
64 | EvtVector4R d1_perp,d1_prime,h1_perp; | |
65 | EvtVector4R D; | |
66 | ||
67 | D=p4_d1p+p4_d2p; | |
68 | ||
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; | |
71 | ||
72 | // orthogonal to both D and d1_perp | |
73 | ||
74 | d1_prime=D.cross(d1_perp); | |
75 | ||
76 | d1_perp= d1_perp/d1_perp.d3mag(); | |
77 | d1_prime= d1_prime/d1_prime.d3mag(); | |
78 | ||
79 | double x,y; | |
80 | ||
81 | x=d1_perp.dot(h1_perp); | |
82 | y=d1_prime.dot(h1_perp); | |
83 | ||
84 | double chi=atan2(y,x); | |
85 | ||
86 | if (chi<0.0) chi+=EvtConst::twoPi; | |
87 | ||
88 | return chi; | |
89 | ||
90 | } | |
91 | ||
92 | ||
93 | ||
94 | double EvtDecayPlaneNormalAngle(const EvtVector4R& p,const EvtVector4R& q, | |
95 | const EvtVector4R& d1,const EvtVector4R& d2){ | |
96 | ||
97 | EvtVector4C lc=dual(directProd(d1,d2)).cont2(q); | |
98 | ||
99 | EvtVector4R l(real(lc.get(0)),real(lc.get(1)), | |
100 | real(lc.get(2)),real(lc.get(3))); | |
101 | ||
102 | double pq=p*q; | |
103 | ||
104 | return q.mass()*(p*l)/sqrt(-(pq*pq-p.mass2()*q.mass2())*l.mass2()); | |
105 | ||
106 | ||
107 | } | |
108 | ||
109 | ||
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 ) | |
113 | { | |
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; | |
122 | ||
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); | |
125 | ||
126 | double phi = atan2(y,x); | |
127 | ||
128 | return phi<0 ? (phi+EvtConst::twoPi) : phi; | |
129 | } | |
130 | ||
131 | EvtComplex wignerD( int j, int m1, int m2, double phi, | |
132 | double theta, double gamma ) | |
133 | { | |
134 | ||
135 | EvtComplex gp(0.0, -phi*m1); | |
136 | EvtComplex gm(0.0, -gamma*m2); | |
137 | ||
138 | return exp( gp ) * EvtdFunction::d(j, m1, m2, theta) * exp( gm ); | |
139 | } | |
140 | ||
141 |