Adding the pt trigger and associated bins in the Psi dependent bf analysis
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 20 May 2012 11:30:52 +0000 (11:30 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 20 May 2012 11:30:52 +0000 (11:30 +0000)
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h
PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C
PWGCF/EBYE/macros/drawBalanceFunctionPsi.C [new file with mode: 0644]
PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C [new file with mode: 0644]

index 9f744e1..1e6a3c8 100644 (file)
@@ -129,10 +129,10 @@ void AliBalancePsi::InitHistograms() {
   axisTitlePair[1]  = "#phi - #Psi_{2} (#circ)"; 
   
   // eta
-  const Int_t kNEtaBins = 20;
+  const Int_t kNEtaBins = 10;
   Double_t etaBins[kNEtaBins+1];
   for(Int_t i = 0; i < kNEtaBins+1; i++)
-    etaBins[i] = -1.0 + i * 0.1;
+    etaBins[i] = -1.0 + i * 0.2;
   iBinSingle[2]       = kNEtaBins;
   dBinsSingle[2]      = etaBins;
   axisTitleSingle[2]  = "#eta"; 
@@ -143,10 +143,10 @@ void AliBalancePsi::InitHistograms() {
   axisTitlePair[2]  = "#eta"; */
 
   // delta eta
-  const Int_t kNDeltaEtaBins = 20;
+  const Int_t kNDeltaEtaBins = 80;
   Double_t deltaEtaBins[kNDeltaEtaBins+1];
   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
-    deltaEtaBins[i] = -2.0 + i * 0.2;
+    deltaEtaBins[i] = -2.0 + i * 0.05;
   iBinPair[2]       = kNDeltaEtaBins;
   dBinsPair[2]      = deltaEtaBins;
   axisTitlePair[2]  = "#Delta #eta"; 
@@ -175,37 +175,15 @@ void AliBalancePsi::InitHistograms() {
   axisTitlePair[8]  = "#Delta y"; */
 
   // phi
-  const Int_t kNPhiBins = 36;
+  const Int_t kNPhiBins = 18;
   Double_t phiBins[kNPhiBins+1];
   for(Int_t i = 0; i < kNPhiBins+1; i++){
-    phiBins[i] = 0.0 + i * 10.;
+    phiBins[i] = 0.0 + i * 20.;
   } 
   iBinSingle[3]       = kNPhiBins;
   dBinsSingle[3]      = phiBins;
   axisTitleSingle[3]  = "#phi (#circ)"; 
-
-  /*iBinPair[4]       = kNPhiBins;
-  dBinsPair[4]      = phiBins;
-  axisTitlePair[4]  = "#phi (#circ)"; */
   
-  // pt(trigger)
-  /*const Int_t kNPtBins = 100;
-  Double_t ptBins[kNPtBins+1];
-  for(Int_t i = 0; i < kNPtBins+1; i++){
-    ptBins[i] = 0.0 + i * 0.2;
-   } 
-  iBinSingle[4]       = kNPtBins;
-  dBinsSingle[4]      = ptBins;
-  axisTitleSingle[4]  = "p_{t}^{trig.} (GeV/c)"; */
-
-  /*iBinPair[5]       = kNPtBins;
-  dBinsPair[5]      = ptBins;
-  axisTitlePair[5]  = "p_{t}^{trig.} (GeV/c)"; 
-
-  iBinPair[6]       = kNPtBins;
-  dBinsPair[6]      = ptBins;
-  axisTitlePair[6]  = "p_{t}^{assoc.} (GeV/c)"; */
-
   // delta phi
   const Int_t kNDeltaPhiBins = 72;
   Double_t deltaPhiBins[kNDeltaPhiBins+1];
@@ -216,6 +194,24 @@ void AliBalancePsi::InitHistograms() {
   dBinsPair[3]      = deltaPhiBins;
   axisTitlePair[3]  = "#Delta #phi (#circ)"; 
 
+  // pt(trigger-associated)
+  const Int_t kNPtBins = 50;
+  Double_t ptBins[kNPtBins+1];
+  for(Int_t i = 0; i < kNPtBins+1; i++){
+    ptBins[i] = 0.0 + i * 0.5;
+   } 
+  iBinSingle[4]       = kNPtBins;
+  dBinsSingle[4]      = ptBins;
+  axisTitleSingle[4]  = "p_{t}^{trig.} (GeV/c)"; 
+
+  iBinPair[4]       = kNPtBins;
+  dBinsPair[4]      = ptBins;
+  axisTitlePair[4]  = "p_{t}^{trig.} (GeV/c)"; 
+
+  iBinPair[5]       = kNPtBins;
+  dBinsPair[5]      = ptBins;
+  axisTitlePair[5]  = "p_{t}^{assoc.} (GeV/c)";   
+
   // qside
   /*const Int_t kNQSideBins = 200;
   Double_t qSideBins[kNQSideBins+1];
@@ -373,7 +369,7 @@ void AliBalancePsi::CalculateBalance(Float_t fCentrality,
     trackVarsSingle[2]    =  eta1;
     //trackVarsSingle[3]    =  rap1;
     trackVarsSingle[3]    =  phi1;  
-    //trackVarsSingle[5]    =  pt1;  
+    trackVarsSingle[4]    =  pt1;  
 
     //Printf("Track(a) %d - phi-Psi: %lf",i+1,trackVarsSingle[1]);
     //fill single particle histograms
@@ -443,11 +439,11 @@ void AliBalancePsi::CalculateBalance(Float_t fCentrality,
       //trackVarsPair[2]    =  eta1;
       //trackVarsPair[3]    =  rap1;
       //trackVarsPair[4]    =  phi1;
-      //trackVarsPair[5]    =  pt1;
-      //trackVarsPair[6]    =  pt2;
       trackVarsPair[2]    =  deta;
       //trackVarsPair[8]    =  dy;
       trackVarsPair[3]    =  dphi;
+      trackVarsPair[4]    =  pt1;
+      trackVarsPair[5]    =  pt2;
       //trackVarsPair[10]   =  qSide;
       //trackVarsPair[11]   =  qOut;
       //trackVarsPair[12]   =  qLong;
@@ -494,6 +490,8 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
   fHistPP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
   fHistNN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
 
+  //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));
+
   // Project into the wanted space (1st: analysis step, 2nd: axis)
   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
@@ -544,9 +542,9 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
     hTemp2->Sumw2();
     hTemp3->Sumw2();
     hTemp4->Sumw2();
-    hTemp1->Add(hTemp3,-2.);
+    hTemp1->Add(hTemp3,-1.);
     hTemp1->Scale(1./hTemp5->GetEntries());
-    hTemp2->Add(hTemp4,-2.);
+    hTemp2->Add(hTemp4,-1.);
     hTemp2->Scale(1./hTemp6->GetEntries());
     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
     gHistBalanceFunctionHistogram->Scale(0.5);
@@ -554,3 +552,99 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
 
   return gHistBalanceFunctionHistogram;
 }
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t centrMin, 
+                                             Double_t centrMax, 
+                                             Double_t psiMin, 
+                                             Double_t psiMax) {
+  //Returns the 2D correlation function for (+-) pairs
+  // centrality: axis 0
+  fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
+  fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
+  
+  // Psi_2: axis 1
+  fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
+  fHistPN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
+    
+  //0:step, 2: Delta eta, 3: Delta phi
+  TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,2,3));
+  //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
+  //Printf("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,2,3)->GetEntries()));
+  if((Double_t)(fHistP->Project(0,2)->GetEntries())!=0)
+    gHist->Scale(1./(Double_t)(fHistP->Project(0,2)->GetEntries()));
+    
+  return gHist;
+}
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t centrMin, 
+                                             Double_t centrMax, 
+                                             Double_t psiMin, 
+                                             Double_t psiMax) {
+  //Returns the 2D correlation function for (+-) pairs
+  // centrality: axis 0
+  fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
+  fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
+  
+  // Psi_2: axis 1
+  fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
+  fHistNP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
+    
+  //0:step, 2: Delta eta, 3: Delta phi
+  TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,2,3));
+  //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
+  //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
+  if((Double_t)(fHistN->Project(0,2)->GetEntries())!=0)
+    gHist->Scale(1./(Double_t)(fHistN->Project(0,2)->GetEntries()));
+    
+  return gHist;
+}
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t centrMin, 
+                                             Double_t centrMax, 
+                                             Double_t psiMin, 
+                                             Double_t psiMax) {
+  //Returns the 2D correlation function for (+-) pairs
+  // centrality: axis 0
+  fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
+  fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
+  
+  // Psi_2: axis 1
+  fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
+  fHistPP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
+      
+  //0:step, 2: Delta eta, 3: Delta phi
+  TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,2,3));
+  //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
+  //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
+  if((Double_t)(fHistP->Project(0,2)->GetEntries())!=0)
+    gHist->Scale(1./(Double_t)(fHistP->Project(0,2)->GetEntries()));
+  
+  return gHist;
+}
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t centrMin, 
+                                             Double_t centrMax, 
+                                             Double_t psiMin, 
+                                             Double_t psiMax) {
+  //Returns the 2D correlation function for (+-) pairs
+  // centrality: axis 0
+  fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
+  fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); 
+  
+  // Psi_2: axis 1
+  fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
+  fHistNN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); 
+    
+  //0:step, 2: Delta eta, 3: Delta phi
+  TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,2,3));
+  //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
+  //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
+  if((Double_t)(fHistN->Project(0,2)->GetEntries())!=0)
+    gHist->Scale(1./(Double_t)(fHistN->Project(0,2)->GetEntries()));
+    
+  return gHist;
+}
index 10580c2..4a09c3b 100644 (file)
 #include <TObject.h>
 #include "TString.h"
 
+#include "AliTHn.h"
+
 #define ANALYSIS_TYPES 7
 #define MAXIMUM_NUMBER_OF_STEPS        1024
 #define MAXIMUM_STEPS_IN_PSI 360
 
-class TGraphErrors;
 class TH1D;
 class TH2D;
 class TH3D;
-class AliTHn;
 
-const Int_t nTrackVariablesSingle = 4;       // track variables in histogram (eta, phi, pTtrig, centrality)
-const Int_t nTrackVariablesPair   = 4;       // track variables in histogram (dEta, dPhi, pT, pTtrig, centrality)
+const Int_t nTrackVariablesSingle = 5;       // track variables in histogram (centrality, phi-Psi2, eta, phi, pTtrig)
+const Int_t nTrackVariablesPair   = 6;       // track variables in histogram (centrality, phi-Psi2, dEta, dPhi, pTtrig, ptAssociated)
 const TString gBFPsiAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
 
 class AliBalancePsi : public TObject {
@@ -60,6 +60,15 @@ class AliBalancePsi : public TObject {
 
   void CalculateBalance(Float_t fCentrality, Double_t gReactionPlane, vector<Double_t> **chargeVector);
   
+  TH2D   *GetCorrelationFunctionPN(Double_t centrMin, Double_t centrMax, 
+                                  Double_t psiMin, Double_t psiMax);
+  TH2D   *GetCorrelationFunctionNP(Double_t centrMin, Double_t centrMax, 
+                                  Double_t psiMin, Double_t psiMax);
+  TH2D   *GetCorrelationFunctionPP(Double_t centrMin, Double_t centrMax, 
+                                  Double_t psiMin, Double_t psiMax);
+  TH2D   *GetCorrelationFunctionNN(Double_t centrMin, Double_t centrMax, 
+                                  Double_t psiMin, Double_t psiMax);
+
   AliTHn *GetHistNp() {return fHistP;}
   AliTHn *GetHistNn() {return fHistN;}
   AliTHn *GetHistNpn() {return fHistPN;}
@@ -67,12 +76,18 @@ class AliBalancePsi : public TObject {
   AliTHn *GetHistNpp() {return fHistPP;}
   AliTHn *GetHistNnn() {return fHistNN;}
 
-  void SetHistNp(AliTHn *gHist) {fHistP = gHist;}
-  void SetHistNn(AliTHn *gHist) {fHistN = gHist;}
-  void SetHistNpn(AliTHn *gHist) {fHistPN = gHist;}
-  void SetHistNnp(AliTHn *gHist) {fHistNP = gHist;}
-  void SetHistNpp(AliTHn *gHist) {fHistPP = gHist;}
-  void SetHistNnn(AliTHn *gHist) {fHistNN = gHist;}
+  void SetHistNp(AliTHn *gHist) {
+    fHistP = gHist; }//fHistP->FillParent(); fHistP->DeleteContainers();}
+  void SetHistNn(AliTHn *gHist) {
+    fHistN = gHist; }//fHistN->FillParent(); fHistN->DeleteContainers();}
+  void SetHistNpn(AliTHn *gHist) {
+    fHistPN = gHist; }//fHistPN->FillParent(); fHistPN->DeleteContainers();}
+  void SetHistNnp(AliTHn *gHist) {
+    fHistNP = gHist; }//fHistNP->FillParent(); fHistNP->DeleteContainers();}
+  void SetHistNpp(AliTHn *gHist) {
+    fHistPP = gHist; }//fHistPP->FillParent(); fHistPP->DeleteContainers();}
+  void SetHistNnn(AliTHn *gHist) {
+    fHistNN = gHist; }//fHistNN->FillParent(); fHistNN->DeleteContainers();}
 
   TH1D *GetBalanceFunctionHistogram(Int_t iVariableSingle,
                                    Int_t iVariablePair,
index 03b5530..2b5ce5b 100644 (file)
@@ -34,38 +34,6 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,
                                                       TString fileNameBase="AnalysisResults") {\r
   // Creates a balance function analysis task and adds it to the analysis manager.\r
   // Get the pointer to the existing analysis manager via the static access method.\r
-  TString centralityName("");\r
-  centralityName+=Form("%.0f",centrMin);\r
-  centralityName+="-";\r
-  centralityName+=Form("%.0f",centrMax);\r
-  centralityName+="_";\r
-  centralityName+=Form("%s",centralityEstimator.Data());\r
-  centralityName+="_";\r
-  centralityName+=Form("vZ%.1f",vertexZ);\r
-  centralityName+="_";\r
-  centralityName+=Form("DCAxy%.1f",DCAxy);\r
-  centralityName+="_";\r
-  centralityName+=Form("DCAz%.1f",DCAz);\r
-  centralityName+="_Pt";\r
-  centralityName+=Form("%.1f",ptMin);\r
-  centralityName+="-";\r
-  centralityName+=Form("%.1f",ptMax);\r
-  centralityName+="_Eta";\r
-  centralityName+=Form("%.1f",etaMin);\r
-  centralityName+="-";\r
-  centralityName+=Form("%.1f",etaMax);\r
-  centralityName+="_Chi";\r
-  centralityName+=Form("%.1f",maxTPCchi2);\r
-  centralityName+="_nClus";\r
-  centralityName+=Form("%d",minNClustersTPC);\r
-  centralityName+="_Bit";\r
-  centralityName+=Form("%d",AODfilterBit);\r
-  if(bCentralTrigger)   centralityName+="_withCentralTrigger";\r
-\r
-\r
-\r
-\r
-\r
   TString outputFileName(fileNameBase);\r
   outputFileName.Append(".root");\r
 \r
@@ -165,11 +133,11 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,
   // Get and connect other common input/output containers via the manager as below\r
   //==============================================================================\r
   TString outputFileName = AliAnalysisManager::GetCommonFileName();\r
-  outputFileName += ":PWGCFEbyE.outputBalanceFunctionAnalysis";\r
-  AliAnalysisDataContainer *coutQA = mgr->CreateContainer(Form("listQA_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
-  AliAnalysisDataContainer *coutBF = mgr->CreateContainer(Form("listBF_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
-  if(gRunShuffling) AliAnalysisDataContainer *coutBFS = mgr->CreateContainer(Form("listBFShuffled_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
-  if(kUsePID) AliAnalysisDataContainer *coutQAPID = mgr->CreateContainer(Form("listQAPID_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  outputFileName += ":PWGCFEbyE.outputBalanceFunctionPsiAnalysis";\r
+  AliAnalysisDataContainer *coutQA = mgr->CreateContainer("listQAPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  AliAnalysisDataContainer *coutBF = mgr->CreateContainer("listBFPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(gRunShuffling) AliAnalysisDataContainer *coutBFS = mgr->CreateContainer("listBFPsiShuffled", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(kUsePID) AliAnalysisDataContainer *coutQAPID = mgr->CreateContainer("listQAPIDPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
 \r
   mgr->ConnectInput(taskBF, 0, mgr->GetCommonInputContainer());\r
   mgr->ConnectOutput(taskBF, 1, coutQA);\r
diff --git a/PWGCF/EBYE/macros/drawBalanceFunctionPsi.C b/PWGCF/EBYE/macros/drawBalanceFunctionPsi.C
new file mode 100644 (file)
index 0000000..00c0bc4
--- /dev/null
@@ -0,0 +1,323 @@
+const Int_t numberOfCentralityBins = 9;
+TString centralityArray[numberOfCentralityBins] = {"0-5","5-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80"};
+Double_t gMinCentrality[numberOfCentralityBins] = {0.,5.,10.,20.,30.,40.,50.,60.,70.};
+Double_t gMaxCentrality[numberOfCentralityBins] = {5.,10.,20.,30.,40.,50.,60.,70.,80.};
+TString gAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
+
+const Int_t gRebin = 1;
+void drawBalanceFunctionPsi(const char* filename = "AnalysisResultsPsi.root", 
+                           Double_t psiMin = 0., Double_t psiMax = 7.5) {
+  //Macro that draws the BF distributions for each centrality bin
+  //for reaction plane dependent analysis
+  //Author: Panos.Christakoglou@nikhef.nl
+  //Load the PWG2ebye library
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libEventMixing.so");
+  gSystem->Load("libCORRFW.so");
+  gSystem->Load("libPWGTools.so");
+  gSystem->Load("libPWGCFebye.so");
+
+  //Prepare the objects and return them
+  TList *listBF = GetListOfObjects(filename);
+  TList *listBFShuffled = GetListOfObjects(filename,kTRUE);
+  if(!listBF) {
+    Printf("The TList object was not created");
+    return;
+  }
+  else 
+    draw(listBF,listBFShuffled,psiMin,psiMax);  
+}
+
+//______________________________________________________//
+TList *GetListOfObjects(const char* filename, Bool_t kShuffling = kFALSE) {
+  //Get the TList objects (QA, bf, bf shuffled)
+  TList *listQA = 0x0;
+  TList *listBF = 0x0;
+  TList *listBFShuffling = 0x0;
+  
+  //Open the file
+  TFile *f = TFile::Open(filename);
+  if((!f)||(!f->IsOpen())) {
+    Printf("The file %s is not found. Aborting...",filename);
+    return listBF;
+  }
+  //f->ls();
+  
+  TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionAnalysis"));
+  if(!dir) {   
+    Printf("The TDirectoryFile is not found. Aborting...",filename);
+    return listBF;
+  }
+  //dir->ls();
+  
+  TString listBFName;
+  if(!kShuffling) 
+    listBFName = "listBF_0-100_V0M_vZ10.0_DCAxy-1.0_DCAz-1.0_Pt0.3-5.0_Eta-0.8-0.8_Chi-1.0_nClus-1_Bit1_withCentralTrigger";
+  else if(kShuffling)
+    listBFName = "listBFShuffled_0-100_V0M_vZ10.0_DCAxy-1.0_DCAz-1.0_Pt0.3-5.0_Eta-0.8-0.8_Chi-1.0_nClus-1_Bit1_withCentralTrigger";
+  listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
+  //listBF->ls();
+
+  //Get the histograms
+  TString histoName;
+  if(!kShuffling)
+    histoName = "fHistPV0M";
+  else if(kShuffling)
+    histoName = "fHistP_shuffleV0M";
+  AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));  
+  if(!fHistP) {
+    Printf("fHistP %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistP->FillParent(); fHistP->DeleteContainers();
+
+  if(!kShuffling)
+    histoName = "fHistNV0M";
+  else if(kShuffling)
+    histoName = "fHistN_shuffleV0M";
+  AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistN) {
+    Printf("fHistN %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistN->FillParent(); fHistN->DeleteContainers();
+    
+  if(!kShuffling)
+    histoName = "fHistPNV0M";
+  else if(kShuffling)
+    histoName = "fHistPN_shuffleV0M";
+  AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistPN) {
+    Printf("fHistPN %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistPN->FillParent(); fHistPN->DeleteContainers();
+  
+  if(!kShuffling)
+    histoName = "fHistNPV0M";
+  else if(kShuffling)
+    histoName = "fHistNP_shuffleV0M";
+  AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistNP) {
+    Printf("fHistNP %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistNP->FillParent(); fHistNP->DeleteContainers();
+
+  if(!kShuffling)
+    histoName = "fHistPPV0M";
+  else if(kShuffling)
+    histoName = "fHistPP_shuffleV0M";
+  AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistPP) {
+    Printf("fHistPP %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistPP->FillParent(); fHistPP->DeleteContainers();
+
+  if(!kShuffling)
+    histoName = "fHistNNV0M";
+  else if(kShuffling)
+    histoName = "fHistNN_shuffleV0M";
+  AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistNN) {
+    Printf("fHistNN %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistNN->FillParent(); fHistNN->DeleteContainers();
+  
+  return listBF;
+}
+
+//______________________________________________________//
+void draw(TList *listBF, TList *listBFShuffled,
+         Double_t psiMin, Double_t psiMax) {
+  gROOT->LoadMacro("~/SetPlotStyle.C");
+  SetPlotStyle();
+  gStyle->SetPalette(1,0);
+  
+  //balance function
+  AliTHn *hP = NULL;
+  AliTHn *hN = NULL;
+  AliTHn *hPN = NULL;
+  AliTHn *hNP = NULL;
+  AliTHn *hPP = NULL;
+  AliTHn *hNN = NULL;
+  //listBF->ls();
+  //Printf("=================");
+  hP = (AliTHn*) listBF->FindObject("fHistPV0M");
+  hN = (AliTHn*) listBF->FindObject("fHistNV0M");
+  hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
+  hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
+  hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
+  hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
+
+  AliBalancePsi *b = new AliBalancePsi();
+  b->SetHistNp(hP);
+  b->SetHistNn(hN);
+  b->SetHistNpn(hPN);
+  b->SetHistNnp(hNP);
+  b->SetHistNpp(hPP);
+  b->SetHistNnn(hNN);
+
+  //balance function shuffling
+  AliTHn *hPShuffled = NULL;
+  AliTHn *hNShuffled = NULL;
+  AliTHn *hPNShuffled = NULL;
+  AliTHn *hNPShuffled = NULL;
+  AliTHn *hPPShuffled = NULL;
+  AliTHn *hNNShuffled = NULL;
+  //listBFShuffled->ls();
+  hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
+  hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
+  hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
+  hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
+  hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
+  hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
+
+  AliBalancePsi *bShuffled = new AliBalancePsi();
+  bShuffled->SetHistNp(hPShuffled);
+  bShuffled->SetHistNn(hNShuffled);
+  bShuffled->SetHistNpn(hPNShuffled);
+  bShuffled->SetHistNnp(hNPShuffled);
+  bShuffled->SetHistNpp(hPPShuffled);
+  bShuffled->SetHistNnn(hNNShuffled);
+
+  TH1D *gHistBalanceFunction[numberOfCentralityBins];
+  TH1D *gHistBalanceFunctionShuffled[numberOfCentralityBins];
+  TCanvas *c1[numberOfCentralityBins];
+  TString histoTitle, pngName;
+  TLegend *legend[numberOfCentralityBins];
+  
+  //loop over the centrality bins
+  for(Int_t iCentralityBin = 0; iCentralityBin < numberOfCentralityBins; iCentralityBin++) {  
+    histoTitle = "Centrality: "; 
+    histoTitle += centralityArray[iCentralityBin]; 
+    histoTitle += "%";
+    histoTitle += " | "; histoTitle += psiMin; 
+    histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax;
+
+    gHistBalanceFunction[iCentralityBin] = b->GetBalanceFunctionHistogram(2,2,gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax);
+    //gHistBalanceFunction[iCentralityBin] = b->GetBalanceFunctionHistogram(3,3,gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax);
+    gHistBalanceFunction[iCentralityBin]->SetMarkerStyle(20);
+    gHistBalanceFunction[iCentralityBin]->SetTitle(histoTitle.Data());
+    gHistBalanceFunction[iCentralityBin]->GetYaxis()->SetTitleOffset(1.3);
+
+    gHistBalanceFunctionShuffled[iCentralityBin] = bShuffled->GetBalanceFunctionHistogram(2,2,gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax);
+    //gHistBalanceFunctionShuffled[iCentralityBin] = b->GetBalanceFunctionHistogram(3,3,gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax);
+    gHistBalanceFunctionShuffled[iCentralityBin]->SetMarkerStyle(24);
+
+    c1[iCentralityBin] = new TCanvas(histoTitle.Data(),"",0,0,600,500);
+    c1[iCentralityBin]->SetFillColor(10); 
+    c1[iCentralityBin]->SetHighLightColor(10);
+    c1[iCentralityBin]->SetLeftMargin(0.15);
+    gHistBalanceFunction[iCentralityBin]->Draw("E");
+    gHistBalanceFunctionShuffled[iCentralityBin]->Draw("ESAME");
+
+    legend[iCentralityBin] = new TLegend(0.18,0.6,0.45,0.82,"","brNDC");
+    legend[iCentralityBin]->SetTextSize(0.045); 
+    legend[iCentralityBin]->SetTextFont(42); 
+    legend[iCentralityBin]->SetBorderSize(0);
+    legend[iCentralityBin]->SetFillStyle(0); 
+    legend[iCentralityBin]->SetFillColor(10);
+    legend[iCentralityBin]->SetMargin(0.25); 
+    legend[iCentralityBin]->SetShadowColor(10);
+    legend[iCentralityBin]->AddEntry(gHistBalanceFunction[iCentralityBin],"Data","lp");
+    legend[iCentralityBin]->AddEntry(gHistBalanceFunctionShuffled[iCentralityBin],"Shuffled data","lp");
+    legend[iCentralityBin]->Draw();
+
+    pngName = "BalanceFunctionDeltaEta.Centrality"; 
+    pngName += centralityArray[iCentralityBin]; 
+    pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
+    pngName += ".png";
+    c1[iCentralityBin]->SaveAs(pngName.Data());
+
+    GetWeightedMean(gHistBalanceFunction[iCentralityBin]);
+    GetWeightedMean(gHistBalanceFunctionShuffled[iCentralityBin]);
+
+    TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
+    meanLatex = "#mu = "; 
+    meanLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetMean());
+    meanLatex += " #pm "; 
+    meanLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetMeanError());
+
+    rmsLatex = "#sigma = "; 
+    rmsLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetRMS());
+    rmsLatex += " #pm "; 
+    rmsLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetRMSError());
+
+    skewnessLatex = "S = "; 
+    skewnessLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetSkewness(1));
+    skewnessLatex += " #pm "; 
+    skewnessLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetSkewness(11));
+    kurtosisLatex = "K = "; 
+    kurtosisLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetKurtosis(1));
+    kurtosisLatex += " #pm "; 
+    kurtosisLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetKurtosis(11));
+
+    TLatex *latex = new TLatex();
+    latex->SetNDC();
+    latex->SetTextSize(0.035);
+    latex->SetTextColor(1);
+    latex->DrawLatex(0.64,0.85,meanLatex.Data());
+    latex->DrawLatex(0.64,0.81,rmsLatex.Data());
+    latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
+    latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
+    Printf("Mean: %lf - Error: %lf",gHistBalanceFunction[iCentralityBin]->GetMean(),gHistBalanceFunction[iCentralityBin]->GetMeanError());
+    Printf("RMS: %lf - Error: %lf",gHistBalanceFunction[iCentralityBin]->GetRMS(),gHistBalanceFunction[iCentralityBin]->GetRMSError());
+    Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunction[iCentralityBin]->GetSkewness(1),gHistBalanceFunction[iCentralityBin]->GetSkewness(11));
+    Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunction[iCentralityBin]->GetKurtosis(1),gHistBalanceFunction[iCentralityBin]->GetKurtosis(11));
+  }
+}
+
+//____________________________________________________________________//
+void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1) {
+  //Prints the calculated width of the BF and its error
+  Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
+  Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
+  Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
+  Double_t deltaBalP2 = 0.0, integral = 0.0;
+  Double_t deltaErrorNew = 0.0;
+
+  //Retrieve this variables from Histogram
+  Int_t fNumberOfBins = gHistBalance->GetNbinsX();
+  Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
+  
+  cout<<"=================================================="<<endl;
+  cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
+  cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
+  cout<<"=================================================="<<endl;
+  for(Int_t i = 1; i <= fNumberOfBins; i++) {
+
+    // this is to simulate |Delta eta| or |Delta phi|
+    if(fNumberOfBins/2 - fStartBin + 1 < i && i < fNumberOfBins/2 + fStartBin ) continue;
+
+    cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<TMath::Abs(gHistBalance->GetBinCenter(i))<<endl;
+
+    gSumXi += TMath::Abs(gHistBalance->GetBinCenter(i)); // this is to simulate |Delta eta| or |Delta phi|
+    gSumBi += gHistBalance->GetBinContent(i); 
+    gSumBiXi += gHistBalance->GetBinContent(i)*TMath::Abs(gHistBalance->GetBinCenter(i));
+    gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
+    gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
+    gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
+    gSumXi2DeltaBi2 += TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2) * TMath::Power(gHistBalance->GetBinError(i),2);
+    
+    deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
+    integral += fP2Step*gHistBalance->GetBinContent(i);
+  }
+  for(Int_t i = fStartBin; i < fNumberOfBins; i++)
+    deltaErrorNew += gHistBalance->GetBinError(i)*(TMath::Abs(gHistBalance->GetBinCenter(i))*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
+  
+  Double_t integralError = TMath::Sqrt(deltaBalP2);
+  
+  Double_t delta = gSumBiXi / gSumBi;
+  Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
+  cout<<"=================================================="<<endl;
+  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
+  cout<<"New error: "<<deltaErrorNew<<endl;
+  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
+  cout<<"=================================================="<<endl;
+  cout<<endl;
+}
diff --git a/PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C b/PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C
new file mode 100644 (file)
index 0000000..a61b3ad
--- /dev/null
@@ -0,0 +1,224 @@
+const Int_t numberOfCentralityBins = 9;
+TString centralityArray[numberOfCentralityBins] = {"0-5","5-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80"};
+Double_t gMinCentrality[numberOfCentralityBins] = {0.,5.,10.,20.,30.,40.,50.,60.,70.};
+Double_t gMaxCentrality[numberOfCentralityBins] = {5.,10.,20.,30.,40.,50.,60.,70.,80.};
+TString gAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
+
+const Int_t gRebin = 1;
+void drawCorrelationFunctionPsi(const char* filename = "AnalysisResults.root", 
+                               Double_t psiMin = 0., Double_t psiMax = 7.5) {
+  //Macro that draws the correlation functions from the balance function
+  //analysis vs the reaction plane
+  //Author: Panos.Christakoglou@nikhef.nl
+  //Load the PWG2ebye library
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libEventMixing.so");
+  gSystem->Load("libCORRFW.so");
+  gSystem->Load("libPWGTools.so");
+  gSystem->Load("libPWGCFebye.so");
+
+  //Prepare the objects and return them
+  TList *list = GetListOfObjects(filename);
+  if(!list) {
+    Printf("The TList object was not created");
+    return;
+  }
+  else 
+    draw(list,psiMin,psiMax);
+}
+
+//______________________________________________________//
+TList *GetListOfObjects(const char* filename) {
+  //Get the TList objects (QA, bf, bf shuffled)
+  TList *listQA = 0x0;
+  TList *listBF = 0x0;
+  TList *listBFShuffling = 0x0;
+  
+  //Open the file
+  TFile *f = TFile::Open(filename);
+  if((!f)||(!f->IsOpen())) {
+    Printf("The file %s is not found. Aborting...",filename);
+    return listBF;
+  }
+  //f->ls();
+  
+  TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionAnalysis"));
+  if(!dir) {   
+    Printf("The TDirectoryFile is not found. Aborting...",filename);
+    return listBF;
+  }
+  //dir->ls();
+  
+  TString listBFName = "listBF_0-100_V0M_vZ10.0_DCAxy-1.0_DCAz-1.0_Pt0.3-5.0_Eta-0.8-0.8_Chi-1.0_nClus-1_Bit1_withCentralTrigger";
+  listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
+  //listBF->ls();
+
+  //Get the histograms
+  TString histoName = "fHistPV0M";
+  AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));  
+  if(!fHistP) {
+    Printf("fHistP %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistP->FillParent(); fHistP->DeleteContainers();
+
+  histoName = "fHistNV0M"; 
+  AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistN) {
+    Printf("fHistN %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistN->FillParent(); fHistN->DeleteContainers();
+    
+  histoName = "fHistPNV0M"; 
+  AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistPN) {
+    Printf("fHistPN %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistPN->FillParent(); fHistPN->DeleteContainers();
+  
+  histoName = "fHistNPV0M";
+  AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistNP) {
+    Printf("fHistNP %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistNP->FillParent(); fHistNP->DeleteContainers();
+
+  histoName = "fHistPPV0M";
+  AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistPP) {
+    Printf("fHistPP %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistPP->FillParent(); fHistPP->DeleteContainers();
+
+  histoName = "fHistNNV0M";
+  AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+  if(!fHistNN) {
+    Printf("fHistNN %s not found!!!",histoName.Data());
+    break;
+  }
+  fHistNN->FillParent(); fHistNN->DeleteContainers();
+  
+  return listBF;
+}
+
+//______________________________________________________//
+void draw(TList *list, Double_t psiMin, Double_t psiMax) {
+  //Draws the correlation functions for every centrality bin
+  //(+-), (-+), (++), (--)
+  gROOT->LoadMacro("~/SetPlotStyle.C");
+  SetPlotStyle();
+  gStyle->SetPalette(1,0);
+  
+  AliTHn *hP = NULL;
+  AliTHn *hN = NULL;
+  AliTHn *hPN = NULL;
+  AliTHn *hNP = NULL;
+  AliTHn *hPP = NULL;
+  AliTHn *hNN = NULL;
+  
+  hP = (AliTHn*) list->FindObject("fHistPV0M");
+  hN = (AliTHn*) list->FindObject("fHistNV0M");
+  hPN = (AliTHn*) list->FindObject("fHistPNV0M");
+  hNP = (AliTHn*) list->FindObject("fHistNPV0M");
+  hPP = (AliTHn*) list->FindObject("fHistPPV0M");
+  hNN = (AliTHn*) list->FindObject("fHistNNV0M");
+
+  //Create the AliBalancePsi object and fill it with the AliTHn objects
+  AliBalancePsi *b = new AliBalancePsi();
+  b->SetHistNp(hP);
+  b->SetHistNn(hN);
+  b->SetHistNpn(hPN);
+  b->SetHistNnp(hNP);
+  b->SetHistNpp(hPP);
+  b->SetHistNnn(hNN);
+  TH2D *gHistPN[numberOfCentralityBins];
+  TH2D *gHistNP[numberOfCentralityBins];
+  TH2D *gHistPP[numberOfCentralityBins];
+  TH2D *gHistNN[numberOfCentralityBins];
+  
+  TCanvas *cPN[numberOfCentralityBins];
+  TCanvas *cNP[numberOfCentralityBins];
+  TCanvas *cPP[numberOfCentralityBins];
+  TCanvas *cNN[numberOfCentralityBins];
+  TString histoTitle, pngName;
+  //loop over the centrality bins
+  for(Int_t iCentralityBin = 0; iCentralityBin < numberOfCentralityBins; iCentralityBin++) {
+    histoTitle = "(+-) | Centrality: "; 
+    histoTitle += centralityArray[iCentralityBin]; 
+    histoTitle += "%";
+    histoTitle += " | "; histoTitle += psiMin; 
+    histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax;
+    gHistPN[iCentralityBin] = b->GetCorrelationFunctionPN(gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax);
+    gHistPN[iCentralityBin]->GetYaxis()->SetTitleOffset(1.5);
+    gHistPN[iCentralityBin]->SetTitle(histoTitle.Data());
+    cPN[iCentralityBin] = new TCanvas(histoTitle.Data(),"",0,0,400,400);
+    cPN[iCentralityBin]->SetFillColor(10); 
+    cPN[iCentralityBin]->SetHighLightColor(10);
+    gHistPN[iCentralityBin]->DrawCopy("lego2");
+    pngName = "DeltaPhiDeltaEta.Centrality"; 
+    pngName += centralityArray[iCentralityBin]; 
+    pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
+    pngName += ".PositiveNegative.png";
+    cPN[iCentralityBin]->SaveAs(pngName.Data());
+
+    histoTitle = "(-+) | Centrality: "; 
+    histoTitle += centralityArray[iCentralityBin]; 
+    histoTitle += "%";
+    histoTitle += " | "; histoTitle += psiMin; 
+    histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax;
+   gHistNP[iCentralityBin] = b->GetCorrelationFunctionNP(gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax);
+    gHistNP[iCentralityBin]->GetYaxis()->SetTitleOffset(1.5);
+    gHistNP[iCentralityBin]->SetTitle(histoTitle.Data());
+    cNP[iCentralityBin] = new TCanvas(histoTitle.Data(),"",400,0,400,400);
+    cNP[iCentralityBin]->SetFillColor(10); 
+    cNP[iCentralityBin]->SetHighLightColor(10);
+    gHistNP[iCentralityBin]->DrawCopy("lego2");
+    pngName = "DeltaPhiDeltaEta.Centrality"; 
+    pngName += centralityArray[iCentralityBin]; 
+    pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
+    pngName += ".NegativePositive.png";
+    cNP[iCentralityBin]->SaveAs(pngName.Data());
+
+    histoTitle = "(++) | Centrality: "; 
+    histoTitle += centralityArray[iCentralityBin]; 
+    histoTitle += "%";
+    histoTitle += " | "; histoTitle += psiMin; 
+    histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax;
+    gHistPP[iCentralityBin] = b->GetCorrelationFunctionPP(gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax);
+    gHistPP[iCentralityBin]->GetYaxis()->SetTitleOffset(1.5);
+    gHistPP[iCentralityBin]->SetTitle(histoTitle.Data());
+    cPP[iCentralityBin] = new TCanvas(histoTitle.Data(),"",0,400,400,400);
+    cPP[iCentralityBin]->SetFillColor(10); 
+    cPP[iCentralityBin]->SetHighLightColor(10);
+    gHistPP[iCentralityBin]->DrawCopy("lego2");
+    pngName = "DeltaPhiDeltaEta.Centrality"; 
+    pngName += centralityArray[iCentralityBin]; 
+    pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
+    pngName += ".PositivePositive.png";
+    cPP[iCentralityBin]->SaveAs(pngName.Data());
+
+    histoTitle = "(--) | Centrality: "; 
+    histoTitle += centralityArray[iCentralityBin]; 
+    histoTitle += "%";
+    histoTitle += " | "; histoTitle += psiMin; 
+    histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax;
+    gHistNN[iCentralityBin] = b->GetCorrelationFunctionNN(gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax);
+    gHistNN[iCentralityBin]->GetYaxis()->SetTitleOffset(1.5);
+    gHistNN[iCentralityBin]->SetTitle(histoTitle.Data());
+    cNN[iCentralityBin] = new TCanvas(histoTitle.Data(),"",400,400,400,400);
+    cNN[iCentralityBin]->SetFillColor(10); 
+    cNN[iCentralityBin]->SetHighLightColor(10);
+    gHistNN[iCentralityBin]->DrawCopy("lego2");
+    pngName = "DeltaPhiDeltaEta.Centrality"; 
+    pngName += centralityArray[iCentralityBin]; 
+    pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
+    pngName += ".NegativeNegative.png";
+    cNN[iCentralityBin]->SaveAs(pngName.Data());
+  }//end of loop over centralities
+}
+