7508a0dd570a340f4005762c76f36e8af7cdf141
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliBalancePsi.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: AliBalancePsi.cxx 54125 2012-01-24 21:07:41Z miweber $ */
15
16 //-----------------------------------------------------------------
17 //           Balance Function class
18 //   This is the class to deal with the Balance Function wrt Psi analysis
19 //   Origin: Panos Christakoglou, Nikhef, Panos.Christakoglou@cern.ch
20 //-----------------------------------------------------------------
21
22
23 //ROOT
24 #include <Riostream.h>
25 #include <TMath.h>
26 #include <TAxis.h>
27 #include <TH2D.h>
28 #include <TH3D.h>
29 #include <TLorentzVector.h>
30 #include <TObjArray.h>
31 #include <TGraphErrors.h>
32 #include <TString.h>
33
34 #include "AliVParticle.h"
35 #include "AliMCParticle.h"
36 #include "AliESDtrack.h"
37 #include "AliAODTrack.h"
38 #include "AliTHn.h"
39
40 #include "AliBalancePsi.h"
41
42 ClassImp(AliBalancePsi)
43
44 //____________________________________________________________________//
45 AliBalancePsi::AliBalancePsi() :
46   TObject(), 
47   bShuffle(kFALSE),
48   fAnalysisLevel("ESD"),
49   fAnalyzedEvents(0) ,
50   fCentralityId(0) ,
51   fCentStart(0.),
52   fCentStop(0.),
53   fHistP(0),
54   fHistN(0),
55   fHistPN(0),
56   fHistNP(0),
57   fHistPP(0),
58   fHistNN(0),
59   fPsiInterval(15.) {
60   // Default constructor
61 }
62
63 //____________________________________________________________________//
64 AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
65   TObject(balance), bShuffle(balance.bShuffle), 
66   fAnalysisLevel(balance.fAnalysisLevel),
67   fAnalyzedEvents(balance.fAnalyzedEvents), 
68   fCentralityId(balance.fCentralityId),
69   fCentStart(balance.fCentStart),
70   fCentStop(balance.fCentStop),
71   fHistP(balance.fHistP),
72   fHistN(balance.fHistN),
73   fHistPN(balance.fHistPN),
74   fHistNP(balance.fHistNP),
75   fHistPP(balance.fHistPP),
76   fHistNN(balance.fHistNN),
77   fPsiInterval(balance.fPsiInterval) {
78   //copy constructor
79 }
80
81 //____________________________________________________________________//
82 AliBalancePsi::~AliBalancePsi() {
83   // Destructor
84   delete fHistP;
85   delete fHistN;
86   delete fHistPN;
87   delete fHistNP;
88   delete fHistPP;
89   delete fHistNN;
90 }
91
92 //____________________________________________________________________//
93 void AliBalancePsi::InitHistograms() {
94   // single particle histograms
95   Int_t anaSteps   = 1;       // analysis steps
96   Int_t iBinSingle[nTrackVariablesSingle];        // binning for track variables
97   Double_t* dBinsSingle[nTrackVariablesSingle];   // bins for track variables  
98   TString axisTitleSingle[nTrackVariablesSingle]; // axis titles for track variables
99   
100   // two particle histograms
101   Int_t iBinPair[nTrackVariablesPair];         // binning for track variables
102   Double_t* dBinsPair[nTrackVariablesPair];    // bins for track variables  
103   TString axisTitlePair[nTrackVariablesPair];  // axis titles for track variables
104
105   //centrality
106   const Int_t kNCentralityBins = 9;
107   Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
108   //const Int_t kNCentralityBins = 200;
109   //Double_t centralityBins[kNCentralityBins+1];
110   //for(Int_t i = 0; i < kNCentralityBins+1; i++)
111   //centralityBins[i] = 0.0 + i * 0.5;
112   iBinSingle[0]       = kNCentralityBins;
113   dBinsSingle[0]      = centralityBins;
114   axisTitleSingle[0]  = "Centrality percentile [%]"; 
115   iBinPair[0]       = kNCentralityBins;
116   dBinsPair[0]      = centralityBins;
117   axisTitlePair[0]  = "Centrality percentile [%]"; 
118
119   //Psi_2
120   const Int_t kNPsi2Bins = 24;
121   Double_t psi2Bins[kNPsi2Bins+1];
122   for(Int_t i = 0; i < kNPsi2Bins+1; i++)
123     psi2Bins[i] = 0.0 + i * 7.5;
124   iBinSingle[1]       = kNPsi2Bins;
125   dBinsSingle[1]      = psi2Bins;
126   axisTitleSingle[1]  = "#phi - #Psi_{2} (#circ)";
127   iBinPair[1]       = kNPsi2Bins;
128   dBinsPair[1]      = psi2Bins;
129   axisTitlePair[1]  = "#phi - #Psi_{2} (#circ)"; 
130   
131   // eta
132   const Int_t kNEtaBins = 10;
133   Double_t etaBins[kNEtaBins+1];
134   for(Int_t i = 0; i < kNEtaBins+1; i++)
135     etaBins[i] = -1.0 + i * 0.2;
136   iBinSingle[2]       = kNEtaBins;
137   dBinsSingle[2]      = etaBins;
138   axisTitleSingle[2]  = "#eta"; 
139   
140   //#eta of triggered particle
141   /*iBinPair[2]       = kNEtaBins;
142   dBinsPair[2]      = etaBins;
143   axisTitlePair[2]  = "#eta"; */
144
145   // delta eta
146   const Int_t kNDeltaEtaBins = 80;
147   Double_t deltaEtaBins[kNDeltaEtaBins+1];
148   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
149     deltaEtaBins[i] = -2.0 + i * 0.05;
150   iBinPair[2]       = kNDeltaEtaBins;
151   dBinsPair[2]      = deltaEtaBins;
152   axisTitlePair[2]  = "#Delta #eta"; 
153
154   // y
155   /*const Int_t kNYBins = 40;
156   Double_t yBins[kNYBins+1];
157   for(Int_t i = 0; i < kNYBins+1; i++)
158     yBins[i] = -1.0 + i * 0.05;
159   iBinSingle[3]       = kNYBins;
160   dBinsSingle[3]      = yBins;
161   axisTitleSingle[3]  = "y"; 
162   
163   //y of triggered particle
164   iBinPair[3]       = kNYBins;
165   dBinsPair[3]      = yBins;
166   axisTitlePair[3]  = "y"; */
167
168   // delta y
169   /*const Int_t kNDeltaYBins = 40;
170   Double_t deltaYBins[kNDeltaYBins+1];
171   for(Int_t i = 0; i < kNDeltaYBins+1; i++)
172     deltaYBins[i] = -2.0 + i * 0.1;
173   iBinPair[8]       = kNDeltaYBins;
174   dBinsPair[8]      = deltaYBins;
175   axisTitlePair[8]  = "#Delta y"; */
176
177   // phi
178   const Int_t kNPhiBins = 18;
179   Double_t phiBins[kNPhiBins+1];
180   for(Int_t i = 0; i < kNPhiBins+1; i++){
181     phiBins[i] = 0.0 + i * 20.;
182   } 
183   iBinSingle[3]       = kNPhiBins;
184   dBinsSingle[3]      = phiBins;
185   axisTitleSingle[3]  = "#phi (#circ)"; 
186   
187   // delta phi
188   const Int_t kNDeltaPhiBins = 72;
189   Double_t deltaPhiBins[kNDeltaPhiBins+1];
190   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
191     deltaPhiBins[i] = -180.0 + i * 5.;
192   } 
193   iBinPair[3]       = kNDeltaPhiBins;
194   dBinsPair[3]      = deltaPhiBins;
195   axisTitlePair[3]  = "#Delta #phi (#circ)"; 
196
197   // pt(trigger-associated)
198   const Int_t kNPtBins = 50;
199   Double_t ptBins[kNPtBins+1];
200   for(Int_t i = 0; i < kNPtBins+1; i++){
201     ptBins[i] = 0.0 + i * 0.5;
202    } 
203   iBinSingle[4]       = kNPtBins;
204   dBinsSingle[4]      = ptBins;
205   axisTitleSingle[4]  = "p_{t}^{trig.} (GeV/c)"; 
206
207   iBinPair[4]       = kNPtBins;
208   dBinsPair[4]      = ptBins;
209   axisTitlePair[4]  = "p_{t}^{trig.} (GeV/c)"; 
210
211   iBinPair[5]       = kNPtBins;
212   dBinsPair[5]      = ptBins;
213   axisTitlePair[5]  = "p_{t}^{assoc.} (GeV/c)";   
214
215   // qside
216   /*const Int_t kNQSideBins = 200;
217   Double_t qSideBins[kNQSideBins+1];
218   for(Int_t i = 0; i < kNQSideBins+1; i++)
219     qSideBins[i] = 0.0 + i * 0.02;
220   iBinPair[10]       = kNQSideBins;
221   dBinsPair[10]      = qSideBins;
222   axisTitlePair[10]  = "q_{side} (GeV/c)";
223
224   // qout
225   const Int_t kNQoutBins = 200;
226   Double_t qoutBins[kNQoutBins+1];
227   for(Int_t i = 0; i < kNQoutBins+1; i++)
228     qoutBins[i] = 0.0 + i * 0.02;
229   iBinPair[11]       = kNQoutBins;
230   dBinsPair[11]      = qoutBins;
231   axisTitlePair[11]  = "q_{out} (GeV/c)";
232
233   // qlong
234   const Int_t kNQlongBins = 200;
235   Double_t qlongBins[kNQlongBins+1];
236   for(Int_t i = 0; i < kNQlongBins+1; i++)
237     qlongBins[i] = 0.0 + i * 0.02;
238   iBinPair[12]       = kNQlongBins;
239   dBinsPair[12]      = qlongBins;
240   axisTitlePair[12]  = "q_{long} (GeV/c)";
241
242   // qinv
243   const Int_t kNQinvBins = 200;
244   Double_t qinvBins[kNQinvBins+1];
245   for(Int_t i = 0; i < kNQinvBins+1; i++)
246     qinvBins[i] = 0.0 + i * 0.02;
247   iBinPair[13]       = kNQinvBins;
248   dBinsPair[13]      = qinvBins;
249   axisTitlePair[13]  = "q_{inv} (GeV/c)";*/
250
251   TString histName;
252   //+ triggered particles
253   histName = "fHistP"; 
254   if(bShuffle) histName.Append("_shuffle");
255   if(fCentralityId) histName += fCentralityId.Data();
256   fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,nTrackVariablesSingle,iBinSingle);
257   for (Int_t j=0; j<nTrackVariablesSingle; j++) {
258     fHistP->SetBinLimits(j, dBinsSingle[j]);
259     fHistP->SetVarTitle(j, axisTitleSingle[j]);
260   }
261
262   //- triggered particles
263   histName = "fHistN"; 
264   if(bShuffle) histName.Append("_shuffle");
265   if(fCentralityId) histName += fCentralityId.Data();
266   fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,nTrackVariablesSingle,iBinSingle);
267   for (Int_t j=0; j<nTrackVariablesSingle; j++) {
268     fHistN->SetBinLimits(j, dBinsSingle[j]);
269     fHistN->SetVarTitle(j, axisTitleSingle[j]);
270   }
271   
272   //+- pairs
273   histName = "fHistPN";
274   if(bShuffle) histName.Append("_shuffle");
275   if(fCentralityId) histName += fCentralityId.Data();
276   fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
277   for (Int_t j=0; j<nTrackVariablesPair; j++) {
278     fHistPN->SetBinLimits(j, dBinsPair[j]);
279     fHistPN->SetVarTitle(j, axisTitlePair[j]);
280   }
281
282   //-+ pairs
283   histName = "fHistNP";
284   if(bShuffle) histName.Append("_shuffle");
285   if(fCentralityId) histName += fCentralityId.Data();
286   fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
287   for (Int_t j=0; j<nTrackVariablesPair; j++) {
288     fHistNP->SetBinLimits(j, dBinsPair[j]);
289     fHistNP->SetVarTitle(j, axisTitlePair[j]);
290   }
291
292   //++ pairs
293   histName = "fHistPP";
294   if(bShuffle) histName.Append("_shuffle");
295   if(fCentralityId) histName += fCentralityId.Data();
296   fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
297   for (Int_t j=0; j<nTrackVariablesPair; j++) {
298     fHistPP->SetBinLimits(j, dBinsPair[j]);
299     fHistPP->SetVarTitle(j, axisTitlePair[j]);
300   }
301
302   //-- pairs
303   histName = "fHistNN";
304   if(bShuffle) histName.Append("_shuffle");
305   if(fCentralityId) histName += fCentralityId.Data();
306   fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
307   for (Int_t j=0; j<nTrackVariablesPair; j++) {
308     fHistNN->SetBinLimits(j, dBinsPair[j]);
309     fHistNN->SetVarTitle(j, axisTitlePair[j]);
310   }
311   Printf("Finished setting up the AliTHn");
312 }
313
314 //____________________________________________________________________//
315 void AliBalancePsi::CalculateBalance(Float_t fCentrality, 
316                                      Double_t gReactionPlane, 
317                                      vector<Double_t> **chargeVector) {
318   // Calculates the balance function
319   fAnalyzedEvents++;
320   Int_t i = 0 , j = 0;
321   
322   // Initialize histograms if not done yet
323   if(!fHistPN){
324     AliWarning("Histograms not yet initialized! --> Will be done now");
325     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
326     InitHistograms();
327   }
328
329   Double_t trackVarsSingle[nTrackVariablesSingle];
330   Double_t trackVarsPair[nTrackVariablesPair];
331
332   Int_t gNtrack = chargeVector[0]->size();
333   //Printf("(AliBalancePsi) Number of tracks: %d",gNtrack);
334   Double_t gPsiMinusPhi = 0.;
335   Double_t dy = 0., deta = 0.;
336   Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
337   Double_t dphi = 0.;
338
339   Double_t charge1  = 0;
340   Double_t eta1 = 0., rap1 = 0.;
341   Double_t px1 = 0., py1 = 0., pz1 = 0.;
342   Double_t pt1 = 0.;
343   Double_t energy1 = 0.;
344   Double_t phi1    = 0.;
345
346   Double_t charge2  = 0;
347   Double_t eta2 = 0., rap2 = 0.;
348   Double_t px2 = 0., py2 = 0., pz2 = 0.;
349   Double_t pt2 = 0.;
350   Double_t energy2 = 0.;
351   Double_t phi2    = 0.;
352   
353   for(i = 0; i < gNtrack; i++) {
354     charge1 = chargeVector[0]->at(i);
355     rap1    = chargeVector[1]->at(i);
356     eta1    = chargeVector[2]->at(i);
357     phi1    = chargeVector[3]->at(i);
358     px1     = chargeVector[4]->at(i);
359     py1     = chargeVector[5]->at(i);
360     pz1     = chargeVector[6]->at(i);
361     pt1     = chargeVector[7]->at(i);
362     energy1 = chargeVector[8]->at(i);
363     gPsiMinusPhi   = TMath::Abs(phi1 - gReactionPlane);
364     if(gPsiMinusPhi > 180.) gPsiMinusPhi = 360. - gPsiMinusPhi;
365     //if(gPsiMinusPhi < -fPsiInterval/2) gPsiMinusPhi = 360. + gPsiMinusPhi;
366     
367     trackVarsSingle[0]    =  fCentrality;
368     trackVarsSingle[1]    =  gPsiMinusPhi;
369     trackVarsSingle[2]    =  eta1;
370     //trackVarsSingle[3]    =  rap1;
371     trackVarsSingle[3]    =  phi1;  
372     trackVarsSingle[4]    =  pt1;  
373
374     //Printf("Track(a) %d - phi-Psi: %lf",i+1,trackVarsSingle[1]);
375     //fill single particle histograms
376     if(charge1 > 0)  
377       fHistP->Fill(trackVarsSingle,0,1.); 
378     else if(charge1 < 0)            
379       fHistN->Fill(trackVarsSingle,0,1.); 
380
381     for(j = 0; j < i; j++) {
382       charge2 = chargeVector[0]->at(j);
383       rap2    = chargeVector[1]->at(j);
384       eta2    = chargeVector[2]->at(j);
385       phi2    = chargeVector[3]->at(j);
386       px2     = chargeVector[4]->at(j);
387       py2     = chargeVector[5]->at(j);
388       pz2     = chargeVector[6]->at(j);
389       pt2     = chargeVector[7]->at(j);
390       energy2 = chargeVector[8]->at(j);
391       //Printf("Track(b) %d - pt: %lf",j+1,pt2);
392       
393       // filling the arrays
394       // RAPIDITY 
395       dy = rap1 - rap2;
396       
397       // Eta
398       deta = eta1 - eta2;
399       
400       //qlong
401       Double_t eTot = energy1 + energy2;
402       Double_t pxTot = px1 + px2;
403       Double_t pyTot = py1 + py2;
404       Double_t pzTot = pz1 + pz2;
405       Double_t q0Tot = energy1 - energy2;
406       Double_t qxTot = px1 - px2;
407       Double_t qyTot = py1 - py2;
408       Double_t qzTot = pz1 - pz2;
409       
410       Double_t eTot2 = eTot*eTot;
411       Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
412       Double_t pzTot2 = pzTot*pzTot;
413       
414       Double_t q0Tot2 = q0Tot*q0Tot;
415       Double_t qTot2  = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;
416       
417       Double_t snn    = eTot2 - pTot2;
418       Double_t ptTot2 = pTot2 - pzTot2 ;
419       Double_t ptTot  = TMath::Sqrt( ptTot2 );
420       
421       qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
422       
423       //qout
424       qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
425       
426       //qside
427       qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
428       
429       //qinv
430       qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
431       
432       //phi
433       dphi = phi2 - phi1;
434       if(dphi < -180.) dphi += 360.;  //dphi should be between -180 and 180!
435       else if(dphi > 180.) dphi -= 360.; //dphi should be between -180 and 180!
436       
437       trackVarsPair[0]    =  fCentrality;             
438       trackVarsPair[1]    =  gPsiMinusPhi;             
439       //trackVarsPair[2]    =  eta1;
440       //trackVarsPair[3]    =  rap1;
441       //trackVarsPair[4]    =  phi1;
442       trackVarsPair[2]    =  deta;
443       //trackVarsPair[8]    =  dy;
444       trackVarsPair[3]    =  dphi;
445       trackVarsPair[4]    =  pt1;
446       trackVarsPair[5]    =  pt2;
447       //trackVarsPair[10]   =  qSide;
448       //trackVarsPair[11]   =  qOut;
449       //trackVarsPair[12]   =  qLong;
450       //trackVarsPair[13]   =  qInv;
451       
452       if( charge1 > 0 && charge2 < 0)  fHistPN->Fill(trackVarsPair,0,1.); 
453       else if( charge1 < 0 && charge2 > 0)  fHistNP->Fill(trackVarsPair,0,1.); 
454       else if( charge1 > 0 && charge2 > 0)  fHistPP->Fill(trackVarsPair,0,1.); 
455       else if( charge1 < 0 && charge2 < 0)  fHistNN->Fill(trackVarsPair,0,1.); 
456       else AliWarning("Wrong charge combination!");
457       
458     }//end of 2nd particle loop
459   }//end of 1st particle loop
460 }  
461
462 //____________________________________________________________________//
463 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
464                                                  Int_t iVariablePair,
465                                                  Double_t centrMin, 
466                                                  Double_t centrMax, 
467                                                  Double_t psiMin, 
468                                                  Double_t psiMax) {
469   //Returns the BF histogram, extracted from the 6 AliTHn objects 
470   //(private members) of the AliBalancePsi class.
471   //iVariableSingle: 2(eta), 3(y), 4(phi) 
472   //iVariablePair: 2(Delta eta), 3(Delta y), 4(Delta phi), 5(qside), 6(qout), 7(qlong) 8(qinv) 
473   TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","phi","qlong","qout","qside","qinv"};
474   TString histName = "gHistBalanceFunctionHistogram";
475   histName += gAnalysisType[iVariablePair];
476
477   // centrality
478   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
479   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
480   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
481   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
482   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
483   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
484
485   // Psi_2
486   fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
487   fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
488   fHistPN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
489   fHistNP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
490   fHistPP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
491   fHistNN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
492
493   //Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0));
494
495   // Project into the wanted space (1st: analysis step, 2nd: axis)
496   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
497   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
498   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
499   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
500   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
501   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
502
503   TH1D *gHistBalanceFunctionHistogram = 0x0;
504   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
505     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
506     gHistBalanceFunctionHistogram->Reset();
507     
508     switch(iVariablePair) {
509     case 2:
510       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
511       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
512       break;
513     case 3:
514       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
515       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
516       break;
517     case 4:
518       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
519       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
520       break;
521     case 5:
522       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
523       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
524       break;
525     case 6:
526       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
527       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
528       break;
529     case 7:
530       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
531       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
532       break;
533     case 8:
534       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
535       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
536       break;
537     default:
538       break;
539     }
540
541     hTemp1->Sumw2();
542     hTemp2->Sumw2();
543     hTemp3->Sumw2();
544     hTemp4->Sumw2();
545     hTemp1->Add(hTemp3,-1.);
546     hTemp1->Scale(1./hTemp5->GetEntries());
547     hTemp2->Add(hTemp4,-1.);
548     hTemp2->Scale(1./hTemp6->GetEntries());
549     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
550     gHistBalanceFunctionHistogram->Scale(0.5);
551   }
552
553   return gHistBalanceFunctionHistogram;
554 }
555
556 //____________________________________________________________________//
557 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t centrMin, 
558                                               Double_t centrMax, 
559                                               Double_t psiMin, 
560                                               Double_t psiMax) {
561   //Returns the 2D correlation function for (+-) pairs
562   // centrality: axis 0
563   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
564   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
565   
566   // Psi_2: axis 1
567   fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
568   fHistPN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
569     
570   //0:step, 2: Delta eta, 3: Delta phi
571   TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,2,3));
572   if(!gHist){
573     AliError("Projection of fHistPN = NULL");
574     return gHist;
575   }
576
577   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
578   //Printf("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,2,3)->GetEntries()));
579   if((Double_t)(fHistP->Project(0,2)->GetEntries())!=0)
580     gHist->Scale(1./(Double_t)(fHistP->Project(0,2)->GetEntries()));
581     
582   return gHist;
583 }
584
585 //____________________________________________________________________//
586 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t centrMin, 
587                                               Double_t centrMax, 
588                                               Double_t psiMin, 
589                                               Double_t psiMax) {
590   //Returns the 2D correlation function for (+-) pairs
591   // centrality: axis 0
592   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
593   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
594   
595   // Psi_2: axis 1
596   fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
597   fHistNP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
598     
599   //0:step, 2: Delta eta, 3: Delta phi
600   TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,2,3));
601   if(!gHist){
602     AliError("Projection of fHistPN = NULL");
603     return gHist;
604   }
605
606   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
607   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
608   if((Double_t)(fHistN->Project(0,2)->GetEntries())!=0)
609     gHist->Scale(1./(Double_t)(fHistN->Project(0,2)->GetEntries()));
610     
611   return gHist;
612 }
613
614 //____________________________________________________________________//
615 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t centrMin, 
616                                               Double_t centrMax, 
617                                               Double_t psiMin, 
618                                               Double_t psiMax) {
619   //Returns the 2D correlation function for (+-) pairs
620   // centrality: axis 0
621   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
622   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
623   
624   // Psi_2: axis 1
625   fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
626   fHistPP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
627       
628   //0:step, 2: Delta eta, 3: Delta phi
629   TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,2,3));
630   if(!gHist){
631     AliError("Projection of fHistPN = NULL");
632     return gHist;
633   }
634
635   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
636   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
637   if((Double_t)(fHistP->Project(0,2)->GetEntries())!=0)
638     gHist->Scale(1./(Double_t)(fHistP->Project(0,2)->GetEntries()));
639   
640   return gHist;
641 }
642
643 //____________________________________________________________________//
644 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t centrMin, 
645                                               Double_t centrMax, 
646                                               Double_t psiMin, 
647                                               Double_t psiMax) {
648   //Returns the 2D correlation function for (+-) pairs
649   // centrality: axis 0
650   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
651   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
652   
653   // Psi_2: axis 1
654   fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
655   fHistNN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
656     
657   //0:step, 2: Delta eta, 3: Delta phi
658   TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,2,3));
659   if(!gHist){
660     AliError("Projection of fHistPN = NULL");
661     return gHist;
662   }
663
664   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
665   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
666   if((Double_t)(fHistN->Project(0,2)->GetEntries())!=0)
667     gHist->Scale(1./(Double_t)(fHistN->Project(0,2)->GetEntries()));
668     
669   return gHist;
670 }