]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGDQ/dielectron/AliDielectronPair.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
update from pr task : sjena
[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>
1750c463 28#include <AliESDEvent.h>
ba15fdfb 29
b2a297fa 30#include "AliDielectronPair.h"
b2a297fa 31
32ClassImp(AliDielectronPair)
33
1750c463 34Double_t AliDielectronPair::fBeamEnergy=-1.;
35
b2a297fa 36AliDielectronPair::AliDielectronPair() :
b2a297fa 37 fType(-1),
a655b716 38 fLabel(-1),
e4339752 39 fPdgCode(0),
b2a297fa 40 fPair(),
572b0139 41 fD1(),
42 fD2(),
b2a297fa 43 fRefD1(),
08b801a6 44 fRefD2(),
45 fKFUsage(kTRUE)
b2a297fa 46{
47 //
48 // Default Constructor
49 //
50
51}
52
53//______________________________________________
54AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
55 AliVTrack * const particle2, Int_t pid2, Char_t type) :
b2a297fa 56 fType(type),
a655b716 57 fLabel(-1),
e4339752 58 fPdgCode(0),
b2a297fa 59 fPair(),
572b0139 60 fD1(),
61 fD2(),
b2a297fa 62 fRefD1(),
08b801a6 63 fRefD2(),
64 fKFUsage(kTRUE)
b2a297fa 65{
66 //
67 // Constructor with tracks
68 //
69 SetTracks(particle1, pid1, particle2, pid2);
70}
71
1201a1a9 72//______________________________________________
73AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
74 const AliKFParticle * const particle2,
75 AliVTrack * const refParticle1,
76 AliVTrack * const refParticle2, Char_t type) :
77 fType(type),
78 fLabel(-1),
e4339752 79 fPdgCode(0),
1201a1a9 80 fPair(),
81 fD1(),
82 fD2(),
83 fRefD1(),
08b801a6 84 fRefD2(),
85 fKFUsage(kTRUE)
1201a1a9 86{
87 //
88 // Constructor with tracks
89 //
90 SetTracks(particle1, particle2,refParticle1,refParticle2);
91}
92
b2a297fa 93//______________________________________________
94AliDielectronPair::~AliDielectronPair()
95{
96 //
97 // Default Destructor
98 //
99
100}
101
102//______________________________________________
103void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
104 AliVTrack * const particle2, Int_t pid2)
105{
106 //
572b0139 107 // Sort particles by pt, first particle larget Pt
108 // set AliKF daughters and pair
1201a1a9 109 // refParticle1 and 2 are the original tracks. In the case of track rotation
110 // they are needed in the framework
b2a297fa 111 //
112 fPair.Initialize();
572b0139 113 fD1.Initialize();
114 fD2.Initialize();
8df8e382 115
b2a297fa 116 AliKFParticle kf1(*particle1,pid1);
117 AliKFParticle kf2(*particle2,pid2);
572b0139 118
b2a297fa 119 fPair.AddDaughter(kf1);
120 fPair.AddDaughter(kf2);
8df8e382 121
4fae8ef9 122 if (particle1->Pt()>particle2->Pt()){
123 fRefD1 = particle1;
124 fRefD2 = particle2;
125 fD1+=kf1;
126 fD2+=kf2;
127 } else {
128 fRefD1 = particle2;
129 fRefD2 = particle1;
130 fD1+=kf2;
131 fD2+=kf1;
132 }
133}
134//______________________________________________
135void AliDielectronPair::SetGammaTracks(AliVTrack * const particle1, Int_t pid1,
136 AliVTrack * const particle2, Int_t pid2)
137{
138 //
139 // Sort particles by pt, first particle larget Pt
140 // set AliKF daughters and a GAMMA pair
141 // refParticle1 and 2 are the original tracks. In the case of track rotation
142 // they are needed in the framework
143 //
144 fD1.Initialize();
145 fD2.Initialize();
146
147 AliKFParticle kf1(*particle1,pid1);
148 AliKFParticle kf2(*particle2,pid2);
149 fPair.ConstructGamma(kf1,kf2);
150
a655b716 151 if (particle1->Pt()>particle2->Pt()){
152 fRefD1 = particle1;
153 fRefD2 = particle2;
572b0139 154 fD1+=kf1;
155 fD2+=kf2;
a655b716 156 } else {
157 fRefD1 = particle2;
158 fRefD2 = particle1;
572b0139 159 fD1+=kf2;
160 fD2+=kf1;
a655b716 161 }
b2a297fa 162}
163
1201a1a9 164//______________________________________________
165void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
166 const AliKFParticle * const particle2,
167 AliVTrack * const refParticle1,
168 AliVTrack * const refParticle2)
169{
170 //
171 // Sort particles by pt, first particle larget Pt
172 // set AliKF daughters and pair
173 // refParticle1 and 2 are the original tracks. In the case of track rotation
174 // they are needed in the framework
175 //
176 fPair.Initialize();
177 fD1.Initialize();
178 fD2.Initialize();
179
180 AliKFParticle kf1(*particle1);
181 AliKFParticle kf2(*particle2);
182
183 fPair.AddDaughter(kf1);
184 fPair.AddDaughter(kf2);
185
186 if (kf1.GetPt()>kf2.GetPt()){
187 fRefD1 = refParticle1;
188 fRefD2 = refParticle2;
189 fD1+=kf1;
190 fD2+=kf2;
191 } else {
192 fRefD1 = refParticle2;
193 fRefD2 = refParticle1;
194 fD1+=kf2;
195 fD2+=kf1;
196 }
197}
198
61d106d3 199//______________________________________________
200void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
201{
202 //
203 // Calculate theta and phi in helicity and Collins-Soper coordinate frame
204 //
1201a1a9 205 Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
206 Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
61d106d3 207 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
208 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
209
1201a1a9 210// AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
211// AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
61d106d3 212
1201a1a9 213// d1->PxPyPz(pxyz1);
214// d2->PxPyPz(pxyz2);
61d106d3 215
1750c463 216 TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
217 TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
61d106d3 218
219 // first & second daughter 4-mom
220 TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
221 TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
222 TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
223 TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
224 // J/Psi 4-momentum vector
225 TLorentzVector motherMom=p1Mom+p2Mom;
226
227 // boost all the 4-mom vectors to the mother rest frame
228 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
229 p1Mom.Boost(beta);
230 p2Mom.Boost(beta);
231 projMom.Boost(beta);
232 targMom.Boost(beta);
233
234 // x,y,z axes
235 TVector3 zAxisHE = (motherMom.Vect()).Unit();
236 TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
237 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
238 TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
239 TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
240
241 // fill theta and phi
1201a1a9 242 if(fD1.GetQ()>0){
61d106d3 243 thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
244 thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
245 phiHE = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
246 phiCS = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
247 } else {
248 thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
249 thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
250 phiHE = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
251 phiCS = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
252 }
253}
254
151ddcbc 255
236e1bda 256//______________________________________________
257Double_t AliDielectronPair::PsiPair(Double_t MagField) const
258{
259 //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
260 //to ID conversions. Adapted from AliTRDv0Info class
99345a64 261 Double_t x, y;//, z;
236e1bda 262 x = fPair.GetX();
263 y = fPair.GetY();
99345a64 264 // z = fPair.GetZ();
236e1bda 265
266 Double_t m1[3] = {0,0,0};
267 Double_t m2[3] = {0,0,0};
268
269 m1[0] = fD1.GetPx();
270 m1[1] = fD1.GetPy();
271 m1[2] = fD1.GetPz();
272
273 m2[0] = fD2.GetPx();
274 m2[1] = fD2.GetPy();
275 m2[2] = fD2.GetPz();
276
277 Double_t deltat = 1.;
278 deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]*m2[0] + m2[1]*m2[1])+1.e-13))-
279 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
280
281 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
282
283 Double_t mom1Prop[3];
284 Double_t mom2Prop[3];
285
286 AliExternalTrackParam *d1 = static_cast<AliExternalTrackParam*>(fRefD1.GetObject());
287 AliExternalTrackParam *d2 = static_cast<AliExternalTrackParam*>(fRefD2.GetObject());
288
289 AliExternalTrackParam nt(*d1), pt(*d2);
236e1bda 290
291 Double_t fPsiPair = 4.;
292 if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
293 fPsiPair = -5.;
294 if(pt.PropagateTo(radiussum,MagField) == 0)
295 fPsiPair = -5.;
296 pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
297 nt.GetPxPyPz(mom2Prop);
298
299
300
301 Double_t pEle =
302 TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
303 Double_t pPos =
304 TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val
305
306 Double_t scalarproduct =
307 mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit
308
309 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
310
311 fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
312
313 return fPsiPair;
314
315}
316
8df8e382 317//______________________________________________
318Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
fce8a399 319 Bool_t isHE, Bool_t isTheta)
61d106d3 320{
321 // The function calculates theta and phi in the mother rest frame with
8df8e382 322 // respect to the helicity coordinate system and Collins-Soper coordinate system
323 // TO DO: generalize for different decays (only J/Psi->e+e- now)
324
325 // Laboratory frame 4-vectors:
326 // projectile beam & target beam 4-mom
61d106d3 327 Double_t px1=d1->Px();
328 Double_t py1=d1->Py();
329 Double_t pz1=d1->Pz();
330 Double_t px2=d2->Px();
331 Double_t py2=d2->Py();
332 Double_t pz2=d2->Pz();
333 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
334 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
b2ad74d0 335// printf(" beam energy %f \n ", fBeamEnergy);
1750c463 336 TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
337 TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
8df8e382 338
339 // first & second daughter 4-mom
61d106d3 340 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
341 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
8df8e382 342 // J/Psi 4-momentum vector
343 TLorentzVector motherMom=p1Mom+p2Mom;
344
345 // boost all the 4-mom vectors to the mother rest frame
346 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
347 p1Mom.Boost(beta);
348 p2Mom.Boost(beta);
349 projMom.Boost(beta);
350 targMom.Boost(beta);
351
352 // x,y,z axes
353 TVector3 zAxis;
354 if(isHE) zAxis = (motherMom.Vect()).Unit();
355 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
356 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
357 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
358
359 // return either theta or phi
360 if(isTheta) {
361 if(d1->Charge()>0)
362 return zAxis.Dot((p1Mom.Vect()).Unit());
363 else
364 return zAxis.Dot((p2Mom.Vect()).Unit());
365
366 }
367 else {
368 if(d1->Charge()>0)
369 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
370 else
371 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
372 }
373}
374
375//______________________________________________
fce8a399 376Double_t AliDielectronPair::ThetaPhiCM(Bool_t isHE, Bool_t isTheta) const {
8df8e382 377 // The function calculates theta and phi in the mother rest frame with
378 // respect to the helicity coordinate system and Collins-Soper coordinate system
379 // TO DO: generalize for different decays (only J/Psi->e+e- now)
380
381 // Laboratory frame 4-vectors:
382 // projectile beam & target beam 4-mom
45b2b1b8 383 AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
384 AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
61d106d3 385
61d106d3 386 Double_t px1=d1->Px();
387 Double_t py1=d1->Py();
388 Double_t pz1=d1->Pz();
389 Double_t px2=d2->Px();
390 Double_t py2=d2->Py();
391 Double_t pz2=d2->Pz();
392 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
393 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
394
1750c463 395 TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
396 TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
61d106d3 397
398 // first & second daughter 4-mom
399 // first & second daughter 4-mom
400 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
401 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
8df8e382 402 // J/Psi 4-momentum vector
403 TLorentzVector motherMom=p1Mom+p2Mom;
404
405 // boost all the 4-mom vectors to the mother rest frame
406 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
407 p1Mom.Boost(beta);
408 p2Mom.Boost(beta);
409 projMom.Boost(beta);
410 targMom.Boost(beta);
411
412 // x,y,z axes
413 TVector3 zAxis;
414 if(isHE) zAxis = (motherMom.Vect()).Unit();
415 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
416 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
417 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
418
419 // return either theta or phi
420 if(isTheta) {
421 if(fD1.GetQ()>0)
422 return zAxis.Dot((p1Mom.Vect()).Unit());
423 else
424 return zAxis.Dot((p2Mom.Vect()).Unit());
425 }
426 else {
427 if(fD1.GetQ()>0)
428 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
429 else
430 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
431 }
432}
2e02dba4 433//______________________________________________
434Double_t AliDielectronPair::GetCosPointingAngle(const AliVVertex *primVtx) const
435{
436 //
437 // Calculate the poiting angle of the pair to the primary vertex and take the cosine
438 //
439 if(!primVtx) return -1.;
440
441 Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
442 deltaPos[0] = fPair.GetX() - primVtx->GetX();
443 deltaPos[1] = fPair.GetY() - primVtx->GetY();
444 deltaPos[2] = fPair.GetZ() - primVtx->GetZ();
445
1598d580 446 Double_t momV02 = Px()*Px() + Py()*Py() + Pz()*Pz();
2e02dba4 447 Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
448
1598d580 449 Double_t cosinePointingAngle = (deltaPos[0]*Px() + deltaPos[1]*Py() + deltaPos[2]*Pz()) / TMath::Sqrt(momV02 * deltaPos2);
2e02dba4 450
451 return TMath::Abs(cosinePointingAngle);
452
453}
ba15fdfb 454
99345a64 455//______________________________________________
456Double_t AliDielectronPair::GetArmAlpha() const
457{
458 //
459 // Calculate the Armenteros-Podolanski Alpha
460 //
461 Int_t qD1 = fD1.GetQ();
462
463 TVector3 momNeg( (qD1<0?fD1.GetPx():fD2.GetPx()),
464 (qD1<0?fD1.GetPy():fD2.GetPy()),
465 (qD1<0?fD1.GetPz():fD2.GetPz()) );
466 TVector3 momPos( (qD1<0?fD2.GetPx():fD1.GetPx()),
467 (qD1<0?fD2.GetPy():fD1.GetPy()),
468 (qD1<0?fD2.GetPz():fD1.GetPz()) );
469 TVector3 momTot(Px(),Py(),Pz());
470
471 Double_t lQlNeg = momNeg.Dot(momTot)/momTot.Mag();
472 Double_t lQlPos = momPos.Dot(momTot)/momTot.Mag();
473
474 return ((lQlPos - lQlNeg)/(lQlPos + lQlNeg));
475}
476
477//______________________________________________
478Double_t AliDielectronPair::GetArmPt() const
479{
480 //
481 // Calculate the Armenteros-Podolanski Pt
482 //
483 Int_t qD1 = fD1.GetQ();
484
485 TVector3 momNeg( (qD1<0?fD1.GetPx():fD2.GetPx()),
486 (qD1<0?fD1.GetPy():fD2.GetPy()),
487 (qD1<0?fD1.GetPz():fD2.GetPz()) );
488 TVector3 momTot(Px(),Py(),Pz());
489
490 return (momNeg.Perp(momTot));
491}
492
1598d580
JB		 493//______________________________________________
494void AliDielectronPair::GetDCA(const AliVVertex *primVtx, Double_t d0z0[2]) const
495{
496 //
497 // Calculate the dca of the mother with respect to the primary vertex
498 //
499 if(!primVtx) return;
500
501 d0z0[0] = TMath::Sqrt(TMath::Power(Xv()-primVtx->GetX(),2) +
502 TMath::Power(Yv()-primVtx->GetY(),2) );
503
504 d0z0[1] = Zv() - primVtx->GetZ();
505 return;
506
507}
508
5720c765 509// //______________________________________________
510// Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
511// {
512// //
513// // Calculate the decay length in XY taking into account the primary vertex position
514// //
515// if(!vtx) return 0;
516// return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt() ;
517// }
ba15fdfb 518
5720c765 519// //______________________________________________
520// Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
521// {
522// //
523// // Calculate the pseudo proper time
524// //
525// Double_t lxy=GetLXY(vtx);
526// Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
527// return psProperDecayLength;
528// }
805bd069 529
530
531//______________________________________________
532Double_t AliDielectronPair::PhivPair(Double_t MagField) const
533{
534 //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
535 //to ID conversions. Angle between ee plane and magnetic field is calculated.
536
537 //Define local buffer variables for leg properties
538 Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
539 Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
540
541 if(MagField>0){
542 if(fD1.GetQ()>0){
543 px1 = fD1.GetPx();
544 py1 = fD1.GetPy();
545 pz1 = fD1.GetPz();
546
547 px2 = fD2.GetPx();
548 py2 = fD2.GetPy();
549 pz2 = fD2.GetPz();
550 }else{
551 px1 = fD2.GetPx();
552 py1 = fD2.GetPy();
553 pz1 = fD2.GetPz();
554
555 px2 = fD1.GetPx();
556 py2 = fD1.GetPy();
557 pz2 = fD1.GetPz();
558 }
559 }else{
560 if(fD1.GetQ()>0){
561 px1 = fD2.GetPx();
562 py1 = fD2.GetPy();
563 pz1 = fD2.GetPz();
564
565 px2 = fD1.GetPx();
566 py2 = fD1.GetPy();
567 pz2 = fD1.GetPz();
568 }else{
569 px1 = fD1.GetPx();
570 py1 = fD1.GetPy();
571 pz1 = fD1.GetPz();
572
573 px2 = fD2.GetPx();
574 py2 = fD2.GetPy();
575 pz2 = fD2.GetPz();
576 }
577 }
578
579 Double_t px = px1+px2;
580 Double_t py = py1+py2;
581 Double_t pz = pz1+pz2;
582 Double_t dppair = TMath::Sqrt(px*px+py*py+pz*pz);
583
584 //unit vector of (pep+pem)
585 Double_t pl = dppair;
586 Double_t ux = px/pl;
587 Double_t uy = py/pl;
588 Double_t uz = pz/pl;
589 Double_t ax = uy/TMath::Sqrt(ux*ux+uy*uy);
590 Double_t ay = -ux/TMath::Sqrt(ux*ux+uy*uy);
591
592 //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
593 //definition.
594 //Double_t ptep = iep->Px()*ax + iep->Py()*ay;
595 //Double_t ptem = iem->Px()*ax + iem->Py()*ay;
596
597 Double_t pxep = px1;
598 Double_t pyep = py1;
599 Double_t pzep = pz1;
600 Double_t pxem = px2;
601 Double_t pyem = py2;
602 Double_t pzem = pz2;
603
604 //vector product of pep X pem
605 Double_t vpx = pyep*pzem - pzep*pyem;
606 Double_t vpy = pzep*pxem - pxep*pzem;
607 Double_t vpz = pxep*pyem - pyep*pxem;
608 Double_t vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
609 //Double_t thev = acos(vpz/vp);
610
611 //unit vector of pep X pem
612 Double_t vx = vpx/vp;
613 Double_t vy = vpy/vp;
614 Double_t vz = vpz/vp;
615
616 //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
617 Double_t wx = uy*vz - uz*vy;
618 Double_t wy = uz*vx - ux*vz;
619 //Double_t wz = ux*vy - uy*vx;
620 //Double_t wl = sqrt(wx*wx+wy*wy+wz*wz);
621 // by construction, (wx,wy,wz) must be a unit vector.
622 // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
623 // should be small if the pair is conversion
624 //
625 Double_t cosPhiV = wx*ax + wy*ay;
626 Double_t phiv = TMath::ACos(cosPhiV);
627
628 return phiv;
629
630}
631
151ddcbc 632//______________________________________________
fce8a399 633Double_t AliDielectronPair::GetPairPlaneAngle(Double_t v0rpH2, Int_t VariNum)const
151ddcbc 634{
635
d6c07bf1 636 // Calculate the angle between electron pair plane and variables
637 // kv0rpH2 is reaction plane angle using V0-A,C,AC,Random
151ddcbc 638
639 Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
640 Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
641
642 px1 = fD1.GetPx();
643 py1 = fD1.GetPy();
644 pz1 = fD1.GetPz();
645
646 px2 = fD2.GetPx();
647 py2 = fD2.GetPy();
648 pz2 = fD2.GetPz();
649
650 //p1+p2
651 Double_t px = px1+px2;
652 Double_t py = py1+py2;
d6c07bf1 653 Double_t pz = pz1+pz2;
151ddcbc 654
655 // normal vector of ee plane
656 Double_t pnorx = py1*pz2 - pz1*py2;
657 Double_t pnory = pz1*px2 - px1*pz2;
658 Double_t pnorz = px1*py2 - py1*px2;
659 Double_t pnor = TMath::Sqrt( pnorx*pnorx + pnory*pnory + pnorz*pnorz );
660
d6c07bf1 661 //unit vector
151ddcbc 662 Double_t upnx = -9999.;
663 Double_t upny = -9999.;
07bc5543 664 Double_t upnz = -9999.;
665
d6c07bf1 666 if (pnor !=0)
667 {
668 upnx= pnorx/pnor;
669 upny= pnory/pnor;
670 upnz= pnorz/pnor;
671 }
151ddcbc 672
151ddcbc 673
d6c07bf1 674 Double_t ax = -9999.;
675 Double_t ay = -9999.;
676 Double_t az = -9999.;
677
678 //variable 1
679 //seeing the angle between ee decay plane and reaction plane by using V0-A,C,AC,Random
680 if(VariNum == 1){
681 ax = TMath::Sin(v0rpH2);
682 ay = -TMath::Cos(v0rpH2);
683 az = 0.0;
684 }
685
686
687 //variable 2
688 //seeing the angle between ee decay plane and (p1+p2) rot ez
689 else if (VariNum == 2 ){
690 ax = py;
691 ay = -px;
692 az = 0.0;
693 }
694
695 //variable 3
696 //seeing the angle between ee decay plane and (p1+p2) rot (p1+p2)x'z
697 else if (VariNum == 3 ){
698 Double_t rotpx = px*TMath::Cos(v0rpH2)+py*TMath::Sin(v0rpH2);
699 //Double_t rotpy = 0.0;
700 // Double_t rotpz = pz;
701
702 ax = py*pz;
703 ay = pz*rotpx-pz*px;
704 az = -rotpx*py;
705 }
706
707 //variable 4
708 //seeing the angle between ee decay plane and (p1+p2) rot ey'
709 else if (VariNum == 4){
710 ax = 0.0;
711 ay = 0.0;
712 az = pz;
713 }
714
715 Double_t denomHelper = ax*ax + ay*ay +az*az;
716 Double_t uax = -9999.;
717 Double_t uay = -9999.;
718 Double_t uaz = -9999.;
719
720 if (denomHelper !=0) {
721 uax = ax/TMath::Sqrt(denomHelper);
722 uay = ay/TMath::Sqrt(denomHelper);
723 uaz = az/TMath::Sqrt(denomHelper);
724 }
725
726 //PM is the angle between Pair plane and a plane decided by using variable 1-4
727
728 Double_t cosPM = upnx*uax + upny*uay + upnz*uaz;
729 Double_t PM = TMath::ACos(cosPM);
730
731 //keep interval [0,pi/2]
732 if(PM > TMath::Pi()/2){
733 PM -= TMath::Pi();
734 PM *= -1.0;
735
736 }
737 return PM;
151ddcbc 738}
b2ad74d0 739
740
36a1e727
JW		 741//_______________________________________________
742Double_t AliDielectronPair::PairPlaneMagInnerProduct(Double_t ZDCrpH1) const
743{
744
745 // Calculate inner product of the strong magnetic field and electron pair plane
746
290227df 747 if(TMath::Abs(ZDCrpH1 - 999.) < 1e-10) return -9999.;
36a1e727
JW		 748
749 Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
750 Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
751
290227df 752 if(fD1.GetQ()<0){
753 px1 = fD1.GetPx();
754 py1 = fD1.GetPy();
755 pz1 = fD1.GetPz();
36a1e727 756
290227df 757 px2 = fD2.GetPx();
758 py2 = fD2.GetPy();
759 pz2 = fD2.GetPz();
760
761 }else{
762 px1 = fD2.GetPx();
763 py1 = fD2.GetPy();
764 pz1 = fD2.GetPz();
765
766 px2 = fD1.GetPx();
767 py2 = fD1.GetPy();
768 pz2 = fD1.GetPz();
769
770 }
36a1e727 771
36a1e727
JW		 772
773 // normal vector of ee plane
290227df 774 Double_t pnorx = py1*pz2 - pz1*py2;
775 Double_t pnory = pz1*px2 - px1*pz2;
776 Double_t pnorz = px1*py2 - py1*px2;
36a1e727
JW		 777 Double_t pnor = TMath::Sqrt( pnorx*pnorx + pnory*pnory + pnorz*pnorz );
778
779 //unit vector
780 Double_t upnx = -9999.;
781 Double_t upny = -9999.;
782 //Double_t upnz = -9999.;
783
290227df 784 if (TMath::Abs(pnor) < 1e-10) return -9999.;
785
36a1e727
JW		 786 upnx= pnorx/pnor;
787 upny= pnory/pnor;
788 //upnz= pnorz/pnor;
789
790 //direction of strong magnetic field
791 Double_t magx = TMath::Cos(ZDCrpH1+(TMath::Pi()/2));
792 Double_t magy = TMath::Sin(ZDCrpH1+(TMath::Pi()/2));
793
794 //inner product of strong magnetic field and ee plane
795 Double_t upnmag = upnx*magx + upny*magy;
796
797 return upnmag;
798}
799
d6c07bf1 800
1750c463 801//______________________________________________
802void AliDielectronPair::SetBeamEnergy(AliVEvent *ev, Double_t beamEbyHand)
803{
804 //
805 // set the beam energy (by hand in case of AODs)
806 //
807 if(ev->IsA()==AliESDEvent::Class())
808 fBeamEnergy = ((AliESDEvent*)ev)->GetBeamEnergy();
809 else
810 fBeamEnergy = beamEbyHand;
811}
151ddcbc 812
696b423a 813