ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Angles/src/Angles.cc
(Generate patch)

Comparing UserCode/MitHzz4l/Angles/src/Angles.cc (file contents):
Revision 1.3 by khahn, Thu May 10 00:14:33 2012 UTC vs.
Revision 1.4 by dkralph, Tue Oct 23 10:10:07 2012 UTC

# Line 2 | Line 2
2   #include "Angles.h"
3   #include "EventData.h"
4  
5 + void fillAngles(Angles &a,
6 +                TLorentzVector lep1, int lep1_q,
7 +                TLorentzVector lep2, int lep2_q,
8 +                TLorentzVector lep3, int lep3_q,
9 +                TLorentzVector lep4, int lep4_q
10 +                ) {
11 +
12 +  LabVectors l;
13 +  
14 +  if( lep1_q > 0 ) {
15 +    l.vecl1p = lep1;  
16 +    l.vecl1m = lep2;  
17 +  } else {
18 +    l.vecl1p = lep2;  
19 +    l.vecl1m = lep1;  
20 +  }
21 +
22 +  if( lep3_q > 0 ) {
23 +    l.vecl2p = lep3;  
24 +    l.vecl2m = lep4;  
25 +  } else {
26 +    l.vecl2p = lep4;  
27 +    l.vecl2m = lep3;  
28 +  }
29 +
30 +  l.vecz1 = (l.vecl1p+l.vecl1m);
31 +  l.vecz2 = (l.vecl2p+l.vecl2m);
32 +  l.vec4l = (l.vecz1+l.vecz2);
33 +  
34 +  getAngles( l, a );
35 + }
36   void fillAngles( EventData &data, Angles &a ) {
37  
38    LabVectors l;
# Line 32 | Line 63 | void fillAngles( EventData &data, Angles
63   void getAngles( LabVectors &l,
64                  Angles &a ) {
65  
66 +  // cout << "ANGLES z1: " << setw(12) << l.vecz1.Pt() << setw(12) << l.vecz1.Eta() << setw(12) << l.vecz1.Phi() << endl;
67 +  // cout << "ANGLES z2: " << setw(12) << l.vecz2.Pt() << setw(12) << l.vecz2.Eta() << setw(12) << l.vecz2.Phi() << endl;
68 +  // cout << "ANGLES 4l: " << setw(12) << l.vec4l.Pt() << setw(12) << l.vec4l.Eta() << setw(12) << l.vec4l.Phi() << endl;
69 +  // cout << "ANGLES l1p:" << setw(12) << l.vecl1p.Pt() << setw(12) << l.vecl1p.Eta() << setw(12) << l.vecl1p.Phi() << endl;
70 +  // cout << "ANGLES l1m:" << setw(12) << l.vecl1m.Pt() << setw(12) << l.vecl1m.Eta() << setw(12) << l.vecl1m.Phi() << endl;
71 +  // cout << "ANGLES l2p:" << setw(12) << l.vecl2p.Pt() << setw(12) << l.vecl2p.Eta() << setw(12) << l.vecl2p.Phi() << endl;
72 +  // cout << "ANGLES l2m:" << setw(12) << l.vecl2m.Pt() << setw(12) << l.vecl2m.Eta() << setw(12) << l.vecl2m.Phi() << endl;
73    //
74    // define the frames
75    // --------------------------------------------------------------------------------------------------
# Line 184 | Line 222 | void getAngles( LabVectors &l,
222   };
223  
224  
225 + void computeMelaAngles(TLorentzVector p4M11, int Z1_lept1Id,
226 +                       TLorentzVector p4M12, int Z1_lept2Id,
227 +                       TLorentzVector p4M21, int Z2_lept1Id,
228 +                       TLorentzVector p4M22, int Z2_lept2Id,
229 +                       float& costhetastar,
230 +                       float& costheta1,
231 +                       float& costheta2,
232 +                       float& Phi,
233 +                       float& Phi1){
234 +  
235 +  //build Z 4-vectors
236 +  TLorentzVector p4Z1 = p4M11 + p4M12;
237 +  TLorentzVector p4Z2 = p4M21 + p4M22;
238 +
239 +  // Sort Z1 leptons so that:
240 +  if ( (Z1_lept1Id*Z1_lept2Id<0 && Z1_lept1Id<0) || // for OS pairs: lep1 must be the negative one
241 +       (Z1_lept1Id*Z1_lept2Id>0 && p4M11.Phi()<=p4M12.Phi()) //for SS pairs: use random deterministic convention
242 +       ) {
243 +    swap(p4M11, p4M12);
244 +  }
245 +  
246 +  // Same for Z2 leptons
247 +  if ( (Z2_lept1Id*Z2_lept2Id<0 && Z2_lept1Id<0) ||
248 +       (Z2_lept1Id*Z2_lept2Id>0 && p4M21.Phi()<=p4M22.Phi())
249 +       ) {
250 +    swap(p4M21, p4M22);
251 +  }
252 +
253 +  
254 +  //build H 4-vectors
255 +  TLorentzVector p4H = p4Z1 + p4Z2;
256 +
257 +  // -----------------------------------
258 +
259 +  //// costhetastar
260 +  TVector3 boostX = -(p4H.BoostVector());
261 +  TLorentzVector thep4Z1inXFrame( p4Z1 );
262 +  TLorentzVector thep4Z2inXFrame( p4Z2 );
263 +  thep4Z1inXFrame.Boost( boostX );
264 +  thep4Z2inXFrame.Boost( boostX );
265 +  TVector3 theZ1X_p3 = TVector3( thep4Z1inXFrame.X(), thep4Z1inXFrame.Y(), thep4Z1inXFrame.Z() );
266 +  TVector3 theZ2X_p3 = TVector3( thep4Z2inXFrame.X(), thep4Z2inXFrame.Y(), thep4Z2inXFrame.Z() );    
267 +  costhetastar = theZ1X_p3.CosTheta();
268 +
269 +  //// --------------------------- costheta1
270 +  TVector3 boostV1 = -(p4Z1.BoostVector());
271 +  TLorentzVector p4M11_BV1( p4M11 );
272 +  TLorentzVector p4M12_BV1( p4M12 );
273 +  TLorentzVector p4M21_BV1( p4M21 );
274 +  TLorentzVector p4M22_BV1( p4M22 );
275 +  p4M11_BV1.Boost( boostV1 );
276 +  p4M12_BV1.Boost( boostV1 );
277 +  p4M21_BV1.Boost( boostV1 );
278 +  p4M22_BV1.Boost( boostV1 );
279 +    
280 +  TLorentzVector p4V2_BV1 = p4M21_BV1 + p4M22_BV1;
281 +  //// costheta1
282 +  costheta1 = -p4V2_BV1.Vect().Dot( p4M11_BV1.Vect() )/p4V2_BV1.Vect().Mag()/p4M11_BV1.Vect().Mag();
283 +
284 +  //// --------------------------- costheta2
285 +  TVector3 boostV2 = -(p4Z2.BoostVector());
286 +  TLorentzVector p4M11_BV2( p4M11 );
287 +  TLorentzVector p4M12_BV2( p4M12 );
288 +  TLorentzVector p4M21_BV2( p4M21 );
289 +  TLorentzVector p4M22_BV2( p4M22 );
290 +  p4M11_BV2.Boost( boostV2 );
291 +  p4M12_BV2.Boost( boostV2 );
292 +  p4M21_BV2.Boost( boostV2 );
293 +  p4M22_BV2.Boost( boostV2 );
294 +    
295 +  TLorentzVector p4V1_BV2 = p4M11_BV2 + p4M12_BV2;
296 +  //// costheta2
297 +  costheta2 = -p4V1_BV2.Vect().Dot( p4M21_BV2.Vect() )/p4V1_BV2.Vect().Mag()/p4M21_BV2.Vect().Mag();
298 +    
299 +  //// --------------------------- Phi and Phi1 (old phistar1 - azimuthal production angle)
300 +  //    TVector3 boostX = -(thep4H.BoostVector());
301 +  TLorentzVector p4M11_BX( p4M11 );
302 +  TLorentzVector p4M12_BX( p4M12 );
303 +  TLorentzVector p4M21_BX( p4M21 );
304 +  TLorentzVector p4M22_BX( p4M22 );
305 +    
306 +  p4M11_BX.Boost( boostX );
307 +  p4M12_BX.Boost( boostX );
308 +  p4M21_BX.Boost( boostX );
309 +  p4M22_BX.Boost( boostX );
310 +    
311 +  TVector3 tmp1 = p4M11_BX.Vect().Cross( p4M12_BX.Vect() );
312 +  TVector3 tmp2 = p4M21_BX.Vect().Cross( p4M22_BX.Vect() );    
313 +    
314 +  TVector3 normal1_BX( tmp1.X()/tmp1.Mag(), tmp1.Y()/tmp1.Mag(), tmp1.Z()/tmp1.Mag() );
315 +  TVector3 normal2_BX( tmp2.X()/tmp2.Mag(), tmp2.Y()/tmp2.Mag(), tmp2.Z()/tmp2.Mag() );
316 +
317 +  //// Phi
318 +  TLorentzVector p4Z1_BX = p4M11_BX + p4M12_BX;    
319 +  float tmpSgnPhi = p4Z1_BX.Vect().Dot( normal1_BX.Cross( normal2_BX) );
320 +  float sgnPhi = tmpSgnPhi/fabs(tmpSgnPhi);
321 +  Phi = sgnPhi * acos( -1.*normal1_BX.Dot( normal2_BX) );
322 +    
323 +    
324 +  //////////////
325 +    
326 +  TVector3 beamAxis(0,0,1);
327 +  TVector3 tmp3 = (p4M11_BX + p4M12_BX).Vect();
328 +    
329 +  TVector3 p3V1_BX( tmp3.X()/tmp3.Mag(), tmp3.Y()/tmp3.Mag(), tmp3.Z()/tmp3.Mag() );
330 +  TVector3 tmp4 = beamAxis.Cross( p3V1_BX );
331 +  TVector3 normalSC_BX( tmp4.X()/tmp4.Mag(), tmp4.Y()/tmp4.Mag(), tmp4.Z()/tmp4.Mag() );
332 +        
333 +  //// Phi1
334 +  float tmpSgnPhi1 = p4Z1_BX.Vect().Dot( normal1_BX.Cross( normalSC_BX) );
335 +  float sgnPhi1 = tmpSgnPhi1/fabs(tmpSgnPhi1);    
336 +  Phi1 = sgnPhi1 * acos( normal1_BX.Dot( normalSC_BX) );    
337 +    
338 + }
339 +
340 + //----------------------------------------------------------------------------------------
341 + void Angles::dump()
342 + {
343 +  cout << "angles: "
344 +     << setw(12) << costheta1 << setw(12) << costheta2 << setw(12) << Phi
345 +     << setw(12) << costhetastar << setw(12) << Phi1
346 +     << endl;
347 + }
348 +  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines