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