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