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