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