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 <TFile.h>
  28 #include <TF1.h>
  29 #include <TH2D.h>
  30 #include <TLorentzVector.h>
  31 #include <TObjArray.h>
  32 #include <TGraphErrors.h>
  33 #include <TString.h>
  34
  35 #include "AliVParticle.h"
  36 #include "AliMCParticle.h"
  37 #include "AliESDtrack.h"
  38 #include "AliAODTrack.h"
  39
  40 #include "AliBalance.h"
  41
  42 using std::cout;
  43 using std::cerr;
  44 using std::endl;
  45
  46 ClassImp(AliBalance)
  47
  48 //____________________________________________________________________//
  49 AliBalance::AliBalance() :
  50   TObject(),
  51   fShuffle(kFALSE),
  52   fHBTcut(kFALSE),
  53   fConversionCut(kFALSE),
  54   fAnalysisLevel("ESD"),
  55   fAnalyzedEvents(0) ,
  56   fCentralityId(0) ,
  57   fCentStart(0.),
  58   fCentStop(0.),
  59   fHistHBTbefore(NULL),
  60   fHistHBTafter(NULL),
  61   fHistConversionbefore(NULL),
  62   fHistConversionafter(NULL)
  63 {
  64   // Default constructor
  65
  66   for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
  67     if(i == 6) {
  68       fNumberOfBins[i] = 180;
  69       fP1Start[i]      = -360.0;
  70       fP1Stop[i]       = 360.0;
  71       fP2Start[i]      = -360.0;
  72       fP2Stop[i]       = 360.0;
  73       fP2Step[i]       = 0.1;
  74     }
  75     else {
  76       fNumberOfBins[i] = 20;
  77       fP1Start[i]      = -1.0;
  78       fP1Stop[i]       = 1.0;
  79       fP2Start[i]      = 0.0;
  80       fP2Stop[i]       = 2.0;
  81     }
  82     fP2Step[i] = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins[i];
  83     fCentStart = 0.;
  84     fCentStop  = 0.;
  85
  86     fNn[i] = 0.0;
  87     fNp[i] = 0.0;
  88
  89     for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
  90       fNpp[i][j] = .0;
  91       fNnn[i][j] = .0;
  92       fNpn[i][j] = .0;
  93       fNnp[i][j] = .0;
  94       fB[i][j] = 0.0;
  95       ferror[i][j] = 0.0;
  96     }
  97
  98     fHistP[i]  = NULL;
  99     fHistN[i]  = NULL;
 100     fHistPP[i] = NULL;
 101     fHistPN[i] = NULL;
 102     fHistNP[i] = NULL;
 103     fHistNN[i] = NULL;
 104
 105   }
 106 }
 107
 108
 109 //____________________________________________________________________//
 110 AliBalance::AliBalance(const AliBalance& balance):
 111   TObject(balance),
 112   fShuffle(balance.fShuffle),
 113   fHBTcut(balance.fHBTcut),
 114   fConversionCut(balance.fConversionCut),
 115   fAnalysisLevel(balance.fAnalysisLevel),
 116   fAnalyzedEvents(balance.fAnalyzedEvents),
 117   fCentralityId(balance.fCentralityId),
 118   fCentStart(balance.fCentStart),
 119   fCentStop(balance.fCentStop),
 120   fHistHBTbefore(balance.fHistHBTbefore),
 121   fHistHBTafter(balance.fHistHBTafter),
 122   fHistConversionbefore(balance.fHistConversionbefore),
 123   fHistConversionafter(balance.fHistConversionafter) {
 124   //copy constructor
 125   for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
 126     fNn[i] = balance.fNn[i];
 127     fNp[i] = balance.fNp[i];
 128
 129     fP1Start[i]      = balance.fP1Start[i];
 130     fP1Stop[i]       = balance.fP1Stop[i];
 131     fNumberOfBins[i] = balance.fNumberOfBins[i];
 132     fP2Start[i]      = balance.fP2Start[i];
 133     fP2Stop[i]       = balance.fP2Stop[i];
 134     fP2Step[i]       = balance.fP2Step[i];
 135     fCentStart       = balance.fCentStart;
 136     fCentStop        = balance.fCentStop;
 137
 138     fHistP[i]        = balance.fHistP[i];
 139     fHistN[i]        = balance.fHistN[i];
 140     fHistPN[i]        = balance.fHistPN[i];
 141     fHistNP[i]        = balance.fHistNP[i];
 142     fHistPP[i]        = balance.fHistPP[i];
 143     fHistNN[i]        = balance.fHistNN[i];
 144
 145     for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
 146       fNpp[i][j] = .0;
 147       fNnn[i][j] = .0;
 148       fNpn[i][j] = .0;
 149       fNnp[i][j] = .0;
 150       fB[i][j] = 0.0;
 151       ferror[i][j] = 0.0;
 152     }
 153   }
 154  }
 155
 156
 157 //____________________________________________________________________//
 158 AliBalance::~AliBalance() {
 159   // Destructor
 160
 161   for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
 162
 163     delete fHistP[i];
 164     delete fHistN[i];
 165     delete fHistPN[i];
 166     delete fHistNP[i];
 167     delete fHistPP[i];
 168     delete fHistNN[i];
 169
 170   }
 171 }
 172
 173 //____________________________________________________________________//
 174 void AliBalance::SetInterval(Int_t iAnalysisType,
 175                              Double_t p1Start, Double_t p1Stop,
 176                              Int_t ibins, Double_t p2Start, Double_t p2Stop) {
 177   // Sets the analyzed interval.
 178   // Set the same Information for all analyses
 179
 180   if(iAnalysisType == -1){
 181     for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
 182       fP1Start[i] = p1Start;
 183       fP1Stop[i] = p1Stop;
 184       fNumberOfBins[i] = ibins;
 185       fP2Start[i] = p2Start;
 186       fP2Stop[i] = p2Stop;
 187       fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
 188     }
 189   }
 190   // Set the Information for one analysis
 191   else if((iAnalysisType > -1) && (iAnalysisType < ANALYSIS_TYPES)) {
 192     fP1Start[iAnalysisType] = p1Start;
 193     fP1Stop[iAnalysisType] = p1Stop;
 194     fNumberOfBins[iAnalysisType] = ibins;
 195     fP2Start[iAnalysisType] = p2Start;
 196     fP2Stop[iAnalysisType] = p2Stop;
 197     fP2Step[iAnalysisType] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[iAnalysisType];
 198   }
 199   else {
 200     AliError("Wrong ANALYSIS number!");
 201   }
 202 }
 203
 204 //____________________________________________________________________//
 205 void AliBalance::InitHistograms() {
 206   //Initialize the histograms
 207
 208   // global switch disabling the reference
 209   // (to avoid "Replacing existing TH1" if several wagons are created in train)
 210   Bool_t oldStatus = TH1::AddDirectoryStatus();
 211   TH1::AddDirectory(kFALSE);
 212
 213   TString histName;
 214   for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
 215     histName = "fHistP"; histName += kBFAnalysisType[iAnalysisType];
 216     if(fShuffle) histName.Append("_shuffle");
 217     if(fCentralityId) histName += fCentralityId.Data();
 218     fHistP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
 219
 220     histName = "fHistN"; histName += kBFAnalysisType[iAnalysisType];
 221     if(fShuffle) histName.Append("_shuffle");
 222     if(fCentralityId) histName += fCentralityId.Data();
 223     fHistN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
 224
 225     histName = "fHistPN"; histName += kBFAnalysisType[iAnalysisType];
 226     if(fShuffle) histName.Append("_shuffle");
 227     if(fCentralityId) histName += fCentralityId.Data();
 228     fHistPN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
 229
 230     histName = "fHistNP"; histName += kBFAnalysisType[iAnalysisType];
 231     if(fShuffle) histName.Append("_shuffle");
 232     if(fCentralityId) histName += fCentralityId.Data();
 233     fHistNP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
 234
 235     histName = "fHistPP"; histName += kBFAnalysisType[iAnalysisType];
 236     if(fShuffle) histName.Append("_shuffle");
 237     if(fCentralityId) histName += fCentralityId.Data();
 238     fHistPP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
 239
 240     histName = "fHistNN"; histName += kBFAnalysisType[iAnalysisType];
 241     if(fShuffle) histName.Append("_shuffle");
 242     if(fCentralityId) histName += fCentralityId.Data();
 243     fHistNN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
 244   }
 245
 246   // QA histograms
 247   fHistHBTbefore        = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,200);
 248   fHistHBTafter         = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,200);
 249   fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,200);
 250   fHistConversionafter  = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,200);
 251
 252   TH1::AddDirectory(oldStatus);
 253
 254 }
 255
 256 //____________________________________________________________________//
 257 void AliBalance::PrintAnalysisSettings() {
 258   //prints the analysis settings
 259
 260   Printf("======================================");
 261   Printf("Analysis level: %s",fAnalysisLevel.Data());
 262   Printf("======================================");
 263   for(Int_t ibin = 0; ibin < ANALYSIS_TYPES; ibin++){
 264     Printf("Interval info for variable %d",ibin);
 265     Printf("Analyzed interval (min.): %lf",fP2Start[ibin]);
 266     Printf("Analyzed interval (max.): %lf",fP2Stop[ibin]);
 267     Printf("Number of bins: %d",fNumberOfBins[ibin]);
 268     Printf("Step: %lf",fP2Step[ibin]);
 269     Printf("          ");
 270   }
 271   Printf("======================================");
 272 }
 273
 274 //____________________________________________________________________//
 275 void AliBalance::CalculateBalance(Float_t fCentrality,vector<Double_t> **chargeVector,Float_t bSign) {
 276   // Calculates the balance function
 277   fAnalyzedEvents++;
 278   Int_t i = 0 , j = 0;
 279   Int_t iBin = 0;
 280
 281   // Initialize histograms if not done yet
 282   if(!fHistPN[0]){
 283     AliWarning("Histograms not yet initialized! --> Will be done now");
 284     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 285     InitHistograms();
 286   }
 287
 288   Int_t gNtrack = chargeVector[0]->size();
 289   //Printf("(AliBalance) Number of tracks: %d",gNtrack);
 290
 291   for(i = 0; i < gNtrack;i++){
 292       Short_t charge          = chargeVector[0]->at(i);
 293       Double_t rapidity       = chargeVector[1]->at(i);
 294       Double_t pseudorapidity = chargeVector[2]->at(i);
 295       Double_t phi            = chargeVector[3]->at(i);
 296
 297       //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
 298       for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
 299         if(iAnalysisType == kEta) {
 300           if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
 301             if(charge > 0) {
 302               fNp[iAnalysisType] += 1.;
 303               fHistP[iAnalysisType]->Fill(fCentrality,pseudorapidity);
 304             }//charge > 0
 305             if(charge < 0) {
 306               fNn[iAnalysisType] += 1.;
 307               fHistN[iAnalysisType]->Fill(fCentrality,pseudorapidity);
 308             }//charge < 0
 309           }//p1 interval check
 310         }//analysis type: eta
 311         else if(iAnalysisType == kPhi) {
 312           if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
 313             if(charge > 0) {
 314               fNp[iAnalysisType] += 1.;
 315               fHistP[iAnalysisType]->Fill(fCentrality,phi);
 316             }//charge > 0
 317             if(charge < 0) {
 318               fNn[iAnalysisType] += 1.;
 319               fHistN[iAnalysisType]->Fill(fCentrality,phi);
 320             }//charge < 0
 321           }//p1 interval check
 322         }//analysis type: phi
 323         else {
 324           if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
 325             if(charge > 0) {
 326               fNp[iAnalysisType] += 1.;
 327               fHistP[iAnalysisType]->Fill(fCentrality,rapidity);
 328             }//charge > 0
 329             if(charge < 0) {
 330               fNn[iAnalysisType] += 1.;
 331               fHistN[iAnalysisType]->Fill(fCentrality,rapidity);
 332             }//charge < 0
 333           }//p1 interval check
 334         }//analysis type: y, qside, qout, qlong, qinv
 335       }//analysis type loop
 336   }
 337
 338   //Printf("Np: %lf - Nn: %lf",fNp[0],fNn[0]);
 339
 340   Double_t dy = 0., deta = 0.;
 341   Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
 342   Double_t dphi = 0.;
 343
 344   Short_t charge1  = 0;
 345   Double_t eta1 = 0., rap1 = 0.;
 346   Double_t px1 = 0., py1 = 0., pz1 = 0.;
 347   Double_t pt1 = 0.;
 348   Double_t energy1 = 0.;
 349   Double_t phi1    = 0.;
 350
 351   Short_t charge2  = 0;
 352   Double_t eta2 = 0., rap2 = 0.;
 353   Double_t px2 = 0., py2 = 0., pz2 = 0.;
 354   Double_t pt2 = 0.;
 355   Double_t energy2 = 0.;
 356   Double_t phi2    = 0.;
 357   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
 358   for(i = 1; i < gNtrack; i++) {
 359
 360       charge1 = chargeVector[0]->at(i);
 361       rap1    = chargeVector[1]->at(i);
 362       eta1    = chargeVector[2]->at(i);
 363       phi1    = chargeVector[3]->at(i);
 364       px1     = chargeVector[4]->at(i);
 365       py1     = chargeVector[5]->at(i);
 366       pz1     = chargeVector[6]->at(i);
 367       pt1     = chargeVector[7]->at(i);
 368       energy1 = chargeVector[8]->at(i);
 369
 370     for(j = 0; j < i; j++) {
 371
 372       charge2 = chargeVector[0]->at(j);
 373       rap2    = chargeVector[1]->at(j);
 374       eta2    = chargeVector[2]->at(j);
 375       phi2    = chargeVector[3]->at(j);
 376       px2     = chargeVector[4]->at(j);
 377       py2     = chargeVector[5]->at(j);
 378       pz2     = chargeVector[6]->at(j);
 379       pt2     = chargeVector[7]->at(j);
 380       energy2 = chargeVector[8]->at(j);
 381
 382         // filling the arrays
 383
 384         // RAPIDITY
 385         dy = TMath::Abs(rap1 - rap2);
 386
 387         // Eta
 388         deta = TMath::Abs(eta1 - eta2);
 389
 390         //qlong
 391         Double_t eTot = energy1 + energy2;
 392         Double_t pxTot = px1 + px2;
 393         Double_t pyTot = py1 + py2;
 394         Double_t pzTot = pz1 + pz2;
 395         Double_t q0Tot = energy1 - energy2;
 396         Double_t qxTot = px1 - px2;
 397         Double_t qyTot = py1 - py2;
 398         Double_t qzTot = pz1 - pz2;
 399
 400         Double_t eTot2 = eTot*eTot;
 401         Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
 402         Double_t pzTot2 = pzTot*pzTot;
 403
 404         Double_t q0Tot2 = q0Tot*q0Tot;
 405         Double_t qTot2  = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;
 406
 407         Double_t snn    = eTot2 - pTot2;
 408         Double_t ptTot2 = pTot2 - pzTot2 ;
 409         Double_t ptTot  = TMath::Sqrt( ptTot2 );
 410
 411         qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
 412
 413         //qout
 414         qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
 415
 416         //qside
 417         qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
 418
 419         //qinv
 420         qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
 421
 422         //phi
 423         dphi = TMath::Abs(phi1 - phi2);
 424         if(dphi>180) dphi = 360 - dphi;  //dphi should be between 0 and 180!
 425
 426         // HBT like cut
 427         if(fHBTcut && charge1 * charge2 > 0){
 428           //if( dphi < 3 || deta < 0.01 ){   // VERSION 1
 429           //  continue;
 430
 431           // VERSION 2 (Taken from DPhiCorrelations)
 432           // the variables & cuthave been developed by the HBT group
 433           // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
 434
 435           fHistHBTbefore->Fill(deta,dphi);
 436
 437           // optimization
 438           if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
 439             {
 440
 441               // phi in rad
 442               Float_t phi1rad = phi1*TMath::DegToRad();
 443               Float_t phi2rad = phi2*TMath::DegToRad();
 444
 445               // check first boundaries to see if is worth to loop and find the minimum
 446               Float_t dphistar1 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 0.8, bSign);
 447               Float_t dphistar2 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 2.5, bSign);
 448
 449               const Float_t kLimit = 0.02 * 3;
 450
 451               Float_t dphistarminabs = 1e5;
 452
 453               if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 )
 454                 {
 455                   for (Double_t rad=0.8; rad<2.51; rad+=0.01)
 456                     {
 457                       Float_t dphistar = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, rad, bSign);
 458                       Float_t dphistarabs = TMath::Abs(dphistar);
 459
 460                       if (dphistarabs < dphistarminabs)
 461                         {
 462                           dphistarminabs = dphistarabs;
 463                         }
 464                     }
 465
 466                   if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02)
 467                     {
 468                       //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));
 469                       continue;
 470                     }
 471                 }
 472             }
 473           fHistHBTafter->Fill(deta,dphi);
 474         }
 475
 476         // conversions
 477         if(fConversionCut){
 478           if (charge1 * charge2 < 0)
 479             {
 480
 481               fHistConversionbefore->Fill(deta,dphi);
 482
 483               Float_t m0 = 0.510e-3;
 484               Float_t tantheta1 = 1e10;
 485
 486               // phi in rad
 487               Float_t phi1rad = phi1*TMath::DegToRad();
 488               Float_t phi2rad = phi2*TMath::DegToRad();
 489
 490               if (eta1 < -1e-10 || eta1 > 1e-10)
 491                 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
 492
 493               Float_t tantheta2 = 1e10;
 494               if (eta2 < -1e-10 || eta2 > 1e-10)
 495                 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
 496
 497               Float_t e1squ = m0 * m0 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
 498               Float_t e2squ = m0 * m0 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
 499
 500               Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
 501
 502               if (masssqu < 0.04*0.04){
 503                 //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));
 504                 continue;
 505               }
 506               fHistConversionafter->Fill(deta,dphi);
 507             }
 508         }
 509
 510
 511         //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
 512         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 513
 514           // rapidity
 515           if( dy > fP2Start[kRapidity] && dy < fP2Stop[kRapidity]){
 516             iBin = Int_t((dy-fP2Start[kRapidity])/fP2Step[kRapidity]);
 517             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 518
 519               if((charge1 > 0)&&(charge2 > 0)) {
 520                 fNpp[kRapidity][iBin] += 1.;
 521                 fHistPP[kRapidity]->Fill(fCentrality,dy);
 522               }
 523               else if((charge1 < 0)&&(charge2 < 0)) {
 524                 fNnn[kRapidity][iBin] += 1.;
 525                 fHistNN[kRapidity]->Fill(fCentrality,dy);
 526               }
 527               else if((charge1 > 0)&&(charge2 < 0)) {
 528                 fNpn[kRapidity][iBin] += 1.;
 529                 fHistPN[kRapidity]->Fill(fCentrality,dy);
 530               }
 531               else if((charge1 < 0)&&(charge2 > 0)) {
 532                 fNpn[kRapidity][iBin] += 1.;
 533                     fHistPN[kRapidity]->Fill(fCentrality,dy);
 534               }
 535             }//BF binning check
 536           }//p2 interval check
 537         }//p1 interval check
 538
 539         // pseudorapidity
 540         if((eta1 >= fP1Start[kEta]) && (eta1 <= fP1Stop[kEta]) && (eta2 >= fP1Start[kEta]) && (eta2 <= fP1Stop[kEta])) {
 541           if( deta > fP2Start[kEta] && deta < fP2Stop[kEta]){
 542             iBin = Int_t((deta-fP2Start[kEta])/fP2Step[kEta]);
 543             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 544               if((charge1 > 0)&&(charge2 > 0)) {
 545                 fNpp[kEta][iBin] += 1.;
 546                 fHistPP[kEta]->Fill(fCentrality,deta);
 547               }
 548               if((charge1 < 0)&&(charge2 < 0)) {
 549                 fNnn[kEta][iBin] += 1.;
 550                     fHistNN[kEta]->Fill(fCentrality,deta);
 551               }
 552               if((charge1 > 0)&&(charge2 < 0)) {
 553                 fNpn[kEta][iBin] += 1.;
 554                 fHistPN[kEta]->Fill(fCentrality,deta);
 555               }
 556               if((charge1 < 0)&&(charge2 > 0)) {
 557                 fNpn[kEta][iBin] += 1.;
 558                     fHistPN[kEta]->Fill(fCentrality,deta);
 559               }
 560             }//BF binning check
 561           }//p2 interval check
 562         }//p1 interval check
 563
 564         // Qlong, out, side, inv
 565         // Check the p1 intervall for rapidity here (like for single tracks above)
 566         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 567           if( qLong > fP2Start[kQlong] && qLong < fP2Stop[kQlong]){
 568             iBin = Int_t((qLong-fP2Start[kQlong])/fP2Step[kQlong]);
 569             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 570               if((charge1 > 0)&&(charge2 > 0)) {
 571                 fNpp[kQlong][iBin] += 1.;
 572                 fHistPP[kQlong]->Fill(fCentrality,qLong);
 573               }
 574               if((charge1 < 0)&&(charge2 < 0)) {
 575                 fNnn[kQlong][iBin] += 1.;
 576                 fHistNN[kQlong]->Fill(fCentrality,qLong);
 577               }
 578               if((charge1 > 0)&&(charge2 < 0)) {
 579                 fNpn[kQlong][iBin] += 1.;
 580                 fHistPN[kQlong]->Fill(fCentrality,qLong);
 581               }
 582               if((charge1 < 0)&&(charge2 > 0)) {
 583                 fNpn[kQlong][iBin] += 1.;
 584                 fHistPN[kQlong]->Fill(fCentrality,qLong);
 585               }
 586             }//BF binning check
 587           }//p2 interval check
 588         }//p1 interval check
 589
 590         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 591           if( qOut > fP2Start[kQout] && qOut < fP2Stop[kQout]){
 592             iBin = Int_t((qOut-fP2Start[kQout])/fP2Step[kQout]);
 593             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 594               if((charge1 > 0)&&(charge2 > 0)) {
 595                 fNpp[kQout][iBin] += 1.;
 596                 fHistPP[kQout]->Fill(fCentrality,qOut);
 597                   }
 598               if((charge1 < 0)&&(charge2 < 0)) {
 599                 fNnn[kQout][iBin] += 1.;
 600                 fHistNN[kQout]->Fill(fCentrality,qOut);
 601               }
 602               if((charge1 > 0)&&(charge2 < 0)) {
 603                 fNpn[kQout][iBin] += 1.;
 604                 fHistPN[kQout]->Fill(fCentrality,qOut);
 605               }
 606               if((charge1 < 0)&&(charge2 > 0)) {
 607                 fNpn[kQout][iBin] += 1.;
 608                 fHistPN[kQout]->Fill(fCentrality,qOut);
 609               }
 610             }//BF binning check
 611           }//p2 interval check
 612         }//p1 interval check
 613
 614         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 615           if( qSide > fP2Start[kQside] && qSide < fP2Stop[kQside]){
 616             iBin = Int_t((qSide-fP2Start[kQside])/fP2Step[kQside]);
 617             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 618               if((charge1 > 0)&&(charge2 > 0)) {
 619                 fNpp[kQside][iBin] += 1.;
 620                 fHistPP[kQside]->Fill(fCentrality,qSide);
 621               }
 622               if((charge1 < 0)&&(charge2 < 0)) {
 623                 fNnn[kQside][iBin] += 1.;
 624                 fHistNN[kQside]->Fill(fCentrality,qSide);
 625               }
 626               if((charge1 > 0)&&(charge2 < 0)) {
 627                 fNpn[kQside][iBin] += 1.;
 628                 fHistPN[kQside]->Fill(fCentrality,qSide);
 629               }
 630               if((charge1 < 0)&&(charge2 > 0)) {
 631                 fNpn[kQside][iBin] += 1.;
 632                 fHistPN[kQside]->Fill(fCentrality,qSide);
 633                           }
 634             }//BF binning check
 635           }//p2 interval check
 636         }//p1 interval check
 637
 638         if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
 639           if( qInv > fP2Start[kQinv] && qInv < fP2Stop[kQinv]){
 640             iBin = Int_t((qInv-fP2Start[kQinv])/fP2Step[kQinv]);
 641             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 642               if((charge1 > 0)&&(charge2 > 0)) {
 643                 fNpp[kQinv][iBin] += 1.;
 644                 fHistPP[kQinv]->Fill(fCentrality,qInv);
 645               }
 646               if((charge1 < 0)&&(charge2 < 0)) {
 647                 fNnn[kQinv][iBin] += 1.;
 648                 fHistNN[kQinv]->Fill(fCentrality,qInv);
 649               }
 650               if((charge1 > 0)&&(charge2 < 0)) {
 651                 fNpn[kQinv][iBin] += 1.;
 652                 fHistPN[kQinv]->Fill(fCentrality,qInv);
 653               }
 654               if((charge1 < 0)&&(charge2 > 0)) {
 655                 fNpn[kQinv][iBin] += 1.;
 656                 fHistPN[kQinv]->Fill(fCentrality,qInv);
 657               }
 658             }//BF binning check
 659           }//p2 interval check
 660         }//p1 interval check
 661
 662         // Phi
 663         if((phi1 >= fP1Start[kPhi]) && (phi1 <= fP1Stop[kPhi]) && (phi2 >= fP1Start[kPhi]) && (phi2 <= fP1Stop[kPhi])) {
 664           if( dphi > fP2Start[kPhi] && dphi < fP2Stop[kPhi]){
 665             iBin = Int_t((dphi-fP2Start[kPhi])/fP2Step[kPhi]);
 666             if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
 667               if((charge1 > 0)&&(charge2 > 0)) {
 668                 fNpp[kPhi][iBin] += 1.;
 669                 fHistPP[kPhi]->Fill(fCentrality,dphi);
 670               }
 671               if((charge1 < 0)&&(charge2 < 0)) {
 672                 fNnn[kPhi][iBin] += 1.;
 673                 fHistNN[kPhi]->Fill(fCentrality,dphi);
 674               }
 675               if((charge1 > 0)&&(charge2 < 0)) {
 676                 fNpn[kPhi][iBin] += 1.;
 677                 fHistPN[kPhi]->Fill(fCentrality,dphi);
 678               }
 679               if((charge1 < 0)&&(charge2 > 0)) {
 680                 fNpn[kPhi][iBin] += 1.;
 681                 fHistPN[kPhi]->Fill(fCentrality,dphi);
 682               }
 683             }//BF binning check
 684           }//p2 interval check
 685         }//p1 interval check
 686     }//end of 2nd particle loop
 687   }//end of 1st particle loop
 688   //Printf("Number of analyzed events: %i",fAnalyzedEvents);
 689   //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]);
 690 }
 691
 692
 693 //____________________________________________________________________//
 694 Double_t AliBalance::GetBalance(Int_t iAnalysisType, Int_t p2) {
 695   // Returns the value of the balance function in bin p2
 696   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];
 697
 698   return fB[iAnalysisType][p2];
 699 }
 700
 701 //____________________________________________________________________//
 702 Double_t AliBalance::GetError(Int_t iAnalysisType, Int_t p2) {
 703   // Returns the error on the BF value for bin p2
 704   // The errors for fNn and fNp are neglected here (0.1 % of total error)
 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])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
 708                               + Double_t(fNnp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
 709                               //+ TMath::Power(fNpn[iAnalysisType][p2]-fNpp[iAnalysisType][p2],2)/TMath::Power(Double_t(fNp[iAnalysisType]),3)
 710                               //+ TMath::Power(fNnp[iAnalysisType][p2]-fNnn[iAnalysisType][p2],2)/TMath::Power(Double_t(fNn[iAnalysisType]),3)
 711                                ) /fP2Step[iAnalysisType];*/
 712
 713   ferror[iAnalysisType][p2] = TMath::Sqrt( Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType])) +
 714                                            Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType])) +
 715                                            Double_t(fNpn[iAnalysisType][p2])*TMath::Power((0.5/Double_t(fNp[iAnalysisType]) + 0.5/Double_t(fNn[iAnalysisType])),2))/fP2Step[iAnalysisType];
 716
 717   return ferror[iAnalysisType][p2];
 718 }
 719 //____________________________________________________________________//
 720 TGraphErrors *AliBalance::DrawBalance(Int_t iAnalysisType) {
 721
 722   // Draws the BF
 723   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
 724   Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
 725   Double_t b[MAXIMUM_NUMBER_OF_STEPS];
 726   Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
 727
 728   if((fNp[iAnalysisType] == 0)||(fNn[iAnalysisType] == 0)) {
 729     cerr<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
 730     return NULL;
 731   }
 732
 733   for(Int_t i = 0; i < fNumberOfBins[iAnalysisType]; i++) {
 734     b[i] = GetBalance(iAnalysisType,i);
 735     ber[i] = GetError(iAnalysisType,i);
 736     x[i] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
 737     xer[i] = 0.0;
 738   }
 739
 740   TGraphErrors *gr = new TGraphErrors(fNumberOfBins[iAnalysisType],x,b,xer,ber);
 741   gr->GetXaxis()->SetTitleColor(1);
 742   if(iAnalysisType==0) {
 743     gr->SetTitle("Balance function B(#Delta y)");
 744     gr->GetXaxis()->SetTitle("#Delta y");
 745     gr->GetYaxis()->SetTitle("B(#Delta y)");
 746   }
 747   if(iAnalysisType==1) {
 748     gr->SetTitle("Balance function B(#Delta #eta)");
 749     gr->GetXaxis()->SetTitle("#Delta #eta");
 750     gr->GetYaxis()->SetTitle("B(#Delta #eta)");
 751   }
 752   if(iAnalysisType==2) {
 753     gr->SetTitle("Balance function B(q_{long})");
 754     gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
 755     gr->GetYaxis()->SetTitle("B(q_{long}) ((GeV/c)^{-1})");
 756   }
 757   if(iAnalysisType==3) {
 758     gr->SetTitle("Balance function B(q_{out})");
 759     gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
 760     gr->GetYaxis()->SetTitle("B(q_{out}) ((GeV/c)^{-1})");
 761   }
 762   if(iAnalysisType==4) {
 763     gr->SetTitle("Balance function B(q_{side})");
 764     gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
 765     gr->GetYaxis()->SetTitle("B(q_{side}) ((GeV/c)^{-1})");
 766   }
 767   if(iAnalysisType==5) {
 768     gr->SetTitle("Balance function B(q_{inv})");
 769     gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 770     gr->GetYaxis()->SetTitle("B(q_{inv}) ((GeV/c)^{-1})");
 771   }
 772   if(iAnalysisType==6) {
 773     gr->SetTitle("Balance function B(#Delta #phi)");
 774     gr->GetXaxis()->SetTitle("#Delta #phi");
 775     gr->GetYaxis()->SetTitle("B(#Delta #phi)");
 776   }
 777
 778   return gr;
 779 }
 780
 781 //____________________________________________________________________//
 782 void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
 783   //Prints the calculated width of the BF and its error
 784   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
 785   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
 786   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
 787   Double_t deltaBalP2 = 0.0, integral = 0.0;
 788   Double_t deltaErrorNew = 0.0;
 789
 790   // cout<<"=================================================="<<endl;
 791   // for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) {
 792   //   x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
 793   //   cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
 794   // }
 795   // cout<<"=================================================="<<endl;
 796   for(Int_t i = 2; i <= fNumberOfBins[iAnalysisType]; i++) {
 797     gSumXi += gHistBalance->GetBinCenter(i);
 798     gSumBi += gHistBalance->GetBinContent(i);
 799     gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
 800     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
 801     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
 802     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
 803     gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
 804
 805     deltaBalP2 += fP2Step[iAnalysisType]*TMath::Power(gHistBalance->GetBinError(i),2);
 806     integral += fP2Step[iAnalysisType]*gHistBalance->GetBinContent(i);
 807   }
 808   for(Int_t i = 1; i < fNumberOfBins[iAnalysisType]; i++)
 809     deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
 810
 811   Double_t integralError = TMath::Sqrt(deltaBalP2);
 812   integralError *= 1.0;
 813
 814   Double_t delta = gSumBiXi / gSumBi; delta *= 1.0;
 815   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
 816   deltaError *= 1.0;
 817   // cout<<"Analysis type: "<<kBFAnalysisType[iAnalysisType].Data()<<endl;
 818   // cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 819   // cout<<"New error: "<<deltaErrorNew<<endl;
 820   // cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
 821   // cout<<"=================================================="<<endl;
 822 }
 823
 824 //____________________________________________________________________//
 825 TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow,Bool_t correctWithEfficiency, Bool_t correctWithAcceptanceOnly, Bool_t correctWithMixed, TH1D *hMixed[4]) {
 826   //Returns the BF histogram, extracted from the 6 TH2D objects
 827   //(private members) of the AliBalance class.
 828   //
 829   // Acceptance correction:
 830   // - only for analysis type = kEta
 831   // - only if etaWindow > 0 (default = -1.)
 832   // - calculated as proposed by STAR
 833   //
 834   TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
 835   TString histName = "gHistBalanceFunctionHistogram";
 836   histName += gAnalysisType[iAnalysisType];
 837
 838   SetInterval(iAnalysisType, fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
 839               fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
 840               fHistPP[iAnalysisType]->GetNbinsY(),
 841               fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),
 842               fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
 843
 844   // determine the projection thresholds
 845   Int_t binMinX, binMinY, binMinZ;
 846   Int_t binMaxX, binMaxY, binMaxZ;
 847
 848   fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMin),binMinX,binMinY,binMinZ);
 849   fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMax),binMaxX,binMaxY,binMaxZ);
 850
 851   TH1D *gHistBalanceFunctionHistogram = new TH1D(histName.Data(),"",fHistPP[iAnalysisType]->GetNbinsY(),fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
 852   switch(iAnalysisType) {
 853   case kRapidity:
 854     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
 855     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
 856     break;
 857   case kEta:
 858     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
 859     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
 860     break;
 861   case kQlong:
 862     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
 863     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
 864     break;
 865   case kQout:
 866     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
 867     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
 868     break;
 869   case kQside:
 870     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
 871     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
 872     break;
 873   case kQinv:
 874     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 875     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
 876     break;
 877   case kPhi:
 878     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
 879     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
 880     break;
 881   default:
 882     break;
 883   }
 884
 885   TH1D *hTemp1 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 886   TH1D *hTemp2 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f_copy",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 887   TH1D *hTemp3 = dynamic_cast<TH1D *>(fHistNN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistNN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 888   TH1D *hTemp4 = dynamic_cast<TH1D *>(fHistPP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 889   TH1D *hTemp5 = dynamic_cast<TH1D *>(fHistN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 890   TH1D *hTemp6 = dynamic_cast<TH1D *>(fHistP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
 891
 892   // get the file with the efficiency matrices
 893   // withAcceptanceOnly: Data single distributions are normalized to 1 (efficiency not taken into account)
 894   // else : Data single distributions are normalized to give single particle efficiency of MC
 895   TFile *fEfficiencyMatrix = NULL;
 896   if(correctWithEfficiency || correctWithMixed){
 897     if(correctWithAcceptanceOnly) fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/accOnlyFromConvolutionAllCent.root");
 898     else  fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/effFromConvolutionAllCent.root");
 899     if(!fEfficiencyMatrix){
 900       AliError("Efficiency histogram file not found");
 901       return NULL;
 902     }
 903   }
 904
 905   // do correction with the efficiency calculated from MC + Data (for single particles and two particle correlations)
 906   // - single particle efficiencies from MC (AliAnalysiTaskEfficiency)
 907   // - two particle efficiencies from convolution of data single particle distributions
 908   //   (normalized to single particle efficiency)
 909   if(iAnalysisType == kEta && etaWindow > 0 && correctWithEfficiency && !correctWithMixed){
 910
 911     TH1F* hEffP  = NULL;
 912     TH1F* hEffN  = NULL;
 913     TH1F* hEffPP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
 914     TH1F* hEffNN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax));
 915     TH1F* hEffPN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax));
 916
 917     // take the data distributions
 918     if(correctWithAcceptanceOnly){
 919       hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_Data",centrMin,centrMax));
 920       hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_Data",centrMin,centrMax));
 921     }
 922     // take the MC distributions
 923     else{
 924       hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_MC",centrMin,centrMax));
 925       hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_MC",centrMin,centrMax));
 926     }
 927
 928     if( !hEffP || !hEffN || !hEffPP || !hEffNN || !hEffPN){
 929       AliError(Form("Efficiency (eta) histograms not found: etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
 930       return NULL;
 931     }
 932
 933     for(Int_t iBin = 0; iBin < hEffP->GetNbinsX(); iBin++){
 934       hTemp5->SetBinError(iBin+1,hTemp5->GetBinError(iBin+1)/hEffN->GetBinContent(hEffN->FindBin(hTemp5->GetBinCenter(iBin+1))));
 935       hTemp5->SetBinContent(iBin+1,hTemp5->GetBinContent(iBin+1)/hEffN->GetBinContent(hEffN->FindBin(hTemp5->GetBinCenter(iBin+1))));
 936
 937       hTemp6->SetBinError(iBin+1,hTemp6->GetBinError(iBin+1)/hEffP->GetBinContent(hEffP->FindBin(hTemp6->GetBinCenter(iBin+1))));
 938       hTemp6->SetBinContent(iBin+1,hTemp6->GetBinContent(iBin+1)/hEffP->GetBinContent(hEffP->FindBin(hTemp6->GetBinCenter(iBin+1))));
 939     }
 940
 941     for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
 942
 943       hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp1->GetBinCenter(iBin+1))));
 944       hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp1->GetBinCenter(iBin+1))));
 945       hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp2->GetBinCenter(iBin+1))));
 946       hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp2->GetBinCenter(iBin+1))));
 947       hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/hEffNN->GetBinContent(hEffNN->FindBin(hTemp3->GetBinCenter(iBin+1))));
 948       hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/hEffNN->GetBinContent(hEffNN->FindBin(hTemp3->GetBinCenter(iBin+1))));
 949       hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/hEffPP->GetBinContent(hEffPP->FindBin(hTemp4->GetBinCenter(iBin+1))));
 950       hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/hEffPP->GetBinContent(hEffPP->FindBin(hTemp4->GetBinCenter(iBin+1))));
 951
 952     }
 953
 954     // TF1 *fPP = new TF1("fPP","pol1",0,1.6);  // phase space factor + efficiency for ++
 955     // fPP->SetParameters(0.736466,-0.461529);
 956     // TF1 *fNN = new TF1("fNN","pol1",0,1.6);  // phase space factor + efficiency for --
 957     // fNN->SetParameters(0.718616,-0.450473);
 958     // TF1 *fPN = new TF1("fPN","pol1",0,1.6);  // phase space factor + efficiency for +-
 959     // fPN->SetParameters(0.727507,-0.455981);
 960
 961     // for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
 962     //   hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1)));
 963     //   hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1)));
 964     //   hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1)));
 965     //   hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1)));
 966     //   hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/fNN->Eval(hTemp1->GetBinCenter(iBin+1)));
 967     //   hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/fNN->Eval(hTemp1->GetBinCenter(iBin+1)));
 968     //   hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/fPP->Eval(hTemp1->GetBinCenter(iBin+1)));
 969     //   hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/fPP->Eval(hTemp1->GetBinCenter(iBin+1)));
 970     // }
 971   }
 972
 973   // do correction with the efficiency calculated from MC + Data (for single particles and two particle correlations)
 974   // - single particle efficiencies from MC (AliAnalysiTaskEfficiency)
 975   // - two particle efficiencies from convolution of data single particle distributions
 976   //   (normalized to single particle efficiency)
 977   if(iAnalysisType == kPhi && correctWithEfficiency && !correctWithMixed){
 978
 979     TH1F* hEffPhiP  = NULL;
 980     TH1F* hEffPhiN  = NULL;
 981     TH1F* hEffPhiPP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
 982     TH1F* hEffPhiNN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax));
 983     TH1F* hEffPhiPN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax));
 984
 985     // take the data distributions
 986     if(correctWithAcceptanceOnly){
 987       hEffPhiP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffP_Cent%.0f-%.0f_Data",centrMin,centrMax));
 988       hEffPhiN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffN_Cent%.0f-%.0f_Data",centrMin,centrMax));
 989     }
 990     // take the MC distributions
 991     else{
 992       hEffPhiP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffP_Cent%.0f-%.0f_MC",centrMin,centrMax));
 993       hEffPhiN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffN_Cent%.0f-%.0f_MC",centrMin,centrMax));
 994     }
 995
 996     if( !hEffPhiP || !hEffPhiN || !hEffPhiPP || !hEffPhiNN || !hEffPhiPN){
 997       AliError("Efficiency (phi) histograms not found");
 998       return NULL;
 999     }
