PWGCF/EBYE/BalanceFunctions/AliBalance.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$ */
  15
  16 //-----------------------------------------------------------------
  17 //           Balance Function class
  18 //   This is the class to deal with the Balance Function analysis
  19 //   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
  20 //-----------------------------------------------------------------
  21
  22
  23 //ROOT
  24 #include <Riostream.h>
  25 #include <TMath.h>
  26 #include <TAxis.h>
  27 #include <TH2D.h>
  28 #include <TLorentzVector.h>
  29 #include <TObjArray.h>
  30 #include <TGraphErrors.h>
  31 #include <TString.h>
  32
  33 #include "AliVParticle.h"
  34 #include "AliMCParticle.h"
  35 #include "AliESDtrack.h"
  36 #include "AliAODTrack.h"
  37
  38 #include "AliBalance.h"
  39
  40 ClassImp(AliBalance)
  41
  42 //____________________________________________________________________//
  43 AliBalance::AliBalance() :
  44   TObject(),
  45   bShuffle(kFALSE),
  46   bHBTcut(kFALSE),
  47   bConversionCut(kFALSE),
  48   fAnalysisLevel("ESD"),
  49   fAnalyzedEvents(0) ,
  50   fCentralityId(0) ,
  51   fCentStart(0.),
  52   fCentStop(0.)
  53 {
  54   // Default constructor
  55
  56   for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
  57     if(i == 6) {
  58       fNumberOfBins[i] = 180;
  59       fP1Start[i]      = -360.0;
  60       fP1Stop[i]       = 360.0;
  61       fP2Start[i]      = -360.0;
  62       fP2Stop[i]       = 360.0;
  63       fP2Step[i]       = 0.1;
  64     }
  65     else {
  66       fNumberOfBins[i] = 20;
  67       fP1Start[i]      = -1.0;
  68       fP1Stop[i]       = 1.0;
  69       fP2Start[i]      = 0.0;
  70       fP2Stop[i]       = 2.0;
  71     }
  72     fP2Step[i] = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins[i];
  73     fCentStart = 0.;
  74     fCentStop  = 0.;
  75
  76     fNn[i] = 0.0;
  77     fNp[i] = 0.0;
  78
  79     for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
  80       fNpp[i][j] = .0;
  81       fNnn[i][j] = .0;
  82       fNpn[i][j] = .0;
  83       fNnp[i][j] = .0;
  84       fB[i][j] = 0.0;
  85       ferror[i][j] = 0.0;
  86     }
  87
  88     fHistP[i]  = NULL;
  89     fHistN[i]  = NULL;
  90     fHistPP[i] = NULL;
  91     fHistPN[i] = NULL;
  92     fHistNP[i] = NULL;
  93     fHistNN[i] = NULL;
  94
  95   }
  96 }
  97
  98
  99 //____________________________________________________________________//
 100 AliBalance::AliBalance(const AliBalance& balance):
 101   TObject(balance),
 102   bShuffle(balance.bShuffle),
 103   bHBTcut(balance.bHBTcut),
 104   bConversionCut(balance.bConversionCut),
 105   fAnalysisLevel(balance.fAnalysisLevel),
 106   fAnalyzedEvents(balance.fAnalyzedEvents),
 107   fCentralityId(balance.fCentralityId),
 108   fCentStart(balance.fCentStart),
 109   fCentStop(balance.fCentStop) {
 110   //copy constructor
 111   for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
 112     fNn[i] = balance.fNn[i];
 113     fNp[i] = balance.fNp[i];
 114
 115     fP1Start[i]      = balance.fP1Start[i];
 116     fP1Stop[i]       = balance.fP1Stop[i];
 117     fNumberOfBins[i] = balance.fNumberOfBins[i];
 118     fP2Start[i]      = balance.fP2Start[i];
 119     fP2Stop[i]       = balance.fP2Stop[i];
 120     fP2Step[i]       = balance.fP2Step[i];
 121     fCentStart       = balance.fCentStart;
 122     fCentStop        = balance.fCentStop;
 123
 124     fHistP[i]        = balance.fHistP[i];
 125     fHistN[i]        = balance.fHistN[i];
 126     fHistPN[i]        = balance.fHistPN[i];
 127     fHistNP[i]        = balance.fHistNP[i];
 128     fHistPP[i]        = balance.fHistPP[i];
 129     fHistNN[i]        = balance.fHistNN[i];
 130
 131     for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
 132       fNpp[i][j] = .0;
 133       fNnn[i][j] = .0;
 134       fNpn[i][j] = .0;
 135       fNnp[i][j] = .0;
 136       fB[i][j] = 0.0;
 137       ferror[i][j] = 0.0;
 138     }
 139   }
 140  }
 141
 142
 143 //____________________________________________________________________//
 144 AliBalance::~AliBalance() {
 145   // Destructor
 146
 147   for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
 148
 149     delete fHistP[i];
 150     delete fHistN[i];
 151     delete fHistPN[i];
 152     delete fHistNP[i];
 153     delete fHistPP[i];
 154     delete fHistNN[i];
 155
 156   }
 157 }
 158
 159 //____________________________________________________________________//
 160 void AliBalance::SetInterval(Int_t iAnalysisType,
 161                              Double_t p1Start, Double_t p1Stop,
 162                              Int_t ibins, Double_t p2Start, Double_t p2Stop) {
 163   // Sets the analyzed interval.
 164   // Set the same Information for all analyses
 165
 166   if(iAnalysisType == -1){
 167     for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
 168       fP1Start[i] = p1Start;
 169       fP1Stop[i] = p1Stop;
 170       fNumberOfBins[i] = ibins;
 171       fP2Start[i] = p2Start;
 172       fP2Stop[i] = p2Stop;
 173       fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
 174     }
 175   }
 176   // Set the Information for one analysis
 177   else if((iAnalysisType > -1) && (iAnalysisType < ANALYSIS_TYPES)) {
 178     fP1Start[iAnalysisType] = p1Start;
 179     fP1Stop[iAnalysisType] = p1Stop;
 180     fNumberOfBins[iAnalysisType] = ibins;
 181     fP2Start[iAnalysisType] = p2Start;
 182     fP2Stop[iAnalysisType] = p2Stop;
 183     fP2Step[iAnalysisType] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[iAnalysisType];
 184   }
 185   else {
 186     AliError("Wrong ANALYSIS number!");
 187   }
 188 }
 189
 190 //____________________________________________________________________//
 191 void AliBalance::InitHistograms() {
 192   //Initialize the histograms
 193   TString histName;
 194   for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
 195     histName = "fHistP"; histName += gBFAnalysisType[iAnalysisType];
 196     if(bShuffle) histName.Append("_shuffle");
 197     if(fCentralityId) histName += fCentralityId.Data();
 198     fHistP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
 199
 200     histName = "fHistN"; histName += gBFAnalysisType[iAnalysisType];
 201     if(bShuffle) histName.Append("_shuffle");
 202     if(fCentralityId) histName += fCentralityId.Data();
 203     fHistN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
 204
 205     histName = "fHistPN"; histName += gBFAnalysisType[iAnalysisType];
 206     if(bShuffle) histName.Append("_shuffle");
 207     if(fCentralityId) histName += fCentralityId.Data();
 208     fHistPN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
 209
 210     histName = "fHistNP"; histName += gBFAnalysisType[iAnalysisType];
 211     if(bShuffle) histName.Append("_shuffle");
 212     if(fCentralityId) histName += fCentralityId.Data();
 213     fHistNP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
 214
 215     histName = "fHistPP"; histName += gBFAnalysisType[iAnalysisType];
 216     if(bShuffle) histName.Append("_shuffle");
 217     if(fCentralityId) histName += fCentralityId.Data();
 218     fHistPP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
 219
 220     histName = "fHistNN"; histName += gBFAnalysisType[iAnalysisType];
 221     if(bShuffle) histName.Append("_shuffle");
 222     if(fCentralityId) histName += fCentralityId.Data();
 223     fHistNN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
 224   }
 225 }
 226
 227 //____________________________________________________________________//
 228 void AliBalance::PrintAnalysisSettings() {
 229
 230   Printf("======================================");
 231   Printf("Analysis level: %s",fAnalysisLevel.Data());
 232   Printf("======================================");
 233   for(Int_t ibin = 0; ibin < ANALYSIS_TYPES; ibin++){
 234     Printf("Interval info for variable %d",ibin);
 235     Printf("Analyzed interval (min.): %lf",fP2Start[ibin]);
 236     Printf("Analyzed interval (max.): %lf",fP2Stop[ibin]);
 237     Printf("Number of bins: %d",fNumberOfBins[ibin]);
 238     Printf("Step: %lf",fP2Step[ibin]);
 239     Printf("          ");
 240   }
 241   Printf("======================================");
 242 }
 243
 244 //____________________________________________________________________//
 245 void AliBalance::CalculateBalance(Float_t fCentrality,vector<Double_t> **chargeVector,Float_t bSign) {
 246   // Calculates the balance function
 247   fAnalyzedEvents++;
 248   Int_t i = 0 , j = 0;
 249   Int_t iBin = 0;
 250
 251   // Initialize histograms if not done yet
 252   if(!fHistPN[0]){
 253     AliWarning("Histograms not yet initialized! --> Will be done now");
 254     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 255     InitHistograms();
 256   }
 257
 258   Int_t gNtrack = chargeVector[0]->size();
 259   //Printf("(AliBalance) Number of tracks: %d",gNtrack);
 260
 261   for(i = 0; i < gNtrack;i++){
 262       Short_t charge          = chargeVector[0]->at(i);
 263       Double_t rapidity       = chargeVector[1]->at(i);
 264       Double_t pseudorapidity = chargeVector[2]->at(i);
 265       Double_t phi            = chargeVector[3]->at(i);
 266
 267       //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
 268       for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
 269         if(iAnalysisType == kEta) {
 270           if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
 271             if(charge > 0) {
 272               fNp[iAnalysisType] += 1.;
 273               fHistP[iAnalysisType]->Fill(fCentrality,pseudorapidity);
 274             }//charge > 0
 275             if(charge < 0) {
 276               fNn[iAnalysisType] += 1.;
 277               fHistN[iAnalysisType]->Fill(fCentrality,pseudorapidity);
 278             }//charge < 0
 279           }//p1 interval check
 280         }//analysis type: eta
 281         else if(iAnalysisType == kPhi) {
 282           if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
 283             if(charge > 0) {
 284               fNp[iAnalysisType] += 1.;
 285               fHistP[iAnalysisType]->Fill(fCentrality,phi);
 286             }//charge > 0
 287             if(charge < 0) {
 288               fNn[iAnalysisType] += 1.;
 289               fHistN[iAnalysisType]->Fill(fCentrality,phi);
 290             }//charge < 0
 291           }//p1 interval check
 292         }//analysis type: phi
 293         else {
 294           if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
 295             if(charge > 0) {
 296               fNp[iAnalysisType] += 1.;
 297               fHistP[iAnalysisType]->Fill(fCentrality,rapidity);
 298             }//charge > 0
 299             if(charge < 0) {
 300               fNn[iAnalysisType] += 1.;
 301               fHistN[iAnalysisType]->Fill(fCentrality,rapidity);
 302             }//charge < 0
 303           }//p1 interval check
 304         }//analysis type: y, qside, qout, qlong, qinv
 305       }//analysis type loop
 306   }
 307
 308   //Printf("Np: %lf - Nn: %lf",fNp[0],fNn[0]);
 309
 310   Double_t dy = 0., deta = 0.;
 311   Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
 312   Double_t dphi = 0.;
 313
 314   Short_t charge1  = 0;
 315   Double_t eta1 = 0., rap1 = 0.;
 316   Double_t px1 = 0., py1 = 0., pz1 = 0.;
 317   Double_t pt1 = 0.;
 318   Double_t energy1 = 0.;
 319   Double_t phi1    = 0.;
 320
 321   Short_t charge2  = 0;
 322   Double_t eta2 = 0., rap2 = 0.;
 323   Double_t px2 = 0., py2 = 0., pz2 = 0.;
 324   Double_t pt2 = 0.;
 325   Double_t energy2 = 0.;
 326   Double_t phi2    = 0.;
 327   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
 328   for(i = 1; i < gNtrack; i++) {
 329
 330       charge1 = chargeVector[0]->at(i);
 331       rap1    = chargeVector[1]->at(i);
 332       eta1    = chargeVector[2]->at(i);
 333       phi1    = chargeVector[3]->at(i);
 334       px1     = chargeVector[4]->at(i);
 335       py1     = chargeVector[5]->at(i);
 336       pz1     = chargeVector[6]->at(i);
 337       pt1     = chargeVector[7]->at(i);
 338       energy1 = chargeVector[8]->at(i);
 339
 340     for(j = 0; j < i; j++) {
 341
 342       charge2 = chargeVector[0]->at(j);
 343       rap2    = chargeVector[1]->at(j);
 344       eta2    = chargeVector[2]->at(j);
 345       phi2    = chargeVector[3]->at(j);
 346       px2     = chargeVector[4]->at(j);
 347       py2     = chargeVector[5]->at(j);
 348       pz2     = chargeVector[6]->at(j);
 349       pt2     = chargeVector[7]->at(j);
 350       energy2 = chargeVector[8]->at(j);
 351
 352         // filling the arrays
 353
 354         // RAPIDITY
 355         dy = TMath::Abs(rap1 - rap2);
 356
 357         // Eta
 358         deta = TMath::Abs(eta1 - eta2);
 359
 360         //qlong
 361         Double_t eTot = energy1 + energy2;
 362         Double_t pxTot = px1 + px2;
 363         Double_t pyTot = py1 + py2;
 364         Double_t pzTot = pz1 + pz2;
 365         Double_t q0Tot = energy1 - energy2;
 366         Double_t qxTot = px1 - px2;
 367         Double_t qyTot = py1 - py2;
 368         Double_t qzTot = pz1 - pz2;
 369
 370         Double_t eTot2 = eTot*eTot;
 371         Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
 372         Double_t pzTot2 = pzTot*pzTot;
 373
 374         Double_t q0Tot2 = q0Tot*q0Tot;
 375         Double_t qTot2  = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;
 376
 377         Double_t snn    = eTot2 - pTot2;
 378         Double_t ptTot2 = pTot2 - pzTot2 ;
 379         Double_t ptTot  = TMath::Sqrt( ptTot2 );
 380
 381         qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
 382
 383         //qout
 384         qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
 385
 386         //qside
 387         qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
 388
 389         //qinv
 390         qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
 391
 392         //phi
 393         dphi = TMath::Abs(phi1 - phi2);
 394         if(dphi>180) dphi = 360 - dphi;  //dphi should be between 0 and 180!
 395
 396         // HBT like cut
 397         if(bHBTcut && charge1 * charge2 > 0){
 398           //if( dphi < 3 || deta < 0.01 ){   // VERSION 1
 399           //  continue;
 400
 401           // VERSION 2 (Taken from DPhiCorrelations)
 402           // the variables & cuthave been developed by the HBT group
 403           // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
 404
 405           // optimization
 406           if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
 407             {
 408
 409               // phi in rad
 410               Float_t phi1rad = phi1*TMath::DegToRad();
 411               Float_t phi2rad = phi2*TMath::DegToRad();
 412
 413               // check first boundaries to see if is worth to loop and find the minimum
 414               Float_t dphistar1 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 0.8, bSign);
 415               Float_t dphistar2 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 2.5, bSign);
 416
 417               const Float_t kLimit = 0.02 * 3;
 418
 419               Float_t dphistarminabs = 1e5;
 420               Float_t dphistarmin = 1e5;
 421
 422               if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 )
 423                 {
 424                   for (Double_t rad=0.8; rad<2.51; rad+=0.01)
 425                     {
 426                       Float_t dphistar = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, rad, bSign);
 427                       Float_t dphistarabs = TMath::Abs(dphistar);
 428
 429                       if (dphistarabs < dphistarminabs)
 430                         {
 431                           dphistarmin = dphistar;
 432                           dphistarminabs = dphistarabs;
 433                         }
 434                     }
 435
 436                   if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02)
 437                     {
 438                       //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));
 439                       continue;
 440                     }
 441                 }
 442             }
 443         }
 444
 445         // conversions
 446         if(bConversionCut){
 447           if (charge1 * charge2 < 0)
 448             {
 449               Float_t m0 = 0.510e-3;
 450               Float_t tantheta1 = 1e10;
 451
 452               // phi in rad
 453               Float_t phi1rad = phi1*TMath::DegToRad();
 454               Float_t phi2rad = phi2*TMath::DegToRad();
 455
 456               if (eta1 < -1e-10 || eta1 > 1e-10)
 457                 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
 458
 459               Float_t tantheta2 = 1e10;
 460               if (eta2 < -1e-10 || eta2 > 1e-10)
 461                 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
 462
 463               Float_t e1squ = m0 * m0 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
 464               Float_t e2squ = m0 * m0 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
 465
 466               Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
 467
 468               if (masssqu < 0.04*0.04){
 469                 //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));
 470                 continue;
 471               }
 472             }
 473         }
 474
 475
 476         //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
 477         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 478
 479           // rapidity
 480           if( dy > fP2Start[kRapidity] && dy < fP2Stop[kRapidity]){
 481             iBin = Int_t((dy-fP2Start[kRapidity])/fP2Step[kRapidity]);
 482             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 483
 484               if((charge1 > 0)&&(charge2 > 0)) {
 485                 fNpp[kRapidity][iBin] += 1.;
 486                 fHistPP[kRapidity]->Fill(fCentrality,dy);
 487               }
 488               else if((charge1 < 0)&&(charge2 < 0)) {
 489                 fNnn[kRapidity][iBin] += 1.;
 490                 fHistNN[kRapidity]->Fill(fCentrality,dy);
 491               }
 492               else if((charge1 > 0)&&(charge2 < 0)) {
 493                 fNpn[kRapidity][iBin] += 1.;
 494                 fHistPN[kRapidity]->Fill(fCentrality,dy);
 495               }
 496               else if((charge1 < 0)&&(charge2 > 0)) {
 497                 fNpn[kRapidity][iBin] += 1.;
 498                     fHistPN[kRapidity]->Fill(fCentrality,dy);
 499               }
 500             }//BF binning check
 501           }//p2 interval check
 502         }//p1 interval check
 503
 504         // pseudorapidity
 505         if((eta1 >= fP1Start[kEta]) && (eta1 <= fP1Stop[kEta]) && (eta2 >= fP1Start[kEta]) && (eta2 <= fP1Stop[kEta])) {
 506           if( deta > fP2Start[kEta] && deta < fP2Stop[kEta]){
 507             iBin = Int_t((deta-fP2Start[kEta])/fP2Step[kEta]);
 508             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 509               if((charge1 > 0)&&(charge2 > 0)) {
 510                 fNpp[kEta][iBin] += 1.;
 511                 fHistPP[kEta]->Fill(fCentrality,deta);
 512               }
 513               if((charge1 < 0)&&(charge2 < 0)) {
 514                 fNnn[kEta][iBin] += 1.;
 515                     fHistNN[kEta]->Fill(fCentrality,deta);
 516               }
 517               if((charge1 > 0)&&(charge2 < 0)) {
 518                 fNpn[kEta][iBin] += 1.;
 519                 fHistPN[kEta]->Fill(fCentrality,deta);
 520               }
 521               if((charge1 < 0)&&(charge2 > 0)) {
 522                 fNpn[kEta][iBin] += 1.;
 523                     fHistPN[kEta]->Fill(fCentrality,deta);
 524               }
 525             }//BF binning check
 526           }//p2 interval check
 527         }//p1 interval check
 528
 529         // Qlong, out, side, inv
 530         // Check the p1 intervall for rapidity here (like for single tracks above)
 531         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 532           if( qLong > fP2Start[kQlong] && qLong < fP2Stop[kQlong]){
 533             iBin = Int_t((qLong-fP2Start[kQlong])/fP2Step[kQlong]);
 534             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 535               if((charge1 > 0)&&(charge2 > 0)) {
 536                 fNpp[kQlong][iBin] += 1.;
 537                 fHistPP[kQlong]->Fill(fCentrality,qLong);
 538               }
 539               if((charge1 < 0)&&(charge2 < 0)) {
 540                 fNnn[kQlong][iBin] += 1.;
 541                 fHistNN[kQlong]->Fill(fCentrality,qLong);
 542               }
 543               if((charge1 > 0)&&(charge2 < 0)) {
 544                 fNpn[kQlong][iBin] += 1.;
 545                 fHistPN[kQlong]->Fill(fCentrality,qLong);
 546               }
 547               if((charge1 < 0)&&(charge2 > 0)) {
 548                 fNpn[kQlong][iBin] += 1.;
 549                 fHistPN[kQlong]->Fill(fCentrality,qLong);
 550               }
 551             }//BF binning check
 552           }//p2 interval check
 553         }//p1 interval check
 554
 555         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 556           if( qOut > fP2Start[kQout] && qOut < fP2Stop[kQout]){
 557             iBin = Int_t((qOut-fP2Start[kQout])/fP2Step[kQout]);
 558             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 559               if((charge1 > 0)&&(charge2 > 0)) {
 560                 fNpp[kQout][iBin] += 1.;
 561                 fHistPP[kQout]->Fill(fCentrality,qOut);
 562                   }
 563               if((charge1 < 0)&&(charge2 < 0)) {
 564                 fNnn[kQout][iBin] += 1.;
 565                 fHistNN[kQout]->Fill(fCentrality,qOut);
 566               }
 567               if((charge1 > 0)&&(charge2 < 0)) {
 568                 fNpn[kQout][iBin] += 1.;
 569                 fHistPN[kQout]->Fill(fCentrality,qOut);
 570               }
 571               if((charge1 < 0)&&(charge2 > 0)) {
 572                 fNpn[kQout][iBin] += 1.;
 573                 fHistPN[kQout]->Fill(fCentrality,qOut);
 574               }
 575             }//BF binning check
 576           }//p2 interval check
 577         }//p1 interval check
 578
 579         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 580           if( qSide > fP2Start[kQside] && qSide < fP2Stop[kQside]){
 581             iBin = Int_t((qSide-fP2Start[kQside])/fP2Step[kQside]);
 582             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 583               if((charge1 > 0)&&(charge2 > 0)) {
 584                 fNpp[kQside][iBin] += 1.;
 585                 fHistPP[kQside]->Fill(fCentrality,qSide);
 586               }
 587               if((charge1 < 0)&&(charge2 < 0)) {
 588                 fNnn[kQside][iBin] += 1.;
 589                 fHistNN[kQside]->Fill(fCentrality,qSide);
 590               }
 591               if((charge1 > 0)&&(charge2 < 0)) {
 592                 fNpn[kQside][iBin] += 1.;
 593                 fHistPN[kQside]->Fill(fCentrality,qSide);
 594               }
 595               if((charge1 < 0)&&(charge2 > 0)) {
 596                 fNpn[kQside][iBin] += 1.;
 597                 fHistPN[kQside]->Fill(fCentrality,qSide);
 598                           }
 599             }//BF binning check
 600           }//p2 interval check
 601         }//p1 interval check
 602
 603         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 604           if( qInv > fP2Start[kQinv] && qInv < fP2Stop[kQinv]){
 605             iBin = Int_t((qInv-fP2Start[kQinv])/fP2Step[kQinv]);
 606             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 607               if((charge1 > 0)&&(charge2 > 0)) {
 608                 fNpp[kQinv][iBin] += 1.;
 609                 fHistPP[kQinv]->Fill(fCentrality,qInv);
 610               }
 611               if((charge1 < 0)&&(charge2 < 0)) {
 612                 fNnn[kQinv][iBin] += 1.;
 613                 fHistNN[kQinv]->Fill(fCentrality,qInv);
 614               }
 615               if((charge1 > 0)&&(charge2 < 0)) {
 616                 fNpn[kQinv][iBin] += 1.;
 617                 fHistPN[kQinv]->Fill(fCentrality,qInv);
 618               }
 619               if((charge1 < 0)&&(charge2 > 0)) {
 620                 fNpn[kQinv][iBin] += 1.;
 621                 fHistPN[kQinv]->Fill(fCentrality,qInv);
 622               }
 623             }//BF binning check
 624           }//p2 interval check
 625         }//p1 interval check
 626
 627         // Phi
 628         if((phi1 >= fP1Start[kPhi]) && (phi1 <= fP1Stop[kPhi]) && (phi2 >= fP1Start[kPhi]) && (phi2 <= fP1Stop[kPhi])) {
 629           if( dphi > fP2Start[kPhi] && dphi < fP2Stop[kPhi]){
 630             iBin = Int_t((dphi-fP2Start[kPhi])/fP2Step[kPhi]);
 631             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 632               if((charge1 > 0)&&(charge2 > 0)) {
 633                 fNpp[kPhi][iBin] += 1.;
 634                 fHistPP[kPhi]->Fill(fCentrality,dphi);
 635               }
 636               if((charge1 < 0)&&(charge2 < 0)) {
 637                 fNnn[kPhi][iBin] += 1.;
 638                 fHistNN[kPhi]->Fill(fCentrality,dphi);
 639               }
 640               if((charge1 > 0)&&(charge2 < 0)) {
 641                 fNpn[kPhi][iBin] += 1.;
 642                 fHistPN[kPhi]->Fill(fCentrality,dphi);
 643               }
 644               if((charge1 < 0)&&(charge2 > 0)) {
 645                 fNpn[kPhi][iBin] += 1.;
 646                 fHistPN[kPhi]->Fill(fCentrality,dphi);
 647               }
 648             }//BF binning check
 649           }//p2 interval check
 650         }//p1 interval check
 651     }//end of 2nd particle loop
 652   }//end of 1st particle loop
 653   //Printf("Number of analyzed events: %i",fAnalyzedEvents);
 654   //Printf("DeltaEta NN[0] = %.0f, PP[0] = %.0f, NP[0] = %.0f, PN[0] = %.0f",fNnn[kEta][0],fNpp[kEta][0],fNnp[kEta][0],fNpn[kEta][0]);
 655 }
 656
 657
 658 //____________________________________________________________________//
 659 Double_t AliBalance::GetBalance(Int_t iAnalysisType, Int_t p2) {
 660   // Returns the value of the balance function in bin p2
 661   fB[iAnalysisType][p2] = 0.5*(((fNpn[iAnalysisType][p2] - 2.*fNnn[iAnalysisType][p2])/fNn[iAnalysisType]) + ((fNpn[iAnalysisType][p2] - 2.*fNpp[iAnalysisType][p2])/fNp[iAnalysisType]))/fP2Step[iAnalysisType];
 662
 663   return fB[iAnalysisType][p2];
 664 }
 665
 666 //____________________________________________________________________//
 667 Double_t AliBalance::GetError(Int_t iAnalysisType, Int_t p2) {
 668   // Returns the error on the BF value for bin p2
 669   // The errors for fNn and fNp are neglected here (0.1 % of total error)
 670   /*ferror[iAnalysisType][p2] = TMath::Sqrt(Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
 671                               + Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType]))
 672                               + Double_t(fNpn[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
 673                               + Double_t(fNnp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
 674                               //+ TMath::Power(fNpn[iAnalysisType][p2]-fNpp[iAnalysisType][p2],2)/TMath::Power(Double_t(fNp[iAnalysisType]),3)
 675                               //+ TMath::Power(fNnp[iAnalysisType][p2]-fNnn[iAnalysisType][p2],2)/TMath::Power(Double_t(fNn[iAnalysisType]),3)
 676                                ) /fP2Step[iAnalysisType];*/
 677
 678   ferror[iAnalysisType][p2] = TMath::Sqrt( Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType])) +
 679                                            Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType])) +
 680                                            Double_t(fNpn[iAnalysisType][p2])*TMath::Power((0.5/Double_t(fNp[iAnalysisType]) + 0.5/Double_t(fNn[iAnalysisType])),2))/fP2Step[iAnalysisType];
 681
 682   return ferror[iAnalysisType][p2];
 683 }
 684 //____________________________________________________________________//
 685 TGraphErrors *AliBalance::DrawBalance(Int_t iAnalysisType) {
 686
 687   // Draws the BF
 688   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
 689   Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
 690   Double_t b[MAXIMUM_NUMBER_OF_STEPS];
 691   Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
 692
 693   if((fNp[iAnalysisType] == 0)||(fNn[iAnalysisType] == 0)) {
 694     cerr<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
 695     return NULL;
 696   }
 697
 698   for(Int_t i = 0; i < fNumberOfBins[iAnalysisType]; i++) {
 699     b[i] = GetBalance(iAnalysisType,i);
 700     ber[i] = GetError(iAnalysisType,i);
 701     x[i] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
 702     xer[i] = 0.0;
 703   }
 704
 705   TGraphErrors *gr = new TGraphErrors(fNumberOfBins[iAnalysisType],x,b,xer,ber);
 706   gr->GetXaxis()->SetTitleColor(1);
 707   if(iAnalysisType==0) {
 708     gr->SetTitle("Balance function B(#Delta y)");
 709     gr->GetXaxis()->SetTitle("#Delta y");
 710     gr->GetYaxis()->SetTitle("B(#Delta y)");
 711   }
 712   if(iAnalysisType==1) {
 713     gr->SetTitle("Balance function B(#Delta #eta)");
 714     gr->GetXaxis()->SetTitle("#Delta #eta");
 715     gr->GetYaxis()->SetTitle("B(#Delta #eta)");
 716   }
 717   if(iAnalysisType==2) {
 718     gr->SetTitle("Balance function B(q_{long})");
 719     gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
 720     gr->GetYaxis()->SetTitle("B(q_{long}) ((GeV/c)^{-1})");
 721   }
 722   if(iAnalysisType==3) {
 723     gr->SetTitle("Balance function B(q_{out})");
 724     gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
 725     gr->GetYaxis()->SetTitle("B(q_{out}) ((GeV/c)^{-1})");
 726   }
 727   if(iAnalysisType==4) {
 728     gr->SetTitle("Balance function B(q_{side})");
 729     gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
 730     gr->GetYaxis()->SetTitle("B(q_{side}) ((GeV/c)^{-1})");
 731   }
 732   if(iAnalysisType==5) {
 733     gr->SetTitle("Balance function B(q_{inv})");
 734     gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 735     gr->GetYaxis()->SetTitle("B(q_{inv}) ((GeV/c)^{-1})");
 736   }
 737   if(iAnalysisType==6) {
 738     gr->SetTitle("Balance function B(#Delta #phi)");
 739     gr->GetXaxis()->SetTitle("#Delta #phi");
 740     gr->GetYaxis()->SetTitle("B(#Delta #phi)");
 741   }
 742
 743   return gr;
 744 }
 745
 746 //____________________________________________________________________//
 747 void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
 748   //Prints the calculated width of the BF and its error
 749   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
 750   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
 751   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
 752   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
 753   Double_t deltaBalP2 = 0.0, integral = 0.0;
 754   Double_t deltaErrorNew = 0.0;
 755
 756   cout<<"=================================================="<<endl;
 757   for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) {
 758     x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
 759     //cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
 760   }
 761   //cout<<"=================================================="<<endl;
 762   for(Int_t i = 2; i <= fNumberOfBins[iAnalysisType]; i++) {
 763     gSumXi += gHistBalance->GetBinCenter(i);
 764     gSumBi += gHistBalance->GetBinContent(i);
 765     gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
 766     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
 767     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
 768     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
 769     gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
 770
 771     deltaBalP2 += fP2Step[iAnalysisType]*TMath::Power(gHistBalance->GetBinError(i),2);
 772     integral += fP2Step[iAnalysisType]*gHistBalance->GetBinContent(i);
 773   }
 774   for(Int_t i = 1; i < fNumberOfBins[iAnalysisType]; i++)
 775     deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
 776
 777   Double_t integralError = TMath::Sqrt(deltaBalP2);
 778
 779   Double_t delta = gSumBiXi / gSumBi;
 780   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
 781   cout<<"Analysis type: "<<gBFAnalysisType[iAnalysisType].Data()<<endl;
 782   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 783   cout<<"New error: "<<deltaErrorNew<<endl;
 784   cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
 785   cout<<"=================================================="<<endl;
 786 }
 787
 788 //____________________________________________________________________//
 789 TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow) {
 790   //Returns the BF histogram, extracted from the 6 TH2D objects
 791   //(private members) of the AliBalance class.
 792   //
 793   // Acceptance correction:
 794   // - only for analysis type = kEta
 795   // - only if etaWindow > 0 (default = -1.)
 796   // - calculated as proposed by STAR
 797   //
 798   TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
 799   TString histName = "gHistBalanceFunctionHistogram";
 800   histName += gAnalysisType[iAnalysisType];
 801
 802   SetInterval(iAnalysisType, fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
 803               fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
 804               fHistPP[iAnalysisType]->GetNbinsY(),
 805               fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),
 806               fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
 807
 808   // determine the projection thresholds
 809   Int_t binMinX, binMinY, binMinZ;
 810   Int_t binMaxX, binMaxY, binMaxZ;
 811
 812   fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMin),binMinX,binMinY,binMinZ);
 813   fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMax),binMaxX,binMaxY,binMaxZ);
 814
 815   TH1D *gHistBalanceFunctionHistogram = new TH1D(histName.Data(),"",fHistPP[iAnalysisType]->GetNbinsY(),fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
 816   switch(iAnalysisType) {
 817   case kRapidity:
 818     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
 819     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
 820     break;
 821   case kEta:
 822     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
 823     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
 824     break;
 825   case kQlong:
 826     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
 827     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
 828     break;
 829   case kQout:
 830     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
 831     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
 832     break;
 833   case kQside:
 834     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
 835     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
 836     break;
 837   case kQinv:
 838     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 839     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
 840     break;
 841   case kPhi:
 842     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
 843     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
 844     break;
 845   default:
 846     break;
 847   }
 848
 849   TH1D *hTemp1 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 850   TH1D *hTemp2 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f_copy",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 851   TH1D *hTemp3 = dynamic_cast<TH1D *>(fHistNN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistNN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 852   TH1D *hTemp4 = dynamic_cast<TH1D *>(fHistPP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 853   TH1D *hTemp5 = dynamic_cast<TH1D *>(fHistN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 854   TH1D *hTemp6 = dynamic_cast<TH1D *>(fHistP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 855
 856   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
 857     hTemp1->Sumw2();
 858     hTemp2->Sumw2();
 859     hTemp3->Sumw2();
 860     hTemp4->Sumw2();
 861     hTemp1->Add(hTemp3,-2.);
 862     hTemp1->Scale(1./hTemp5->GetEntries());
 863     hTemp2->Add(hTemp4,-2.);
 864     hTemp2->Scale(1./hTemp6->GetEntries());
 865     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 866     gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]);
 867   }
 868
 869   // do the acceptance correction (only for Eta and etaWindow > 0)
 870   if(iAnalysisType == kEta && etaWindow > 0){
 871     for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
 872
 873       Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
 874       Double_t corrected    = notCorrected / (1 - (gHistBalanceFunctionHistogram->GetBinCenter(iBin+1))/ etaWindow );
 875       gHistBalanceFunctionHistogram->SetBinContent(iBin+1, corrected);
 876
 877     }
 878   }
 879
 880   PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);
 881
 882   return gHistBalanceFunctionHistogram;
 883 }