31de69be0aaa4e9b905983835351d0b6d5a2f408
[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  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
73 AliAODPair::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 /************************************************************************/
125 AliAODPair::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
175 AliAODPair& AliAODPair::operator=(const AliAODPair& in)
176 {
177  //Assigment operator
178  in.Copy(*this);
179  return *this;
180 }
181 /************************************************************************/
182
183 Double_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
199 Double_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
212 Double_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
286 Double_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
301 Double_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
314 Double_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
363 Double_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
377 Double_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
391 Double_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
404 Double_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
418 Double_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 }