1 /*************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////
18 // Dielectron Pair class. Internally it makes use of AliKFParticle. //
20 ///////////////////////////////////////////////////////////////////////////
23 #include <TDatabasePDG.h>
24 #include <AliVTrack.h>
25 #include <AliVVertex.h>
27 #include <AliExternalTrackParam.h>
28 #include <AliESDEvent.h>
30 #include "AliDielectronPair.h"
32 ClassImp(AliDielectronPair)
34 Double_t AliDielectronPair::fBeamEnergy=-1.;
36 AliDielectronPair::AliDielectronPair() :
48 // Default Constructor
53 //______________________________________________
54 AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
55 AliVTrack * const particle2, Int_t pid2, Char_t type) :
67 // Constructor with tracks
69 SetTracks(particle1, pid1, particle2, pid2);
72 //______________________________________________
73 AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
74 const AliKFParticle * const particle2,
75 AliVTrack * const refParticle1,
76 AliVTrack * const refParticle2, Char_t type) :
88 // Constructor with tracks
90 SetTracks(particle1, particle2,refParticle1,refParticle2);
93 //______________________________________________
94 AliDielectronPair::~AliDielectronPair()
102 //______________________________________________
103 void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
104 AliVTrack * const particle2, Int_t pid2)
107 // Sort particles by pt, first particle larget Pt
108 // set AliKF daughters and pair
109 // refParticle1 and 2 are the original tracks. In the case of track rotation
110 // they are needed in the framework
116 AliKFParticle kf1(*particle1,pid1);
117 AliKFParticle kf2(*particle2,pid2);
119 fPair.AddDaughter(kf1);
120 fPair.AddDaughter(kf2);
122 if (particle1->Pt()>particle2->Pt()){
134 //______________________________________________
135 void AliDielectronPair::SetGammaTracks(AliVTrack * const particle1, Int_t pid1,
136 AliVTrack * const particle2, Int_t pid2)
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
147 AliKFParticle kf1(*particle1,pid1);
148 AliKFParticle kf2(*particle2,pid2);
149 fPair.ConstructGamma(kf1,kf2);
151 if (particle1->Pt()>particle2->Pt()){
164 //______________________________________________
165 void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
166 const AliKFParticle * const particle2,
167 AliVTrack * const refParticle1,
168 AliVTrack * const refParticle2)
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
180 AliKFParticle kf1(*particle1);
181 AliKFParticle kf2(*particle2);
183 fPair.AddDaughter(kf1);
184 fPair.AddDaughter(kf2);
186 if (kf1.GetPt()>kf2.GetPt()){
187 fRefD1 = refParticle1;
188 fRefD2 = refParticle2;
192 fRefD1 = refParticle2;
193 fRefD2 = refParticle1;
199 //______________________________________________
200 void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
203 // Calculate theta and phi in helicity and Collins-Soper coordinate frame
205 Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
206 Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
207 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
208 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
210 // AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
211 // AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
213 // d1->PxPyPz(pxyz1);
214 // d2->PxPyPz(pxyz2);
216 TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
217 TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
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;
227 // boost all the 4-mom vectors to the mother rest frame
228 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
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();
241 // fill theta and phi
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));
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));
255 //______________________________________________
256 void AliDielectronPair::GetRotPair(Double_t &RotPairx, Double_t &RotPairy, Double_t &RotPairz) const
258 // calculation of rotation p1 p2
259 Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
260 Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
270 // normal vector of ee plane
271 Double_t pnorx = py1*pz2 - pz1*py2;
272 Double_t pnory = pz1*px2 - px1*pz2;
273 Double_t pnorz = px1*py2 - py1*px2;
274 Double_t pnor = TMath::Sqrt( pnorx*pnorx + pnory*pnory + pnorz*pnorz );
277 Double_t upnx = -9999.;
278 Double_t upny = -9999.;
279 Double_t upnz = -9999.;
295 //______________________________________________
296 Double_t AliDielectronPair::PsiPair(Double_t MagField) const
298 //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
299 //to ID conversions. Adapted from AliTRDv0Info class
305 Double_t m1[3] = {0,0,0};
306 Double_t m2[3] = {0,0,0};
316 Double_t deltat = 1.;
317 deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]*m2[0] + m2[1]*m2[1])+1.e-13))-
318 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
320 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
322 Double_t mom1Prop[3];
323 Double_t mom2Prop[3];
325 AliExternalTrackParam *d1 = static_cast<AliExternalTrackParam*>(fRefD1.GetObject());
326 AliExternalTrackParam *d2 = static_cast<AliExternalTrackParam*>(fRefD2.GetObject());
328 AliExternalTrackParam nt(*d1), pt(*d2);
330 Double_t fPsiPair = 4.;
331 if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
333 if(pt.PropagateTo(radiussum,MagField) == 0)
335 pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
336 nt.GetPxPyPz(mom2Prop);
341 TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
343 TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val
345 Double_t scalarproduct =
346 mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit
348 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
350 fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
356 //______________________________________________
357 Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
358 const Bool_t isHE, const Bool_t isTheta)
360 // The function calculates theta and phi in the mother rest frame with
361 // respect to the helicity coordinate system and Collins-Soper coordinate system
362 // TO DO: generalize for different decays (only J/Psi->e+e- now)
364 // Laboratory frame 4-vectors:
365 // projectile beam & target beam 4-mom
366 Double_t px1=d1->Px();
367 Double_t py1=d1->Py();
368 Double_t pz1=d1->Pz();
369 Double_t px2=d2->Px();
370 Double_t py2=d2->Py();
371 Double_t pz2=d2->Pz();
372 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
373 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
374 // printf(" beam energy %f \n ", fBeamEnergy);
375 TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
376 TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
378 // first & second daughter 4-mom
379 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
380 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
381 // J/Psi 4-momentum vector
382 TLorentzVector motherMom=p1Mom+p2Mom;
384 // boost all the 4-mom vectors to the mother rest frame
385 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
393 if(isHE) zAxis = (motherMom.Vect()).Unit();
394 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
395 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
396 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
398 // return either theta or phi
401 return zAxis.Dot((p1Mom.Vect()).Unit());
403 return zAxis.Dot((p2Mom.Vect()).Unit());
408 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
410 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
414 //______________________________________________
415 Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
416 // The function calculates theta and phi in the mother rest frame with
417 // respect to the helicity coordinate system and Collins-Soper coordinate system
418 // TO DO: generalize for different decays (only J/Psi->e+e- now)
420 // Laboratory frame 4-vectors:
421 // projectile beam & target beam 4-mom
422 AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
423 AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
425 Double_t px1=d1->Px();
426 Double_t py1=d1->Py();
427 Double_t pz1=d1->Pz();
428 Double_t px2=d2->Px();
429 Double_t py2=d2->Py();
430 Double_t pz2=d2->Pz();
431 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
432 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
434 TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
435 TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
437 // first & second daughter 4-mom
438 // first & second daughter 4-mom
439 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
440 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
441 // J/Psi 4-momentum vector
442 TLorentzVector motherMom=p1Mom+p2Mom;
444 // boost all the 4-mom vectors to the mother rest frame
445 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
453 if(isHE) zAxis = (motherMom.Vect()).Unit();
454 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
455 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
456 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
458 // return either theta or phi
461 return zAxis.Dot((p1Mom.Vect()).Unit());
463 return zAxis.Dot((p2Mom.Vect()).Unit());
467 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
469 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
472 //______________________________________________
473 Double_t AliDielectronPair::GetCosPointingAngle(const AliVVertex *primVtx) const
476 // Calculate the poiting angle of the pair to the primary vertex and take the cosine
478 if(!primVtx) return -1.;
480 Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
481 deltaPos[0] = fPair.GetX() - primVtx->GetX();
482 deltaPos[1] = fPair.GetY() - primVtx->GetY();
483 deltaPos[2] = fPair.GetZ() - primVtx->GetZ();
485 Double_t momV02 = fPair.GetPx()*fPair.GetPx() + fPair.GetPy()*fPair.GetPy() + fPair.GetPz()*fPair.GetPz();
486 Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
488 Double_t cosinePointingAngle = (deltaPos[0]*fPair.GetPx() + deltaPos[1]*fPair.GetPy() + deltaPos[2]*fPair.GetPz()) / TMath::Sqrt(momV02 * deltaPos2);
490 return TMath::Abs(cosinePointingAngle);
494 //______________________________________________
495 Double_t AliDielectronPair::GetArmAlpha() const
498 // Calculate the Armenteros-Podolanski Alpha
500 Int_t qD1 = fD1.GetQ();
502 TVector3 momNeg( (qD1<0?fD1.GetPx():fD2.GetPx()),
503 (qD1<0?fD1.GetPy():fD2.GetPy()),
504 (qD1<0?fD1.GetPz():fD2.GetPz()) );
505 TVector3 momPos( (qD1<0?fD2.GetPx():fD1.GetPx()),
506 (qD1<0?fD2.GetPy():fD1.GetPy()),
507 (qD1<0?fD2.GetPz():fD1.GetPz()) );
508 TVector3 momTot(Px(),Py(),Pz());
510 Double_t lQlNeg = momNeg.Dot(momTot)/momTot.Mag();
511 Double_t lQlPos = momPos.Dot(momTot)/momTot.Mag();
513 return ((lQlPos - lQlNeg)/(lQlPos + lQlNeg));
516 //______________________________________________
517 Double_t AliDielectronPair::GetArmPt() const
520 // Calculate the Armenteros-Podolanski Pt
522 Int_t qD1 = fD1.GetQ();
524 TVector3 momNeg( (qD1<0?fD1.GetPx():fD2.GetPx()),
525 (qD1<0?fD1.GetPy():fD2.GetPy()),
526 (qD1<0?fD1.GetPz():fD2.GetPz()) );
527 TVector3 momTot(Px(),Py(),Pz());
529 return (momNeg.Perp(momTot));
532 // //______________________________________________
533 // Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
536 // // Calculate the decay length in XY taking into account the primary vertex position
538 // if(!vtx) return 0;
539 // return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt() ;
542 // //______________________________________________
543 // Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
546 // // Calculate the pseudo proper time
548 // Double_t lxy=GetLXY(vtx);
549 // Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
550 // return psProperDecayLength;
554 //______________________________________________
555 Double_t AliDielectronPair::PhivPair(Double_t MagField) const
557 //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
558 //to ID conversions. Angle between ee plane and magnetic field is calculated.
560 //Define local buffer variables for leg properties
561 Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
562 Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
602 Double_t px = px1+px2;
603 Double_t py = py1+py2;
604 Double_t pz = pz1+pz2;
605 Double_t dppair = TMath::Sqrt(px*px+py*py+pz*pz);
607 //unit vector of (pep+pem)
608 Double_t pl = dppair;
612 Double_t ax = uy/TMath::Sqrt(ux*ux+uy*uy);
613 Double_t ay = -ux/TMath::Sqrt(ux*ux+uy*uy);
615 //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
617 //Double_t ptep = iep->Px()*ax + iep->Py()*ay;
618 //Double_t ptem = iem->Px()*ax + iem->Py()*ay;
627 //vector product of pep X pem
628 Double_t vpx = pyep*pzem - pzep*pyem;
629 Double_t vpy = pzep*pxem - pxep*pzem;
630 Double_t vpz = pxep*pyem - pyep*pxem;
631 Double_t vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
632 //Double_t thev = acos(vpz/vp);
634 //unit vector of pep X pem
635 Double_t vx = vpx/vp;
636 Double_t vy = vpy/vp;
637 Double_t vz = vpz/vp;
639 //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
640 Double_t wx = uy*vz - uz*vy;
641 Double_t wy = uz*vx - ux*vz;
642 //Double_t wz = ux*vy - uy*vx;
643 //Double_t wl = sqrt(wx*wx+wy*wy+wz*wz);
644 // by construction, (wx,wy,wz) must be a unit vector.
645 // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
646 // should be small if the pair is conversion
648 Double_t cosPhiV = wx*ax + wy*ay;
649 Double_t phiv = TMath::ACos(cosPhiV);
655 //______________________________________________
656 Double_t AliDielectronPair::PairPlaneAngle() const
659 // Calculate the angle between electron pair plane and VZERO-C reaction plane for 2nd harmonic
662 Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
663 Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
674 Double_t px = px1+px2;
675 Double_t py = py1+py2;
676 //Double_t pz = pz1+pz2;
678 // normal vector of ee plane
679 Double_t pnorx = py1*pz2 - pz1*py2;
680 Double_t pnory = pz1*px2 - px1*pz2;
681 Double_t pnorz = px1*py2 - py1*px2;
682 Double_t pnor = TMath::Sqrt( pnorx*pnorx + pnory*pnory + pnorz*pnorz );
685 Double_t upnx = -9999.;
686 Double_t upny = -9999.;
687 Double_t upnz = -9999.;
696 // normal vector of strong magnetic field plane
697 //rotation coordinates (x,y,z)->(x',y',z')
703 Double_t denomHelper = ax*ax + ay*ay +az*az;
704 Double_t uax = -9999.;
705 Double_t uay = -9999.;
706 Double_t uaz = -9999.;
707 if (denomHelper !=0) {
708 uax = ax/TMath::Sqrt(denomHelper);
709 uay = ay/TMath::Sqrt(denomHelper);
710 uaz = az/TMath::Sqrt(denomHelper);
712 //PM is the angle between Pair plane and Magnetic field plane
713 Double_t cosPM = upnx*uax + upny*uay + upnz*uaz;
714 Double_t PM = TMath::ACos(cosPM);
716 //keep interval [0,pi/2]
717 if(PM > TMath::Pi()/2){
726 //______________________________________________
727 void AliDielectronPair::SetBeamEnergy(AliVEvent *ev, Double_t beamEbyHand)
730 // set the beam energy (by hand in case of AODs)
732 if(ev->IsA()==AliESDEvent::Class())
733 fBeamEnergy = ((AliESDEvent*)ev)->GetBeamEnergy();
735 fBeamEnergy = beamEbyHand;