PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx

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