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 <TCanvas.h>
  26 #include <TMath.h>
  27 #include <TAxis.h>
  28 #include <TH2D.h>
  29 #include <TH3D.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 #include "AliTHn.h"
  40
  41 #include "AliBalancePsi.h"
  42
  43 ClassImp(AliBalancePsi)
  44
  45 //____________________________________________________________________//
  46 AliBalancePsi::AliBalancePsi() :
  47   TObject(),
  48   bShuffle(kFALSE),
  49   fAnalysisLevel("ESD"),
  50   fAnalyzedEvents(0) ,
  51   fCentralityId(0) ,
  52   fCentStart(0.),
  53   fCentStop(0.),
  54   fHistP(0),
  55   fHistN(0),
  56   fHistPN(0),
  57   fHistNP(0),
  58   fHistPP(0),
  59   fHistNN(0),
  60   fPsiInterval(15.) {
  61   // Default constructor
  62 }
  63
  64 //____________________________________________________________________//
  65 AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
  66   TObject(balance), bShuffle(balance.bShuffle),
  67   fAnalysisLevel(balance.fAnalysisLevel),
  68   fAnalyzedEvents(balance.fAnalyzedEvents),
  69   fCentralityId(balance.fCentralityId),
  70   fCentStart(balance.fCentStart),
  71   fCentStop(balance.fCentStop),
  72   fHistP(balance.fHistP),
  73   fHistN(balance.fHistN),
  74   fHistPN(balance.fHistPN),
  75   fHistNP(balance.fHistNP),
  76   fHistPP(balance.fHistPP),
  77   fHistNN(balance.fHistNN),
  78   fPsiInterval(balance.fPsiInterval) {
  79   //copy constructor
  80 }
  81
  82 //____________________________________________________________________//
  83 AliBalancePsi::~AliBalancePsi() {
  84   // Destructor
  85   delete fHistP;
  86   delete fHistN;
  87   delete fHistPN;
  88   delete fHistNP;
  89   delete fHistPP;
  90   delete fHistNN;
  91 }
  92
  93 //____________________________________________________________________//
  94 void AliBalancePsi::InitHistograms() {
  95   // single particle histograms
  96   Int_t anaSteps   = 1;       // analysis steps
  97   Int_t iBinSingle[nTrackVariablesSingle];        // binning for track variables
  98   Double_t* dBinsSingle[nTrackVariablesSingle];   // bins for track variables
  99   TString axisTitleSingle[nTrackVariablesSingle]; // axis titles for track variables
 100
 101   // two particle histograms
 102   Int_t iBinPair[nTrackVariablesPair];         // binning for track variables
 103   Double_t* dBinsPair[nTrackVariablesPair];    // bins for track variables
 104   TString axisTitlePair[nTrackVariablesPair];  // axis titles for track variables
 105
 106   //centrality
 107   /*const Int_t kNCentralityBins = 9;
 108   Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
 109   iBinSingle[0]       = kNCentralityBins;
 110   dBinsSingle[0]      = centralityBins;
 111   axisTitleSingle[0]  = "Centrality percentile [%]";
 112   iBinPair[0]       = kNCentralityBins;
 113   dBinsPair[0]      = centralityBins;
 114   axisTitlePair[0]  = "Centrality percentile [%]"; */
 115
 116   //Psi_2
 117   const Int_t kNPsi2Bins = 3;
 118   Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5};
 119   iBinSingle[0]       = kNPsi2Bins;
 120   dBinsSingle[0]      = psi2Bins;
 121   axisTitleSingle[0]  = "#phi - #Psi_{2} (a.u.)";
 122   iBinPair[0]       = kNPsi2Bins;
 123   dBinsPair[0]      = psi2Bins;
 124   axisTitlePair[0]  = "#phi - #Psi_{2} (a.u.)";
 125
 126    // delta eta
 127   const Int_t kNDeltaEtaBins = 80;
 128   Double_t deltaEtaBins[kNDeltaEtaBins+1];
 129   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
 130     deltaEtaBins[i] = -2.0 + i * 0.05;
 131   iBinPair[1]       = kNDeltaEtaBins;
 132   dBinsPair[1]      = deltaEtaBins;
 133   axisTitlePair[1]  = "#Delta #eta";
 134
 135    // delta phi
 136   const Int_t kNDeltaPhiBins = 72;
 137   Double_t deltaPhiBins[kNDeltaPhiBins+1];
 138   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
 139     deltaPhiBins[i] = -180.0 + i * 5.;
 140   }
 141   iBinPair[2]       = kNDeltaPhiBins;
 142   dBinsPair[2]      = deltaPhiBins;
 143   axisTitlePair[2]  = "#Delta #phi (#circ)";
 144
 145   // pt(trigger-associated)
 146   const Int_t kNPtBins = 40;
 147   Double_t ptBins[kNPtBins+1];
 148   for(Int_t i = 0; i < kNPtBins+1; i++){
 149     ptBins[i] = 0.0 + i * 0.5;
 150    }
 151   iBinSingle[1]       = kNPtBins;
 152   dBinsSingle[1]      = ptBins;
 153   axisTitleSingle[1]  = "p_{t}^{trig.} (GeV/c)";
 154
 155   iBinPair[3]       = kNPtBins;
 156   dBinsPair[3]      = ptBins;
 157   axisTitlePair[3]  = "p_{t}^{trig.} (GeV/c)";
 158
 159   iBinPair[4]       = kNPtBins;
 160   dBinsPair[4]      = ptBins;
 161   axisTitlePair[4]  = "p_{t}^{assoc.} (GeV/c)";
 162
 163   TString histName;
 164   //+ triggered particles
 165   histName = "fHistP";
 166   if(bShuffle) histName.Append("_shuffle");
 167   if(fCentralityId) histName += fCentralityId.Data();
 168   fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,nTrackVariablesSingle,iBinSingle);
 169   for (Int_t j=0; j<nTrackVariablesSingle; j++) {
 170     fHistP->SetBinLimits(j, dBinsSingle[j]);
 171     fHistP->SetVarTitle(j, axisTitleSingle[j]);
 172   }
 173
 174   //- triggered particles
 175   histName = "fHistN";
 176   if(bShuffle) histName.Append("_shuffle");
 177   if(fCentralityId) histName += fCentralityId.Data();
 178   fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,nTrackVariablesSingle,iBinSingle);
 179   for (Int_t j=0; j<nTrackVariablesSingle; j++) {
 180     fHistN->SetBinLimits(j, dBinsSingle[j]);
 181     fHistN->SetVarTitle(j, axisTitleSingle[j]);
 182   }
 183
 184   //+- pairs
 185   histName = "fHistPN";
 186   if(bShuffle) histName.Append("_shuffle");
 187   if(fCentralityId) histName += fCentralityId.Data();
 188   fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
 189   for (Int_t j=0; j<nTrackVariablesPair; j++) {
 190     fHistPN->SetBinLimits(j, dBinsPair[j]);
 191     fHistPN->SetVarTitle(j, axisTitlePair[j]);
 192   }
 193
 194   //-+ pairs
 195   histName = "fHistNP";
 196   if(bShuffle) histName.Append("_shuffle");
 197   if(fCentralityId) histName += fCentralityId.Data();
 198   fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
 199   for (Int_t j=0; j<nTrackVariablesPair; j++) {
 200     fHistNP->SetBinLimits(j, dBinsPair[j]);
 201     fHistNP->SetVarTitle(j, axisTitlePair[j]);
 202   }
 203
 204   //++ pairs
 205   histName = "fHistPP";
 206   if(bShuffle) histName.Append("_shuffle");
 207   if(fCentralityId) histName += fCentralityId.Data();
 208   fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
 209   for (Int_t j=0; j<nTrackVariablesPair; j++) {
 210     fHistPP->SetBinLimits(j, dBinsPair[j]);
 211     fHistPP->SetVarTitle(j, axisTitlePair[j]);
 212   }
 213
 214   //-- pairs
 215   histName = "fHistNN";
 216   if(bShuffle) histName.Append("_shuffle");
 217   if(fCentralityId) histName += fCentralityId.Data();
 218   fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, nTrackVariablesPair, iBinPair);
 219   for (Int_t j=0; j<nTrackVariablesPair; j++) {
 220     fHistNN->SetBinLimits(j, dBinsPair[j]);
 221     fHistNN->SetVarTitle(j, axisTitlePair[j]);
 222   }
 223   AliInfo("Finished setting up the AliTHn");
 224 }
 225
 226 //____________________________________________________________________//
 227 void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
 228                                      TObjArray *particles,
 229                                      TObjArray *particlesMixed ) {
 230   // Calculates the balance function
 231   fAnalyzedEvents++;
 232
 233   // Initialize histograms if not done yet
 234   if(!fHistPN){
 235     AliWarning("Histograms not yet initialized! --> Will be done now");
 236     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 237     InitHistograms();
 238   }
 239
 240   Double_t trackVariablesSingle[nTrackVariablesSingle];
 241   Double_t trackVariablesPair[nTrackVariablesPair];
 242
 243   if (!particles){
 244     AliWarning("particles TObjArray is NULL pointer --> return");
 245     return;
 246   }
 247
 248   // define end of particle loops
 249   Int_t iMax = particles->GetEntriesFast();
 250   Int_t jMax = iMax;
 251   if (particlesMixed)
 252     jMax = particlesMixed->GetEntriesFast();
 253
 254   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
 255   TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
 256
 257   TArrayF secondEta(jMax);
 258   TArrayF secondPhi(jMax);
 259   TArrayF secondPt(jMax);
 260
 261   for (Int_t i=0; i<jMax; i++){
 262     secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
 263     secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
 264     secondPt[i]  = ((AliVParticle*) particlesSecond->At(i))->Pt();
 265   }
 266
 267   // 1st particle loop
 268   for (Int_t i=0; i<iMax; i++)
 269     {
 270
 271       AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
 272
 273       // some optimization
 274       Float_t firstEta = firstParticle->Eta();
 275       Float_t firstPhi = firstParticle->Phi();
 276       Float_t firstPt  = firstParticle->Pt();
 277
 278       // Event plane (determine psi bin)
 279       Double_t gPsiMinusPhi    =   0.;
 280       Double_t gPsiMinusPhiBin = -10.;
 281       gPsiMinusPhi   = TMath::Abs(firstPhi - gReactionPlane);
 282       if((gPsiMinusPhi <= 7.5)||
 283          ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
 284         gPsiMinusPhiBin = 0.0;
 285       else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
 286               ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
 287               ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
 288               ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
 289         gPsiMinusPhiBin = 1.0;
 290       else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
 291               ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
 292         gPsiMinusPhiBin = 2.0;
 293       else continue;
 294
 295
 296       Short_t  charge = (Short_t) firstParticle->Charge();
 297
 298       trackVariablesSingle[0]    =  gPsiMinusPhiBin;
 299       trackVariablesSingle[1]    =  firstPt;
 300
 301       //fill single particle histograms
 302       if(charge > 0)      fHistP->Fill(trackVariablesSingle,0,1.);
 303       else if(charge < 0) fHistN->Fill(trackVariablesSingle,0,1.);
 304
 305
 306
 307       // 2nd particle loop (only for j < i for non double counting in the same pT region)
 308       // --> SAME pT region for trigger and assoc: NO double counting with this
 309       // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i
 310       //                          --> can be handled afterwards by using assoc. as trigger as well ?!
 311       for(Int_t j = 0; j < i; j++) {   // or go to full here (everything prepared)?
 312
 313         if (particlesMixed && jMax < i)  // if the mixed track number is smaller than the main event one (could be done better if one loops over all tracks)
 314           break;
 315
 316         AliVParticle* secondParticle = (AliVParticle*) particlesSecond->At(j);
 317
 318         Short_t charge2 = (Short_t) secondParticle->Charge();
 319
 320         trackVariablesPair[0]    =  gPsiMinusPhiBin;
 321         trackVariablesPair[1]    =  firstEta - secondEta[j];  // delta eta
 322         trackVariablesPair[2]    =  firstPhi - secondPhi[j];  // delta phi
 323         if (trackVariablesPair[2] > 180.)   // delta phi between -180 and 180
 324           trackVariablesPair[2] -= 360.;
 325         if (trackVariablesPair[2] <  - 180.)
 326           trackVariablesPair[2] += 360.;
 327
 328         trackVariablesPair[3]    =  firstPt;      // pt trigger
 329         trackVariablesPair[4]    =  secondPt[j];  // pt
 330         //      trackVariablesPair[5]    =  fCentrality;  // centrality
 331
 332         if( charge > 0 && charge2 < 0)  fHistPN->Fill(trackVariablesPair,0,1.);
 333         else if( charge < 0 && charge2 > 0)  fHistNP->Fill(trackVariablesPair,0,1.);
 334         else if( charge > 0 && charge2 > 0)  fHistPP->Fill(trackVariablesPair,0,1.);
 335         else if( charge < 0 && charge2 < 0)  fHistNN->Fill(trackVariablesPair,0,1.);
 336         else AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
 337
 338       }//end of 2nd particle loop
 339     }//end of 1st particle loop
 340 }
 341
 342 //____________________________________________________________________//
 343 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
 344                                                  Int_t iVariablePair,
 345                                                  Double_t psiMin,
 346                                                  Double_t psiMax) {
 347   //Returns the BF histogram, extracted from the 6 AliTHn objects
 348   //(private members) of the AliBalancePsi class.
 349   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 350   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 351   TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","phi","qlong","qout","qside","qinv"};
 352   TString histName = "gHistBalanceFunctionHistogram";
 353   histName += gAnalysisType[iVariablePair];
 354
 355   // Psi_2
 356   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 357   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 358   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 359   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 360   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 361   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 362
 363   //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));
 364
 365   // Project into the wanted space (1st: analysis step, 2nd: axis)
 366   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
 367   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
 368   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
 369   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
 370   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
 371   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
 372
 373   TH1D *gHistBalanceFunctionHistogram = 0x0;
 374   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 375     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
 376     gHistBalanceFunctionHistogram->Reset();
 377
 378     switch(iVariablePair) {
 379     case 1:
 380       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
 381       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
 382       break;
 383     case 2:
 384       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
 385       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
 386       break;
 387     default:
 388       break;
 389     }
 390
 391     hTemp1->Sumw2();
 392     hTemp2->Sumw2();
 393     hTemp3->Sumw2();
 394     hTemp4->Sumw2();
 395     hTemp1->Add(hTemp3,-1.);
 396     hTemp1->Scale(1./hTemp5->GetEntries());
 397     hTemp2->Add(hTemp4,-1.);
 398     hTemp2->Scale(1./hTemp6->GetEntries());
 399     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 400     gHistBalanceFunctionHistogram->Scale(0.5);
 401   }
 402
 403   return gHistBalanceFunctionHistogram;
 404 }
 405
 406 //____________________________________________________________________//
 407 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
 408                                                         Double_t psiMax) {
 409   //Returns the BF histogram in Delta eta vs Delta phi,
 410   //extracted from the 6 AliTHn objects
 411   //(private members) of the AliBalancePsi class.
 412   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 413   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 414   TString histName = "gHistBalanceFunctionHistogram2D";
 415
 416   // Psi_2
 417   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 418   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 419   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 420   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 421   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 422   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 423
 424   AliInfo(Form("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)));
 425
 426   // Project into the wanted space (1st: analysis step, 2nd: axis)
 427   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
 428   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
 429   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
 430   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
 431   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
 432   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
 433
 434   TH2D *gHistBalanceFunctionHistogram = 0x0;
 435   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 436     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
 437     gHistBalanceFunctionHistogram->Reset();
 438     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
 439     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta #phi (deg.)");
 440     gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta #eta,#Delta #phi)");
 441
 442     hTemp1->Sumw2();
 443     hTemp2->Sumw2();
 444     hTemp3->Sumw2();
 445     hTemp4->Sumw2();
 446     hTemp1->Add(hTemp3,-1.);
 447     hTemp1->Scale(1./hTemp5->GetEntries());
 448     hTemp2->Add(hTemp4,-1.);
 449     hTemp2->Scale(1./hTemp6->GetEntries());
 450     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 451     gHistBalanceFunctionHistogram->Scale(0.5);
 452   }
 453
 454   return gHistBalanceFunctionHistogram;
 455 }
 456
 457 //____________________________________________________________________//
 458 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
 459                                               Double_t psiMax) {
 460   //Returns the 2D correlation function for (+-) pairs
 461   // Psi_2: axis 0
 462   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 463   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 464   //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
 465   //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
 466
 467   TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
 468   TCanvas *c1 = new TCanvas("c1","");
 469   c1->cd();
 470   if(!gHistTest){
 471     AliError("Projection of fHistP = NULL");
 472     return gHistTest;
 473   }
 474   else{
 475     gHistTest->DrawCopy("colz");
 476   }
 477
 478   //0:step, 2: Delta eta, 3: Delta phi
 479   TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
 480   if(!gHist){
 481     AliError("Projection of fHistPN = NULL");
 482     return gHist;
 483   }
 484
 485   AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
 486   AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
 487   AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
 488
 489   TCanvas *c2 = new TCanvas("c2","");
 490   c2->cd();
 491   fHistPN->Project(0,1,2)->DrawCopy("colz");
 492
 493   if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
 494     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
 495
 496   return gHist;
 497 }
 498
 499 //____________________________________________________________________//
 500 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
 501                                               Double_t psiMax) {
 502   //Returns the 2D correlation function for (+-) pairs
 503   // Psi_2: axis 0
 504   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 505   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 506
 507   //0:step, 2: Delta eta, 3: Delta phi
 508   TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
 509   if(!gHist){
 510     AliError("Projection of fHistPN = NULL");
 511     return gHist;
 512   }
 513
 514   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
 515   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
 516   if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
 517     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
 518
 519   return gHist;
 520 }
 521
 522 //____________________________________________________________________//
 523 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
 524                                               Double_t psiMax) {
 525   //Returns the 2D correlation function for (+-) pairs
 526   // Psi_2: axis 0
 527   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 528   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 529
 530   //0:step, 2: Delta eta, 3: Delta phi
 531   TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
 532   if(!gHist){
 533     AliError("Projection of fHistPN = NULL");
 534     return gHist;
 535   }
 536
 537   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
 538   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
 539   if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
 540     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
 541
 542   return gHist;
 543 }
 544
 545 //____________________________________________________________________//
 546 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
 547                                               Double_t psiMax) {
 548   //Returns the 2D correlation function for (+-) pairs
 549   // Psi_2: axis 0
 550   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 551   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 552
 553   //0:step, 2: Delta eta, 3: Delta phi
 554   TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
 555   if(!gHist){
 556     AliError("Projection of fHistPN = NULL");
 557     return gHist;
 558   }
 559
 560   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
 561   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
 562   if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
 563     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
 564
 565   return gHist;
 566 }
 567