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   fHistPsiMinusPhi(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   fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
  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 fHistPsiMinusPhi;
 111 }
 112
 113 //____________________________________________________________________//
 114 void AliBalancePsi::InitHistograms() {
 115   // single particle histograms
 116
 117   // global switch disabling the reference
 118   // (to avoid "Replacing existing TH1" if several wagons are created in train)
 119   Bool_t oldStatus = TH1::AddDirectoryStatus();
 120   TH1::AddDirectory(kFALSE);
 121
 122   Int_t anaSteps   = 1;       // analysis steps
 123   Int_t iBinSingle[kTrackVariablesSingle];        // binning for track variables
 124   Double_t* dBinsSingle[kTrackVariablesSingle];   // bins for track variables
 125   TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
 126
 127   // two particle histograms
 128   Int_t iBinPair[kTrackVariablesPair];         // binning for track variables
 129   Double_t* dBinsPair[kTrackVariablesPair];    // bins for track variables
 130   TString axisTitlePair[kTrackVariablesPair];  // axis titles for track variables
 131
 132   //centrality
 133   /*const Int_t kNCentralityBins = 9;
 134   Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
 135   iBinSingle[0]       = kNCentralityBins;
 136   dBinsSingle[0]      = centralityBins;
 137   axisTitleSingle[0]  = "Centrality percentile [%]";
 138   iBinPair[0]       = kNCentralityBins;
 139   dBinsPair[0]      = centralityBins;
 140   axisTitlePair[0]  = "Centrality percentile [%]"; */
 141
 142   //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)
 143   const Int_t kNPsi2Bins = 4;
 144   Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
 145   iBinSingle[0]       = kNPsi2Bins;
 146   dBinsSingle[0]      = psi2Bins;
 147   axisTitleSingle[0]  = "#varphi - #Psi_{2} (a.u.)";
 148   iBinPair[0]       = kNPsi2Bins;
 149   dBinsPair[0]      = psi2Bins;
 150   axisTitlePair[0]  = "#varphi - #Psi_{2} (a.u.)";
 151
 152    // delta eta
 153   const Int_t kNDeltaEtaBins = 80;
 154   Double_t deltaEtaBins[kNDeltaEtaBins+1];
 155   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
 156     deltaEtaBins[i] = -2.0 + i * 0.05;
 157   iBinPair[1]       = kNDeltaEtaBins;
 158   dBinsPair[1]      = deltaEtaBins;
 159   axisTitlePair[1]  = "#Delta#eta";
 160
 161    // delta phi
 162   const Int_t kNDeltaPhiBins = 72;
 163   Double_t deltaPhiBins[kNDeltaPhiBins+1];
 164   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
 165     //deltaPhiBins[i] = -180.0 + i * 5.;
 166     deltaPhiBins[i] = -TMath::Pi()/2. + i * 5.*TMath::Pi()/180.;
 167   }
 168   iBinPair[2]       = kNDeltaPhiBins;
 169   dBinsPair[2]      = deltaPhiBins;
 170   axisTitlePair[2]  = "#Delta#varphi (rad)";
 171
 172   // pt(trigger-associated)
 173   const Int_t kNPtBins = 16;
 174   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.};
 175   //for(Int_t i = 0; i < kNPtBins+1; i++){
 176   //ptBins[i] = 0.2 + i * 0.5;
 177   //}
 178   iBinSingle[1]       = kNPtBins;
 179   dBinsSingle[1]      = ptBins;
 180   axisTitleSingle[1]  = "p_{T,trig.} (GeV/c)";
 181
 182   iBinPair[3]       = kNPtBins;
 183   dBinsPair[3]      = ptBins;
 184   axisTitlePair[3]  = "p_{T,trig.} (GeV/c)";
 185
 186   iBinPair[4]       = kNPtBins;
 187   dBinsPair[4]      = ptBins;
 188   axisTitlePair[4]  = "p_{T,assoc.} (GeV/c)";
 189
 190   TString histName;
 191   //+ triggered particles
 192   histName = "fHistP";
 193   if(fShuffle) histName.Append("_shuffle");
 194   if(fCentralityId) histName += fCentralityId.Data();
 195   fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
 196   for (Int_t j=0; j<kTrackVariablesSingle; j++) {
 197     fHistP->SetBinLimits(j, dBinsSingle[j]);
 198     fHistP->SetVarTitle(j, axisTitleSingle[j]);
 199   }
 200
 201   //- triggered particles
 202   histName = "fHistN";
 203   if(fShuffle) histName.Append("_shuffle");
 204   if(fCentralityId) histName += fCentralityId.Data();
 205   fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
 206   for (Int_t j=0; j<kTrackVariablesSingle; j++) {
 207     fHistN->SetBinLimits(j, dBinsSingle[j]);
 208     fHistN->SetVarTitle(j, axisTitleSingle[j]);
 209   }
 210
 211   //+- pairs
 212   histName = "fHistPN";
 213   if(fShuffle) histName.Append("_shuffle");
 214   if(fCentralityId) histName += fCentralityId.Data();
 215   fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 216   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 217     fHistPN->SetBinLimits(j, dBinsPair[j]);
 218     fHistPN->SetVarTitle(j, axisTitlePair[j]);
 219   }
 220
 221   //-+ pairs
 222   histName = "fHistNP";
 223   if(fShuffle) histName.Append("_shuffle");
 224   if(fCentralityId) histName += fCentralityId.Data();
 225   fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 226   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 227     fHistNP->SetBinLimits(j, dBinsPair[j]);
 228     fHistNP->SetVarTitle(j, axisTitlePair[j]);
 229   }
 230
 231   //++ pairs
 232   histName = "fHistPP";
 233   if(fShuffle) histName.Append("_shuffle");
 234   if(fCentralityId) histName += fCentralityId.Data();
 235   fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 236   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 237     fHistPP->SetBinLimits(j, dBinsPair[j]);
 238     fHistPP->SetVarTitle(j, axisTitlePair[j]);
 239   }
 240
 241   //-- pairs
 242   histName = "fHistNN";
 243   if(fShuffle) histName.Append("_shuffle");
 244   if(fCentralityId) histName += fCentralityId.Data();
 245   fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
 246   for (Int_t j=0; j<kTrackVariablesPair; j++) {
 247     fHistNN->SetBinLimits(j, dBinsPair[j]);
 248     fHistNN->SetVarTitle(j, axisTitlePair[j]);
 249   }
 250   AliInfo("Finished setting up the AliTHn");
 251
 252   // QA histograms
 253   fHistHBTbefore        = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
 254   fHistHBTafter         = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
 255   fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,2.*TMath::Pi());
 256   fHistConversionafter  = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,2.*TMath::Pi());
 257   fHistPsiMinusPhi     = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
 258
 259   TH1::AddDirectory(oldStatus);
 260
 261 }
 262
 263 //____________________________________________________________________//
 264 void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
 265                                      TObjArray *particles,
 266                                      TObjArray *particlesMixed,
 267                                      Float_t bSign) {
 268   // Calculates the balance function
 269   fAnalyzedEvents++;
 270
 271   // Initialize histograms if not done yet
 272   if(!fHistPN){
 273     AliWarning("Histograms not yet initialized! --> Will be done now");
 274     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 275     InitHistograms();
 276   }
 277
 278   Double_t trackVariablesSingle[kTrackVariablesSingle];
 279   Double_t trackVariablesPair[kTrackVariablesPair];
 280
 281   if (!particles){
 282     AliWarning("particles TObjArray is NULL pointer --> return");
 283     return;
 284   }
 285
 286   // define end of particle loops
 287   Int_t iMax = particles->GetEntriesFast();
 288   Int_t jMax = iMax;
 289   if (particlesMixed)
 290     jMax = particlesMixed->GetEntriesFast();
 291
 292   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
 293   TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
 294
 295   TArrayF secondEta(jMax);
 296   TArrayF secondPhi(jMax);
 297   TArrayF secondPt(jMax);
 298
 299   for (Int_t i=0; i<jMax; i++){
 300     secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
 301     secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
 302     secondPt[i]  = ((AliVParticle*) particlesSecond->At(i))->Pt();
 303   }
 304
 305   // 1st particle loop
 306   for (Int_t i = 0; i < iMax; i++) {
 307     AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
 308
 309     // some optimization
 310     Float_t firstEta = firstParticle->Eta();
 311     Float_t firstPhi = firstParticle->Phi();
 312     Float_t firstPt  = firstParticle->Pt();
 313
 314     // Event plane (determine psi bin)
 315     Double_t gPsiMinusPhi    =   0.;
 316     Double_t gPsiMinusPhiBin = -10.;
 317     gPsiMinusPhi   = TMath::Abs(firstPhi - gReactionPlane);
 318     //in-plane
 319     if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
 320        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
 321       gPsiMinusPhiBin = 0.0;
 322     //intermediate
 323     else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
 324             ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
 325             ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
 326             ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
 327       gPsiMinusPhiBin = 1.0;
 328     //out of plane
 329     else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
 330             ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
 331       gPsiMinusPhiBin = 2.0;
 332     //everything else
 333     else
 334       gPsiMinusPhiBin = 3.0;
 335
 336     fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
 337
 338     Short_t  charge1 = (Short_t) firstParticle->Charge();
 339
 340     trackVariablesSingle[0]    =  gPsiMinusPhiBin;
 341     trackVariablesSingle[1]    =  firstPt;
 342
 343     //fill single particle histograms
 344     if(charge1 > 0)      fHistP->Fill(trackVariablesSingle,0,1.);
 345     else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,1.);
 346
 347     // 2nd particle loop (only for j < i for non double counting in the same pT region)
 348     // --> SAME pT region for trigger and assoc: NO double counting with this
 349     // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i
 350     //                          --> can be handled afterwards by using assoc. as trigger as well ?!
 351     for(Int_t j = 0; j < iMax; j++) {
 352       if (particlesMixed && j == jMax-1 )  // if the mixed track number is smaller than the main event one
 353         break;
 354
 355       if(j == i) continue; // no auto correlations
 356
 357       AliVParticle* secondParticle = (AliVParticle*) particlesSecond->At(j);
 358       Short_t charge2 = (Short_t) secondParticle->Charge();
 359
 360       trackVariablesPair[0]    =  gPsiMinusPhiBin;
 361       trackVariablesPair[1]    =  firstEta - secondEta[j];  // delta eta
 362       trackVariablesPair[2]    =  firstPhi - secondPhi[j];  // delta phi
 363       //if (trackVariablesPair[2] > 180.)   // delta phi between -180 and 180
 364       //trackVariablesPair[2] -= 360.;
 365       //if (trackVariablesPair[2] <  - 180.)
 366       //trackVariablesPair[2] += 360.;
 367       if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi
 368         trackVariablesPair[2] -= 2.*TMath::Pi();
 369       if (trackVariablesPair[2] <  - TMath::Pi())
 370         trackVariablesPair[2] += 2.*TMath::Pi();
 371       if (trackVariablesPair[2] <  - TMath::Pi()/2.)
 372       trackVariablesPair[2] += 2.*TMath::Pi();
 373
 374       trackVariablesPair[3]    =  firstPt;      // pt trigger
 375       trackVariablesPair[4]    =  secondPt[j];  // pt
 376       //        trackVariablesPair[5]    =  fCentrality;  // centrality
 377
 378       // HBT like cut
 379       if(fHBTCut && charge1 * charge2 > 0){
 380         //if( dphi < 3 || deta < 0.01 ){   // VERSION 1
 381         //  continue;
 382
 383         Double_t deta = firstEta - secondEta[j];
 384         Double_t dphi = firstPhi - secondPhi[j];
 385         // VERSION 2 (Taken from DPhiCorrelations)
 386         // the variables & cuthave been developed by the HBT group
 387         // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
 388         fHistHBTbefore->Fill(deta,dphi);
 389
 390         // optimization
 391         if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
 392           {
 393             // phi in rad
 394             //Float_t phi1rad = firstPhi*TMath::DegToRad();
 395             //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
 396             Float_t phi1rad = firstPhi;
 397             Float_t phi2rad = secondPhi[j];
 398
 399             // check first boundaries to see if is worth to loop and find the minimum
 400             Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
 401             Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
 402
 403             const Float_t kLimit = 0.02 * 3;
 404
 405             Float_t dphistarminabs = 1e5;
 406             Float_t dphistarmin = 1e5;
 407
 408             if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
 409               for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
 410                 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
 411                 Float_t dphistarabs = TMath::Abs(dphistar);
 412
 413                 if (dphistarabs < dphistarminabs) {
 414                   dphistarmin = dphistar;
 415                   dphistarminabs = dphistarabs;
 416                 }
 417               }
 418
 419               if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
 420                 //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));
 421                 continue;
 422               }
 423             }
 424           }
 425         fHistHBTafter->Fill(deta,dphi);
 426       }//HBT cut
 427
 428       // conversions
 429       if(fConversionCut) {
 430         if (charge1 * charge2 < 0) {
 431           Double_t deta = firstEta - secondEta[j];
 432           Double_t dphi = firstPhi - secondPhi[j];
 433           fHistConversionbefore->Fill(deta,dphi);
 434
 435           Float_t m0 = 0.510e-3;
 436           Float_t tantheta1 = 1e10;
 437
 438           // phi in rad
 439           //Float_t phi1rad = firstPhi*TMath::DegToRad();
 440           //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
 441           Float_t phi1rad = firstPhi;
 442           Float_t phi2rad = secondPhi[j];
 443
 444           if (firstEta < -1e-10 || firstEta > 1e-10)
 445             tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
 446
 447           Float_t tantheta2 = 1e10;
 448           if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
 449             tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
 450
 451           Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
 452           Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
 453
 454           Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
 455
 456           if (masssqu < 0.04*0.04){
 457             //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));
 458             continue;
 459           }
 460           fHistConversionafter->Fill(deta,dphi);
 461         }
 462       }//conversion cut
 463
 464       if( charge1 > 0 && charge2 < 0)  fHistPN->Fill(trackVariablesPair,0,1.);
 465       else if( charge1 < 0 && charge2 > 0)  fHistNP->Fill(trackVariablesPair,0,1.);
 466       else if( charge1 > 0 && charge2 > 0)  fHistPP->Fill(trackVariablesPair,0,1.);
 467       else if( charge1 < 0 && charge2 < 0)  fHistNN->Fill(trackVariablesPair,0,1.);
 468       else {
 469         //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
 470         continue;
 471       }
 472     }//end of 2nd particle loop
 473   }//end of 1st particle loop
 474 }
 475
 476 //____________________________________________________________________//
 477 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
 478                                                  Int_t iVariablePair,
 479                                                  Double_t psiMin,
 480                                                  Double_t psiMax,
 481                                                  Double_t ptTriggerMin,
 482                                                  Double_t ptTriggerMax,
 483                                                  Double_t ptAssociatedMin,
 484                                                  Double_t ptAssociatedMax) {
 485   //Returns the BF histogram, extracted from the 6 AliTHn objects
 486   //(private members) of the AliBalancePsi class.
 487   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 488   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 489   // Psi_2
 490   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 491   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 492   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 493   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 494   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 495   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 496
 497   // pt trigger
 498   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 499     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 500     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 501     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 502     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 503     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 504     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 505   }
 506
 507   // pt associated
 508   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 509     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 510     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 511     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 512     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 513   }
 514
 515   //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));
 516
 517   // Project into the wanted space (1st: analysis step, 2nd: axis)
 518   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
 519   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
 520   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
 521   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
 522   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
 523   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
 524
 525   TH1D *gHistBalanceFunctionHistogram = 0x0;
 526   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 527     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
 528     gHistBalanceFunctionHistogram->Reset();
 529
 530     switch(iVariablePair) {
 531     case 1:
 532       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
 533       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
 534       break;
 535     case 2:
 536       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
 537       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
 538       break;
 539     default:
 540       break;
 541     }
 542
 543     hTemp1->Sumw2();
 544     hTemp2->Sumw2();
 545     hTemp3->Sumw2();
 546     hTemp4->Sumw2();
 547     hTemp1->Add(hTemp3,-1.);
 548     hTemp1->Scale(1./hTemp5->GetEntries());
 549     hTemp2->Add(hTemp4,-1.);
 550     hTemp2->Scale(1./hTemp6->GetEntries());
 551     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 552     gHistBalanceFunctionHistogram->Scale(0.5);
 553   }
 554
 555   return gHistBalanceFunctionHistogram;
 556 }
 557
 558 //____________________________________________________________________//
 559 TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
 560                                                          Int_t iVariablePair,
 561                                                          Double_t psiMin,
 562                                                          Double_t psiMax,
 563                                                          Double_t ptTriggerMin,
 564                                                          Double_t ptTriggerMax,
 565                                                          Double_t ptAssociatedMin,
 566                                                          Double_t ptAssociatedMax,
 567                                                          AliBalancePsi *bfMix) {
 568   //Returns the BF histogram, extracted from the 6 AliTHn objects
 569   //after dividing each correlation function by the Event Mixing one
 570   //(private members) of the AliBalancePsi class.
 571   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 572   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 573   // Psi_2
 574   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 575   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 576   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 577   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 578   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 579   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 580
 581   // pt trigger
 582   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 583     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 584     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 585     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 586     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 587     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 588     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 589   }
 590
 591   // pt associated
 592   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 593     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 594     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 595     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 596     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 597   }
 598
 599   // ============================================================================================
 600   // the same for event mixing
 601   AliTHn *fHistPMix = bfMix->GetHistNp();
 602   AliTHn *fHistNMix = bfMix->GetHistNn();
 603   AliTHn *fHistPNMix = bfMix->GetHistNpn();
 604   AliTHn *fHistNPMix = bfMix->GetHistNnp();
 605   AliTHn *fHistPPMix = bfMix->GetHistNpp();
 606   AliTHn *fHistNNMix = bfMix->GetHistNnn();
 607
 608   // Psi_2
 609   fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 610   fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 611   fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 612   fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 613   fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 614   fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 615
 616   // pt trigger
 617   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 618     fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 619     fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 620     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 621     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 622     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 623     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 624   }
 625
 626   // pt associated
 627   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 628     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 629     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 630     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 631     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 632   }
 633   // ============================================================================================
 634
 635   //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));
 636
 637   // Project into the wanted space (1st: analysis step, 2nd: axis)
 638   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
 639   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
 640   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
 641   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
 642   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
 643   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
 644
 645   // ============================================================================================
 646   // the same for event mixing
 647   TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,iVariablePair);
 648   TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,iVariablePair);
 649   TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,iVariablePair);
 650   TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,iVariablePair);
 651   TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
 652   TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
 653   // ============================================================================================
 654
 655   TH1D *gHistBalanceFunctionHistogram = 0x0;
 656   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 657     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
 658     gHistBalanceFunctionHistogram->Reset();
 659
 660     switch(iVariablePair) {
 661     case 1:
 662       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
 663       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
 664       break;
 665     case 2:
 666       gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
 667       gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
 668       break;
 669     default:
 670       break;
 671     }
 672
 673     hTemp1->Sumw2();
 674     hTemp2->Sumw2();
 675     hTemp3->Sumw2();
 676     hTemp4->Sumw2();
 677     hTemp1Mix->Sumw2();
 678     hTemp2Mix->Sumw2();
 679     hTemp3Mix->Sumw2();
 680     hTemp4Mix->Sumw2();
 681
 682     hTemp1->Scale(1./hTemp5->GetEntries());
 683     hTemp3->Scale(1./hTemp5->GetEntries());
 684     hTemp2->Scale(1./hTemp6->GetEntries());
 685     hTemp4->Scale(1./hTemp6->GetEntries());
 686     hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
 687     hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
 688     hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
 689     hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
 690
 691     hTemp1->Divide(hTemp1Mix);
 692     hTemp2->Divide(hTemp2Mix);
 693     hTemp3->Divide(hTemp3Mix);
 694     hTemp4->Divide(hTemp4Mix);
 695
 696     hTemp1->Add(hTemp3,-1.);
 697     hTemp2->Add(hTemp4,-1.);
 698
 699     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 700     gHistBalanceFunctionHistogram->Scale(0.5);
 701   }
 702
 703   return gHistBalanceFunctionHistogram;
 704 }
 705
 706 //____________________________________________________________________//
 707 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
 708                                                         Double_t psiMax,
 709                                                         Double_t ptTriggerMin,
 710                                                         Double_t ptTriggerMax,
 711                                                         Double_t ptAssociatedMin,
 712                                                         Double_t ptAssociatedMax) {
 713   //Returns the BF histogram in Delta eta vs Delta phi,
 714   //extracted from the 6 AliTHn objects
 715   //(private members) of the AliBalancePsi class.
 716   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 717   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 718   TString histName = "gHistBalanceFunctionHistogram2D";
 719
 720   // Psi_2
 721   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 722   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 723   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 724   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 725   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 726   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 727
 728   // pt trigger
 729   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 730     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 731     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 732     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 733     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 734     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 735     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 736   }
 737
 738   // pt associated
 739   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 740     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 741     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 742     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 743     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 744   }
 745
 746   //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)));
 747
 748   // Project into the wanted space (1st: analysis step, 2nd: axis)
 749   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
 750   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
 751   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
 752   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
 753   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
 754   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
 755
 756   TH2D *gHistBalanceFunctionHistogram = 0x0;
 757   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 758     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
 759     gHistBalanceFunctionHistogram->Reset();
 760     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
 761     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
 762     gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
 763
 764     hTemp1->Sumw2();
 765     hTemp2->Sumw2();
 766     hTemp3->Sumw2();
 767     hTemp4->Sumw2();
 768     hTemp1->Add(hTemp3,-1.);
 769     hTemp1->Scale(1./hTemp5->GetEntries());
 770     hTemp2->Add(hTemp4,-1.);
 771     hTemp2->Scale(1./hTemp6->GetEntries());
 772     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 773     gHistBalanceFunctionHistogram->Scale(0.5);
 774   }
 775
 776   return gHistBalanceFunctionHistogram;
 777 }
 778
 779 //____________________________________________________________________//
 780 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
 781                                                                 Double_t psiMax,
 782                                                                 Double_t ptTriggerMin,
 783                                                                 Double_t ptTriggerMax,
 784                                                                 Double_t ptAssociatedMin,
 785                                                                 Double_t ptAssociatedMax,
 786                                                                 AliBalancePsi *bfMix) {
 787   //Returns the BF histogram in Delta eta vs Delta phi,
 788   //after dividing each correlation function by the Event Mixing one
 789   //extracted from the 6 AliTHn objects
 790   //(private members) of the AliBalancePsi class.
 791   //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
 792   //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
 793   TString histName = "gHistBalanceFunctionHistogram2D";
 794
 795   if(!bfMix){
 796     AliError("balance function object for event mixing not available");
 797     return NULL;
 798   }
 799
 800   // Psi_2
 801   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 802   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 803   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 804   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 805   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 806   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 807
 808   // pt trigger
 809   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 810     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 811     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 812     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 813     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 814     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 815     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 816   }
 817
 818   // pt associated
 819   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 820     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 821     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 822     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 823     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 824   }
 825
 826
 827   // ============================================================================================
 828   // the same for event mixing
 829   AliTHn *fHistPMix = bfMix->GetHistNp();
 830   AliTHn *fHistNMix = bfMix->GetHistNn();
 831   AliTHn *fHistPNMix = bfMix->GetHistNpn();
 832   AliTHn *fHistNPMix = bfMix->GetHistNnp();
 833   AliTHn *fHistPPMix = bfMix->GetHistNpp();
 834   AliTHn *fHistNNMix = bfMix->GetHistNnn();
 835
 836   // Psi_2
 837   fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 838   fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 839   fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 840   fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 841   fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 842   fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 843
 844   // pt trigger
 845   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 846     fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 847     fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 848     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 849     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 850     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 851     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 852   }
 853
 854   // pt associated
 855   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
 856     fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 857     fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 858     fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 859     fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 860   }
 861   // ============================================================================================
 862
 863
 864   //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)));
 865
 866   // Project into the wanted space (1st: analysis step, 2nd: axis)
 867   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
 868   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
 869   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
 870   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
 871   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
 872   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
 873
 874   // ============================================================================================
 875   // the same for event mixing
 876   TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
 877   TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
 878   TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
 879   TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
 880   TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
 881   TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
 882   // ============================================================================================
 883
 884   TH2D *gHistBalanceFunctionHistogram = 0x0;
 885   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
 886     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
 887     gHistBalanceFunctionHistogram->Reset();
 888     gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
 889     gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
 890     gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
 891
 892     hTemp1->Sumw2();
 893     hTemp2->Sumw2();
 894     hTemp3->Sumw2();
 895     hTemp4->Sumw2();
 896     hTemp1Mix->Sumw2();
 897     hTemp2Mix->Sumw2();
 898     hTemp3Mix->Sumw2();
 899     hTemp4Mix->Sumw2();
 900
 901     hTemp1->Scale(1./hTemp5->GetEntries());
 902     hTemp3->Scale(1./hTemp5->GetEntries());
 903     hTemp2->Scale(1./hTemp6->GetEntries());
 904     hTemp4->Scale(1./hTemp6->GetEntries());
 905     hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
 906     hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
 907     hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
 908     hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
 909
 910     hTemp1->Divide(hTemp1Mix);
 911     hTemp2->Divide(hTemp2Mix);
 912     hTemp3->Divide(hTemp3Mix);
 913     hTemp4->Divide(hTemp4Mix);
 914
 915     hTemp1->Add(hTemp3,-1.);
 916     hTemp2->Add(hTemp4,-1.);
 917
 918     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
 919     gHistBalanceFunctionHistogram->Scale(0.5);
 920   }
 921
 922   return gHistBalanceFunctionHistogram;
 923 }
 924
 925
 926 //____________________________________________________________________//
 927 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
 928                                               Double_t psiMax,
 929                                               Double_t ptTriggerMin,
 930                                               Double_t ptTriggerMax,
 931                                               Double_t ptAssociatedMin,
 932                                               Double_t ptAssociatedMax,
 933                                               Bool_t   normToTrig) {
 934   //Returns the 2D correlation function for (+-) pairs
 935   // Psi_2: axis 0
 936   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 937   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 938
 939   // pt trigger
 940   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 941     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 942     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
 943   }
 944
 945   // pt associated
 946   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
 947     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
 948
 949   //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
 950   //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
 951
 952   //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
 953   //TCanvas *c1 = new TCanvas("c1","");
 954   //c1->cd();
 955   //if(!gHistTest){
 956   //AliError("Projection of fHistP = NULL");
 957   //return gHistTest;
 958   //}
 959   //else{
 960   //gHistTest->DrawCopy("colz");
 961   //}
 962
 963   //0:step, 1: Delta eta, 2: Delta phi
 964   TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
 965   if(!gHist){
 966     AliError("Projection of fHistPN = NULL");
 967     return gHist;
 968   }
 969
 970   //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
 971   //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
 972   //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
 973
 974   //TCanvas *c2 = new TCanvas("c2","");
 975   //c2->cd();
 976   //fHistPN->Project(0,1,2)->DrawCopy("colz");
 977
 978   if(normToTrig && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
 979     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
 980
 981   return gHist;
 982 }
 983
 984 //____________________________________________________________________//
 985 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
 986                                               Double_t psiMax,
 987                                               Double_t ptTriggerMin,
 988                                               Double_t ptTriggerMax,
 989                                               Double_t ptAssociatedMin,
 990                                               Double_t ptAssociatedMax,
 991                                               Bool_t   normToTrig) {
 992   //Returns the 2D correlation function for (+-) pairs
 993   // Psi_2: axis 0
 994   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 995   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
 996
 997   // pt trigger
 998   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
 999     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1000     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1001   }
