a385fb6edee094a26caefd3d8112b4985f1deac9
[u/mrichter/AliRoot.git] / HBTAN / AliHBTPair.cxx
1 #include "AliHBTPair.h"
2 #include "AliHBTParticle.h"
3 #include "AliHBTLLWeights.h"
4
5 ClassImp(AliHBTPair)
6
7 /************************************************************************/
8 AliHBTPair::AliHBTPair(Bool_t rev):
9  fPart1(0x0),
10  fPart2(0x0),
11  fSwapedPair(0x0),
12  fQSideCMSLC(0.0),
13  fQSideCMSLCNotCalc(kTRUE),
14  fQOutCMSLC(0.0),
15  fQOutCMSLCNotCalc(kTRUE),
16  fQLongCMSLC(0.0),
17  fQLongCMSLCNotCalc(kTRUE),
18  fQInv(0.0),
19  fQInvNotCalc(kTRUE),
20  fInvMass(0.0),
21  fInvMassNotCalc(kTRUE),
22  fKt(0.0),
23  fKtNotCalc(kTRUE),
24  fKStar(0.0),
25  fKStarNotCalc(kTRUE),
26  fPInv(0.0),
27  fQSide(0.0),
28  fOut(0.0),
29  fQLong(0.0),
30  fMt(0.0),
31  fMtNotCalc(kTRUE),
32  fInvMassSqr(0.0),
33  fMassSqrNotCalc(kTRUE),
34  fQInvL(0.0),
35  fQInvLNotCalc(kTRUE),
36  fLLWeight(0.0),
37  ffLLWeightNotCalc(kTRUE),
38  fPxSum(0.0),
39  fPySum(0.0),
40  fPzSum(0.0),
41  fESum(0.0),
42  fSumsNotCalc(kTRUE),
43  fPxDiff(0.0),
44  fPyDiff(0.0),
45  fPzDiff(0.0),
46  fEDiff(0.0),
47  fDiffsNotCalc(kTRUE),
48  fGammaCMSLC(0.0),
49  fGammaCMSLCNotCalc(kTRUE),
50  fChanged(kTRUE)
51  {
52 //value of rev defines if it is Swaped
53 //if you pass kTRUE swpaped pair will NOT be created
54 //though you wont be able to get the swaped pair from this pair
55
56   if(!rev) fSwapedPair = new AliHBTPair(kTRUE); //if false create swaped pair
57   
58  }
59 /************************************************************************/
60
61 AliHBTPair::AliHBTPair(AliHBTParticle* part1, AliHBTParticle* part2, Bool_t rev):
62  fPart1(part1),
63  fPart2(part2),
64  fSwapedPair(0x0),
65  fQSideCMSLC(0.0),
66  fQSideCMSLCNotCalc(kTRUE),
67  fQOutCMSLC(0.0),
68  fQOutCMSLCNotCalc(kTRUE),
69  fQLongCMSLC(0.0),
70  fQLongCMSLCNotCalc(kTRUE),
71  fQInv(0.0),
72  fQInvNotCalc(kTRUE),
73  fInvMass(0.0),
74  fInvMassNotCalc(kTRUE),
75  fKt(0.0),
76  fKtNotCalc(kTRUE),
77  fKStar(0.0),
78  fKStarNotCalc(kTRUE),
79  fPInv(0.0),
80  fQSide(0.0),
81  fOut(0.0),
82  fQLong(0.0),
83  fMt(0.0),
84  fMtNotCalc(kTRUE),
85  fInvMassSqr(0.0),
86  fMassSqrNotCalc(kTRUE),
87  fQInvL(0.0),
88  fQInvLNotCalc(kTRUE),
89  fLLWeight(0.0),
90  ffLLWeightNotCalc(kTRUE),
91  fPxSum(0.0),
92  fPySum(0.0),
93  fPzSum(0.0),
94  fESum(0.0),
95  fSumsNotCalc(kTRUE),
96  fPxDiff(0.0),
97  fPyDiff(0.0),
98  fPzDiff(0.0),
99  fEDiff(0.0),
100  fDiffsNotCalc(kTRUE),
101  fGammaCMSLC(0.0),
102  fGammaCMSLCNotCalc(kTRUE),
103  fChanged(kTRUE)
104  {
105 //value of rev defines if it is Swaped
106 //if you pass kTRUE swpaped pair will NOT be created
107 //though you wont be able to get the swaped pair from this pair
108
109   if(!rev) fSwapedPair = new AliHBTPair(part2,part1,kTRUE); //if false create swaped pair
110   
111  }
112 /************************************************************************/
113
114 Double_t AliHBTPair::GetInvMass()
115 {
116 //Returns qinv value for a pair
117   if(fInvMassNotCalc)
118    {
119      CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method 
120      
121      if(fInvMassSqr<0)  fInvMass = TMath::Sqrt(-fInvMassSqr);
122      else fInvMass = TMath::Sqrt(fInvMassSqr); 
123      
124      fInvMassNotCalc = kFALSE;
125    }
126   return fInvMass;
127 }
128 /************************************************************************/
129 Double_t AliHBTPair::GetQSideCMSLC()
130 {
131   //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
132  
133   if (fQSideCMSLCNotCalc)
134    {
135     fQSideCMSLC = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
136     fQSideCMSLCNotCalc = kFALSE;
137    }
138   return fQSideCMSLC;
139 }
140 /************************************************************************/
141 Double_t AliHBTPair::GetQOutCMSLC()
142 {
143  if(fQOutCMSLCNotCalc)
144   {
145    CalculateSums();
146    CalculateDiffs();
147
148    if (fPart1->GetMass() != fPart2->GetMass())
149     {
150 /*    
151       //STAR algorithm
152       Double_t beta  = fPzSum/fESum;
153       Double_t gamma = GetGammaToCMSLC();
154       Double_t el = gamma * (fPart1->Energy() - beta * fPart1->Pz());
155       Double_t x  = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
156       beta  = 2.0*GetKt()/GetMt();
157       gamma = GetMt()/GetQInv();
158       fQOutCMSLC = gamma * (x - beta * el);
159 */
160
161       //beta=fPzSum/fESum;    // Longit. V == beta
162       Double_t beta=fPzSum/fESum;
163       Double_t gamma = GetGammaToCMSLC();
164       
165       Double_t cosphi=fPxSum/(2.0*GetKt());  // cos(phi)
166       Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
167       
168 //      ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
169 //      ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
170       Double_t tmp;
171       tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
172       Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
173       Double_t part1Px = tmp;
174
175       tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
176       Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
177       Double_t part2Px = tmp;
178       
179       
180 //      LTR(part1Pz,E1,beta,GetGammaToCMSLC(),part1Pz,E1a);
181 //      LTR(part2Pz,E2,beta,GetGammaToCMSLC(),part2Pz,E2a);
182       Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->Energy());
183       Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->Energy());
184
185       Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
186       Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
187       Double_t part1E=TMath::Sqrt(fPart1->GetMass()*fPart1->GetMass()+part1P2);
188       Double_t part2E=TMath::Sqrt(fPart2->GetMass()*fPart2->GetMass()+part2P2);
189       Double_t sumE=part1E+part2E;
190       Double_t sumPx=part1Px+part2Px;
191       Double_t sumPy=part1Py+part2Py;
192       Double_t sumPZ=part1Pz+part2Pz;
193       Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
194
195       Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
196       Double_t hf = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/(relmass*relmass);
197       fQOutCMSLC=(part1Px-part2Px);//== id
198       fQOutCMSLC=fQOutCMSLC-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
199     }
200    else
201     {
202       Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
203       fQOutCMSLC = 0.5*k2/GetKt();
204    // if (non-id)  fQOutCMSLC=fQOutCMSLC - sumPx*HF;
205     }
206
207     
208    fQOutCMSLCNotCalc = kFALSE;
209   }
210  return fQOutCMSLC;
211 }
212 /************************************************************************/
213 Double_t AliHBTPair::GetQLongCMSLC()
214 {
215  if (fQLongCMSLCNotCalc)
216   {
217     CalculateSums();
218     CalculateDiffs();
219     Double_t beta = fPzSum/fESum;
220     fQLongCMSLC = GetGammaToCMSLC() * ( fPzDiff - beta*fEDiff );
221     fQLongCMSLCNotCalc = kFALSE;
222   }
223  return fQLongCMSLC; 
224 }
225 /************************************************************************/
226 Double_t AliHBTPair::GetKt()
227 {
228   if(fKtNotCalc)
229    { 
230      CalculateSums();
231      fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
232      fKtNotCalc = kFALSE;
233    }
234   return fKt;
235 }
236 /************************************************************************/
237
238 Double_t AliHBTPair::GetKStar()
239 {
240   if (fKStarNotCalc)
241    { 
242     CalculateSums();
243
244     Double_t Ptrans = fPxSum*fPxSum + fPySum*fPySum;
245     Double_t Mtrans = fESum*fESum - fPzSum*fPzSum;
246     Double_t Pinv =   TMath::Sqrt(Mtrans - Ptrans);
247
248     Double_t Q = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/Pinv;
249     
250     CalculateQInvL();
251     
252     Q = Q*Q - fQInvL;
253     if ( Q < 0)
254      {
255         Info("GetKStar","Q = %f",Q);
256         fPart1->Print();
257         fPart2->Print();
258         Q = TMath::Abs(Q);
259      }
260      
261     Q = TMath::Sqrt(Q);
262     fKStar = Q/2.;
263     fKStarNotCalc = kFALSE;
264    }
265   return fKStar;
266 }
267 /************************************************************************/
268
269 Double_t AliHBTPair::GetQInv()
270 {
271 //returns Qinv 
272 //warning for non-id particles you want to use 2*KStar
273   if(fQInvNotCalc)
274    {
275     CalculateQInvL();
276     fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
277     fQInvNotCalc = kFALSE;
278    }
279   return fQInv;
280 }
281 /************************************************************************/
282
283 Double_t AliHBTPair::GetGammaToCMSLC()
284 {
285   if(fGammaCMSLCNotCalc)
286    {
287      CalculateSums();
288      Double_t beta = fPzSum/fESum;
289      fGammaCMSLC = 1.0/TMath::Sqrt(1.0 - beta*beta);
290      fGammaCMSLCNotCalc = kFALSE;
291    }
292   return fGammaCMSLC;
293 }
294 /************************************************************************/
295
296 Double_t AliHBTPair::GetMt()
297 {
298   if (fMtNotCalc)
299    {
300      CalculateSums();
301      fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
302      fMtNotCalc = kFALSE;
303    }
304    return fMt;
305 }
306 /************************************************************************/
307
308 Double_t AliHBTPair::GetLLWeight()
309 {
310   if (ffLLWeightNotCalc)
311    {
312       fLLWeight = AliHBTLLWeights::Instance()->GetWeight(this);
313       ffLLWeightNotCalc = kFALSE;
314    }
315   return fLLWeight; 
316 }