PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx

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