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->Integral());
 785     hTemp2->Add(hTemp4,-1.);
 786     hTemp2->Scale(1./hTemp6->Integral());
 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->Integral());
1144       hTemp3->Scale(1./hTemp5->Integral());
1145       hTemp2->Scale(1./hTemp6->Integral());
1146       hTemp4->Scale(1./hTemp6->Integral());
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->Integral());
1152       hTemp3Mix->Scale(1./hTemp5Mix->Integral());
1153       hTemp2Mix->Scale(1./hTemp6Mix->Integral());
1154       hTemp4Mix->Scale(1./hTemp6Mix->Integral());
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->Integral());
1328     hTemp2->Add(hTemp4,-1.);
1329     hTemp2->Scale(1./hTemp6->Integral());
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->Integral());
1534       hTemp3->Scale(1./hTemp5->Integral());
1535       hTemp2->Scale(1./hTemp6->Integral());
1536       hTemp4->Scale(1./hTemp6->Integral());
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->Integral());
1793       hTemp3->Scale(1./hTemp5->Integral());
1794       hTemp2->Scale(1./hTemp6->Integral());
1795       hTemp4->Scale(1./hTemp6->Integral());
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                                             Bool_t normToTrig,
1897                                             TH2D *hVertexCentrality) {
1898
1899   // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
1900   // does the division by event mixing inside,
1901   // and averaging over several vertexZ and centrality bins
1902
1903   // security checks
1904   if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
1905     AliError("Only types allowed: PN,NP,PP,NN,ALL");
1906     return NULL;
1907   }
1908   if(!bMixed){
1909     AliError("No Event Mixing AliTHn");
1910     return NULL;
1911   }
1912   if(normToTrig && !hVertexCentrality){
1913     AliError("Per-trigger yield option chosen, but QA histogram for event statistics not available");
1914     return NULL;
1915   }
1916
1917   TH2D *gHist  = NULL;
1918   TH2D *fSame  = NULL;
1919   TH2D *fMixed = NULL;
1920
1921   // ranges (in bins) for vertexZ and centrality (psi)
1922   Int_t binPsiMin    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin+0.00001);
1923   Int_t binPsiMax    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1924   Int_t binVertexMin = 0;
1925   Int_t binVertexMax = 0;
1926   if(fVertexBinning){
1927     binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin+0.00001);
1928     binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1929   }
1930   Double_t binPsiLowEdge    = 0.;
1931   Double_t binPsiUpEdge     = 1000.;
1932   Double_t binVertexLowEdge = 0.;
1933   Double_t binVertexUpEdge  = 1000.;
1934
1935   // loop over all bins
1936   for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1937     for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1938
1939       cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin)  "<<endl;
1940
1941       // determine the bin edges for this bin
1942       binPsiLowEdge    = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi) + 0.00001;
1943       binPsiUpEdge     = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi) - 0.00001;
1944       if(fVertexBinning){
1945         binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex) + 0.00001;
1946         binVertexUpEdge  = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex) - 0.00001;
1947       }
1948
1949       // get the 2D histograms for the correct type
1950       if(type=="PN"){
1951         fSame  = GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1952         fMixed = bMixed->GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1953       }
1954       else if(type=="NP"){
1955         fSame  = GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1956         fMixed = bMixed->GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1957       }
1958       else if(type=="PP"){
1959         fSame  = GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1960         fMixed = bMixed->GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1961       }
1962       else if(type=="NN"){
1963         fSame  = GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1964         fMixed = bMixed->GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1965       }
1966       else if(type=="ALL"){
1967         fSame  = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1968         fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1969       }
1970
1971       if(fMixed && normToTrig && fMixed->Integral()>0){
1972
1973         // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1974         // do it only on away-side (due to two-track cuts): pi +- pi/3.
1975
1976         Int_t binXmin  = fMixed->GetXaxis()->FindBin(0-10e-5);
1977         Int_t binXmax  = fMixed->GetXaxis()->FindBin(0+10e-5);
1978         Double_t binsX = (Double_t)(binXmax - binXmin + 1);
1979         Int_t binYmin  = fMixed->GetYaxis()->FindBin(2./3.*TMath::Pi());
1980         Int_t binYmax  = fMixed->GetYaxis()->FindBin(4./3.*TMath::Pi());
1981         Double_t binsY = (Double_t)(binYmax - binYmin + 1);
1982
1983         Double_t mixedNorm = fMixed->Integral(binXmin,binXmax,binYmin,binYmax);
1984         mixedNorm /= binsX * binsY;
1985
1986         // finite bin correction
1987         Double_t binWidthEta = fMixed->GetXaxis()->GetBinWidth(fMixed->GetNbinsX());
1988         Double_t maxEta      = fMixed->GetXaxis()->GetBinUpEdge(fMixed->GetNbinsX());
1989
1990         Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
1991         //Printf("Finite bin correction: %f", finiteBinCorrection);
1992         mixedNorm /= finiteBinCorrection;
1993
1994         fMixed->Scale(1./mixedNorm);
1995       }
1996
1997       if(fSame && fMixed){
1998         // then get the correlation function (divide fSame/fmixed)
1999         fSame->Divide(fMixed);
2000
2001         // // averaging with number of triggers:
2002         // // average over number of triggers in each sub-bin
2003         // Double_t NTrigSubBin = 0;
2004         // if(type=="PN" || type=="PP")
2005         //   NTrigSubBin = (Double_t)(fHistP->Project(0,1)->Integral());
2006         // else if(type=="NP" || type=="NN")
2007         //   NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral());
2008         // fSame->Scale(NTrigSubBin);
2009
2010         // averaging with number of events:
2011         // average over number of events in each sub-bin
2012         Int_t binStatsVertexLowEdge = hVertexCentrality->GetXaxis()->FindBin(binVertexLowEdge + 0.00001) + 0.00001;
2013         Int_t binStatsVertexUpEdge  = hVertexCentrality->GetXaxis()->FindBin(binVertexUpEdge - 0.00001) - 0.00001;
2014         Int_t binStatsPsiLowEdge    = hVertexCentrality->GetYaxis()->FindBin(binPsiLowEdge + 0.00001) + 0.00001;
2015         Int_t binStatsPsiUpEdge     = hVertexCentrality->GetYaxis()->FindBin(binPsiUpEdge - 0.00001) - 0.00001;
2016
2017         Double_t NEventsSubBin = (Double_t)hVertexCentrality->Integral(binStatsVertexLowEdge,binStatsVertexUpEdge,binStatsPsiLowEdge,binStatsPsiUpEdge);
2018
2019         Printf("Averaging from %d < z < %d and %d < cent < %d ",binStatsVertexLowEdge,binStatsVertexUpEdge,binStatsPsiLowEdge,binStatsPsiUpEdge);
2020         Printf("Averaging from %.2f < z < %.2f and %.2f < cent < %.2f --> %.2f ",binVertexLowEdge,binVertexUpEdge,binPsiLowEdge,binPsiUpEdge,NEventsSubBin);
2021         fSame->Scale(NEventsSubBin);
2022
2023         // OLD and NEW averaging:
2024         // for the first: clone
2025         if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
2026           gHist = (TH2D*)fSame->Clone();
2027         }
2028         else{  // otherwise: add for averaging
2029           gHist->Add(fSame);
2030         }
2031       }
2032     }
2033   }
2034
2035   if(gHist){
2036
2037     //// OLD averaging:
2038     //// average over number of bins nbinsVertex * nbinsPsi
2039     // gHist->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
2040
2041     // // averaging with number of triggers:
2042     // // first set to full range and then obtain number of all triggers
2043     // Double_t NTrigAll = 0;
2044     // if(type=="PN" || type=="PP"){
2045     //   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2046     //   fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2047     //   fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2048     //   NTrigAll = (Double_t)(fHistP->Project(0,1)->Integral());
2049     // }
2050     // else if(type=="NP" || type=="NN"){
2051     //   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2052     //   fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2053     //   fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2054     //   NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral());
2055     // }
2056     // gHist->Scale(1./NTrigAll);
2057
2058     // averaging with number of events:
2059     // first set to full range and then obtain number of all events
2060     Int_t binStatsAllVertexLowEdge = hVertexCentrality->GetXaxis()->FindBin(vertexZMin + 0.00001) + 0.00001;
2061     Int_t binStatsAllVertexUpEdge  = hVertexCentrality->GetXaxis()->FindBin(vertexZMax - 0.00001) - 0.00001;
2062     Int_t binStatsAllPsiLowEdge    = hVertexCentrality->GetYaxis()->FindBin(psiMin + 0.00001) + 0.00001;
2063     Int_t binStatsAllPsiUpEdge     = hVertexCentrality->GetYaxis()->FindBin(psiMax - 0.00001) - 0.00001;
2064
2065     Double_t NEventsAll = (Double_t)hVertexCentrality->Integral(binStatsAllVertexLowEdge,binStatsAllVertexUpEdge,binStatsAllPsiLowEdge,binStatsAllPsiUpEdge);
2066
2067     //Printf("All from %d < z < %d and %d < cent < %d ",binStatsAllVertexLowEdge,binStatsAllVertexUpEdge,binStatsAllPsiLowEdge,binStatsAllPsiUpEdge);
2068     //Printf("All from %.2f < z < %.2f and %.2f < cent < %.2f --> %.2f ",vertexZMin,vertexZMax,psiMin,psiMax,NEventsAll);
2069
2070     gHist->Scale(1./NEventsAll);
2071
2072   }
2073
2074   return gHist;
2075 }
2076
2077
2078 //____________________________________________________________________//
2079 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
2080                                               Double_t psiMax,
2081                                               Double_t vertexZMin,
2082                                               Double_t vertexZMax,
2083                                               Double_t ptTriggerMin,
2084                                               Double_t ptTriggerMax,
2085                                               Double_t ptAssociatedMin,
2086                                               Double_t ptAssociatedMax) {
2087   //Returns the 2D correlation function for (+-) pairs
2088
2089   // security checks
2090   if(psiMin > psiMax-0.00001){
2091     AliError("psiMax <= psiMin");
2092     return NULL;
2093   }
2094   if(vertexZMin > vertexZMax-0.00001){
2095     AliError("vertexZMax <= vertexZMin");
2096     return NULL;
2097   }
2098   if(ptTriggerMin > ptTriggerMax-0.00001){
2099     AliError("ptTriggerMax <= ptTriggerMin");
2100     return NULL;
2101   }
2102   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2103     AliError("ptAssociatedMax <= ptAssociatedMin");
2104     return NULL;
2105   }
2106
2107   // Psi_2: axis 0
2108   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2109   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2110
2111   // Vz
2112   if(fVertexBinning){
2113     //Printf("Cutting on Vz...");
2114     fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2115     fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2116   }
2117
2118   // pt trigger
2119   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2120     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2121     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2122   }
2123
2124   // pt associated
2125   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2126     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2127
2128   //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
2129   //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
2130
2131   //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
2132   //TCanvas *c1 = new TCanvas("c1","");
2133   //c1->cd();
2134   //if(!gHistTest){
2135   //AliError("Projection of fHistP = NULL");
2136   //return gHistTest;
2137   //}
2138   //else{
2139   //gHistTest->DrawCopy("colz");
2140   //}
2141
2142   //0:step, 1: Delta eta, 2: Delta phi
2143   TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2144   if(!gHist){
2145     AliError("Projection of fHistPN = NULL");
2146     return gHist;
2147   }
2148
2149   //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
2150   //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
2151   //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
2152   //AliInfo(Form("Integral (1D): %lf",(Double_t)(fHistP->Project(0,1)->Integral())));
2153   //AliInfo(Form("Integral (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->Integral())));
2154
2155   //TCanvas *c2 = new TCanvas("c2","");
2156   //c2->cd();
2157   //fHistPN->Project(0,1,2)->DrawCopy("colz");
2158
2159   if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2160     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
2161
2162   //normalize to bin width
2163   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2164
2165   return gHist;
2166 }
2167
2168 //____________________________________________________________________//
2169 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
2170                                               Double_t psiMax,
2171                                               Double_t vertexZMin,
2172                                               Double_t vertexZMax,
2173                                               Double_t ptTriggerMin,
2174                                               Double_t ptTriggerMax,
2175                                               Double_t ptAssociatedMin,
2176                                               Double_t ptAssociatedMax) {
2177   //Returns the 2D correlation function for (+-) pairs
2178
2179   // security checks
2180   if(psiMin > psiMax-0.00001){
2181     AliError("psiMax <= psiMin");
2182     return NULL;
2183   }
2184   if(vertexZMin > vertexZMax-0.00001){
2185     AliError("vertexZMax <= vertexZMin");
2186     return NULL;
2187   }
2188   if(ptTriggerMin > ptTriggerMax-0.00001){
2189     AliError("ptTriggerMax <= ptTriggerMin");
2190     return NULL;
2191   }
2192   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2193     AliError("ptAssociatedMax <= ptAssociatedMin");
2194     return NULL;
2195   }
2196
2197   // Psi_2: axis 0
2198   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2199   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2200
2201   // Vz
2202   if(fVertexBinning){
2203     fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2204     fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2205   }
2206
2207   // pt trigger
2208   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2209     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2210     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2211   }
2212
2213   // pt associated
2214   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2215     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2216
2217   //0:step, 1: Delta eta, 2: Delta phi
2218   TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2219   if(!gHist){
2220     AliError("Projection of fHistPN = NULL");
2221     return gHist;
2222   }
2223
2224   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2225   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
2226   if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2227     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
2228
2229   //normalize to bin width
2230   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2231
2232   return gHist;
2233 }
2234
2235 //____________________________________________________________________//
2236 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
2237                                               Double_t psiMax,
2238                                               Double_t vertexZMin,
2239                                               Double_t vertexZMax,
2240                                               Double_t ptTriggerMin,
2241                                               Double_t ptTriggerMax,
2242                                               Double_t ptAssociatedMin,
2243                                               Double_t ptAssociatedMax) {
2244   //Returns the 2D correlation function for (+-) pairs
2245
2246   // security checks
2247   if(psiMin > psiMax-0.00001){
2248     AliError("psiMax <= psiMin");
2249     return NULL;
2250   }
2251   if(vertexZMin > vertexZMax-0.00001){
2252     AliError("vertexZMax <= vertexZMin");
2253     return NULL;
2254   }
2255   if(ptTriggerMin > ptTriggerMax-0.00001){
2256     AliError("ptTriggerMax <= ptTriggerMin");
2257     return NULL;
2258   }
2259   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2260     AliError("ptAssociatedMax <= ptAssociatedMin");
2261     return NULL;
2262   }
2263
2264   // Psi_2: axis 0
2265   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2266   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2267
2268   // Vz
2269   if(fVertexBinning){
2270     fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2271     fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2272   }
2273
2274   // pt trigger
2275   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2276     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2277     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2278   }
2279
2280   // pt associated
2281   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2282     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2283
2284   //0:step, 1: Delta eta, 2: Delta phi
2285   TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2286   if(!gHist){
2287     AliError("Projection of fHistPN = NULL");
2288     return gHist;
2289   }
2290
2291   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
2292   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
2293   if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2294     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
2295
2296   //normalize to bin width
2297   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2298
2299   return gHist;
2300 }
2301
2302 //____________________________________________________________________//
2303 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
2304                                               Double_t psiMax,
2305                                               Double_t vertexZMin,
2306                                               Double_t vertexZMax,
2307                                               Double_t ptTriggerMin,
2308                                               Double_t ptTriggerMax,
2309                                               Double_t ptAssociatedMin,
2310                                               Double_t ptAssociatedMax) {
2311   //Returns the 2D correlation function for (+-) pairs
2312
2313   // security checks
2314   if(psiMin > psiMax-0.00001){
2315     AliError("psiMax <= psiMin");
2316     return NULL;
2317   }
2318   if(vertexZMin > vertexZMax-0.00001){
2319     AliError("vertexZMax <= vertexZMin");
2320     return NULL;
2321   }
2322   if(ptTriggerMin > ptTriggerMax-0.00001){
2323     AliError("ptTriggerMax <= ptTriggerMin");
2324     return NULL;
2325   }
2326   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2327     AliError("ptAssociatedMax <= ptAssociatedMin");
2328     return NULL;
2329   }
2330
2331   // Psi_2: axis 0
2332   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2333   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2334
2335   // Vz
2336   if(fVertexBinning){
2337     fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2338     fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2339   }
2340
2341   // pt trigger
2342   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2343     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2344     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2345   }
2346
2347   // pt associated
2348   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2349     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2350
2351   //0:step, 1: Delta eta, 2: Delta phi
2352   TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2353   if(!gHist){
2354     AliError("Projection of fHistPN = NULL");
2355     return gHist;
2356   }
2357
2358   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2359   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
2360   if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2361     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
2362
2363   //normalize to bin width
2364   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2365
2366   return gHist;
2367 }
2368
2369 //____________________________________________________________________//
2370 TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
2371                                                              Double_t psiMax,
2372                                                              Double_t vertexZMin,
2373                                                              Double_t vertexZMax,
2374                                                              Double_t ptTriggerMin,
2375                                                              Double_t ptTriggerMax,
2376                                                              Double_t ptAssociatedMin,
2377                                                              Double_t ptAssociatedMax) {
2378   //Returns the 2D correlation function for the sum of all charge combination pairs
2379
2380   // security checks
2381   if(psiMin > psiMax-0.00001){
2382     AliError("psiMax <= psiMin");
2383     return NULL;
2384   }
2385   if(vertexZMin > vertexZMax-0.00001){
2386     AliError("vertexZMax <= vertexZMin");
2387     return NULL;
2388   }
2389   if(ptTriggerMin > ptTriggerMax-0.00001){
2390     AliError("ptTriggerMax <= ptTriggerMin");
2391     return NULL;
2392   }
2393   if(ptAssociatedMin > ptAssociatedMax-0.00001){
2394     AliError("ptAssociatedMax <= ptAssociatedMin");
2395     return NULL;
2396   }
2397
2398   // Psi_2: axis 0
2399   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2400   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2401   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2402   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2403   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2404   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2405
2406   // Vz
2407   if(fVertexBinning){
2408     fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2409     fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2410     fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2411     fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2412     fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2413     fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2414   }
2415
2416   // pt trigger
2417   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2418     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2419     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2420     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2421     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2422     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2423     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2424   }
2425
2426   // pt associated
2427   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2428     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2429     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2430     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2431     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2432
2433   //0:step, 1: Delta eta, 2: Delta phi
2434   TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2435   if(!gHistNN){
2436     AliError("Projection of fHistNN = NULL");
2437     return gHistNN;
2438   }
2439   TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2440   if(!gHistPP){
2441     AliError("Projection of fHistPP = NULL");
2442     return gHistPP;
2443   }
2444   TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2445   if(!gHistNP){
2446     AliError("Projection of fHistNP = NULL");
2447     return gHistNP;
2448   }
2449   TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2450   if(!gHistPN){
2451     AliError("Projection of fHistPN = NULL");
2452     return gHistPN;
2453   }
2454
2455   // sum all 2 particle histograms
2456   gHistNN->Add(gHistPP);
2457   gHistNN->Add(gHistNP);
2458   gHistNN->Add(gHistPN);
2459
2460   // divide by sum of + and - triggers
2461   if((Double_t)(fHistN->Project(0,1)->Integral())>0 && (Double_t)(fHistP->Project(0,1)->Integral())>0)
2462     gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral() + fHistN->Project(0,1)->Integral()));
2463
2464   //normalize to bin width
2465   gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
2466
2467   return gHistNN;
2468 }
2469
2470 //____________________________________________________________________//
2471
2472 Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
2473                                            Double_t &mean, Double_t &meanError,
2474                                            Double_t &sigma, Double_t &sigmaError,
2475                                            Double_t &skewness, Double_t &skewnessError,
2476                                            Double_t &kurtosis, Double_t &kurtosisError) {
2477   //
2478   // helper method to calculate the moments and errors of a TH1D anlytically
2479   // fVariable = 1 for Delta eta (moments in full range)
2480   //           = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2)
2481   //
2482
2483   Bool_t success = kFALSE;
2484   mean          = -1.;
2485   meanError     = -1.;
2486   sigma         = -1.;
2487   sigmaError    = -1.;
2488   skewness      = -1.;
2489   skewnessError = -1.;
2490   kurtosis      = -1.;
2491   kurtosisError = -1.;
2492
2493   if(gHist){
2494
2495     // ----------------------------------------------------------------------
2496     // basic parameters of histogram
2497
2498     Int_t fNumberOfBins = gHist->GetNbinsX();
2499     //    Int_t fBinWidth     = gHist->GetBinWidth(1); // assume equal binning
2500
2501
2502     // ----------------------------------------------------------------------
2503     // ZYAM (for partially negative distributions)
2504     // --> we subtract always the minimum value
2505     Double_t zeroYield    = 0.;
2506     Double_t zeroYieldCur = -FLT_MAX;
2507     if(kUseZYAM){
2508       for(Int_t iMin = 0; iMin<2; iMin++){
2509         zeroYieldCur = gHist->GetMinimum(zeroYieldCur);
2510         zeroYield   += zeroYieldCur;
2511       }
2512       zeroYield /= 2.;
2513       //zeroYield = gHist->GetMinimum();
2514     }
2515     // ----------------------------------------------------------------------
2516     // first calculate the mean
2517
2518     Double_t fWeightedAverage   = 0.;
2519     Double_t fNormalization     = 0.;
2520
2521     for(Int_t i = 1; i <= fNumberOfBins; i++) {
2522
2523       // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2524       if(fVariable == 2 &&
2525          (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2526           gHist->GetBinCenter(i) > TMath::Pi()/2)){
2527         continue;
2528       }
2529
2530       fWeightedAverage   += (gHist->GetBinContent(i)-zeroYield) * gHist->GetBinCenter(i);
2531       fNormalization     += (gHist->GetBinContent(i)-zeroYield);
2532     }
2533
2534     mean = fWeightedAverage / fNormalization;
2535
2536     // ----------------------------------------------------------------------
2537     // then calculate the higher moments
2538
2539     Double_t fMu  = 0.;
2540     Double_t fMu2 = 0.;
2541     Double_t fMu3 = 0.;
2542     Double_t fMu4 = 0.;
2543     Double_t fMu5 = 0.;
2544     Double_t fMu6 = 0.;
2545     Double_t fMu7 = 0.;
2546     Double_t fMu8 = 0.;
2547
2548     for(Int_t i = 1; i <= fNumberOfBins; i++) {
2549
2550       // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2551       if(fVariable == 2 &&
2552          (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2553           gHist->GetBinCenter(i) > TMath::Pi()/2)){
2554         continue;
2555       }
2556
2557       fMu  += (gHist->GetBinContent(i)-zeroYield) * (gHist->GetBinCenter(i) - mean);
2558       fMu2 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
2559       fMu3 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
2560       fMu4 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
2561       fMu5 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),5);
2562       fMu6 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),6);
2563       fMu7 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),7);
2564       fMu8 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),8);
2565     }
2566
2567     // normalize to bin entries!
2568     fMu  /= fNormalization;
2569     fMu2 /= fNormalization;
2570     fMu3 /= fNormalization;
2571     fMu4 /= fNormalization;
2572     fMu5 /= fNormalization;
2573     fMu6 /= fNormalization;
2574     fMu7 /= fNormalization;
2575     fMu8 /= fNormalization;
2576
2577     sigma    = TMath::Sqrt(fMu2);
2578     skewness = fMu3 / TMath::Power(sigma,3);
2579     kurtosis = fMu4 / TMath::Power(sigma,4) - 3;
2580
2581     // normalize with sigma^r
2582     fMu  /= TMath::Power(sigma,1);
2583     fMu2 /= TMath::Power(sigma,2);
2584     fMu3 /= TMath::Power(sigma,3);
2585     fMu4 /= TMath::Power(sigma,4);
2586     fMu5 /= TMath::Power(sigma,5);
2587     fMu6 /= TMath::Power(sigma,6);
2588     fMu7 /= TMath::Power(sigma,7);
2589     fMu8 /= TMath::Power(sigma,8);
2590
2591     // ----------------------------------------------------------------------
2592     // then calculate the higher moment errors
2593     // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
2594     // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl;
2595     // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl;
2596
2597     Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError
2598
2599     // // assuming normal distribution (as done in ROOT)
2600     // meanError     = sigma / TMath::Sqrt(normError);
2601     // sigmaError    = TMath::Sqrt(sigma*sigma/(2*normError));
2602     // skewnessError = TMath::Sqrt(6./(normError));
2603     // kurtosisError = TMath::Sqrt(24./(normError));
2604
2605     // use delta theorem paper (Luo - arXiv:1109.0593v1)
2606     Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError));
2607     Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError);
2608     Double_t Lambda33 = TMath::Abs((-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError);
2609     //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
2610     //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
2611     //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
2612
2613     // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
2614     // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
2615
2616     if (TMath::Sqrt(normError) != 0){
2617       meanError        = sigma / TMath::Sqrt(normError);
2618       sigmaError       = TMath::Sqrt(Lambda11);
2619       skewnessError    = TMath::Sqrt(Lambda22);
2620       kurtosisError    = TMath::Sqrt(Lambda33);
2621
2622       success = kTRUE;
2623
2624     }
2625     else success = kFALSE;
2626
2627   }
2628   return success;
2629 }
2630
2631
2632 //____________________________________________________________________//
2633 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) {
2634   //
2635   // calculates dphistar
2636   //
2637   Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2638
2639   static const Double_t kPi = TMath::Pi();
2640
2641   // circularity
2642 //   if (dphistar > 2 * kPi)
2643 //     dphistar -= 2 * kPi;
2644 //   if (dphistar < -2 * kPi)
2645 //     dphistar += 2 * kPi;
2646
2647   if (dphistar > kPi)
2648     dphistar = kPi * 2 - dphistar;
2649   if (dphistar < -kPi)
2650     dphistar = -kPi * 2 - dphistar;
2651   if (dphistar > kPi) // might look funny but is needed
2652     dphistar = kPi * 2 - dphistar;
2653
2654   return dphistar;
2655 }
2656
2657 //____________________________________________________________________//
2658 Double_t* AliBalancePsi::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
2659 {
2660   // This method is a copy from AliUEHist::GetBinning
2661   // takes the binning from <configuration> identified by <tag>
2662   // configuration syntax example:
2663   // 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
2664   // phi: .....
2665   //
2666   // returns bin edges which have to be deleted by the caller
2667
2668   TString config(configuration);
2669   TObjArray* lines = config.Tokenize("\n");
2670   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2671   {
2672     TString line(lines->At(i)->GetName());
2673     if (line.BeginsWith(TString(tag) + ":"))
2674     {
2675       line.Remove(0, strlen(tag) + 1);
2676       line.ReplaceAll(" ", "");
2677       TObjArray* binning = line.Tokenize(",");
2678       Double_t* bins = new Double_t[binning->GetEntriesFast()];
2679       for (Int_t j=0; j<binning->GetEntriesFast(); j++)
2680         bins[j] = TString(binning->At(j)->GetName()).Atof();
2681
2682       nBins = binning->GetEntriesFast() - 1;
2683
2684       delete binning;
2685       delete lines;
2686       return bins;
2687     }
2688   }
2689
2690   delete lines;
2691   AliFatal(Form("Tag %s not found in %s", tag, configuration));
2692   return 0;
2693 }