1002
1003   // pt associated
1004   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1005     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1006
1007   //0:step, 1: Delta eta, 2: Delta phi
1008   TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
1009   if(!gHist){
1010     AliError("Projection of fHistPN = NULL");
1011     return gHist;
1012   }
1013
1014   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1015   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
1016   if(normToTrig && (Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
1017     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
1018
1019   return gHist;
1020 }
1021
1022 //____________________________________________________________________//
1023 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
1024                                               Double_t psiMax,
1025                                               Double_t ptTriggerMin,
1026                                               Double_t ptTriggerMax,
1027                                               Double_t ptAssociatedMin,
1028                                               Double_t ptAssociatedMax,
1029                                               Bool_t   normToTrig) {
1030   //Returns the 2D correlation function for (+-) pairs
1031   // Psi_2: axis 0
1032   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1033   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1034
1035   // pt trigger
1036   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1037     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1038     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1039   }
1040
1041   // pt associated
1042   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1043     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1044
1045   //0:step, 1: Delta eta, 2: Delta phi
1046   TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
1047   if(!gHist){
1048     AliError("Projection of fHistPN = NULL");
1049     return gHist;
1050   }
1051
1052   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
1053   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
1054   if(normToTrig && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
1055     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
1056
1057   return gHist;
1058 }
1059
1060 //____________________________________________________________________//
1061 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
1062                                               Double_t psiMax,
1063                                               Double_t ptTriggerMin,
1064                                               Double_t ptTriggerMax,
1065                                               Double_t ptAssociatedMin,
1066                                               Double_t ptAssociatedMax,
1067                                               Bool_t   normToTrig) {
1068   //Returns the 2D correlation function for (+-) pairs
1069   // Psi_2: axis 0
1070   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1071   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1072
1073   // pt trigger
1074   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1075     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1076     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1077   }
1078
1079   // pt associated
1080   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1081     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1082
1083   //0:step, 1: Delta eta, 2: Delta phi
1084   TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
1085   if(!gHist){
1086     AliError("Projection of fHistPN = NULL");
1087     return gHist;
1088   }
1089
1090   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1091   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
1092   if(normToTrig && (Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
1093     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
1094
1095   return gHist;
1096 }
1097
1098 //____________________________________________________________________//
1099 TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
1100                                                              Double_t psiMax,
1101                                                              Double_t ptTriggerMin,
1102                                                              Double_t ptTriggerMax,
1103                                                              Double_t ptAssociatedMin,
1104                                                              Double_t ptAssociatedMax,
1105                                                              Bool_t   normToTrig) {
1106   //Returns the 2D correlation function for the sum of all charge combination pairs
1107   // Psi_2: axis 0
1108   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1109   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1110   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1111   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1112   fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1113   fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1114
1115   // pt trigger
1116   if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1117     fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1118     fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1119     fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1120     fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1121     fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1122     fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1123   }
1124
1125   // pt associated
1126   if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1127     fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1128     fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1129     fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1130     fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
1131
1132   //0:step, 1: Delta eta, 2: Delta phi
1133   TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
1134   if(!gHistNN){
1135     AliError("Projection of fHistNN = NULL");
1136     return gHistNN;
1137   }
1138   TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
1139   if(!gHistPP){
1140     AliError("Projection of fHistPP = NULL");
1141     return gHistPP;
1142   }
1143   TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
1144   if(!gHistNP){
1145     AliError("Projection of fHistNP = NULL");
1146     return gHistNP;
1147   }
1148   TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
1149   if(!gHistPN){
1150     AliError("Projection of fHistPN = NULL");
1151     return gHistPN;
1152   }
1153
1154   // sum all 2 particle histograms
1155   gHistNN->Add(gHistPP);
1156   gHistNN->Add(gHistNP);
1157   gHistNN->Add(gHistPN);
1158
1159   // if normalization to trigger particle, divide by sum of + and - triggers
1160   if(normToTrig && (Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
1161     gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
1162
1163   return gHistNN;
1164 }
1165
1166 //____________________________________________________________________//
1167 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) {
1168   //
1169   // calculates dphistar
1170   //
1171   Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1172
1173   static const Double_t kPi = TMath::Pi();
1174
1175   // circularity
1176 //   if (dphistar > 2 * kPi)
1177 //     dphistar -= 2 * kPi;
1178 //   if (dphistar < -2 * kPi)
1179 //     dphistar += 2 * kPi;
1180
1181   if (dphistar > kPi)
1182     dphistar = kPi * 2 - dphistar;
1183   if (dphistar < -kPi)
1184     dphistar = -kPi * 2 - dphistar;
1185   if (dphistar > kPi) // might look funny but is needed
1186     dphistar = kPi * 2 - dphistar;
1187
1188   return dphistar;
1189 }
1190