]> 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
Distance between creation points CFs
[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();
525 return TMath::Sqrt( fPxDiff*fPxDiff + fPyDiff*fPyDiff + fPzDiff*fPzDiff);
526
527}
528/************************************************************************/
529
530Double_t AliAODPair::GetRStar()
531{
532//Returns distance between particles vertexes in thir CMS
533
534
535 CalculateSums();
536
537 Double_t klen = fPxSum*fPxSum + fPySum*fPySum + fPzSum*fPzSum;
538 klen = TMath::Sqrt(klen);
539
540 Double_t aBeta = klen/fESum;
541 Double_t aGamma = 1.0/TMath::Sqrt(1.0 - aBeta*aBeta);
542
543
544 Double_t alpha = -TMath::ATan2(fPySum,fPzSum);
545 Double_t beta = TMath::ATan2(fPxSum,TMath::Hypot(fPySum,fPzSum));
546
547 Double_t sinalpha = TMath::Sin(alpha);
548 Double_t cosalpha = TMath::Cos(alpha);
549 Double_t sinbeta = TMath::Sin(beta);
550 Double_t cosbeta = TMath::Cos(beta);
551
552 Double_t v1xP = fPart1->Vx();
553 Double_t v2xP = fPart2->Vx();
554 Double_t v1yP = fPart1->Vy()*cosalpha + fPart1->Vz()*sinalpha;
555 Double_t v2yP = fPart2->Vy()*cosalpha + fPart2->Vz()*sinalpha;
556 Double_t v1zP =-fPart1->Vy()*sinalpha + fPart1->Vz()*cosalpha;
557 Double_t v2zP =-fPart2->Vy()*sinalpha + fPart2->Vz()*cosalpha;
558
559
560///////////////////////////////////////////////////
561
562// Double_t p1yP = fPart1->Py()*cosalpha + fPart1->Pz()*sinalpha;
563// Double_t p2yP = fPart2->Py()*cosalpha + fPart2->Pz()*sinalpha;
564//
565// Double_t p1zP =-fPart1->Py()*sinalpha + fPart1->Pz()*cosalpha;
566// Double_t p2zP =-fPart2->Py()*sinalpha + fPart2->Pz()*cosalpha;
567//
568//
569// Double_t p1x = fPart1->Px()*cosbeta - p1zP*sinbeta;
570// Double_t p2x = fPart2->Px()*cosbeta - p2zP*sinbeta;
571// Double_t p1z = fPart1->Px()*sinbeta + p1zP*cosbeta;
572// Double_t p2z = fPart2->Px()*sinbeta + p2zP*cosbeta;
573
574// Info("","%f %f %f",p1yP,p2yP,p1yP+p2yP);
575// Info("","%f %f %f",p1x,p2x,p1x+p2x);
576
577// Info("","%f %f ",p1x+p2x,p1yP+p2yP);
578
579///////////////////////////////////////////////////
580
581
582 Double_t v1x = v1xP*cosbeta - v1zP*sinbeta;
583 Double_t v2x = v2xP*cosbeta - v2zP*sinbeta;
584 Double_t v1y = v1yP;
585 Double_t v2y = v2yP;
586 Double_t v1z = v1xP*sinbeta + v1zP*cosbeta;
587 Double_t v2z = v2xP*sinbeta + v2zP*cosbeta;
588
589
590 Double_t v1zB=aGamma*(v1z-aBeta*fPart1->T());
591 Double_t v2zB=aGamma*(v2z-aBeta*fPart2->T());
592
593
594
595 Double_t dx = v1x - v2x;
596 Double_t dy = v1y - v2y;
597 Double_t dz = v1zB - v2zB;
598
599 Double_t rstar = TMath::Sqrt( dx*dx + dy*dy + dz*dz);
600
601 return rstar;
602}