Upadte from Taku
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronPair.cxx
CommitLineData
b2a297fa 1/*************************************************************************
2* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
15
16///////////////////////////////////////////////////////////////////////////
17// //
18// Dielectron Pair class. Internally it makes use of AliKFParticle. //
19// //
20///////////////////////////////////////////////////////////////////////////
21
22
ba15fdfb 23#include <TDatabasePDG.h>
24#include <AliVTrack.h>
25#include <AliVVertex.h>
26#include <AliPID.h>
236e1bda 27#include <AliExternalTrackParam.h>
ba15fdfb 28
b2a297fa 29#include "AliDielectronPair.h"
b2a297fa 30
31ClassImp(AliDielectronPair)
32
33AliDielectronPair::AliDielectronPair() :
b2a297fa 34 fType(-1),
a655b716 35 fLabel(-1),
e4339752 36 fPdgCode(0),
b2a297fa 37 fPair(),
572b0139 38 fD1(),
39 fD2(),
b2a297fa 40 fRefD1(),
08b801a6 41 fRefD2(),
42 fKFUsage(kTRUE)
b2a297fa 43{
44 //
45 // Default Constructor
46 //
47
48}
49
50//______________________________________________
51AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
52 AliVTrack * const particle2, Int_t pid2, Char_t type) :
b2a297fa 53 fType(type),
a655b716 54 fLabel(-1),
e4339752 55 fPdgCode(0),
b2a297fa 56 fPair(),
572b0139 57 fD1(),
58 fD2(),
b2a297fa 59 fRefD1(),
08b801a6 60 fRefD2(),
61 fKFUsage(kTRUE)
b2a297fa 62{
63 //
64 // Constructor with tracks
65 //
66 SetTracks(particle1, pid1, particle2, pid2);
67}
68
69//______________________________________________
1201a1a9 70AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
71 const AliKFParticle * const particle2,
72 AliVTrack * const refParticle1,
73 AliVTrack * const refParticle2, Char_t type) :
74 fType(type),
75 fLabel(-1),
e4339752 76 fPdgCode(0),
1201a1a9 77 fPair(),
78 fD1(),
79 fD2(),
80 fRefD1(),
08b801a6 81 fRefD2(),
82 fKFUsage(kTRUE)
1201a1a9 83{
84 //
85 // Constructor with tracks
86 //
87 SetTracks(particle1, particle2,refParticle1,refParticle2);
88}
89
90//______________________________________________
b2a297fa 91AliDielectronPair::~AliDielectronPair()
92{
93 //
94 // Default Destructor
95 //
96
97}
98
99//______________________________________________
100void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
101 AliVTrack * const particle2, Int_t pid2)
102{
103 //
572b0139 104 // Sort particles by pt, first particle larget Pt
105 // set AliKF daughters and pair
1201a1a9 106 // refParticle1 and 2 are the original tracks. In the case of track rotation
107 // they are needed in the framework
b2a297fa 108 //
109 fPair.Initialize();
572b0139 110 fD1.Initialize();
111 fD2.Initialize();
8df8e382 112
b2a297fa 113 AliKFParticle kf1(*particle1,pid1);
114 AliKFParticle kf2(*particle2,pid2);
572b0139 115
b2a297fa 116 fPair.AddDaughter(kf1);
117 fPair.AddDaughter(kf2);
8df8e382 118
a655b716 119 if (particle1->Pt()>particle2->Pt()){
120 fRefD1 = particle1;
121 fRefD2 = particle2;
572b0139 122 fD1+=kf1;
123 fD2+=kf2;
a655b716 124 } else {
125 fRefD1 = particle2;
126 fRefD2 = particle1;
572b0139 127 fD1+=kf2;
128 fD2+=kf1;
a655b716 129 }
b2a297fa 130}
131
8df8e382 132//______________________________________________
1201a1a9 133void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
134 const AliKFParticle * const particle2,
135 AliVTrack * const refParticle1,
136 AliVTrack * const refParticle2)
137{
138 //
139 // Sort particles by pt, first particle larget Pt
140 // set AliKF daughters and pair
141 // refParticle1 and 2 are the original tracks. In the case of track rotation
142 // they are needed in the framework
143 //
144 fPair.Initialize();
145 fD1.Initialize();
146 fD2.Initialize();
147
148 AliKFParticle kf1(*particle1);
149 AliKFParticle kf2(*particle2);
150
151 fPair.AddDaughter(kf1);
152 fPair.AddDaughter(kf2);
153
154 if (kf1.GetPt()>kf2.GetPt()){
155 fRefD1 = refParticle1;
156 fRefD2 = refParticle2;
157 fD1+=kf1;
158 fD2+=kf2;
159 } else {
160 fRefD1 = refParticle2;
161 fRefD2 = refParticle1;
162 fD1+=kf2;
163 fD2+=kf1;
164 }
165}
166
167//______________________________________________
61d106d3 168void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
169{
170 //
171 // Calculate theta and phi in helicity and Collins-Soper coordinate frame
172 //
173 const Double_t kBeamEnergy = 3500.;
1201a1a9 174 Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
175 Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
61d106d3 176 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
177 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
178
1201a1a9 179// AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
180// AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
61d106d3 181
1201a1a9 182// d1->PxPyPz(pxyz1);
183// d2->PxPyPz(pxyz2);
61d106d3 184
185 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
186 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
187
188 // first & second daughter 4-mom
189 TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
190 TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
191 TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
192 TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
193 // J/Psi 4-momentum vector
194 TLorentzVector motherMom=p1Mom+p2Mom;
195
196 // boost all the 4-mom vectors to the mother rest frame
197 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
198 p1Mom.Boost(beta);
199 p2Mom.Boost(beta);
200 projMom.Boost(beta);
201 targMom.Boost(beta);
202
203 // x,y,z axes
204 TVector3 zAxisHE = (motherMom.Vect()).Unit();
205 TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
206 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
207 TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
208 TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
209
210 // fill theta and phi
1201a1a9 211 if(fD1.GetQ()>0){
61d106d3 212 thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
213 thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
214 phiHE = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
215 phiCS = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
216 } else {
217 thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
218 thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
219 phiHE = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
220 phiCS = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
221 }
222}
223
224//______________________________________________
236e1bda 225Double_t AliDielectronPair::PsiPair(Double_t MagField) const
226{
227 //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
228 //to ID conversions. Adapted from AliTRDv0Info class
229 Double_t x, y, z;
230 x = fPair.GetX();
231 y = fPair.GetY();
232 z = fPair.GetZ();
233
234 Double_t m1[3] = {0,0,0};
235 Double_t m2[3] = {0,0,0};
236
237 m1[0] = fD1.GetPx();
238 m1[1] = fD1.GetPy();
239 m1[2] = fD1.GetPz();
240
241 m2[0] = fD2.GetPx();
242 m2[1] = fD2.GetPy();
243 m2[2] = fD2.GetPz();
244
245 Double_t deltat = 1.;
246 deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]*m2[0] + m2[1]*m2[1])+1.e-13))-
247 TMath::ATan(m1[2]/(TMath::Sqrt(m1[0]*m1[0] + m1[1]*m1[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
248
249 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
250
251 Double_t mom1Prop[3];
252 Double_t mom2Prop[3];
253
254 AliExternalTrackParam *d1 = static_cast<AliExternalTrackParam*>(fRefD1.GetObject());
255 AliExternalTrackParam *d2 = static_cast<AliExternalTrackParam*>(fRefD2.GetObject());
256
257 AliExternalTrackParam nt(*d1), pt(*d2);
236e1bda 258
259 Double_t fPsiPair = 4.;
260 if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
261 fPsiPair = -5.;
262 if(pt.PropagateTo(radiussum,MagField) == 0)
263 fPsiPair = -5.;
264 pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
265 nt.GetPxPyPz(mom2Prop);
266
267
268
269 Double_t pEle =
270 TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
271 Double_t pPos =
272 TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val
273
274 Double_t scalarproduct =
275 mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit
276
277 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
278
279 fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
280
281 return fPsiPair;
282
283}
284
285//______________________________________________
8df8e382 286Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
61d106d3 287 const Bool_t isHE, const Bool_t isTheta)
288{
289 // The function calculates theta and phi in the mother rest frame with
8df8e382 290 // respect to the helicity coordinate system and Collins-Soper coordinate system
291 // TO DO: generalize for different decays (only J/Psi->e+e- now)
292
293 // Laboratory frame 4-vectors:
294 // projectile beam & target beam 4-mom
61d106d3 295 // TODO: need to retrieve the beam energy from somewhere
296 const Double_t kBeamEnergy = 3500.;
297 Double_t px1=d1->Px();
298 Double_t py1=d1->Py();
299 Double_t pz1=d1->Pz();
300 Double_t px2=d2->Px();
301 Double_t py2=d2->Py();
302 Double_t pz2=d2->Pz();
303 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
304 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
305
306 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
307 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
8df8e382 308
309 // first & second daughter 4-mom
61d106d3 310 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
311 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
8df8e382 312 // J/Psi 4-momentum vector
313 TLorentzVector motherMom=p1Mom+p2Mom;
314
315 // boost all the 4-mom vectors to the mother rest frame
316 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
317 p1Mom.Boost(beta);
318 p2Mom.Boost(beta);
319 projMom.Boost(beta);
320 targMom.Boost(beta);
321
322 // x,y,z axes
323 TVector3 zAxis;
324 if(isHE) zAxis = (motherMom.Vect()).Unit();
325 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
326 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
327 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
328
329 // return either theta or phi
330 if(isTheta) {
331 if(d1->Charge()>0)
332 return zAxis.Dot((p1Mom.Vect()).Unit());
333 else
334 return zAxis.Dot((p2Mom.Vect()).Unit());
335
336 }
337 else {
338 if(d1->Charge()>0)
339 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
340 else
341 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
342 }
343}
344
345//______________________________________________
346Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
347 // The function calculates theta and phi in the mother rest frame with
348 // respect to the helicity coordinate system and Collins-Soper coordinate system
349 // TO DO: generalize for different decays (only J/Psi->e+e- now)
350
351 // Laboratory frame 4-vectors:
352 // projectile beam & target beam 4-mom
45b2b1b8 353 AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
354 AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
61d106d3 355
356 const Double_t kBeamEnergy = 3500.;
357 Double_t px1=d1->Px();
358 Double_t py1=d1->Py();
359 Double_t pz1=d1->Pz();
360 Double_t px2=d2->Px();
361 Double_t py2=d2->Py();
362 Double_t pz2=d2->Pz();
363 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
364 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
365
366 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
367 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
368
369 // first & second daughter 4-mom
370 // first & second daughter 4-mom
371 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
372 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
8df8e382 373 // J/Psi 4-momentum vector
374 TLorentzVector motherMom=p1Mom+p2Mom;
375
376 // boost all the 4-mom vectors to the mother rest frame
377 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
378 p1Mom.Boost(beta);
379 p2Mom.Boost(beta);
380 projMom.Boost(beta);
381 targMom.Boost(beta);
382
383 // x,y,z axes
384 TVector3 zAxis;
385 if(isHE) zAxis = (motherMom.Vect()).Unit();
386 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
387 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
388 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
389
390 // return either theta or phi
391 if(isTheta) {
392 if(fD1.GetQ()>0)
393 return zAxis.Dot((p1Mom.Vect()).Unit());
394 else
395 return zAxis.Dot((p2Mom.Vect()).Unit());
396 }
397 else {
398 if(fD1.GetQ()>0)
399 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
400 else
401 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
402 }
403}
2e02dba4 404//______________________________________________
405Double_t AliDielectronPair::GetCosPointingAngle(const AliVVertex *primVtx) const
406{
407 //
408 // Calculate the poiting angle of the pair to the primary vertex and take the cosine
409 //
410 if(!primVtx) return -1.;
411
412 Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
413 deltaPos[0] = fPair.GetX() - primVtx->GetX();
414 deltaPos[1] = fPair.GetY() - primVtx->GetY();
415 deltaPos[2] = fPair.GetZ() - primVtx->GetZ();
416
417 Double_t momV02 = fPair.GetPx()*fPair.GetPx() + fPair.GetPy()*fPair.GetPy() + fPair.GetPz()*fPair.GetPz();
418 Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
419
420 Double_t cosinePointingAngle = (deltaPos[0]*fPair.GetPx() + deltaPos[1]*fPair.GetPy() + deltaPos[2]*fPair.GetPz()) / TMath::Sqrt(momV02 * deltaPos2);
421
422 return TMath::Abs(cosinePointingAngle);
423
424}
ba15fdfb 425
5720c765 426// //______________________________________________
427// Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
428// {
429// //
430// // Calculate the decay length in XY taking into account the primary vertex position
431// //
432// if(!vtx) return 0;
433// return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt() ;
434// }
ba15fdfb 435
5720c765 436// //______________________________________________
437// Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
438// {
439// //
440// // Calculate the pseudo proper time
441// //
442// Double_t lxy=GetLXY(vtx);
443// Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
444// return psProperDecayLength;
445// }
805bd069 446
447
448//______________________________________________
449Double_t AliDielectronPair::PhivPair(Double_t MagField) const
450{
451 //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
452 //to ID conversions. Angle between ee plane and magnetic field is calculated.
453
454 //Define local buffer variables for leg properties
455 Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
456 Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
457
458 if(MagField>0){
459 if(fD1.GetQ()>0){
460 px1 = fD1.GetPx();
461 py1 = fD1.GetPy();
462 pz1 = fD1.GetPz();
463
464 px2 = fD2.GetPx();
465 py2 = fD2.GetPy();
466 pz2 = fD2.GetPz();
467 }else{
468 px1 = fD2.GetPx();
469 py1 = fD2.GetPy();
470 pz1 = fD2.GetPz();
471
472 px2 = fD1.GetPx();
473 py2 = fD1.GetPy();
474 pz2 = fD1.GetPz();
475 }
476 }else{
477 if(fD1.GetQ()>0){
478 px1 = fD2.GetPx();
479 py1 = fD2.GetPy();
480 pz1 = fD2.GetPz();
481
482 px2 = fD1.GetPx();
483 py2 = fD1.GetPy();
484 pz2 = fD1.GetPz();
485 }else{
486 px1 = fD1.GetPx();
487 py1 = fD1.GetPy();
488 pz1 = fD1.GetPz();
489
490 px2 = fD2.GetPx();
491 py2 = fD2.GetPy();
492 pz2 = fD2.GetPz();
493 }
494 }
495
496 Double_t px = px1+px2;
497 Double_t py = py1+py2;
498 Double_t pz = pz1+pz2;
499 Double_t dppair = TMath::Sqrt(px*px+py*py+pz*pz);
500
501 //unit vector of (pep+pem)
502 Double_t pl = dppair;
503 Double_t ux = px/pl;
504 Double_t uy = py/pl;
505 Double_t uz = pz/pl;
506 Double_t ax = uy/TMath::Sqrt(ux*ux+uy*uy);
507 Double_t ay = -ux/TMath::Sqrt(ux*ux+uy*uy);
508
509 //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
510 //definition.
511 //Double_t ptep = iep->Px()*ax + iep->Py()*ay;
512 //Double_t ptem = iem->Px()*ax + iem->Py()*ay;
513
514 Double_t pxep = px1;
515 Double_t pyep = py1;
516 Double_t pzep = pz1;
517 Double_t pxem = px2;
518 Double_t pyem = py2;
519 Double_t pzem = pz2;
520
521 //vector product of pep X pem
522 Double_t vpx = pyep*pzem - pzep*pyem;
523 Double_t vpy = pzep*pxem - pxep*pzem;
524 Double_t vpz = pxep*pyem - pyep*pxem;
525 Double_t vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
526 //Double_t thev = acos(vpz/vp);
527
528 //unit vector of pep X pem
529 Double_t vx = vpx/vp;
530 Double_t vy = vpy/vp;
531 Double_t vz = vpz/vp;
532
533 //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
534 Double_t wx = uy*vz - uz*vy;
535 Double_t wy = uz*vx - ux*vz;
536 //Double_t wz = ux*vy - uy*vx;
537 //Double_t wl = sqrt(wx*wx+wy*wy+wz*wz);
538 // by construction, (wx,wy,wz) must be a unit vector.
539 // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
540 // should be small if the pair is conversion
541 //
542 Double_t cosPhiV = wx*ax + wy*ay;
543 Double_t phiv = TMath::ACos(cosPhiV);
544
545 return phiv;
546
547}
548
549
550