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>
29 #include "AliDielectronPair.h"
31 ClassImp(AliDielectronPair)
33 AliDielectronPair::AliDielectronPair() :
43 // Default Constructor
48 //______________________________________________
49 AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
50 AliVTrack * const particle2, Int_t pid2, Char_t type) :
60 // Constructor with tracks
62 SetTracks(particle1, pid1, particle2, pid2);
65 //______________________________________________
66 AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
67 const AliKFParticle * const particle2,
68 AliVTrack * const refParticle1,
69 AliVTrack * const refParticle2, Char_t type) :
79 // Constructor with tracks
81 SetTracks(particle1, particle2,refParticle1,refParticle2);
84 //______________________________________________
85 AliDielectronPair::~AliDielectronPair()
93 //______________________________________________
94 void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
95 AliVTrack * const particle2, Int_t pid2)
98 // Sort particles by pt, first particle larget Pt
99 // set AliKF daughters and pair
100 // refParticle1 and 2 are the original tracks. In the case of track rotation
101 // they are needed in the framework
107 AliKFParticle kf1(*particle1,pid1);
108 AliKFParticle kf2(*particle2,pid2);
110 fPair.AddDaughter(kf1);
111 fPair.AddDaughter(kf2);
113 if (particle1->Pt()>particle2->Pt()){
126 //______________________________________________
127 void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
128 const AliKFParticle * const particle2,
129 AliVTrack * const refParticle1,
130 AliVTrack * const refParticle2)
133 // Sort particles by pt, first particle larget Pt
134 // set AliKF daughters and pair
135 // refParticle1 and 2 are the original tracks. In the case of track rotation
136 // they are needed in the framework
142 AliKFParticle kf1(*particle1);
143 AliKFParticle kf2(*particle2);
145 fPair.AddDaughter(kf1);
146 fPair.AddDaughter(kf2);
148 if (kf1.GetPt()>kf2.GetPt()){
149 fRefD1 = refParticle1;
150 fRefD2 = refParticle2;
154 fRefD1 = refParticle2;
155 fRefD2 = refParticle1;
161 //______________________________________________
162 void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
165 // Calculate theta and phi in helicity and Collins-Soper coordinate frame
167 const Double_t kBeamEnergy = 3500.;
168 Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
169 Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
170 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
171 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
173 // AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
174 // AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
176 // d1->PxPyPz(pxyz1);
177 // d2->PxPyPz(pxyz2);
179 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
180 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
182 // first & second daughter 4-mom
183 TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
184 TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
185 TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
186 TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
187 // J/Psi 4-momentum vector
188 TLorentzVector motherMom=p1Mom+p2Mom;
190 // boost all the 4-mom vectors to the mother rest frame
191 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
198 TVector3 zAxisHE = (motherMom.Vect()).Unit();
199 TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
200 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
201 TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
202 TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
204 // fill theta and phi
206 thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
207 thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
208 phiHE = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
209 phiCS = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
211 thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
212 thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
213 phiHE = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
214 phiCS = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
218 //______________________________________________
219 Double_t AliDielectronPair::PsiPair(Double_t MagField) const
221 //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
222 //to ID conversions. Adapted from AliTRDv0Info class
228 Double_t m1[3] = {0,0,0};
229 Double_t m2[3] = {0,0,0};
239 Double_t deltat = 1.;
240 deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]*m2[0] + m2[1]*m2[1])+1.e-13))-
241 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
243 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
245 Double_t mom1Prop[3];
246 Double_t mom2Prop[3];
248 AliExternalTrackParam *d1 = static_cast<AliExternalTrackParam*>(fRefD1.GetObject());
249 AliExternalTrackParam *d2 = static_cast<AliExternalTrackParam*>(fRefD2.GetObject());
251 AliExternalTrackParam nt(*d1), pt(*d2);
252 //AliExternalTrackParam nt(), pt();
253 nt.CopyFromVTrack(d1);
254 pt.CopyFromVTrack(d2);
256 Double_t fPsiPair = 4.;
257 if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
259 if(pt.PropagateTo(radiussum,MagField) == 0)
261 pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
262 nt.GetPxPyPz(mom2Prop);
267 TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
269 TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val
271 Double_t scalarproduct =
272 mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit
274 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
276 fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
282 //______________________________________________
283 Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
284 const Bool_t isHE, const Bool_t isTheta)
286 // The function calculates theta and phi in the mother rest frame with
287 // respect to the helicity coordinate system and Collins-Soper coordinate system
288 // TO DO: generalize for different decays (only J/Psi->e+e- now)
290 // Laboratory frame 4-vectors:
291 // projectile beam & target beam 4-mom
292 // TODO: need to retrieve the beam energy from somewhere
293 const Double_t kBeamEnergy = 3500.;
294 Double_t px1=d1->Px();
295 Double_t py1=d1->Py();
296 Double_t pz1=d1->Pz();
297 Double_t px2=d2->Px();
298 Double_t py2=d2->Py();
299 Double_t pz2=d2->Pz();
300 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
301 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
303 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
304 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
306 // first & second daughter 4-mom
307 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
308 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
309 // J/Psi 4-momentum vector
310 TLorentzVector motherMom=p1Mom+p2Mom;
312 // boost all the 4-mom vectors to the mother rest frame
313 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
321 if(isHE) zAxis = (motherMom.Vect()).Unit();
322 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
323 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
324 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
326 // return either theta or phi
329 return zAxis.Dot((p1Mom.Vect()).Unit());
331 return zAxis.Dot((p2Mom.Vect()).Unit());
336 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
338 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
342 //______________________________________________
343 Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
344 // The function calculates theta and phi in the mother rest frame with
345 // respect to the helicity coordinate system and Collins-Soper coordinate system
346 // TO DO: generalize for different decays (only J/Psi->e+e- now)
348 // Laboratory frame 4-vectors:
349 // projectile beam & target beam 4-mom
350 AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
351 AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
353 const Double_t kBeamEnergy = 3500.;
354 Double_t px1=d1->Px();
355 Double_t py1=d1->Py();
356 Double_t pz1=d1->Pz();
357 Double_t px2=d2->Px();
358 Double_t py2=d2->Py();
359 Double_t pz2=d2->Pz();
360 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
361 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
363 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
364 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
366 // first & second daughter 4-mom
367 // first & second daughter 4-mom
368 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
369 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
370 // J/Psi 4-momentum vector
371 TLorentzVector motherMom=p1Mom+p2Mom;
373 // boost all the 4-mom vectors to the mother rest frame
374 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
382 if(isHE) zAxis = (motherMom.Vect()).Unit();
383 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
384 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
385 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
387 // return either theta or phi
390 return zAxis.Dot((p1Mom.Vect()).Unit());
392 return zAxis.Dot((p2Mom.Vect()).Unit());
396 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
398 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
402 // //______________________________________________
403 // Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
406 // // Calculate the decay length in XY taking into account the primary vertex position
408 // if(!vtx) return 0;
409 // return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt() ;
412 // //______________________________________________
413 // Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
416 // // Calculate the pseudo proper time
418 // Double_t lxy=GetLXY(vtx);
419 // Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
420 // return psProperDecayLength;