Adding the SetOwner
[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
339     if(track) {
340       Short_t charge = track->Charge();
341       if(charge > 0) fNp += 1.;
342       if(charge < 0) fNn += 1.;
343     }
344     else continue;
345   }
346   //Printf("Np: %lf - Nn: %lf",fNp,fNn);
347
348   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
349   if(fAnalysisType==0) {
350     for(i = 1; i < gNtrack; i++) {
351       if(fAnalysisLevel == "ESD")
352         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
353       else if(fAnalysisLevel == "AOD")
354         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
355       else if(fAnalysisLevel == "MC")
356         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
357
358       Short_t charge1 = 0;
359       Double_t pZ1 = 0., energy1 = 0.;
360       if(track1) {
361         charge1 = track1->Charge();
362         pZ1 = track1->Pz();
363         energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
364                               TMath::Power(track1->M(),2));
365       }
366       else continue;
367
368       for(j = 0; j < i; j++) {
369         if(fAnalysisLevel == "ESD")
370           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
371         else if(fAnalysisLevel == "AOD")
372           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
373         else if(fAnalysisLevel == "MC")
374           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
375
376         if(track2) {
377           Short_t charge2 = track2->Charge();
378           Double_t pZ2 = track2->Pz();
379           Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
380                                          TMath::Power(track2->M(),2));
381           
382           Double_t rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1)); 
383           Double_t rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2)); 
384           Double_t dy = TMath::Abs(rap1 - rap2);
385           ibin = Int_t(dy/fP2Step);
386           if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
387           if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
388           if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
389           if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
390         }
391         else continue;
392       }
393     }
394   }//case 0
395   if(fAnalysisType==1) {
396     for(i = 1; i < gNtrack; i++) {
397       if(fAnalysisLevel == "ESD")
398         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
399       else if(fAnalysisLevel == "AOD")
400         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
401       else if(fAnalysisLevel == "MC")
402         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
403
404       Short_t charge1 = 0;
405       Double_t pZ1 = 0., p1 = 0.;
406       if(track1) {
407         charge1 = track1->Charge();
408         pZ1 = track1->Pz();
409         p1 = track1->P();
410       }
411       else continue;
412
413       for(j = 0; j < i; j++) {
414         if(fAnalysisLevel == "ESD")
415           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
416         else if(fAnalysisLevel == "AOD")
417           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
418         else if(fAnalysisLevel == "MC")
419           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
420
421         if(track2) {
422           Short_t charge2 = track2->Charge();
423           Double_t pZ2 = track2->Pz();
424           Double_t p2 = track2->P();
425           Double_t eta1 = 0.5*log((p1 + pZ1)/(p1 - pZ1)); 
426           Double_t eta2 = 0.5*log((p2 + pZ2)/(p2 - pZ2)); 
427           Double_t deta = TMath::Abs(eta1 - eta2);
428           ibin = Int_t(deta/fP2Step);
429           
430           if((charge1 > 0.)&&(charge2 > 0.)) fNpp[ibin] += 1.;
431           if((charge1 < 0.)&&(charge2 < 0.)) fNnn[ibin] += 1.;
432           if((charge1 > 0.)&&(charge2 < 0.)) fNpn[ibin] += 1.;
433           if((charge1 < 0.)&&(charge2 > 0.)) fNpn[ibin] += 1.;
434           //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]);      
435         }
436         else continue;
437       }
438     }
439   }//case 1
440   if(fAnalysisType==2) {
441     for(i = 1; i < gNtrack; i++) {
442       if(fAnalysisLevel == "ESD")
443         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
444       else if(fAnalysisLevel == "AOD")
445         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
446       else if(fAnalysisLevel == "MC")
447         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
448
449       Short_t charge1 = 0;
450       Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
451       if(track1) {
452         charge1 = track1->Charge();
453         pX1 = track1->Px();
454         pY1 = track1->Py();
455         pZ1 = track1->Pz();
456         energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
457                                        TMath::Power(track1->M(),2));
458       }
459       else continue;
460
461       for(j = 0; j < i; j++) {
462         if(fAnalysisLevel == "ESD")
463           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
464         else if(fAnalysisLevel == "AOD")
465           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
466         else if(fAnalysisLevel == "MC")
467           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
468
469         if(track2) {
470           Short_t charge2 = track2->Charge();
471           Double_t pX2 = track2->Px();
472           Double_t pY2 = track2->Py();
473           Double_t pZ2 = track2->Pz();
474           Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
475                                          TMath::Power(track2->M(),2));
476           Double_t eTot = energy1 + energy2;
477           Double_t pxTot = pX1 + pX2;
478           Double_t pyTot = pY1 + pY2;
479           Double_t pzTot = pZ1 + pZ2;
480           Double_t q0Tot = energy1 - energy2;
481           Double_t qzTot = pZ1 - pZ2;
482           Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
483           Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
484           Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
485           ibin = Int_t(qLong/fP2Step);
486           //cout<<ibin<<endl;
487           if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
488           if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
489           if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
490           if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
491         }
492         else continue;
493       }
494     }
495   }//case 2
496   if(fAnalysisType==3) {
497     for(i = 1; i < gNtrack; i++) {
498       if(fAnalysisLevel == "ESD")
499         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
500       else if(fAnalysisLevel == "AOD")
501         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
502       else if(fAnalysisLevel == "MC")
503         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
504
505       Short_t charge1 = 0;
506       Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
507       if(track1) {
508         charge1 = track1->Charge();
509         pX1 = track1->Px();
510         pY1 = track1->Py();
511         pZ1 = track1->Pz();
512         energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
513                               TMath::Power(track1->M(),2));
514       }
515       else continue;
516
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
525         if(track2) {
526           Short_t charge2 = track2->Charge();
527           Double_t pX2 = track2->Px();
528           Double_t pY2 = track2->Py();
529           Double_t pZ2 = track2->Pz();
530           Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
531                                          TMath::Power(track2->M(),2));
532           Double_t eTot = energy1 + energy2;
533           Double_t pxTot = pX1 + pX2;
534           Double_t pyTot = pY1 + pY2;
535           Double_t pzTot = pZ1 + pZ2;
536           Double_t qxTot = pX1 - pX2;
537           Double_t qyTot = pY1 - pY2;
538           Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
539           Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
540           Double_t qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
541           ibin = Int_t(qOut/fP2Step);
542           //cout<<ibin<<endl;
543           if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
544           if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
545           if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
546           if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
547         }
548         else continue;
549       }
550     }
551   }//case 3
552   if(fAnalysisType==4) {
553     for(i = 1; i < gNtrack; i++) {
554       if(fAnalysisLevel == "ESD")
555         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
556       else if(fAnalysisLevel == "AOD")
557         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
558       else if(fAnalysisLevel == "MC")
559         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
560
561       Short_t charge1 = 0;
562       Double_t pX1 = 0., pY1 = 0.;
563       if(track1) {
564         charge1 = track1->Charge();
565         pX1 = track1->Px();
566         pY1 = track1->Py();
567         //Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
568         //TMath::Power(track1->M(),2));
569       }
570       else continue;
571
572       for(j = 0; j < i; j++) {
573         if(fAnalysisLevel == "ESD")
574           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
575         else if(fAnalysisLevel == "AOD")
576           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
577         else if(fAnalysisLevel == "MC")
578           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
579
580         if(track2) {
581           Short_t charge2 = track2->Charge();
582           Double_t pX2 = track2->Px();
583           Double_t pY2 = track2->Py();
584           //Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
585           //TMath::Power(track2->M(),2));
586           //Double_t eTot = energy1 + energy2;
587           Double_t pxTot = pX1 + pX2;
588           Double_t pyTot = pY1 + pY2;
589           Double_t qxTot = pX1 - pX2;
590           Double_t qyTot = pY1 - pY2;
591           Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
592           Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
593           ibin = Int_t(qSide/fP2Step);
594           //cout<<ibin<<endl;
595           if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
596           if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
597           if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
598           if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
599         }
600         else continue;
601       }
602     }
603   }//case 4
604   if(fAnalysisType==5) {
605     for(i = 1; i < gNtrack; i++) {
606       if(fAnalysisLevel == "ESD")
607         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
608       else if(fAnalysisLevel == "AOD")
609         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
610       else if(fAnalysisLevel == "MC")
611         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
612
613       Short_t charge1 = 0;
614       Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
615       if(track1) {
616         charge1 = track1->Charge();
617         pX1 = track1->Px();
618         pY1 = track1->Py();
619         pZ1 = track1->Pz();
620         energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
621                               TMath::Power(track1->M(),2));
622       }
623       else continue;
624
625       for(j = 0; j < i; j++) {
626         if(fAnalysisLevel == "ESD")
627           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
628         else if(fAnalysisLevel == "AOD")
629           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
630         else if(fAnalysisLevel == "MC")
631           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
632
633         if(track2) {
634           Short_t charge2 = track2->Charge();
635           Double_t pX2 = track2->Px();
636           Double_t pY2 = track2->Py();
637           Double_t pZ2 = track2->Pz();
638           Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
639                                          TMath::Power(track2->M(),2));
640           Double_t q0Tot = energy1 - energy2;
641           Double_t qxTot = pX1 - pX2;
642           Double_t qyTot = pY1 - pY2;
643           Double_t qzTot = pZ1 - pZ2;
644           Double_t qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
645           ibin = Int_t(qInv/fP2Step);
646           //cout<<ibin<<endl;
647           if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
648           if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
649           if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
650           if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
651         }
652         else continue;
653       }
654     }
655   }//case 5     
656   if(fAnalysisType==6) {
657     for(i = 1; i < gNtrack; i++) {
658       if(fAnalysisLevel == "ESD")
659         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
660       else if(fAnalysisLevel == "AOD")
661         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
662       else if(fAnalysisLevel == "MC")
663         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
664
665       Short_t charge1 = 0;
666       Double_t phi1 = 0.;
667       if(track1) {
668         charge1 = track1->Charge();
669         phi1 = track1->Phi();
670       }
671       else continue;
672
673       for(j = 0; j < i; j++) {
674         if(fAnalysisLevel == "ESD")
675           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
676         else if(fAnalysisLevel == "AOD")
677           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
678         else if(fAnalysisLevel == "MC")
679           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
680
681         if(track2) {
682           Short_t charge2 = track2->Charge();
683           Double_t phi2 = track2->Phi();
684           Double_t dphi = TMath::Abs(phi1 - phi2);
685           ibin = Int_t(dphi/fP2Step);
686           if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
687           if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
688           if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
689           if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
690         }
691         else continue;
692       }
693     }
694   }//case 6
695
696   /*for(Int_t i = 0; i < fNumberOfBins; i++) 
697     Printf("bin: %d - Npp: %lf - Nnn: %lf - Nnp: %lf - Npn: %lf",i,fNpp[i],fNnn[i],fNpn[i],fNpn[i]);*/
698 }
699
700 //____________________________________________________________________//
701 TH1F *AliBalance::GetHistNnn() {
702   //Return the control histogram of the -- component
703   for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
704     fHistfNnn->SetBinContent(iBin,GetNnn(iBin-1));
705
706   return fHistfNnn;
707 }
708
709 //____________________________________________________________________//
710 TH1F *AliBalance::GetHistNpp() {
711   //Return the control histogram of the ++ component
712   for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
713     fHistfNpp->SetBinContent(iBin,GetNpp(iBin-1));
714
715   return fHistfNpp;
716 }
717
718 //____________________________________________________________________//
719 TH1F *AliBalance::GetHistNpn() {
720   //Return the control histogram of the +- component
721   for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
722     fHistfNpn->SetBinContent(iBin,GetNpn(iBin-1));
723
724   return fHistfNpn;
725 }
726
727 //____________________________________________________________________//
728 Double_t AliBalance::GetBalance(Int_t p2) {
729   // Returns the value of the balance function in bin p2
730   fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
731   
732   return fB[p2];
733 }
734     
735 //____________________________________________________________________//
736 Double_t AliBalance::GetError(Int_t p2) {
737   // Returns the error on the BF value for bin p2
738   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;
739
740   return ferror[p2];
741 }
742
743 //____________________________________________________________________//
744 void AliBalance::PrintResults() {
745   // Prints the results
746   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
747   Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
748   Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
749   Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
750   Double_t deltaBalP2 = 0.0, integral = 0.0;
751   Double_t deltaErrorNew = 0.0;
752
753   cout<<"=================================================="<<endl;
754   for(Int_t i = 0; i < fNumberOfBins; i++) { 
755     //x[i] = fP2Start + fP2Step*i + fP2Step/2;
756     x[i] = fP2Step*i + fP2Step/2;
757     cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
758   } 
759   cout<<"=================================================="<<endl;
760   for(Int_t i = 1; i < fNumberOfBins; i++) {
761     fSumXi += x[i];
762     fSumBi += fB[i];
763     fSumBiXi += fB[i]*x[i];
764     fSumBiXi2 += fB[i]*TMath::Power(x[i],2);
765     fSumBi2Xi2 += TMath::Power(fB[i],2)*TMath::Power(x[i],2);
766     fSumDeltaBi2 +=  TMath::Power(ferror[i],2);
767     fSumXi2DeltaBi2 += TMath::Power(x[i],2) * TMath::Power(ferror[i],2);
768     
769     deltaBalP2 += fP2Step*TMath::Power(ferror[i],2);
770     integral += fP2Step*fB[i];
771   }
772   for(Int_t i = 1; i < fNumberOfBins; i++) deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/TMath::Power(fSumBi,2);
773    
774   Double_t integralError = TMath::Sqrt(deltaBalP2);
775   
776   Double_t delta = fSumBiXi / fSumBi;
777   Double_t deltaError = (fSumBiXi / fSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + TMath::Power((fSumDeltaBi2/fSumBi),2) );
778  
779   cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
780   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
781   cout<<"New error: "<<deltaErrorNew<<endl;
782   cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
783   cout<<"=================================================="<<endl;
784 }
785   
786 //____________________________________________________________________//
787 TGraphErrors *AliBalance::DrawBalance() {
788   // Draws the BF
789   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
790   Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
791   Double_t b[MAXIMUM_NUMBER_OF_STEPS];
792   Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
793
794   if((fNp == 0)||(fNn == 0)) {
795     cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
796     cout<<"Aborting....."<<endl;
797     abort();
798   }
799   
800   for(Int_t i = 0; i < fNumberOfBins; i++) {
801     b[i] = GetBalance(i);
802     ber[i] = GetError(i);
803     //x[i] = fP2Start + fP2Step*i + fP2Step/2;
804     x[i] = fP2Step*i + fP2Step/2;
805     xer[i] = 0.0;
806   }
807   
808   TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber);
809   gr->SetMarkerStyle(25);
810   gr->GetXaxis()->SetTitleColor(1);
811   if(fAnalysisType==0) {
812     gr->GetXaxis()->SetTitle("#Delta y");
813     gr->GetYaxis()->SetTitle("B(#Delta y)");
814   }
815   if(fAnalysisType==1) {
816     gr->GetXaxis()->SetTitle("#Delta #eta");
817     gr->GetYaxis()->SetTitle("B(#Delta #eta)");
818   }
819   if(fAnalysisType==2) {
820     gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
821     gr->GetYaxis()->SetTitle("B(q_{long}) [(GeV/c)^{-1}]");
822   }
823   if(fAnalysisType==3) {
824     gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
825     gr->GetYaxis()->SetTitle("B(q_{out}) [(GeV/c)^{-1}]");
826   }
827   if(fAnalysisType==4) {
828     gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
829     gr->GetYaxis()->SetTitle("B(q_{side}) [(GeV/c)^{-1}]");
830   }
831   if(fAnalysisType==5) {
832     gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
833     gr->GetYaxis()->SetTitle("B(q_{inv}) [(GeV/c)^{-1}]");
834   }
835   if(fAnalysisType==6) {
836     gr->GetXaxis()->SetTitle("#Delta #phi");
837     gr->GetYaxis()->SetTitle("B(#Delta #phi)");
838   }
839
840   return gr;
841 }
842
843 //____________________________________________________________________//
844 void AliBalance::Merge(AliBalance *b) {
845   //Merging function to be used for proof and grid
846   fNp += b->GetNp();
847   fNn += b->GetNn();
848   for(Int_t i = 0; i < fNumberOfBins; i++) {
849     fNnn[i] += b->GetNnn(i);
850     fNpp[i] += b->GetNpp(i);
851     fNpn[i] += b->GetNpn(i);
852   }
853 }