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