Line data Source code
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 0 : double pd=p*d;
37 0 : double pq=p*q;
38 0 : double qd=q*d;
39 0 : double mp2=p.mass2();
40 0 : double mq2=q.mass2();
41 0 : double md2=d.mass2();
42 :
43 0 : double cost=(pd*mq2-pq*qd)/sqrt((pq*pq-mq2*mp2)*(qd*qd-mq2*md2));
44 :
45 0 : 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 0 : EvtVector4R p4_d1p,p4_h1p,p4_h2p,p4_d2p;
54 :
55 :
56 : // boost all vectors parent restframe
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 0 : p4_d1p=boostTo(p4_d1,p4_p,true);
60 0 : p4_d2p=boostTo(p4_d2,p4_p,true);
61 0 : p4_h1p=boostTo(p4_h1,p4_p,true);
62 0 : p4_h2p=boostTo(p4_h2,p4_p,true);
63 :
64 :
65 0 : EvtVector4R d1_perp,d1_prime,h1_perp;
66 0 : EvtVector4R D;
67 :
68 0 : D=p4_d1p+p4_d2p;
69 :
70 0 : d1_perp=p4_d1p-(D.dot(p4_d1p)/D.dot(D))*D;
71 0 : h1_perp=p4_h1p-(D.dot(p4_h1p)/D.dot(D))*D;
72 :
73 : // orthogonal to both D and d1_perp
74 :
75 0 : d1_prime=D.cross(d1_perp);
76 :
77 0 : d1_perp= d1_perp/d1_perp.d3mag();
78 0 : d1_prime= d1_prime/d1_prime.d3mag();
79 :
80 : double x,y;
81 :
82 0 : x=d1_perp.dot(h1_perp);
83 0 : y=d1_prime.dot(h1_perp);
84 :
85 0 : double chi=atan2(y,x);
86 :
87 0 : if (chi<0.0) chi+=EvtConst::twoPi;
88 :
89 0 : return chi;
90 :
91 0 : }
92 :
93 :
94 :
95 : double EvtDecayPlaneNormalAngle(const EvtVector4R& p,const EvtVector4R& q,
96 : const EvtVector4R& d1,const EvtVector4R& d2){
97 :
98 0 : EvtVector4C lc=dual(EvtGenFunctions::directProd(d1,d2)).cont2(q);
99 :
100 0 : EvtVector4R l(real(lc.get(0)),real(lc.get(1)),
101 0 : real(lc.get(2)),real(lc.get(3)));
102 :
103 0 : double pq=p*q;
104 :
105 0 : return q.mass()*(p*l)/sqrt(-(pq*pq-p.mass2()*q.mass2())*l.mass2());
106 :
107 :
108 0 : }
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 0 : double eq = (p * q) / p.mass();
116 0 : double ed = (p * d) / p.mass();
117 0 : double mq = q.mass();
118 0 : double q2 = p.mag2r3(q);
119 0 : double qd = p.dotr3(q,d);
120 0 : double zq = p.dotr3(z,q);
121 0 : double zd = p.dotr3(z,d);
122 0 : double alpha = (eq - mq)/(q2 * mq) * qd - ed/mq;
123 :
124 0 : double y = p.scalartripler3(z,q,d) + alpha * p.scalartripler3(z,q,q);
125 0 : double x = (zq * (qd + alpha * q2) - q2 * (zd + alpha * zq)) / sqrt(q2);
126 :
127 0 : double phi = atan2(y,x);
128 :
129 0 : 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 0 : EvtComplex gp(0.0, -phi*m1);
137 0 : EvtComplex gm(0.0, -gamma*m2);
138 :
139 0 : return exp( gp ) * EvtdFunction::d(j, m1, m2, theta) * exp( gm );
140 0 : }
141 :
142 :
|
