Catching up to NewIO -> Particle stores all passible PID and their probabilities
[u/mrichter/AliRoot.git] / HBTAN / AliHBTPair.cxx
1 #include "AliHBTPair.h"
2 #include "AliHBTParticle.h"
3
4 ClassImp(AliHBTPair)
5
6 /************************************************************************/
7 AliHBTPair::AliHBTPair(Bool_t rev):
8  fPart1(0x0),
9  fPart2(0x0)
10  {
11 //value of rev defines if it is Swaped
12 //if you pass kTRUE swpaped pair will NOT be created
13 //though you wont be able to get the swaped pair from this pair
14
15   if(!rev) fSwapedPair = new AliHBTPair(kTRUE); //if false create swaped pair
16   else fSwapedPair = 0x0; //if true set the pointer to NULL
17   
18  }
19 /************************************************************************/
20
21 AliHBTPair::AliHBTPair(AliHBTParticle* part1, AliHBTParticle* part2, Bool_t rev):
22  fPart1(part1),
23  fPart2(part2)
24  {
25 //value of rev defines if it is Swaped
26 //if you pass kTRUE swpaped pair will NOT be created
27 //though you wont be able to get the swaped pair from this pair
28
29   if(!rev) fSwapedPair = new AliHBTPair(part2,part1,kTRUE); //if false create swaped pair
30   else fSwapedPair = 0x0; //if true set the pointer to NULL
31   
32  }
33 /************************************************************************/
34
35 Double_t AliHBTPair::GetInvMass()
36 {
37 //Returns qinv value for a pair
38   if(fInvMassNotCalc)
39    {
40      CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method 
41      
42      if(fInvMassSqr<0)  fInvMass = TMath::Sqrt(-fInvMassSqr);
43      else fInvMass = TMath::Sqrt(fInvMassSqr); 
44      
45      fInvMassNotCalc = kFALSE;
46    }
47   return fInvMass;
48 }
49 /************************************************************************/
50 Double_t AliHBTPair::GetQSideCMSLC()
51 {
52   //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
53  
54   if (fQSideCMSLCNotCalc)
55    {
56     fQSideCMSLC = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
57     fQSideCMSLCNotCalc = kFALSE;
58    }
59   return fQSideCMSLC;
60 }
61 /************************************************************************/
62 Double_t AliHBTPair::GetQOutCMSLC()
63 {
64  if(fQOutCMSLCNotCalc)
65   {
66    CalculateSums();
67    CalculateDiffs();
68
69    if (fPart1->GetMass() != fPart2->GetMass())
70     {
71 /*    
72       //STAR algorithm
73       Double_t beta  = fPzSum/fESum;
74       Double_t gamma = GetGammaToCMSLC();
75       Double_t el = gamma * (fPart1->Energy() - beta * fPart1->Pz());
76       Double_t x  = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
77       beta  = 2.0*GetKt()/GetMt();
78       gamma = GetMt()/GetQInv();
79       fQOutCMSLC = gamma * (x - beta * el);
80 */
81
82       //beta=fPzSum/fESum;    // Longit. V == beta
83       Double_t beta=fPzSum/fESum;
84       Double_t gamma = GetGammaToCMSLC();
85       
86       Double_t cosphi=fPxSum/(2.0*GetKt());  // cos(phi)
87       Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
88       
89 //      ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
90 //      ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
91       Double_t tmp;
92       tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
93       Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
94       Double_t part1Px = tmp;
95
96       tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
97       Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
98       Double_t part2Px = tmp;
99       
100       
101 //      LTR(part1Pz,E1,beta,GetGammaToCMSLC(),part1Pz,E1a);
102 //      LTR(part2Pz,E2,beta,GetGammaToCMSLC(),part2Pz,E2a);
103       Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->Energy());
104       Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->Energy());
105
106       Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
107       Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
108       Double_t part1E=TMath::Sqrt(fPart1->GetMass()*fPart1->GetMass()+part1P2);
109       Double_t part2E=TMath::Sqrt(fPart2->GetMass()*fPart2->GetMass()+part2P2);
110       Double_t sumE=part1E+part2E;
111       Double_t sumPx=part1Px+part2Px;
112       Double_t sumPy=part1Py+part2Py;
113       Double_t sumPZ=part1Pz+part2Pz;
114       Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
115
116       Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
117       Double_t hf = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/(relmass*relmass);
118       fQOutCMSLC=(part1Px-part2Px);//== id
119       fQOutCMSLC=fQOutCMSLC-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
120     }
121    else
122     {
123       Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
124       fQOutCMSLC = 0.5*k2/GetKt();
125    // if (non-id)  fQOutCMSLC=fQOutCMSLC - sumPx*HF;
126     }
127
128     
129    fQOutCMSLCNotCalc = kFALSE;
130   }
131  return fQOutCMSLC;
132 }
133 /************************************************************************/
134 Double_t AliHBTPair::GetQLongCMSLC()
135 {
136  if (fQLongCMSLCNotCalc)
137   {
138     CalculateSums();
139     CalculateDiffs();
140     Double_t beta = fPzSum/fESum;
141     fQLongCMSLC = GetGammaToCMSLC() * ( fPzDiff - beta*fEDiff );
142     fQLongCMSLCNotCalc = kFALSE;
143   }
144  return fQLongCMSLC; 
145 }
146 /************************************************************************/
147 Double_t AliHBTPair::GetKt()
148 {
149   if(fKtNotCalc)
150    { 
151      CalculateSums();
152      fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
153      fKtNotCalc = kFALSE;
154    }
155   return fKt;
156 }
157 /************************************************************************/
158
159 Double_t AliHBTPair::GetKStar()
160 {
161   if (fKStarNotCalc)
162    { 
163     CalculateSums();
164
165     Double_t Ptrans = fPxSum*fPxSum + fPySum*fPySum;
166     Double_t Mtrans = fESum*fESum - fPzSum*fPzSum;
167     Double_t Pinv =   TMath::Sqrt(Mtrans - Ptrans);
168
169     Double_t Q = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/Pinv;
170     
171     CalculateQInvL();
172     
173     Q = TMath::Sqrt( Q*Q - fQInvL);
174     fKStar = Q/2.;
175     fKStarNotCalc = kFALSE;
176    }
177   return fKStar;
178 }
179 /************************************************************************/
180
181 Double_t AliHBTPair::GetQInv()
182 {
183 //returns Qinv 
184 //warning for non-id particles you want to use 2*KStar
185   if(fQInvNotCalc)
186    {
187     CalculateQInvL();
188     fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
189     fQInvNotCalc = kFALSE;
190    }
191   return fQInv;
192 }
193 /************************************************************************/
194
195 Double_t AliHBTPair::GetGammaToCMSLC()
196 {
197   if(fGammaCMSLCNotCalc)
198    {
199      CalculateSums();
200      Double_t beta = fPzSum/fESum;
201      fGammaCMSLC = 1.0/TMath::Sqrt(1.0 - beta*beta);
202      fGammaCMSLCNotCalc = kFALSE;
203    }
204   return fGammaCMSLC;
205 }
206 /************************************************************************/
207
208 Double_t AliHBTPair::GetMt()
209 {
210   if (fMtNotCalc)
211    {
212      CalculateSums();
213      fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
214      fMtNotCalc = kFALSE;
215    }
216    return fMt;
217 }
218 /************************************************************************/
219
220
221
222
223
224
225