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