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