Q transverse added
[u/mrichter/AliRoot.git] / ANALYSIS / AliAODPair.cxx
CommitLineData
78d7c6d3 1#include "AliAODPair.h"
2//_________________________________________________________________________
3///////////////////////////////////////////////////////////////////////////
4//
5// class AliAODPair
6//
7// class implements pair of particles and taking care of caluclation (almost)
8// all of pair properties (Qinv, InvMass,...)
9//
10// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
11//
12////////////////////////////////////////////////////////////////////////////
13
14#include "AliVAODParticle.h"
15#include "AliTrackPoints.h"
16
17ClassImp(AliAODPair)
18
19/************************************************************************/
20AliAODPair::AliAODPair(Bool_t rev):
21 fPart1(0x0),
22 fPart2(0x0),
23 fSwappedPair(0x0),
24 fQSideLCMS(0.0),
25 fQSideLCMSNotCalc(kTRUE),
26 fQOutLCMS(0.0),
27 fQOutLCMSNotCalc(kTRUE),
28 fQLongLCMS(0.0),
29 fQLongLCMSNotCalc(kTRUE),
f422a7da 30 fQtLCMS(0.0),
31 fQtLCMSNotCalc(kTRUE),
78d7c6d3 32 fQInv(0.0),
33 fQInvNotCalc(kTRUE),
34 fInvMass(0.0),
35 fInvMassNotCalc(kTRUE),
36 fKt(0.0),
37 fKtNotCalc(kTRUE),
38 fKStar(0.0),
39 fKStarNotCalc(kTRUE),
40 fPInv(0.0),
41 fQSide(0.0),
42 fOut(0.0),
43 fQLong(0.0),
44 fMt(0.0),
45 fMtNotCalc(kTRUE),
46 fInvMassSqr(0.0),
47 fMassSqrNotCalc(kTRUE),
48 fQInvL(0.0),
49 fQInvLNotCalc(kTRUE),
50 fAvarageDistance(0.0),
51 fAvarageDistanceNotCalc(kTRUE),
52 fPxSum(0.0),
53 fPySum(0.0),
54 fPzSum(0.0),
55 fESum(0.0),
56 fSumsNotCalc(kTRUE),
57 fPxDiff(0.0),
58 fPyDiff(0.0),
59 fPzDiff(0.0),
60 fEDiff(0.0),
61 fDiffsNotCalc(kTRUE),
62 fGammaLCMS(0.0),
63 fGammaLCMSNotCalc(kTRUE),
64 fChanged(kTRUE)
65 {
66//value of rev defines if it is Swapped
67//if you pass kTRUE swpaped pair will NOT be created
68//though you wont be able to get the swaped pair from this pair
69
70 if(!rev) fSwappedPair = new AliAODPair(kTRUE); //if false create swaped pair
71
72 }
73/************************************************************************/
74
75AliAODPair::AliAODPair(AliVAODParticle* part1, AliVAODParticle* part2, Bool_t rev):
76 fPart1(part1),
77 fPart2(part2),
78 fSwappedPair(0x0),
79 fQSideLCMS(0.0),
80 fQSideLCMSNotCalc(kTRUE),
81 fQOutLCMS(0.0),
82 fQOutLCMSNotCalc(kTRUE),
83 fQLongLCMS(0.0),
84 fQLongLCMSNotCalc(kTRUE),
f422a7da 85 fQtLCMS(0.0),
86 fQtLCMSNotCalc(kTRUE),
78d7c6d3 87 fQInv(0.0),
88 fQInvNotCalc(kTRUE),
89 fInvMass(0.0),
90 fInvMassNotCalc(kTRUE),
91 fKt(0.0),
92 fKtNotCalc(kTRUE),
93 fKStar(0.0),
94 fKStarNotCalc(kTRUE),
95 fPInv(0.0),
96 fQSide(0.0),
97 fOut(0.0),
98 fQLong(0.0),
99 fMt(0.0),
100 fMtNotCalc(kTRUE),
101 fInvMassSqr(0.0),
102 fMassSqrNotCalc(kTRUE),
103 fQInvL(0.0),
104 fQInvLNotCalc(kTRUE),
105 fAvarageDistance(0.0),
106 fAvarageDistanceNotCalc(kTRUE),
107 fPxSum(0.0),
108 fPySum(0.0),
109 fPzSum(0.0),
110 fESum(0.0),
111 fSumsNotCalc(kTRUE),
112 fPxDiff(0.0),
113 fPyDiff(0.0),
114 fPzDiff(0.0),
115 fEDiff(0.0),
116 fDiffsNotCalc(kTRUE),
117 fGammaLCMS(0.0),
118 fGammaLCMSNotCalc(kTRUE),
119 fChanged(kTRUE)
120 {
121//value of rev defines if it is Swapped
122//if you pass kTRUE swpaped pair will NOT be created
123//though you wont be able to get the swaped pair from this pair
124
125 if(!rev) fSwappedPair = new AliAODPair(part2,part1,kTRUE); //if false create swaped pair
126
127 }
128/************************************************************************/
129AliAODPair::AliAODPair(const AliAODPair& in):
130 TObject(in),
131 fPart1(0x0),
132 fPart2(0x0),
133 fSwappedPair(0x0),
134 fQSideLCMS(0.0),
135 fQSideLCMSNotCalc(kTRUE),
136 fQOutLCMS(0.0),
137 fQOutLCMSNotCalc(kTRUE),
138 fQLongLCMS(0.0),
139 fQLongLCMSNotCalc(kTRUE),
f422a7da 140 fQtLCMS(0.0),
141 fQtLCMSNotCalc(kTRUE),
78d7c6d3 142 fQInv(0.0),
143 fQInvNotCalc(kTRUE),
144 fInvMass(0.0),
145 fInvMassNotCalc(kTRUE),
146 fKt(0.0),
147 fKtNotCalc(kTRUE),
148 fKStar(0.0),
149 fKStarNotCalc(kTRUE),
150 fPInv(0.0),
151 fQSide(0.0),
152 fOut(0.0),
153 fQLong(0.0),
154 fMt(0.0),
155 fMtNotCalc(kTRUE),
156 fInvMassSqr(0.0),
157 fMassSqrNotCalc(kTRUE),
158 fQInvL(0.0),
159 fQInvLNotCalc(kTRUE),
160 fAvarageDistance(0.0),
161 fAvarageDistanceNotCalc(kTRUE),
162 fPxSum(0.0),
163 fPySum(0.0),
164 fPzSum(0.0),
165 fESum(0.0),
166 fSumsNotCalc(kTRUE),
167 fPxDiff(0.0),
168 fPyDiff(0.0),
169 fPzDiff(0.0),
170 fEDiff(0.0),
171 fDiffsNotCalc(kTRUE),
172 fGammaLCMS(0.0),
173 fGammaLCMSNotCalc(kTRUE),
174 fChanged(kTRUE)
175{
176 //cpy constructor
177 in.Copy(*this);
178}
179/************************************************************************/
180
181AliAODPair& AliAODPair::operator=(const AliAODPair& in)
182{
183 //Assigment operator
184 in.Copy(*this);
185 return *this;
186}
187/************************************************************************/
188
189Double_t AliAODPair::GetInvMass()
190{
191//Returns qinv value for a pair
192 if(fInvMassNotCalc)
193 {
194 CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
195
196 if(fInvMassSqr<0) fInvMass = TMath::Sqrt(-fInvMassSqr);
197 else fInvMass = TMath::Sqrt(fInvMassSqr);
198
199 fInvMassNotCalc = kFALSE;
200 }
201 return fInvMass;
202}
203/************************************************************************/
204
205Double_t AliAODPair::GetQSideLCMS()
206{
207//return Q Side in Central Of Mass System in Longitudialy Comoving Frame
208
209 if (fQSideLCMSNotCalc)
210 {
211 fQSideLCMS = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
212 fQSideLCMSNotCalc = kFALSE;
213 }
214 return fQSideLCMS;
215}
216/************************************************************************/
217
218Double_t AliAODPair::GetQOutLCMS()
219{
220 //caculates Qout in Center Of Mass Longitudionally Co-Moving
221 if(fQOutLCMSNotCalc)
222 {
223 CalculateSums();
224 CalculateDiffs();
225
226 if (fPart1->Mass() != fPart2->Mass())
227 {
228/*
229 //STAR algorithm
230 Double_t beta = fPzSum/fESum;
231 Double_t gamma = GetGammaToLCMS();
232 Double_t el = gamma * (fPart1->E() - beta * fPart1->Pz());
233 Double_t x = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
234 beta = 2.0*GetKt()/GetMt();
235 gamma = GetMt()/GetQInv();
236 fQOutLCMS = gamma * (x - beta * el);
237*/
238
239 //beta=fPzSum/fESum; // Longit. V == beta
240 Double_t beta=fPzSum/fESum;
241 Double_t gamma = GetGammaToLCMS();
242
243 Double_t cosphi=fPxSum/(2.0*GetKt()); // cos(phi)
244 Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
245
246// ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
247// ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
248 Double_t tmp;
249 tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
250 Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
251 Double_t part1Px = tmp;
252
253 tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
254 Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
255 Double_t part2Px = tmp;
256
257
258// LTR(part1Pz,E1,beta,GetGammaToLCMS(),part1Pz,E1a);
259// LTR(part2Pz,E2,beta,GetGammaToLCMS(),part2Pz,E2a);
260 Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->E());
261 Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->E());
262
263 Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
264 Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
265 Double_t part1E=TMath::Sqrt(fPart1->Mass()*fPart1->Mass()+part1P2);
266 Double_t part2E=TMath::Sqrt(fPart2->Mass()*fPart2->Mass()+part2P2);
267 Double_t sumE=part1E+part2E;
268 Double_t sumPx=part1Px+part2Px;
269 Double_t sumPy=part1Py+part2Py;
270 Double_t sumPZ=part1Pz+part2Pz;
271 Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
272
273 Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
274 Double_t hf = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/(relmass*relmass);
275 fQOutLCMS=(part1Px-part2Px);//== id
276 fQOutLCMS=fQOutLCMS-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
277 }
278 else
279 {
280 Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
281 fQOutLCMS = 0.5*k2/GetKt();
282 // if (non-id) fQOutLCMS=fQOutLCMS - sumPx*HF;
283 }
284
285
286 fQOutLCMSNotCalc = kFALSE;
287 }
288 return fQOutLCMS;
289}
290/************************************************************************/
291
292Double_t AliAODPair::GetQLongLCMS()
293{
294 //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
295 if (fQLongLCMSNotCalc)
296 {
297 CalculateSums();
298 CalculateDiffs();
299 Double_t beta = fPzSum/fESum;
300 fQLongLCMS = GetGammaToLCMS() * ( fPzDiff - beta*fEDiff );
301 fQLongLCMSNotCalc = kFALSE;
302 }
303 return fQLongLCMS;
304}
305/************************************************************************/
306
f422a7da 307Double_t AliAODPair::GetQtLCMS()
308{
309 //returns Q transverse CMS longitudionally co-moving
310 if (fQtLCMSNotCalc)
311 {
312 fQtLCMS = TMath::Hypot(GetQOutLCMS(),GetQSideLCMS());
313 fQtLCMSNotCalc = kFALSE;
314 }
315 return fQtLCMS;
316}
317
78d7c6d3 318Double_t AliAODPair::GetKt()
319{
320 //calculates the evarage momentum of the pair
321 if(fKtNotCalc)
322 {
323 CalculateSums();
324 fKt = 0.5*TMath::Hypot(fPxSum,fPySum);
325 fKtNotCalc = kFALSE;
326 }
327 return fKt;
328}
329/************************************************************************/
330
331Double_t AliAODPair::GetKStar()
332{
333 //calculates invariant velocity difference
334 if (fKStarNotCalc)
335 {
336 CalculateSums();
337
338 Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
339 Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
340 if (ptrans > mtrans)
341 {
342 Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
343 Error("GetKStar","Particle1:");
344 fPart1->Print();
345 Error("GetKStar","Particle2:");
346 fPart2->Print();
347 Error("GetKStar","");
348
349 fKStar = 10e5;
350 fKStarNotCalc = kFALSE;
351 return fKStar;
352 }
353 Double_t pinv = TMath::Sqrt(mtrans - ptrans);
354
355 Double_t q = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/pinv;
356
357 CalculateQInvL();
358
359 q = q*q - fQInvL;
360 if ( q < 0)
361 {
362 Info("GetKStar","Sqrt of negative number q = %f",q);
363 Error("GetKStar","Particle1:");
364 fPart1->Print();
365 Error("GetKStar","Particle2:");
366 fPart2->Print();
367 fKStar = 10e5;
368 fKStarNotCalc = kFALSE;
369 return fKStar;
370 }
371
372 q = TMath::Sqrt(q);
373 fKStar = q/2.;
374 fKStarNotCalc = kFALSE;
375 }
376 return fKStar;
377}
378/************************************************************************/
379
380Double_t AliAODPair::GetQInv()
381{
382//returns Qinv
383//warning for non-id particles you want to use 2*KStar
384 if(fQInvNotCalc)
385 {
386 CalculateQInvL();
387 fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
388 fQInvNotCalc = kFALSE;
389 }
390 return fQInv;
391}
392/************************************************************************/
393
394Double_t AliAODPair::GetGammaToLCMS()
395{
396 //calculates gamma factor of the boost to LCMS
397 if(fGammaLCMSNotCalc)
398 {
399 CalculateSums();
400 Double_t beta = fPzSum/fESum;
401 fGammaLCMS = 1.0/TMath::Sqrt(1.0 - beta*beta);
402 fGammaLCMSNotCalc = kFALSE;
403 }
404 return fGammaLCMS;
405}
406/************************************************************************/
407
408Double_t AliAODPair::GetMt()
409{
410 //Calculates transverse mass of the pair
411 if (fMtNotCalc)
412 {
413 CalculateSums();
414 fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
415 fMtNotCalc = kFALSE;
416 }
417 return fMt;
418}
419/************************************************************************/
420
421Double_t AliAODPair::GetAvarageDistance()
422{
423//returns and buffers avarage distance between two tracks calculated
424// out of track points (see AliAODTrackPoints class)
425
426 if (fAvarageDistanceNotCalc)
427 {
428 fAvarageDistance = AvDistance();
429 fAvarageDistanceNotCalc = kFALSE;
430 }
431 return fAvarageDistance;
432}
433/************************************************************************/
434
435Double_t AliAODPair::AvDistance()
436{
437 //returns avarage distance between two tracks in range
438 //as defined in Track-Points of AliVAODParticle
439 //returns negative value if error uccured f.g. tracks do not have track-points
440 AliTrackPoints* tpts1 = fPart1->GetTPCTrackPoints();
441 if ( tpts1 == 0x0)
442 {//it could be simulated pair
443// Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
444 return -1.0;
445 }
446
447 AliTrackPoints* tpts2 = fPart2->GetTPCTrackPoints();
448 if ( tpts2 == 0x0)
449 {
450// Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
451 return -1.0;
452 }
453
454 return tpts1->AvarageDistance(*tpts2);
455}