1000
1001     for(Int_t iBin = 0; iBin < hEffPhiP->GetNbinsX(); iBin++){
1002       hTemp5->SetBinError(iBin+1,hTemp5->GetBinError(iBin+1)/hEffPhiN->GetBinContent(hEffPhiN->FindBin(hTemp5->GetBinCenter(iBin+1))));
1003       hTemp5->SetBinContent(iBin+1,hTemp5->GetBinContent(iBin+1)/hEffPhiN->GetBinContent(hEffPhiN->FindBin(hTemp5->GetBinCenter(iBin+1))));
1004
1005       hTemp6->SetBinError(iBin+1,hTemp6->GetBinError(iBin+1)/hEffPhiP->GetBinContent(hEffPhiP->FindBin(hTemp6->GetBinCenter(iBin+1))));
1006       hTemp6->SetBinContent(iBin+1,hTemp6->GetBinContent(iBin+1)/hEffPhiP->GetBinContent(hEffPhiP->FindBin(hTemp6->GetBinCenter(iBin+1))));
1007     }
1008
1009     for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
1010
1011       hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp1->GetBinCenter(iBin+1))));
1012       hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp1->GetBinCenter(iBin+1))));
1013       hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp2->GetBinCenter(iBin+1))));
1014       hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp2->GetBinCenter(iBin+1))));
1015       hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/hEffPhiNN->GetBinContent(hEffPhiNN->FindBin(hTemp3->GetBinCenter(iBin+1))));
1016       hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/hEffPhiNN->GetBinContent(hEffPhiNN->FindBin(hTemp3->GetBinCenter(iBin+1))));
1017       hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/hEffPhiPP->GetBinContent(hEffPhiPP->FindBin(hTemp4->GetBinCenter(iBin+1))));
1018       hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/hEffPhiPP->GetBinContent(hEffPhiPP->FindBin(hTemp4->GetBinCenter(iBin+1))));
1019
1020     }
1021   }
1022
1023   // do the correction with the event mixing directly!
1024   if(correctWithMixed){
1025
1026     // take the MC distributions (for average efficiency)
1027     TH1F* hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_MC",centrMin,centrMax));
1028     TH1F* hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_MC",centrMin,centrMax));
1029
1030     TH1F* hEffPP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
1031     TH1F* hEffNN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax));
1032     TH1F* hEffPN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax));
1033
1034     if( !hEffP || !hEffN){
1035       AliError(Form("Efficiency (eta) histograms not found: etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
1036       return NULL;
1037     }
1038
1039     if(hMixed[0] && hMixed[1] && hMixed[2] && hMixed[3]){
1040
1041       // scale to average efficiency in the pt region (0.3-1.5) and |eta| < 0.8
1042       // by multiplying the average single particle efficiencies from HIJING
1043       // here we assume that the distributions are 1:
1044       // - in the integral for dphi (for averaging over sector structure)
1045       // - in the maximum for deta
1046       Double_t normPMC = (Double_t)hEffP->Integral()/(Double_t)hEffP->GetNbinsX();
1047       Double_t normNMC = (Double_t)hEffN->Integral()/(Double_t)hEffN->GetNbinsX();
1048       Double_t normPPMC = (Double_t)hEffPP->Integral()/(Double_t)hEffPP->GetNbinsX();
1049       Double_t normNNMC = (Double_t)hEffNN->Integral()/(Double_t)hEffNN->GetNbinsX();
1050       Double_t normPNMC = (Double_t)hEffPN->Integral()/(Double_t)hEffPN->GetNbinsX();
1051
1052       hMixed[0]->Scale(normPNMC);
1053       hMixed[1]->Scale(normPNMC);
1054       hMixed[2]->Scale(normNNMC);
1055       hMixed[3]->Scale(normPPMC);
1056
1057       // divide by event mixing
1058       hTemp1->Divide(hMixed[0]);
1059       hTemp2->Divide(hMixed[1]);
1060       hTemp3->Divide(hMixed[2]);
1061       hTemp4->Divide(hMixed[3]);
1062
1063       // scale also single histograms with average efficiency
1064       hTemp5->Scale(1./normNMC);
1065       hTemp6->Scale(1./normPMC);
1066
1067     }
1068     else{
1069       AliError("Correction with EventMixing requested, but not all Histograms there!");
1070       return NULL;
1071     }
1072   }
1073
1074
1075   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
1076     hTemp1->Sumw2();
1077     hTemp2->Sumw2();
1078     hTemp3->Sumw2();
1079     hTemp4->Sumw2();
1080     hTemp1->Add(hTemp3,-2.);
1081     hTemp1->Scale(1./hTemp5->Integral());
1082     hTemp2->Add(hTemp4,-2.);
1083     hTemp2->Scale(1./hTemp6->Integral());
1084     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1085     gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]);
1086   }
1087
1088   // do the acceptance correction (only for Eta and etaWindow > 0)
1089   if(iAnalysisType == kEta && etaWindow > 0 && !correctWithEfficiency && !correctWithMixed){
1090     for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
1091
1092       Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
1093       Double_t corrected    = notCorrected / (1 - (gHistBalanceFunctionHistogram->GetBinCenter(iBin+1))/ etaWindow );
1094       gHistBalanceFunctionHistogram->SetBinContent(iBin+1, corrected);
1095       gHistBalanceFunctionHistogram->SetBinError(iBin+1,corrected/notCorrected*gHistBalanceFunctionHistogram->GetBinError(iBin+1));
1096
1097     }
1098   }
1099
1100   if(fEfficiencyMatrix)   fEfficiencyMatrix->Close();
1101
1102   PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);
1103
1104   return gHistBalanceFunctionHistogram;
1105 }