Fix in the print of the results (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   fNp = 0.0;
50   fNn = 0.0;
51   fP2_Start = 0.0;
52   fP2_Stop = 0.0;
53   fP2_Step = 0.0; 
54   fAnalysisType = 0;
55   fNumberOfBins = 0;
56   fNtrack = 0;
57   fAnalyzedEvents = 0;
58 }
59
60 //----------------------------------------//
61 AliBalance::AliBalance(Double_t P2_Start, Double_t P2_Stop, Int_t P2_bins)
62 {
63   this->fP2_Start = P2_Start;
64   this->fP2_Stop = P2_Stop;
65   this->fNumberOfBins = P2_bins;
66   this->fP2_Step = fabs(fP2_Start - fP2_Stop) / (Double_t)fNumberOfBins;
67 }
68
69 //----------------------------------------//
70 AliBalance::~AliBalance()
71 {
72 }
73
74 //----------------------------------------//
75 void AliBalance::SetNumberOfBins(Int_t ibins)
76 {
77   this->fNumberOfBins = ibins ;
78 }
79
80 //----------------------------------------//
81 void AliBalance::SetInterval(Double_t P2_Start, Double_t P2_Stop)
82 {
83   this->fP2_Start = P2_Start;
84   this->fP2_Stop = P2_Stop;
85   this->fP2_Step = fabs(P2_Start - P2_Stop) / (Double_t)fNumberOfBins;
86 }
87
88 //----------------------------------------//
89 void AliBalance::SetAnalysisType(Int_t iType)
90 {
91   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
92   this->fAnalysisType = iType; 
93   if(fAnalysisType==0)
94     {
95       cout<<" ====================== "<<endl;
96       cout<<"||Analysis selected: y||"<<endl;
97       cout<<" ====================== "<<endl;
98     } 
99   else if(fAnalysisType==1)
100     {
101       cout<<" ======================== "<<endl;
102       cout<<"||Analysis selected: eta||"<<endl;
103       cout<<" ======================== "<<endl;
104     }
105   else if(fAnalysisType==2)
106     {
107       cout<<" ========================== "<<endl;
108       cout<<"||Analysis selected: Qlong||"<<endl;
109       cout<<" ========================== "<<endl;
110     }
111   else if(fAnalysisType==3)
112     {
113       cout<<" ========================= "<<endl;
114       cout<<"||Analysis selected: Qout||"<<endl;
115       cout<<" ========================= "<<endl;
116     }
117   else if(fAnalysisType==4)
118     {
119       cout<<" ========================== "<<endl;
120       cout<<"||Analysis selected: Qside||"<<endl;
121       cout<<" ========================== "<<endl;
122     }
123   else if(fAnalysisType==5)
124     {
125       cout<<" ========================= "<<endl;
126       cout<<"||Analysis selected: Qinv||"<<endl;
127       cout<<" ========================= "<<endl;
128     }
129   else if(fAnalysisType==6)
130     {
131       cout<<" ======================== "<<endl;
132       cout<<"||Analysis selected: phi||"<<endl;
133       cout<<" ======================== "<<endl;
134     }
135   else
136     {
137       cout<<"Selection of analysis mode failed!!!"<<endl;
138       cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
139       abort();
140     }
141 }
142
143 //----------------------------------------//
144 void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim)
145 {
146   this->fV = P;
147   this->fCharge = charge;
148   fNtrack = dim;
149 }
150
151
152 //----------------------------------------//
153 void AliBalance::CalculateBalance()
154 {
155   fAnalyzedEvents++;
156   Int_t i = 0 , j = 0;
157   Int_t ibin = 0;
158
159   for(i = 0; i < fNtrack; i++)
160     {
161       if(fCharge[i] > 0)
162         fNp += 1.;
163       if(fCharge[i] < 0)
164         fNn += 1.;
165     }
166
167   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
168    if(fAnalysisType==0)
169     {
170       for(i = 1; i < fNtrack; i++)
171         {
172           for(j = 0; j < i; j++)
173             {
174               Double_t rap1 = 0.5*log((fV[i].E() - fV[i].Pz())/(fV[i].E() + fV[i].Pz())); 
175               Double_t rap2 = 0.5*log((fV[j].E() - fV[j].Pz())/(fV[j].E() + fV[j].Pz())); 
176               Double_t dy = TMath::Abs(rap1 - rap2);
177               ibin = Int_t(dy/fP2_Step);
178               if((fCharge[i] > 0)&&(fCharge[j] > 0))
179                 fNpp[ibin] += 1.;
180               if((fCharge[i] < 0)&&(fCharge[j] < 0))
181                 fNnn[ibin] += 1.;
182               if((fCharge[i] > 0)&&(fCharge[j] < 0))
183                 fNpn[ibin] += 1.;
184               if((fCharge[i] < 0)&&(fCharge[j] > 0))
185                 fNpn[ibin] += 1.;
186             }
187         }
188     }//case 0
189    if(fAnalysisType==1)
190     {
191       for(i = 1; i < fNtrack; i++)
192         {
193           for(j = 0; j < i; j++)
194             {
195               Double_t P1 = sqrt(pow(fV[i].Px(),2) + pow(fV[i].Py(),2) + pow(fV[i].Pz(),2)); 
196               Double_t P2 = sqrt(pow(fV[j].Px(),2) + pow(fV[j].Py(),2) + pow(fV[j].Pz(),2));
197               Double_t eta1 = 0.5*log((P1 - fV[i].Pz())/(P1 + fV[i].Pz())); 
198               Double_t eta2 = 0.5*log((P2 - fV[j].Pz())/(P2 + fV[j].Pz())); 
199               Double_t deta = TMath::Abs(eta1 - eta2);
200               ibin = Int_t(deta/fP2_Step);
201               if((fCharge[i] > 0)&&(fCharge[j] > 0))
202                 fNpp[ibin] += 1.;
203               if((fCharge[i] < 0)&&(fCharge[j] < 0))
204                 fNnn[ibin] += 1.;
205               if((fCharge[i] > 0)&&(fCharge[j] < 0))
206                 fNpn[ibin] += 1.;
207               if((fCharge[i] < 0)&&(fCharge[j] > 0))
208                 fNpn[ibin] += 1.;
209             }
210         }
211     }//case 1
212    if(fAnalysisType==2)
213     {
214       for(i = 1; i < fNtrack; i++)
215         {
216           for(j = 0; j < i; j++)
217             {
218               Double_t ETot = fV[i].E() + fV[j].E();
219               Double_t PxTot = fV[i].Px() + fV[j].Px();
220               Double_t PyTot = fV[i].Py() + fV[j].Py();
221               Double_t PzTot = fV[i].Pz() + fV[j].Pz();
222               Double_t Q0Tot = fV[i].E() - fV[j].E();
223               Double_t QzTot = fV[i].Pz() - fV[j].Pz();
224               Double_t Snn = pow(ETot,2) - pow(PxTot,2) - pow(PyTot,2) - pow(PzTot,2);
225               Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
226               Double_t Qlong = TMath::Abs(ETot*QzTot - PzTot*Q0Tot)/sqrt(Snn + pow(PtTot,2));
227               ibin = Int_t(Qlong/fP2_Step);
228               //cout<<ibin<<endl;
229               if((fCharge[i] > 0)&&(fCharge[j] > 0))
230                 fNpp[ibin] += 1.;
231               if((fCharge[i] < 0)&&(fCharge[j] < 0))
232                 fNnn[ibin] += 1.;
233               if((fCharge[i] > 0)&&(fCharge[j] < 0))
234                 fNpn[ibin] += 1.;
235               if((fCharge[i] < 0)&&(fCharge[j] > 0))
236                 fNpn[ibin] += 1.;
237             }
238         }
239     }//case 2
240   if(fAnalysisType==3)
241     {
242       for(i = 1; i < fNtrack; i++)
243         {
244           for(j = 0; j < i; j++)
245             {
246               Double_t ETot = fV[i].E() + fV[j].E();
247               Double_t PxTot = fV[i].Px() + fV[j].Px();
248               Double_t PyTot = fV[i].Py() + fV[j].Py();
249               Double_t PzTot = fV[i].Pz() + fV[j].Pz();
250               Double_t QxTot = fV[i].Px() - fV[j].Px();
251               Double_t QyTot = fV[i].Py() - fV[j].Py();
252               Double_t Snn = pow(ETot,2) - pow(PxTot,2) - pow(PyTot,2) - pow(PzTot,2);
253               Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
254               Double_t Qout = sqrt(Snn/(Snn + pow(PtTot,2))) * TMath::Abs(PxTot*QxTot + PyTot*QyTot)/PtTot;
255               ibin = Int_t(Qout/fP2_Step);
256               //cout<<ibin<<endl;
257               if((fCharge[i] > 0)&&(fCharge[j] > 0))
258                 fNpp[ibin] += 1.;
259               if((fCharge[i] < 0)&&(fCharge[j] < 0))
260                 fNnn[ibin] += 1.;
261               if((fCharge[i] > 0)&&(fCharge[j] < 0))
262                 fNpn[ibin] += 1.;
263               if((fCharge[i] < 0)&&(fCharge[j] > 0))
264                 fNpn[ibin] += 1.;
265             }
266         }
267     }//case 3
268   if(fAnalysisType==4)
269     {
270       for(i = 1; i < fNtrack; i++)
271         {
272           for(j = 0; j < i; j++)
273             {
274               Double_t PxTot = fV[i].Px() + fV[j].Px();
275               Double_t PyTot = fV[i].Py() + fV[j].Py();
276               Double_t QxTot = fV[i].Px() - fV[j].Px();
277               Double_t QyTot = fV[i].Py() - fV[j].Py();
278               Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
279               Double_t Qside = TMath::Abs(PxTot*QyTot - PyTot*QxTot)/PtTot;
280               ibin = Int_t(Qside/fP2_Step);
281               //cout<<ibin<<endl;
282               if((fCharge[i] > 0)&&(fCharge[j] > 0))
283                 fNpp[ibin] += 1.;
284               if((fCharge[i] < 0)&&(fCharge[j] < 0))
285                 fNnn[ibin] += 1.;
286               if((fCharge[i] > 0)&&(fCharge[j] < 0))
287                 fNpn[ibin] += 1.;
288               if((fCharge[i] < 0)&&(fCharge[j] > 0))
289                 fNpn[ibin] += 1.;
290             }
291         }
292     }//case 4
293    if(fAnalysisType==5)
294     {
295       for(i = 1; i < fNtrack; i++)
296         {
297           for(j = 0; j < i; j++)
298             {
299               Double_t Q0Tot = fV[i].E() - fV[j].E();
300               Double_t QxTot = fV[i].Px() - fV[j].Px();
301               Double_t QyTot = fV[i].Py() - fV[j].Py();
302               Double_t QzTot = fV[i].Pz() - fV[j].Pz();
303               Double_t Qinv = sqrt(TMath::Abs(-pow(Q0Tot,2) +pow(QxTot,2) +pow(QyTot,2) +pow(QzTot,2)));
304               ibin = Int_t(Qinv/fP2_Step);
305               //cout<<ibin<<endl;
306               if((fCharge[i] > 0)&&(fCharge[j] > 0))
307                 fNpp[ibin] += 1.;
308               if((fCharge[i] < 0)&&(fCharge[j] < 0))
309                 fNnn[ibin] += 1.;
310               if((fCharge[i] > 0)&&(fCharge[j] < 0))
311                 fNpn[ibin] += 1.;
312               if((fCharge[i] < 0)&&(fCharge[j] > 0))
313                 fNpn[ibin] += 1.;
314             }
315         }
316     }//case 5   
317   if(fAnalysisType==6)
318     {
319       for(i = 1; i < fNtrack; i++)
320         {
321           for(j = 0; j < i; j++)
322             {
323               Double_t phi1 = TMath::ATan(fV[i].Py()/fV[i].Px())*180.0/TMath::Pi();
324               Double_t phi2 = TMath::ATan(fV[j].Py()/fV[j].Px())*180.0/TMath::Pi();
325               Double_t dphi = TMath::Abs(phi1 - phi2);
326               ibin = Int_t(dphi/fP2_Step);
327               if((fCharge[i] > 0)&&(fCharge[j] > 0))
328                 fNpp[ibin] += 1.;
329               if((fCharge[i] < 0)&&(fCharge[j] < 0))
330                 fNnn[ibin] += 1.;
331               if((fCharge[i] > 0)&&(fCharge[j] < 0))
332                 fNpn[ibin] += 1.;
333               if((fCharge[i] < 0)&&(fCharge[j] > 0))
334                 fNpn[ibin] += 1.;
335             }
336         }
337     }//case 6
338 }
339
340 //----------------------------------------//
341 Double_t AliBalance::GetBalance(Int_t p2)
342 {
343   fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2_Step;
344
345   return fB[p2];
346 }
347     
348 //----------------------------------------//
349 Double_t AliBalance::GetError(Int_t p2)
350 {
351   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;
352
353   return ferror[p2];
354 }
355
356 //----------------------------------------//
357 void AliBalance::PrintResults()
358 {
359   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
360   Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
361   Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
362   Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
363   Double_t delta_bal_P2 = 0.0, Integral = 0.0;
364   Double_t DeltaErrorNew = 0.0;
365
366   cout<<"=================================================="<<endl;
367   for(Int_t i = 0; i < fNumberOfBins; i++)
368     { 
369       x[i] = fP2_Start + fP2_Step*i + fP2_Step/2 ;
370       cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
371     } 
372   cout<<"=================================================="<<endl;
373   for(Int_t i = 1; i < fNumberOfBins; i++)
374     {
375       fSumXi += x[i];
376       fSumBi += fB[i];
377       fSumBiXi += fB[i]*x[i];
378       fSumBiXi2 += fB[i]*pow(x[i],2);
379       fSumBi2Xi2 += pow(fB[i],2)*pow(x[i],2);
380       fSumDeltaBi2 +=  pow(ferror[i],2) ;
381       fSumXi2DeltaBi2 += pow(x[i],2) * pow(ferror[i],2) ;
382       
383       delta_bal_P2 += fP2_Step*pow(ferror[i],2) ;
384       Integral += fP2_Step*fB[i] ;
385     }
386   for(Int_t i = 1; i < fNumberOfBins; i++)
387     {
388       DeltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/pow(fSumBi,2);
389     }
390   Double_t IntegralError = sqrt(delta_bal_P2) ;
391   
392   Double_t Delta = fSumBiXi / fSumBi ;
393   Double_t DeltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) ) ;
394  
395   cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
396   cout<<"Width: "<<Delta<<"\t Error: "<<DeltaError<<endl;
397   cout<<"New error: "<<DeltaErrorNew<<endl;
398   cout<<"Interval: "<<Integral<<"\t Error: "<<IntegralError<<endl;
399   cout<<"=================================================="<<endl;
400 }
401