PWGCF/EBYE/BalanceFunctions/AliBalanceTriggered.cxx

   1 /**************************************************************************
   2  * Author: Panos Christakoglou.                                           *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 //-----------------------------------------------------------------
  15 //           Balance Function class
  16 //   This is the class to deal with the Balance Function analysis
  17 //   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
  18 //   Modified: Michael Weber, m.weber@cern.ch
  19 //-----------------------------------------------------------------
  20
  21
  22 //ROOT
  23 #include <Riostream.h>
  24 #include <TMath.h>
  25 #include <TAxis.h>
  26 #include <TH1D.h>
  27 #include <TH2D.h>
  28
  29
  30 #include <AliTHn.h>
  31 #include <TLorentzVector.h>
  32 #include <TObjArray.h>
  33 #include <TGraphErrors.h>
  34 #include <TString.h>
  35
  36 #include "AliVParticle.h"
  37 #include "AliMCParticle.h"
  38 #include "AliESDtrack.h"
  39 #include "AliAODTrack.h"
  40
  41 #include "AliBalanceTriggered.h"
  42
  43 ClassImp(AliBalanceTriggered)
  44
  45 //____________________________________________________________________//
  46 AliBalanceTriggered::AliBalanceTriggered() :
  47   TObject(),
  48   fAnalysisLevel("AOD"),
  49   fHistP(0x0),
  50   fHistN(0x0),
  51   fHistPN(0x0),
  52   fHistNP(0x0),
  53   fHistPP(0x0),
  54   fHistNN(0x0)
  55 {
  56   // Default constructor
  57
  58 }
  59
  60
  61 //____________________________________________________________________//
  62 AliBalanceTriggered::AliBalanceTriggered(const AliBalanceTriggered& balance):
  63   TObject(balance),
  64   fAnalysisLevel(balance.fAnalysisLevel),
  65   fHistP(balance.fHistP),
  66   fHistN(balance.fHistN),
  67   fHistPN(balance.fHistPN),
  68   fHistNP(balance.fHistNP),
  69   fHistPP(balance.fHistPP),
  70   fHistNN(balance.fHistNN)
  71 {
  72   //copy constructor
  73
  74 }
  75
  76 //____________________________________________________________________//
  77 AliBalanceTriggered::~AliBalanceTriggered() {
  78   // Destructor
  79
  80   delete fHistP;
  81   delete fHistN;
  82   delete fHistPN;
  83   delete fHistNP;
  84   delete fHistPP;
  85   delete fHistNN;
  86
  87
  88 }
  89
  90 //____________________________________________________________________//
  91 void AliBalanceTriggered::InitHistograms() {
  92   //Initialize the histograms
  93
  94   // global switch disabling the reference
  95   // (to avoid "Replacing existing TH1" if several wagons are created in train)
  96   Bool_t oldStatus = TH1::AddDirectoryStatus();
  97   TH1::AddDirectory(kFALSE);
  98
  99   TString title    = "";      // histogram title
 100   Int_t anaSteps   = 1;       // analysis steps
 101
 102   // single particle histograms
 103   Int_t iBinSingle[kTrackVarsSingle];         // binning for track variables
 104   Double_t* dBinsSingle[kTrackVarsSingle];    // bins for track variables
 105   TString axisTitleSingle[kTrackVarsSingle];  // axis titles for track variables
 106
 107   // two particle histograms
 108   Int_t iBinPair[kTrackVarsPair];         // binning for track variables
 109   Double_t* dBinsPair[kTrackVarsPair];    // bins for track variables
 110   TString axisTitlePair[kTrackVarsPair];  // axis titles for track variables
 111
 112
 113   //-----------------------------------------------------------
 114   // histogram settings (hard coded!)
 115   //-----------------------------------------------------------
 116
 117   // eta
 118   const Int_t kNEtaBins = 40;
 119   Double_t etaBins[kNEtaBins+1];
 120   for(Int_t i = 0; i < kNEtaBins+1; i++){
 121     etaBins[i] = -1.0 + i * 0.05;
 122   }
 123   iBinSingle[0]       = kNEtaBins;
 124   dBinsSingle[0]      = etaBins;
 125   axisTitleSingle[0]  = "#eta";
 126
 127   // delta eta
 128   const Int_t kNDeltaEtaBins = 40;
 129   Double_t deltaEtaBins[kNDeltaEtaBins+1];
 130   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++){
 131     deltaEtaBins[i] = -2.0 + i * 0.1;
 132   }
 133   iBinPair[0]       = kNDeltaEtaBins;
 134   dBinsPair[0]      = deltaEtaBins;
 135   axisTitlePair[0]  = "#Delta #eta";
 136
 137   // phi
 138   const Int_t kNPhiBins = 72;
 139   Double_t phiBins[kNPhiBins+1];
 140   for(Int_t i = 0; i < kNPhiBins+1; i++){
 141     phiBins[i] = 0.0 + i * 5.;
 142   }
 143   iBinSingle[1]       = kNPhiBins;
 144   dBinsSingle[1]      = phiBins;
 145   axisTitleSingle[1]  = "#phi (#circ)";
 146
 147   // delta phi
 148   const Int_t kNDeltaPhiBins = 72;
 149   Double_t deltaPhiBins[kNDeltaPhiBins+1];
 150   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
 151     deltaPhiBins[i] = -180.0 + i * 5.;
 152   }
 153   iBinPair[1]       = kNDeltaPhiBins;
 154   dBinsPair[1]      = deltaPhiBins;
 155   axisTitlePair[1]  = "#Delta #phi (#circ)";
 156
 157   // pT
 158   const Int_t kNPtBins = 22;
 159   Double_t pTBins[kNPtBins+1] = {0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 15.0, 20.0};
 160   iBinSingle[2]       = kNPtBins;
 161   dBinsSingle[2]      = pTBins;
 162   axisTitleSingle[2]  = "p_{T,trig} (GeV/c)";
 163
 164   iBinPair[2]       = kNPtBins;
 165   dBinsPair[2]      = pTBins;
 166   axisTitlePair[2]  = "p_{T} (GeV/c)";
 167
 168   iBinPair[3]       = kNPtBins;
 169   dBinsPair[3]      = pTBins;
 170   axisTitlePair[3]  = "p_{T,trig} (GeV/c)";
 171
 172   // centrality
 173   const Int_t kNCentBins = 9;
 174   Double_t centBins[kNCentBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
 175
 176   iBinSingle[3]       = kNCentBins;
 177   dBinsSingle[3]      = centBins;
 178   axisTitleSingle[3]  = "centrality";
 179
 180   iBinPair[4]       = kNCentBins;
 181   dBinsPair[4]      = centBins;
 182   axisTitlePair[4]  = "centrality";
 183
 184   //-----------------------------------------------------------
 185
 186
 187
 188   // User settings (depending on the analysis settings [only one at the moment])
 189   if(fAnalysisLevel == "AOD"){
 190     title = "hdEtaVsdPhi";
 191   }
 192
 193   //-----------------------------------------------------------
 194   //-----------------------------------------------------------
 195   // Histogram creation
 196
 197   // histogram for negative particles
 198   fHistN = new AliTHn(Form("fHistN"), Form("%s_N",title.Data()), anaSteps, kTrackVarsSingle, iBinSingle);
 199   for (Int_t j=0; j<kTrackVarsSingle; j++)
 200     {
 201       fHistN->SetBinLimits(j, dBinsSingle[j]);
 202       fHistN->SetVarTitle(j, axisTitleSingle[j]);
 203     }
 204
 205   // histogram for positive particles
 206   fHistP = new AliTHn(Form("fHistP"), Form("%s_P",title.Data()), anaSteps, kTrackVarsSingle, iBinSingle);
 207   for (Int_t j=0; j<kTrackVarsSingle; j++)
 208     {
 209       fHistP->SetBinLimits(j, dBinsSingle[j]);
 210       fHistP->SetVarTitle(j, axisTitleSingle[j]);
 211     }
 212
 213   // histogram for +- pairs
 214   fHistPN = new AliTHn(Form("fHistPN"), Form("%s_PN",title.Data()), anaSteps, kTrackVarsPair, iBinPair);
 215   for (Int_t j=0; j<kTrackVarsPair; j++)
 216     {
 217       fHistPN->SetBinLimits(j, dBinsPair[j]);
 218       fHistPN->SetVarTitle(j, axisTitlePair[j]);
 219     }
 220
 221   // histogram for -+ pairs
 222   fHistNP = new AliTHn(Form("fHistNP"), Form("%s_NP",title.Data()), anaSteps, kTrackVarsPair, iBinPair);
 223   for (Int_t j=0; j<kTrackVarsPair; j++)
 224     {
 225       fHistNP->SetBinLimits(j, dBinsPair[j]);
 226       fHistNP->SetVarTitle(j, axisTitlePair[j]);
 227     }
 228
 229   // histogram for ++ pairs
 230   fHistPP = new AliTHn(Form("fHistPP"), Form("%s_PP",title.Data()), anaSteps, kTrackVarsPair, iBinPair);
 231   for (Int_t j=0; j<kTrackVarsPair; j++)
 232     {
 233       fHistPP->SetBinLimits(j, dBinsPair[j]);
 234       fHistPP->SetVarTitle(j, axisTitlePair[j]);
 235     }
 236
 237   // histogram for -- pairs
 238   fHistNN = new AliTHn(Form("fHistNN"), Form("%s_NN",title.Data()), anaSteps, kTrackVarsPair, iBinPair);
 239   for (Int_t j=0; j<kTrackVarsPair; j++)
 240     {
 241       fHistNN->SetBinLimits(j, dBinsPair[j]);
 242       fHistNN->SetVarTitle(j, axisTitlePair[j]);
 243     }
 244
 245   //-----------------------------------------------------------
 246   //-----------------------------------------------------------
 247
 248   TH1::AddDirectory(oldStatus);
 249
 250 }
 251
 252 //____________________________________________________________________//
 253 void AliBalanceTriggered::FillBalance(Float_t fCentrality,TObjArray *particles, TObjArray *particlesMixed) {
 254   // Calculates the balance function
 255
 256
 257   // Initialize histograms if not done yet
 258   if(!fHistPN){
 259     AliWarning("Histograms not yet initialized! --> Will be done now");
 260     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 261     InitHistograms();
 262   }
 263
 264   Double_t trackVarsSingle[kTrackVarsSingle];
 265   Double_t trackVarsPair[kTrackVarsPair];
 266
 267   if (!particles){
 268     AliWarning("particles TObjArray is NULL pointer --> return");
 269     return;
 270   }
 271
 272   // define end of particle loops
 273   Int_t iMax = particles->GetEntriesFast();
 274   Int_t jMax = iMax;
 275   if (particlesMixed)
 276     jMax = particlesMixed->GetEntriesFast();
 277
 278   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
 279   TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
 280
 281   TArrayF secondEta(jMax);
 282   TArrayF secondPhi(jMax);
 283   TArrayF secondPt(jMax);
 284
 285   for (Int_t i=0; i<jMax; i++){
 286     secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
 287     secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
 288     secondPt[i]  = ((AliVParticle*) particlesSecond->At(i))->Pt();
 289   }
 290
 291   // 1st particle loop
 292   for (Int_t i=0; i<iMax; i++)
 293     {
 294
 295       AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
 296
 297       // some optimization
 298       Float_t firstEta = firstParticle->Eta();
 299       Float_t firstPhi = firstParticle->Phi();
 300       Float_t firstPt  = firstParticle->Pt();
 301
 302
 303       Short_t  charge = (Short_t) firstParticle->Charge();
 304       trackVarsSingle[0]    =  firstEta;    //eta
 305       trackVarsSingle[1]    =  firstPhi;    //phi
 306       trackVarsSingle[2]    =  firstPt;     //pt trigger
 307       trackVarsSingle[3]    =  fCentrality; //centrality
 308
 309       //fill single particle histograms
 310       if(charge > 0)  fHistP->Fill(trackVarsSingle,0,1.);
 311       else            fHistN->Fill(trackVarsSingle,0,1.);
 312
 313
 314       // 2nd particle loop (only for j < i for non double counting in the same pT region)
 315       // --> SAME pT region for trigger and assoc: NO double counting with this
 316       // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i
 317       //                          --> can be handled afterwards by using assoc. as trigger as well ?!
 318       for(Int_t j = 0; j < i; j++) {   // or go to full here (everything prepared)?
 319
 320         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)
 321           break;
 322
 323         AliVParticle* secondParticle = (AliVParticle*) particlesSecond->At(j);
 324
 325         Short_t charge2 = (Short_t) secondParticle->Charge();
 326         trackVarsPair[0]    =  firstEta - secondEta[j];  // delta eta
 327         trackVarsPair[1]    =  firstPhi - secondPhi[j];  // delta phi
 328         if (trackVarsPair[1] > 180)   // delta phi between -180 and 180
 329           trackVarsPair[1] -= 360;
 330         if (trackVarsPair[1] <  - 180)
 331           trackVarsPair[1] += 360;
 332
 333         trackVarsPair[2]    =  firstPt;      // pt trigger
 334         trackVarsPair[3]    =  secondPt[j];  // pt
 335         trackVarsPair[4]    =  fCentrality;  // centrality
 336
 337         if( charge > 0 && charge2 < 0)  fHistPN->Fill(trackVarsPair,0,1.);
 338         else if( charge < 0 && charge2 > 0)  fHistNP->Fill(trackVarsPair,0,1.);
 339         else if( charge > 0 && charge2 > 0)  fHistPP->Fill(trackVarsPair,0,1.);
 340         else if( charge < 0 && charge2 < 0)  fHistNN->Fill(trackVarsPair,0,1.);
 341         else AliWarning("Wrong charge combination!");
 342
 343       }//end of 2nd particle loop
 344     }//end of 1st particle loop
 345 }
 346
 347
 348 TH1D* AliBalanceTriggered::GetBalanceFunctionHistogram1D(Int_t var, Double_t pTMinTrigger, Double_t pTMaxTrigger, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax, Double_t etaGap){
 349
 350   // check which variable should be analyzed
 351   // 0 = Delta eta
 352   // 1 = Delta phi
 353
 354   if( var < 0 || var > 1){
 355     AliError("Only Variable 0 (= Delta eta) or 1 (= Delta phi) allowed");
 356     return NULL;
 357   }
 358
 359
 360   // Choose region to analyze
 361   // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
 362   // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
 363
 364   // pT trigger
 365   fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 366   fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 367   fHistPN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 368   fHistNP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 369   fHistPP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 370   fHistNN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 371
 372   // pT
 373   fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 374   fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 375   fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 376   fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 377
 378   // centrality
 379   fHistP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax);
 380   fHistN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax);
 381   fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 382   fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 383   fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 384   fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 385
 386   if (etaGap > 0){
 387     fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(20 + etaGap, 9999);
 388     fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(20 + etaGap, 9999);
 389     fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(20 + etaGap, 9999);
 390     fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(20 + etaGap, 9999);
 391     fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(20 + etaGap, 9999);
 392     fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(20 + etaGap, 9999);
 393   }
 394
 395   // Project into the wanted space (1st: analysis step, 2nd: axis)
 396   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,var);
 397   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,var);
 398   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,var);
 399   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,var);
 400   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,var);
 401   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,var);
 402
 403   if (etaGap > 0){
 404     fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(-9999, 9999);
 405     fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(-9999, 9999);
 406     fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(-9999, 9999);
 407     fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(-9999, 9999);
 408     fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(-9999, 9999);
 409     fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(-9999, 9999);
 410   }
 411
 412   TH1D* gHistBalanceFunctionHistogram = NULL;
 413
 414   // Calculate BF
 415   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 416
 417     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
 418     gHistBalanceFunctionHistogram->Reset();
 419
 420     // Calculate BF
 421     hTemp1->Add(hTemp3,-1.);
 422     hTemp1->Scale(1./hTemp5->GetEntries());
 423     hTemp2->Add(hTemp4,-1.);
 424     hTemp2->Scale(1./hTemp6->GetEntries());
 425     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 426     gHistBalanceFunctionHistogram->Scale(0.5/1.);
 427   }
 428
 429   return gHistBalanceFunctionHistogram;
 430 }
 431
 432 TH2D* AliBalanceTriggered::GetBalanceFunctionHistogram2D(Int_t var1, Int_t var2, Double_t pTMinTrigger, Double_t pTMaxTrigger, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax){
 433
 434   // check which variable should be analyzed
 435   // 0 = Delta eta
 436   // 1 = Delta phi
 437
 438   if( var1 < 0 || var1 > 1 || var2 < 0 || var2 > 1){
 439     AliError("Only Variable 0 (= Delta eta) or 1 (= Delta phi) allowed");
 440     return NULL;
 441   }
 442
 443
 444   // Choose region to analyze
 445   // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
 446   // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
 447
 448   // pT trigger
 449   fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 450   fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 451   fHistPN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 452   fHistNP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 453   fHistPP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 454   fHistNN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 455
 456   // pT
 457   fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 458   fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 459   fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 460   fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 461
 462   // centrality
 463   fHistP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax);
 464   fHistN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax);
 465   fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 466   fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 467   fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 468   fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 469
 470
 471   // Project into the wanted space (1st: analysis step, 2nd: axis1, 3rd: axis2)
 472   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,var1,var2);
 473   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,var1,var2);
 474   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,var1,var2);
 475   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,var1,var2);
 476   TH2D* hTemp5 = (TH2D*)fHistP->Project(0,var1,var2);
 477   TH2D* hTemp6 = (TH2D*)fHistN->Project(0,var1,var2);
 478
 479   TH2D* gHistBalanceFunctionHistogram = NULL;
 480
 481   // Calculate BF
 482   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 483
 484     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
 485     gHistBalanceFunctionHistogram->Reset();
 486
 487     // Calculate BF
 488     hTemp1->Add(hTemp3,-1.);
 489     hTemp1->Scale(1./hTemp5->GetEntries());
 490     hTemp2->Add(hTemp4,-1.);
 491     hTemp2->Scale(1./hTemp6->GetEntries());
 492     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 493   }
 494
 495   return gHistBalanceFunctionHistogram;
 496 }
 497
 498
 499 TH1D* AliBalanceTriggered::GetHistogram1D(Int_t histo, Int_t var, Double_t pTMinTrigger, Double_t pTMaxTrigger, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax, Double_t etaGap){
 500
 501   // check which variable should be analyzed
 502   //
 503   // single histograms:
 504   // 0 = eta
 505   // 1 = phi
 506   // 2 = pT, trigger
 507   // 3 = centrality
 508   //
 509   // pair histograms:
 510   // 0 = Delta eta
 511   // 1 = Delta phi
 512   // 2 = pT
 513   // 3 = pT, trigger
 514   // 4 = centrality
 515
 516   if(histo < 0 || histo > 5){
 517     AliError("Only 6 histograms available: 0(P), 1(N), 2(PN), 3(NP), 4(PP), 5(NN)");
 518     return NULL;
 519   }
 520
 521   if( histo > 1 && (var < 0 || var > 5)){
 522     AliError("Only Variable 0 to 4 allowed for pair histograms (histo > 1)");
 523     return NULL;
 524   }
 525   if( histo < 2 && (var < 0 || var > 4)){
 526     AliError("Only Variable 0 to 3 allowed for single histograms (histo < 2)");
 527     return NULL;
 528   }
 529
 530   // get the histogram
 531   AliTHn *gTHn = NULL;
 532   switch(histo){
 533
 534   case 0:
 535     gTHn = fHistP;
 536     break;
 537
 538   case 1:
 539     gTHn = fHistN;
 540     break;
 541
 542   case 2:
 543     gTHn = fHistPN;
 544     break;
 545
 546   case 3:
 547     gTHn = fHistNP;
 548     break;
 549
 550   case 4:
 551     gTHn = fHistPP;
 552     break;
 553
 554   case 5:
 555     gTHn = fHistNN;
 556     break;
 557
 558   }
 559
 560   if(!gTHn){
 561     AliError(Form("AliTHn number %d = NULL",histo));
 562     return NULL;
 563   }
 564
 565   // Choose region to analyze
 566   // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
 567   // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
 568
 569   // pT trigger
 570   gTHn->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 571
 572   // pT
 573   if(histo > 1) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 574
 575   // centrality
 576   if(histo < 2) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax);
 577   else          gTHn->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 578
 579   if (etaGap > 0) gTHn->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(20 + etaGap, 9999);
 580
 581   // Project into the wanted space (1st: analysis step, 2nd: axis)
 582   TH1D* gHisto      = (TH1D*)gTHn->Project(0,var);
 583
 584
 585   gTHn->GetGrid(0)->GetGrid()->GetAxis(0)->SetRange(-9999, 9999);
 586
 587
 588   return gHisto;
 589 }
 590
 591
 592 TH2D* AliBalanceTriggered::GetHistogram2D(Int_t histo, Int_t var1, Int_t var2, Double_t pTMinTrigger, Double_t pTMaxTrigger, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax){
 593
 594   // check which variable should be analyzed
 595   //
 596   // single histograms:
 597   // 0 = eta
 598   // 1 = phi
 599   // 2 = pT, trigger
 600   // 3 = centrality
 601   //
 602   // pair histograms:
 603   // 0 = Delta eta
 604   // 1 = Delta phi
 605   // 2 = pT
 606   // 3 = pT, trigger
 607   // 4 = centrality
 608
 609   if(histo < 0 || histo > 5){
 610     AliError("Only 6 histograms available: 0(P), 1(N), 2(PN), 3(NP), 4(PP), 5(NN)");
 611     return NULL;
 612   }
 613
 614   if( histo > 1 && (var1 < 0 || var1 > 5 || var2 < 0 || var2 > 5)){
 615     AliError("Only Variable 0 to 4 allowed for pair histograms (histo > 1)");
 616     return NULL;
 617   }
 618   if( histo < 2 && (var1 < 0 || var1 > 4 || var2 < 0 || var2 > 4)){
 619     AliError("Only Variable 0 to 3 allowed for single histograms (histo < 2)");
 620     return NULL;
 621   }
 622
 623   // get the histogram
 624   AliTHn *gTHn = NULL;
 625   switch(histo){
 626
 627   case 0:
 628     gTHn = fHistP;
 629     break;
 630
 631   case 1:
 632     gTHn = fHistN;
 633     break;
 634
 635   case 2:
 636     gTHn = fHistPN;
 637     break;
 638
 639   case 3:
 640     gTHn = fHistNP;
 641     break;
 642
 643   case 4:
 644     gTHn = fHistPP;
 645     break;
 646
 647   case 5:
 648     gTHn = fHistNN;
 649     break;
 650
 651   }
 652
 653   if(!gTHn){
 654     AliError(Form("AliTHn number %d = NULL",histo));
 655     return NULL;
 656   }
 657
 658   // Choose region to analyze
 659   // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
 660   // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
 661
 662   // pT trigger
 663   gTHn->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger);
 664
 665   // pT
 666   if(histo > 1) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax);
 667
 668   // centrality
 669   if(histo < 2) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax);
 670   else          gTHn->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax);
 671
 672
 673   // Project into the wanted space (1st: analysis step, 2nd: axis)
 674   TH2D* gHisto = (TH2D*)gTHn->Project(0,var1,var2);
 675
 676   return gHisto;
 677 }