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