5e4a6441190bd80fcfd0cb4a9ca6541140be6c15
[u/mrichter/AliRoot.git] / PWG2 / EBYE / 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 <TAxis.h>
27 #include <TLorentzVector.h>
28 #include <TObjArray.h>
29 #include <TGraphErrors.h>
30 #include <TString.h>
31 #include <TH1F.h>
32
33 #include "AliVParticle.h"
34 #include "AliMCParticle.h"
35 #include "AliESDtrack.h"
36 #include "AliAODTrack.h"
37
38 #include "AliBalance.h"
39
40 ClassImp(AliBalance)
41
42 //____________________________________________________________________//
43 AliBalance::AliBalance() :
44   TObject(), 
45   fAnalysisLevel("ESD"), fNumberOfBins(0),
46   fAnalysisType(0), fAnalyzedEvents(0), fP2Start(0),
47   fP2Stop(0), fP2Step(0), fNn(0), fNp(0),
48   fHistfNnn(new TH1F("fHistfNnn","(--) component;;Entries",
49                      fNumberOfBins,fP2Start,fP2Stop)),
50   fHistfNpp(new TH1F("fHistfNpp","(++) component;;Entries",
51                      fNumberOfBins,fP2Start,fP2Stop)),
52   fHistfNpn(new TH1F("fHistfNpn","(+-) component;;Entries",
53                      fNumberOfBins,fP2Start,fP2Stop)) {
54   // Default constructor
55   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
56     fNpp[i] = .0;
57     fNnn[i] = .0;
58     fNpn[i] = .0;
59     fB[i] = 0.0;
60     ferror[i] = 0.0;
61   } 
62
63   switch(fAnalysisType) {
64   case 0:
65     fHistfNnn->GetXaxis()->SetTitle("#Delta y");
66     fHistfNpp->GetXaxis()->SetTitle("#Delta y");
67     fHistfNpn->GetXaxis()->SetTitle("#Delta y");
68     break;
69   case 1:
70     fHistfNnn->GetXaxis()->SetTitle("#Delta #eta");
71     fHistfNpp->GetXaxis()->SetTitle("#Delta #eta");
72     fHistfNpn->GetXaxis()->SetTitle("#Delta #eta");
73     break;
74   case 2:
75     fHistfNnn->GetXaxis()->SetTitle("q_{long} (GeV/c)");
76     fHistfNpp->GetXaxis()->SetTitle("q_{long} (GeV/c)");
77     fHistfNpn->GetXaxis()->SetTitle("q_{long} (GeV/c)");
78     break;
79   case 3:
80     fHistfNnn->GetXaxis()->SetTitle("q_{out} (GeV/c)");
81     fHistfNpp->GetXaxis()->SetTitle("q_{out} (GeV/c)");
82     fHistfNpn->GetXaxis()->SetTitle("q_{out} (GeV/c)");
83     break;
84   case 4:
85     fHistfNnn->GetXaxis()->SetTitle("q_{side} (GeV/c)");
86     fHistfNpp->GetXaxis()->SetTitle("q_{side} (GeV/c)");
87     fHistfNpn->GetXaxis()->SetTitle("q_{side} (GeV/c)");
88     break;
89   case 5:
90     fHistfNnn->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
91     fHistfNpp->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
92     fHistfNpn->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
93     break;
94   case 6:
95     fHistfNnn->GetXaxis()->SetTitle("#Delta #phi");
96     fHistfNpp->GetXaxis()->SetTitle("#Delta #phi");
97     fHistfNpn->GetXaxis()->SetTitle("#Delta #phi");
98     break;
99   default:
100     break;
101   }
102 }
103
104 //____________________________________________________________________//
105 AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins) :
106   TObject(), fAnalysisLevel("ESD"),
107   fNumberOfBins(p2Bins), fAnalysisType(0), 
108   fAnalyzedEvents(0), fP2Start(p2Start), fP2Stop(p2Stop), 
109   fP2Step(TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins), 
110   fNn(0), fNp(0),
111   fHistfNnn(new TH1F("fHistfNnn","(--) component;;Entries",
112                      fNumberOfBins,fP2Start,fP2Stop)),
113   fHistfNpp(new TH1F("fHistfNpp","(++) component;;Entries",
114                      fNumberOfBins,fP2Start,fP2Stop)),
115   fHistfNpn(new TH1F("fHistfNpn","(+-) component;;Entries",
116                      fNumberOfBins,fP2Start,fP2Stop)) {
117   // Constructor
118   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
119     fNpp[i] = .0;
120     fNnn[i] = .0;
121     fNpn[i] = .0;
122     fB[i] = 0.0;
123     ferror[i] = 0.0;
124   } 
125
126   switch(fAnalysisType) {
127   case 0:
128     fHistfNnn->GetXaxis()->SetTitle("#Delta y");
129     fHistfNpp->GetXaxis()->SetTitle("#Delta y");
130     fHistfNpn->GetXaxis()->SetTitle("#Delta y");
131     break;
132   case 1:
133     fHistfNnn->GetXaxis()->SetTitle("#Delta #eta");
134     fHistfNpp->GetXaxis()->SetTitle("#Delta #eta");
135     fHistfNpn->GetXaxis()->SetTitle("#Delta #eta");
136     break;
137   case 2:
138     fHistfNnn->GetXaxis()->SetTitle("q_{long} (GeV/c)");
139     fHistfNpp->GetXaxis()->SetTitle("q_{long} (GeV/c)");
140     fHistfNpn->GetXaxis()->SetTitle("q_{long} (GeV/c)");
141     break;
142   case 3:
143     fHistfNnn->GetXaxis()->SetTitle("q_{out} (GeV/c)");
144     fHistfNpp->GetXaxis()->SetTitle("q_{out} (GeV/c)");
145     fHistfNpn->GetXaxis()->SetTitle("q_{out} (GeV/c)");
146     break;
147   case 4:
148     fHistfNnn->GetXaxis()->SetTitle("q_{side} (GeV/c)");
149     fHistfNpp->GetXaxis()->SetTitle("q_{side} (GeV/c)");
150     fHistfNpn->GetXaxis()->SetTitle("q_{side} (GeV/c)");
151     break;
152   case 5:
153     fHistfNnn->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
154     fHistfNpp->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
155     fHistfNpn->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
156     break;
157   case 6:
158     fHistfNnn->GetXaxis()->SetTitle("#Delta #phi");
159     fHistfNpp->GetXaxis()->SetTitle("#Delta #phi");
160     fHistfNpn->GetXaxis()->SetTitle("#Delta #phi");
161     break;
162   default:
163     break;
164   }
165 }
166
167 //____________________________________________________________________//
168 AliBalance::AliBalance(const AliBalance& balance):
169   TObject(balance), fAnalysisLevel(balance.fAnalysisLevel),
170   fNumberOfBins(balance.fNumberOfBins),
171   fAnalysisType(balance.fAnalysisType),
172   fAnalyzedEvents(balance.fAnalyzedEvents),
173   fP2Start(balance.fP2Start),
174   fP2Stop(balance.fP2Stop),
175   fP2Step(balance.fP2Step),
176   fNn(balance.fNn),
177   fNp(balance.fNp),
178   fHistfNnn(balance.fHistfNnn), 
179   fHistfNpp(balance.fHistfNpp), 
180   fHistfNpn(balance.fHistfNpn) {
181   //copy constructor
182   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
183     fNpp[i] = .0;
184     fNnn[i] = .0;
185     fNpn[i] = .0;
186     fB[i] = 0.0;
187     ferror[i] = 0.0;
188   } 
189 }
190
191 //____________________________________________________________________//
192 AliBalance::~AliBalance() {
193   // Destructor
194   if(fHistfNnn) delete fHistfNnn;
195   if(fHistfNpp) delete fHistfNpp;
196   if(fHistfNpn) delete fHistfNpn;
197 }
198
199 //____________________________________________________________________//
200 void AliBalance::SetNnn(Double_t *nn) {
201   // Setter of the Nnn term
202   for(Int_t i = 0; i < fNumberOfBins; i++) fNnn[i] = nn[i];
203 }
204
205 //____________________________________________________________________//
206 void AliBalance::SetNpp(Double_t *pp) {
207   // Setter of the Npp term
208   for(Int_t i = 0; i < fNumberOfBins; i++) fNpp[i] = pp[i];
209 }
210
211 //____________________________________________________________________//
212 void AliBalance::SetNpn(Double_t *pn) {
213   // Setter of the Npn term
214   for(Int_t i = 0; i < fNumberOfBins; i++) fNpn[i] = pn[i];
215 }
216
217 //____________________________________________________________________//
218 void AliBalance::SetNumberOfBins(Int_t ibins) {
219   // Sets the number of bins for the analyzed interval
220   fNumberOfBins = ibins;
221 }
222
223 //____________________________________________________________________//
224 void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop) {
225   // Sets the analyzed interval. 
226   // The analysis variable is set by SetAnalysisType
227   fP2Start = p2Start;
228   fP2Stop = p2Stop;
229   fP2Step = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins;
230 }
231
232 //____________________________________________________________________//
233 void AliBalance::SetAnalysisType(Int_t iType) {
234   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
235   this->fAnalysisType = iType; 
236   if(fAnalysisType==0) {
237     cout<<" ====================== "<<endl;
238     cout<<"||Analysis selected: y||"<<endl;
239     cout<<" ====================== "<<endl;
240   } 
241   else if(fAnalysisType==1) {
242     cout<<" ======================== "<<endl;
243     cout<<"||Analysis selected: eta||"<<endl;
244     cout<<" ======================== "<<endl;
245   }
246   else if(fAnalysisType==2) {
247     cout<<" ========================== "<<endl;
248     cout<<"||Analysis selected: Qlong||"<<endl;
249     cout<<" ========================== "<<endl;
250   }
251   else if(fAnalysisType==3) {
252     cout<<" ========================= "<<endl;
253     cout<<"||Analysis selected: Qout||"<<endl;
254     cout<<" ========================= "<<endl;
255   }
256   else if(fAnalysisType==4) {
257     cout<<" ========================== "<<endl;
258     cout<<"||Analysis selected: Qside||"<<endl;
259     cout<<" ========================== "<<endl;
260   }
261   else if(fAnalysisType==5) {
262     cout<<" ========================= "<<endl;
263     cout<<"||Analysis selected: Qinv||"<<endl;
264     cout<<" ========================= "<<endl;
265   }
266   else if(fAnalysisType==6) {
267     cout<<" ======================== "<<endl;
268     cout<<"||Analysis selected: phi||"<<endl;
269     cout<<" ======================== "<<endl;
270   }
271   else {
272     cout<<"Selection of analysis mode failed!!!"<<endl;
273     cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
274     abort();
275   }
276 }
277
278 //____________________________________________________________________//
279 void AliBalance::PrintAnalysisSettings() {
280   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
281   TString analysisType;
282   switch(fAnalysisType) {
283   case 0:
284     analysisType = "Rapidity"; 
285     break;
286   case 1:
287     analysisType = "Pseudo-rapidity"; 
288     break;
289   case 2:
290     analysisType = "Qlong"; 
291     break;
292   case 3:
293     analysisType = "Qout"; 
294     break;
295   case 4:
296     analysisType = "Qside"; 
297     break;
298   case 5:
299     analysisType = "Qinv"; 
300     break;
301   case 6:
302     analysisType = "Phi"; 
303     break;
304   default:
305     break;
306   }
307   
308   Printf("======================================");
309   Printf("Analysis level: %s",fAnalysisLevel.Data());
310   Printf("Analysis type: %s",analysisType.Data());
311   Printf("Analyzed interval (min.): %lf",fP2Start);
312   Printf("Analyzed interval (max.): %lf",fP2Stop);
313   Printf("Number of bins: %d",fNumberOfBins);
314   Printf("Step: %lf",fP2Step);
315   Printf("======================================");
316 }
317
318 //____________________________________________________________________//
319 void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
320   // Calculates the balance function
321   fAnalyzedEvents++;
322   Int_t i = 0 , j = 0;
323   Int_t ibin = 0;
324   
325   AliVParticle* track = 0;
326   AliVParticle* track1 = 0;
327   AliVParticle* track2 = 0;
328     
329   //Printf("(AliBalance) Number of tracks: %d",gTrackArray->GetEntries());
330   Int_t gNtrack = gTrackArray->GetEntries();
331   for(i = 0; i < gNtrack; i++) {
332     if(fAnalysisLevel == "ESD")
333       track = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
334     else if(fAnalysisLevel == "AOD")
335       track = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
336     else if(fAnalysisLevel == "MC")
337       track = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
338     Short_t charge = track->Charge();
339     if(charge > 0) fNp += 1.;
340     if(charge < 0) fNn += 1.;
341   }
342   //Printf("Np: %lf - Nn: %lf",fNp,fNn);
343
344   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
345   if(fAnalysisType==0) {
346     for(i = 1; i < gNtrack; i++) {
347       if(fAnalysisLevel == "ESD")
348         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
349       else if(fAnalysisLevel == "AOD")
350         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
351       else if(fAnalysisLevel == "MC")
352         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
353       Short_t charge1 = track1->Charge();
354       Double_t pZ1 = track1->Pz();
355       Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
356                                      TMath::Power(track1->M(),2));
357       for(j = 0; j < i; j++) {
358         if(fAnalysisLevel == "ESD")
359           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
360         else if(fAnalysisLevel == "AOD")
361           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
362         else if(fAnalysisLevel == "MC")
363           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
364         Short_t charge2 = track2->Charge();
365         Double_t pZ2 = track2->Pz();
366         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
367                                        TMath::Power(track2->M(),2));
368
369         Double_t rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1)); 
370         Double_t rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2)); 
371         Double_t dy = TMath::Abs(rap1 - rap2);
372         ibin = Int_t(dy/fP2Step);
373         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
374         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
375         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
376         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
377       }
378     }
379   }//case 0
380   if(fAnalysisType==1) {
381     for(i = 1; i < gNtrack; i++) {
382       if(fAnalysisLevel == "ESD")
383         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
384       else if(fAnalysisLevel == "AOD")
385         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
386       else if(fAnalysisLevel == "MC")
387         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
388       Short_t charge1 = track1->Charge();
389       Double_t pZ1 = track1->Pz();
390       Double_t p1 = track1->P();
391       for(j = 0; j < i; j++) {
392         if(fAnalysisLevel == "ESD")
393           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
394         else if(fAnalysisLevel == "AOD")
395           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
396         else if(fAnalysisLevel == "MC")
397           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
398         Short_t charge2 = track2->Charge();
399         Double_t pZ2 = track2->Pz();
400         Double_t p2 = track2->P();
401         Double_t eta1 = 0.5*log((p1 + pZ1)/(p1 - pZ1)); 
402         Double_t eta2 = 0.5*log((p2 + pZ2)/(p2 - pZ2)); 
403         Double_t deta = TMath::Abs(eta1 - eta2);
404         ibin = Int_t(deta/fP2Step);
405         
406         if((charge1 > 0.)&&(charge2 > 0.)) fNpp[ibin] += 1.;
407         if((charge1 < 0.)&&(charge2 < 0.)) fNnn[ibin] += 1.;
408         if((charge1 > 0.)&&(charge2 < 0.)) fNpn[ibin] += 1.;
409         if((charge1 < 0.)&&(charge2 > 0.)) fNpn[ibin] += 1.;
410         //Printf("charge1: %d - eta1: %lf - charge2: %d - eta2: %lf - deta: %lf - ibin: %d - fNpp: %lf - fNnn: %lf - fNpn: %lf",charge1,eta1,charge2,eta2,deta,ibin,fNpp[ibin],fNnn[ibin],fNpn[ibin]);      
411       }
412     }
413   }//case 1
414   if(fAnalysisType==2) {
415     for(i = 1; i < gNtrack; i++) {
416       if(fAnalysisLevel == "ESD")
417         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
418       else if(fAnalysisLevel == "AOD")
419         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
420       else if(fAnalysisLevel == "MC")
421         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
422       Short_t charge1 = track1->Charge();
423       Double_t pX1 = track1->Px();
424       Double_t pY1 = track1->Py();
425       Double_t pZ1 = track1->Pz();
426       Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
427                                      TMath::Power(track1->M(),2));
428       for(j = 0; j < i; j++) {
429         if(fAnalysisLevel == "ESD")
430           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
431         else if(fAnalysisLevel == "AOD")
432           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
433         else if(fAnalysisLevel == "MC")
434           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
435         Short_t charge2 = track2->Charge();
436         Double_t pX2 = track2->Px();
437         Double_t pY2 = track2->Py();
438         Double_t pZ2 = track2->Pz();
439         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
440                                        TMath::Power(track2->M(),2));
441         Double_t eTot = energy1 + energy2;
442         Double_t pxTot = pX1 + pX2;
443         Double_t pyTot = pY1 + pY2;
444         Double_t pzTot = pZ1 + pZ2;
445         Double_t q0Tot = energy1 - energy2;
446         Double_t qzTot = pZ1 - pZ2;
447         Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
448         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
449         Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
450         ibin = Int_t(qLong/fP2Step);
451         //cout<<ibin<<endl;
452         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
453         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
454         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
455         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
456       }
457     }
458   }//case 2
459   if(fAnalysisType==3) {
460     for(i = 1; i < gNtrack; i++) {
461       if(fAnalysisLevel == "ESD")
462         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
463       else if(fAnalysisLevel == "AOD")
464         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
465       else if(fAnalysisLevel == "MC")
466         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
467       Short_t charge1 = track1->Charge();
468       Double_t pX1 = track1->Px();
469       Double_t pY1 = track1->Py();
470       Double_t pZ1 = track1->Pz();
471       Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
472                                      TMath::Power(track1->M(),2));
473       for(j = 0; j < i; j++) {
474         if(fAnalysisLevel == "ESD")
475           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
476         else if(fAnalysisLevel == "AOD")
477           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
478         else if(fAnalysisLevel == "MC")
479           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
480         Short_t charge2 = track2->Charge();
481         Double_t pX2 = track2->Px();
482         Double_t pY2 = track2->Py();
483         Double_t pZ2 = track2->Pz();
484         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
485                                        TMath::Power(track2->M(),2));
486         Double_t eTot = energy1 + energy2;
487         Double_t pxTot = pX1 + pX2;
488         Double_t pyTot = pY1 + pY2;
489         Double_t pzTot = pZ1 + pZ2;
490         Double_t qxTot = pX1 - pX2;
491         Double_t qyTot = pY1 - pY2;
492         Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
493         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
494         Double_t qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
495         ibin = Int_t(qOut/fP2Step);
496         //cout<<ibin<<endl;
497         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
498         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
499         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
500         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
501       }
502     }
503   }//case 3
504   if(fAnalysisType==4) {
505     for(i = 1; i < gNtrack; i++) {
506       if(fAnalysisLevel == "ESD")
507         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
508       else if(fAnalysisLevel == "AOD")
509         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
510       else if(fAnalysisLevel == "MC")
511         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
512       Short_t charge1 = track1->Charge();
513       Double_t pX1 = track1->Px();
514       Double_t pY1 = track1->Py();
515       //Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
516       //TMath::Power(track1->M(),2));
517       for(j = 0; j < i; j++) {
518         if(fAnalysisLevel == "ESD")
519           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
520         else if(fAnalysisLevel == "AOD")
521           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
522         else if(fAnalysisLevel == "MC")
523           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
524         Short_t charge2 = track2->Charge();
525         Double_t pX2 = track2->Px();
526         Double_t pY2 = track2->Py();
527         //Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
528         //TMath::Power(track2->M(),2));
529         //Double_t eTot = energy1 + energy2;
530         Double_t pxTot = pX1 + pX2;
531         Double_t pyTot = pY1 + pY2;
532         Double_t qxTot = pX1 - pX2;
533         Double_t qyTot = pY1 - pY2;
534         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
535         Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
536         ibin = Int_t(qSide/fP2Step);
537         //cout<<ibin<<endl;
538         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
539         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
540         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
541         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
542       }
543     }
544   }//case 4
545   if(fAnalysisType==5) {
546     for(i = 1; i < gNtrack; i++) {
547       if(fAnalysisLevel == "ESD")
548         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
549       else if(fAnalysisLevel == "AOD")
550         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
551       else if(fAnalysisLevel == "MC")
552         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
553       Short_t charge1 = track1->Charge();
554       Double_t pX1 = track1->Px();
555       Double_t pY1 = track1->Py();
556       Double_t pZ1 = track1->Pz();
557       Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
558                                      TMath::Power(track1->M(),2));
559       for(j = 0; j < i; j++) {
560         if(fAnalysisLevel == "ESD")
561           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
562         else if(fAnalysisLevel == "AOD")
563           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
564         else if(fAnalysisLevel == "MC")
565           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
566         Short_t charge2 = track2->Charge();
567         Double_t pX2 = track2->Px();
568         Double_t pY2 = track2->Py();
569         Double_t pZ2 = track2->Pz();
570         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
571                                        TMath::Power(track2->M(),2));
572         Double_t q0Tot = energy1 - energy2;
573         Double_t qxTot = pX1 - pX2;
574         Double_t qyTot = pY1 - pY2;
575         Double_t qzTot = pZ1 - pZ2;
576         Double_t qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
577         ibin = Int_t(qInv/fP2Step);
578         //cout<<ibin<<endl;
579         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
580         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
581         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
582         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
583       }
584     }
585   }//case 5     
586   if(fAnalysisType==6) {
587     for(i = 1; i < gNtrack; i++) {
588       if(fAnalysisLevel == "ESD")
589         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
590       else if(fAnalysisLevel == "AOD")
591         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
592       else if(fAnalysisLevel == "MC")
593         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
594       Short_t charge1 = track1->Charge();
595       Double_t phi1 = track1->Phi();
596       for(j = 0; j < i; j++) {
597         if(fAnalysisLevel == "ESD")
598           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
599         else if(fAnalysisLevel == "AOD")
600           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
601         else if(fAnalysisLevel == "MC")
602           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
603         Short_t charge2 = track2->Charge();
604         Double_t phi2 = track2->Phi();
605         Double_t dphi = TMath::Abs(phi1 - phi2);
606         ibin = Int_t(dphi/fP2Step);
607         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
608         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
609         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
610         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
611       }
612     }
613   }//case 6
614
615   /*for(Int_t i = 0; i < fNumberOfBins; i++) 
616     Printf("bin: %d - Npp: %lf - Nnn: %lf - Nnp: %lf - Npn: %lf",i,fNpp[i],fNnn[i],fNpn[i],fNpn[i]);*/
617 }
618
619 //____________________________________________________________________//
620 Double_t AliBalance::GetBalance(Int_t p2) {
621   // Returns the value of the balance function in bin p2
622   fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
623   
624   return fB[p2];
625 }
626     
627 //____________________________________________________________________//
628 Double_t AliBalance::GetError(Int_t p2) {
629   // Returns the error on the BF value for bin p2
630   ferror[p2] = TMath::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])*TMath::Power((0.5/Double_t(fNp) + 0.5/Double_t(fNn)),2))/fP2Step;
631
632   return ferror[p2];
633 }
634
635 //____________________________________________________________________//
636 void AliBalance::PrintResults() {
637   // Prints the results
638   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
639   Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
640   Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
641   Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
642   Double_t deltaBalP2 = 0.0, integral = 0.0;
643   Double_t deltaErrorNew = 0.0;
644
645   cout<<"=================================================="<<endl;
646   for(Int_t i = 0; i < fNumberOfBins; i++) { 
647     //x[i] = fP2Start + fP2Step*i + fP2Step/2;
648     x[i] = fP2Step*i + fP2Step/2;
649     cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
650   } 
651   cout<<"=================================================="<<endl;
652   for(Int_t i = 1; i < fNumberOfBins; i++) {
653     fSumXi += x[i];
654     fSumBi += fB[i];
655     fSumBiXi += fB[i]*x[i];
656     fSumBiXi2 += fB[i]*TMath::Power(x[i],2);
657     fSumBi2Xi2 += TMath::Power(fB[i],2)*TMath::Power(x[i],2);
658     fSumDeltaBi2 +=  TMath::Power(ferror[i],2);
659     fSumXi2DeltaBi2 += TMath::Power(x[i],2) * TMath::Power(ferror[i],2);
660     
661     deltaBalP2 += fP2Step*TMath::Power(ferror[i],2);
662     integral += fP2Step*fB[i];
663   }
664   for(Int_t i = 1; i < fNumberOfBins; i++) deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/TMath::Power(fSumBi,2);
665    
666   Double_t integralError = TMath::Sqrt(deltaBalP2);
667   
668   Double_t delta = fSumBiXi / fSumBi;
669   Double_t deltaError = (fSumBiXi / fSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + TMath::Power((fSumDeltaBi2/fSumBi),2) );
670  
671   cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
672   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
673   cout<<"New error: "<<deltaErrorNew<<endl;
674   cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
675   cout<<"=================================================="<<endl;
676 }
677   
678 //____________________________________________________________________//
679 TGraphErrors *AliBalance::DrawBalance() {
680   // Draws the BF
681   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
682   Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
683   Double_t b[MAXIMUM_NUMBER_OF_STEPS];
684   Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
685
686   if((fNp == 0)||(fNn == 0)) {
687     cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
688     cout<<"Aborting....."<<endl;
689     abort();
690   }
691   
692   for(Int_t i = 0; i < fNumberOfBins; i++) {
693     b[i] = GetBalance(i);
694     ber[i] = GetError(i);
695     //x[i] = fP2Start + fP2Step*i + fP2Step/2;
696     x[i] = fP2Step*i + fP2Step/2;
697     xer[i] = 0.0;
698   }
699   
700   TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber);
701   gr->SetMarkerStyle(25);
702   gr->GetXaxis()->SetTitleColor(1);
703   if(fAnalysisType==0) {
704     gr->GetXaxis()->SetTitle("#Delta y");
705     gr->GetYaxis()->SetTitle("B(#Delta y)");
706   }
707   if(fAnalysisType==1) {
708     gr->GetXaxis()->SetTitle("#Delta #eta");
709     gr->GetYaxis()->SetTitle("B(#Delta #eta)");
710   }
711   if(fAnalysisType==2) {
712     gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
713     gr->GetYaxis()->SetTitle("B(q_{long}) [(GeV/c)^{-1}]");
714   }
715   if(fAnalysisType==3) {
716     gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
717     gr->GetYaxis()->SetTitle("B(q_{out}) [(GeV/c)^{-1}]");
718   }
719   if(fAnalysisType==4) {
720     gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
721     gr->GetYaxis()->SetTitle("B(q_{side}) [(GeV/c)^{-1}]");
722   }
723   if(fAnalysisType==5) {
724     gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
725     gr->GetYaxis()->SetTitle("B(q_{inv}) [(GeV/c)^{-1}]");
726   }
727   if(fAnalysisType==6) {
728     gr->GetXaxis()->SetTitle("#Delta #phi");
729     gr->GetYaxis()->SetTitle("B(#Delta #phi)");
730   }
731
732   return gr;
733 }
734
735 //____________________________________________________________________//
736 void AliBalance::Merge(AliBalance *b) {
737   //Merging function to be used for proof and grid
738   fNp += b->GetNp();
739   fNn += b->GetNn();
740   for(Int_t i = 0; i < fNumberOfBins; i++) {
741     fNnn[i] += b->GetNnn(i);
742     fNpp[i] += b->GetNpp(i);
743     fNpn[i] += b->GetNpn(i);
744   }
745 }