major dielectron update (included also the data and plotting macros for paper)
[u/mrichter/AliRoot.git] / PWG3 / 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>
27
b2a297fa 28#include "AliDielectronPair.h"
b2a297fa 29
30ClassImp(AliDielectronPair)
31
32AliDielectronPair::AliDielectronPair() :
b2a297fa 33 fType(-1),
a655b716 34 fLabel(-1),
b2a297fa 35 fPair(),
572b0139 36 fD1(),
37 fD2(),
b2a297fa 38 fRefD1(),
39 fRefD2()
40{
41 //
42 // Default Constructor
43 //
44
45}
46
47//______________________________________________
48AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
49 AliVTrack * const particle2, Int_t pid2, Char_t type) :
b2a297fa 50 fType(type),
a655b716 51 fLabel(-1),
b2a297fa 52 fPair(),
572b0139 53 fD1(),
54 fD2(),
b2a297fa 55 fRefD1(),
56 fRefD2()
57{
58 //
59 // Constructor with tracks
60 //
61 SetTracks(particle1, pid1, particle2, pid2);
62}
63
64//______________________________________________
1201a1a9 65AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
66 const AliKFParticle * const particle2,
67 AliVTrack * const refParticle1,
68 AliVTrack * const refParticle2, Char_t type) :
69 fType(type),
70 fLabel(-1),
71 fPair(),
72 fD1(),
73 fD2(),
74 fRefD1(),
75 fRefD2()
76{
77 //
78 // Constructor with tracks
79 //
80 SetTracks(particle1, particle2,refParticle1,refParticle2);
81}
82
83//______________________________________________
b2a297fa 84AliDielectronPair::~AliDielectronPair()
85{
86 //
87 // Default Destructor
88 //
89
90}
91
92//______________________________________________
93void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
94 AliVTrack * const particle2, Int_t pid2)
95{
96 //
572b0139 97 // Sort particles by pt, first particle larget Pt
98 // set AliKF daughters and pair
1201a1a9 99 // refParticle1 and 2 are the original tracks. In the case of track rotation
100 // they are needed in the framework
b2a297fa 101 //
102 fPair.Initialize();
572b0139 103 fD1.Initialize();
104 fD2.Initialize();
8df8e382 105
b2a297fa 106 AliKFParticle kf1(*particle1,pid1);
107 AliKFParticle kf2(*particle2,pid2);
572b0139 108
b2a297fa 109 fPair.AddDaughter(kf1);
110 fPair.AddDaughter(kf2);
8df8e382 111
a655b716 112 if (particle1->Pt()>particle2->Pt()){
113 fRefD1 = particle1;
114 fRefD2 = particle2;
572b0139 115 fD1+=kf1;
116 fD2+=kf2;
a655b716 117 } else {
118 fRefD1 = particle2;
119 fRefD2 = particle1;
572b0139 120 fD1+=kf2;
121 fD2+=kf1;
a655b716 122 }
b2a297fa 123}
124
8df8e382 125//______________________________________________
1201a1a9 126void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
127 const AliKFParticle * const particle2,
128 AliVTrack * const refParticle1,
129 AliVTrack * const refParticle2)
130{
131 //
132 // Sort particles by pt, first particle larget Pt
133 // set AliKF daughters and pair
134 // refParticle1 and 2 are the original tracks. In the case of track rotation
135 // they are needed in the framework
136 //
137 fPair.Initialize();
138 fD1.Initialize();
139 fD2.Initialize();
140
141 AliKFParticle kf1(*particle1);
142 AliKFParticle kf2(*particle2);
143
144 fPair.AddDaughter(kf1);
145 fPair.AddDaughter(kf2);
146
147 if (kf1.GetPt()>kf2.GetPt()){
148 fRefD1 = refParticle1;
149 fRefD2 = refParticle2;
150 fD1+=kf1;
151 fD2+=kf2;
152 } else {
153 fRefD1 = refParticle2;
154 fRefD2 = refParticle1;
155 fD1+=kf2;
156 fD2+=kf1;
157 }
158}
159
160//______________________________________________
61d106d3 161void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
162{
163 //
164 // Calculate theta and phi in helicity and Collins-Soper coordinate frame
165 //
166 const Double_t kBeamEnergy = 3500.;
1201a1a9 167 Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
168 Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
61d106d3 169 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
170 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
171
1201a1a9 172// AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
173// AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
61d106d3 174
1201a1a9 175// d1->PxPyPz(pxyz1);
176// d2->PxPyPz(pxyz2);
61d106d3 177
178 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
179 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
180
181 // first & second daughter 4-mom
182 TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
183 TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
184 TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
185 TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
186 // J/Psi 4-momentum vector
187 TLorentzVector motherMom=p1Mom+p2Mom;
188
189 // boost all the 4-mom vectors to the mother rest frame
190 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
191 p1Mom.Boost(beta);
192 p2Mom.Boost(beta);
193 projMom.Boost(beta);
194 targMom.Boost(beta);
195
196 // x,y,z axes
197 TVector3 zAxisHE = (motherMom.Vect()).Unit();
198 TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
199 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
200 TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
201 TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
202
203 // fill theta and phi
1201a1a9 204 if(fD1.GetQ()>0){
61d106d3 205 thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
206 thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
207 phiHE = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
208 phiCS = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
209 } else {
210 thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
211 thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
212 phiHE = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
213 phiCS = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
214 }
215}
216
217//______________________________________________
8df8e382 218Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
61d106d3 219 const Bool_t isHE, const Bool_t isTheta)
220{
221 // The function calculates theta and phi in the mother rest frame with
8df8e382 222 // respect to the helicity coordinate system and Collins-Soper coordinate system
223 // TO DO: generalize for different decays (only J/Psi->e+e- now)
224
225 // Laboratory frame 4-vectors:
226 // projectile beam & target beam 4-mom
61d106d3 227 // TODO: need to retrieve the beam energy from somewhere
228 const Double_t kBeamEnergy = 3500.;
229 Double_t px1=d1->Px();
230 Double_t py1=d1->Py();
231 Double_t pz1=d1->Pz();
232 Double_t px2=d2->Px();
233 Double_t py2=d2->Py();
234 Double_t pz2=d2->Pz();
235 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
236 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
237
238 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
239 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
8df8e382 240
241 // first & second daughter 4-mom
61d106d3 242 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
243 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
8df8e382 244 // J/Psi 4-momentum vector
245 TLorentzVector motherMom=p1Mom+p2Mom;
246
247 // boost all the 4-mom vectors to the mother rest frame
248 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
249 p1Mom.Boost(beta);
250 p2Mom.Boost(beta);
251 projMom.Boost(beta);
252 targMom.Boost(beta);
253
254 // x,y,z axes
255 TVector3 zAxis;
256 if(isHE) zAxis = (motherMom.Vect()).Unit();
257 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
258 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
259 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
260
261 // return either theta or phi
262 if(isTheta) {
263 if(d1->Charge()>0)
264 return zAxis.Dot((p1Mom.Vect()).Unit());
265 else
266 return zAxis.Dot((p2Mom.Vect()).Unit());
267
268 }
269 else {
270 if(d1->Charge()>0)
271 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
272 else
273 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
274 }
275}
276
277//______________________________________________
278Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
279 // The function calculates theta and phi in the mother rest frame with
280 // respect to the helicity coordinate system and Collins-Soper coordinate system
281 // TO DO: generalize for different decays (only J/Psi->e+e- now)
282
283 // Laboratory frame 4-vectors:
284 // projectile beam & target beam 4-mom
45b2b1b8 285 AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
286 AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
61d106d3 287
288 const Double_t kBeamEnergy = 3500.;
289 Double_t px1=d1->Px();
290 Double_t py1=d1->Py();
291 Double_t pz1=d1->Pz();
292 Double_t px2=d2->Px();
293 Double_t py2=d2->Py();
294 Double_t pz2=d2->Pz();
295 Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
296 Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
297
298 TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
299 TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
300
301 // first & second daughter 4-mom
302 // first & second daughter 4-mom
303 TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
304 TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
8df8e382 305 // J/Psi 4-momentum vector
306 TLorentzVector motherMom=p1Mom+p2Mom;
307
308 // boost all the 4-mom vectors to the mother rest frame
309 TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
310 p1Mom.Boost(beta);
311 p2Mom.Boost(beta);
312 projMom.Boost(beta);
313 targMom.Boost(beta);
314
315 // x,y,z axes
316 TVector3 zAxis;
317 if(isHE) zAxis = (motherMom.Vect()).Unit();
318 else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
319 TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
320 TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
321
322 // return either theta or phi
323 if(isTheta) {
324 if(fD1.GetQ()>0)
325 return zAxis.Dot((p1Mom.Vect()).Unit());
326 else
327 return zAxis.Dot((p2Mom.Vect()).Unit());
328 }
329 else {
330 if(fD1.GetQ()>0)
331 return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
332 else
333 return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
334 }
335}
ba15fdfb 336
337//______________________________________________
338Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
339{
340 //
341 // Calculate the decay length in XY taking into account the primary vertex position
342 //
343 if(!vtx) return 0;
344 return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt() ;
345}
346
347//______________________________________________
348Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
349{
350 //
351 // Calculate the pseudo proper time
352 //
353 Double_t lxy=GetLXY(vtx);
354 Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
355 return psProperDecayLength;
356}