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 #include "AliAnalysisTaskTriggeredBF.h"
  41
  42 #include "AliBalancePsi.h"
  43
  44 ClassImp(AliBalancePsi)
  45
  46 //____________________________________________________________________//
  47 AliBalancePsi::AliBalancePsi() :
  48   TObject(),
  49   fShuffle(kFALSE),
  50   fAnalysisLevel("ESD"),
  51   fAnalyzedEvents(0) ,
  52   fCentralityId(0) ,
  53   fCentStart(0.),
  54   fCentStop(0.),
  55   fHistP(0),
  56   fHistN(0),
  57   fHistPN(0),
  58   fHistNP(0),
  59   fHistPP(0),
  60   fHistNN(0),
  61   fHistHBTbefore(0),
  62   fHistHBTafter(0),
  63   fHistConversionbefore(0),
  64   fHistConversionafter(0),
  65   fHistPsiMinusPhi(0),
  66   fPsiInterval(15.),
  67   fDeltaEtaMax(2.0),
  68   fHBTCut(kFALSE),
  69   fConversionCut(kFALSE),
  70   fEventClass("EventPlane"){
  71   // Default constructor
  72 }
  73
  74 //____________________________________________________________________//
  75 AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
  76   TObject(balance), fShuffle(balance.fShuffle),
  77   fAnalysisLevel(balance.fAnalysisLevel),
  78   fAnalyzedEvents(balance.fAnalyzedEvents),
  79   fCentralityId(balance.fCentralityId),
  80   fCentStart(balance.fCentStart),
  81   fCentStop(balance.fCentStop),
  82   fHistP(balance.fHistP),
  83   fHistN(balance.fHistN),
  84   fHistPN(balance.fHistPN),
  85   fHistNP(balance.fHistNP),
  86   fHistPP(balance.fHistPP),
  87   fHistNN(balance.fHistNN),
  88   fHistHBTbefore(balance.fHistHBTbefore),
  89   fHistHBTafter(balance.fHistHBTafter),
  90   fHistConversionbefore(balance.fHistConversionbefore),
  91   fHistConversionafter(balance.fHistConversionafter),
  92   fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
  93   fPsiInterval(balance.fPsiInterval),
  94   fDeltaEtaMax(balance.fDeltaEtaMax),
  95   fHBTCut(balance.fHBTCut),
  96   fConversionCut(balance.fConversionCut),
  97   fEventClass("EventPlane"){
  98   //copy constructor
  99 }
 100
 101 //____________________________________________________________________//
 102 AliBalancePsi::~AliBalancePsi() {
 103   // Destructor
 104   delete fHistP;
 105   delete fHistN;
 106   delete fHistPN;
 107   delete fHistNP;
 108   delete fHistPP;
 109   delete fHistNN;
 110
 111   delete fHistHBTbefore;
 112   delete fHistHBTafter;
 113   delete fHistConversionbefore;
 114   delete fHistConversionafter;
 115   delete fHistPsiMinusPhi;
 116
 117 }
 118
 119 //____________________________________________________________________//
 120 void AliBalancePsi::InitHistograms() {
 121   // single particle histograms
 122
 123   // global switch disabling the reference
 124   // (to avoid "Replacing existing TH1" if several wagons are created in train)
 125   Bool_t oldStatus = TH1::AddDirectoryStatus();
 126   TH1::AddDirectory(kFALSE);
 127
 128   Int_t anaSteps   = 1;       // analysis steps
 129   Int_t iBinSingle[kTrackVariablesSingle];        // binning for track variables
 130   Double_t* dBinsSingle[kTrackVariablesSingle];   // bins for track variables
 131   TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
 132
 133   // two particle histograms
 134   Int_t iBinPair[kTrackVariablesPair];         // binning for track variables
 135   Double_t* dBinsPair[kTrackVariablesPair];    // bins for track variables
 136   TString axisTitlePair[kTrackVariablesPair];  // axis titles for track variables
 137   /**********************************************************
 138
 139   ======> Modification: Change Event Classification Scheme
 140
 141   ---> fEventClass == "EventPlane"
 142
 143    Default operation with Event Plane
 144
 145   ---> fEventClass == "Multiplicity"
 146
 147    Work with kTPCITStracklet multiplicity (from GetReferenceMultiplicity)
 148
 149   ---> fEventClass == "Centrality"
 150
 151    Work with Centrality Bins
 152
 153   ***********************************************************/
 154
 155   //--- Multiplicity Bins ------------------------------------
 156     const Int_t kMultBins = 8;
 157     //A first rough attempt at four bins
 158     Double_t kMultBinLimits[kMultBins+1]={0,10,20,30,40,50,60,70,80};
 159   //----------------------------------------------------------
 160
 161   //--- Centrality Bins --------------------------------------
 162     const Int_t kNCentralityBins = 9;
 163     Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
 164   //----------------------------------------------------------
 165
 166   //--- Event Plane Bins -------------------------------------
 167     //Psi_2: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest)
 168     const Int_t kNPsi2Bins = 4;
 169     Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
 170   //----------------------------------------------------------
 171
 172   //Depending on fEventClass Variable, do one thing or the other...
 173     if(fEventClass == "Multiplicity"){
 174         iBinSingle[0]       = kMultBins;
 175         dBinsSingle[0]      = kMultBinLimits;
 176         axisTitleSingle[0]  = "kTPCITStracklet multiplicity";
 177         iBinPair[0]       = kMultBins;
 178         dBinsPair[0]      = kMultBinLimits;
 179         axisTitlePair[0]  = "kTPCITStracklet multiplicity";
 180     }
 181     if(fEventClass == "Centrality"){
 182         iBinSingle[0]       = kNCentralityBins;
 183         dBinsSingle[0]      = centralityBins;
 184         axisTitleSingle[0]  = "Centrality percentile [%]";
 185         iBinPair[0]       = kNCentralityBins;
 186         dBinsPair[0]      = centralityBins;
 187         axisTitlePair[0]  = "Centrality percentile [%]";
 188     }
 189     if(fEventClass == "EventPlane"){
 190         iBinSingle[0]       = kNPsi2Bins;
 191         dBinsSingle[0]      = psi2Bins;
 192         axisTitleSingle[0]  = "#varphi - #Psi_{2} (a.u.)";
 193         iBinPair[0]       = kNPsi2Bins;
 194         dBinsPair[0]      = psi2Bins;
 195         axisTitlePair[0]  = "#varphi - #Psi_{2} (a.u.)";
 196     }
 197
 198    // delta eta
 199   const Int_t kNDeltaEtaBins = 80;
 200   Double_t deltaEtaBins[kNDeltaEtaBins+1];
 201   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
 202     deltaEtaBins[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBins;
 203   iBinPair[1]       = kNDeltaEtaBins;
 204   dBinsPair[1]      = deltaEtaBins;
 205   axisTitlePair[1]  = "#Delta#eta";
 206
 207    // delta phi
 208   const Int_t kNDeltaPhiBins = 72;
 209   Double_t deltaPhiBins[kNDeltaPhiBins+1];
 210   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
 211     //deltaPhiBins[i] = -180.0 + i * 5.;
 212     deltaPhiBins[i] = -TMath::Pi()/2. + i * 5.*TMath::Pi()/180.;
 213   }
 214   iBinPair[2]       = kNDeltaPhiBins;
 215   dBinsPair[2]      = deltaPhiBins;
 216   axisTitlePair[2]  = "#Delta#varphi (rad)";
 217
 218   // pt(trigger-associated)
 219   const Int_t kNPtBins = 16;
 220   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.};
 221   //for(Int_t i = 0; i < kNPtBins+1; i++){
 222   //ptBins[i] = 0.2 + i * 0.5;
 223   //}
 224   iBinSingle[1]       = kNPtBins;
 225   dBinsSingle[1]      = ptBins;
 226   axisTitleSingle[1]  = "p_{T,trig.} (GeV/c)";
 227
 228   iBinPair[3]       = kNPtBins;
 229   dBinsPair[3]      = ptBins;
 230   axisTitlePair[3]  = "p_{T,trig.} (GeV/c)";
 231
 232   iBinPair[4]       = kNPtBins;
 233   dBinsPair[4]      = ptBins;
 234   axisTitlePair[4]  = "p_{T,assoc.} (GeV/c)";
 235
 236   TString histName;
 237   //+ triggered particles
 238   histName = "fHistP";
 239   if(fShuffle) histName.Append("_shuffle");
 240   if(fCentralityId) histName += fCentralityId.Data();
 241   fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
 242   for (Int_t j=0; j<kTrackVariablesSingle; j++) {
 243     fHistP->SetBinLimits(j, dBinsSingle[j]);
 244     fHistP->SetVarTitle(j, axisTitleSingle[j]);
 245   }
 246
 247   //- triggered particles
 248   histName = "fHistN";
 249   if(fShuffle) histName.Append("_shuffle");
 250   if(fCentralityId) histName += fCentralityId.Data();
 251   fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
 252   for (Int_t j=0; j<kTrackVariablesSingle; j++) {
 253     fHistN->SetBinLimits(j, dBinsSingle[j]);
 254     fHistN->SetVarTitle(j, axisTitleSingle[j]);
 255   }
 256
 257   //+- pairs
 258   histName = "fHistPN";
 259   if(fShuffle) histName.Append("_shuffle");
 260   if(fCentralityId) histName += fCentralityId.Data();
 261   fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 262   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 263     fHistPN->SetBinLimits(j, dBinsPair[j]);
 264     fHistPN->SetVarTitle(j, axisTitlePair[j]);
 265   }
 266
 267   //-+ pairs
 268   histName = "fHistNP";
 269   if(fShuffle) histName.Append("_shuffle");
 270   if(fCentralityId) histName += fCentralityId.Data();
 271   fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 272   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 273     fHistNP->SetBinLimits(j, dBinsPair[j]);
 274     fHistNP->SetVarTitle(j, axisTitlePair[j]);
 275   }
 276
 277   //++ pairs
 278   histName = "fHistPP";
 279   if(fShuffle) histName.Append("_shuffle");
 280   if(fCentralityId) histName += fCentralityId.Data();
 281   fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 282   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 283     fHistPP->SetBinLimits(j, dBinsPair[j]);
 284     fHistPP->SetVarTitle(j, axisTitlePair[j]);
 285   }
 286
 287   //-- pairs
 288   histName = "fHistNN";
 289   if(fShuffle) histName.Append("_shuffle");
 290   if(fCentralityId) histName += fCentralityId.Data();
 291   fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 292   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 293     fHistNN->SetBinLimits(j, dBinsPair[j]);
 294     fHistNN->SetVarTitle(j, axisTitlePair[j]);
 295   }
 296   AliInfo("Finished setting up the AliTHn");
 297
 298   // QA histograms
 299   fHistHBTbefore        = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
 300   fHistHBTafter         = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
 301   fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,2.*TMath::Pi());
 302   fHistConversionafter  = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,2.*TMath::Pi());
 303   fHistPsiMinusPhi     = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
 304
 305   TH1::AddDirectory(oldStatus);
 306
 307 }
 308
 309 //____________________________________________________________________//
 310 void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
 311                                      TObjArray *particles,
 312                                      TObjArray *particlesMixed,
 313                      Float_t bSign,
 314                      Double_t kMultorCent) {
 315   // Calculates the balance function
 316   fAnalyzedEvents++;
 317
 318   // Initialize histograms if not done yet
 319   if(!fHistPN){
 320     AliWarning("Histograms not yet initialized! --> Will be done now");
 321     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 322     InitHistograms();
 323   }
 324
 325   Double_t trackVariablesSingle[kTrackVariablesSingle];
 326   Double_t trackVariablesPair[kTrackVariablesPair];
 327
 328   if (!particles){
 329     AliWarning("particles TObjArray is NULL pointer --> return");
 330     return;
 331   }
 332
 333   // define end of particle loops
 334   Int_t iMax = particles->GetEntriesFast();
 335   Int_t jMax = iMax;
 336   if (particlesMixed)
 337     jMax = particlesMixed->GetEntriesFast();
 338
 339   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
 340   TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
 341
 342   TArrayF secondEta(jMax);
 343   TArrayF secondPhi(jMax);
 344   TArrayF secondPt(jMax);
 345   TArrayS secondCharge(jMax);
 346   TArrayD secondCorrection(jMax);
 347
 348   for (Int_t i=0; i<jMax; i++){
 349     secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
 350     secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
 351     secondPt[i]  = ((AliVParticle*) particlesSecond->At(i))->Pt();
 352     secondCharge[i]  = (Short_t)((AliVParticle*) particlesSecond->At(i))->Charge();
 353     secondCorrection[i]  = (Double_t)((AliBFBasicParticle*) particlesSecond->At(i))->Correction();   //==========================correction
 354   }
 355
 356   // 1st particle loop
 357   for (Int_t i = 0; i < iMax; i++) {
 358     //AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
 359     AliBFBasicParticle* firstParticle = (AliBFBasicParticle*) particles->At(i); //==========================correction
 360
 361     // some optimization
 362     Float_t firstEta = firstParticle->Eta();
 363     Float_t firstPhi = firstParticle->Phi();
 364     Float_t firstPt  = firstParticle->Pt();
 365     Float_t firstCorrection  = firstParticle->Correction();//==========================correction
 366
 367     // Event plane (determine psi bin)
 368     Double_t gPsiMinusPhi    =   0.;
 369     Double_t gPsiMinusPhiBin = -10.;
 370     gPsiMinusPhi   = TMath::Abs(firstPhi - gReactionPlane);
 371     //in-plane
 372     if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
 373        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
 374       gPsiMinusPhiBin = 0.0;
 375     //intermediate
 376     else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
 377             ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
 378             ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
 379             ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
 380       gPsiMinusPhiBin = 1.0;
 381     //out of plane
 382     else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
 383             ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
 384       gPsiMinusPhiBin = 2.0;
 385     //everything else
 386     else
 387       gPsiMinusPhiBin = 3.0;
 388
 389     fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
 390
 391     Short_t  charge1 = (Short_t) firstParticle->Charge();
 392
 393     trackVariablesSingle[0]    =  gPsiMinusPhiBin;
 394     trackVariablesSingle[1]    =  firstPt;
 395       if(fEventClass=="Multiplicity" || fEventClass == "Centrality" ) trackVariablesSingle[0] = kMultorCent;
 396
 397     //fill single particle histograms
 398     if(charge1 > 0)      fHistP->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
 399     else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,firstCorrection);  //==========================correction
 400
 401     // 2nd particle loop
 402     for(Int_t j = 0; j < jMax; j++) {
 403
 404       if(!particlesMixed && j == i) continue; // no auto correlations (only for non mixing)
 405
 406       // pT,Assoc < pT,Trig
 407       if(firstPt < secondPt[j]) continue;
 408
 409       Short_t charge2 = secondCharge[j];
 410
 411       trackVariablesPair[0]    =  trackVariablesSingle[0];
 412       trackVariablesPair[1]    =  firstEta - secondEta[j];  // delta eta
 413       trackVariablesPair[2]    =  firstPhi - secondPhi[j];  // delta phi
 414       //if (trackVariablesPair[2] > 180.)   // delta phi between -180 and 180
 415       //trackVariablesPair[2] -= 360.;
 416       //if (trackVariablesPair[2] <  - 180.)
 417       //trackVariablesPair[2] += 360.;
 418       if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi
 419         trackVariablesPair[2] -= 2.*TMath::Pi();
 420       if (trackVariablesPair[2] <  - TMath::Pi())
 421         trackVariablesPair[2] += 2.*TMath::Pi();
 422       if (trackVariablesPair[2] <  - TMath::Pi()/2.)
 423       trackVariablesPair[2] += 2.*TMath::Pi();
 424
 425       trackVariablesPair[3]    =  firstPt;      // pt trigger
 426       trackVariablesPair[4]    =  secondPt[j];  // pt
 427       //        trackVariablesPair[5]    =  fCentrality;  // centrality
 428
 429       // HBT like cut
 430       if(fHBTCut){ // VERSION 3 (all pairs)
 431         //if(fHBTCut && charge1 * charge2 > 0){  // VERSION 2 (only for LS)
 432         //if( dphi < 3 || deta < 0.01 ){   // VERSION 1
 433         //  continue;
 434
 435         Double_t deta = firstEta - secondEta[j];
 436         Double_t dphi = firstPhi - secondPhi[j];
 437         // VERSION 2 (Taken from DPhiCorrelations)
 438         // the variables & cuthave been developed by the HBT group
 439         // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
 440         fHistHBTbefore->Fill(deta,dphi);
 441
 442         // optimization
 443         if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
 444           {
 445             // phi in rad
 446             //Float_t phi1rad = firstPhi*TMath::DegToRad();
 447             //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
 448             Float_t phi1rad = firstPhi;
 449             Float_t phi2rad = secondPhi[j];
 450
 451             // check first boundaries to see if is worth to loop and find the minimum
 452             Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
 453             Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
 454
 455             const Float_t kLimit = 0.02 * 3;
 456
 457             Float_t dphistarminabs = 1e5;
 458             Float_t dphistarmin = 1e5;
 459
 460             if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
 461               for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
 462                 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
 463                 Float_t dphistarabs = TMath::Abs(dphistar);
 464
 465                 if (dphistarabs < dphistarminabs) {
 466                   dphistarmin = dphistar;
 467                   dphistarminabs = dphistarabs;
 468                 }
 469               }
 470
 471               if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
 472                 //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));
 473                 continue;
 474               }
 475             }
 476           }
 477         fHistHBTafter->Fill(deta,dphi);
 478       }//HBT cut
 479
 480       // conversions
 481       if(fConversionCut) {
 482         if (charge1 * charge2 < 0) {
 483           Double_t deta = firstEta - secondEta[j];
 484           Double_t dphi = firstPhi - secondPhi[j];
 485           fHistConversionbefore->Fill(deta,dphi);
 486
 487           Float_t m0 = 0.510e-3;
 488           Float_t tantheta1 = 1e10;
 489
 490           // phi in rad
 491           //Float_t phi1rad = firstPhi*TMath::DegToRad();
 492           //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
 493           Float_t phi1rad = firstPhi;
 494           Float_t phi2rad = secondPhi[j];
 495
 496           if (firstEta < -1e-10 || firstEta > 1e-10)
 497             tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
 498
 499           Float_t tantheta2 = 1e10;
 500           if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
 501             tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
 502
 503           Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
 504           Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
 505
 506           Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
 507
 508           if (masssqu < 0.04*0.04){
 509             //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));
 510             continue;
 511           }
 512           fHistConversionafter->Fill(deta,dphi);
 513         }
 514       }//conversion cut
 515
 516       if( charge1 > 0 && charge2 < 0)  fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
 517       else if( charge1 < 0 && charge2 > 0)  fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
 518       else if( charge1 > 0 && charge2 > 0)  fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
 519       else if( charge1 < 0 && charge2 < 0)  fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
 520       else {
 521         //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
 522         continue;
 523       }
 524     }//end of 2nd particle loop
 525   }//end of 1st particle loop
 526 }
 527
 528 //____________________________________________________________________//
 529 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
 530                                                  Int_t iVariablePair,
 531                                                  Double_t psiMin,
 532                                                  Double_t psiMax,
 533                                                  Double_t ptTriggerMin,
 534                                                  Double_t ptTriggerMax,
 535                                                  Double_t ptAssociatedMin,
 536                                                  Double_t ptAssociatedMax) {
 537   //Returns the BF histogram, extracted from the 6 AliTHn objects
 538   //(private members) of the AliBalancePsi class.
 539   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 540   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 541   // Psi_2
 542   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 543   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 544   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 545   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 546   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 547   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 548
 549   // pt trigger
 550   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 551     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 552     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 553     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 554     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 555     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 556     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 557   }
 558
 559   // pt associated
 560   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 561     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 562     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 563     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 564     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 565   }
 566
 567   //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));
 568
 569   // Project into the wanted space (1st: analysis step, 2nd: axis)
 570   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
 571   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
 572   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
 573   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
 574   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
 575   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
 576
 577   TH1D *gHistBalanceFunctionHistogram = 0x0;
 578   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 579     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
 580     gHistBalanceFunctionHistogram->Reset();
 581
 582     switch(iVariablePair) {
 583     case 1:
 584       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
 585       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
 586       break;
 587     case 2:
 588       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
 589       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
 590       break;
 591     default:
 592       break;
 593     }
 594
 595     hTemp1->Sumw2();
 596     hTemp2->Sumw2();
 597     hTemp3->Sumw2();
 598     hTemp4->Sumw2();
 599     hTemp1->Add(hTemp3,-1.);
 600     hTemp1->Scale(1./hTemp5->GetEntries());
 601     hTemp2->Add(hTemp4,-1.);
 602     hTemp2->Scale(1./hTemp6->GetEntries());
 603     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 604     gHistBalanceFunctionHistogram->Scale(0.5);
 605
 606     //normalize to bin width
 607     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
 608   }
 609
 610   return gHistBalanceFunctionHistogram;
 611 }
 612
 613 //____________________________________________________________________//
 614 TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
 615                                                          Int_t iVariablePair,
 616                                                          Double_t psiMin,
 617                                                          Double_t psiMax,
 618                                                          Double_t ptTriggerMin,
 619                                                          Double_t ptTriggerMax,
 620                                                          Double_t ptAssociatedMin,
 621                                                          Double_t ptAssociatedMax,
 622                                                          AliBalancePsi *bfMix) {
 623   //Returns the BF histogram, extracted from the 6 AliTHn objects
 624   //after dividing each correlation function by the Event Mixing one
 625   //(private members) of the AliBalancePsi class.
 626   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 627   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 628   // Psi_2
 629   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 630   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 631   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 632   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 633   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 634   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 635
 636   // pt trigger
 637   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 638     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 639     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 640     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 641     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 642     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 643     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 644   }
 645
 646   // pt associated
 647   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 648     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 649     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 650     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 651     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 652   }
 653
 654   // ============================================================================================
 655   // the same for event mixing
 656   AliTHn *fHistPMix = bfMix->GetHistNp();
 657   AliTHn *fHistNMix = bfMix->GetHistNn();
 658   AliTHn *fHistPNMix = bfMix->GetHistNpn();
 659   AliTHn *fHistNPMix = bfMix->GetHistNnp();
 660   AliTHn *fHistPPMix = bfMix->GetHistNpp();
 661   AliTHn *fHistNNMix = bfMix->GetHistNnn();
 662
 663   // Psi_2
 664   fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 665   fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 666   fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 667   fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 668   fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 669   fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 670
 671   // pt trigger
 672   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 673     fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 674     fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 675     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 676     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 677     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 678     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 679   }
 680
 681   // pt associated
 682   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 683     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 684     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 685     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 686     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 687   }
 688   // ============================================================================================
 689
 690   //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));
 691
 692   // Project into the wanted space (1st: analysis step, 2nd: axis)
 693   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
 694   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
 695   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
 696   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
 697   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
 698   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
 699
 700   // ============================================================================================
 701   // the same for event mixing
 702   TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,iVariablePair);
 703   TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,iVariablePair);
 704   TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,iVariablePair);
 705   TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,iVariablePair);
 706   TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
 707   TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
 708   // ============================================================================================
 709
 710   TH1D *gHistBalanceFunctionHistogram = 0x0;
 711   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 712     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
 713     gHistBalanceFunctionHistogram->Reset();
 714
 715     switch(iVariablePair) {
 716     case 1:
 717       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
 718       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
 719       break;
 720     case 2:
 721       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
 722       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
 723       break;
 724     default:
 725       break;
 726     }
 727
 728     hTemp1->Sumw2();
 729     hTemp2->Sumw2();
 730     hTemp3->Sumw2();
 731     hTemp4->Sumw2();
 732     hTemp1Mix->Sumw2();
 733     hTemp2Mix->Sumw2();
 734     hTemp3Mix->Sumw2();
 735     hTemp4Mix->Sumw2();
 736
 737     hTemp1->Scale(1./hTemp5->GetEntries());
 738     hTemp3->Scale(1./hTemp5->GetEntries());
 739     hTemp2->Scale(1./hTemp6->GetEntries());
 740     hTemp4->Scale(1./hTemp6->GetEntries());
 741     hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
 742     hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
 743     hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
 744     hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
 745
 746     hTemp1->Divide(hTemp1Mix);
 747     hTemp2->Divide(hTemp2Mix);
 748     hTemp3->Divide(hTemp3Mix);
 749     hTemp4->Divide(hTemp4Mix);
 750
 751     hTemp1->Add(hTemp3,-1.);
 752     hTemp2->Add(hTemp4,-1.);
 753
 754     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 755     gHistBalanceFunctionHistogram->Scale(0.5);
 756
 757     //normalize to bin width
 758     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
 759   }
 760
 761   return gHistBalanceFunctionHistogram;
 762 }
 763
 764 //____________________________________________________________________//
 765 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
 766                                                         Double_t psiMax,
 767                                                         Double_t ptTriggerMin,
 768                                                         Double_t ptTriggerMax,
 769                                                         Double_t ptAssociatedMin,
 770                                                         Double_t ptAssociatedMax) {
 771   //Returns the BF histogram in Delta eta vs Delta phi,
 772   //extracted from the 6 AliTHn objects
 773   //(private members) of the AliBalancePsi class.
 774   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 775   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 776   TString histName = "gHistBalanceFunctionHistogram2D";
 777
 778   // Psi_2
 779   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 780   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 781   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 782   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 783   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 784   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 785
 786   // pt trigger
 787   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 788     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 789     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 790     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 791     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 792     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 793     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 794   }
 795
 796   // pt associated
 797   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 798     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 799     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 800     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 801     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 802   }
 803
 804   //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)));
 805
 806   // Project into the wanted space (1st: analysis step, 2nd: axis)
 807   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
 808   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
 809   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
 810   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
 811   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
 812   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
 813
 814   TH2D *gHistBalanceFunctionHistogram = 0x0;
 815   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 816     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
 817     gHistBalanceFunctionHistogram->Reset();
 818     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
 819     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
 820     gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
 821
 822     hTemp1->Sumw2();
 823     hTemp2->Sumw2();
 824     hTemp3->Sumw2();
 825     hTemp4->Sumw2();
 826     hTemp1->Add(hTemp3,-1.);
 827     hTemp1->Scale(1./hTemp5->GetEntries());
 828     hTemp2->Add(hTemp4,-1.);
 829     hTemp2->Scale(1./hTemp6->GetEntries());
 830     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 831     gHistBalanceFunctionHistogram->Scale(0.5);
 832
 833     //normalize to bin width
 834     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
 835   }
 836
 837   return gHistBalanceFunctionHistogram;
 838 }
 839
 840 //____________________________________________________________________//
 841 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
 842                                                                 Double_t psiMax,
 843                                                                 Double_t ptTriggerMin,
 844                                                                 Double_t ptTriggerMax,
 845                                                                 Double_t ptAssociatedMin,
 846                                                                 Double_t ptAssociatedMax,
 847                                                                 AliBalancePsi *bfMix) {
 848   //Returns the BF histogram in Delta eta vs Delta phi,
 849   //after dividing each correlation function by the Event Mixing one
 850   //extracted from the 6 AliTHn objects
 851   //(private members) of the AliBalancePsi class.
 852   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 853   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 854   TString histName = "gHistBalanceFunctionHistogram2D";
 855
 856   if(!bfMix){
 857     AliError("balance function object for event mixing not available");
 858     return NULL;
 859   }
 860
 861   // Psi_2
 862   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 863   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 864   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 865   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 866   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 867   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 868
 869   // pt trigger
 870   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 871     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 872     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 873     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 874     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 875     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 876     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 877   }
 878
 879   // pt associated
 880   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 881     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 882     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 883     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 884     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 885   }
 886
 887
 888   // ============================================================================================
 889   // the same for event mixing
 890   AliTHn *fHistPMix = bfMix->GetHistNp();
 891   AliTHn *fHistNMix = bfMix->GetHistNn();
 892   AliTHn *fHistPNMix = bfMix->GetHistNpn();
 893   AliTHn *fHistNPMix = bfMix->GetHistNnp();
 894   AliTHn *fHistPPMix = bfMix->GetHistNpp();
 895   AliTHn *fHistNNMix = bfMix->GetHistNnn();
 896
 897   // Psi_2
 898   fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 899   fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 900   fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 901   fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 902   fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 903   fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 904
 905   // pt trigger
 906   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 907     fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 908     fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 909     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 910     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 911     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 912     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 913   }
 914
 915   // pt associated
 916   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 917     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 918     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 919     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 920     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 921   }
 922   // ============================================================================================
 923
 924
 925   //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)));
 926
 927   // Project into the wanted space (1st: analysis step, 2nd: axis)
 928   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
 929   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
 930   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
 931   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
 932   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
 933   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
 934
 935   // ============================================================================================
 936   // the same for event mixing
 937   TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
 938   TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
 939   TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
 940   TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
 941   TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
 942   TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
 943   // ============================================================================================
 944
 945   TH2D *gHistBalanceFunctionHistogram = 0x0;
 946   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 947     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
 948     gHistBalanceFunctionHistogram->Reset();
 949     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
 950     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
 951     gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
 952
 953     hTemp1->Sumw2();
 954     hTemp2->Sumw2();
 955     hTemp3->Sumw2();
 956     hTemp4->Sumw2();
 957     hTemp1Mix->Sumw2();
 958     hTemp2Mix->Sumw2();
 959     hTemp3Mix->Sumw2();
 960     hTemp4Mix->Sumw2();
 961
 962     hTemp1->Scale(1./hTemp5->GetEntries());
 963     hTemp3->Scale(1./hTemp5->GetEntries());
 964     hTemp2->Scale(1./hTemp6->GetEntries());
 965     hTemp4->Scale(1./hTemp6->GetEntries());
 966     hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
 967     hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
 968     hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
 969     hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
 970
 971     hTemp1->Divide(hTemp1Mix);
 972     hTemp2->Divide(hTemp2Mix);
 973     hTemp3->Divide(hTemp3Mix);
 974     hTemp4->Divide(hTemp4Mix);
 975
 976     hTemp1->Add(hTemp3,-1.);
 977     hTemp2->Add(hTemp4,-1.);
 978
 979     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 980     gHistBalanceFunctionHistogram->Scale(0.5);
 981
 982     //normalize to bin width
 983     gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
 984   }
 985
 986   return gHistBalanceFunctionHistogram;
 987 }
 988
 989 //____________________________________________________________________//
 990 TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
 991                                                         Double_t psiMin,
 992                                                         Double_t psiMax,
 993                                                         Double_t ptTriggerMin,
 994                                                         Double_t ptTriggerMax,
 995                                                         Double_t ptAssociatedMin,
 996                                                         Double_t ptAssociatedMax,
 997                                                         AliBalancePsi *bfMix) {
 998   //Returns the BF histogram in Delta eta OR Delta phi,
 999   //after dividing each correlation function by the Event Mixing one
1000   // (But the division is done here in 2D, this was basically done to check the Integral)
1001   //extracted from the 6 AliTHn objects
1002   //(private members) of the AliBalancePsi class.
1003   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1004   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1005   TString histName = "gHistBalanceFunctionHistogram1D";
1006
1007   if(!bfMix){
1008     AliError("balance function object for event mixing not available");
1009     return NULL;
1010   }
1011
1012   // Psi_2
1013   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1014   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1015   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1016   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1017   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1018   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1019
1020   // pt trigger
1021   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1022     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1023     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1024     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1025     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1026     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1027     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1028   }
1029
1030   // pt associated
1031   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1032     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1033     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1034     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1035     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1036   }
1037
1038
1039   // ============================================================================================
1040   // the same for event mixing
1041   AliTHn *fHistPMix = bfMix->GetHistNp();
1042   AliTHn *fHistNMix = bfMix->GetHistNn();
1043   AliTHn *fHistPNMix = bfMix->GetHistNpn();
1044   AliTHn *fHistNPMix = bfMix->GetHistNnp();
1045   AliTHn *fHistPPMix = bfMix->GetHistNpp();
1046   AliTHn *fHistNNMix = bfMix->GetHistNnn();
1047
1048   // Psi_2
1049   fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1050   fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1051   fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1052   fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1053   fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1054   fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1055
1056   // pt trigger
1057   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1058     fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1059     fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1060     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1061     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1062     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1063     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1064   }
1065
1066   // pt associated
1067   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1068     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1069     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1070     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1071     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1072   }
1073   // ============================================================================================
1074
1075
1076   //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)));
1077
1078   // Project into the wanted space (1st: analysis step, 2nd: axis)
1079   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1080   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1081   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1082   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1083   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1084   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1085
1086   // ============================================================================================
1087   // the same for event mixing
1088   TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1089   TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1090   TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1091   TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1092   TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1093   TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1094   // ============================================================================================
1095
1096   TH1D *gHistBalanceFunctionHistogram = 0x0;
1097   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1098
1099     hTemp1->Sumw2();
1100     hTemp2->Sumw2();
1101     hTemp3->Sumw2();
1102     hTemp4->Sumw2();
1103     hTemp1Mix->Sumw2();
1104     hTemp2Mix->Sumw2();
1105     hTemp3Mix->Sumw2();
1106     hTemp4Mix->Sumw2();
1107
1108     hTemp1->Scale(1./hTemp5->GetEntries());
1109     hTemp3->Scale(1./hTemp5->GetEntries());
1110     hTemp2->Scale(1./hTemp6->GetEntries());
1111     hTemp4->Scale(1./hTemp6->GetEntries());
1112
1113     hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1114     hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1115     hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1116     hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
1117
1118     hTemp1->Divide(hTemp1Mix);
1119     hTemp2->Divide(hTemp2Mix);
1120     hTemp3->Divide(hTemp3Mix);
1121     hTemp4->Divide(hTemp4Mix);
1122
1123     // now only project on one axis
1124     TH1D *h1DTemp1 = NULL;
1125     TH1D *h1DTemp2 = NULL;
1126     TH1D *h1DTemp3 = NULL;
1127     TH1D *h1DTemp4 = NULL;
1128     if(!bPhi){
1129       h1DTemp1 = (TH1D*)hTemp1->ProjectionX(Form("%s_projX",hTemp1->GetName()));
1130       h1DTemp2 = (TH1D*)hTemp2->ProjectionX(Form("%s_projX",hTemp2->GetName()));
1131       h1DTemp3 = (TH1D*)hTemp3->ProjectionX(Form("%s_projX",hTemp3->GetName()));
1132       h1DTemp4 = (TH1D*)hTemp4->ProjectionX(Form("%s_projX",hTemp4->GetName()));
1133     }
1134     else{
1135       h1DTemp1 = (TH1D*)hTemp1->ProjectionY(Form("%s_projX",hTemp1->GetName()));
1136       h1DTemp2 = (TH1D*)hTemp2->ProjectionY(Form("%s_projX",hTemp2->GetName()));
1137       h1DTemp3 = (TH1D*)hTemp3->ProjectionY(Form("%s_projX",hTemp3->GetName()));
1138       h1DTemp4 = (TH1D*)hTemp4->ProjectionY(Form("%s_projX",hTemp4->GetName()));
1139     }
1140
1141     gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1142     gHistBalanceFunctionHistogram->Reset();
1143     if(!bPhi){
1144       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1145       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1146     }
1147     else{
1148       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1149       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1150     }
1151
1152     h1DTemp1->Add(h1DTemp3,-1.);
1153     h1DTemp2->Add(h1DTemp4,-1.);
1154
1155     gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1156     gHistBalanceFunctionHistogram->Scale(0.5);
1157   }
1158
1159   return gHistBalanceFunctionHistogram;
1160 }
1161
1162
1163 //____________________________________________________________________//
1164 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
1165                                               Double_t psiMax,
1166                                               Double_t ptTriggerMin,
1167                                               Double_t ptTriggerMax,
1168                                               Double_t ptAssociatedMin,
1169                                               Double_t ptAssociatedMax) {
1170   //Returns the 2D correlation function for (+-) pairs
1171   // Psi_2: axis 0
1172   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1173   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1174
1175   // pt trigger
1176   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1177     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1178     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1179   }
1180
1181   // pt associated
1182   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1183     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1184
1185   //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1186   //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1187
1188   //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
1189   //TCanvas *c1 = new TCanvas("c1","");
1190   //c1->cd();
1191   //if(!gHistTest){
1192   //AliError("Projection of fHistP = NULL");
1193   //return gHistTest;
1194   //}
1195   //else{
1196   //gHistTest->DrawCopy("colz");
1197   //}
1198
1199   //0:step, 1: Delta eta, 2: Delta phi
1200   TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
1201   if(!gHist){
1202     AliError("Projection of fHistPN = NULL");
1203     return gHist;
1204   }
1205
1206   //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
1207   //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
1208   //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
1209
1210   //TCanvas *c2 = new TCanvas("c2","");
1211   //c2->cd();
1212   //fHistPN->Project(0,1,2)->DrawCopy("colz");
1213
1214   if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
1215     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
1216
1217   //normalize to bin width
1218   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
1219
1220   return gHist;
1221 }
1222
1223 //____________________________________________________________________//
1224 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
1225                                               Double_t psiMax,
1226                                               Double_t ptTriggerMin,
1227                                               Double_t ptTriggerMax,
1228                                               Double_t ptAssociatedMin,
1229                                               Double_t ptAssociatedMax) {
1230   //Returns the 2D correlation function for (+-) pairs
1231   // Psi_2: axis 0
1232   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1233   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1234
1235   // pt trigger
1236   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1237     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1238     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1239   }
1240
1241   // pt associated
1242   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1243     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1244
1245   //0:step, 1: Delta eta, 2: Delta phi
1246   TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
1247   if(!gHist){
1248     AliError("Projection of fHistPN = NULL");
1249     return gHist;
1250   }
1251
1252   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1253   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
1254   if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
1255     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
1256
1257   //normalize to bin width
1258   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
1259
1260   return gHist;
1261 }
1262
1263 //____________________________________________________________________//
1264 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
1265                                               Double_t psiMax,
1266                                               Double_t ptTriggerMin,
1267                                               Double_t ptTriggerMax,
1268                                               Double_t ptAssociatedMin,
1269                                               Double_t ptAssociatedMax) {
1270   //Returns the 2D correlation function for (+-) pairs
1271   // Psi_2: axis 0
1272   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1273   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1274
1275   // pt trigger
1276   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1277     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1278     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1279   }
1280
1281   // pt associated
1282   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1283     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1284
1285   //0:step, 1: Delta eta, 2: Delta phi
1286   TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
1287   if(!gHist){
1288     AliError("Projection of fHistPN = NULL");
1289     return gHist;
1290   }
1291
1292   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
1293   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
1294   if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
1295     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
1296
1297   //normalize to bin width
1298   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
1299
1300   return gHist;
1301 }
1302
1303 //____________________________________________________________________//
1304 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
1305                                               Double_t psiMax,
1306                                               Double_t ptTriggerMin,
1307                                               Double_t ptTriggerMax,
1308                                               Double_t ptAssociatedMin,
1309                                               Double_t ptAssociatedMax) {
1310   //Returns the 2D correlation function for (+-) pairs
1311   // Psi_2: axis 0
1312   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1313   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1314
1315   // pt trigger
1316   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1317     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1318     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1319   }
1320
1321   // pt associated
1322   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1323     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1324
1325   //0:step, 1: Delta eta, 2: Delta phi
1326   TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
1327   if(!gHist){
1328     AliError("Projection of fHistPN = NULL");
1329     return gHist;
1330   }
1331
1332   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1333   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
1334   if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
1335     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
1336
1337   //normalize to bin width
1338   gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
1339
1340   return gHist;
1341 }
1342
1343 //____________________________________________________________________//
1344 TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
1345                                                              Double_t psiMax,
1346                                                              Double_t ptTriggerMin,
1347                                                              Double_t ptTriggerMax,
1348                                                              Double_t ptAssociatedMin,
1349                                                              Double_t ptAssociatedMax) {
1350   //Returns the 2D correlation function for the sum of all charge combination pairs
1351   // Psi_2: axis 0
1352   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1353   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1354   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1355   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1356   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1357   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1358
1359   // pt trigger
1360   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1361     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1362     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1363     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1364     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1365     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1366     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1367   }
1368
1369   // pt associated
1370   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1371     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1372     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1373     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1374     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1375
1376   //0:step, 1: Delta eta, 2: Delta phi
1377   TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
1378   if(!gHistNN){
1379     AliError("Projection of fHistNN = NULL");
1380     return gHistNN;
1381   }
1382   TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
1383   if(!gHistPP){
1384     AliError("Projection of fHistPP = NULL");
1385     return gHistPP;
1386   }
1387   TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
1388   if(!gHistNP){
1389     AliError("Projection of fHistNP = NULL");
1390     return gHistNP;
1391   }
1392   TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
1393   if(!gHistPN){
1394     AliError("Projection of fHistPN = NULL");
1395     return gHistPN;
1396   }
1397
1398   // sum all 2 particle histograms
1399   gHistNN->Add(gHistPP);
1400   gHistNN->Add(gHistNP);
1401   gHistNN->Add(gHistPN);
1402
1403   // divide by sum of + and - triggers
1404   if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
1405     gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
1406
1407   //normalize to bin width
1408   gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
1409
1410   return gHistNN;
1411 }
1412
1413 //____________________________________________________________________//
1414
1415 Bool_t AliBalancePsi::GetMomentsAnalytical(TH1D* gHist,
1416                                            Double_t &mean, Double_t &meanError,
1417                                            Double_t &sigma, Double_t &sigmaError,
1418                                            Double_t &skewness, Double_t &skewnessError,
1419                                            Double_t &kurtosis, Double_t &kurtosisError) {
1420   //
1421   // helper method to calculate the moments and errors of a TH1D anlytically
1422   //
1423
1424   Bool_t success = kFALSE;
1425   mean          = -1.;
1426   meanError     = -1.;
1427   sigma         = -1.;
1428   sigmaError    = -1.;
1429   skewness      = -1.;
1430   skewnessError = -1.;
1431   kurtosis      = -1.;
1432   kurtosisError = -1.;
1433
1434   if(gHist){
1435
1436     // ----------------------------------------------------------------------
1437     // basic parameters of histogram
1438
1439     Int_t fNumberOfBins = gHist->GetNbinsX();
1440     //    Int_t fBinWidth     = gHist->GetBinWidth(1); // assume equal binning
1441
1442
1443     // ----------------------------------------------------------------------
1444     // first calculate the mean
1445
1446     Double_t fWeightedAverage   = 0.;
1447     Double_t fNormalization = 0.;
1448
1449     for(Int_t i = 1; i <= fNumberOfBins; i++) {
1450
1451       fWeightedAverage   += gHist->GetBinContent(i) * gHist->GetBinCenter(i);
1452       fNormalization     += gHist->GetBinContent(i);
1453
1454     }
1455
1456     mean = fWeightedAverage / fNormalization;
1457
1458
1459     // ----------------------------------------------------------------------
1460     // then calculate the higher moments
1461
1462     Double_t fDelta  = 0.;
1463     Double_t fDelta2 = 0.;
1464     Double_t fDelta3 = 0.;
1465     Double_t fDelta4 = 0.;
1466
1467     for(Int_t i = 1; i <= fNumberOfBins; i++) {
1468       fDelta  += gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean;
1469       fDelta2 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),2);
1470       fDelta3 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),3);
1471       fDelta4 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),4);
1472     }
1473
1474     sigma    = fDelta2 / fNormalization;
1475     skewness = fDelta3 / fNormalization / TMath::Power(sigma,3);
1476     kurtosis = fDelta4 / fNormalization / TMath::Power(sigma,4) - 3;
1477
1478     success = kTRUE;
1479   }
1480
1481
1482   return success;
1483 }
1484
1485
1486 //____________________________________________________________________//
1487 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) {
1488   //
1489   // calculates dphistar
1490   //
1491   Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1492
1493   static const Double_t kPi = TMath::Pi();
1494
1495   // circularity
1496 //   if (dphistar > 2 * kPi)
1497 //     dphistar -= 2 * kPi;
1498 //   if (dphistar < -2 * kPi)
1499 //     dphistar += 2 * kPi;
1500
1501   if (dphistar > kPi)
1502     dphistar = kPi * 2 - dphistar;
1503   if (dphistar < -kPi)
1504     dphistar = -kPi * 2 - dphistar;
1505   if (dphistar > kPi) // might look funny but is needed
1506     dphistar = kPi * 2 - dphistar;
1507
1508   return dphistar;
1509 }
1510