PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx

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