1 |
– |
#include "Angles.h" |
2 |
– |
#include <iostream> |
3 |
– |
|
4 |
– |
void getAngles( LabVectors &l, |
5 |
– |
Angles &a ) { |
6 |
– |
|
7 |
– |
// |
8 |
– |
// define the frames |
9 |
– |
// -------------------------------------------------------------------------------------------------- |
10 |
– |
TVector3 Xframe = l.vec4l.BoostVector(); |
11 |
– |
TVector3 Z1frame = l.vecz1.BoostVector(); |
12 |
– |
TVector3 Z2frame = l.vecz2.BoostVector(); |
13 |
– |
|
14 |
– |
// |
15 |
– |
// costheta(l1m,X) in Z1 frame |
16 |
– |
// -------------------------------------------------------------------------------------------------- |
17 |
– |
TLorentzVector vecl1m_in_Z1frame( l.vecl1m ); // make a copies of the vectors |
18 |
– |
TLorentzVector vec4l_in_Z1frame( l.vec4l ); // |
19 |
– |
vecl1m_in_Z1frame.Boost( -1*Z1frame ); // now go to the Z1 rest frame |
20 |
– |
vec4l_in_Z1frame.Boost( -1*Z1frame ); // |
21 |
– |
a.costheta1 |
22 |
– |
= vecl1m_in_Z1frame.Vect().Dot(-1*vec4l_in_Z1frame.Vect()) / |
23 |
– |
(vecl1m_in_Z1frame.Vect().Mag()*vec4l_in_Z1frame.Vect().Mag()); |
24 |
– |
|
25 |
– |
// |
26 |
– |
// costheta(l3,X) in Z2 frame |
27 |
– |
// -------------------------------------------------------------------------------------------------- |
28 |
– |
TLorentzVector vecl2m_in_Z2frame( l.vecl2m ); // make copies of the vectors |
29 |
– |
TLorentzVector vec4l_in_Z2frame( l.vec4l ); // |
30 |
– |
vecl2m_in_Z2frame.Boost( -1*Z2frame ); // now go to the Z2 rest frame |
31 |
– |
vec4l_in_Z2frame.Boost( -1*Z2frame ); |
32 |
– |
a.costheta2 = |
33 |
– |
vecl2m_in_Z2frame.Vect().Dot(-1*vec4l_in_Z2frame.Vect()) / |
34 |
– |
(vecl2m_in_Z2frame.Vect().Mag()*vec4l_in_Z2frame.Vect().Mag()); |
35 |
– |
|
36 |
– |
|
37 |
– |
// |
38 |
– |
// Phi(Z1,Z2) in X frame |
39 |
– |
// -------------------------------------------------------------------------------------------------- |
40 |
– |
|
41 |
– |
// this is essentially what they do and it reproduces their numbers, |
42 |
– |
// it's the angle between l1 and the z2 decay plane in the Z1 frame |
43 |
– |
TLorentzVector vecz2_in_Z1frame( l.vecz2 ); // |
44 |
– |
TLorentzVector vecl2m_in_Z1frame( l.vecl2m ); // |
45 |
– |
TLorentzVector vecl2p_in_Z1frame( l.vecl2p ); // |
46 |
– |
vecz2_in_Z1frame.Boost( -1*Z1frame ); // |
47 |
– |
vecl2m_in_Z1frame.Boost( -1*Z1frame ); // |
48 |
– |
vecl2p_in_Z1frame.Boost( -1*Z1frame ); // |
49 |
– |
|
50 |
– |
// |
51 |
– |
// define the Z2 coordinate system in the Z1 frame |
52 |
– |
// x : opposite the Z2 dir, lies in the Z2 decay plane |
53 |
– |
// z : orthogonal to the Z2 decay plane |
54 |
– |
// y : orthogonal to x&z, lies in the Z2 decay plane |
55 |
– |
TVector3 nnewX = ( -1*vecz2_in_Z1frame.Vect()).Unit(); |
56 |
– |
TVector3 nnewZ = (vecl2m_in_Z1frame.Vect().Cross(vecl2p_in_Z1frame.Vect())).Unit(); |
57 |
– |
TVector3 nnewY = nnewZ.Cross(nnewX).Unit(); |
58 |
– |
TRotation rotation2; |
59 |
– |
rotation2 = (rotation2.RotateAxes( nnewX, nnewY, nnewZ )).Inverse(); |
60 |
– |
TVector3 vecl1m_in_Z2coords_Z1frame = vecl1m_in_Z1frame.Vect().Transform(rotation2); |
61 |
– |
|
62 |
– |
// |
63 |
– |
// axes are redfined so that Phi is the angle between the Z1 orthogonal |
64 |
– |
// and the Z2 decay plane as given by y |
65 |
– |
// |
66 |
– |
TRotation rotation3; |
67 |
– |
TVector3 nnnewX( 0,1,0 ); |
68 |
– |
TVector3 nnnewY( 0,0,1 ); |
69 |
– |
TVector3 nnnewZ( 1,0,0 ); |
70 |
– |
rotation3 = (rotation3.RotateAxes( nnnewX, nnnewY, nnnewZ )).Inverse(); |
71 |
– |
TVector3 vecl1m_in_Z2coords_redefined_Z1frame = vecl1m_in_Z2coords_Z1frame.Transform(rotation3); |
72 |
– |
a.Phi = vecl1m_in_Z2coords_redefined_Z1frame.Phi(); |
73 |
– |
|
74 |
– |
|
75 |
– |
// |
76 |
– |
// rather, I'd expect something like this ... |
77 |
– |
// |
78 |
– |
TLorentzVector vecl1p_in_Xframe( l.vecl1p ); // make a copies of the vectors |
79 |
– |
TLorentzVector vecl1m_in_Xframe( l.vecl1m ); // |
80 |
– |
TLorentzVector vecl2p_in_Xframe( l.vecl2p ); // |
81 |
– |
TLorentzVector vecl2m_in_Xframe( l.vecl2m ); // |
82 |
– |
|
83 |
– |
vecl1p_in_Xframe.Boost( -1*Xframe ); // go to the X rest frame |
84 |
– |
vecl1m_in_Xframe.Boost( -1*Xframe ); // |
85 |
– |
vecl2p_in_Xframe.Boost( -1*Xframe ); // |
86 |
– |
vecl2m_in_Xframe.Boost( -1*Xframe ); // |
87 |
– |
|
88 |
– |
/* |
89 |
– |
// sanity check ... |
90 |
– |
TLorentzVector z1inX = vecl1p_in_Xframe+vecl1m_in_Xframe; |
91 |
– |
TLorentzVector z2inX = vecl2p_in_Xframe+vecl2m_in_Xframe; |
92 |
– |
float dotz1z2 = z1inX.Vect().Dot(z2inX.Vect()) / (z1inX.Vect().Mag()*z2inX.Vect().Mag()); |
93 |
– |
std::cout << "dotz1z2: " << dotz1z2 << std::endl; |
94 |
– |
*/ |
95 |
– |
|
96 |
– |
TVector3 Z1ortho_in_Xframe = // define the decay planes via the orthogonals |
97 |
– |
vecl1m_in_Xframe.Vect().Cross(vecl1p_in_Xframe.Vect()); |
98 |
– |
TVector3 Z2ortho_in_Xframe = |
99 |
– |
vecl2m_in_Xframe.Vect().Cross(vecl2p_in_Xframe.Vect()); |
100 |
– |
|
101 |
– |
/* |
102 |
– |
// sanity check ... why are these small but not exactly zero ??? |
103 |
– |
float z1odot = Z1ortho_in_Xframe.Dot(z1inX.Vect())/(Z1ortho_in_Xframe.Mag()*z1inX.Vect().Mag()); |
104 |
– |
float z2odot = Z2ortho_in_Xframe.Dot(z2inX.Vect())/(Z2ortho_in_Xframe.Mag()*z2inX.Vect().Mag()); |
105 |
– |
float l1podot = Z1ortho_in_Xframe.Dot(vecl1p_in_Xframe.Vect())/(Z1ortho_in_Xframe.Mag()*vecl1p_in_Xframe.Vect().Mag()); |
106 |
– |
float l1modot = Z1ortho_in_Xframe.Dot(vecl1m_in_Xframe.Vect())/(Z1ortho_in_Xframe.Mag()*vecl1m_in_Xframe.Vect().Mag()); |
107 |
– |
std::cout << "Z1ortho_in_Xframe.Dot(z1inX): " << z1odot << "\t" |
108 |
– |
<< "Z2ortho_in_Xframe.Dot(z2inX): " << z2odot << "\t" |
109 |
– |
<< "Z1ortho_in_Xframe.Dot(veclp1_in_Xframe) : " << l1podot << "\t" |
110 |
– |
<< "Z1ortho_in_Xframe.Dot(veclm1_in_Xframe) : " << l1modot << "\t" |
111 |
– |
<< std::endl; |
112 |
– |
*/ |
113 |
– |
|
114 |
– |
|
115 |
– |
// a.Phi = Z1ortho_in_Xframe.DeltaPhi(Z2ortho_in_Xframe); |
116 |
– |
//a.Phi = Z1ortho_in_Xframe.Angle(Z2ortho_in_Xframe); |
117 |
– |
|
118 |
– |
|
119 |
– |
|
120 |
– |
// |
121 |
– |
// costheta*(Z1,pp) in X frame |
122 |
– |
// -------------------------------------------------------------------------------------------------- |
123 |
– |
TLorentzVector vecz1_in_Xframe( l.vecz1 ); // make a copies of the vectors |
124 |
– |
TLorentzVector vecpp_in_Xframe( 0., 0., 1., 0. ); // define pp axis |
125 |
– |
|
126 |
– |
vecz1_in_Xframe.Boost( -1*Xframe ); // go to the X rest frame |
127 |
– |
vecpp_in_Xframe.Boost( -1*Xframe ); // |
128 |
– |
|
129 |
– |
a.costhetastar = |
130 |
– |
vecz1_in_Xframe.Vect().Dot(vecpp_in_Xframe.Vect()) / |
131 |
– |
(vecz1_in_Xframe.Vect().Mag()*vecpp_in_Xframe.Vect().Mag()); |
132 |
– |
|
133 |
– |
|
134 |
– |
// |
135 |
– |
// Phi1(Z1,pp) in X frame |
136 |
– |
// -------------------------------------------------------------------------------------------------- |
137 |
– |
|
138 |
– |
// define the Z1 decay plane in the Xframe |
139 |
– |
// y: orthogonal to the Z1 decay plane |
140 |
– |
// x: lies in the Z1 decay plane, but orthogonal to Z1 |
141 |
– |
// z: z1 in Xframe |
142 |
– |
TVector3 newY = (vecl1p_in_Xframe.Vect().Unit().Cross(vecl1m_in_Xframe.Vect().Unit())).Unit(); |
143 |
– |
TVector3 newX = (newY.Cross(vecz1_in_Xframe.Vect().Unit())).Unit(); |
144 |
– |
TVector3 newZ( (vecz1_in_Xframe.Vect()).Unit() ); |
145 |
– |
|
146 |
– |
TRotation rotation; |
147 |
– |
rotation = (rotation.RotateAxes( newX, newY, newZ )).Inverse(); |
148 |
– |
TVector3 vecpp( 0., 0., 1. ); |
149 |
– |
TVector3 vecpp_in_Z1coords_Xframe = vecpp.Transform(rotation); |
150 |
– |
a.Phi1 = vecpp_in_Z1coords_Xframe.Phi(); |
151 |
– |
|
152 |
– |
a.mz1 = (l.vecl1m+l.vecl1p).M(); |
153 |
– |
a.mz2 = (l.vecl2m+l.vecl2p).M(); |
154 |
– |
a.m4l = (l.vec4l).M(); |
155 |
– |
|
156 |
– |
}; |
157 |
– |
|
158 |
– |
|