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