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