]>
Commit | Line | Data |
---|---|---|
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 | |
32 | ClassImp(AliDielectronPair) | |
33 | ||
1750c463 | 34 | Double_t AliDielectronPair::fBeamEnergy=-1.; |
35 | ||
b2a297fa | 36 | AliDielectronPair::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 | //______________________________________________ | |
54 | AliDielectronPair::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 | //______________________________________________ |
73 | AliDielectronPair::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 | //______________________________________________ |
94 | AliDielectronPair::~AliDielectronPair() | |
95 | { | |
96 | // | |
97 | // Default Destructor | |
98 | // | |
99 | ||
100 | } | |
101 | ||
102 | //______________________________________________ | |
103 | void 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 | //______________________________________________ | |
135 | void 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 | //______________________________________________ |
165 | void 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 | //______________________________________________ |
200 | void 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 | //______________________________________________ |
256 | void AliDielectronPair::GetRotPair(Double_t &RotPairx, Double_t &RotPairy, Double_t &RotPairz) const | |
257 | { | |
258 | // calculation of rotation p1 p2 | |
259 | Double_t px1=-9999.,py1=-9999.,pz1=-9999.; | |
260 | Double_t px2=-9999.,py2=-9999.,pz2=-9999.; | |
261 | ||
262 | px1 = fD1.GetPx(); | |
263 | py1 = fD1.GetPy(); | |
264 | pz1 = fD1.GetPz(); | |
265 | ||
266 | px2 = fD2.GetPx(); | |
267 | py2 = fD2.GetPy(); | |
268 | pz2 = fD2.GetPz(); | |
269 | ||
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 ); | |
275 | ||
276 | //unit vector | |
277 | Double_t upnx = -9999.; | |
278 | Double_t upny = -9999.; | |
279 | Double_t upnz = -9999.; | |
280 | if (pnor !=0) | |
281 | { | |
282 | upnx= pnorx/pnor; | |
283 | upny= pnory/pnor; | |
284 | upnz= pnorz/pnor; | |
285 | } | |
286 | ||
287 | ||
288 | RotPairx = upnx; | |
289 | RotPairy = upny; | |
290 | RotPairz = upnz; | |
291 | ||
292 | } | |
293 | ||
294 | ||
236e1bda | 295 | //______________________________________________ |
296 | Double_t AliDielectronPair::PsiPair(Double_t MagField) const | |
297 | { | |
298 | //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX | |
299 | //to ID conversions. Adapted from AliTRDv0Info class | |
99345a64 | 300 | Double_t x, y;//, z; |
236e1bda | 301 | x = fPair.GetX(); |
302 | y = fPair.GetY(); | |
99345a64 | 303 | // z = fPair.GetZ(); |
236e1bda | 304 | |
305 | Double_t m1[3] = {0,0,0}; | |
306 | Double_t m2[3] = {0,0,0}; | |
307 | ||
308 | m1[0] = fD1.GetPx(); | |
309 | m1[1] = fD1.GetPy(); | |
310 | m1[2] = fD1.GetPz(); | |
311 | ||
312 | m2[0] = fD2.GetPx(); | |
313 | m2[1] = fD2.GetPy(); | |
314 | m2[2] = fD2.GetPz(); | |
315 | ||
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 | |
319 | ||
320 | Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated | |
321 | ||
322 | Double_t mom1Prop[3]; | |
323 | Double_t mom2Prop[3]; | |
324 | ||
325 | AliExternalTrackParam *d1 = static_cast<AliExternalTrackParam*>(fRefD1.GetObject()); | |
326 | AliExternalTrackParam *d2 = static_cast<AliExternalTrackParam*>(fRefD2.GetObject()); | |
327 | ||
328 | AliExternalTrackParam nt(*d1), pt(*d2); | |
236e1bda | 329 | |
330 | Double_t fPsiPair = 4.; | |
331 | if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside | |
332 | fPsiPair = -5.; | |
333 | if(pt.PropagateTo(radiussum,MagField) == 0) | |
334 | fPsiPair = -5.; | |
335 | pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation | |
336 | nt.GetPxPyPz(mom2Prop); | |
337 | ||
338 | ||
339 | ||
340 | Double_t pEle = | |
341 | TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val | |
342 | Double_t pPos = | |
343 | TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val | |
344 | ||
345 | Double_t scalarproduct = | |
346 | mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit | |
347 | ||
348 | Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks | |
349 | ||
350 | fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair)); | |
351 | ||
352 | return fPsiPair; | |
353 | ||
354 | } | |
355 | ||
8df8e382 | 356 | //______________________________________________ |
357 | Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2, | |
61d106d3 | 358 | const Bool_t isHE, const Bool_t isTheta) |
359 | { | |
360 | // The function calculates theta and phi in the mother rest frame with | |
8df8e382 | 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) | |
363 | ||
364 | // Laboratory frame 4-vectors: | |
365 | // projectile beam & target beam 4-mom | |
61d106d3 | 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); | |
b2ad74d0 | 374 | // printf(" beam energy %f \n ", fBeamEnergy); |
1750c463 | 375 | TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass)); |
376 | TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass)); | |
8df8e382 | 377 | |
378 | // first & second daughter 4-mom | |
61d106d3 | 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)); | |
8df8e382 | 381 | // J/Psi 4-momentum vector |
382 | TLorentzVector motherMom=p1Mom+p2Mom; | |
383 | ||
384 | // boost all the 4-mom vectors to the mother rest frame | |
385 | TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect(); | |
386 | p1Mom.Boost(beta); | |
387 | p2Mom.Boost(beta); | |
388 | projMom.Boost(beta); | |
389 | targMom.Boost(beta); | |
390 | ||
391 | // x,y,z axes | |
392 | TVector3 zAxis; | |
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(); | |
397 | ||
398 | // return either theta or phi | |
399 | if(isTheta) { | |
400 | if(d1->Charge()>0) | |
401 | return zAxis.Dot((p1Mom.Vect()).Unit()); | |
402 | else | |
403 | return zAxis.Dot((p2Mom.Vect()).Unit()); | |
404 | ||
405 | } | |
406 | else { | |
407 | if(d1->Charge()>0) | |
408 | return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis)); | |
409 | else | |
410 | return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis)); | |
411 | } | |
412 | } | |
413 | ||
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) | |
419 | ||
420 | // Laboratory frame 4-vectors: | |
421 | // projectile beam & target beam 4-mom | |
45b2b1b8 | 422 | AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject()); |
423 | AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject()); | |
61d106d3 | 424 | |
61d106d3 | 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); | |
433 | ||
1750c463 | 434 | TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass)); |
435 | TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass)); | |
61d106d3 | 436 | |
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)); | |
8df8e382 | 441 | // J/Psi 4-momentum vector |
442 | TLorentzVector motherMom=p1Mom+p2Mom; | |
443 | ||
444 | // boost all the 4-mom vectors to the mother rest frame | |
445 | TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect(); | |
446 | p1Mom.Boost(beta); | |
447 | p2Mom.Boost(beta); | |
448 | projMom.Boost(beta); | |
449 | targMom.Boost(beta); | |
450 | ||
451 | // x,y,z axes | |
452 | TVector3 zAxis; | |
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(); | |
457 | ||
458 | // return either theta or phi | |
459 | if(isTheta) { | |
460 | if(fD1.GetQ()>0) | |
461 | return zAxis.Dot((p1Mom.Vect()).Unit()); | |
462 | else | |
463 | return zAxis.Dot((p2Mom.Vect()).Unit()); | |
464 | } | |
465 | else { | |
466 | if(fD1.GetQ()>0) | |
467 | return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis)); | |
468 | else | |
469 | return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis)); | |
470 | } | |
471 | } | |
2e02dba4 | 472 | //______________________________________________ |
473 | Double_t AliDielectronPair::GetCosPointingAngle(const AliVVertex *primVtx) const | |
474 | { | |
475 | // | |
476 | // Calculate the poiting angle of the pair to the primary vertex and take the cosine | |
477 | // | |
478 | if(!primVtx) return -1.; | |
479 | ||
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(); | |
484 | ||
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]; | |
487 | ||
488 | Double_t cosinePointingAngle = (deltaPos[0]*fPair.GetPx() + deltaPos[1]*fPair.GetPy() + deltaPos[2]*fPair.GetPz()) / TMath::Sqrt(momV02 * deltaPos2); | |
489 | ||
490 | return TMath::Abs(cosinePointingAngle); | |
491 | ||
492 | } | |
ba15fdfb | 493 | |
99345a64 | 494 | //______________________________________________ |
495 | Double_t AliDielectronPair::GetArmAlpha() const | |
496 | { | |
497 | // | |
498 | // Calculate the Armenteros-Podolanski Alpha | |
499 | // | |
500 | Int_t qD1 = fD1.GetQ(); | |
501 | ||
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()); | |
509 | ||
510 | Double_t lQlNeg = momNeg.Dot(momTot)/momTot.Mag(); | |
511 | Double_t lQlPos = momPos.Dot(momTot)/momTot.Mag(); | |
512 | ||
513 | return ((lQlPos - lQlNeg)/(lQlPos + lQlNeg)); | |
514 | } | |
515 | ||
516 | //______________________________________________ | |
517 | Double_t AliDielectronPair::GetArmPt() const | |
518 | { | |
519 | // | |
520 | // Calculate the Armenteros-Podolanski Pt | |
521 | // | |
522 | Int_t qD1 = fD1.GetQ(); | |
523 | ||
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()); | |
528 | ||
529 | return (momNeg.Perp(momTot)); | |
530 | } | |
531 | ||
5720c765 | 532 | // //______________________________________________ |
533 | // Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const | |
534 | // { | |
535 | // // | |
536 | // // Calculate the decay length in XY taking into account the primary vertex position | |
537 | // // | |
538 | // if(!vtx) return 0; | |
539 | // return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt() ; | |
540 | // } | |
ba15fdfb | 541 | |
5720c765 | 542 | // //______________________________________________ |
543 | // Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const | |
544 | // { | |
545 | // // | |
546 | // // Calculate the pseudo proper time | |
547 | // // | |
548 | // Double_t lxy=GetLXY(vtx); | |
549 | // Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt(); | |
550 | // return psProperDecayLength; | |
551 | // } | |
805bd069 | 552 | |
553 | ||
554 | //______________________________________________ | |
555 | Double_t AliDielectronPair::PhivPair(Double_t MagField) const | |
556 | { | |
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. | |
559 | ||
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.; | |
563 | ||
564 | if(MagField>0){ | |
565 | if(fD1.GetQ()>0){ | |
566 | px1 = fD1.GetPx(); | |
567 | py1 = fD1.GetPy(); | |
568 | pz1 = fD1.GetPz(); | |
569 | ||
570 | px2 = fD2.GetPx(); | |
571 | py2 = fD2.GetPy(); | |
572 | pz2 = fD2.GetPz(); | |
573 | }else{ | |
574 | px1 = fD2.GetPx(); | |
575 | py1 = fD2.GetPy(); | |
576 | pz1 = fD2.GetPz(); | |
577 | ||
578 | px2 = fD1.GetPx(); | |
579 | py2 = fD1.GetPy(); | |
580 | pz2 = fD1.GetPz(); | |
581 | } | |
582 | }else{ | |
583 | if(fD1.GetQ()>0){ | |
584 | px1 = fD2.GetPx(); | |
585 | py1 = fD2.GetPy(); | |
586 | pz1 = fD2.GetPz(); | |
587 | ||
588 | px2 = fD1.GetPx(); | |
589 | py2 = fD1.GetPy(); | |
590 | pz2 = fD1.GetPz(); | |
591 | }else{ | |
592 | px1 = fD1.GetPx(); | |
593 | py1 = fD1.GetPy(); | |
594 | pz1 = fD1.GetPz(); | |
595 | ||
596 | px2 = fD2.GetPx(); | |
597 | py2 = fD2.GetPy(); | |
598 | pz2 = fD2.GetPz(); | |
599 | } | |
600 | } | |
601 | ||
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); | |
606 | ||
607 | //unit vector of (pep+pem) | |
608 | Double_t pl = dppair; | |
609 | Double_t ux = px/pl; | |
610 | Double_t uy = py/pl; | |
611 | Double_t uz = pz/pl; | |
612 | Double_t ax = uy/TMath::Sqrt(ux*ux+uy*uy); | |
613 | Double_t ay = -ux/TMath::Sqrt(ux*ux+uy*uy); | |
614 | ||
615 | //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by | |
616 | //definition. | |
617 | //Double_t ptep = iep->Px()*ax + iep->Py()*ay; | |
618 | //Double_t ptem = iem->Px()*ax + iem->Py()*ay; | |
619 | ||
620 | Double_t pxep = px1; | |
621 | Double_t pyep = py1; | |
622 | Double_t pzep = pz1; | |
623 | Double_t pxem = px2; | |
624 | Double_t pyem = py2; | |
625 | Double_t pzem = pz2; | |
626 | ||
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); | |
633 | ||
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; | |
638 | ||
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 | |
647 | // | |
648 | Double_t cosPhiV = wx*ax + wy*ay; | |
649 | Double_t phiv = TMath::ACos(cosPhiV); | |
650 | ||
651 | return phiv; | |
652 | ||
653 | } | |
654 | ||
151ddcbc | 655 | //______________________________________________ |
f3f5888a | 656 | Double_t AliDielectronPair::PairPlaneAngle() const |
151ddcbc | 657 | { |
658 | ||
659 | // Calculate the angle between electron pair plane and VZERO-C reaction plane for 2nd harmonic | |
f3f5888a | 660 | |
151ddcbc | 661 | |
662 | Double_t px1=-9999.,py1=-9999.,pz1=-9999.; | |
663 | Double_t px2=-9999.,py2=-9999.,pz2=-9999.; | |
664 | ||
665 | px1 = fD1.GetPx(); | |
666 | py1 = fD1.GetPy(); | |
667 | pz1 = fD1.GetPz(); | |
668 | ||
669 | px2 = fD2.GetPx(); | |
670 | py2 = fD2.GetPy(); | |
671 | pz2 = fD2.GetPz(); | |
672 | ||
673 | //p1+p2 | |
674 | Double_t px = px1+px2; | |
675 | Double_t py = py1+py2; | |
f3f5888a | 676 | //Double_t pz = pz1+pz2; |
151ddcbc | 677 | |
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 ); | |
683 | ||
684 | //unit vector | |
685 | Double_t upnx = -9999.; | |
686 | Double_t upny = -9999.; | |
07bc5543 | 687 | Double_t upnz = -9999.; |
688 | ||
151ddcbc | 689 | if (pnor !=0) |
690 | { | |
691 | upnx= pnorx/pnor; | |
692 | upny= pnory/pnor; | |
07bc5543 | 693 | upnz= pnorz/pnor; |
151ddcbc | 694 | } |
151ddcbc | 695 | |
696 | // normal vector of strong magnetic field plane | |
697 | //rotation coordinates (x,y,z)->(x',y',z') | |
f3f5888a | 698 | //(p1+p2),(0,0,1) |
699 | Double_t ax = py; | |
700 | Double_t ay = -px; | |
701 | Double_t az = 0.0; | |
702 | ||
696b423a | 703 | Double_t denomHelper = ax*ax + ay*ay +az*az; |
704 | Double_t uax = -9999.; | |
705 | Double_t uay = -9999.; | |
07bc5543 | 706 | Double_t uaz = -9999.; |
696b423a | 707 | if (denomHelper !=0) { |
708 | uax = ax/TMath::Sqrt(denomHelper); | |
709 | uay = ay/TMath::Sqrt(denomHelper); | |
07bc5543 | 710 | uaz = az/TMath::Sqrt(denomHelper); |
696b423a | 711 | } |
151ddcbc | 712 | //PM is the angle between Pair plane and Magnetic field plane |
07bc5543 | 713 | Double_t cosPM = upnx*uax + upny*uay + upnz*uaz; |
151ddcbc | 714 | Double_t PM = TMath::ACos(cosPM); |
715 | ||
716 | //keep interval [0,pi/2] | |
717 | if(PM > TMath::Pi()/2){ | |
718 | PM -= TMath::Pi(); | |
719 | PM *= -1.0; | |
720 | ||
721 | } | |
722 | return PM; | |
723 | } | |
b2ad74d0 | 724 | |
725 | ||
1750c463 | 726 | //______________________________________________ |
727 | void AliDielectronPair::SetBeamEnergy(AliVEvent *ev, Double_t beamEbyHand) | |
728 | { | |
729 | // | |
730 | // set the beam energy (by hand in case of AODs) | |
731 | // | |
732 | if(ev->IsA()==AliESDEvent::Class()) | |
733 | fBeamEnergy = ((AliESDEvent*)ev)->GetBeamEnergy(); | |
734 | else | |
735 | fBeamEnergy = beamEbyHand; | |
736 | } | |
151ddcbc | 737 | |
696b423a | 738 |