]>
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 | |
0ca57c2f | 57 | // This does not boost particle to parent rest frame !!! |
58 | // It goes from parents rest frame to frame where parent has given momentum. | |
59 | p4_d1p=boostTo(p4_d1,p4_p,true); | |
60 | p4_d2p=boostTo(p4_d2,p4_p,true); | |
61 | p4_h1p=boostTo(p4_h1,p4_p,true); | |
62 | p4_h2p=boostTo(p4_h2,p4_p,true); | |
da0e9ce3 | 63 | |
64 | ||
65 | EvtVector4R d1_perp,d1_prime,h1_perp; | |
66 | EvtVector4R D; | |
67 | ||
68 | D=p4_d1p+p4_d2p; | |
69 | ||
70 | d1_perp=p4_d1p-(D.dot(p4_d1p)/D.dot(D))*D; | |
71 | h1_perp=p4_h1p-(D.dot(p4_h1p)/D.dot(D))*D; | |
72 | ||
73 | // orthogonal to both D and d1_perp | |
74 | ||
75 | d1_prime=D.cross(d1_perp); | |
76 | ||
77 | d1_perp= d1_perp/d1_perp.d3mag(); | |
78 | d1_prime= d1_prime/d1_prime.d3mag(); | |
79 | ||
80 | double x,y; | |
81 | ||
82 | x=d1_perp.dot(h1_perp); | |
83 | y=d1_prime.dot(h1_perp); | |
84 | ||
85 | double chi=atan2(y,x); | |
86 | ||
87 | if (chi<0.0) chi+=EvtConst::twoPi; | |
88 | ||
89 | return chi; | |
90 | ||
91 | } | |
92 | ||
93 | ||
94 | ||
95 | double EvtDecayPlaneNormalAngle(const EvtVector4R& p,const EvtVector4R& q, | |
96 | const EvtVector4R& d1,const EvtVector4R& d2){ | |
97 | ||
0ca57c2f | 98 | EvtVector4C lc=dual(EvtGenFunctions::directProd(d1,d2)).cont2(q); |
da0e9ce3 | 99 | |
100 | EvtVector4R l(real(lc.get(0)),real(lc.get(1)), | |
101 | real(lc.get(2)),real(lc.get(3))); | |
102 | ||
103 | double pq=p*q; | |
104 | ||
105 | return q.mass()*(p*l)/sqrt(-(pq*pq-p.mass2()*q.mass2())*l.mass2()); | |
106 | ||
107 | ||
108 | } | |
109 | ||
110 | ||
111 | // Calculate phi using the given 4 vectors (all in the same frame) | |
112 | double EvtDecayAnglePhi( const EvtVector4R& z, const EvtVector4R& p, const | |
113 | EvtVector4R& q, const EvtVector4R& d ) | |
114 | { | |
115 | double eq = (p * q) / p.mass(); | |
116 | double ed = (p * d) / p.mass(); | |
117 | double mq = q.mass(); | |
118 | double q2 = p.mag2r3(q); | |
119 | double qd = p.dotr3(q,d); | |
120 | double zq = p.dotr3(z,q); | |
121 | double zd = p.dotr3(z,d); | |
122 | double alpha = (eq - mq)/(q2 * mq) * qd - ed/mq; | |
123 | ||
124 | double y = p.scalartripler3(z,q,d) + alpha * p.scalartripler3(z,q,q); | |
125 | double x = (zq * (qd + alpha * q2) - q2 * (zd + alpha * zq)) / sqrt(q2); | |
126 | ||
127 | double phi = atan2(y,x); | |
128 | ||
129 | return phi<0 ? (phi+EvtConst::twoPi) : phi; | |
130 | } | |
131 | ||
132 | EvtComplex wignerD( int j, int m1, int m2, double phi, | |
133 | double theta, double gamma ) | |
134 | { | |
135 | ||
136 | EvtComplex gp(0.0, -phi*m1); | |
137 | EvtComplex gm(0.0, -gamma*m2); | |
138 | ||
139 | return exp( gp ) * EvtdFunction::d(j, m1, m2, theta) * exp( gm ); | |
140 | } | |
141 | ||
142 |