Coding convention
[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 //ROOT
24 #include <Riostream.h>
25 #include <TMath.h>
26 #include <TLorentzVector.h>
27
28 #include "AliBalance.h"
29
30 ClassImp(AliBalance)
31
32 //----------------------------------------//
33 AliBalance::AliBalance()
34 {
35   // Default constructor
36   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++)
37     {
38       fNpp[i] = .0;
39       fNnn[i] = .0;
40       fNpn[i] = .0;
41       fB[i] = 0.0;
42       ferror[i] = 0.0;
43     } 
44   fNp = 0.0;
45   fNn = 0.0;
46   fP2Start = 0.0;
47   fP2Stop = 0.0;
48   fP2Step = 0.0; 
49   fAnalysisType = 0;
50   fNumberOfBins = 0;
51   fNtrack = 0;
52   fAnalyzedEvents = 0;
53 }
54
55 //----------------------------------------//
56 AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins)
57 {
58   // Constructor
59   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++)
60     {
61       fNpp[i] = .0;
62       fNnn[i] = .0;
63       fNpn[i] = .0;
64       fB[i] = 0.0;
65       ferror[i] = 0.0;
66     } 
67   fNp = 0.0;
68   fNn = 0.0;
69   fAnalysisType = 0;
70   fNtrack = 0;
71   fAnalyzedEvents = 0;
72
73   fP2Start = p2Start;
74   fP2Stop = p2Stop;
75   fNumberOfBins = p2Bins;
76   fP2Step = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins;
77 }
78
79 //----------------------------------------//
80 AliBalance::~AliBalance()
81 {
82   // Destructor
83 }
84
85 //----------------------------------------//
86 void AliBalance::SetNumberOfBins(Int_t ibins)
87 {
88   // Sets the number of bins for the analyzed interval
89   fNumberOfBins = ibins ;
90 }
91
92 //----------------------------------------//
93 void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop)
94 {
95   // Sets the analyzed interval. 
96   // The analysis variable is set by SetAnalysisType
97   fP2Start = p2Start;
98   fP2Stop = p2Stop;
99   fP2Step = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins;
100 }
101
102 //----------------------------------------//
103 void AliBalance::SetAnalysisType(Int_t iType)
104 {
105   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
106   this->fAnalysisType = iType; 
107   if(fAnalysisType==0)
108     {
109       cout<<" ====================== "<<endl;
110       cout<<"||Analysis selected: y||"<<endl;
111       cout<<" ====================== "<<endl;
112     } 
113   else if(fAnalysisType==1)
114     {
115       cout<<" ======================== "<<endl;
116       cout<<"||Analysis selected: eta||"<<endl;
117       cout<<" ======================== "<<endl;
118     }
119   else if(fAnalysisType==2)
120     {
121       cout<<" ========================== "<<endl;
122       cout<<"||Analysis selected: Qlong||"<<endl;
123       cout<<" ========================== "<<endl;
124     }
125   else if(fAnalysisType==3)
126     {
127       cout<<" ========================= "<<endl;
128       cout<<"||Analysis selected: Qout||"<<endl;
129       cout<<" ========================= "<<endl;
130     }
131   else if(fAnalysisType==4)
132     {
133       cout<<" ========================== "<<endl;
134       cout<<"||Analysis selected: Qside||"<<endl;
135       cout<<" ========================== "<<endl;
136     }
137   else if(fAnalysisType==5)
138     {
139       cout<<" ========================= "<<endl;
140       cout<<"||Analysis selected: Qinv||"<<endl;
141       cout<<" ========================= "<<endl;
142     }
143   else if(fAnalysisType==6)
144     {
145       cout<<" ======================== "<<endl;
146       cout<<"||Analysis selected: phi||"<<endl;
147       cout<<" ======================== "<<endl;
148     }
149   else
150     {
151       cout<<"Selection of analysis mode failed!!!"<<endl;
152       cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
153       abort();
154     }
155 }
156
157 //----------------------------------------//
158 void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim)
159 {
160   // Sets a new particle with given 4-momentum and charge.
161   // dim is the size of the array of charges and corresponds
162   // to the number of selected tracks.
163   this->fV = P;
164   this->fCharge = charge;
165   fNtrack = dim;
166 }
167
168
169 //----------------------------------------//
170 void AliBalance::CalculateBalance()
171 {
172   // Calculates the balance function
173   fAnalyzedEvents++;
174   Int_t i = 0 , j = 0;
175   Int_t ibin = 0;
176
177   for(i = 0; i < fNtrack; i++)
178     {
179       if(fCharge[i] > 0)
180         fNp += 1.;
181       if(fCharge[i] < 0)
182         fNn += 1.;
183     }
184
185   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
186    if(fAnalysisType==0)
187     {
188       for(i = 1; i < fNtrack; i++)
189         {
190           for(j = 0; j < i; j++)
191             {
192               Double_t rap1 = 0.5*log((fV[i].E() - fV[i].Pz())/(fV[i].E() + fV[i].Pz())); 
193               Double_t rap2 = 0.5*log((fV[j].E() - fV[j].Pz())/(fV[j].E() + fV[j].Pz())); 
194               Double_t dy = TMath::Abs(rap1 - rap2);
195               ibin = Int_t(dy/fP2Step);
196               if((fCharge[i] > 0)&&(fCharge[j] > 0))
197                 fNpp[ibin] += 1.;
198               if((fCharge[i] < 0)&&(fCharge[j] < 0))
199                 fNnn[ibin] += 1.;
200               if((fCharge[i] > 0)&&(fCharge[j] < 0))
201                 fNpn[ibin] += 1.;
202               if((fCharge[i] < 0)&&(fCharge[j] > 0))
203                 fNpn[ibin] += 1.;
204             }
205         }
206     }//case 0
207    if(fAnalysisType==1)
208     {
209       for(i = 1; i < fNtrack; i++)
210         {
211           for(j = 0; j < i; j++)
212             {
213               Double_t p1 = sqrt(pow(fV[i].Px(),2) + pow(fV[i].Py(),2) + pow(fV[i].Pz(),2)); 
214               Double_t p2 = sqrt(pow(fV[j].Px(),2) + pow(fV[j].Py(),2) + pow(fV[j].Pz(),2));
215               Double_t eta1 = 0.5*log((p1 - fV[i].Pz())/(p1 + fV[i].Pz())); 
216               Double_t eta2 = 0.5*log((p2 - fV[j].Pz())/(p2 + fV[j].Pz())); 
217               Double_t deta = TMath::Abs(eta1 - eta2);
218               ibin = Int_t(deta/fP2Step);
219               if((fCharge[i] > 0)&&(fCharge[j] > 0))
220                 fNpp[ibin] += 1.;
221               if((fCharge[i] < 0)&&(fCharge[j] < 0))
222                 fNnn[ibin] += 1.;
223               if((fCharge[i] > 0)&&(fCharge[j] < 0))
224                 fNpn[ibin] += 1.;
225               if((fCharge[i] < 0)&&(fCharge[j] > 0))
226                 fNpn[ibin] += 1.;
227             }
228         }
229     }//case 1
230    if(fAnalysisType==2)
231     {
232       for(i = 1; i < fNtrack; i++)
233         {
234           for(j = 0; j < i; j++)
235             {
236               Double_t eTot = fV[i].E() + fV[j].E();
237               Double_t pxTot = fV[i].Px() + fV[j].Px();
238               Double_t pyTot = fV[i].Py() + fV[j].Py();
239               Double_t pzTot = fV[i].Pz() + fV[j].Pz();
240               Double_t q0Tot = fV[i].E() - fV[j].E();
241               Double_t qzTot = fV[i].Pz() - fV[j].Pz();
242               Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
243               Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
244               Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/sqrt(snn + pow(ptTot,2));
245               ibin = Int_t(qLong/fP2Step);
246               //cout<<ibin<<endl;
247               if((fCharge[i] > 0)&&(fCharge[j] > 0))
248                 fNpp[ibin] += 1.;
249               if((fCharge[i] < 0)&&(fCharge[j] < 0))
250                 fNnn[ibin] += 1.;
251               if((fCharge[i] > 0)&&(fCharge[j] < 0))
252                 fNpn[ibin] += 1.;
253               if((fCharge[i] < 0)&&(fCharge[j] > 0))
254                 fNpn[ibin] += 1.;
255             }
256         }
257     }//case 2
258   if(fAnalysisType==3)
259     {
260       for(i = 1; i < fNtrack; i++)
261         {
262           for(j = 0; j < i; j++)
263             {
264               Double_t eTot = fV[i].E() + fV[j].E();
265               Double_t pxTot = fV[i].Px() + fV[j].Px();
266               Double_t pyTot = fV[i].Py() + fV[j].Py();
267               Double_t pzTot = fV[i].Pz() + fV[j].Pz();
268               Double_t qxTot = fV[i].Px() - fV[j].Px();
269               Double_t qyTot = fV[i].Py() - fV[j].Py();
270               Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
271               Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
272               Double_t qOut = sqrt(snn/(snn + pow(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
273               ibin = Int_t(qOut/fP2Step);
274               //cout<<ibin<<endl;
275               if((fCharge[i] > 0)&&(fCharge[j] > 0))
276                 fNpp[ibin] += 1.;
277               if((fCharge[i] < 0)&&(fCharge[j] < 0))
278                 fNnn[ibin] += 1.;
279               if((fCharge[i] > 0)&&(fCharge[j] < 0))
280                 fNpn[ibin] += 1.;
281               if((fCharge[i] < 0)&&(fCharge[j] > 0))
282                 fNpn[ibin] += 1.;
283             }
284         }
285     }//case 3
286   if(fAnalysisType==4)
287     {
288       for(i = 1; i < fNtrack; i++)
289         {
290           for(j = 0; j < i; j++)
291             {
292               Double_t pxTot = fV[i].Px() + fV[j].Px();
293               Double_t pyTot = fV[i].Py() + fV[j].Py();
294               Double_t qxTot = fV[i].Px() - fV[j].Px();
295               Double_t qyTot = fV[i].Py() - fV[j].Py();
296               Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
297               Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
298               ibin = Int_t(qSide/fP2Step);
299               //cout<<ibin<<endl;
300               if((fCharge[i] > 0)&&(fCharge[j] > 0))
301                 fNpp[ibin] += 1.;
302               if((fCharge[i] < 0)&&(fCharge[j] < 0))
303                 fNnn[ibin] += 1.;
304               if((fCharge[i] > 0)&&(fCharge[j] < 0))
305                 fNpn[ibin] += 1.;
306               if((fCharge[i] < 0)&&(fCharge[j] > 0))
307                 fNpn[ibin] += 1.;
308             }
309         }
310     }//case 4
311    if(fAnalysisType==5)
312     {
313       for(i = 1; i < fNtrack; i++)
314         {
315           for(j = 0; j < i; j++)
316             {
317               Double_t q0Tot = fV[i].E() - fV[j].E();
318               Double_t qxTot = fV[i].Px() - fV[j].Px();
319               Double_t qyTot = fV[i].Py() - fV[j].Py();
320               Double_t qzTot = fV[i].Pz() - fV[j].Pz();
321               Double_t qInv = sqrt(TMath::Abs(-pow(q0Tot,2) +pow(qxTot,2) +pow(qyTot,2) +pow(qzTot,2)));
322               ibin = Int_t(qInv/fP2Step);
323               //cout<<ibin<<endl;
324               if((fCharge[i] > 0)&&(fCharge[j] > 0))
325                 fNpp[ibin] += 1.;
326               if((fCharge[i] < 0)&&(fCharge[j] < 0))
327                 fNnn[ibin] += 1.;
328               if((fCharge[i] > 0)&&(fCharge[j] < 0))
329                 fNpn[ibin] += 1.;
330               if((fCharge[i] < 0)&&(fCharge[j] > 0))
331                 fNpn[ibin] += 1.;
332             }
333         }
334     }//case 5   
335   if(fAnalysisType==6)
336     {
337       for(i = 1; i < fNtrack; i++)
338         {
339           for(j = 0; j < i; j++)
340             {
341               Double_t phi1 = TMath::ATan(fV[i].Py()/fV[i].Px())*180.0/TMath::Pi();
342               Double_t phi2 = TMath::ATan(fV[j].Py()/fV[j].Px())*180.0/TMath::Pi();
343               Double_t dphi = TMath::Abs(phi1 - phi2);
344               ibin = Int_t(dphi/fP2Step);
345               if((fCharge[i] > 0)&&(fCharge[j] > 0))
346                 fNpp[ibin] += 1.;
347               if((fCharge[i] < 0)&&(fCharge[j] < 0))
348                 fNnn[ibin] += 1.;
349               if((fCharge[i] > 0)&&(fCharge[j] < 0))
350                 fNpn[ibin] += 1.;
351               if((fCharge[i] < 0)&&(fCharge[j] > 0))
352                 fNpn[ibin] += 1.;
353             }
354         }
355     }//case 6
356 }
357
358 //----------------------------------------//
359 Double_t AliBalance::GetBalance(Int_t p2)
360 {
361   // Returns the value of the balance function in bin p2
362   fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
363
364   return fB[p2];
365 }
366     
367 //----------------------------------------//
368 Double_t AliBalance::GetError(Int_t p2)
369 {
370   // Returns the error on the BF value for bin p2
371   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))/fP2Step;
372
373   return ferror[p2];
374 }
375
376 //----------------------------------------//
377 void AliBalance::PrintResults()
378 {
379   // Prints the results
380   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
381   Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
382   Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
383   Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
384   Double_t deltaBalP2 = 0.0, integral = 0.0;
385   Double_t deltaErrorNew = 0.0;
386
387   cout<<"=================================================="<<endl;
388   for(Int_t i = 0; i < fNumberOfBins; i++)
389     { 
390       x[i] = fP2Start + fP2Step*i + fP2Step/2 ;
391       cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
392     } 
393   cout<<"=================================================="<<endl;
394   for(Int_t i = 1; i < fNumberOfBins; i++)
395     {
396       fSumXi += x[i];
397       fSumBi += fB[i];
398       fSumBiXi += fB[i]*x[i];
399       fSumBiXi2 += fB[i]*pow(x[i],2);
400       fSumBi2Xi2 += pow(fB[i],2)*pow(x[i],2);
401       fSumDeltaBi2 +=  pow(ferror[i],2) ;
402       fSumXi2DeltaBi2 += pow(x[i],2) * pow(ferror[i],2) ;
403       
404       deltaBalP2 += fP2Step*pow(ferror[i],2) ;
405       integral += fP2Step*fB[i] ;
406     }
407   for(Int_t i = 1; i < fNumberOfBins; i++)
408     {
409       deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/pow(fSumBi,2);
410     }
411   Double_t integralError = sqrt(deltaBalP2) ;
412   
413   Double_t delta = fSumBiXi / fSumBi ;
414   Double_t deltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) ) ;
415  
416   cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
417   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
418   cout<<"New error: "<<deltaErrorNew<<endl;
419   cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
420   cout<<"=================================================="<<endl;
421 }
422