PWGCF/EBYE/BalanceFunctions/AliBalancePsi.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: AliBalancePsi.cxx 54125 2012-01-24 21:07:41Z miweber $ */
  15
  16 //-----------------------------------------------------------------
  17 //           Balance Function class
  18 //   This is the class to deal with the Balance Function wrt Psi analysis
  19 //   Origin: Panos Christakoglou, Nikhef, Panos.Christakoglou@cern.ch
  20 //-----------------------------------------------------------------
  21
  22
  23 //ROOT
  24 #include <Riostream.h>
  25 #include <TMath.h>
  26 #include <TAxis.h>
  27 #include <TH2D.h>
  28 #include <TH3D.h>
  29 #include <TLorentzVector.h>
  30 #include <TObjArray.h>
  31 #include <TGraphErrors.h>
  32 #include <TString.h>
  33
  34 #include "AliVParticle.h"
  35 #include "AliMCParticle.h"
  36 #include "AliESDtrack.h"
  37 #include "AliAODTrack.h"
  38 #include "AliTHn.h"
  39
  40 #include "AliBalancePsi.h"
  41
  42 ClassImp(AliBalancePsi)
  43
  44 //____________________________________________________________________//
  45 AliBalancePsi::AliBalancePsi() :
  46   TObject(),
  47   bShuffle(kFALSE),
  48   fAnalysisLevel("ESD"),
  49   fAnalyzedEvents(0) ,
  50   fCentralityId(0) ,
  51   fCentStart(0.),
  52   fCentStop(0.),
  53   fHistP(0),
  54   fHistN(0),
  55   fHistPN(0),
  56   fHistNP(0),
  57   fHistPP(0),
  58   fHistNN(0),
  59   fPsiInterval(15.) {
  60   // Default constructor
  61 }
  62
  63 //____________________________________________________________________//
  64 AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
  65   TObject(balance), bShuffle(balance.bShuffle),
  66   fAnalysisLevel(balance.fAnalysisLevel),
  67   fAnalyzedEvents(balance.fAnalyzedEvents),
  68   fCentralityId(balance.fCentralityId),
  69   fCentStart(balance.fCentStart),
  70   fCentStop(balance.fCentStop),
  71   fHistP(balance.fHistP),
  72   fHistN(balance.fHistN),
  73   fHistPN(balance.fHistPN),
  74   fHistNP(balance.fHistNP),
  75   fHistPP(balance.fHistPP),
  76   fHistNN(balance.fHistNN),
  77   fPsiInterval(balance.fPsiInterval) {
  78   //copy constructor
  79 }
  80
  81 //____________________________________________________________________//
  82 AliBalancePsi::~AliBalancePsi() {
  83   // Destructor
  84   delete fHistP;
  85   delete fHistN;
  86   delete fHistPN;
  87   delete fHistNP;
  88   delete fHistPP;
  89   delete fHistNN;
  90 }
  91
  92 //____________________________________________________________________//
  93 void AliBalancePsi::InitHistograms() {
  94   // single particle histograms
  95   Int_t anaSteps   = 1;       // analysis steps
  96   Int_t iBinSingle[nTrackVariablesSingle];        // binning for track variables
  97   Double_t* dBinsSingle[nTrackVariablesSingle];   // bins for track variables
  98   TString axisTitleSingle[nTrackVariablesSingle]; // axis titles for track variables
  99
 100   // two particle histograms
 101   Int_t iBinPair[nTrackVariablesPair];         // binning for track variables
 102   Double_t* dBinsPair[nTrackVariablesPair];    // bins for track variables
 103   TString axisTitlePair[nTrackVariablesPair];  // axis titles for track variables
 104
 105   //centrality
 106   const Int_t kNCentralityBins = 9;
 107   Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
 108   //const Int_t kNCentralityBins = 200;
 109   //Double_t centralityBins[kNCentralityBins+1];
 110   //for(Int_t i = 0; i < kNCentralityBins+1; i++)
 111   //centralityBins[i] = 0.0 + i * 0.5;
 112   iBinSingle[0]       = kNCentralityBins;
 113   dBinsSingle[0]      = centralityBins;
 114   axisTitleSingle[0]  = "Centrality percentile [%]";
 115   iBinPair[0]       = kNCentralityBins;
 116   dBinsPair[0]      = centralityBins;
 117   axisTitlePair[0]  = "Centrality percentile [%]";
 118
 119   //Psi_2
 120   const Int_t kNPsi2Bins = 24;
 121   Double_t psi2Bins[kNPsi2Bins+1];
 122   for(Int_t i = 0; i < kNPsi2Bins+1; i++)
 123     psi2Bins[i] = 0.0 + i * 7.5;
 124   iBinSingle[1]       = kNPsi2Bins;
 125   dBinsSingle[1]      = psi2Bins;
 126   axisTitleSingle[1]  = "#phi - #Psi_{2} (#circ)";
 127   iBinPair[1]       = kNPsi2Bins;
 128   dBinsPair[1]      = psi2Bins;
 129   axisTitlePair[1]  = "#phi - #Psi_{2} (#circ)";
 130
 131   // eta
 132   const Int_t kNEtaBins = 10;
 133   Double_t etaBins[kNEtaBins+1];
 134   for(Int_t i = 0; i < kNEtaBins+1; i++)
 135     etaBins[i] = -1.0 + i * 0.2;
 136   iBinSingle[2]       = kNEtaBins;
 137   dBinsSingle[2]      = etaBins;
 138   axisTitleSingle[2]  = "#eta";
 139
 140   //#eta of triggered particle
 141   /*iBinPair[2]       = kNEtaBins;
 142   dBinsPair[2]      = etaBins;
 143   axisTitlePair[2]  = "#eta"; */
 144
 145   // delta eta
 146   const Int_t kNDeltaEtaBins = 80;
 147   Double_t deltaEtaBins[kNDeltaEtaBins+1];
 148   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
 149     deltaEtaBins[i] = -2.0 + i * 0.05;
 150   iBinPair[2]       = kNDeltaEtaBins;
 151   dBinsPair[2]      = deltaEtaBins;
 152   axisTitlePair[2]  = "#Delta #eta";
 153
 154   // y
 155   /*const Int_t kNYBins = 40;
 156   Double_t yBins[kNYBins+1];
 157   for(Int_t i = 0; i < kNYBins+1; i++)
 158     yBins[i] = -1.0 + i * 0.05;
 159   iBinSingle[3]       = kNYBins;
 160   dBinsSingle[3]      = yBins;
 161   axisTitleSingle[3]  = "y";
 162
 163   //y of triggered particle
 164   iBinPair[3]       = kNYBins;
 165   dBinsPair[3]      = yBins;
 166   axisTitlePair[3]  = "y"; */
 167
 168   // delta y
 169   /*const Int_t kNDeltaYBins = 40;
 170   Double_t deltaYBins[kNDeltaYBins+1];
 171   for(Int_t i = 0; i < kNDeltaYBins+1; i++)
 172     deltaYBins[i] = -2.0 + i * 0.1;
 173   iBinPair[8]       = kNDeltaYBins;
 174   dBinsPair[8]      = deltaYBins;
 175   axisTitlePair[8]  = "#Delta y"; */
 176
 177   // phi
 178   const Int_t kNPhiBins = 18;
 179   Double_t phiBins[kNPhiBins+1];
 180   for(Int_t i = 0; i < kNPhiBins+1; i++){
 181     phiBins[i] = 0.0 + i * 20.;
 182   }
 183   iBinSingle[3]       = kNPhiBins;
 184   dBinsSingle[3]      = phiBins;
 185   axisTitleSingle[3]  = "#phi (#circ)";
 186
 187   // delta phi
 188   const Int_t kNDeltaPhiBins = 72;
 189   Double_t deltaPhiBins[kNDeltaPhiBins+1];
 190   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
 191     deltaPhiBins[i] = -180.0 + i * 5.;
 192   }
 193   iBinPair[3]       = kNDeltaPhiBins;
 194   dBinsPair[3]      = deltaPhiBins;
 195   axisTitlePair[3]  = "#Delta #phi (#circ)";
 196
 197   // pt(trigger-associated)
 198   const Int_t kNPtBins = 50;
 199   Double_t ptBins[kNPtBins+1];
 200   for(Int_t i = 0; i < kNPtBins+1; i++){
 201     ptBins[i] = 0.0 + i * 0.5;
 202    }
 203   iBinSingle[4]       = kNPtBins;
 204   dBinsSingle[4]      = ptBins;
 205   axisTitleSingle[4]  = "p_{t}^{trig.} (GeV/c)";
 206
 207   iBinPair[4]       = kNPtBins;
 208   dBinsPair[4]      = ptBins;
 209   axisTitlePair[4]  = "p_{t}^{trig.} (GeV/c)";
 210
 211   iBinPair[5]       = kNPtBins;
 212   dBinsPair[5]      = ptBins;
 213   axisTitlePair[5]  = "p_{t}^{assoc.} (GeV/c)";
 214
 215   // qside
 216   /*const Int_t kNQSideBins = 200;
 217   Double_t qSideBins[kNQSideBins+1];
 218   for(Int_t i = 0; i < kNQSideBins+1; i++)
 219     qSideBins[i] = 0.0 + i * 0.02;
 220   iBinPair[10]       = kNQSideBins;
 221   dBinsPair[10]      = qSideBins;
 222   axisTitlePair[10]  = "q_{side} (GeV/c)";
 223
 224   // qout
 225   const Int_t kNQoutBins = 200;
 226   Double_t qoutBins[kNQoutBins+1];
 227   for(Int_t i = 0; i < kNQoutBins+1; i++)
 228     qoutBins[i] = 0.0 + i * 0.02;
 229   iBinPair[11]       = kNQoutBins;
 230   dBinsPair[11]      = qoutBins;
 231   axisTitlePair[11]  = "q_{out} (GeV/c)";
 232
 233   // qlong
 234   const Int_t kNQlongBins = 200;
 235   Double_t qlongBins[kNQlongBins+1];
 236   for(Int_t i = 0; i < kNQlongBins+1; i++)
 237     qlongBins[i] = 0.0 + i * 0.02;
 238   iBinPair[12]       = kNQlongBins;
 239   dBinsPair[12]      = qlongBins;
 240   axisTitlePair[12]  = "q_{long} (GeV/c)";
 241
 242   // qinv
 243   const Int_t kNQinvBins = 200;
 244   Double_t qinvBins[kNQinvBins+1];
 245   for(Int_t i = 0; i < kNQinvBins+1; i++)
 246     qinvBins[i] = 0.0 + i * 0.02;
 247   iBinPair[13]       = kNQinvBins;
 248   dBinsPair[13]      = qinvBins;
 249   axisTitlePair[13]  = "q_{inv} (GeV/c)";*/
 250
 251   TString histName;
 252   //+ triggered particles
 253   histName = "fHistP";
 254   if(bShuffle) histName.Append("_shuffle");
 255   if(fCentralityId) histName += fCentralityId.Data();
 256   fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,nTrackVariablesSingle,iBinSingle);
 257   for (Int_t j=0; j<nTrackVariablesSingle; j++) {
 258     fHistP->SetBinLimits(j, dBinsSingle[j]);
 259     fHistP->SetVarTitle(j, axisTitleSingle[j]);
 260   }
 261
 262   //- triggered particles
 263   histName = "fHistN";
 264   if(bShuffle) histName.Append("_shuffle");
 265   if(fCentralityId) histName += fCentralityId.Data();
 266   fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,nTrackVariablesSingle,iBinSingle);
 267   for (Int_t j=0; j<nTrackVariablesSingle; j++) {
 268     fHistN->SetBinLimits(j, dBinsSingle[j]);
 269     fHistN->SetVarTitle(j, axisTitleSingle[j]);
 270   }
 271
 272   //+- pairs
 273   histName = "fHistPN";
 274   if(bShuffle) histName.Append("_shuffle");
 275   if(fCentralityId) histName += fCentralityId.Data();
 276   fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
 277   for (Int_t j=0; j<nTrackVariablesPair; j++) {
 278     fHistPN->SetBinLimits(j, dBinsPair[j]);
 279     fHistPN->SetVarTitle(j, axisTitlePair[j]);
 280   }
 281
 282   //-+ pairs
 283   histName = "fHistNP";
 284   if(bShuffle) histName.Append("_shuffle");
 285   if(fCentralityId) histName += fCentralityId.Data();
 286   fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
 287   for (Int_t j=0; j<nTrackVariablesPair; j++) {
 288     fHistNP->SetBinLimits(j, dBinsPair[j]);
 289     fHistNP->SetVarTitle(j, axisTitlePair[j]);
 290   }
 291
 292   //++ pairs
 293   histName = "fHistPP";
 294   if(bShuffle) histName.Append("_shuffle");
 295   if(fCentralityId) histName += fCentralityId.Data();
 296   fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
 297   for (Int_t j=0; j<nTrackVariablesPair; j++) {
 298     fHistPP->SetBinLimits(j, dBinsPair[j]);
 299     fHistPP->SetVarTitle(j, axisTitlePair[j]);
 300   }
 301
 302   //-- pairs
 303   histName = "fHistNN";
 304   if(bShuffle) histName.Append("_shuffle");
 305   if(fCentralityId) histName += fCentralityId.Data();
 306   fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
 307   for (Int_t j=0; j<nTrackVariablesPair; j++) {
 308     fHistNN->SetBinLimits(j, dBinsPair[j]);
 309     fHistNN->SetVarTitle(j, axisTitlePair[j]);
 310   }
 311   Printf("Finished setting up the AliTHn");
 312 }
 313
 314 //____________________________________________________________________//
 315 void AliBalancePsi::CalculateBalance(Float_t fCentrality,
 316                                      Double_t gReactionPlane,
 317                                      vector<Double_t> **chargeVector) {
 318   // Calculates the balance function
 319   fAnalyzedEvents++;
 320   Int_t i = 0 , j = 0;
 321
 322   // Initialize histograms if not done yet
 323   if(!fHistPN){
 324     AliWarning("Histograms not yet initialized! --> Will be done now");
 325     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 326     InitHistograms();
 327   }
 328
 329   Double_t trackVarsSingle[nTrackVariablesSingle];
 330   Double_t trackVarsPair[nTrackVariablesPair];
 331
 332   Int_t gNtrack = chargeVector[0]->size();
 333   //Printf("(AliBalancePsi) Number of tracks: %d",gNtrack);
 334   Double_t gPsiMinusPhi = 0.;
 335   Double_t dy = 0., deta = 0.;
 336   Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
 337   Double_t dphi = 0.;
 338
 339   Double_t charge1  = 0;
 340   Double_t eta1 = 0., rap1 = 0.;
 341   Double_t px1 = 0., py1 = 0., pz1 = 0.;
 342   Double_t pt1 = 0.;
 343   Double_t energy1 = 0.;
 344   Double_t phi1    = 0.;
 345
 346   Double_t charge2  = 0;
 347   Double_t eta2 = 0., rap2 = 0.;
 348   Double_t px2 = 0., py2 = 0., pz2 = 0.;
 349   Double_t pt2 = 0.;
 350   Double_t energy2 = 0.;
 351   Double_t phi2    = 0.;
 352
 353   for(i = 0; i < gNtrack; i++) {
 354     charge1 = chargeVector[0]->at(i);
 355     rap1    = chargeVector[1]->at(i);
 356     eta1    = chargeVector[2]->at(i);
 357     phi1    = chargeVector[3]->at(i);
 358     px1     = chargeVector[4]->at(i);
 359     py1     = chargeVector[5]->at(i);
 360     pz1     = chargeVector[6]->at(i);
 361     pt1     = chargeVector[7]->at(i);
 362     energy1 = chargeVector[8]->at(i);
 363     gPsiMinusPhi   = TMath::Abs(phi1 - gReactionPlane);
 364     if(gPsiMinusPhi > 180.) gPsiMinusPhi = 360. - gPsiMinusPhi;
 365     //if(gPsiMinusPhi < -fPsiInterval/2) gPsiMinusPhi = 360. + gPsiMinusPhi;
 366
 367     trackVarsSingle[0]    =  fCentrality;
 368     trackVarsSingle[1]    =  gPsiMinusPhi;
 369     trackVarsSingle[2]    =  eta1;
 370     //trackVarsSingle[3]    =  rap1;
 371     trackVarsSingle[3]    =  phi1;
 372     trackVarsSingle[4]    =  pt1;
 373
 374     //Printf("Track(a) %d - phi-Psi: %lf",i+1,trackVarsSingle[1]);
 375     //fill single particle histograms
 376     if(charge1 > 0)
 377       fHistP->Fill(trackVarsSingle,0,1.);
 378     else if(charge1 < 0)
 379       fHistN->Fill(trackVarsSingle,0,1.);
 380
 381     for(j = 0; j < i; j++) {
 382       charge2 = chargeVector[0]->at(j);
 383       rap2    = chargeVector[1]->at(j);
 384       eta2    = chargeVector[2]->at(j);
 385       phi2    = chargeVector[3]->at(j);
 386       px2     = chargeVector[4]->at(j);
 387       py2     = chargeVector[5]->at(j);
 388       pz2     = chargeVector[6]->at(j);
 389       pt2     = chargeVector[7]->at(j);
 390       energy2 = chargeVector[8]->at(j);
 391       //Printf("Track(b) %d - pt: %lf",j+1,pt2);
 392
 393       // filling the arrays
 394       // RAPIDITY
 395       dy = rap1 - rap2;
 396
 397       // Eta
 398       deta = eta1 - eta2;
 399
 400       //qlong
 401       Double_t eTot = energy1 + energy2;
 402       Double_t pxTot = px1 + px2;
 403       Double_t pyTot = py1 + py2;
 404       Double_t pzTot = pz1 + pz2;
 405       Double_t q0Tot = energy1 - energy2;
 406       Double_t qxTot = px1 - px2;
 407       Double_t qyTot = py1 - py2;
 408       Double_t qzTot = pz1 - pz2;
 409
 410       Double_t eTot2 = eTot*eTot;
 411       Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
 412       Double_t pzTot2 = pzTot*pzTot;
 413
 414       Double_t q0Tot2 = q0Tot*q0Tot;
 415       Double_t qTot2  = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;
 416
 417       Double_t snn    = eTot2 - pTot2;
 418       Double_t ptTot2 = pTot2 - pzTot2 ;
 419       Double_t ptTot  = TMath::Sqrt( ptTot2 );
 420
 421       qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
 422
 423       //qout
 424       qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
 425
 426       //qside
 427       qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
 428
 429       //qinv
 430       qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
 431
 432       //phi
 433       dphi = phi2 - phi1;
 434       if(dphi < -180.) dphi += 360.;  //dphi should be between -180 and 180!
 435       else if(dphi > 180.) dphi -= 360.; //dphi should be between -180 and 180!
 436
 437       trackVarsPair[0]    =  fCentrality;
 438       trackVarsPair[1]    =  gPsiMinusPhi;
 439       //trackVarsPair[2]    =  eta1;
 440       //trackVarsPair[3]    =  rap1;
 441       //trackVarsPair[4]    =  phi1;
 442       trackVarsPair[2]    =  deta;
 443       //trackVarsPair[8]    =  dy;
 444       trackVarsPair[3]    =  dphi;
 445       trackVarsPair[4]    =  pt1;
 446       trackVarsPair[5]    =  pt2;
 447       //trackVarsPair[10]   =  qSide;
 448       //trackVarsPair[11]   =  qOut;
 449       //trackVarsPair[12]   =  qLong;
 450       //trackVarsPair[13]   =  qInv;
 451
 452       if( charge1 > 0 && charge2 < 0)  fHistPN->Fill(trackVarsPair,0,1.);
 453       else if( charge1 < 0 && charge2 > 0)  fHistNP->Fill(trackVarsPair,0,1.);
 454       else if( charge1 > 0 && charge2 > 0)  fHistPP->Fill(trackVarsPair,0,1.);
 455       else if( charge1 < 0 && charge2 < 0)  fHistNN->Fill(trackVarsPair,0,1.);
 456       else AliWarning("Wrong charge combination!");
 457
 458     }//end of 2nd particle loop
 459   }//end of 1st particle loop
 460 }
 461
 462 //____________________________________________________________________//
 463 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
 464                                                  Int_t iVariablePair,
 465                                                  Double_t centrMin,
 466                                                  Double_t centrMax,
 467                                                  Double_t psiMin,
 468                                                  Double_t psiMax) {
 469   //Returns the BF histogram, extracted from the 6 AliTHn objects
 470   //(private members) of the AliBalancePsi class.
 471   //iVariableSingle: 2(eta), 3(y), 4(phi)
 472   //iVariablePair: 2(Delta eta), 3(Delta y), 4(Delta phi), 5(qside), 6(qout), 7(qlong) 8(qinv)
 473   TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","phi","qlong","qout","qside","qinv"};
 474   TString histName = "gHistBalanceFunctionHistogram";
 475   histName += gAnalysisType[iVariablePair];
 476
 477   // centrality
 478   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 479   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 480   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 481   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 482   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 483   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 484
 485   // Psi_2
 486   fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 487   fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 488   fHistPN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 489   fHistNP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 490   fHistPP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 491   fHistNN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 492
 493   //Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0));
 494
 495   // Project into the wanted space (1st: analysis step, 2nd: axis)
 496   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
 497   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
 498   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
 499   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
 500   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
 501   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
 502
 503   TH1D *gHistBalanceFunctionHistogram = 0x0;
 504   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 505     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
 506     gHistBalanceFunctionHistogram->Reset();
 507
 508     switch(iVariablePair) {
 509     case 2:
 510       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
 511       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
 512       break;
 513     case 3:
 514       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
 515       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
 516       break;
 517     case 4:
 518       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
 519       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
 520       break;
 521     case 5:
 522       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
 523       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
 524       break;
 525     case 6:
 526       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
 527       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
 528       break;
 529     case 7:
 530       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
 531       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
 532       break;
 533     case 8:
 534       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 535       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
 536       break;
 537     default:
 538       break;
 539     }
 540
 541     hTemp1->Sumw2();
 542     hTemp2->Sumw2();
 543     hTemp3->Sumw2();
 544     hTemp4->Sumw2();
 545     hTemp1->Add(hTemp3,-1.);
 546     hTemp1->Scale(1./hTemp5->GetEntries());
 547     hTemp2->Add(hTemp4,-1.);
 548     hTemp2->Scale(1./hTemp6->GetEntries());
 549     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 550     gHistBalanceFunctionHistogram->Scale(0.5);
 551   }
 552
 553   return gHistBalanceFunctionHistogram;
 554 }
 555
 556 //____________________________________________________________________//
 557 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t centrMin,
 558                                               Double_t centrMax,
 559                                               Double_t psiMin,
 560                                               Double_t psiMax) {
 561   //Returns the 2D correlation function for (+-) pairs
 562   // centrality: axis 0
 563   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 564   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 565
 566   // Psi_2: axis 1
 567   fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 568   fHistPN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 569
 570   //0:step, 2: Delta eta, 3: Delta phi
 571   TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,2,3));
 572   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
 573   //Printf("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,2,3)->GetEntries()));
 574   if((Double_t)(fHistP->Project(0,2)->GetEntries())!=0)
 575     gHist->Scale(1./(Double_t)(fHistP->Project(0,2)->GetEntries()));
 576
 577   return gHist;
 578 }
 579
 580 //____________________________________________________________________//
 581 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t centrMin,
 582                                               Double_t centrMax,
 583                                               Double_t psiMin,
 584                                               Double_t psiMax) {
 585   //Returns the 2D correlation function for (+-) pairs
 586   // centrality: axis 0
 587   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 588   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 589
 590   // Psi_2: axis 1
 591   fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 592   fHistNP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 593
 594   //0:step, 2: Delta eta, 3: Delta phi
 595   TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,2,3));
 596   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
 597   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
 598   if((Double_t)(fHistN->Project(0,2)->GetEntries())!=0)
 599     gHist->Scale(1./(Double_t)(fHistN->Project(0,2)->GetEntries()));
 600
 601   return gHist;
 602 }
 603
 604 //____________________________________________________________________//
 605 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t centrMin,
 606                                               Double_t centrMax,
 607                                               Double_t psiMin,
 608                                               Double_t psiMax) {
 609   //Returns the 2D correlation function for (+-) pairs
 610   // centrality: axis 0
 611   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 612   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 613
 614   // Psi_2: axis 1
 615   fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 616   fHistPP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 617
 618   //0:step, 2: Delta eta, 3: Delta phi
 619   TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,2,3));
 620   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
 621   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
 622   if((Double_t)(fHistP->Project(0,2)->GetEntries())!=0)
 623     gHist->Scale(1./(Double_t)(fHistP->Project(0,2)->GetEntries()));
 624
 625   return gHist;
 626 }
 627
 628 //____________________________________________________________________//
 629 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t centrMin,
 630                                               Double_t centrMax,
 631                                               Double_t psiMin,
 632                                               Double_t psiMax) {
 633   //Returns the 2D correlation function for (+-) pairs
 634   // centrality: axis 0
 635   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 636   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax);
 637
 638   // Psi_2: axis 1
 639   fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 640   fHistNN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax);
 641
 642   //0:step, 2: Delta eta, 3: Delta phi
 643   TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,2,3));
 644   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
 645   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
 646   if((Double_t)(fHistN->Project(0,2)->GetEntries())!=0)
 647     gHist->Scale(1./(Double_t)(fHistN->Project(0,2)->GetEntries()));
 648
 649   return gHist;
 650 }