PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx

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