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