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