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   fShuffle(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   fHistHBTbefore(0),
  61   fHistHBTafter(0),
  62   fHistConversionbefore(0),
  63   fHistConversionafter(0),
  64   fHistPshiMinusPhi(0),
  65   fPsiInterval(15.),
  66   fHBTCut(kFALSE),
  67   fConversionCut(kFALSE) {
  68   // Default constructor
  69 }
  70
  71 //____________________________________________________________________//
  72 AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
  73   TObject(balance), fShuffle(balance.fShuffle),
  74   fAnalysisLevel(balance.fAnalysisLevel),
  75   fAnalyzedEvents(balance.fAnalyzedEvents),
  76   fCentralityId(balance.fCentralityId),
  77   fCentStart(balance.fCentStart),
  78   fCentStop(balance.fCentStop),
  79   fHistP(balance.fHistP),
  80   fHistN(balance.fHistN),
  81   fHistPN(balance.fHistPN),
  82   fHistNP(balance.fHistNP),
  83   fHistPP(balance.fHistPP),
  84   fHistNN(balance.fHistNN),
  85   fHistHBTbefore(balance.fHistHBTbefore),
  86   fHistHBTafter(balance.fHistHBTafter),
  87   fHistConversionbefore(balance.fHistConversionbefore),
  88   fHistConversionafter(balance.fHistConversionafter),
  89   fHistPshiMinusPhi(balance.fHistPshiMinusPhi),
  90   fPsiInterval(balance.fPsiInterval),
  91   fHBTCut(balance.fHBTCut),
  92   fConversionCut(balance.fConversionCut) {
  93   //copy constructor
  94 }
  95
  96 //____________________________________________________________________//
  97 AliBalancePsi::~AliBalancePsi() {
  98   // Destructor
  99   delete fHistP;
 100   delete fHistN;
 101   delete fHistPN;
 102   delete fHistNP;
 103   delete fHistPP;
 104   delete fHistNN;
 105
 106   delete fHistHBTbefore;
 107   delete fHistHBTafter;
 108   delete fHistConversionbefore;
 109   delete fHistConversionafter;
 110   delete fHistPshiMinusPhi;
 111 }
 112
 113 //____________________________________________________________________//
 114 void AliBalancePsi::InitHistograms() {
 115   // single particle histograms
 116   Int_t anaSteps   = 1;       // analysis steps
 117   Int_t iBinSingle[kTrackVariablesSingle];        // binning for track variables
 118   Double_t* dBinsSingle[kTrackVariablesSingle];   // bins for track variables
 119   TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
 120
 121   // two particle histograms
 122   Int_t iBinPair[kTrackVariablesPair];         // binning for track variables
 123   Double_t* dBinsPair[kTrackVariablesPair];    // bins for track variables
 124   TString axisTitlePair[kTrackVariablesPair];  // axis titles for track variables
 125
 126   //centrality
 127   /*const Int_t kNCentralityBins = 9;
 128   Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
 129   iBinSingle[0]       = kNCentralityBins;
 130   dBinsSingle[0]      = centralityBins;
 131   axisTitleSingle[0]  = "Centrality percentile [%]";
 132   iBinPair[0]       = kNCentralityBins;
 133   dBinsPair[0]      = centralityBins;
 134   axisTitlePair[0]  = "Centrality percentile [%]"; */
 135
 136   //Psi_2: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (all)
 137   const Int_t kNPsi2Bins = 4;
 138   Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
 139   iBinSingle[0]       = kNPsi2Bins;
 140   dBinsSingle[0]      = psi2Bins;
 141   axisTitleSingle[0]  = "#phi - #Psi_{2} (a.u.)";
 142   iBinPair[0]       = kNPsi2Bins;
 143   dBinsPair[0]      = psi2Bins;
 144   axisTitlePair[0]  = "#phi - #Psi_{2} (a.u.)";
 145
 146    // delta eta
 147   const Int_t kNDeltaEtaBins = 80;
 148   Double_t deltaEtaBins[kNDeltaEtaBins+1];
 149   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
 150     deltaEtaBins[i] = -2.0 + i * 0.05;
 151   iBinPair[1]       = kNDeltaEtaBins;
 152   dBinsPair[1]      = deltaEtaBins;
 153   axisTitlePair[1]  = "#Delta #eta";
 154
 155    // delta phi
 156   const Int_t kNDeltaPhiBins = 72;
 157   Double_t deltaPhiBins[kNDeltaPhiBins+1];
 158   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
 159     deltaPhiBins[i] = -180.0 + i * 5.;
 160   }
 161   iBinPair[2]       = kNDeltaPhiBins;
 162   dBinsPair[2]      = deltaPhiBins;
 163   axisTitlePair[2]  = "#Delta #phi (#circ)";
 164
 165   // pt(trigger-associated)
 166   const Int_t kNPtBins = 16;
 167   Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
 168   //for(Int_t i = 0; i < kNPtBins+1; i++){
 169   //ptBins[i] = 0.2 + i * 0.5;
 170   //}
 171   iBinSingle[1]       = kNPtBins;
 172   dBinsSingle[1]      = ptBins;
 173   axisTitleSingle[1]  = "p_{t}^{trig.} (GeV/c)";
 174
 175   iBinPair[3]       = kNPtBins;
 176   dBinsPair[3]      = ptBins;
 177   axisTitlePair[3]  = "p_{t}^{trig.} (GeV/c)";
 178
 179   iBinPair[4]       = kNPtBins;
 180   dBinsPair[4]      = ptBins;
 181   axisTitlePair[4]  = "p_{t}^{assoc.} (GeV/c)";
 182
 183   TString histName;
 184   //+ triggered particles
 185   histName = "fHistP";
 186   if(fShuffle) histName.Append("_shuffle");
 187   if(fCentralityId) histName += fCentralityId.Data();
 188   fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
 189   for (Int_t j=0; j<kTrackVariablesSingle; j++) {
 190     fHistP->SetBinLimits(j, dBinsSingle[j]);
 191     fHistP->SetVarTitle(j, axisTitleSingle[j]);
 192   }
 193
 194   //- triggered particles
 195   histName = "fHistN";
 196   if(fShuffle) histName.Append("_shuffle");
 197   if(fCentralityId) histName += fCentralityId.Data();
 198   fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
 199   for (Int_t j=0; j<kTrackVariablesSingle; j++) {
 200     fHistN->SetBinLimits(j, dBinsSingle[j]);
 201     fHistN->SetVarTitle(j, axisTitleSingle[j]);
 202   }
 203
 204   //+- pairs
 205   histName = "fHistPN";
 206   if(fShuffle) histName.Append("_shuffle");
 207   if(fCentralityId) histName += fCentralityId.Data();
 208   fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 209   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 210     fHistPN->SetBinLimits(j, dBinsPair[j]);
 211     fHistPN->SetVarTitle(j, axisTitlePair[j]);
 212   }
 213
 214   //-+ pairs
 215   histName = "fHistNP";
 216   if(fShuffle) histName.Append("_shuffle");
 217   if(fCentralityId) histName += fCentralityId.Data();
 218   fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 219   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 220     fHistNP->SetBinLimits(j, dBinsPair[j]);
 221     fHistNP->SetVarTitle(j, axisTitlePair[j]);
 222   }
 223
 224   //++ pairs
 225   histName = "fHistPP";
 226   if(fShuffle) histName.Append("_shuffle");
 227   if(fCentralityId) histName += fCentralityId.Data();
 228   fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 229   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 230     fHistPP->SetBinLimits(j, dBinsPair[j]);
 231     fHistPP->SetVarTitle(j, axisTitlePair[j]);
 232   }
 233
 234   //-- pairs
 235   histName = "fHistNN";
 236   if(fShuffle) histName.Append("_shuffle");
 237   if(fCentralityId) histName += fCentralityId.Data();
 238   fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 239   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 240     fHistNN->SetBinLimits(j, dBinsPair[j]);
 241     fHistNN->SetVarTitle(j, axisTitlePair[j]);
 242   }
 243   AliInfo("Finished setting up the AliTHn");
 244
 245   // QA histograms
 246   fHistHBTbefore        = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,200);
 247   fHistHBTafter         = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,200);
 248   fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,200);
 249   fHistConversionafter  = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,200);
 250   fHistPshiMinusPhi     = new TH2D("fHistPshiMinusPhi","",4,-0.5,3.5,100,0,360.);
 251 }
 252
 253 //____________________________________________________________________//
 254 void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
 255                                      TObjArray *particles,
 256                                      TObjArray *particlesMixed,
 257                                      Float_t bSign) {
 258   // Calculates the balance function
 259   fAnalyzedEvents++;
 260
 261   // Initialize histograms if not done yet
 262   if(!fHistPN){
 263     AliWarning("Histograms not yet initialized! --> Will be done now");
 264     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 265     InitHistograms();
 266   }
 267
 268   Double_t trackVariablesSingle[kTrackVariablesSingle];
 269   Double_t trackVariablesPair[kTrackVariablesPair];
 270
 271   if (!particles){
 272     AliWarning("particles TObjArray is NULL pointer --> return");
 273     return;
 274   }
 275
 276   // define end of particle loops
 277   Int_t iMax = particles->GetEntriesFast();
 278   Int_t jMax = iMax;
 279   if (particlesMixed)
 280     jMax = particlesMixed->GetEntriesFast();
 281
 282   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
 283   TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
 284
 285   TArrayF secondEta(jMax);
 286   TArrayF secondPhi(jMax);
 287   TArrayF secondPt(jMax);
 288
 289   for (Int_t i=0; i<jMax; i++){
 290     secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
 291     secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
 292     secondPt[i]  = ((AliVParticle*) particlesSecond->At(i))->Pt();
 293   }
 294
 295   // 1st particle loop
 296   for (Int_t i = 0; i < iMax; i++) {
 297     AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
 298
 299     // some optimization
 300     Float_t firstEta = firstParticle->Eta();
 301     Float_t firstPhi = firstParticle->Phi();
 302     Float_t firstPt  = firstParticle->Pt();
 303
 304     // Event plane (determine psi bin)
 305     Double_t gPsiMinusPhi    =   0.;
 306     Double_t gPsiMinusPhiBin = -10.;
 307     gPsiMinusPhi   = TMath::Abs(firstPhi - gReactionPlane);
 308     //in-plane
 309     if((gPsiMinusPhi <= 7.5)||
 310        ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
 311       gPsiMinusPhiBin = 0.0;
 312     //intermediate
 313     else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
 314             ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
 315             ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
 316             ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
 317       gPsiMinusPhiBin = 1.0;
 318     //out of plane
 319     else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
 320             ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
 321       gPsiMinusPhiBin = 2.0;
 322     //everything else
 323     else
 324       gPsiMinusPhiBin = 3.0;
 325
 326     fHistPshiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
 327
 328     Short_t  charge1 = (Short_t) firstParticle->Charge();
 329
 330     trackVariablesSingle[0]    =  gPsiMinusPhiBin;
 331     trackVariablesSingle[1]    =  firstPt;
 332
 333     //fill single particle histograms
 334     if(charge1 > 0)      fHistP->Fill(trackVariablesSingle,0,1.);
 335     else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,1.);
 336
 337     // 2nd particle loop (only for j < i for non double counting in the same pT region)
 338     // --> SAME pT region for trigger and assoc: NO double counting with this
 339     // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i
 340     //                          --> can be handled afterwards by using assoc. as trigger as well ?!
 341     for(Int_t j = 0; j < iMax; j++) {
 342
 343       if (particlesMixed && j == jMax-1 )  // if the mixed track number is smaller than the main event one
 344         break;
 345
 346       if(j == i) continue; // no auto correlations
 347
 348       AliVParticle* secondParticle = (AliVParticle*) particlesSecond->At(j);
 349
 350       Short_t charge2 = (Short_t) secondParticle->Charge();
 351
 352       trackVariablesPair[0]    =  gPsiMinusPhiBin;
 353       trackVariablesPair[1]    =  firstEta - secondEta[j];  // delta eta
 354       trackVariablesPair[2]    =  firstPhi - secondPhi[j];  // delta phi
 355       if (trackVariablesPair[2] > 180.)   // delta phi between -180 and 180
 356         trackVariablesPair[2] -= 360.;
 357       if (trackVariablesPair[2] <  - 180.)
 358         trackVariablesPair[2] += 360.;
 359
 360       trackVariablesPair[3]    =  firstPt;      // pt trigger
 361       trackVariablesPair[4]    =  secondPt[j];  // pt
 362       //        trackVariablesPair[5]    =  fCentrality;  // centrality
 363
 364       // HBT like cut
 365       if(fHBTCut && charge1 * charge2 > 0){
 366         //if( dphi < 3 || deta < 0.01 ){   // VERSION 1
 367         //  continue;
 368
 369         Double_t deta = firstEta - secondEta[j];
 370         Double_t dphi = firstPhi - secondPhi[j];
 371         // VERSION 2 (Taken from DPhiCorrelations)
 372         // the variables & cuthave been developed by the HBT group
 373         // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
 374         fHistHBTbefore->Fill(deta,dphi);
 375
 376         // optimization
 377         if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
 378           {
 379             // phi in rad
 380             Float_t phi1rad = firstPhi*TMath::DegToRad();
 381             Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
 382
 383             // check first boundaries to see if is worth to loop and find the minimum
 384             Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
 385             Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
 386
 387             const Float_t kLimit = 0.02 * 3;
 388
 389             Float_t dphistarminabs = 1e5;
 390             Float_t dphistarmin = 1e5;
 391
 392             if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
 393               for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
 394                 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
 395                 Float_t dphistarabs = TMath::Abs(dphistar);
 396
 397                 if (dphistarabs < dphistarminabs) {
 398                   dphistarmin = dphistar;
 399                   dphistarminabs = dphistarabs;
 400                 }
 401               }
 402
 403               if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
 404                 //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));
 405                 continue;
 406               }
 407             }
 408           }
 409         fHistHBTafter->Fill(deta,dphi);
 410       }//HBT cut
 411
 412       // conversions
 413       if(fConversionCut) {
 414         if (charge1 * charge2 < 0) {
 415           Double_t deta = firstEta - secondEta[j];
 416           Double_t dphi = firstPhi - secondPhi[j];
 417           fHistConversionbefore->Fill(deta,dphi);
 418
 419           Float_t m0 = 0.510e-3;
 420           Float_t tantheta1 = 1e10;
 421
 422           // phi in rad
 423           Float_t phi1rad = firstPhi*TMath::DegToRad();
 424           Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
 425
 426           if (firstEta < -1e-10 || firstEta > 1e-10)
 427             tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
 428
 429           Float_t tantheta2 = 1e10;
 430           if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
 431             tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
 432
 433           Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
 434           Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
 435
 436           Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
 437
 438           if (masssqu < 0.04*0.04){
 439             //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));
 440             continue;
 441           }
 442           fHistConversionafter->Fill(deta,dphi);
 443         }
 444       }//conversion cut
 445
 446       if( charge1 > 0 && charge2 < 0)  fHistPN->Fill(trackVariablesPair,0,1.);
 447       else if( charge1 < 0 && charge2 > 0)  fHistNP->Fill(trackVariablesPair,0,1.);
 448       else if( charge1 > 0 && charge2 > 0)  fHistPP->Fill(trackVariablesPair,0,1.);
 449       else if( charge1 < 0 && charge2 < 0)  fHistNN->Fill(trackVariablesPair,0,1.);
 450       else {
 451         //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
 452         continue;
 453       }
 454     }//end of 2nd particle loop
 455   }//end of 1st particle loop
 456 }
 457
 458 //____________________________________________________________________//
 459 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
 460                                                  Int_t iVariablePair,
 461                                                  Double_t psiMin,
 462                                                  Double_t psiMax,
 463                                                  Double_t ptTriggerMin,
 464                                                  Double_t ptTriggerMax,
 465                                                  Double_t ptAssociatedMin,
 466                                                  Double_t ptAssociatedMax) {
 467   //Returns the BF histogram, extracted from the 6 AliTHn objects
 468   //(private members) of the AliBalancePsi class.
 469   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 470   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 471   // Psi_2
 472   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 473   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 474   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 475   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 476   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 477   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 478
 479   // pt trigger
 480   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 481     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 482     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 483     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 484     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 485     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 486     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 487   }
 488
 489   // pt associated
 490   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 491     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 492     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 493     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 494     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 495   }
 496
 497   //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));
 498
 499   // Project into the wanted space (1st: analysis step, 2nd: axis)
 500   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
 501   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
 502   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
 503   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
 504   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
 505   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
 506
 507   TH1D *gHistBalanceFunctionHistogram = 0x0;
 508   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 509     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
 510     gHistBalanceFunctionHistogram->Reset();
 511
 512     switch(iVariablePair) {
 513     case 1:
 514       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
 515       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
 516       break;
 517     case 2:
 518       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
 519       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
 520       break;
 521     default:
 522       break;
 523     }
 524
 525     hTemp1->Sumw2();
 526     hTemp2->Sumw2();
 527     hTemp3->Sumw2();
 528     hTemp4->Sumw2();
 529     hTemp1->Add(hTemp3,-1.);
 530     hTemp1->Scale(1./hTemp5->GetEntries());
 531     hTemp2->Add(hTemp4,-1.);
 532     hTemp2->Scale(1./hTemp6->GetEntries());
 533     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 534     gHistBalanceFunctionHistogram->Scale(0.5);
 535   }
 536
 537   return gHistBalanceFunctionHistogram;
 538 }
 539
 540 //____________________________________________________________________//
 541 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
 542                                                         Double_t psiMax,
 543                                                         Double_t ptTriggerMin,
 544                                                         Double_t ptTriggerMax,
 545                                                         Double_t ptAssociatedMin,
 546                                                         Double_t ptAssociatedMax) {
 547   //Returns the BF histogram in Delta eta vs Delta phi,
 548   //extracted from the 6 AliTHn objects
 549   //(private members) of the AliBalancePsi class.
 550   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 551   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 552   TString histName = "gHistBalanceFunctionHistogram2D";
 553
 554   // Psi_2
 555   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 556   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 557   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 558   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 559   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 560   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 561
 562   // pt trigger
 563   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 564     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 565     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 566     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 567     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 568     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 569     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 570   }
 571
 572   // pt associated
 573   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 574     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 575     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 576     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 577     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 578   }
 579
 580   //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)));
 581
 582   // Project into the wanted space (1st: analysis step, 2nd: axis)
 583   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
 584   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
 585   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
 586   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
 587   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
 588   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
 589
 590   TH2D *gHistBalanceFunctionHistogram = 0x0;
 591   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 592     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
 593     gHistBalanceFunctionHistogram->Reset();
 594     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
 595     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta #phi (deg.)");
 596     gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta #eta,#Delta #phi)");
 597
 598     hTemp1->Sumw2();
 599     hTemp2->Sumw2();
 600     hTemp3->Sumw2();
 601     hTemp4->Sumw2();
 602     hTemp1->Add(hTemp3,-1.);
 603     hTemp1->Scale(1./hTemp5->GetEntries());
 604     hTemp2->Add(hTemp4,-1.);
 605     hTemp2->Scale(1./hTemp6->GetEntries());
 606     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 607     gHistBalanceFunctionHistogram->Scale(0.5);
 608   }
 609
 610   return gHistBalanceFunctionHistogram;
 611 }
 612
 613 //____________________________________________________________________//
 614 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
 615                                               Double_t psiMax,
 616                                               Double_t ptTriggerMin,
 617                                               Double_t ptTriggerMax,
 618                                               Double_t ptAssociatedMin,
 619                                               Double_t ptAssociatedMax) {
 620   //Returns the 2D correlation function for (+-) pairs
 621   // Psi_2: axis 0
 622   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 623   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 624
 625   // pt trigger
 626   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 627     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 628     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 629   }
 630
 631   // pt associated
 632   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
 633     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 634
 635   //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
 636   //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
 637
 638   //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
 639   //TCanvas *c1 = new TCanvas("c1","");
 640   //c1->cd();
 641   //if(!gHistTest){
 642   //AliError("Projection of fHistP = NULL");
 643   //return gHistTest;
 644   //}
 645   //else{
 646   //gHistTest->DrawCopy("colz");
 647   //}
 648
 649   //0:step, 1: Delta eta, 2: Delta phi
 650   TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
 651   if(!gHist){
 652     AliError("Projection of fHistPN = NULL");
 653     return gHist;
 654   }
 655
 656   //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
 657   //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
 658   //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
 659
 660   //TCanvas *c2 = new TCanvas("c2","");
 661   //c2->cd();
 662   //fHistPN->Project(0,1,2)->DrawCopy("colz");
 663
 664   if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
 665     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
 666
 667   return gHist;
 668 }
 669
 670 //____________________________________________________________________//
 671 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
 672                                               Double_t psiMax,
 673                                               Double_t ptTriggerMin,
 674                                               Double_t ptTriggerMax,
 675                                               Double_t ptAssociatedMin,
 676                                               Double_t ptAssociatedMax) {
 677   //Returns the 2D correlation function for (+-) pairs
 678   // Psi_2: axis 0
 679   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 680   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 681
 682   // pt trigger
 683   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 684     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 685     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 686   }
 687
 688   // pt associated
 689   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
 690     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 691
 692   //0:step, 1: Delta eta, 2: Delta phi
 693   TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
 694   if(!gHist){
 695     AliError("Projection of fHistPN = NULL");
 696     return gHist;
 697   }
 698
 699   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
 700   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
 701   if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
 702     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
 703
 704   return gHist;
 705 }
 706
 707 //____________________________________________________________________//
 708 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
 709                                               Double_t psiMax,
 710                                               Double_t ptTriggerMin,
 711                                               Double_t ptTriggerMax,
 712                                               Double_t ptAssociatedMin,
 713                                               Double_t ptAssociatedMax) {
 714   //Returns the 2D correlation function for (+-) pairs
 715   // Psi_2: axis 0
 716   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 717   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 718
 719   // pt trigger
 720   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 721     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 722     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 723   }
 724
 725   // pt associated
 726   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
 727     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 728
 729   //0:step, 1: Delta eta, 2: Delta phi
 730   TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
 731   if(!gHist){
 732     AliError("Projection of fHistPN = NULL");
 733     return gHist;
 734   }
 735
 736   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
 737   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
 738   if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
 739     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
 740
 741   return gHist;
 742 }
 743
 744 //____________________________________________________________________//
 745 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
 746                                               Double_t psiMax,
 747                                               Double_t ptTriggerMin,
 748                                               Double_t ptTriggerMax,
 749                                               Double_t ptAssociatedMin,
 750                                               Double_t ptAssociatedMax) {
 751   //Returns the 2D correlation function for (+-) pairs
 752   // Psi_2: axis 0
 753   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 754   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 755
 756   // pt trigger
 757   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 758     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 759     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 760   }
 761
 762   // pt associated
 763   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
 764     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 765
 766   //0:step, 1: Delta eta, 2: Delta phi
 767   TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
 768   if(!gHist){
 769     AliError("Projection of fHistPN = NULL");
 770     return gHist;
 771   }
 772
 773   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
 774   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
 775   if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
 776     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
 777
 778   return gHist;
 779 }
 780
 781 //____________________________________________________________________//
 782 Float_t AliBalancePsi::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign) {
 783   //
 784   // calculates dphistar
 785   //
 786   Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
 787
 788   static const Double_t kPi = TMath::Pi();
 789
 790   // circularity
 791 //   if (dphistar > 2 * kPi)
 792 //     dphistar -= 2 * kPi;
 793 //   if (dphistar < -2 * kPi)
 794 //     dphistar += 2 * kPi;
 795
 796   if (dphistar > kPi)
 797     dphistar = kPi * 2 - dphistar;
 798   if (dphistar < -kPi)
 799     dphistar = -kPi * 2 - dphistar;
 800   if (dphistar > kPi) // might look funny but is needed
 801     dphistar = kPi * 2 - dphistar;
 802
 803   return dphistar;
 804 }
 805