]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/AliAODPair.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Access to the ESD objects stored in a tree. Coding conventions (T.Kuhr)
[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"
78d7c6d3 16ClassImp(AliAODPair)
17
18/************************************************************************/
19AliAODPair::AliAODPair(Bool_t rev):
20 fPart1(0x0),
21 fPart2(0x0),
22 fSwappedPair(0x0),
23 fQSideLCMS(0.0),
24 fQSideLCMSNotCalc(kTRUE),
25 fQOutLCMS(0.0),
26 fQOutLCMSNotCalc(kTRUE),
27 fQLongLCMS(0.0),
28 fQLongLCMSNotCalc(kTRUE),
f422a7da 29 fQtLCMS(0.0),
30 fQtLCMSNotCalc(kTRUE),
fb629e2d 31 fQt(0.0),
32 fQtNotCalc(kTRUE),
78d7c6d3 33 fQInv(0.0),
34 fQInvNotCalc(kTRUE),
35 fInvMass(0.0),
36 fInvMassNotCalc(kTRUE),
37 fKt(0.0),
38 fKtNotCalc(kTRUE),
39 fKStar(0.0),
40 fKStarNotCalc(kTRUE),
a40c0433 41 fKStarOut(0.0),
42 fKStarSide(0.0),
43 fKStarLong(0.0),
44 fKStarCompNotCalc(kTRUE),
78d7c6d3 45 fPInv(0.0),
46 fQSide(0.0),
47 fOut(0.0),
48 fQLong(0.0),
49 fMt(0.0),
50 fMtNotCalc(kTRUE),
51 fInvMassSqr(0.0),
52 fMassSqrNotCalc(kTRUE),
53 fQInvL(0.0),
54 fQInvLNotCalc(kTRUE),
55 fAvarageDistance(0.0),
56 fAvarageDistanceNotCalc(kTRUE),
57 fPxSum(0.0),
58 fPySum(0.0),
59 fPzSum(0.0),
60 fESum(0.0),
61 fSumsNotCalc(kTRUE),
62 fPxDiff(0.0),
63 fPyDiff(0.0),
64 fPzDiff(0.0),
65 fEDiff(0.0),
66 fDiffsNotCalc(kTRUE),
67 fGammaLCMS(0.0),
68 fGammaLCMSNotCalc(kTRUE),
69 fChanged(kTRUE)
70 {
71//value of rev defines if it is Swapped
72//if you pass kTRUE swpaped pair will NOT be created
73//though you wont be able to get the swaped pair from this pair
74
75 if(!rev) fSwappedPair = new AliAODPair(kTRUE); //if false create swaped pair
76
77 }
78/************************************************************************/
79
80AliAODPair::AliAODPair(AliVAODParticle* part1, AliVAODParticle* part2, Bool_t rev):
81 fPart1(part1),
82 fPart2(part2),
83 fSwappedPair(0x0),
84 fQSideLCMS(0.0),
85 fQSideLCMSNotCalc(kTRUE),
86 fQOutLCMS(0.0),
87 fQOutLCMSNotCalc(kTRUE),
88 fQLongLCMS(0.0),
89 fQLongLCMSNotCalc(kTRUE),
f422a7da 90 fQtLCMS(0.0),
91 fQtLCMSNotCalc(kTRUE),
78d7c6d3 92 fQInv(0.0),
93 fQInvNotCalc(kTRUE),
94 fInvMass(0.0),
95 fInvMassNotCalc(kTRUE),
96 fKt(0.0),
97 fKtNotCalc(kTRUE),
98 fKStar(0.0),
99 fKStarNotCalc(kTRUE),
a40c0433 100 fKStarOut(0.0),
101 fKStarSide(0.0),
102 fKStarLong(0.0),
103 fKStarCompNotCalc(kTRUE),
78d7c6d3 104 fPInv(0.0),
105 fQSide(0.0),
106 fOut(0.0),
107 fQLong(0.0),
108 fMt(0.0),
109 fMtNotCalc(kTRUE),
110 fInvMassSqr(0.0),
111 fMassSqrNotCalc(kTRUE),
112 fQInvL(0.0),
113 fQInvLNotCalc(kTRUE),
114 fAvarageDistance(0.0),
115 fAvarageDistanceNotCalc(kTRUE),
116 fPxSum(0.0),
117 fPySum(0.0),
118 fPzSum(0.0),
119 fESum(0.0),
120 fSumsNotCalc(kTRUE),
121 fPxDiff(0.0),
122 fPyDiff(0.0),
123 fPzDiff(0.0),
124 fEDiff(0.0),
125 fDiffsNotCalc(kTRUE),
126 fGammaLCMS(0.0),
127 fGammaLCMSNotCalc(kTRUE),
128 fChanged(kTRUE)
129 {
130//value of rev defines if it is Swapped
131//if you pass kTRUE swpaped pair will NOT be created
132//though you wont be able to get the swaped pair from this pair
133
134 if(!rev) fSwappedPair = new AliAODPair(part2,part1,kTRUE); //if false create swaped pair
135
136 }
137/************************************************************************/
138AliAODPair::AliAODPair(const AliAODPair& in):
139 TObject(in),
140 fPart1(0x0),
141 fPart2(0x0),
142 fSwappedPair(0x0),
143 fQSideLCMS(0.0),
144 fQSideLCMSNotCalc(kTRUE),
145 fQOutLCMS(0.0),
146 fQOutLCMSNotCalc(kTRUE),
147 fQLongLCMS(0.0),
148 fQLongLCMSNotCalc(kTRUE),
f422a7da 149 fQtLCMS(0.0),
150 fQtLCMSNotCalc(kTRUE),
78d7c6d3 151 fQInv(0.0),
152 fQInvNotCalc(kTRUE),
153 fInvMass(0.0),
154 fInvMassNotCalc(kTRUE),
155 fKt(0.0),
156 fKtNotCalc(kTRUE),
157 fKStar(0.0),
158 fKStarNotCalc(kTRUE),
a40c0433 159 fKStarOut(0.0),
160 fKStarSide(0.0),
161 fKStarLong(0.0),
162 fKStarCompNotCalc(kTRUE),
78d7c6d3 163 fPInv(0.0),
164 fQSide(0.0),
165 fOut(0.0),
166 fQLong(0.0),
167 fMt(0.0),
168 fMtNotCalc(kTRUE),
169 fInvMassSqr(0.0),
170 fMassSqrNotCalc(kTRUE),
171 fQInvL(0.0),
172 fQInvLNotCalc(kTRUE),
173 fAvarageDistance(0.0),
174 fAvarageDistanceNotCalc(kTRUE),
175 fPxSum(0.0),
176 fPySum(0.0),
177 fPzSum(0.0),
178 fESum(0.0),
179 fSumsNotCalc(kTRUE),
180 fPxDiff(0.0),
181 fPyDiff(0.0),
182 fPzDiff(0.0),
183 fEDiff(0.0),
184 fDiffsNotCalc(kTRUE),
185 fGammaLCMS(0.0),
186 fGammaLCMSNotCalc(kTRUE),
187 fChanged(kTRUE)
188{
189 //cpy constructor
190 in.Copy(*this);
191}
192/************************************************************************/
193
194AliAODPair& AliAODPair::operator=(const AliAODPair& in)
195{
196 //Assigment operator
197 in.Copy(*this);
198 return *this;
199}
200/************************************************************************/
201
202Double_t AliAODPair::GetInvMass()
203{
204//Returns qinv value for a pair
205 if(fInvMassNotCalc)
206 {
207 CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
208
209 if(fInvMassSqr<0) fInvMass = TMath::Sqrt(-fInvMassSqr);
210 else fInvMass = TMath::Sqrt(fInvMassSqr);
211
212 fInvMassNotCalc = kFALSE;
213 }
214 return fInvMass;
215}
216/************************************************************************/
217
218Double_t AliAODPair::GetQSideLCMS()
219{
220//return Q Side in Central Of Mass System in Longitudialy Comoving Frame
221
222 if (fQSideLCMSNotCalc)
223 {
224 fQSideLCMS = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
225 fQSideLCMSNotCalc = kFALSE;
226 }
227 return fQSideLCMS;
228}
229/************************************************************************/
230
231Double_t AliAODPair::GetQOutLCMS()
232{
233 //caculates Qout in Center Of Mass Longitudionally Co-Moving
234 if(fQOutLCMSNotCalc)
235 {
236 CalculateSums();
237 CalculateDiffs();
238
239 if (fPart1->Mass() != fPart2->Mass())
240 {
241/*
242 //STAR algorithm
243 Double_t beta = fPzSum/fESum;
244 Double_t gamma = GetGammaToLCMS();
245 Double_t el = gamma * (fPart1->E() - beta * fPart1->Pz());
246 Double_t x = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
247 beta = 2.0*GetKt()/GetMt();
248 gamma = GetMt()/GetQInv();
249 fQOutLCMS = gamma * (x - beta * el);
250*/
251
252 //beta=fPzSum/fESum; // Longit. V == beta
253 Double_t beta=fPzSum/fESum;
254 Double_t gamma = GetGammaToLCMS();
255
256 Double_t cosphi=fPxSum/(2.0*GetKt()); // cos(phi)
257 Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
258
259// ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
260// ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
261 Double_t tmp;
262 tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
263 Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
264 Double_t part1Px = tmp;
265
266 tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
267 Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
268 Double_t part2Px = tmp;
269
270
271// LTR(part1Pz,E1,beta,GetGammaToLCMS(),part1Pz,E1a);
272// LTR(part2Pz,E2,beta,GetGammaToLCMS(),part2Pz,E2a);
273 Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->E());
274 Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->E());
275
276 Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
277 Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
278 Double_t part1E=TMath::Sqrt(fPart1->Mass()*fPart1->Mass()+part1P2);
279 Double_t part2E=TMath::Sqrt(fPart2->Mass()*fPart2->Mass()+part2P2);
280 Double_t sumE=part1E+part2E;
281 Double_t sumPx=part1Px+part2Px;
282 Double_t sumPy=part1Py+part2Py;
283 Double_t sumPZ=part1Pz+part2Pz;
284 Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
285
286 Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
287 Double_t hf = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/(relmass*relmass);
288 fQOutLCMS=(part1Px-part2Px);//== id
289 fQOutLCMS=fQOutLCMS-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
290 }
291 else
292 {
293 Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
294 fQOutLCMS = 0.5*k2/GetKt();
295 // if (non-id) fQOutLCMS=fQOutLCMS - sumPx*HF;
296 }
297
298
299 fQOutLCMSNotCalc = kFALSE;
300 }
301 return fQOutLCMS;
302}
303/************************************************************************/
304
305Double_t AliAODPair::GetQLongLCMS()
306{
307 //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
308 if (fQLongLCMSNotCalc)
309 {
310 CalculateSums();
311 CalculateDiffs();
312 Double_t beta = fPzSum/fESum;
313 fQLongLCMS = GetGammaToLCMS() * ( fPzDiff - beta*fEDiff );
314 fQLongLCMSNotCalc = kFALSE;
315 }
316 return fQLongLCMS;
317}
318/************************************************************************/
319
f422a7da 320Double_t AliAODPair::GetQtLCMS()
321{
322 //returns Q transverse CMS longitudionally co-moving
323 if (fQtLCMSNotCalc)
324 {
325 fQtLCMS = TMath::Hypot(GetQOutLCMS(),GetQSideLCMS());
326 fQtLCMSNotCalc = kFALSE;
327 }
328 return fQtLCMS;
329}
fb629e2d 330/************************************************************************/
331
332Double_t AliAODPair::GetQt()
333{
334 //returns Q transverse CMS longitudionally co-moving
335 if (fQtNotCalc)
336 {
05eb9603 337 CalculateSums();
338 CalculateDiffs();
339
fb629e2d 340 Double_t dotprod = fPxSum*fPxDiff + fPySum*fPyDiff + fPzSum*fPzDiff;
341 Double_t klen = fPxSum*fPxSum + fPySum*fPySum + fPzSum*fPzSum;
342 klen = TMath::Sqrt(klen);
343 Double_t qlen = fPxDiff*fPxDiff + fPyDiff*fPyDiff + fPzDiff*fPzDiff;
344 qlen = TMath::Sqrt(qlen);
05eb9603 345 Double_t tmp = klen*qlen;
346 if (tmp == 0.0)
347 {
348 fQt = 10e5;
349 fQtNotCalc = kFALSE;
350 return fQt;
351 }
352 Double_t cosopenangle = dotprod/tmp;
fb629e2d 353 Double_t sinopenangle = TMath::Sqrt(1.0 - cosopenangle*cosopenangle);
354
355 fQt = sinopenangle*qlen;
356 fQtNotCalc = kFALSE;
357 }
358 return fQt;
359}
360/************************************************************************/
f422a7da 361
78d7c6d3 362Double_t AliAODPair::GetKt()
363{
364 //calculates the evarage momentum of the pair
365 if(fKtNotCalc)
366 {
367 CalculateSums();
368 fKt = 0.5*TMath::Hypot(fPxSum,fPySum);
369 fKtNotCalc = kFALSE;
370 }
371 return fKt;
372}
373/************************************************************************/
374
375Double_t AliAODPair::GetKStar()
376{
377 //calculates invariant velocity difference
378 if (fKStarNotCalc)
379 {
380 CalculateSums();
381
382 Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
383 Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
384 if (ptrans > mtrans)
385 {
386 Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
387 Error("GetKStar","Particle1:");
388 fPart1->Print();
389 Error("GetKStar","Particle2:");
390 fPart2->Print();
391 Error("GetKStar","");
392
393 fKStar = 10e5;
394 fKStarNotCalc = kFALSE;
395 return fKStar;
396 }
397 Double_t pinv = TMath::Sqrt(mtrans - ptrans);
398
399 Double_t q = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/pinv;
400
401 CalculateQInvL();
402
403 q = q*q - fQInvL;
404 if ( q < 0)
405 {
406 Info("GetKStar","Sqrt of negative number q = %f",q);
407 Error("GetKStar","Particle1:");
408 fPart1->Print();
409 Error("GetKStar","Particle2:");
410 fPart2->Print();
411 fKStar = 10e5;
412 fKStarNotCalc = kFALSE;
413 return fKStar;
414 }
415
416 q = TMath::Sqrt(q);
417 fKStar = q/2.;
418 fKStarNotCalc = kFALSE;
419 }
420 return fKStar;
421}
422/************************************************************************/
a40c0433 423Double_t AliAODPair::GetKStarOut()
424{
425 CalculateKStarComp();
426 return fKStarOut;
427}
428/************************************************************************/
429Double_t AliAODPair::GetKStarSide()
430{
431 CalculateKStarComp();
432 return fKStarSide;
433}
434/************************************************************************/
435Double_t AliAODPair::GetKStarLong()
436{
437 CalculateKStarComp();
438 return fKStarLong;
439}
440/************************************************************************/
78d7c6d3 441
442Double_t AliAODPair::GetQInv()
443{
444//returns Qinv
445//warning for non-id particles you want to use 2*KStar
446 if(fQInvNotCalc)
447 {
448 CalculateQInvL();
449 fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
450 fQInvNotCalc = kFALSE;
451 }
452 return fQInv;
453}
454/************************************************************************/
455
456Double_t AliAODPair::GetGammaToLCMS()
457{
458 //calculates gamma factor of the boost to LCMS
459 if(fGammaLCMSNotCalc)
460 {
461 CalculateSums();
462 Double_t beta = fPzSum/fESum;
463 fGammaLCMS = 1.0/TMath::Sqrt(1.0 - beta*beta);
464 fGammaLCMSNotCalc = kFALSE;
465 }
466 return fGammaLCMS;
467}
468/************************************************************************/
469
470Double_t AliAODPair::GetMt()
471{
472 //Calculates transverse mass of the pair
473 if (fMtNotCalc)
474 {
475 CalculateSums();
476 fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
477 fMtNotCalc = kFALSE;
478 }
479 return fMt;
480}
481/************************************************************************/
482
483Double_t AliAODPair::GetAvarageDistance()
484{
485//returns and buffers avarage distance between two tracks calculated
486// out of track points (see AliAODTrackPoints class)
487
488 if (fAvarageDistanceNotCalc)
489 {
490 fAvarageDistance = AvDistance();
491 fAvarageDistanceNotCalc = kFALSE;
492 }
493 return fAvarageDistance;
494}
495/************************************************************************/
496
497Double_t AliAODPair::AvDistance()
498{
499 //returns avarage distance between two tracks in range
500 //as defined in Track-Points of AliVAODParticle
501 //returns negative value if error uccured f.g. tracks do not have track-points
502 AliTrackPoints* tpts1 = fPart1->GetTPCTrackPoints();
503 if ( tpts1 == 0x0)
504 {//it could be simulated pair
505// Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
506 return -1.0;
507 }
508
509 AliTrackPoints* tpts2 = fPart2->GetTPCTrackPoints();
510 if ( tpts2 == 0x0)
511 {
512// Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
513 return -1.0;
514 }
515
516 return tpts1->AvarageDistance(*tpts2);
517}
8ea30edc 518/************************************************************************/
519
520Double_t AliAODPair::GetR()
521{
522//Returns distance between particles vertexes in thir CMS
523
524 CalculateDiffs();
5c03edbe 525
526 Double_t vxDiff = fPart1->Vx() - fPart2->Vx();
527 Double_t vyDiff = fPart1->Vy() - fPart2->Vy();
528 Double_t vzDiff = fPart1->Vz() - fPart2->Vz();
529
530 Double_t r = TMath::Sqrt( vxDiff*vxDiff + vyDiff*vyDiff + vzDiff*vzDiff );
531 return r;
8ea30edc 532
533}
534/************************************************************************/
535
536Double_t AliAODPair::GetRStar()
537{
538//Returns distance between particles vertexes in thir CMS
539
540
541 CalculateSums();
542
543 Double_t klen = fPxSum*fPxSum + fPySum*fPySum + fPzSum*fPzSum;
544 klen = TMath::Sqrt(klen);
545
546 Double_t aBeta = klen/fESum;
547 Double_t aGamma = 1.0/TMath::Sqrt(1.0 - aBeta*aBeta);
548
549
550 Double_t alpha = -TMath::ATan2(fPySum,fPzSum);
551 Double_t beta = TMath::ATan2(fPxSum,TMath::Hypot(fPySum,fPzSum));
552
553 Double_t sinalpha = TMath::Sin(alpha);
554 Double_t cosalpha = TMath::Cos(alpha);
555 Double_t sinbeta = TMath::Sin(beta);
556 Double_t cosbeta = TMath::Cos(beta);
557
558 Double_t v1xP = fPart1->Vx();
559 Double_t v2xP = fPart2->Vx();
560 Double_t v1yP = fPart1->Vy()*cosalpha + fPart1->Vz()*sinalpha;
561 Double_t v2yP = fPart2->Vy()*cosalpha + fPart2->Vz()*sinalpha;
562 Double_t v1zP =-fPart1->Vy()*sinalpha + fPart1->Vz()*cosalpha;
563 Double_t v2zP =-fPart2->Vy()*sinalpha + fPart2->Vz()*cosalpha;
564
565
566///////////////////////////////////////////////////
567
568// Double_t p1yP = fPart1->Py()*cosalpha + fPart1->Pz()*sinalpha;
569// Double_t p2yP = fPart2->Py()*cosalpha + fPart2->Pz()*sinalpha;
570//
571// Double_t p1zP =-fPart1->Py()*sinalpha + fPart1->Pz()*cosalpha;
572// Double_t p2zP =-fPart2->Py()*sinalpha + fPart2->Pz()*cosalpha;
573//
574//
575// Double_t p1x = fPart1->Px()*cosbeta - p1zP*sinbeta;
576// Double_t p2x = fPart2->Px()*cosbeta - p2zP*sinbeta;
577// Double_t p1z = fPart1->Px()*sinbeta + p1zP*cosbeta;
578// Double_t p2z = fPart2->Px()*sinbeta + p2zP*cosbeta;
579
580// Info("","%f %f %f",p1yP,p2yP,p1yP+p2yP);
581// Info("","%f %f %f",p1x,p2x,p1x+p2x);
582
583// Info("","%f %f ",p1x+p2x,p1yP+p2yP);
584
585///////////////////////////////////////////////////
586
587
588 Double_t v1x = v1xP*cosbeta - v1zP*sinbeta;
589 Double_t v2x = v2xP*cosbeta - v2zP*sinbeta;
590 Double_t v1y = v1yP;
591 Double_t v2y = v2yP;
592 Double_t v1z = v1xP*sinbeta + v1zP*cosbeta;
593 Double_t v2z = v2xP*sinbeta + v2zP*cosbeta;
594
595
596 Double_t v1zB=aGamma*(v1z-aBeta*fPart1->T());
597 Double_t v2zB=aGamma*(v2z-aBeta*fPart2->T());
598
599
600
601 Double_t dx = v1x - v2x;
602 Double_t dy = v1y - v2y;
603 Double_t dz = v1zB - v2zB;
604
605 Double_t rstar = TMath::Sqrt( dx*dx + dy*dy + dz*dz);
606
607 return rstar;
608}
6808c13f 609/************************************************************************/
af2872da 610
611void AliAODPair::MirrorSecond()
612{
613//makes local copy of the second particle and mirrors their momenta
614//for its deletion is responsible who calls this method
615 fPart2 = (AliVAODParticle*)fPart2->Clone();
616 fPart2->SetMomentum(-fPart2->Px(),-fPart2->Py(),-fPart2->Pz(),fPart2->E());
6808c13f 617 Changed();
af2872da 618}
6808c13f 619/************************************************************************/
af2872da 620
621void AliAODPair::DeleteSecond()
622{
623//Deletes second particle
624 delete fPart2;
625 fPart2 = 0x0;
626}