ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Angles/src/Angles_JHU.cc
Revision: 1.2
Committed: Mon Feb 13 10:06:01 2012 UTC (13 years, 3 months ago) by khahn
Content type: text/plain
Branch: MAIN
CVS Tags: synced_FSR_2, synced_FSR, synched2, synched
Changes since 1.1: +14 -14 lines
Log Message: 
*** empty log message ***

File Contents

# User Rev Content
1 khahn 1.1 #include "Angles.h"
2     #include <iostream>
3    
4     void getAngles_JHU( LabVectors &v, Angles &a ) {
5    
6     TLorentzVector thep4H = v.vec4l;
7     TLorentzVector thep4Z1 = v.vecz1;
8     TLorentzVector thep4M11 = v.vecl1m;
9     TLorentzVector thep4M12 = v.vecl1p;
10     TLorentzVector thep4Z2 = v.vecz2;
11     TLorentzVector thep4M21 = v.vecl2m;
12     TLorentzVector thep4M22 = v.vecl2p;
13    
14 khahn 1.2 float& costheta1 = a.costheta1;
15     float& costheta2 = a.costheta2;
16     float& phi = a.Phi;
17     float& costhetastar = a.costhetastar;
18     float& phistar1 = a.Phi1;
19     float& phistar2 = a.Phi2;
20     float& phi1 = a.phi1;
21     float& phi2 = a.phi2;
22 khahn 1.1
23     //std::cout << "In calculate angles..." << std::endl;
24    
25 khahn 1.2 float norm;
26 khahn 1.1
27     TVector3 boostX = -(thep4H.BoostVector());
28     TLorentzVector thep4Z1inXFrame( thep4Z1 );
29     TLorentzVector thep4Z2inXFrame( thep4Z2 );
30     thep4Z1inXFrame.Boost( boostX );
31     thep4Z2inXFrame.Boost( boostX );
32     TVector3 theZ1X_p3 = TVector3( thep4Z1inXFrame.X(), thep4Z1inXFrame.Y(), thep4Z1inXFrame.Z() );
33     TVector3 theZ2X_p3 = TVector3( thep4Z2inXFrame.X(), thep4Z2inXFrame.Y(), thep4Z2inXFrame.Z() );
34    
35     // calculate phi1, phi2, costhetastar
36     phi1 = theZ1X_p3.Phi();
37     phi2 = theZ2X_p3.Phi();
38     ///////////////////////////////////////////////
39     // check for z1/z2 convention
40     ///////////////////////////////////////////////
41     TLorentzVector p4H, p4Z1, p4M11, p4M12, p4Z2, p4M21, p4M22;
42     p4Z1 = thep4Z1; p4M11 = thep4M11; p4M12 = thep4M12;
43     p4Z2 = thep4Z2; p4M21 = thep4M21; p4M22 = thep4M22;
44     costhetastar = theZ1X_p3.CosTheta();
45    
46     // now helicity angles................................
47     TVector3 boostZ1 = -(p4Z1.BoostVector());
48     TLorentzVector p4Z2Z1(p4Z2);
49     p4Z2Z1.Boost(boostZ1);
50     //find the decay axis
51     /////TVector3 unitx_1 = -Hep3Vector(p4Z2Z1);
52     TVector3 unitx_1( -p4Z2Z1.X(), -p4Z2Z1.Y(), -p4Z2Z1.Z() );
53     norm = 1/(unitx_1.Mag());
54     unitx_1*=norm;
55     //boost daughters of z2
56     TLorentzVector p4M21Z1(p4M21);
57     TLorentzVector p4M22Z1(p4M22);
58     p4M21Z1.Boost(boostZ1);
59     p4M22Z1.Boost(boostZ1);
60     //create z and y axes
61     /////TVector3 unitz_1 = Hep3Vector(p4M21Z1).cross(Hep3Vector(p4M22Z1));
62     TVector3 p4M21Z1_p3( p4M21Z1.X(), p4M21Z1.Y(), p4M21Z1.Z() );
63     TVector3 p4M22Z1_p3( p4M22Z1.X(), p4M22Z1.Y(), p4M22Z1.Z() );
64     TVector3 unitz_1 = p4M21Z1_p3.Cross( p4M22Z1_p3 );
65     norm = 1/(unitz_1.Mag());
66     unitz_1 *= norm;
67     TVector3 unity_1 = unitz_1.Cross(unitx_1);
68    
69     //caculate theta1
70     TLorentzVector p4M11Z1(p4M11);
71     p4M11Z1.Boost(boostZ1);
72     TVector3 p3M11( p4M11Z1.X(), p4M11Z1.Y(), p4M11Z1.Z() );
73     TVector3 unitM11 = p3M11.Unit();
74     std::cout << "orig.X(): " << unitM11.X() << "\t"
75     << "orig.Y(): " << unitM11.Y() << "\t"
76     << "orig.Z(): " << unitM11.Z() << "\t"
77     << std::endl;
78 khahn 1.2 float x_m11 = unitM11.Dot(unitx_1); float y_m11 = unitM11.Dot(unity_1); float z_m11 = unitM11.Dot(unitz_1);
79 khahn 1.1 std::cout << "dotx: " << x_m11 << "\t"
80     << "doty: " << y_m11 << "\t"
81     << "dotx: " << z_m11 << "\t"
82     << std::endl;
83     TVector3 M11_Z1frame(y_m11, z_m11, x_m11);
84     std::cout << "M11_Z1frame.X(): " << M11_Z1frame.X() << "\t"
85     << "M11_Z1frame.Y(): " << M11_Z1frame.Y() << "\t"
86     << "M11_Z1frame.Z(): " << M11_Z1frame.Z() << "\t"
87     << std::endl;
88     costheta1 = M11_Z1frame.CosTheta();
89     //std::cout << "theta1: " << M11_Z1frame.Theta() << std::endl;
90     //////-----------------------old way of calculating phi---------------/////////
91     phi = M11_Z1frame.Phi();
92    
93     //set axes for other system
94     TVector3 boostZ2 = -(p4Z2.BoostVector());
95     TLorentzVector p4Z1Z2(p4Z1);
96     p4Z1Z2.Boost(boostZ2);
97     TVector3 unitx_2( -p4Z1Z2.X(), -p4Z1Z2.Y(), -p4Z1Z2.Z() );
98     norm = 1/(unitx_2.Mag());
99     unitx_2*=norm;
100     //boost daughters of z2
101     TLorentzVector p4M11Z2(p4M11);
102     TLorentzVector p4M12Z2(p4M12);
103     p4M11Z2.Boost(boostZ2);
104     p4M12Z2.Boost(boostZ2);
105     TVector3 p4M11Z2_p3( p4M11Z2.X(), p4M11Z2.Y(), p4M11Z2.Z() );
106     TVector3 p4M12Z2_p3( p4M12Z2.X(), p4M12Z2.Y(), p4M12Z2.Z() );
107     TVector3 unitz_2 = p4M11Z2_p3.Cross( p4M12Z2_p3 );
108     norm = 1/(unitz_2.Mag());
109     unitz_2*=norm;
110     TVector3 unity_2 = unitz_2.Cross(unitx_2);
111     //calcuate theta2
112     TLorentzVector p4M21Z2(p4M21);
113     p4M21Z2.Boost(boostZ2);
114     TVector3 p3M21( p4M21Z2.X(), p4M21Z2.Y(), p4M21Z2.Z() );
115     TVector3 unitM21 = p3M21.Unit();
116 khahn 1.2 float x_m21 = unitM21.Dot(unitx_2); float y_m21 = unitM21.Dot(unity_2); float z_m21 = unitM21.Dot(unitz_2);
117 khahn 1.1 TVector3 M21_Z2frame(y_m21, z_m21, x_m21);
118     costheta2 = M21_Z2frame.CosTheta();
119    
120     // calculate phi
121     //calculating phi_n
122     TLorentzVector n_p4Z1inXFrame( p4Z1 );
123     TLorentzVector n_p4M11inXFrame( p4M11 );
124     n_p4Z1inXFrame.Boost( boostX );
125     n_p4M11inXFrame.Boost( boostX );
126     TVector3 n_p4Z1inXFrame_unit = n_p4Z1inXFrame.Vect().Unit();
127     TVector3 n_p4M11inXFrame_unit = n_p4M11inXFrame.Vect().Unit();
128     TVector3 n_unitz_1( n_p4Z1inXFrame_unit );
129     //// y-axis is defined by neg lepton cross z-axis
130     //// the subtle part is here...
131     //////////TVector3 n_unity_1 = n_p4M11inXFrame_unit.Cross( n_unitz_1 );
132     TVector3 n_unity_1 = n_unitz_1.Cross( n_p4M11inXFrame_unit );
133     TVector3 n_unitx_1 = n_unity_1.Cross( n_unitz_1 );
134    
135     TLorentzVector n_p4M21inXFrame( p4M21 );
136     n_p4M21inXFrame.Boost( boostX );
137     TVector3 n_p4M21inXFrame_unit = n_p4M21inXFrame.Vect().Unit();
138     //rotate into other plane
139     TVector3 n_p4M21inXFrame_unitprime( n_p4M21inXFrame_unit.Dot(n_unitx_1), n_p4M21inXFrame_unit.Dot(n_unity_1), n_p4M21inXFrame_unit.Dot(n_unitz_1) );
140    
141     /// and then calculate phistar1
142     TVector3 n_p4PartoninXFrame_unit( 0.0, 0.0, 1.0 );
143     TVector3 n_p4PartoninXFrame_unitprime( n_p4PartoninXFrame_unit.Dot(n_unitx_1), n_p4PartoninXFrame_unit.Dot(n_unity_1), n_p4PartoninXFrame_unit.Dot(n_unitz_1) );
144     // negative sign is for arrow convention in paper
145     phistar1 = (n_p4PartoninXFrame_unitprime.Phi());
146    
147     // and the calculate phistar2
148     TLorentzVector n_p4Z2inXFrame( p4Z2 );
149     n_p4Z2inXFrame.Boost( boostX );
150     TVector3 n_p4Z2inXFrame_unit = n_p4Z2inXFrame.Vect().Unit();
151     ///////TLorentzVector n_p4M21inXFrame( p4M21 );
152     //////n_p4M21inXFrame.Boost( boostX );
153     ////TVector3 n_p4M21inXFrame_unit = n_p4M21inXFrame.Vect().Unit();
154     TVector3 n_unitz_2( n_p4Z2inXFrame_unit );
155     //// y-axis is defined by neg lepton cross z-axis
156     //// the subtle part is here...
157     //////TVector3 n_unity_2 = n_p4M21inXFrame_unit.Cross( n_unitz_2 );
158     TVector3 n_unity_2 = n_unitz_2.Cross( n_p4M21inXFrame_unit );
159     TVector3 n_unitx_2 = n_unity_2.Cross( n_unitz_2 );
160     TVector3 n_p4PartoninZ2PlaneFrame_unitprime( n_p4PartoninXFrame_unit.Dot(n_unitx_2), n_p4PartoninXFrame_unit.Dot(n_unity_2), n_p4PartoninXFrame_unit.Dot(n_unitz_2) );
161     phistar2 = (n_p4PartoninZ2PlaneFrame_unitprime.Phi());
162    
163 khahn 1.2 // a.mz1 = (v.vecl1m+v.vecl1p).M();
164     // a.mz2 = (v.vecl2m+v.vecl2p).M();
165     // a.m4l = (v.vec4l).M();
166 khahn 1.1
167     };
168    
169