Balance function analysis (P.Christakoglou)
[u/mrichter/AliRoot.git] / ANALYSIS / AliBalance.cxx
1 /**************************************************************************
2  * Author: Panos Christakoglou.                                           *
3  * Contributors are mentioned in the code where appropriate.              *
4  *                                                                        *
5  * Permission to use, copy, modify and distribute this software and its   *
6  * documentation strictly for non-commercial purposes is hereby granted   *
7  * without fee, provided that the above copyright notice appears in all   *
8  * copies and that both the copyright notice and this permission notice   *
9  * appear in the supporting documentation. The authors make no claims     *
10  * about the suitability of this software for any purpose. It is          *
11  * provided "as is" without express or implied warranty.                  *
12  **************************************************************************/
13
14 /* $Id$ */
15
16 //-----------------------------------------------------------------
17 //           Balance Function class
18 //   This is the class to deal with the Balance Function analysis
19 //   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
20 //-----------------------------------------------------------------
21
22
23 #include <stdlib.h>
24
25 //ROOT
26 #include <Riostream.h>
27 #include <TROOT.h>
28 #include <TObject.h>
29 #include <TSystem.h>
30 #include <TObject.h>
31 #include <TMath.h>
32 #include <TLorentzVector.h>
33
34 #include "AliBalance.h"
35
36 ClassImp(AliBalance)
37
38 //----------------------------------------//
39 AliBalance::AliBalance()
40 {
41   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++)
42     {
43       fNpp[i] = .0;
44       fNnn[i] = .0;
45       fNpn[i] = .0;
46       fB[i] = 0.0;
47       ferror[i] = 0.0;
48     } 
49   fAnalyzedEvents = 0;
50 }
51
52 //----------------------------------------//
53 AliBalance::AliBalance(Double_t P2_Start, Double_t P2_Stop, Int_t P2_bins)
54 {
55   this->fP2_Start = P2_Start;
56   this->fP2_Stop = P2_Stop;
57   this->fNumberOfBins = P2_bins;
58   this->fP2_Step = fabs(fP2_Start - fP2_Stop) / (Double_t)fNumberOfBins;
59 }
60
61 //----------------------------------------//
62 AliBalance::~AliBalance()
63 {
64 }
65
66 //----------------------------------------//
67 void AliBalance::SetNumberOfBins(Int_t ibins)
68 {
69   this->fNumberOfBins = ibins ;
70 }
71
72 //----------------------------------------//
73 void AliBalance::SetInterval(Double_t P2_Start, Double_t P2_Stop)
74 {
75   this->fP2_Start = P2_Start;
76   this->fP2_Stop = P2_Stop;
77   this->fP2_Step = fabs(P2_Start - P2_Stop) / (Double_t)fNumberOfBins;
78 }
79
80 //----------------------------------------//
81 void AliBalance::SetAnalysisType(Int_t iType)
82 {
83   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
84   this->fAnalysisType = iType; 
85   if(fAnalysisType==0)
86     {
87       cout<<" ====================== "<<endl;
88       cout<<"||Analysis selected: y||"<<endl;
89       cout<<" ====================== "<<endl;
90     } 
91   else if(fAnalysisType==1)
92     {
93       cout<<" ======================== "<<endl;
94       cout<<"||Analysis selected: eta||"<<endl;
95       cout<<" ======================== "<<endl;
96     }
97   else if(fAnalysisType==2)
98     {
99       cout<<" ========================== "<<endl;
100       cout<<"||Analysis selected: Qlong||"<<endl;
101       cout<<" ========================== "<<endl;
102     }
103   else if(fAnalysisType==3)
104     {
105       cout<<" ========================= "<<endl;
106       cout<<"||Analysis selected: Qout||"<<endl;
107       cout<<" ========================= "<<endl;
108     }
109   else if(fAnalysisType==4)
110     {
111       cout<<" ========================== "<<endl;
112       cout<<"||Analysis selected: Qside||"<<endl;
113       cout<<" ========================== "<<endl;
114     }
115   else if(fAnalysisType==5)
116     {
117       cout<<" ========================= "<<endl;
118       cout<<"||Analysis selected: Qinv||"<<endl;
119       cout<<" ========================= "<<endl;
120     }
121   else if(fAnalysisType==6)
122     {
123       cout<<" ======================== "<<endl;
124       cout<<"||Analysis selected: phi||"<<endl;
125       cout<<" ======================== "<<endl;
126     }
127   else
128     {
129       cout<<"Selection of analysis mode failed!!!"<<endl;
130       cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
131       abort();
132     }
133 }
134
135 //----------------------------------------//
136 void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim)
137 {
138   this->fV = P;
139   this->fCharge = charge;
140   fNtrack = dim;
141 }
142
143
144 //----------------------------------------//
145 void AliBalance::CalculateBalance()
146 {
147   fAnalyzedEvents++;
148   Int_t i = 0 , j = 0;
149   Int_t ibin = 0;
150
151   for(i = 0; i < fNtrack; i++)
152     {
153       if(fCharge[i] > 0)
154         fNp += 1.;
155       if(fCharge[i] < 0)
156         fNn += 1.;
157     }
158
159   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
160    if(fAnalysisType==0)
161     {
162       for(i = 1; i < fNtrack; i++)
163         {
164           for(j = 0; j < i; j++)
165             {
166               Double_t rap1 = 0.5*log((fV[i].E() - fV[i].Pz())/(fV[i].E() + fV[i].Pz())); 
167               Double_t rap2 = 0.5*log((fV[j].E() - fV[j].Pz())/(fV[j].E() + fV[j].Pz())); 
168               Double_t dy = TMath::Abs(rap1 - rap2);
169               ibin = Int_t(dy/fP2_Step);
170               if((fCharge[i] > 0)&&(fCharge[j] > 0))
171                 fNpp[ibin] += 1.;
172               if((fCharge[i] < 0)&&(fCharge[j] < 0))
173                 fNnn[ibin] += 1.;
174               if((fCharge[i] > 0)&&(fCharge[j] < 0))
175                 fNpn[ibin] += 1.;
176               if((fCharge[i] < 0)&&(fCharge[j] > 0))
177                 fNpn[ibin] += 1.;
178             }
179         }
180     }//case 0
181    if(fAnalysisType==1)
182     {
183       for(i = 1; i < fNtrack; i++)
184         {
185           for(j = 0; j < i; j++)
186             {
187               Double_t P1 = sqrt(pow(fV[i].Px(),2) + pow(fV[i].Py(),2) + pow(fV[i].Pz(),2)); 
188               Double_t P2 = sqrt(pow(fV[j].Px(),2) + pow(fV[j].Py(),2) + pow(fV[j].Pz(),2));
189               Double_t eta1 = 0.5*log((P1 - fV[i].Pz())/(P1 + fV[i].Pz())); 
190               Double_t eta2 = 0.5*log((P2 - fV[j].Pz())/(P2 + fV[j].Pz())); 
191               Double_t deta = TMath::Abs(eta1 - eta2);
192               ibin = Int_t(deta/fP2_Step);
193               if((fCharge[i] > 0)&&(fCharge[j] > 0))
194                 fNpp[ibin] += 1.;
195               if((fCharge[i] < 0)&&(fCharge[j] < 0))
196                 fNnn[ibin] += 1.;
197               if((fCharge[i] > 0)&&(fCharge[j] < 0))
198                 fNpn[ibin] += 1.;
199               if((fCharge[i] < 0)&&(fCharge[j] > 0))
200                 fNpn[ibin] += 1.;
201             }
202         }
203     }//case 1
204    if(fAnalysisType==2)
205     {
206       for(i = 1; i < fNtrack; i++)
207         {
208           for(j = 0; j < i; j++)
209             {
210               Double_t ETot = fV[i].E() + fV[j].E();
211               Double_t PxTot = fV[i].Px() + fV[j].Px();
212               Double_t PyTot = fV[i].Py() + fV[j].Py();
213               Double_t PzTot = fV[i].Pz() + fV[j].Pz();
214               Double_t Q0Tot = fV[i].E() - fV[j].E();
215               Double_t QzTot = fV[i].Pz() - fV[j].Pz();
216               Double_t Snn = pow(ETot,2) - pow(PxTot,2) - pow(PyTot,2) - pow(PzTot,2);
217               Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
218               Double_t Qlong = TMath::Abs(ETot*QzTot - PzTot*Q0Tot)/sqrt(Snn + pow(PtTot,2));
219               ibin = Int_t(Qlong/fP2_Step);
220               //cout<<ibin<<endl;
221               if((fCharge[i] > 0)&&(fCharge[j] > 0))
222                 fNpp[ibin] += 1.;
223               if((fCharge[i] < 0)&&(fCharge[j] < 0))
224                 fNnn[ibin] += 1.;
225               if((fCharge[i] > 0)&&(fCharge[j] < 0))
226                 fNpn[ibin] += 1.;
227               if((fCharge[i] < 0)&&(fCharge[j] > 0))
228                 fNpn[ibin] += 1.;
229             }
230         }
231     }//case 2
232   if(fAnalysisType==3)
233     {
234       for(i = 1; i < fNtrack; i++)
235         {
236           for(j = 0; j < i; j++)
237             {
238               Double_t ETot = fV[i].E() + fV[j].E();
239               Double_t PxTot = fV[i].Px() + fV[j].Px();
240               Double_t PyTot = fV[i].Py() + fV[j].Py();
241               Double_t PzTot = fV[i].Pz() + fV[j].Pz();
242               Double_t QxTot = fV[i].Px() - fV[j].Px();
243               Double_t QyTot = fV[i].Py() - fV[j].Py();
244               Double_t Snn = pow(ETot,2) - pow(PxTot,2) - pow(PyTot,2) - pow(PzTot,2);
245               Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
246               Double_t Qout = sqrt(Snn/(Snn + pow(PtTot,2))) * TMath::Abs(PxTot*QxTot + PyTot*QyTot)/PtTot;
247               ibin = Int_t(Qout/fP2_Step);
248               //cout<<ibin<<endl;
249               if((fCharge[i] > 0)&&(fCharge[j] > 0))
250                 fNpp[ibin] += 1.;
251               if((fCharge[i] < 0)&&(fCharge[j] < 0))
252                 fNnn[ibin] += 1.;
253               if((fCharge[i] > 0)&&(fCharge[j] < 0))
254                 fNpn[ibin] += 1.;
255               if((fCharge[i] < 0)&&(fCharge[j] > 0))
256                 fNpn[ibin] += 1.;
257             }
258         }
259     }//case 3
260   if(fAnalysisType==4)
261     {
262       for(i = 1; i < fNtrack; i++)
263         {
264           for(j = 0; j < i; j++)
265             {
266               Double_t PxTot = fV[i].Px() + fV[j].Px();
267               Double_t PyTot = fV[i].Py() + fV[j].Py();
268               Double_t QxTot = fV[i].Px() - fV[j].Px();
269               Double_t QyTot = fV[i].Py() - fV[j].Py();
270               Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
271               Double_t Qside = TMath::Abs(PxTot*QyTot - PyTot*QxTot)/PtTot;
272               ibin = Int_t(Qside/fP2_Step);
273               //cout<<ibin<<endl;
274               if((fCharge[i] > 0)&&(fCharge[j] > 0))
275                 fNpp[ibin] += 1.;
276               if((fCharge[i] < 0)&&(fCharge[j] < 0))
277                 fNnn[ibin] += 1.;
278               if((fCharge[i] > 0)&&(fCharge[j] < 0))
279                 fNpn[ibin] += 1.;
280               if((fCharge[i] < 0)&&(fCharge[j] > 0))
281                 fNpn[ibin] += 1.;
282             }
283         }
284     }//case 4
285    if(fAnalysisType==5)
286     {
287       for(i = 1; i < fNtrack; i++)
288         {
289           for(j = 0; j < i; j++)
290             {
291               Double_t Q0Tot = fV[i].E() - fV[j].E();
292               Double_t QxTot = fV[i].Px() - fV[j].Px();
293               Double_t QyTot = fV[i].Py() - fV[j].Py();
294               Double_t QzTot = fV[i].Pz() - fV[j].Pz();
295               Double_t Qinv = sqrt(TMath::Abs(-pow(Q0Tot,2) +pow(QxTot,2) +pow(QyTot,2) +pow(QzTot,2)));
296               ibin = Int_t(Qinv/fP2_Step);
297               //cout<<ibin<<endl;
298               if((fCharge[i] > 0)&&(fCharge[j] > 0))
299                 fNpp[ibin] += 1.;
300               if((fCharge[i] < 0)&&(fCharge[j] < 0))
301                 fNnn[ibin] += 1.;
302               if((fCharge[i] > 0)&&(fCharge[j] < 0))
303                 fNpn[ibin] += 1.;
304               if((fCharge[i] < 0)&&(fCharge[j] > 0))
305                 fNpn[ibin] += 1.;
306             }
307         }
308     }//case 5   
309   if(fAnalysisType==6)
310     {
311       for(i = 1; i < fNtrack; i++)
312         {
313           for(j = 0; j < i; j++)
314             {
315               Double_t phi1 = TMath::ATan(fV[i].Py()/fV[i].Px())*180.0/TMath::Pi();
316               Double_t phi2 = TMath::ATan(fV[j].Py()/fV[j].Px())*180.0/TMath::Pi();
317               Double_t dphi = TMath::Abs(phi1 - phi2);
318               ibin = Int_t(dphi/fP2_Step);
319               if((fCharge[i] > 0)&&(fCharge[j] > 0))
320                 fNpp[ibin] += 1.;
321               if((fCharge[i] < 0)&&(fCharge[j] < 0))
322                 fNnn[ibin] += 1.;
323               if((fCharge[i] > 0)&&(fCharge[j] < 0))
324                 fNpn[ibin] += 1.;
325               if((fCharge[i] < 0)&&(fCharge[j] > 0))
326                 fNpn[ibin] += 1.;
327             }
328         }
329     }//case 6
330 }
331
332 //----------------------------------------//
333 Double_t AliBalance::GetBalance(Int_t p2)
334 {
335   fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2_Step;
336
337   return fB[p2];
338 }
339     
340 //----------------------------------------//
341 Double_t AliBalance::GetError(Int_t p2)
342 {
343   ferror[p2] = sqrt( Double_t(fNpp[p2])/(Double_t(fNp)*Double_t(fNp)) + Double_t(fNnn[p2])/(Double_t(fNn)*Double_t(fNn)) + Double_t(fNpn[p2])*pow((0.5/Double_t(fNp) + 0.5/Double_t(fNn)),2))/fP2_Step;
344
345   return ferror[p2];
346 }
347
348 //----------------------------------------//
349 void AliBalance::PrintResults()
350 {
351   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
352   Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
353   Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
354   Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
355   Double_t delta_bal_P2 = 0.0, Integral = 0.0;
356   Double_t DeltaErrorNew = 0.0;
357
358   cout<<"=================================================="<<endl;
359   for(Int_t i = 0; i < fNumberOfBins; i++)
360     { 
361       Double_t x = fP2_Start + fP2_Step*i + fP2_Step/2 ;
362       cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x<<endl;
363     } 
364   cout<<"=================================================="<<endl;
365   for(Int_t i = 1; i < fNumberOfBins; i++)
366     {
367       fSumXi += x[i];
368       fSumBi += fB[i];
369       fSumBiXi += fB[i]*x[i];
370       fSumBiXi2 += fB[i]*pow(x[i],2);
371       fSumBi2Xi2 += pow(fB[i],2)*pow(x[i],2);
372       fSumDeltaBi2 +=  pow(ferror[i],2) ;
373       fSumXi2DeltaBi2 += pow(x[i],2) * pow(ferror[i],2) ;
374       
375       delta_bal_P2 += fP2_Step*pow(ferror[i],2) ;
376       Integral += fP2_Step*fB[i] ;
377     }
378   for(Int_t i = 1; i < fNumberOfBins; i++)
379     {
380       DeltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/pow(fSumBi,2);
381     }
382   Double_t IntegralError = sqrt(delta_bal_P2) ;
383   
384   Double_t Delta = fSumBiXi / fSumBi ;
385   Double_t DeltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) ) ;
386  
387   cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
388   cout<<"Width: "<<Delta<<"\t Error: "<<DeltaError<<endl;
389   cout<<"New error: "<<DeltaErrorNew<<endl;
390   cout<<"Interval: "<<Integral<<"\t Error: "<<IntegralError<<endl;
391   cout<<"=================================================="<<endl;
392 }
393