| 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 |
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 |
|
| 23 |
//std::cout << "In calculate angles..." << std::endl;
|
| 24 |
|
| 25 |
float norm;
|
| 26 |
|
| 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 |
float x_m11 = unitM11.Dot(unitx_1); float y_m11 = unitM11.Dot(unity_1); float z_m11 = unitM11.Dot(unitz_1);
|
| 79 |
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 |
float x_m21 = unitM21.Dot(unitx_2); float y_m21 = unitM21.Dot(unity_2); float z_m21 = unitM21.Dot(unitz_2);
|
| 117 |
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 |
// a.mz1 = (v.vecl1m+v.vecl1p).M();
|
| 164 |
// a.mz2 = (v.vecl2m+v.vecl2p).M();
|
| 165 |
// a.m4l = (v.vec4l).M();
|
| 166 |
|
| 167 |
};
|
| 168 |
|
| 169 |
|
