Revert change on Pair pT cut:
[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 <TCanvas.h>
26 #include <TMath.h>
27 #include <TAxis.h>
28 #include <TH2D.h>
29 #include <TH3D.h>
30 #include <TLorentzVector.h>
31 #include <TObjArray.h>
32 #include <TGraphErrors.h>
33 #include <TString.h>
34 #include <TSpline.h>
35 #include <TRandom3.h>
36
37 #include "AliVParticle.h"
38 #include "AliMCParticle.h"
39 #include "AliESDtrack.h"
40 #include "AliAODTrack.h"
41 #include "AliTHn.h"
42 #include "AliAnalysisTaskTriggeredBF.h"
43
44 #include "AliBalancePsi.h"
45 using std::cout;
46 using std::endl;
47 using std::cerr;
48
49 ClassImp(AliBalancePsi)
50
51 //____________________________________________________________________//
52 AliBalancePsi::AliBalancePsi() :
53   TObject(), 
54   fShuffle(kFALSE),
55   fAnalysisLevel("ESD"),
56   fAnalyzedEvents(0) ,
57   fCentralityId(0) ,
58   fCentStart(0.),
59   fCentStop(0.),
60   fHistP(0),
61   fHistN(0),
62   fHistPN(0),
63   fHistNP(0),
64   fHistPP(0),
65   fHistNN(0),
66   fHistHBTbefore(0),
67   fHistHBTafter(0),
68   fHistConversionbefore(0),
69   fHistConversionafter(0),
70   fHistPsiMinusPhi(0),
71   fHistResonancesBefore(0),
72   fHistResonancesRho(0),
73   fHistResonancesK0(0),
74   fHistResonancesLambda(0),
75   fHistQbefore(0),
76   fHistQafter(0),
77   fPsiInterval(15.),
78   fDeltaEtaMax(2.0),
79   fResonancesCut(kFALSE),
80   fHBTCut(kFALSE),
81   fHBTCutValue(0.02),
82   fConversionCut(kFALSE),
83   fInvMassCutConversion(0.04),
84   fQCut(kFALSE),
85   fDeltaPtMin(0.0),
86   fVertexBinning(kFALSE),
87   fCustomBinning(""),
88   fBinningString(""),
89   fEventClass("EventPlane"){
90   // Default constructor
91 }
92
93 //____________________________________________________________________//
94 AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
95   TObject(balance), fShuffle(balance.fShuffle), 
96   fAnalysisLevel(balance.fAnalysisLevel),
97   fAnalyzedEvents(balance.fAnalyzedEvents), 
98   fCentralityId(balance.fCentralityId),
99   fCentStart(balance.fCentStart),
100   fCentStop(balance.fCentStop),
101   fHistP(balance.fHistP),
102   fHistN(balance.fHistN),
103   fHistPN(balance.fHistPN),
104   fHistNP(balance.fHistNP),
105   fHistPP(balance.fHistPP),
106   fHistNN(balance.fHistNN),
107   fHistHBTbefore(balance.fHistHBTbefore),
108   fHistHBTafter(balance.fHistHBTafter),
109   fHistConversionbefore(balance.fHistConversionbefore),
110   fHistConversionafter(balance.fHistConversionafter),
111   fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
112   fHistResonancesBefore(balance.fHistResonancesBefore),
113   fHistResonancesRho(balance.fHistResonancesRho),
114   fHistResonancesK0(balance.fHistResonancesK0),
115   fHistResonancesLambda(balance.fHistResonancesLambda),
116   fHistQbefore(balance.fHistQbefore),
117   fHistQafter(balance.fHistQafter),
118   fPsiInterval(balance.fPsiInterval),
119   fDeltaEtaMax(balance.fDeltaEtaMax),
120   fResonancesCut(balance.fResonancesCut),
121   fHBTCut(balance.fHBTCut),
122   fHBTCutValue(balance.fHBTCutValue),
123   fConversionCut(balance.fConversionCut),
124   fInvMassCutConversion(balance.fInvMassCutConversion),
125   fQCut(balance.fQCut),
126   fDeltaPtMin(balance.fDeltaPtMin),
127   fVertexBinning(balance.fVertexBinning),
128   fCustomBinning(balance.fCustomBinning),
129   fBinningString(balance.fBinningString),
130   fEventClass("EventPlane"){
131   //copy constructor
132 }
133
134 //____________________________________________________________________//
135 AliBalancePsi::~AliBalancePsi() {
136   // Destructor
137   delete fHistP;
138   delete fHistN;
139   delete fHistPN;
140   delete fHistNP;
141   delete fHistPP;
142   delete fHistNN;
143
144   delete fHistHBTbefore;
145   delete fHistHBTafter;
146   delete fHistConversionbefore;
147   delete fHistConversionafter;
148   delete fHistPsiMinusPhi;
149   delete fHistResonancesBefore;
150   delete fHistResonancesRho;
151   delete fHistResonancesK0;
152   delete fHistResonancesLambda;
153   delete fHistQbefore;
154   delete fHistQafter;
155     
156 }
157
158 //____________________________________________________________________//
159 void AliBalancePsi::InitHistograms() {
160   // single particle histograms
161
162   // global switch disabling the reference 
163   // (to avoid "Replacing existing TH1" if several wagons are created in train)
164   Bool_t oldStatus = TH1::AddDirectoryStatus();
165   TH1::AddDirectory(kFALSE);
166
167   Int_t anaSteps   = 1;       // analysis steps
168   Int_t iBinSingle[kTrackVariablesSingle];        // binning for track variables
169   Double_t* dBinsSingle[kTrackVariablesSingle];   // bins for track variables  
170   TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
171   
172   // two particle histograms
173   Int_t iBinPair[kTrackVariablesPair];         // binning for track variables
174   Double_t* dBinsPair[kTrackVariablesPair];    // bins for track variables  
175   TString axisTitlePair[kTrackVariablesPair];  // axis titles for track variables
176
177
178
179   // =========================================================
180   // The default string (from older versions of AliBalancePsi)
181   // =========================================================
182   TString defaultBinningStr;
183   defaultBinningStr = "multiplicity: 0,10,20,30,40,50,60,70,80,100,100000\n"                // Multiplicity Bins
184     "centrality: 0.,5.,10.,20.,30.,40.,50.,60.,70.,80.\n"                                   // Centrality Bins
185     "centralityVertex: 0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.\n" // Centrality Bins (Vertex Binning)
186     "eventPlane: -0.5,0.5,1.5,2.5,3.5\n"                                                    // Event Plane Bins (Psi: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest))
187     "deltaEta: -1.6, -1.56, -1.52, -1.48, -1.44, -1.4, -1.36, -1.32, -1.28, -1.24, -1.2, -1.16, -1.12, -1.08, -1.04, -1, -0.96, -0.92, -0.88, -0.84, -0.8, -0.76, -0.72, -0.68, -0.64, -0.6, -0.56, -0.52, -0.48, -0.44, -0.4, -0.36, -0.32, -0.28, -0.24, -0.2, -0.16, -0.12, -0.08, -0.04, 0, 0.04, 0.08, 0.12, 0.16, 0.2, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48, 0.52, 0.56, 0.6, 0.64, 0.68, 0.72, 0.76, 0.8, 0.84, 0.88, 0.92, 0.96, 1, 1.04, 1.08, 1.12, 1.16, 1.2, 1.24, 1.28, 1.32, 1.36, 1.4, 1.44, 1.48, 1.52, 1.56, 1.6\n" // Delta Eta Bins
188     "deltaEtaVertex: -1.6, -1.52, -1.44, -1.36, -1.28, -1.2, -1.12, -1.04, -0.96, -0.88, -0.8, -0.72, -0.64, -0.56, -0.48, -0.4, -0.32, -0.24, -0.16, -0.08, 0, 0.08, 0.16, 0.24, 0.32, 0.4, 0.48, 0.56, 0.64, 0.72, 0.8, 0.88, 0.96, 1.04, 1.12, 1.2, 1.28, 1.36, 1.44, 1.52, 1.6\n" // Delta Eta Bins (Vertex Binning)
189     "deltaPhi: -1.5708, -1.48353, -1.39626, -1.309, -1.22173, -1.13446, -1.0472, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.0872665, 0, 0.0872665, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.0472, 1.13446, 1.22173, 1.309, 1.39626, 1.48353, 1.5708, 1.65806, 1.74533, 1.8326, 1.91986, 2.00713, 2.0944, 2.18166, 2.26893, 2.35619, 2.44346, 2.53073, 2.61799, 2.70526, 2.79253, 2.87979, 2.96706, 3.05433, 3.14159, 3.22886, 3.31613, 3.40339, 3.49066, 3.57792, 3.66519, 3.75246, 3.83972, 3.92699, 4.01426, 4.10152, 4.18879, 4.27606, 4.36332, 4.45059, 4.53786, 4.62512, 4.71239\n" // Delta Phi Bins
190     "pT: 0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.\n"             // pT Bins
191     "pTVertex: 0.2,1.0,2.0,3.0,4.0,8.0,15.0\n"                                              // pT Bins (Vertex Binning)
192     "vertex: -10., 10.\n"                                                                   // Vertex Bins
193     "vertexVertex: -10., -7., -5., -3., -1., 1., 3., 5., 7., 10.\n"                         // Vertex Bins (Vertex Binning)
194     ;
195   
196   
197   // =========================================================
198   // Customization (adopted from AliUEHistograms)
199   // =========================================================
200
201   TObjArray* lines = defaultBinningStr.Tokenize("\n");
202   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
203   {
204     TString line(lines->At(i)->GetName());
205     TString tag = line(0, line.Index(":")+1);
206     if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
207       fBinningString += line + "\n";
208     else
209       AliInfo(Form("Using custom binning for %s", tag.Data()));
210   }
211   delete lines;
212   fBinningString += fCustomBinning;
213   
214   AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
215
216
217   // =========================================================
218   // Now set the bins
219   // =========================================================
220
221   //Depending on fEventClass Variable, do one thing or the other...
222   /**********************************************************
223    
224   ======> Modification: Change Event Classification Scheme
225     
226   ---> fEventClass == "EventPlane"
227    
228    Default operation with Event Plane 
229    
230   ---> fEventClass == "Multiplicity"
231    
232    Work with reference multiplicity (from GetReferenceMultiplicity, which one is decided in the configuration)
233    
234   ---> fEventClass == "Centrality" 
235    
236    Work with Centrality Bins
237
238   ***********************************************************/
239   if(fEventClass == "Multiplicity"){
240     dBinsSingle[0]     = GetBinning(fBinningString, "multiplicity", iBinSingle[0]);
241     dBinsPair[0]       = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
242     axisTitleSingle[0] = "reference multiplicity";
243     axisTitlePair[0]   = "reference multiplicity";
244   }
245   if(fEventClass == "Centrality"){
246     // fine binning in case of vertex Z binning
247     if(fVertexBinning){
248       dBinsSingle[0]     = GetBinning(fBinningString, "centralityVertex", iBinSingle[0]);
249       dBinsPair[0]       = GetBinning(fBinningString, "centralityVertex", iBinPair[0]);
250     }
251     else{
252       dBinsSingle[0]     = GetBinning(fBinningString, "centrality", iBinSingle[0]);
253       dBinsPair[0]       = GetBinning(fBinningString, "centrality", iBinPair[0]);
254     }
255     axisTitleSingle[0] = "Centrality percentile [%]";
256     axisTitlePair[0]   = "Centrality percentile [%]";
257   }
258   if(fEventClass == "EventPlane"){
259     dBinsSingle[0]     = GetBinning(fBinningString, "eventPlane", iBinSingle[0]);
260     dBinsPair[0]       = GetBinning(fBinningString, "eventPlane", iBinPair[0]);
261     axisTitleSingle[0] = "#varphi - #Psi_{2} (a.u.)";
262     axisTitlePair[0]   = "#varphi - #Psi_{2} (a.u.)";
263   }
264   
265
266   // Delta Eta and Delta Phi
267   // (coarse binning in case of vertex Z binning)
268   if(fVertexBinning){
269     dBinsPair[1]       = GetBinning(fBinningString, "deltaEtaVertex", iBinPair[1]);
270   }
271   else{
272     dBinsPair[1]       = GetBinning(fBinningString, "deltaEta", iBinPair[1]);
273   }
274   axisTitlePair[1]  = "#Delta#eta"; 
275   
276   dBinsPair[2]       = GetBinning(fBinningString, "deltaPhi", iBinPair[2]);
277   axisTitlePair[2]   = "#Delta#varphi (rad)";  
278   
279
280   // pT Trig and pT Assoc
281   // (coarse binning in case of vertex Z binning)
282   if(fVertexBinning){
283     dBinsSingle[1]   = GetBinning(fBinningString, "pTVertex", iBinSingle[1]);
284     dBinsPair[3]     = GetBinning(fBinningString, "pTVertex", iBinPair[3]);
285     dBinsPair[4]     = GetBinning(fBinningString, "pTVertex", iBinPair[4]);
286   }
287   else{
288     dBinsSingle[1]   = GetBinning(fBinningString, "pT", iBinSingle[1]);
289     dBinsPair[3]     = GetBinning(fBinningString, "pT", iBinPair[3]);
290     dBinsPair[4]     = GetBinning(fBinningString, "pT", iBinPair[4]);
291   }
292   
293   axisTitleSingle[1]  = "p_{T,trig.} (GeV/c)"; 
294   axisTitlePair[3]    = "p_{T,trig.} (GeV/c)"; 
295   axisTitlePair[4]    = "p_{T,assoc.} (GeV/c)";  
296  
297
298   // vertex Z binning or not
299   if(fVertexBinning){
300     dBinsSingle[2]   = GetBinning(fBinningString, "vertexVertex", iBinSingle[2]);
301     dBinsPair[5]     = GetBinning(fBinningString, "vertexVertex", iBinPair[5]);
302   }
303   else{
304     dBinsSingle[2]   = GetBinning(fBinningString, "vertex", iBinSingle[2]);
305     dBinsPair[5]     = GetBinning(fBinningString, "vertex", iBinPair[5]);
306   }
307
308   axisTitleSingle[2]  = "v_{Z} (cm)"; 
309   axisTitlePair[5]    = "v_{Z} (cm)"; 
310
311
312
313   // =========================================================
314   // Create the Output objects (AliTHn)
315   // =========================================================
316
317   TString histName;
318   //+ triggered particles
319   histName = "fHistP"; 
320   if(fShuffle) histName.Append("_shuffle");
321   if(fCentralityId) histName += fCentralityId.Data();
322   fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
323   for (Int_t j=0; j<kTrackVariablesSingle; j++) {
324     fHistP->SetBinLimits(j, dBinsSingle[j]);
325     fHistP->SetVarTitle(j, axisTitleSingle[j]);
326   }
327
328   //- triggered particles
329   histName = "fHistN"; 
330   if(fShuffle) histName.Append("_shuffle");
331   if(fCentralityId) histName += fCentralityId.Data();
332   fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
333   for (Int_t j=0; j<kTrackVariablesSingle; j++) {
334     fHistN->SetBinLimits(j, dBinsSingle[j]);
335     fHistN->SetVarTitle(j, axisTitleSingle[j]);
336   }
337   
338   //+- pairs
339   histName = "fHistPN";
340   if(fShuffle) histName.Append("_shuffle");
341   if(fCentralityId) histName += fCentralityId.Data();
342   fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
343   for (Int_t j=0; j<kTrackVariablesPair; j++) {
344     fHistPN->SetBinLimits(j, dBinsPair[j]);
345     fHistPN->SetVarTitle(j, axisTitlePair[j]);
346   }
347
348   //-+ pairs
349   histName = "fHistNP";
350   if(fShuffle) histName.Append("_shuffle");
351   if(fCentralityId) histName += fCentralityId.Data();
352   fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
353   for (Int_t j=0; j<kTrackVariablesPair; j++) {
354     fHistNP->SetBinLimits(j, dBinsPair[j]);
355     fHistNP->SetVarTitle(j, axisTitlePair[j]);
356   }
357
358   //++ pairs
359   histName = "fHistPP";
360   if(fShuffle) histName.Append("_shuffle");
361   if(fCentralityId) histName += fCentralityId.Data();
362   fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
363   for (Int_t j=0; j<kTrackVariablesPair; j++) {
364     fHistPP->SetBinLimits(j, dBinsPair[j]);
365     fHistPP->SetVarTitle(j, axisTitlePair[j]);
366   }
367
368   //-- pairs
369   histName = "fHistNN";
370   if(fShuffle) histName.Append("_shuffle");
371   if(fCentralityId) histName += fCentralityId.Data();
372   fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
373   for (Int_t j=0; j<kTrackVariablesPair; j++) {
374     fHistNN->SetBinLimits(j, dBinsPair[j]);
375     fHistNN->SetVarTitle(j, axisTitlePair[j]);
376   }
377
378   AliInfo("Finished setting up the AliTHn");
379
380   // QA histograms
381   fHistHBTbefore        = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
382   fHistHBTafter         = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
383   fHistConversionbefore = new TH3D("fHistConversionbefore","before Conversion cut;#Delta#eta;#Delta#phi;M_{inv}^{2}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
384   fHistConversionafter  = new TH3D("fHistConversionafter","after Conversion cut;#Delta#eta;#Delta#phi;M_{inv}^{2}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
385   fHistPsiMinusPhi      = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
386   fHistResonancesBefore = new TH3D("fHistResonancesBefore","before resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
387   fHistResonancesRho    = new TH3D("fHistResonancesRho","after #rho resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
388   fHistResonancesK0     = new TH3D("fHistResonancesK0","after #rho, K0 resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
389   fHistResonancesLambda = new TH3D("fHistResonancesLambda","after #rho, K0, Lambda resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
390   fHistQbefore          = new TH3D("fHistQbefore","before momentum difference cut;#Delta#eta;#Delta#phi;|#Delta p_{T}| (GeV/c)",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
391   fHistQafter           = new TH3D("fHistQafter","after momentum difference cut;#Delta#eta;#Delta#phi;|#Delta p_{T}| (GeV/c)",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
392
393   TH1::AddDirectory(oldStatus);
394
395 }
396
397 //____________________________________________________________________//
398 void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
399                                      TObjArray *particles, 
400                                      TObjArray *particlesMixed,
401                                      Float_t bSign,
402                                      Double_t kMultorCent,
403                                      Double_t vertexZ) {
404   // Calculates the balance function
405   fAnalyzedEvents++;
406     
407   // Initialize histograms if not done yet
408   if(!fHistPN){
409     AliWarning("Histograms not yet initialized! --> Will be done now");
410     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
411     InitHistograms();
412   }
413
414   Double_t trackVariablesSingle[kTrackVariablesSingle];
415   Double_t trackVariablesPair[kTrackVariablesPair];
416
417   if (!particles){
418     AliWarning("particles TObjArray is NULL pointer --> return");
419     return;
420   }
421   
422   // define end of particle loops
423   Int_t iMax = particles->GetEntriesFast();
424   Int_t jMax = iMax;
425   if (particlesMixed)
426     jMax = particlesMixed->GetEntriesFast();
427
428   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
429   TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
430
431   TArrayF secondEta(jMax);
432   TArrayF secondPhi(jMax);
433   TArrayF secondPt(jMax);
434   TArrayS secondCharge(jMax);
435   TArrayD secondCorrection(jMax);
436
437   for (Int_t i=0; i<jMax; i++){
438     secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
439     secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
440     secondPt[i]  = ((AliVParticle*) particlesSecond->At(i))->Pt();
441     secondCharge[i]  = (Short_t)((AliVParticle*) particlesSecond->At(i))->Charge();
442     secondCorrection[i]  = (Double_t)((AliBFBasicParticle*) particlesSecond->At(i))->Correction();   //==========================correction
443   }
444   
445   //TLorenzVector implementation for resonances
446   TLorentzVector vectorMother, vectorDaughter[2];
447   TParticle pPion, pProton, pRho0, pK0s, pLambda;
448   pPion.SetPdgCode(211); //pion
449   pRho0.SetPdgCode(113); //rho0
450   pK0s.SetPdgCode(310); //K0s
451   pProton.SetPdgCode(2212); //proton
452   pLambda.SetPdgCode(3122); //Lambda
453   Double_t gWidthForRho0 = 0.01;
454   Double_t gWidthForK0s = 0.01;
455   Double_t gWidthForLambda = 0.006;
456   Double_t nSigmaRejection = 3.0;
457
458   // 1st particle loop
459   for (Int_t i = 0; i < iMax; i++) {
460     //AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
461     AliBFBasicParticle* firstParticle = (AliBFBasicParticle*) particles->At(i); //==========================correction
462     
463     // some optimization
464     Float_t firstEta = firstParticle->Eta();
465     Float_t firstPhi = firstParticle->Phi();
466     Float_t firstPt  = firstParticle->Pt();
467     Float_t firstCorrection  = firstParticle->Correction();//==========================correction
468
469     // Event plane (determine psi bin)
470     Double_t gPsiMinusPhi    =   0.;
471     Double_t gPsiMinusPhiBin = -10.;
472     gPsiMinusPhi   = TMath::Abs(firstPhi - gReactionPlane);
473     //in-plane
474     if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
475        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
476       gPsiMinusPhiBin = 0.0;
477     //intermediate
478     else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
479             ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
480             ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
481             ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
482       gPsiMinusPhiBin = 1.0;
483     //out of plane
484     else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
485             ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
486       gPsiMinusPhiBin = 2.0;
487     //everything else
488     else 
489       gPsiMinusPhiBin = 3.0;
490     
491     fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
492
493     Short_t  charge1 = (Short_t) firstParticle->Charge();
494     
495     trackVariablesSingle[0]    =  gPsiMinusPhiBin;
496     trackVariablesSingle[1]    =  firstPt;
497       if(fEventClass=="Multiplicity" || fEventClass == "Centrality" ) trackVariablesSingle[0] = kMultorCent;
498     trackVariablesSingle[2]    =  vertexZ;
499
500     
501     //fill single particle histograms
502     if(charge1 > 0)      fHistP->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
503     else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,firstCorrection);  //==========================correction
504     
505     // 2nd particle loop
506     for(Int_t j = 0; j < jMax; j++) {   
507
508       if(!particlesMixed && j == i) continue; // no auto correlations (only for non mixing)
509
510       // pT,Assoc < pT,Trig
511       if(firstPt < secondPt[j]) continue;
512
513       Short_t charge2 = secondCharge[j];
514       
515       trackVariablesPair[0]    =  trackVariablesSingle[0];
516       trackVariablesPair[1]    =  firstEta - secondEta[j];  // delta eta
517       trackVariablesPair[2]    =  firstPhi - secondPhi[j];  // delta phi
518       //if (trackVariablesPair[2] > 180.)   // delta phi between -180 and 180 
519       //trackVariablesPair[2] -= 360.;
520       //if (trackVariablesPair[2] <  - 180.) 
521       //trackVariablesPair[2] += 360.;
522       if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi 
523         trackVariablesPair[2] -= 2.*TMath::Pi();
524       if (trackVariablesPair[2] <  - TMath::Pi()) 
525         trackVariablesPair[2] += 2.*TMath::Pi();
526       if (trackVariablesPair[2] <  - TMath::Pi()/2.) 
527       trackVariablesPair[2] += 2.*TMath::Pi();
528       
529       trackVariablesPair[3]    =  firstPt;      // pt trigger
530       trackVariablesPair[4]    =  secondPt[j];  // pt
531       trackVariablesPair[5]    =  vertexZ;      // z of the primary vertex
532       
533       //Exclude resonances for the calculation of pairs by looking 
534       //at the invariant mass and not considering the pairs that 
535       //fall within 3sigma from the mass peak of: rho0, K0s, Lambda
536       if(fResonancesCut) {
537         if (charge1 * charge2 < 0) {
538
539           //rho0
540           vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
541           vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
542           vectorMother = vectorDaughter[0] + vectorDaughter[1];
543           fHistResonancesBefore->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
544           if(TMath::Abs(vectorMother.M() - pRho0.GetMass()) <= nSigmaRejection*gWidthForRho0)
545             continue;
546           fHistResonancesRho->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
547           
548           //K0s
549           if(TMath::Abs(vectorMother.M() - pK0s.GetMass()) <= nSigmaRejection*gWidthForK0s)
550             continue;
551           fHistResonancesK0->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
552           
553           
554           //Lambda
555           vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
556           vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pProton.GetMass());
557           vectorMother = vectorDaughter[0] + vectorDaughter[1];
558           if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
559             continue;
560           
561           vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pProton.GetMass());
562           vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
563           vectorMother = vectorDaughter[0] + vectorDaughter[1];
564           if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
565             continue;
566           fHistResonancesLambda->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
567         
568         }//unlike-sign only
569       }//resonance cut
570
571       // HBT like cut
572       //if(fHBTCut){ // VERSION 3 (all pairs)
573       if(fHBTCut && charge1 * charge2 > 0){  // VERSION 2 (only for LS)
574         //if( dphi < 3 || deta < 0.01 ){   // VERSION 1
575         //  continue;
576         
577         Double_t deta = firstEta - secondEta[j];
578         Double_t dphi = firstPhi - secondPhi[j];
579         // VERSION 2 (Taken from DPhiCorrelations)
580         // the variables & cuthave been developed by the HBT group 
581         // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
582         fHistHBTbefore->Fill(deta,dphi);
583         
584         // optimization
585         if (TMath::Abs(deta) < fHBTCutValue * 2.5 * 3) //fHBTCutValue = 0.02 [default for dphicorrelations]
586           {
587             // phi in rad
588             //Float_t phi1rad = firstPhi*TMath::DegToRad();
589             //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
590             Float_t phi1rad = firstPhi;
591             Float_t phi2rad = secondPhi[j];
592             
593             // check first boundaries to see if is worth to loop and find the minimum
594             Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
595             Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
596             
597             const Float_t kLimit = fHBTCutValue * 3;
598             
599             Float_t dphistarminabs = 1e5;
600             Float_t dphistarmin = 1e5;
601             
602             if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
603               for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
604                 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
605                 Float_t dphistarabs = TMath::Abs(dphistar);
606                 
607                 if (dphistarabs < dphistarminabs) {
608                   dphistarmin = dphistar;
609                   dphistarminabs = dphistarabs;
610                 }
611               }
612               
613               if (dphistarminabs < fHBTCutValue && TMath::Abs(deta) < fHBTCutValue) {
614                 //AliInfo(Form("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f", i, j, deta, dphi, dphistarminabs, dphistar1, dphistar2, phi1rad, pt1, charge1, phi2rad, pt2, charge2, bSign));
615                 continue;
616               }
617             }
618           }
619         fHistHBTafter->Fill(deta,dphi);
620       }//HBT cut
621         
622       // conversions
623       if(fConversionCut) {
624         if (charge1 * charge2 < 0) {
625           Double_t deta = firstEta - secondEta[j];
626           Double_t dphi = firstPhi - secondPhi[j];
627           
628           Float_t m0 = 0.510e-3;
629           Float_t tantheta1 = 1e10;
630           
631           // phi in rad
632           //Float_t phi1rad = firstPhi*TMath::DegToRad();
633           //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
634           Float_t phi1rad = firstPhi;
635           Float_t phi2rad = secondPhi[j];
636           
637           if (firstEta < -1e-10 || firstEta > 1e-10)
638             tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
639           
640           Float_t tantheta2 = 1e10;
641           if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
642             tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
643           
644           Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
645           Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
646           
647           Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
648
649           fHistConversionbefore->Fill(deta,dphi,masssqu);
650           
651           if (masssqu < fInvMassCutConversion*fInvMassCutConversion){
652             //AliInfo(Form("Conversion: Removed track pair %d %d with [[%f %f] %f %f] %d %d <- %f %f  %f %f   %f %f ", i, j, deta, dphi, masssqu, charge1, charge2,eta1,eta2,phi1,phi2,pt1,pt2));
653             continue;
654           }
655           fHistConversionafter->Fill(deta,dphi,masssqu);
656         }
657       }//conversion cut
658
659       // momentum difference cut - suppress femtoscopic effects
660       if(fQCut){ 
661
662         //Double_t ptMin        = 0.1; //const for the time being (should be changeable later on)
663         Double_t ptDifference = TMath::Abs( firstPt - secondPt[j]);
664
665         fHistQbefore->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
666         if(ptDifference < fDeltaPtMin) continue;
667         fHistQafter->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
668
669       }
670
671       if( charge1 > 0 && charge2 < 0)  fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
672       else if( charge1 < 0 && charge2 > 0)  fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction 
673       else if( charge1 > 0 && charge2 > 0)  fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction 
674       else if( charge1 < 0 && charge2 < 0)  fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction 
675       else {
676         //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
677         continue;
678       }
679     }//end of 2nd particle loop
680   }//end of 1st particle loop
681 }  
682
683 //____________________________________________________________________//
684 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
685                                                  Int_t iVariablePair,
686                                                  Double_t psiMin, 
687                                                  Double_t psiMax,
688                                                  Double_t vertexZMin,
689                                                  Double_t vertexZMax,
690                                                  Double_t ptTriggerMin,
691                                                  Double_t ptTriggerMax,
692                                                  Double_t ptAssociatedMin,
693                                                  Double_t ptAssociatedMax) {
694   //Returns the BF histogram, extracted from the 6 AliTHn objects 
695   //(private members) of the AliBalancePsi class.
696   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
697   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
698
699   // security checks
700   if(psiMin > psiMax-0.00001){
701     AliError("psiMax <= psiMin");
702     return NULL;
703   }
704   if(vertexZMin > vertexZMax-0.00001){
705     AliError("vertexZMax <= vertexZMin");
706     return NULL;
707   }
708   if(ptTriggerMin > ptTriggerMax-0.00001){
709     AliError("ptTriggerMax <= ptTriggerMin");
710     return NULL;
711   }
712   if(ptAssociatedMin > ptAssociatedMax-0.00001){
713     AliError("ptAssociatedMax <= ptAssociatedMin");
714     return NULL;
715   }
716
717   // Psi_2
718   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
719   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
720   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
721   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
722   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
723   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
724
725   // Vz
726  if(fVertexBinning){
727    fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
728    fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
729    fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
730    fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
731    fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
732    fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
733  }
734
735   // pt trigger
736   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
737     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
738     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
739     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
740     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
741     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
742     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
743   }
744
745   // pt associated
746   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
747     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
748     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
749     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
750     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
751   }
752
753   //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));
754
755   // Project into the wanted space (1st: analysis step, 2nd: axis)
756   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair); //
757   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair); //
758   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair); //
759   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair); //
760   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle); //
761   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle); //
762
763   TH1D *gHistBalanceFunctionHistogram = 0x0;
764   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
765     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
766     gHistBalanceFunctionHistogram->Reset();
767     
768     switch(iVariablePair) {
769     case 1:
770       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
771       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
772       break;
773     case 2:
774       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
775       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
776       break;
777     default:
778       break;
779     }
780
781     hTemp1->Sumw2();
782     hTemp2->Sumw2();
783     hTemp3->Sumw2();
784     hTemp4->Sumw2();
785     hTemp1->Add(hTemp3,-1.);
786     hTemp1->Scale(1./hTemp5->Integral());
787     hTemp2->Add(hTemp4,-1.);
788     hTemp2->Scale(1./hTemp6->Integral());
789     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
790     gHistBalanceFunctionHistogram->Scale(0.5);
791
792     //normalize to bin width
793     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
794   }
795
796   return gHistBalanceFunctionHistogram;
797 }
798
799 //____________________________________________________________________//
800 TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
801                                                          Int_t iVariablePair,
802                                                          Double_t psiMin, 
803                                                          Double_t psiMax,
804                                                          Double_t vertexZMin,
805                                                          Double_t vertexZMax,
806                                                          Double_t ptTriggerMin,
807                                                          Double_t ptTriggerMax,
808                                                          Double_t ptAssociatedMin,
809                                                          Double_t ptAssociatedMax,
810                                                          AliBalancePsi *bfMix) {
811   //Returns the BF histogram, extracted from the 6 AliTHn objects 
812   //after dividing each correlation function by the Event Mixing one 
813   //(private members) of the AliBalancePsi class.
814   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
815   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
816
817   // security checks
818   if(psiMin > psiMax-0.00001){
819     AliError("psiMax <= psiMin");
820     return NULL;
821   }
822   if(vertexZMin > vertexZMax-0.00001){
823     AliError("vertexZMax <= vertexZMin");
824     return NULL;
825   }
826   if(ptTriggerMin > ptTriggerMax-0.00001){
827     AliError("ptTriggerMax <= ptTriggerMin");
828     return NULL;
829   }
830   if(ptAssociatedMin > ptAssociatedMax-0.00001){
831     AliError("ptAssociatedMax <= ptAssociatedMin");
832     return NULL;
833   }
834
835   // pt trigger (fixed for all vertices/psi/centralities)
836   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
837     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
838     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
839     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
840     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
841     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
842     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
843   }
844
845   // pt associated (fixed for all vertices/psi/centralities)
846   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
847     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
848     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
849     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
850     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
851   }
852
853   // ============================================================================================
854   // the same for event mixing
855   AliTHn *fHistPMix = bfMix->GetHistNp();
856   AliTHn *fHistNMix = bfMix->GetHistNn();
857   AliTHn *fHistPNMix = bfMix->GetHistNpn();
858   AliTHn *fHistNPMix = bfMix->GetHistNnp();
859   AliTHn *fHistPPMix = bfMix->GetHistNpp();
860   AliTHn *fHistNNMix = bfMix->GetHistNnn();
861
862   // pt trigger (fixed for all vertices/psi/centralities)
863   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
864     fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
865     fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
866     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
867     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
868     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
869     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
870   }
871
872   // pt associated (fixed for all vertices/psi/centralities)
873   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
874     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
875     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
876     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
877     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
878   }
879
880   // ============================================================================================
881   // ranges (in bins) for vertexZ and centrality (psi)
882   Int_t binPsiMin    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
883   Int_t binPsiMax    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
884   Int_t binVertexMin = 0;
885   Int_t binVertexMax = 0;
886   if(fVertexBinning){
887     binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
888     binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
889   }  
890   Double_t binPsiLowEdge    = 0.;
891   Double_t binPsiUpEdge     = 1000.;
892   Double_t binVertexLowEdge = 0.;
893   Double_t binVertexUpEdge  = 1000.;
894
895   TH1D* h1 = NULL;
896   TH1D* h2 = NULL;
897   TH1D* h3 = NULL;
898   TH1D* h4 = NULL;
899
900   // loop over all bins
901   for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
902     for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
903   
904       cout<<"In the balance function (1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin)  "<<endl;
905
906       // determine the bin edges for this bin
907       binPsiLowEdge    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
908       binPsiUpEdge     = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
909       if(fVertexBinning){
910         binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
911         binVertexUpEdge  = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
912       }
913       
914       // Psi_2
915       fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
916       fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
917       fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
918       fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
919       fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
920       fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
921
922       // Vz
923       if(fVertexBinning){
924         fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
925         fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
926         fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
927         fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
928         fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
929         fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
930       }
931
932       // ============================================================================================
933       // the same for event mixing
934       
935       // Psi_2
936       fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
937       fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
938       fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
939       fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
940       fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
941       fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
942       
943       // Vz
944       if(fVertexBinning){
945         fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
946         fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
947         fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
948         fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
949         fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
950         fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
951       }
952       
953       // ============================================================================================
954
955       //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));
956       
957       // Project into the wanted space (1st: analysis step, 2nd: axis)
958       TH1D* hTempHelper1 = (TH1D*)fHistPN->Project(0,iVariablePair);
959       TH1D* hTempHelper2 = (TH1D*)fHistNP->Project(0,iVariablePair);
960       TH1D* hTempHelper3 = (TH1D*)fHistPP->Project(0,iVariablePair);
961       TH1D* hTempHelper4 = (TH1D*)fHistNN->Project(0,iVariablePair);
962       TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
963       TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
964       
965       // ============================================================================================
966       // the same for event mixing
967       TH1D* hTempHelper1Mix = (TH1D*)fHistPNMix->Project(0,iVariablePair);
968       TH1D* hTempHelper2Mix = (TH1D*)fHistNPMix->Project(0,iVariablePair);
969       TH1D* hTempHelper3Mix = (TH1D*)fHistPPMix->Project(0,iVariablePair);
970       TH1D* hTempHelper4Mix = (TH1D*)fHistNNMix->Project(0,iVariablePair);
971       TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
972       TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
973       // ============================================================================================
974
975       hTempHelper1->Sumw2();
976       hTempHelper2->Sumw2();
977       hTempHelper3->Sumw2();
978       hTempHelper4->Sumw2();
979       hTemp5->Sumw2();
980       hTemp6->Sumw2();
981       hTempHelper1Mix->Sumw2();
982       hTempHelper2Mix->Sumw2();
983       hTempHelper3Mix->Sumw2();
984       hTempHelper4Mix->Sumw2();
985       hTemp5Mix->Sumw2();
986       hTemp6Mix->Sumw2();
987
988       // first put everything on positive x - axis (to be comparable with published data) --> ONLY IN THIS OPTION!
989
990       Double_t helperEndBin = 1.6;
991       if(iVariablePair==2) helperEndBin = TMath::Pi();
992
993       TH1D* hTempPos1 = new TH1D(Form("hTempPos1_%d_%d",iBinPsi,iBinVertex),Form("hTempPos1_%d_%d",iBinPsi,iBinVertex),hTempHelper1->GetNbinsX()/2,0,helperEndBin);
994       TH1D* hTempPos2 = new TH1D(Form("hTempPos2_%d_%d",iBinPsi,iBinVertex),Form("hTempPos2_%d_%d",iBinPsi,iBinVertex),hTempHelper2->GetNbinsX()/2,0,helperEndBin);
995       TH1D* hTempPos3 = new TH1D(Form("hTempPos3_%d_%d",iBinPsi,iBinVertex),Form("hTempPos3_%d_%d",iBinPsi,iBinVertex),hTempHelper3->GetNbinsX()/2,0,helperEndBin);
996       TH1D* hTempPos4 = new TH1D(Form("hTempPos4_%d_%d",iBinPsi,iBinVertex),Form("hTempPos4_%d_%d",iBinPsi,iBinVertex),hTempHelper4->GetNbinsX()/2,0,helperEndBin);
997       TH1D* hTempPos1Mix = new TH1D(Form("hTempPos1Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos1Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper1Mix->GetNbinsX()/2,0,helperEndBin);
998       TH1D* hTempPos2Mix = new TH1D(Form("hTempPos2Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos2Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper2Mix->GetNbinsX()/2,0,helperEndBin);
999       TH1D* hTempPos3Mix = new TH1D(Form("hTempPos3Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos3Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper3Mix->GetNbinsX()/2,0,helperEndBin);
1000       TH1D* hTempPos4Mix = new TH1D(Form("hTempPos4Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos4Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper4Mix->GetNbinsX()/2,0,helperEndBin);
1001
1002       TH1D* hTemp1 = new TH1D(Form("hTemp1_%d_%d",iBinPsi,iBinVertex),Form("hTemp1_%d_%d",iBinPsi,iBinVertex),hTempHelper1->GetNbinsX()/2,0,helperEndBin);
1003       TH1D* hTemp2 = new TH1D(Form("hTemp2_%d_%d",iBinPsi,iBinVertex),Form("hTemp2_%d_%d",iBinPsi,iBinVertex),hTempHelper2->GetNbinsX()/2,0,helperEndBin);
1004       TH1D* hTemp3 = new TH1D(Form("hTemp3_%d_%d",iBinPsi,iBinVertex),Form("hTemp3_%d_%d",iBinPsi,iBinVertex),hTempHelper3->GetNbinsX()/2,0,helperEndBin);
1005       TH1D* hTemp4 = new TH1D(Form("hTemp4_%d_%d",iBinPsi,iBinVertex),Form("hTemp4_%d_%d",iBinPsi,iBinVertex),hTempHelper4->GetNbinsX()/2,0,helperEndBin);
1006       TH1D* hTemp1Mix = new TH1D(Form("hTemp1Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp1Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper1Mix->GetNbinsX()/2,0,helperEndBin);
1007       TH1D* hTemp2Mix = new TH1D(Form("hTemp2Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp2Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper2Mix->GetNbinsX()/2,0,helperEndBin);
1008       TH1D* hTemp3Mix = new TH1D(Form("hTemp3Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp3Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper3Mix->GetNbinsX()/2,0,helperEndBin);
1009       TH1D* hTemp4Mix = new TH1D(Form("hTemp4Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp4Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper4Mix->GetNbinsX()/2,0,helperEndBin);
1010
1011
1012       hTempPos1->Sumw2();
1013       hTempPos2->Sumw2();
1014       hTempPos3->Sumw2();
1015       hTempPos4->Sumw2();
1016       hTempPos1Mix->Sumw2();
1017       hTempPos2Mix->Sumw2();
1018       hTempPos3Mix->Sumw2();
1019       hTempPos4Mix->Sumw2();
1020
1021       hTemp1->Sumw2();
1022       hTemp2->Sumw2();
1023       hTemp3->Sumw2();
1024       hTemp4->Sumw2();
1025       hTemp1Mix->Sumw2();
1026       hTemp2Mix->Sumw2();
1027       hTemp3Mix->Sumw2();
1028       hTemp4Mix->Sumw2();
1029
1030
1031       for(Int_t i=0;i<hTempHelper1->GetNbinsX();i++){
1032
1033         Double_t binCenter  = hTempHelper1->GetXaxis()->GetBinCenter(i+1);
1034
1035         if(iVariablePair==1){
1036           if(binCenter>0){
1037             hTempPos1->SetBinContent(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinContent(i+1));
1038             hTempPos2->SetBinContent(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinContent(i+1));
1039             hTempPos3->SetBinContent(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinContent(i+1));
1040             hTempPos4->SetBinContent(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinContent(i+1));
1041             hTempPos1Mix->SetBinContent(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinContent(i+1));
1042             hTempPos2Mix->SetBinContent(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinContent(i+1));
1043             hTempPos3Mix->SetBinContent(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinContent(i+1));
1044             hTempPos4Mix->SetBinContent(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinContent(i+1));
1045
1046             hTempPos1->SetBinError(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinError(i+1));
1047             hTempPos2->SetBinError(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinError(i+1));
1048             hTempPos3->SetBinError(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinError(i+1));
1049             hTempPos4->SetBinError(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinError(i+1));
1050             hTempPos1Mix->SetBinError(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinError(i+1));
1051             hTempPos2Mix->SetBinError(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinError(i+1));
1052             hTempPos3Mix->SetBinError(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinError(i+1));
1053             hTempPos4Mix->SetBinError(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinError(i+1));
1054           }
1055           else{
1056             hTemp1->SetBinContent(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinContent(i+1));
1057             hTemp2->SetBinContent(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinContent(i+1));
1058             hTemp3->SetBinContent(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinContent(i+1));
1059             hTemp4->SetBinContent(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinContent(i+1));
1060             hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1061             hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1062             hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1063             hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1064
1065             hTemp1->SetBinError(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinError(i+1));
1066             hTemp2->SetBinError(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinError(i+1));
1067             hTemp3->SetBinError(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinError(i+1));
1068             hTemp4->SetBinError(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinError(i+1));
1069             hTemp1Mix->SetBinError(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinError(i+1));
1070             hTemp2Mix->SetBinError(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinError(i+1));
1071             hTemp3Mix->SetBinError(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinError(i+1));
1072             hTemp4Mix->SetBinError(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinError(i+1));
1073           }
1074         }
1075         else if(iVariablePair==2){
1076           if(binCenter>0 && binCenter<TMath::Pi()){
1077             hTempPos1->SetBinContent(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinContent(i+1));
1078             hTempPos2->SetBinContent(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinContent(i+1));
1079             hTempPos3->SetBinContent(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinContent(i+1));
1080             hTempPos4->SetBinContent(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinContent(i+1));
1081             hTempPos1Mix->SetBinContent(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinContent(i+1));
1082             hTempPos2Mix->SetBinContent(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinContent(i+1));
1083             hTempPos3Mix->SetBinContent(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinContent(i+1));
1084             hTempPos4Mix->SetBinContent(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinContent(i+1));
1085
1086             hTempPos1->SetBinError(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinError(i+1));
1087             hTempPos2->SetBinError(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinError(i+1));
1088             hTempPos3->SetBinError(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinError(i+1));
1089             hTempPos4->SetBinError(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinError(i+1));
1090             hTempPos1Mix->SetBinError(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinError(i+1));
1091             hTempPos2Mix->SetBinError(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinError(i+1));
1092             hTempPos3Mix->SetBinError(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinError(i+1));
1093             hTempPos4Mix->SetBinError(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinError(i+1));
1094           }
1095           else if(binCenter<0){
1096             hTemp1->SetBinContent(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinContent(i+1));
1097             hTemp2->SetBinContent(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinContent(i+1));
1098             hTemp3->SetBinContent(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinContent(i+1));
1099             hTemp4->SetBinContent(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinContent(i+1));
1100             hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1101             hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1102             hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1103             hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1104
1105             hTemp1->SetBinError(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinError(i+1));
1106             hTemp2->SetBinError(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinError(i+1));
1107             hTemp3->SetBinError(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinError(i+1));
1108             hTemp4->SetBinError(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinError(i+1));
1109             hTemp1Mix->SetBinError(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinError(i+1));
1110             hTemp2Mix->SetBinError(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinError(i+1));
1111             hTemp3Mix->SetBinError(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinError(i+1));
1112             hTemp4Mix->SetBinError(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinError(i+1));
1113           }
1114           else{
1115             hTemp1->SetBinContent(hTemp1->FindBin(TMath::Pi()-binCenter),hTempHelper1->GetBinContent(i+1));
1116             hTemp2->SetBinContent(hTemp2->FindBin(TMath::Pi()-binCenter),hTempHelper2->GetBinContent(i+1));
1117             hTemp3->SetBinContent(hTemp3->FindBin(TMath::Pi()-binCenter),hTempHelper3->GetBinContent(i+1));
1118             hTemp4->SetBinContent(hTemp4->FindBin(TMath::Pi()-binCenter),hTempHelper4->GetBinContent(i+1));
1119             hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(TMath::Pi()-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1120             hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(TMath::Pi()-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1121             hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(TMath::Pi()-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1122             hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(TMath::Pi()-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1123
1124             hTemp1->SetBinError(hTemp1->FindBin(TMath::Pi()-binCenter),hTempHelper1->GetBinError(i+1));
1125             hTemp2->SetBinError(hTemp2->FindBin(TMath::Pi()-binCenter),hTempHelper2->GetBinError(i+1));
1126             hTemp3->SetBinError(hTemp3->FindBin(TMath::Pi()-binCenter),hTempHelper3->GetBinError(i+1));
1127             hTemp4->SetBinError(hTemp4->FindBin(TMath::Pi()-binCenter),hTempHelper4->GetBinError(i+1));
1128             hTemp1Mix->SetBinError(hTemp1Mix->FindBin(TMath::Pi()-binCenter),hTempHelper1Mix->GetBinError(i+1));
1129             hTemp2Mix->SetBinError(hTemp2Mix->FindBin(TMath::Pi()-binCenter),hTempHelper2Mix->GetBinError(i+1));
1130             hTemp3Mix->SetBinError(hTemp3Mix->FindBin(TMath::Pi()-binCenter),hTempHelper3Mix->GetBinError(i+1));
1131             hTemp4Mix->SetBinError(hTemp4Mix->FindBin(TMath::Pi()-binCenter),hTempHelper4Mix->GetBinError(i+1));
1132           }
1133         }
1134       }
1135
1136       hTemp1->Add(hTempPos1);
1137       hTemp2->Add(hTempPos2);
1138       hTemp3->Add(hTempPos3);
1139       hTemp4->Add(hTempPos4);
1140       hTemp1Mix->Add(hTempPos1Mix);
1141       hTemp2Mix->Add(hTempPos2Mix);
1142       hTemp3Mix->Add(hTempPos3Mix);
1143       hTemp4Mix->Add(hTempPos4Mix);
1144
1145       hTemp1->Scale(1./hTemp5->Integral());
1146       hTemp3->Scale(1./hTemp5->Integral());
1147       hTemp2->Scale(1./hTemp6->Integral());
1148       hTemp4->Scale(1./hTemp6->Integral());
1149
1150       // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1151       // does not work here, so normalize also to trigger particles 
1152       // --> careful: gives different integrals then as with full 2D method 
1153       hTemp1Mix->Scale(1./hTemp5Mix->Integral());
1154       hTemp3Mix->Scale(1./hTemp5Mix->Integral());
1155       hTemp2Mix->Scale(1./hTemp6Mix->Integral());
1156       hTemp4Mix->Scale(1./hTemp6Mix->Integral());
1157
1158       hTemp1->Divide(hTemp1Mix);
1159       hTemp2->Divide(hTemp2Mix);
1160       hTemp3->Divide(hTemp3Mix);
1161       hTemp4->Divide(hTemp4Mix);
1162
1163       // for the first: clone
1164       if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1165         h1 = (TH1D*)hTemp1->Clone();
1166         h2 = (TH1D*)hTemp2->Clone();
1167         h3 = (TH1D*)hTemp3->Clone();
1168         h4 = (TH1D*)hTemp4->Clone();
1169       }
1170       else{  // otherwise: add for averaging
1171         h1->Add(hTemp1);
1172         h2->Add(hTemp2);
1173         h3->Add(hTemp3);
1174         h4->Add(hTemp4);
1175       }
1176
1177       delete hTemp1;
1178       delete hTemp2;
1179       delete hTemp3;
1180       delete hTemp4;
1181       delete hTemp1Mix;
1182       delete hTemp2Mix;
1183       delete hTemp3Mix;
1184       delete hTemp4Mix;
1185
1186       delete hTempPos1;
1187       delete hTempPos2;
1188       delete hTempPos3;
1189       delete hTempPos4;
1190       delete hTempPos1Mix;
1191       delete hTempPos2Mix;
1192       delete hTempPos3Mix;
1193       delete hTempPos4Mix;
1194     }
1195   }
1196
1197   TH1D *gHistBalanceFunctionHistogram = 0x0;
1198   if((h1)&&(h2)&&(h3)&&(h4)) {
1199
1200     // average over number of bins nbinsVertex * nbinsPsi
1201     h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1202     h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1203     h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1204     h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1205
1206     gHistBalanceFunctionHistogram = (TH1D*)h1->Clone();
1207     gHistBalanceFunctionHistogram->Reset();
1208     
1209     switch(iVariablePair) {
1210     case 1:
1211       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1212       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1213       break;
1214     case 2:
1215       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1216       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1217       break;
1218     default:
1219       break;
1220     }
1221
1222     h1->Add(h3,-1.);
1223     h2->Add(h4,-1.);
1224
1225     gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.);
1226     gHistBalanceFunctionHistogram->Scale(0.5);
1227
1228     //normalize to bin width
1229     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
1230   }
1231
1232   return gHistBalanceFunctionHistogram;
1233 }
1234
1235 //____________________________________________________________________//
1236 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin, 
1237                                                         Double_t psiMax,
1238                                                         Double_t vertexZMin,
1239                                                         Double_t vertexZMax,
1240                                                         Double_t ptTriggerMin,
1241                                                         Double_t ptTriggerMax,
1242                                                         Double_t ptAssociatedMin,
1243                                                         Double_t ptAssociatedMax) {
1244   //Returns the BF histogram in Delta eta vs Delta phi, 
1245   //extracted from the 6 AliTHn objects 
1246   //(private members) of the AliBalancePsi class.
1247   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1248   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1249
1250   // security checks
1251   if(psiMin > psiMax-0.00001){
1252     AliError("psiMax <= psiMin");
1253     return NULL;
1254   }
1255   if(vertexZMin > vertexZMax-0.00001){
1256     AliError("vertexZMax <= vertexZMin");
1257     return NULL;
1258   }
1259   if(ptTriggerMin > ptTriggerMax-0.00001){
1260     AliError("ptTriggerMax <= ptTriggerMin");
1261     return NULL;
1262   }
1263   if(ptAssociatedMin > ptAssociatedMax-0.00001){
1264     AliError("ptAssociatedMax <= ptAssociatedMin");
1265     return NULL;
1266   }
1267
1268   TString histName = "gHistBalanceFunctionHistogram2D";
1269
1270   // Psi_2
1271   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
1272   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
1273   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
1274   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
1275   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
1276   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
1277
1278   // Vz
1279   if(fVertexBinning){
1280     fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1281     fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1282     fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1283     fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1284     fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1285     fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1286   }
1287
1288   // pt trigger
1289   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1290     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1291     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1292     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1293     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1294     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1295     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1296   }
1297
1298   // pt associated
1299   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1300     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1301     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1302     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1303     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1304   }
1305
1306   //AliInfo(Form("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)));
1307
1308   // Project into the wanted space (1st: analysis step, 2nd: axis)
1309   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1310   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1311   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1312   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1313   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1314   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1315
1316   TH2D *gHistBalanceFunctionHistogram = 0x0;
1317   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1318     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1319     gHistBalanceFunctionHistogram->Reset();
1320     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");   
1321     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1322     gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");   
1323     
1324     hTemp1->Sumw2();
1325     hTemp2->Sumw2();
1326     hTemp3->Sumw2();
1327     hTemp4->Sumw2();
1328     hTemp1->Add(hTemp3,-1.);
1329     hTemp1->Scale(1./hTemp5->Integral());
1330     hTemp2->Add(hTemp4,-1.);
1331     hTemp2->Scale(1./hTemp6->Integral());
1332     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1333     gHistBalanceFunctionHistogram->Scale(0.5);
1334
1335     //normalize to bin width
1336     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
1337   }
1338
1339   return gHistBalanceFunctionHistogram;
1340 }
1341
1342 //____________________________________________________________________//
1343 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin, 
1344                                                                 Double_t psiMax,
1345                                                                 Double_t vertexZMin,
1346                                                                 Double_t vertexZMax,
1347                                                                 Double_t ptTriggerMin,
1348                                                                 Double_t ptTriggerMax,
1349                                                                 Double_t ptAssociatedMin,
1350                                                                 Double_t ptAssociatedMax,
1351                                                                 AliBalancePsi *bfMix) {
1352   //Returns the BF histogram in Delta eta vs Delta phi,
1353   //after dividing each correlation function by the Event Mixing one 
1354   //extracted from the 6 AliTHn objects 
1355   //(private members) of the AliBalancePsi class.
1356   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1357   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1358   TString histName = "gHistBalanceFunctionHistogram2D";
1359
1360   if(!bfMix){
1361     AliError("balance function object for event mixing not available");
1362     return NULL;
1363   }
1364
1365   // security checks
1366   if(psiMin > psiMax-0.00001){
1367     AliError("psiMax <= psiMin");
1368     return NULL;
1369   }
1370   if(vertexZMin > vertexZMax-0.00001){
1371     AliError("vertexZMax <= vertexZMin");
1372     return NULL;
1373   }
1374   if(ptTriggerMin > ptTriggerMax-0.00001){
1375     AliError("ptTriggerMax <= ptTriggerMin");
1376     return NULL;
1377   }
1378   if(ptAssociatedMin > ptAssociatedMax-0.00001){
1379     AliError("ptAssociatedMax <= ptAssociatedMin");
1380     return NULL;
1381   }
1382
1383   // pt trigger (fixed for all vertices/psi/centralities)
1384   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1385     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1386     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1387     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1388     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1389     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1390     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1391   }
1392
1393   // pt associated (fixed for all vertices/psi/centralities)
1394   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1395     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1396     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1397     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1398     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1399   }
1400
1401   // ============================================================================================
1402   // the same for event mixing
1403   AliTHn *fHistPMix = bfMix->GetHistNp();
1404   AliTHn *fHistNMix = bfMix->GetHistNn();
1405   AliTHn *fHistPNMix = bfMix->GetHistNpn();
1406   AliTHn *fHistNPMix = bfMix->GetHistNnp();
1407   AliTHn *fHistPPMix = bfMix->GetHistNpp();
1408   AliTHn *fHistNNMix = bfMix->GetHistNnn();
1409
1410   // pt trigger (fixed for all vertices/psi/centralities)
1411   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1412     fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1413     fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1414     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1415     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1416     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1417     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1418   }
1419
1420   // pt associated (fixed for all vertices/psi/centralities)
1421   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1422     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1423     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1424     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1425     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1426   }
1427
1428   // ============================================================================================
1429   // ranges (in bins) for vertexZ and centrality (psi)
1430   Int_t binPsiMin    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1431   Int_t binPsiMax    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1432   Int_t binVertexMin = 0;
1433   Int_t binVertexMax = 0;
1434   if(fVertexBinning){
1435     binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1436     binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1437   }  
1438   Double_t binPsiLowEdge    = 0.;
1439   Double_t binPsiUpEdge     = 0.;
1440   Double_t binVertexLowEdge = 0.;
1441   Double_t binVertexUpEdge  = 0.;
1442
1443   TH2D* h1 = NULL;
1444   TH2D* h2 = NULL;
1445   TH2D* h3 = NULL;
1446   TH2D* h4 = NULL;
1447
1448   // loop over all bins
1449   for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1450     for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1451   
1452       cout<<"In the balance function (2D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin)  "<<endl;
1453
1454       // determine the bin edges for this bin
1455       binPsiLowEdge    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1456       binPsiUpEdge     = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
1457       if(fVertexBinning){
1458         binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1459         binVertexUpEdge  = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1460       }
1461
1462       // Psi_2
1463       fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1464       fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1465       fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1466       fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1467       fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1468       fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1469   
1470       // Vz
1471       if(fVertexBinning){
1472         fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1473         fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1474         fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1475         fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1476         fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1477         fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1478       }
1479
1480
1481
1482       // ============================================================================================
1483       // the same for event mixing
1484       
1485       // Psi_2
1486       fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1487       fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1488       fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1489       fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1490       fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1491       fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); 
1492       
1493       // Vz
1494       if(fVertexBinning){
1495         fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1496         fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1497         fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1498         fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1499         fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1500         fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); 
1501       }
1502       
1503       // ============================================================================================
1504       
1505       
1506       //AliInfo(Form("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)));
1507
1508       // Project into the wanted space (1st: analysis step, 2nd: axis)
1509       TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1510       TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1511       TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1512       TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1513       TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1514       TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1515       
1516       // ============================================================================================
1517       // the same for event mixing
1518       TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1519       TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1520       TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1521       TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1522       // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1523       // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1524       // ============================================================================================
1525       
1526       hTemp1->Sumw2();
1527       hTemp2->Sumw2();
1528       hTemp3->Sumw2();
1529       hTemp4->Sumw2();
1530       hTemp1Mix->Sumw2();
1531       hTemp2Mix->Sumw2();
1532       hTemp3Mix->Sumw2();
1533       hTemp4Mix->Sumw2();
1534       
1535       hTemp1->Scale(1./hTemp5->Integral());
1536       hTemp3->Scale(1./hTemp5->Integral());
1537       hTemp2->Scale(1./hTemp6->Integral());
1538       hTemp4->Scale(1./hTemp6->Integral());
1539
1540       // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1541       Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1542       mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1543       hTemp1Mix->Scale(1./mixedNorm1);
1544       Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1545       mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1546       hTemp2Mix->Scale(1./mixedNorm2);
1547       Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1548       mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1549       hTemp3Mix->Scale(1./mixedNorm3);
1550       Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1551       mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1552       hTemp4Mix->Scale(1./mixedNorm4);
1553       
1554       hTemp1->Divide(hTemp1Mix);
1555       hTemp2->Divide(hTemp2Mix);
1556       hTemp3->Divide(hTemp3Mix);
1557       hTemp4->Divide(hTemp4Mix);
1558       
1559       // for the first: clone
1560       if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1561         h1 = (TH2D*)hTemp1->Clone();
1562         h2 = (TH2D*)hTemp2->Clone();
1563         h3 = (TH2D*)hTemp3->Clone();
1564         h4 = (TH2D*)hTemp4->Clone();
1565       }
1566       else{  // otherwise: add for averaging
1567         h1->Add(hTemp1);
1568         h2->Add(hTemp2);
1569         h3->Add(hTemp3);
1570         h4->Add(hTemp4);
1571       } 
1572     }
1573   }
1574    
1575   TH2D *gHistBalanceFunctionHistogram = 0x0;
1576   if((h1)&&(h2)&&(h3)&&(h4)) {
1577
1578     // average over number of bins nbinsVertex * nbinsPsi
1579     h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1580     h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1581     h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1582     h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1583
1584     gHistBalanceFunctionHistogram = (TH2D*)h1->Clone();
1585     gHistBalanceFunctionHistogram->Reset();
1586     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");   
1587     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1588     gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");   
1589
1590     h1->Add(h3,-1.);
1591     h2->Add(h4,-1.);
1592     
1593     gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.);
1594     gHistBalanceFunctionHistogram->Scale(0.5);
1595     
1596     //normalize to bin width
1597     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
1598   }
1599   
1600   return gHistBalanceFunctionHistogram;
1601 }
1602
1603 //____________________________________________________________________//
1604 TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
1605                                                         Double_t psiMin,
1606                                                         Double_t psiMax,
1607                                                         Double_t vertexZMin,
1608                                                         Double_t vertexZMax,
1609                                                         Double_t ptTriggerMin,
1610                                                         Double_t ptTriggerMax,
1611                                                         Double_t ptAssociatedMin,
1612                                                         Double_t ptAssociatedMax,
1613                                                         AliBalancePsi *bfMix) {
1614   //Returns the BF histogram in Delta eta OR Delta phi,
1615   //after dividing each correlation function by the Event Mixing one
1616   // (But the division is done here in 2D, this was basically done to check the Integral)
1617   //extracted from the 6 AliTHn objects
1618   //(private members) of the AliBalancePsi class.
1619   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1620   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1621   TString histName = "gHistBalanceFunctionHistogram1D";
1622
1623   if(!bfMix){
1624     AliError("balance function object for event mixing not available");
1625     return NULL;
1626   }
1627
1628   // security checks
1629   if(psiMin > psiMax-0.00001){
1630     AliError("psiMax <= psiMin");
1631     return NULL;
1632   }
1633   if(vertexZMin > vertexZMax-0.00001){
1634     AliError("vertexZMax <= vertexZMin");
1635     return NULL;
1636   }
1637   if(ptTriggerMin > ptTriggerMax-0.00001){
1638     AliError("ptTriggerMax <= ptTriggerMin");
1639     return NULL;
1640   }
1641   if(ptAssociatedMin > ptAssociatedMax-0.00001){
1642     AliError("ptAssociatedMax <= ptAssociatedMin");
1643     return NULL;
1644   }
1645
1646   // pt trigger (fixed for all vertices/psi/centralities)
1647   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1648     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1649     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1650     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1651     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1652     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1653     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1654   }
1655
1656   // pt associated (fixed for all vertices/psi/centralities)
1657   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1658     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1659     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1660     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1661     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1662   }
1663
1664   // ============================================================================================
1665   // the same for event mixing
1666   AliTHn *fHistPMix = bfMix->GetHistNp();
1667   AliTHn *fHistNMix = bfMix->GetHistNn();
1668   AliTHn *fHistPNMix = bfMix->GetHistNpn();
1669   AliTHn *fHistNPMix = bfMix->GetHistNnp();
1670   AliTHn *fHistPPMix = bfMix->GetHistNpp();
1671   AliTHn *fHistNNMix = bfMix->GetHistNnn();
1672
1673   // pt trigger (fixed for all vertices/psi/centralities)
1674   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1675     fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1676     fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1677     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1678     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1679     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1680     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1681   }
1682
1683   // pt associated (fixed for all vertices/psi/centralities)
1684   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1685     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1686     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1687     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1688     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1689   }
1690
1691   // ============================================================================================
1692   // ranges (in bins) for vertexZ and centrality (psi)
1693   Int_t binPsiMin    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1694   Int_t binPsiMax    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1695   Int_t binVertexMin = 0;
1696   Int_t binVertexMax = 0;
1697   if(fVertexBinning){
1698     binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1699     binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1700   }  
1701
1702   TH2D* h1 = NULL;
1703   TH2D* h2 = NULL;
1704   TH2D* h3 = NULL;
1705   TH2D* h4 = NULL;
1706
1707   // loop over all bins
1708   for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1709     for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1710   
1711       cout<<"In the balance function (2D->1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin)  "<<endl;
1712
1713       // Psi_2
1714       fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1715       fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1716       fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1717       fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1718       fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1719       fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1720   
1721       // Vz
1722       if(fVertexBinning){
1723         fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1724         fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1725         fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1726         fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1727         fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1728         fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1729       }
1730
1731       // ============================================================================================
1732       // the same for event mixing
1733       
1734       // Psi_2
1735       fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1736       fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1737       fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1738       fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1739       fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1740       fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1741
1742       // Vz
1743       if(fVertexBinning){
1744         fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1745         fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1746         fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1747         fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1748         fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1749         fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
1750       }
1751       // ============================================================================================
1752
1753       //AliInfo(Form("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)));
1754       
1755       // Project into the wanted space (1st: analysis step, 2nd: axis)
1756       TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1757       TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1758       TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1759       TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1760       TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1761       TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1762
1763       // ============================================================================================
1764       // the same for event mixing
1765       TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1766       TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1767       TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1768       TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1769       // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1770       // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1771       // ============================================================================================
1772       
1773       hTemp1->Sumw2();
1774       hTemp2->Sumw2();
1775       hTemp3->Sumw2();
1776       hTemp4->Sumw2();
1777       hTemp1Mix->Sumw2();
1778       hTemp2Mix->Sumw2();
1779       hTemp3Mix->Sumw2();
1780       hTemp4Mix->Sumw2();
1781       
1782       hTemp1->Scale(1./hTemp5->Integral());
1783       hTemp3->Scale(1./hTemp5->Integral());
1784       hTemp2->Scale(1./hTemp6->Integral());
1785       hTemp4->Scale(1./hTemp6->Integral());
1786
1787       // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1788       Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1789       mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1790       hTemp1Mix->Scale(1./mixedNorm1);
1791       Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1792       mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1793       hTemp2Mix->Scale(1./mixedNorm2);
1794       Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1795       mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1796       hTemp3Mix->Scale(1./mixedNorm3);
1797       Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1798       mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1799       hTemp4Mix->Scale(1./mixedNorm4);
1800
1801       hTemp1->Divide(hTemp1Mix);
1802       hTemp2->Divide(hTemp2Mix);
1803       hTemp3->Divide(hTemp3Mix);
1804       hTemp4->Divide(hTemp4Mix);
1805
1806       // for the first: clone
1807       if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1808         h1 = (TH2D*)hTemp1->Clone();
1809         h2 = (TH2D*)hTemp2->Clone();
1810         h3 = (TH2D*)hTemp3->Clone();
1811         h4 = (TH2D*)hTemp4->Clone();
1812       }
1813       else{  // otherwise: add for averaging
1814         h1->Add(hTemp1);
1815         h2->Add(hTemp2);
1816         h3->Add(hTemp3);
1817         h4->Add(hTemp4);
1818       }
1819
1820     }
1821   }
1822
1823   TH1D *gHistBalanceFunctionHistogram = 0x0;
1824   if((h1)&&(h2)&&(h3)&&(h4)) {
1825
1826     // average over number of bins nbinsVertex * nbinsPsi
1827     h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1828     h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1829     h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1830     h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1831
1832     // now only project on one axis
1833     TH1D *h1DTemp1 = NULL;
1834     TH1D *h1DTemp2 = NULL;
1835     TH1D *h1DTemp3 = NULL;
1836     TH1D *h1DTemp4 = NULL;
1837     if(!bPhi){
1838       h1DTemp1 = (TH1D*)h1->ProjectionX(Form("%s_projX",h1->GetName()));
1839       h1DTemp2 = (TH1D*)h2->ProjectionX(Form("%s_projX",h2->GetName()));
1840       h1DTemp3 = (TH1D*)h3->ProjectionX(Form("%s_projX",h3->GetName()));
1841       h1DTemp4 = (TH1D*)h4->ProjectionX(Form("%s_projX",h4->GetName()));
1842     }
1843     else{
1844       h1DTemp1 = (TH1D*)h1->ProjectionY(Form("%s_projY",h1->GetName()));
1845       h1DTemp2 = (TH1D*)h2->ProjectionY(Form("%s_projY",h2->GetName()));
1846       h1DTemp3 = (TH1D*)h3->ProjectionY(Form("%s_projY",h3->GetName()));
1847       h1DTemp4 = (TH1D*)h4->ProjectionY(Form("%s_projY",h4->GetName()));
1848     }
1849
1850     gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1851     gHistBalanceFunctionHistogram->Reset();
1852     if(!bPhi){
1853       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");  
1854       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");  
1855     }
1856     else{
1857       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1858       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");  
1859     }
1860
1861     h1DTemp1->Add(h1DTemp3,-1.);
1862     h1DTemp2->Add(h1DTemp4,-1.);
1863
1864     gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1865     gHistBalanceFunctionHistogram->Scale(0.5);
1866
1867     //normalize to bin width
1868     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
1869   }
1870
1871   return gHistBalanceFunctionHistogram;
1872 }
1873
1874
1875 //____________________________________________________________________//
1876 TH2D *AliBalancePsi::GetCorrelationFunction(TString type,
1877                                             Double_t psiMin, 
1878                                             Double_t psiMax,
1879                                             Double_t vertexZMin,
1880                                             Double_t vertexZMax,
1881                                             Double_t ptTriggerMin,
1882                                             Double_t ptTriggerMax,
1883                                             Double_t ptAssociatedMin,
1884                                             Double_t ptAssociatedMax,
1885                                             AliBalancePsi *bMixed,
1886                                             Bool_t normToTrig,
1887                                             Double_t normalizationRangePhi) {
1888
1889   // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
1890   // does the division by event mixing inside,
1891   // and averaging over several vertexZ and centrality bins
1892
1893   // security checks
1894   if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
1895     AliError("Only types allowed: PN,NP,PP,NN,ALL");
1896     return NULL;
1897   }
1898   if(!bMixed){
1899     AliError("No Event Mixing AliTHn");
1900     return NULL;
1901   }
1902
1903   TH2D *gHist  = NULL;
1904   TH2D *fSame  = NULL;
1905   TH2D *fMixed = NULL;
1906
1907   // ranges (in bins) for vertexZ and centrality (psi)
1908   Int_t binPsiMin    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin+0.00001);
1909   Int_t binPsiMax    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1910   Int_t binVertexMin = 0;
1911   Int_t binVertexMax = 0;
1912   if(fVertexBinning){
1913     binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin+0.00001);
1914     binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1915   }  
1916   Double_t binPsiLowEdge    = 0.;
1917   Double_t binPsiUpEdge     = 1000.;
1918   Double_t binVertexLowEdge = 0.;
1919   Double_t binVertexUpEdge  = 1000.;
1920
1921   // loop over all bins
1922   for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1923     for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1924
1925       //cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin)  "<<endl;
1926
1927       // determine the bin edges for this bin
1928       binPsiLowEdge    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi) + 0.00001;
1929       binPsiUpEdge     = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi) - 0.00001;
1930       if(fVertexBinning){
1931         binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex) + 0.00001;
1932         binVertexUpEdge  = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex) - 0.00001;
1933       }
1934
1935       // get the 2D histograms for the correct type
1936       if(type=="PN"){
1937         fSame  = GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1938         fMixed = bMixed->GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1939       }
1940       else if(type=="NP"){
1941         fSame  = GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1942         fMixed = bMixed->GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1943       }
1944       else if(type=="PP"){
1945         fSame  = GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1946         fMixed = bMixed->GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1947       }
1948       else if(type=="NN"){
1949         fSame  = GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1950         fMixed = bMixed->GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1951       }
1952       else if(type=="ALL"){
1953         fSame  = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1954         fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1955       }
1956     
1957       if(fMixed && normToTrig && fMixed->Integral()>0){
1958         
1959         // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1960         // do it only on away-side (due to two-track cuts): pi +- pi/6.
1961         Int_t binXmin  = fMixed->GetXaxis()->FindBin(0-10e-5);
1962         Int_t binXmax  = fMixed->GetXaxis()->FindBin(0+10e-5);
1963         Double_t binsX = (Double_t)(binXmax - binXmin + 1);
1964         Int_t binYmin  = fMixed->GetYaxis()->FindBin(TMath::Pi() - normalizationRangePhi);
1965         Int_t binYmax  = fMixed->GetYaxis()->FindBin(TMath::Pi() + normalizationRangePhi - 0.00001);
1966         Double_t binsY = (Double_t)(binYmax - binYmin + 1);
1967         
1968         Double_t mixedNorm = fMixed->Integral(binXmin,binXmax,binYmin,binYmax);
1969         mixedNorm /= binsX * binsY;
1970
1971         // finite bin correction
1972         Double_t binWidthEta = fMixed->GetXaxis()->GetBinWidth(fMixed->GetNbinsX());
1973         Double_t maxEta      = fMixed->GetXaxis()->GetBinUpEdge(fMixed->GetNbinsX());
1974         
1975         Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
1976         //Printf("Finite bin correction: %f", finiteBinCorrection);
1977         mixedNorm /= finiteBinCorrection;
1978         
1979         fMixed->Scale(1./mixedNorm);
1980       }
1981
1982       if(fSame && fMixed){
1983         // then get the correlation function (divide fSame/fmixed)
1984         fSame->Divide(fMixed);
1985         
1986         // averaging with number of triggers:
1987         // average over number of triggers in each sub-bin
1988         Double_t NTrigSubBin = 0;
1989         if(type=="PN" || type=="PP")
1990           NTrigSubBin = (Double_t)(fHistP->Project(0,1)->Integral());
1991         else if(type=="NP" || type=="NN")
1992           NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral());
1993         else if(type=="ALL")
1994           NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
1995         fSame->Scale(NTrigSubBin);
1996         
1997         // for the first: clone
1998         if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
1999           gHist = (TH2D*)fSame->Clone();
2000         }
2001         else{  // otherwise: add for averaging
2002           gHist->Add(fSame);
2003         }
2004       }
2005     }
2006   }
2007
2008   if(gHist){
2009     
2010     // averaging with number of triggers:
2011     // first set to full range and then obtain number of all triggers 
2012     Double_t NTrigAll = 0;
2013     if(type=="PN" || type=="PP"){
2014       fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2015       fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2016       fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2017       NTrigAll = (Double_t)(fHistP->Project(0,1)->Integral());
2018     }
2019     else if(type=="NP" || type=="NN"){
2020       fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2021       fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2022       fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2023       NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral());
2024     }
2025     else if(type=="ALL"){
2026       fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2027       fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2028       fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2029       fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2030       fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2031       fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2032       NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
2033     }
2034     gHist->Scale(1./NTrigAll);
2035    
2036   }
2037   
2038   return gHist;
2039 }
2040
2041
2042 //____________________________________________________________________//
2043 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin, 
2044                                               Double_t psiMax,
2045                                               Double_t vertexZMin,
2046                                               Double_t vertexZMax,
2047                                               Double_t ptTriggerMin,
2048                                               Double_t ptTriggerMax,
2049                                               Double_t ptAssociatedMin,
2050                                               Double_t ptAssociatedMax) {
2051   //Returns the 2D correlation function for (+-) pairs
2052
2053   // security checks
2054   if(psiMin > psiMax-0.00001){
2055     AliError("psiMax <= psiMin");
2056     return NULL;
2057   }
2058   if(vertexZMin > vertexZMax-0.00001){
2059     AliError("vertexZMax <= vertexZMin");
2060     return NULL;
2061   }
2062   if(ptTriggerMin > ptTriggerMax-0.00001){
2063     AliError("ptTriggerMax <= ptTriggerMin");
2064     return NULL;
2065   }
2066   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2067     AliError("ptAssociatedMax <= ptAssociatedMin");
2068     return NULL;
2069   }
2070
2071   // Psi_2: axis 0
2072   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2073   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2074
2075   // Vz
2076   if(fVertexBinning){
2077     //Printf("Cutting on Vz...");
2078     fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2079     fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2080   }
2081
2082   // pt trigger
2083   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2084     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2085     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2086   }
2087
2088   // pt associated
2089   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2090     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2091
2092   //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5); 
2093   //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5); 
2094
2095   //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
2096   //TCanvas *c1 = new TCanvas("c1","");
2097   //c1->cd();
2098   //if(!gHistTest){
2099   //AliError("Projection of fHistP = NULL");
2100   //return gHistTest;
2101   //}
2102   //else{
2103   //gHistTest->DrawCopy("colz");
2104   //}
2105
2106   //0:step, 1: Delta eta, 2: Delta phi
2107   TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2108   if(!gHist){
2109     AliError("Projection of fHistPN = NULL");
2110     return gHist;
2111   }
2112
2113   //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
2114   //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
2115   //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
2116   //AliInfo(Form("Integral (1D): %lf",(Double_t)(fHistP->Project(0,1)->Integral())));
2117   //AliInfo(Form("Integral (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->Integral())));
2118   
2119   //TCanvas *c2 = new TCanvas("c2","");
2120   //c2->cd();
2121   //fHistPN->Project(0,1,2)->DrawCopy("colz");
2122
2123   if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2124     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
2125
2126   //normalize to bin width
2127   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2128     
2129   return gHist;
2130 }
2131
2132 //____________________________________________________________________//
2133 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin, 
2134                                               Double_t psiMax,
2135                                               Double_t vertexZMin,
2136                                               Double_t vertexZMax,
2137                                               Double_t ptTriggerMin,
2138                                               Double_t ptTriggerMax,
2139                                               Double_t ptAssociatedMin,
2140                                               Double_t ptAssociatedMax) {
2141   //Returns the 2D correlation function for (+-) pairs
2142
2143   // security checks
2144   if(psiMin > psiMax-0.00001){
2145     AliError("psiMax <= psiMin");
2146     return NULL;
2147   }
2148   if(vertexZMin > vertexZMax-0.00001){
2149     AliError("vertexZMax <= vertexZMin");
2150     return NULL;
2151   }
2152   if(ptTriggerMin > ptTriggerMax-0.00001){
2153     AliError("ptTriggerMax <= ptTriggerMin");
2154     return NULL;
2155   }
2156   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2157     AliError("ptAssociatedMax <= ptAssociatedMin");
2158     return NULL;
2159   }
2160
2161   // Psi_2: axis 0
2162   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2163   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2164
2165   // Vz
2166   if(fVertexBinning){
2167     fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2168     fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2169   }
2170     
2171   // pt trigger
2172   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2173     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2174     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2175   }
2176
2177   // pt associated
2178   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2179     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2180
2181   //0:step, 1: Delta eta, 2: Delta phi
2182   TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2183   if(!gHist){
2184     AliError("Projection of fHistPN = NULL");
2185     return gHist;
2186   }
2187
2188   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2189   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
2190   if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2191     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
2192
2193   //normalize to bin width
2194   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2195     
2196   return gHist;
2197 }
2198
2199 //____________________________________________________________________//
2200 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin, 
2201                                               Double_t psiMax,
2202                                               Double_t vertexZMin,
2203                                               Double_t vertexZMax,
2204                                               Double_t ptTriggerMin,
2205                                               Double_t ptTriggerMax,
2206                                               Double_t ptAssociatedMin,
2207                                               Double_t ptAssociatedMax) {
2208   //Returns the 2D correlation function for (+-) pairs
2209
2210   // security checks
2211   if(psiMin > psiMax-0.00001){
2212     AliError("psiMax <= psiMin");
2213     return NULL;
2214   }
2215   if(vertexZMin > vertexZMax-0.00001){
2216     AliError("vertexZMax <= vertexZMin");
2217     return NULL;
2218   }
2219   if(ptTriggerMin > ptTriggerMax-0.00001){
2220     AliError("ptTriggerMax <= ptTriggerMin");
2221     return NULL;
2222   }
2223   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2224     AliError("ptAssociatedMax <= ptAssociatedMin");
2225     return NULL;
2226   }
2227
2228   // Psi_2: axis 0
2229   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2230   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2231
2232   // Vz
2233   if(fVertexBinning){
2234     fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2235     fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2236   }
2237
2238   // pt trigger
2239   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2240     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2241     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2242   }
2243
2244   // pt associated
2245   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2246     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2247       
2248   //0:step, 1: Delta eta, 2: Delta phi
2249   TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2250   if(!gHist){
2251     AliError("Projection of fHistPN = NULL");
2252     return gHist;
2253   }
2254
2255   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
2256   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
2257   if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2258     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
2259
2260   //normalize to bin width
2261   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2262   
2263   return gHist;
2264 }
2265
2266 //____________________________________________________________________//
2267 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin, 
2268                                               Double_t psiMax,
2269                                               Double_t vertexZMin,
2270                                               Double_t vertexZMax,
2271                                               Double_t ptTriggerMin,
2272                                               Double_t ptTriggerMax,
2273                                               Double_t ptAssociatedMin,
2274                                               Double_t ptAssociatedMax) {
2275   //Returns the 2D correlation function for (+-) pairs
2276
2277   // security checks
2278   if(psiMin > psiMax-0.00001){
2279     AliError("psiMax <= psiMin");
2280     return NULL;
2281   }
2282   if(vertexZMin > vertexZMax-0.00001){
2283     AliError("vertexZMax <= vertexZMin");
2284     return NULL;
2285   }
2286   if(ptTriggerMin > ptTriggerMax-0.00001){
2287     AliError("ptTriggerMax <= ptTriggerMin");
2288     return NULL;
2289   }
2290   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2291     AliError("ptAssociatedMax <= ptAssociatedMin");
2292     return NULL;
2293   }
2294
2295   // Psi_2: axis 0
2296   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2297   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2298
2299   // Vz
2300   if(fVertexBinning){
2301     fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2302     fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2303   }
2304
2305   // pt trigger
2306   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2307     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2308     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2309   }
2310
2311   // pt associated
2312   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2313     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2314     
2315   //0:step, 1: Delta eta, 2: Delta phi
2316   TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2317   if(!gHist){
2318     AliError("Projection of fHistPN = NULL");
2319     return gHist;
2320   }
2321
2322   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2323   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
2324   if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2325     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
2326
2327   //normalize to bin width
2328   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2329     
2330   return gHist;
2331 }
2332
2333 //____________________________________________________________________//
2334 TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin, 
2335                                                              Double_t psiMax,
2336                                                              Double_t vertexZMin,
2337                                                              Double_t vertexZMax,
2338                                                              Double_t ptTriggerMin,
2339                                                              Double_t ptTriggerMax,
2340                                                              Double_t ptAssociatedMin,
2341                                                              Double_t ptAssociatedMax) {
2342   //Returns the 2D correlation function for the sum of all charge combination pairs
2343
2344   // security checks
2345   if(psiMin > psiMax-0.00001){
2346     AliError("psiMax <= psiMin");
2347     return NULL;
2348   }
2349   if(vertexZMin > vertexZMax-0.00001){
2350     AliError("vertexZMax <= vertexZMin");
2351     return NULL;
2352   }
2353   if(ptTriggerMin > ptTriggerMax-0.00001){
2354     AliError("ptTriggerMax <= ptTriggerMin");
2355     return NULL;
2356   }
2357   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2358     AliError("ptAssociatedMax <= ptAssociatedMin");
2359     return NULL;
2360   }
2361
2362   // Psi_2: axis 0
2363   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2364   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2365   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2366   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2367   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2368   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); 
2369  
2370   // Vz
2371   if(fVertexBinning){
2372     fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2373     fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2374     fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2375     fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2376     fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2377     fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); 
2378   }
2379
2380   // pt trigger
2381   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2382     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2383     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2384     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2385     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2386     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2387     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2388   }
2389
2390   // pt associated
2391   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)){
2392     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2393     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2394     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2395     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2396   }
2397
2398   //0:step, 1: Delta eta, 2: Delta phi
2399   TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2400   if(!gHistNN){
2401     AliError("Projection of fHistNN = NULL");
2402     return gHistNN;
2403   }
2404   TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2405   if(!gHistPP){
2406     AliError("Projection of fHistPP = NULL");
2407     return gHistPP;
2408   }
2409   TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2410   if(!gHistNP){
2411     AliError("Projection of fHistNP = NULL");
2412     return gHistNP;
2413   }
2414   TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2415   if(!gHistPN){
2416     AliError("Projection of fHistPN = NULL");
2417     return gHistPN;
2418   }
2419
2420   // sum all 2 particle histograms
2421   gHistNN->Add(gHistPP);
2422   gHistNN->Add(gHistNP);
2423   gHistNN->Add(gHistPN);
2424
2425   // divide by sum of + and - triggers
2426   if((Double_t)(fHistN->Project(0,1)->Integral())>0 && (Double_t)(fHistP->Project(0,1)->Integral())>0)
2427     gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral()));
2428
2429   //normalize to bin width
2430   gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
2431   
2432   return gHistNN;
2433 }
2434
2435 //____________________________________________________________________//
2436 Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
2437                                            Double_t &mean, Double_t &meanError,
2438                                            Double_t &sigma, Double_t &sigmaError,
2439                                            Double_t &skewness, Double_t &skewnessError,
2440                                            Double_t &kurtosis, Double_t &kurtosisError) {
2441   //
2442   // helper method to calculate the moments and errors of a TH1D anlytically
2443   // fVariable = 1 for Delta eta (moments in full range)
2444   //           = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2)
2445   //
2446   
2447   Bool_t success = kFALSE;
2448   mean          = -1.;
2449   meanError     = -1.;
2450   sigma         = -1.;
2451   sigmaError    = -1.;
2452   skewness      = -1.;
2453   skewnessError = -1.;
2454   kurtosis      = -1.;
2455   kurtosisError = -1.;
2456
2457   if(gHist){
2458
2459     // ----------------------------------------------------------------------
2460     // basic parameters of histogram
2461
2462     Int_t fNumberOfBins = gHist->GetNbinsX();
2463     //    Int_t fBinWidth     = gHist->GetBinWidth(1); // assume equal binning
2464
2465
2466     // ----------------------------------------------------------------------
2467     // ZYAM (for partially negative distributions)
2468     // --> we subtract always the minimum value
2469     Double_t zeroYield    = 0.;
2470     Double_t zeroYieldCur = -FLT_MAX;
2471     if(kUseZYAM){
2472       for(Int_t iMin = 0; iMin<2; iMin++){
2473         zeroYieldCur = gHist->GetMinimum(zeroYieldCur);
2474         zeroYield   += zeroYieldCur;
2475       }
2476       zeroYield /= 2.;
2477       //zeroYield = gHist->GetMinimum();
2478     }
2479     // ----------------------------------------------------------------------
2480     // first calculate the mean
2481
2482     Double_t fWeightedAverage   = 0.;
2483     Double_t fNormalization     = 0.;
2484
2485     for(Int_t i = 1; i <= fNumberOfBins; i++) {
2486
2487       // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2488       if(fVariable == 2 && 
2489          (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2490           gHist->GetBinCenter(i) > TMath::Pi()/2)){
2491         continue;
2492       }
2493
2494       fWeightedAverage   += (gHist->GetBinContent(i)-zeroYield) * gHist->GetBinCenter(i);
2495       fNormalization     += (gHist->GetBinContent(i)-zeroYield);
2496     }  
2497     
2498     mean = fWeightedAverage / fNormalization;
2499
2500     // ----------------------------------------------------------------------
2501     // then calculate the higher moments
2502
2503     Double_t fMu  = 0.;
2504     Double_t fMu2 = 0.;
2505     Double_t fMu3 = 0.;
2506     Double_t fMu4 = 0.;
2507     Double_t fMu5 = 0.;
2508     Double_t fMu6 = 0.;
2509     Double_t fMu7 = 0.;
2510     Double_t fMu8 = 0.;
2511
2512     for(Int_t i = 1; i <= fNumberOfBins; i++) {
2513
2514       // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2515       if(fVariable == 2 && 
2516          (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2517           gHist->GetBinCenter(i) > TMath::Pi()/2)){
2518         continue;
2519       }
2520
2521       fMu  += (gHist->GetBinContent(i)-zeroYield) * (gHist->GetBinCenter(i) - mean);
2522       fMu2 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
2523       fMu3 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
2524       fMu4 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
2525       fMu5 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),5);
2526       fMu6 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),6);
2527       fMu7 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),7);
2528       fMu8 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),8);
2529     }
2530
2531     // normalize to bin entries!
2532     fMu  /= fNormalization;    
2533     fMu2 /= fNormalization;    
2534     fMu3 /= fNormalization;    
2535     fMu4 /= fNormalization;    
2536     fMu5 /= fNormalization;    
2537     fMu6 /= fNormalization;    
2538     fMu7 /= fNormalization;    
2539     fMu8 /= fNormalization;    
2540
2541     sigma    = TMath::Sqrt(fMu2);
2542     skewness = fMu3 / TMath::Power(sigma,3);
2543     kurtosis = fMu4 / TMath::Power(sigma,4) - 3;
2544
2545     // normalize with sigma^r
2546     fMu  /= TMath::Power(sigma,1);    
2547     fMu2 /= TMath::Power(sigma,2);    
2548     fMu3 /= TMath::Power(sigma,3);    
2549     fMu4 /= TMath::Power(sigma,4);    
2550     fMu5 /= TMath::Power(sigma,5);    
2551     fMu6 /= TMath::Power(sigma,6);    
2552     fMu7 /= TMath::Power(sigma,7);    
2553     fMu8 /= TMath::Power(sigma,8);    
2554   
2555     // ----------------------------------------------------------------------
2556     // then calculate the higher moment errors
2557     // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
2558     // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl;
2559     // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl;
2560
2561     Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError
2562
2563     // // assuming normal distribution (as done in ROOT)
2564     // meanError     = sigma / TMath::Sqrt(normError); 
2565     // sigmaError    = TMath::Sqrt(sigma*sigma/(2*normError));
2566     // skewnessError = TMath::Sqrt(6./(normError));
2567     // kurtosisError = TMath::Sqrt(24./(normError));
2568
2569     // use delta theorem paper (Luo - arXiv:1109.0593v1)
2570     Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError));
2571     Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError);
2572     Double_t Lambda33 = TMath::Abs((-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError);
2573     //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
2574     //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
2575     //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
2576
2577     // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
2578     // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
2579
2580     if (TMath::Sqrt(normError) != 0){
2581       meanError        = sigma / TMath::Sqrt(normError); 
2582       sigmaError       = TMath::Sqrt(Lambda11);
2583       skewnessError    = TMath::Sqrt(Lambda22);
2584       kurtosisError    = TMath::Sqrt(Lambda33);
2585
2586       success = kTRUE;    
2587           
2588     }
2589     else success = kFALSE;
2590   
2591   }
2592   return success;
2593 }
2594
2595
2596 //____________________________________________________________________//
2597 Float_t AliBalancePsi::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign) { 
2598   //
2599   // calculates dphistar
2600   //
2601   Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2602   
2603   static const Double_t kPi = TMath::Pi();
2604   
2605   // circularity
2606 //   if (dphistar > 2 * kPi)
2607 //     dphistar -= 2 * kPi;
2608 //   if (dphistar < -2 * kPi)
2609 //     dphistar += 2 * kPi;
2610   
2611   if (dphistar > kPi)
2612     dphistar = kPi * 2 - dphistar;
2613   if (dphistar < -kPi)
2614     dphistar = -kPi * 2 - dphistar;
2615   if (dphistar > kPi) // might look funny but is needed
2616     dphistar = kPi * 2 - dphistar;
2617   
2618   return dphistar;
2619 }
2620
2621 //____________________________________________________________________//
2622 Double_t* AliBalancePsi::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
2623 {
2624   // This method is a copy from AliUEHist::GetBinning
2625   // takes the binning from <configuration> identified by <tag>
2626   // configuration syntax example:
2627   // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
2628   // phi: .....
2629   //
2630   // returns bin edges which have to be deleted by the caller
2631   
2632   TString config(configuration);
2633   TObjArray* lines = config.Tokenize("\n");
2634   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2635   {
2636     TString line(lines->At(i)->GetName());
2637     if (line.BeginsWith(TString(tag) + ":"))
2638     {
2639       line.Remove(0, strlen(tag) + 1);
2640       line.ReplaceAll(" ", "");
2641       TObjArray* binning = line.Tokenize(",");
2642       Double_t* bins = new Double_t[binning->GetEntriesFast()];
2643       for (Int_t j=0; j<binning->GetEntriesFast(); j++)
2644         bins[j] = TString(binning->At(j)->GetName()).Atof();
2645       
2646       nBins = binning->GetEntriesFast() - 1;
2647
2648       delete binning;
2649       delete lines;
2650       return bins;
2651     }
2652   }
2653   
2654   delete lines;
2655   AliFatal(Form("Tag %s not found in %s", tag, configuration));
2656   return 0;
2657 }
2658
2659 //____________________________________________________________________//
2660 Double_t AliBalancePsi::GetFWHM(Int_t gDeltaEtaPhi, TH1D* gHist,
2661                                 Double_t &fwhm_spline, Double_t &fwhmError) {
2662  
2663   // helper method to calculate the fwhm and errors of a TH1D 
2664   if(gHist){
2665     Int_t repeat = 10000;
2666     
2667     if (gDeltaEtaPhi == 1){ 
2668       fwhm_spline  = 0.;
2669       fwhmError    = 0.;
2670       gHist->Smooth(4); 
2671
2672       Double_t xmin = gHist->GetXaxis()->GetXmin();
2673       Double_t xmax = gHist->GetXaxis()->GetXmax();
2674       
2675       //+++++++FWHM TSpline+++++++++++++++++++++++++++//
2676       TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0); 
2677       spline3->SetLineColor(kGreen+2); 
2678       spline3->SetLineWidth(2);
2679       //spline3->Draw("SAME");
2680       
2681       Int_t Nbin = gHist->GetNbinsX();
2682       Int_t bin_max_hist_y = gHist->GetMaximumBin();
2683       Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y);
2684       Double_t max_spline = spline3->Eval(bins_center_x);      
2685       Double_t y_spline = -999;
2686       Double_t y_spline_r = -999;
2687       Double_t x_spline_l = -999;
2688       Double_t x_spline_r = -999;
2689       
2690       for (Int_t i= 0; i < bin_max_hist_y*1000; i++){
2691         y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000));
2692         if (y_spline > max_spline/2){
2693           x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000);
2694           break;
2695         }
2696       } 
2697       for (Int_t j= 0; j < bin_max_hist_y*1000; j++){
2698         y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000));
2699         if (y_spline_r < max_spline/2){
2700           x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000);
2701           break;
2702         }
2703       }   
2704       fwhm_spline = x_spline_r - x_spline_l;
2705       //+++++++FWHM TSpline++++++++++++++++++++++++//
2706     
2707       //+++++++Error Calculation SPLINE++++++++++++//
2708       Double_t dmeans,dsigmas;
2709       TSpline3 *spline1;  
2710       Int_t bin_max_hist_y1;
2711       Double_t bins_center_x1;
2712       Double_t max_spline1;
2713       Double_t y_spline1 = - 999.; 
2714       Double_t y_spline_r1 = - 999.; 
2715       Double_t x_spline_l1 = -999.; 
2716       Double_t x_spline_r1 = -999.;
2717       Double_t fwhm_spline1 = 0.;
2718       Double_t fwhm_T = 0.;  
2719  
2720       TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax);
2721       TRandom3 *rgauss = new TRandom3(0);      
2722       TH1F *hists = new TH1F("hists","hists",1000,1.,3);
2723       
2724       for (Int_t f=0; f<repeat; f++){ //100
2725         for (Int_t h=0; h<Nbin+1; h++){   
2726           dmeans = gHist->GetBinContent(h);
2727           dsigmas = gHist->GetBinError(h);
2728           Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas);
2729           hEmpty2->SetBinContent(h,rgauss_point);
2730           
2731           //++++++++++++//  
2732           hEmpty2->Smooth(4);   
2733           //++++++++++++//
2734         }           
2735         spline1 = new TSpline3(hEmpty2,"b2e2",0,0); 
2736         spline1->SetLineColor(kMagenta+1); 
2737         spline1->SetLineWidth(1);
2738         //spline1->Draw("SAME");  
2739         
2740         bin_max_hist_y1 = hEmpty2->GetMaximumBin();
2741         bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1);
2742         max_spline1 = spline1->Eval(bins_center_x1);
2743         
2744         for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){
2745           y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000));
2746           if (y_spline1 > max_spline1/2){
2747             x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000);
2748             break;
2749           }
2750         } 
2751         for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){
2752           y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000));
2753           if (y_spline_r1 < max_spline1/2){
2754             x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000);
2755             break;
2756           }
2757         }     
2758         
2759         fwhm_spline1 = x_spline_r1 - x_spline_l1;
2760         hists->Fill(fwhm_spline1);   
2761         fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1); 
2762       }
2763      
2764       fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1)));
2765       //+++++++Error Calculation SPLINE+++++++++++//
2766     }
2767     else{  
2768       fwhm_spline  = 0.;
2769       fwhmError    = 0.;
2770       gHist->Smooth(4);
2771       
2772       Double_t xmin = gHist->GetXaxis()->GetXmin();
2773       Double_t xmax = gHist->GetXaxis()->GetXmax();
2774       
2775       //+++++++FWHM TSpline+++++++++++++++++++++++++++//
2776       TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0); 
2777       spline3->SetLineColor(kGreen+2); 
2778       spline3->SetLineWidth(2);
2779       //spline3->Draw("SAME");  
2780       
2781       Int_t Nbin = gHist->GetNbinsX();
2782       Int_t bin_max_hist_y = gHist->GetMaximumBin();
2783       Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y);
2784       Double_t max_spline = spline3->Eval(bins_center_x);
2785       
2786       Int_t bin_min_hist_y = gHist->GetMinimumBin();
2787       Double_t bins_center_xmin = gHist->GetBinCenter(bin_min_hist_y);
2788       Double_t min_spline = spline3->Eval(bins_center_xmin);
2789       
2790       Double_t y_spline = -999.;
2791       Double_t y_spline_r = -999.;
2792       Double_t x_spline_l = -999.;
2793       Double_t x_spline_r = -999.;
2794       
2795       for (Int_t i= 0; i < bin_max_hist_y*1000; i++){
2796         y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000));
2797         if (y_spline > (max_spline+min_spline)/2){
2798           x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000);
2799           break;
2800         }
2801       } 
2802       for (Int_t j= 0; j < bin_max_hist_y*1000; j++){
2803         y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000));
2804         if (y_spline_r < (max_spline+min_spline)/2){
2805           x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000);
2806           break;
2807         }
2808       }   
2809       fwhm_spline = x_spline_r - x_spline_l;
2810       //+++++++FWHM TSpline++++++++++++++++++++++++//
2811       
2812       //+++++++Error Calculation SPLINE++++++++++++//
2813       Double_t dmeans,dsigmas;
2814       TSpline3 *spline1;  
2815       Int_t bin_max_hist_y1;
2816       Double_t bins_center_x1;
2817       Double_t max_spline1;
2818       Double_t y_spline1 = -999.;
2819       Double_t y_spline_r1 = -999.;
2820       Double_t x_spline_l1 = -999.;
2821       Double_t x_spline_r1 = -999.;
2822       Double_t fwhm_spline1 = 0.;
2823       Double_t fwhm_T = 0.;   
2824
2825       TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax);
2826       TRandom3 *rgauss = new TRandom3(0);      
2827       Int_t bin_min_hist_y1; 
2828       Double_t bins_center_xmin1,min_spline1;
2829       
2830       for (Int_t f=0; f<repeat; f++){ //100
2831         for (Int_t h=0; h<Nbin+1; h++){   
2832           dmeans = gHist->GetBinContent(h);
2833           dsigmas = gHist->GetBinError(h);
2834           Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas);
2835           hEmpty2->SetBinContent(h,rgauss_point);
2836           
2837           //++++++++++++//  
2838           hEmpty2->Smooth(4);   
2839           //++++++++++++//
2840         }           
2841         spline1 = new TSpline3(hEmpty2,"b2e2",0,0); 
2842         spline1->SetLineColor(kMagenta+1); 
2843         spline1->SetLineWidth(1);
2844         //spline1->Draw("SAME");  
2845         
2846         bin_max_hist_y1 = hEmpty2->GetMaximumBin();
2847         bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1);
2848         max_spline1 = spline1->Eval(bins_center_x1);
2849         
2850         bin_min_hist_y1 = hEmpty2->GetMinimumBin();
2851         bins_center_xmin1 = hEmpty2->GetBinCenter(bin_min_hist_y1);
2852         min_spline1 = spline1->Eval(bins_center_xmin1);
2853         
2854         for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){
2855           y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000));
2856           if (y_spline1 > (max_spline1+min_spline1)/2){
2857             x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000);
2858             break;
2859           }
2860         } 
2861         for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){
2862           y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000));
2863           if (y_spline_r1 < (max_spline1+min_spline1)/2){
2864             x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000);
2865             break;
2866           }
2867         }     
2868         
2869         fwhm_spline1 = x_spline_r1 - x_spline_l1; 
2870         fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1);    
2871       }           
2872       fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1)));
2873     }
2874   }
2875   return fwhm_spline;
2876 }